R_Bayes_jags
授業用.
_ 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店
応用編を,BRugs パッケージではなく,JAGSを使って実行してみる.
p.134の分散分析.
課題は,大教室で受けた講義の期末試験を同じ教室で受ける場合と,試験だけ小教室で受ける二つのパターンと,さらには,その逆の二つのパターンで,成績に有意に差があるかどうかを調べるという分散分析である.
rjags パッケージの読み込み
library(rjags)
まずデータは次のように用意.これを jagsData.txt として用意しておく.もっとも『マルコフ連鎖モンテカルロ法』付録5-12_分散分析フォルダにあるdata.txtをそのまま読み込んでも良い(ただし「 '21' 行目には,4 個の要素がありません」 と警告が表示される.)
y Lect Test LT 22 0 0 1 15 0 0 1 20 0 0 1 17 0 0 1 16 0 0 1 5 1 0 0 以下略
このデータを次のように読み込む.
d <- read.table("jagsData.txt", header = T)
初期値 list の準備
inits <-list(mu=1, a=0, b=0, c=0, tau.e=1)
モデルファイルjagsModel.txt
model { for(n in 1:20){ y[n] ~ dnorm(theta[n], tau.e) theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n] } mu~dnorm(0.0,1.0E-4) a~dnorm(0.0,1.0E-4) b~dnorm(0.0,1.0E-4) c~dnorm(0.0,1.0E-4) tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e) }
モデルの定義
m <- jags.model( file = "jagsModel.txt", data = d, inits = list(inits, inits, inits), nchain = 3 )
burn-in を行う
update(m, 1000)
MCMC 計算で事後分布からサンプリング,その結果をうけとる
x <- coda.samples( m, c("mu", "a", "b", "c", "tau.e"), thin = 100, n.iter = 20000 )
class(x) [1] "mcmc.list"
結果を見てみる
summary(x) Iterations = 1100:21000 Thinning interval = 100 Number of chains = 3 Sample size per chain = 200 1. Empirical mean and standard deviation for each variable, plus standard error of the mean: Mean SD Naive SE Time-series SE a -0.9706 1.51672 0.061920 0.062622 b -1.0420 1.50942 0.061622 0.072638 c 13.0763 1.48156 0.060485 0.058367 mu 5.0077 1.52229 0.062147 0.062471 tau.e 0.1066 0.03818 0.001559 0.001535 2. Quantiles for each variable: 2.5% 25% 50% 75% 97.5% a -3.79705 -2.03311 -0.9784 0.0934 1.8935 b -4.02259 -2.00015 -1.0789 -0.1027 2.0586 c 10.27393 12.07503 13.0831 14.0117 15.8806 mu 1.87623 3.99396 5.1173 5.9961 7.8435 tau.e 0.04677 0.07727 0.1008 0.1297 0.2024
図を描く
par(ask=TRUE) # パラメータが5個あるので一枚では収まらない plot(x)
上の実行例ではモデルファイルにデータ長を 20 をハードコーティングしているが,これを変数,例えば N に変更するなら,例えばデータ list を準備して
list.data <- list( y = d$y, Lect = d$Lect, Test = d$Test, LT = d$LT, N = nrow(d) )
自明だが,モデルファイルは次のように修正
model { for(n in 1:N){ y[n] ~ dnorm(theta[n], tau.e) theta[n] <- mu + a*Lect[n] + b*Test[n] + c*LT[n] } mu~dnorm(0.0,1.0E-4) a~dnorm(0.0,1.0E-4) b~dnorm(0.0,1.0E-4) c~dnorm(0.0,1.0E-4) tau.e~dgamma(1.0E-3,1.0E-3); sigma.e <- 1.0/sqrt(tau.e) }
このリストをdataに指定して実行
m <- jags.model( file = "jagsModel.txt", data = list.data, inits = list(inits, inits, inits), nchain = 3 )
_ データブロックとモデルブロック
JAGSのモデルブロックで,例えばデータブロックで初期化した X を以下のように右辺に置くと
X[i, j, k, 1:L] ~ dmulti( p[i, j, k, 1:L], n[i, j, k] ) n[i, j, k] <- sum(X[i, j, k, ] )
Unable to resolve parameter n[1,1,1] (one of its ancestors may be undefined)
と怒られる.2行目の処理はデータブロックに移さないといけないらしい. JAGSマニュアルver1.0.3の7.0.5節 Data transformations にもそうある.
_ varは予約語
モデルにうっかり使うとエラーになるので,var0 などとする.
syntax error, unexpected VAR, expecting NAME or FUNC or FOR or '}' 以下にエラー jags.model...
_ DIC
こんな感じで算出すればよい
x <- dic.samples( m, n.iter = 10000, type = "pD" ) ## x$deviance ## x$penalty x1 <- as.mcmc(x) summary(x1)
_ WinBUGSのstructure()関数について
添付ファイル: jagsData.txt 1843件
[詳細]
jagsModel.txt 1815件
[詳細]
Last-modified: 2009-08-13 (木) 16:10:24 (5532d)