勉強などのメモ

勉強用のメモ

Chapter4 練習問題

t検定に相当することをStanで行う。以下のRコードで生成されたデータを用いて、グループ1(Y1)の平均パラメーターμ1とグループ2(Y2)の平均パラメーターμ2に差があるか、Prob[μ1<μ2](μ1<μ2である確率)を求めることで判断したい

set.seed(123)
N1 <- 30
N2 <- 20
Y1 <- rnorm(n=N1, mean=0, sd=5)
Y2 <- rnorm(n=N2, mean=1, sd=4)

(1)各グループの値に差が認められるか、おおよそ把握するために図を描け

hist(Y1, breaks=seq(-10,10,2))

f:id:prism0081:20170902015501j:plain

hist(Y2, breaks=seq(-10,10,2))

f:id:prism0081:20170902015508j:plain

とりあえず簡単にヒストグラムで確認。まあ違うっぽい。

(2)各グループで標準偏差が等しいと仮定してモデル式を書け

Y_{1}[n_{1}] \sim Normal(\mu_{1}, \sigma)   (n_{1}=1,\cdots,30)
Y_{2}[n_{2}] \sim Normal(\mu_{2}, \sigma)   (n_{2}=1,\cdots,20)

(3)Stanで(2)のモデルファイルを作成して実行せよ

Stanのコードは以下のとおり。

data{
  int N1;
  int N2;
  real Y1[N1];
  real Y2[N2];
}

parameters{
  real myu1;
  real myu2;
  real<lower=0> sig;
}

model{
  for(n in 1:N1)
    Y1[n] ~ normal(myu1, sig);
  for(n in 1:N2)
    Y2[n] ~ normal(myu2, sig);
}

これをRで実行する。

library(rstan)
model1 <- stan(file="Rstan4-4-1.stan", data=list(N1=N1,N2=N2,Y1=Y1,Y2=Y2))

トレースプロットをみてみる。

traceplot(model1)

f:id:prism0081:20170902021821j:plain

問題なさそう。

(4)得られたMCMCサンプルからR側でProb[μ1<μ2]を計算せよ

ベタ打ちだけどこれでいいのかな。

myu1 <- rstan::extract(model1)$myu1

myu2 <- rstan::extract(model1)$myu2

prob <- sum(ifelse(myu1<myu2,1,0)) / nrow(myu1)

prob

[1] 0.92175

(5)グループごとの標準偏差が異なる場合をモデル式で表現せよ。同様にProb[μ1<μ2]を計算せよ

モデル式は以下。

Y_{1}[n_{1}] \sim Normal(\mu_{1}, \sigma_{1})   (n_{1}=1,\cdots,30)
Y_{2}[n_{2}] \sim Normal(\mu_{2}, \sigma_{2})   (n_{2}=1,\cdots,20)

σに添え字がついただけ。

Stanのコードも微修正の範囲。

data{
  int N1;
  int N2;
  real Y1[N1];
  real Y2[N2];
}

parameters{
  real myu1;
  real myu2;
  real<lower=0> sig1;
  real<lower=0> sig2;
}

model{
  for(n in 1:N1)
    Y1[n] ~ normal(myu1, sig1);
  for(n in 1:N2)
    Y2[n] ~ normal(myu2, sig2);
}

これを実行する。

model2 <- stan(file="Rstan4-4-2.stan", data=list(N1=N1, N2=N2, Y1=Y1,Y2=Y2))
traceplot(model2)

f:id:prism0081:20170902022911j:plain

myu12 <- rstan::extract(model2)$myu1

myu22 <- rstan::extract(model2)$myu2

prob2 <- sum(ifelse(myu12<myu22,1,0)) / nrow(myu12)

変数名適当すぎるけど、結果は・・・

prob2

[1] 0.93625

これであっているのかわからないけど、とりあえず問題なく処理はできた。