cmdstanrで推定後にWAICを計算
結論
この関数を使おう
※stanファイルの対数尤度の変数名がlog_lik
として推定している事を仮定
waic2 <- function(fit){ # log_lik <- rstan::extract(fit)$log_lik # rtan用 log_lik <- fit$draws("log_lik") %>% as_draws_df() %>% select(contains("log_lik")) # cmdstanr用 lppd <- sum(log(colMeans(exp(log_lik)))) p_waic <- sum(apply(log_lik,2,var)) waic <- -2*lppd+2*p_waic return(list(waic=waic,lppd=lppd,p_waic=p_waic)) }
小噺
Rのcmdstanr
パッケージを使ってMCMCをしたときの推定結果は、rstan
パッケージのそれとは異なります。rstan
で推定した場合はstanfit
オブジェクト、cmdstanr
で推定した場合はCmdStanFit
オブジェクトとなります。
対数尤度からWAIC((Watanabe-Akaike) Widely Applicable Information Criterion)を計算するとき、loo
パッケージを使うこともできますが、loo
パッケージのwaic()
は引数がstanfit
なので、cmdstanrで推定した場合、stanfit
に直してあげる必要があります。
なので、自分はこちらの記事にある自作関数で計算します。
ただ、この記事のwaic2()
もコードの通り、stanfit
を引数としているので、CmdStanFit
用に個人的に直しました。
補足
stanfitオブジェクトから特定の推定結果を抜き出すときrstan::extract()
でパラメータを指定します。この時の出力は2次元配列(パラメータ×mcmcサンプル総数)です。
同様の処理をcmdstanrの時は、fit$draws("log_lik")で行えます。
しかし、cmdstanrだと3次元配列(パラメータ×チェイン×iteration)で出力されます。
なので、as_draws_ds()
でdata.frame型に変換します。この時、2次元データに落とし込めますが、log_likの他に'.chain', '.iteration', '.draw'という変数もついてきます。これが便利かどうかはおいておいて、select()で除外してあげれば、あとはstanfit
のときと同じ処理でいける。
以上で終わりだ それだけ