試しに計算してみたのですが,まだ細部を理解できていないので,検証用にメモを残しておきます.
豊田本2017(『実践ベイズモデリング』)にStan
からWAIC
を計算する手順が解説されていたので,それを参考にして
JAGS
から作ったサンプルでWAIC
を計算してみます.
例によってJAGS
をR
から呼び出すのにパピー4匹本のDBDA2E-utilities.R
を使います.
まずJAGS
用のmodel{ }
内に
for (i in 1:N) {
loglikelihood[i] <- log(dnorm(y[i], mu[i], 1/ (s * s)));
}
などを追加して,推定したパラメータとデータを用いた対数尤度を計算します.
この場合はパラメータはmu
とs
です.mu
は説明変数に回帰しているので添え字[i]
がついています.
確率密度関数は応答関数にあわせて適宜修正しましょう.
coda.samples
の結果を格納します.
mcmcCoda = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nPerChain , thin=thinSteps )
ここから対数尤度を整形して取り出し,looパッケージのwaic関数に渡します.
こちらの記事を参考にさせていただきました.
ito-hi.blog.so-net.ne.jp
library(loo)
loglik1 <- sapply(1:N, function(i)
unlist(mcmcCoda[, paste("loglikelihood[", i, "]", sep = "")]))
waic(loglik1)
これでStanから計算した結果とほぼ同じWAIC
が出力されました.
loo
パッケージのextract_log_lik
関数を使えばもっと簡潔になりそう.