『岩波データサイエンス Vol.1』を読んでいます。
前からベイズ統計に興味があったので、その勉強内容をまとめます。
今回は、『階層ベイズ最初の一歩』の内容をトレースしてみます。
問題設定
日本国内の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()
以下、パラメーターの事後分布。