ベイズ統計モデルを勉強している

『岩波データサイエンス Vol.1』を読んでいます。

https://www.amazon.co.jp/%E5%B2%A9%E6%B3%A2%E3%83%87%E3%83%BC%E3%82%BF%E3%82%B5%E3%82%A4%E3%82%A8%E3%83%B3%E3%82%B9-Vol-1-%E5%B2%A9%E6%B3%A2%E3%83%87%E3%83%BC%E3%82%BF%E3%82%B5%E3%82%A4%E3%82%A8%E3%83%B3%E3%82%B9%E5%88%8A%E8%A1%8C%E5%A7%94%E5%93%A1%E4%BC%9A/dp/4000298518

前からベイズ統計に興味があったので、その勉強内容をまとめます。

今回は、『階層ベイズ最初の一歩』の内容をトレースしてみます。

問題設定

日本国内の10県で、県ごとに異なる2種類の給食タイプが身長ののびに影響を与えるか調査。

半分の5県では普通の給食、もう半分が異なるタイプの給食を与えられる。

現在の身長と、1年後の身長を測定し、特殊なタイプの給食の効果を推定する。

正解としては、特殊なタイプの給食が身長に与える効果は無い。

以下に解析対象のデータを示す。

測定回目の0が現在の身長測定結果を、1が1年後の測定結果を表す。

給食タイプの0が普通の給食を、1が特殊な給食を表す。

つまり、県A,B,D,E,Iが特殊な給食を与えられる県。

モデル1 直線あてはめ

後日まとめます。。。

モデル2 ベイズモデル

後日まとめます。。。

以下、コードと結果なぶり書き。

library(rjags) # R と JAGS をつなげる package
library(R2WinBUGS) # write.model() 使うため
model.bugs <- function()
{
  for (i in 1:N) {
    mean.Y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta[1] + (beta[2] + beta[3] * X[i]) * Age[i]
  }
  for (k in 1:N.beta) {
    beta[k] ~ dunif(-10000, 10000)
  }
  tau <- 1/(sd * sd)
  sd ~ dunif(0.00000E+00, 10000)
}
file.model <- "model.bug.txt"
write.model(model.bugs, file.model)

load("data.RData") # データ d の読みこみ
N.beta <- 3 # beta[1], beta[2], beta[3]
list.data <- list( # データ
  mean.Y = d$mean.Y,
  X = d$X, Age = d$Age,
  N = nrow(d), N.beta = N.beta
)
# パラメーターの初期値
list.inits <- list(beta = rep(0, N.beta), sd = 1)

# 事後分布からのサンプリングの詳細
n.burnin <- 1000
n.chain <- 3
n.thin <- 10
n.iter <- n.thin * 3000

model <- jags.model(
  file = file.model, data = list.data,
  inits = list.inits, n.chain = n.chain
)
update(model, n.burnin) # burn in

# 推定結果を post.mcmc.list に格納する
post.mcmc.list <- coda.samples(
  model = model,
  variable.names = c("mu", names(list.inits)),
  n.iter = n.iter,
  thin = n.thin
)

pdf("plot1.pdf") 
plot(post.mcmc.list) 
dev.off()

以下、パラメーターの事後分布。