to be continued...

Rや心理統計学の備忘録的な

There are three kinds of lies: lies, damned lies, and statistics.

- 嘘には三種類ある。嘘、大嘘、そして統計だ -
Benjamn Disraeli(19世紀のイギリス首相)

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のときと同じ処理でいける。

以上で終わりだ それだけ

to be continued ...