The One with ...

研究用のメモなどを公開しています

JAGSでWAICを計算するためのメモ

試しに計算してみたのですが,まだ細部を理解できていないので,検証用にメモを残しておきます.

豊田本2017(『実践ベイズモデリング』)にStanからWAICを計算する手順が解説されていたので,それを参考にして JAGSから作ったサンプルでWAICを計算してみます.

例によってJAGSRから呼び出すのにパピー4匹本のDBDA2E-utilities.Rを使います.

まずJAGS用のmodel{ }内に

 # Generated quantities
    for (i in 1:N) {
        loglikelihood[i] <- log(dnorm(y[i], mu[i], 1/ (s * s)));
    }

などを追加して,推定したパラメータとデータを用いた対数尤度を計算します. この場合はパラメータはmusです.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関数を使えばもっと簡潔になりそう.