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))
hist(Y2, breaks=seq(-10,10,2))
とりあえず簡単にヒストグラムで確認。まあ違うっぽい。
(2)各グループで標準偏差が等しいと仮定してモデル式を書け
(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)
問題なさそう。
(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]を計算せよ
モデル式は以下。
σに添え字がついただけ。
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)
myu12 <- rstan::extract(model2)$myu1 myu22 <- rstan::extract(model2)$myu2 prob2 <- sum(ifelse(myu12<myu22,1,0)) / nrow(myu12)
変数名適当すぎるけど、結果は・・・
prob2 [1] 0.93625
これであっているのかわからないけど、とりあえず問題なく処理はできた。