データ分析のお勉強

メモ代わりにいろいろ書いていきます

【StanとRでベイズ統計モデリング】chapter4練習問題

練習問題

  • 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))

f:id:ritsujiro:20161030111229p:plain

うん。見た感じ差はありそうだぞ。

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

Y1, Y2それぞれモデル式で表すと、

Y_1[N_1]  = μ_1 + ε[N_1]
ε[N_1] \sim N(0, σ)

より、

Y_1[N_1]  \sim N(μ_1, σ)

Y2も同様で

Y_2[N_2]  \sim N(μ_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 にするだけなので、
Y_1[N_1]  \sim N(μ_1, σ_1)
Y_2[N_2]  \sim N(μ_2, σ_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検定すら十分に理解できてないのかと落ち込みました。