R_Bayes_jags のバックアップ(No.5) - アールメカブ

アールメカブ


R_Bayes_jags のバックアップ(No.5)


Rの備忘録

授業用.

_ 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店

応用編を,BRugs パッケージではなく,JAGSを使って実行してみる.

p.134の分散分析.

課題は,大教室で受けた講義の期末試験を同じ教室で受ける場合と,試験だけ小教室で受ける二つのパターンと,さらには,その逆の二つのパターンで,成績に有意に差があるかどうかを調べるという分散分析である.

rjags パッケージの読み込み

library(rjags)

まずデータは次のように用意.これを filejagsData.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)

モデルファイルfilejagsModel.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 にもそうある.

_ WinBUGSのstructure()関数について