緑本メモ: MCMCサンプルをJAGSとRで実行 - ランダム効果1個版 (10.3.1節)

緑本サイトからdata7a.csv と model.bug.txt をダウンロードする。
9章の例題をJAGSでやるメモの記載の通り、ダウンロードした model.bug.txt の文末にセミコロン; を追加しておく。

このファイル(データとモデル)をもとに下のようにして、パラメータの推定・可視化と、MCMCが収束しているかの評価値をそれぞれ出力する。

パラメータの推定と可視化

下のRコードで、切片βと個体差のばらつきsのパラメータをJAGSによるMCMCによって推定する。

library('rjags')

d <- read.csv("data7a.csv")
filename <- 'model.bug.txt'
N <- nrow(d)

list.data <- list(
  N = N,
  Y = d$y
)

inits <- list(
  beta  = c(0),
  r     = rnorm(N, 0, 0.1),
  s     = 1
)

m <- jags.model(
  file = filename,
  data = list.data,
  inits = list(inits, inits, inits),
  n.chain = 3
)

update(m, 100)

post.list <- coda.samples(
  m,
  c("beta", "s"),
  thin = 10, n.iter = 10100
)
summary(post.list)

plot(post.list)

summary(..) で出力されるβ,sの推定値：

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean     SD Naive SE Time-series SE
beta 0.05688 0.3395 0.006168       0.011131
s    3.03808 0.3723 0.006764       0.007663

plot で出力されるサンプリングのトレースと、周辺事後分布の近似：

MCMCの収束の評価値の出力

9章の例題をJAGSでやるメモと同様に R2jags で収束評価値 R-hat を出力する：

library('R2jags')

post.jags <- jags(
  data = list.data,
  inits = list(inits, inits, inits),
  parameters.to.save = c("beta", "s"),
  n.iter = 10100,
  model.file = filename,
  n.chains = 3,
  n.thin = 10,
  n.burnin = 100)

print(post.jags) # R-hat

出力結果（右から二番目の列に R-hat）：

Inference for Bugs model at "model.bug.txt", fit using jags,
 3 chains, each with 10100 iterations (first 100 discarded), n.thin = 10
 n.sims = 3000 iterations saved
         mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
beta       0.045   0.332  -0.607  -0.176   0.039   0.262   0.709 1.001  2500
s          3.019   0.364   2.393   2.758   2.996   3.247   3.806 1.001  3000
deviance 223.184  13.582 197.821 213.551 222.686 232.133 251.500 1.001  3000

R-hatが1.0付近で1.1より小さければ収束していると判断。