RJagsでMCMCやってみた (緑本9章)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

9章の例題をやってみたメモ。
MacなのでWinBUGSは動かないからRJagsで。

インストール

Jags本体と、Rから使うためのRJagsを入れる必要がある。

Jags

Mac用のパッケージ(dmg)を入れた。
homebrewにもJagsはあるようだが、PATHがRJagsが想定しているデフォルトとは違うらしく面倒そう。特にこだわりが無ければパッケージを入れた方が楽。

RJags

普通に install.packages("rjags") するだけ

Rコード

使用するデータは書籍のサポートサイトからダウンロードできる。
同じ人の別記事にあるコードをちょっと変えただけで動いた。

library(rjags)

# 本のサポートサイトからダウンロードしたデータ読み込み
load("d.RData")

# データの準備
list.data <- list(
    N = nrow(d),
    Y = d$y,
    X = d$x,
    Mean.X = mean(d$x)
)

# 初期値の準備
inits <- list(
    beta1 = 0,
    beta2 = 0
)

# R内でのモデルの定義
m <- jags.model(
    file = "model.txt",
    data = list.data,
    inits = list(inits, inits, inits),
    n.chain = 3
)

# burn-in のためのカラまわしMCMC計算
update(m, 100)

# MCMC計算で事後分布からサンプリング
x <- coda.samples(
    m,
    c("beta1", "beta2"),
    thin = 3, n.iter = 1500, n.burnin = 100
)

plot(x)

• 本の例題(WinBUGS)と違うのは、RJags では update(m, 100) がburn-inのための空回しらしい点。n.burnin=100 でも指定できているようだからどちらかが要らないと思う。

BUGSのコード (model.txt)

model
{
    for (i in 1:N) {
        Y[i] ~ dpois(lambda[i]);
        log(lambda[i]) <- beta1 + beta2 * (X[i] - Mean.X);
    }
    beta1 ~ dnorm(0, 1.0E-4);
    beta2 ~ dnorm(0, 1.0E-4);
}

• Jagsではセミコロン必須らしい

プロット結果

参考サイト

• 生態学データ解析 - 本/データ解析のための統計モデリング入門
  • 緑本のサポートサイト。データセットのダウンロードもできる。
• 生態学データ解析 - R で JAGS
  • 同じ人のRJags解説。