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

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

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

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

まずJAGS用のmodel{ }内に

 # Generated quantities
    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関数に渡します．

こちらの記事を参考にさせていただきました．

library(loo)
loglik1 <- sapply(1:N, function(i)
    unlist(mcmcCoda[, paste("loglikelihood[", i, "]", sep = "")]))
waic(loglik1)

これでStanから計算した結果とほぼ同じWAICが出力されました． looパッケージのextract_log_lik関数を使えばもっと簡潔になりそう．

The One with ...