【StanとRでベイズ統計モデリング】chapter4練習問題
参考にしたサイト
- GitHub - MatsuuraKentaro/RStanBook
- p値を計算したくなる検定の数々を試しにStanによるベイジアンモデリングで代替してみた - 六本木で働くデータサイエンティストのブログ
- Stan超初心者入門
練習問題
- t検定に相当することをStanで行う
- グループY1の平均とグループY2の平均に差があるかを判断せよ
元データの生成
set.seed(123) N1 <- 30 N2 <- 20 Y1 <- rnorm(n=N1, mean=0, sd=5) Y2 <- rnorm(n=N2, mean=1, sd=4)
中身の確認
> Y1 [1] -2.8023782 -1.1508874 7.7935416 0.3525420 0.6464387 8.5753249 [7] 2.3045810 -6.3253062 -3.4342643 -2.2283099 6.1204090 1.7990691 [13] 2.0038573 0.5534136 -2.7792057 8.9345657 2.4892524 -9.8330858 [19] 3.5067795 -2.3639570 -5.3391185 -1.0898746 -5.1300222 -3.6444561 [25] -3.1251963 -8.4334666 4.1889352 0.7668656 -5.6906847 6.2690746 > Y2 [1] 2.7058569 -0.1802859 4.5805026 4.5125340 4.2863243 3.7545610 [7] 3.2156706 0.7523532 -0.2238507 -0.5218840 -1.7788279 0.1683309 [13] -4.0615854 9.6758239 5.8318480 -3.4924343 -0.6115393 -0.8666214 [19] 4.1198605 0.6665237
はい。
練習問題を解く前に普通にRでt検定してみましょう。
> t.test(Y1, Y2, var.equal = T) #Student Two Sample t-test data: Y1 and Y2 t = -1.4772, df = 48, p-value = 0.1461 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.3967325 0.6723788 sample estimates: mean of x mean of y -0.2355188 1.6266580 > t.test(Y1, Y2, var.equal = F) #Welch Welch Two Sample t-test data: Y1 and Y2 t = -1.5884, df = 47.914, p-value = 0.1188 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.2195205 0.4951669 sample estimates: mean of x mean of y -0.2355188 1.6266580
等分散を仮定した場合とそうでない場合の両方で、差があるとは認められないという結果でした(95%有意水準)。
p値は0.14と0.12。これは練習問題で解いた結果の確認用です。
それでは練習問題解いていきます。もう眠くなってきたので淡々といきますよ。
Q1 各グループに差が認められるか、おおよそ把握するために図を書け
par(mfrow=c(2,1)) hist(Y1, breaks=seq(-15,15,2)) hist(Y2, breaks=seq(-15,15,2))
うん。見た感じ差はありそうだぞ。
Q3 StanでQ2のモデルファイルを作成して実行せよ(Q4でR側でMCMCサンプルの操作をするので、generated quantitiesブロックは使わない)
これは普通に、N1, N2, Y1[N1], Y2[N2] を渡して、σ,μ1,μ2 を推定させるので↓こんな感じ。
data { int N1; int N2; real Y1[N1]; real Y2[N2]; } parameters { real<lower=0> sigma; real mu1; real mu2; } model { for (n1 in 1:N1) { Y1[n1] ~ normal(mu1, sigma); } for (n2 in 1:N2) { Y2[n2] ~ normal(mu2, sigma); } }
※最終行が空白じゃないとエラーになる。なぜに?
これを chap4_practice_1.stan という名前にして保存。
↓のRコードで実行する。
data <- list(N1=length(Y1), N2=length(Y2), Y1=Y1, Y2=Y2) fit <- stan(file = "chap4_practice_1.stan", data=data)
Q4 得られたMCMCサンプルからR側でProb[μ1 < μ2] を計算せよ(ヒント:MCMCサンプルの中でμ1<μ2となる回数をカウントして、全サンプル数で割る)
ヒントのとおりにやるのみ
ms <- rstan::extract(fit) df <- data.frame( mu1=ms$mu1, mu2=ms$mu2, diff=ms$mu2-ms$mu1 ) prob <- sum(ifelse(df$diff>0, 1, 0)) / nrow(ms$mu1) prob
結果は
> prob [1] 0.93225
Q5 グループごとの標準偏差が異なる場合をモデル式で表現せよ。また同様にProb[μ1 < μ2]を計算せよ
Q2で共通でおいていたσを、σ1,σ2 にするだけなので、
Stanコードもsigma を sigma1, sigma2 にするだけなので
data { int N1; int N2; real Y1[N1]; real Y2[N2]; } parameters { real<lower=0> sigma1; real<lower=0> sigma2; real mu1; real mu2; } model { for (n1 in 1:N1) { Y1[n1] ~ normal(mu1, sigma1); } for (n2 in 1:N2) { Y2[n2] ~ normal(mu2, sigma2); } }
これを chap4_practice_2.stan という名前にして保存。
↓のRコードで実行する。
data <- list(N1=length(Y1), N2=length(Y2), Y1=Y1, Y2=Y2) fit <- stan(file = "chap4_practice_2.stan", data=data) ms <- rstan::extract(fit) df <- data.frame( mu1=ms$mu1, mu2=ms$mu2, diff=ms$mu2-ms$mu1 ) prob <- sum(ifelse(df$diff>0, 1, 0)) / nrow(ms$mu1) prob
結果は
> prob [1] 0.93175
練習問題はこれで終了。
モデル式書くところでわりと手間取った・・・。まだt検定すら十分に理解できてないのかと落ち込みました。