アールメカブ

• 追加された行はこの色です。
• 削除された行はこの色です。
• R_Bayes_jags へ行く。
• R_Bayes_jags の差分を削除

[[Rの備忘録]]

授業用．

#contents

* 豊田秀樹『マルコフ連鎖モンテカルロ法』朝倉書店 [#f02738ac]
応用編を，BRugs パッケージではなく，[[JAGS:http://www-fis.iarc.fr/~martyn/software/jags/]]を使って実行してみる． 

p.134の分散分析.

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


rjags パッケージの読み込み
 library(rjags)

まずデータは次のように用意．これを &ref(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)
モデルファイル&ref(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
 )
//## リストとして読み込んだら data = list.data,
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
 )

* データブロックとモデルブロック [#ze9cc2a7]
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 にもそうある．
[[JAGSマニュアル:http://www-fis.iarc.fr/~martyn/software/jags/]]ver1.0.3の7.0.5節 Data transformations にもそうある．
//-- JAGSでは値を固定する母数の初期値として NA を指定するとエラーになる．

* varは予約語 [#d5840e35]
モデルにうっかり使うとエラーになるので，var0 などとする．

 syntax error, unexpected VAR, expecting NAME or FUNC or FOR or '}'
  以下にエラー jags.model...

* DIC [#t798c6d2]
こんな感じで算出すればよい

 x <- dic.samples(
	m,
	n.iter = 10000,
    type = "pD"
 )
 ## x$deviance
 ## x$penalty
 x1 <- as.mcmc(x)
 summary(x1)

* WinBUGSの[[structure()関数について>R_winBUGS_structure]] [#tcc9f4f4]