[stan] [jags] ggmcmc でMCMCの事後診断

2019/12/15 追記: 現在は ggmcmcよりもbayesplotのほうがおすすめです
ill-identified.hatenablog.com

概要

これまで, stan などのサンプリング結果を R で処理するのが面倒だと思っていたのだが, いまさら ggmcmc パッケージという便利なものに気づいた. rstan, rjags などの R と連携できるパッケージと組み合わせるとトレースプロットやコレログラム, 事後密度やヒストグラムを簡単に出力してくれるが, 若干使いにくい*1ところもあるので使い方と合わせて改変したものについても言及する.

参考: http://xavier-fim.net/packages/ggmcmc/

内容

ggmcmc の構文は簡単で, 最短で

library(ggmcmc)
ggmcmc(ggs(stan/jagsの出力オブジェクト), file=ファイル名)

だけで出力できる.グラフは ggplot2 で描画されるので実際きれいである.上記の参考ページにはもっと多くの画像がある.

f:id:ill-identified:20161007231714p:plain

ちなみに ggmcmc は ggplot2 のほか, dplyr と tidyr に依存しているので呼び出した時点でこれらも読み込まれる.

ggs() 関数は stanfit や rjag オブジェクトを ggmcmc が読み取れるよう変換するものだが, 使いにくいところもある. 具体的には, そのままだと全てのパラメータに対して計算してしまうので, 例えば情報量規準の評価のために WAICを計算させてみた：Taglibro de H：SSブログのように尤度を計算していた場合, 各観測値ごとに計算した大量の尤度についても律儀にヒストグラムなどを作成してしまう*2.ggs() には計算対象とするパラメータを指定する引数があるのだが, ERE を書くことでしか指定できないので, いまいち使いにくい.

また入力データが大きいと遅くなる*3 のも気になった.

これらの問題を解決するために, 自分なりに改変した関数 ggs2 を書いておく. rstan パッケージの stanfit, rjags の mcmc.list オブジェクト, runjags の runjags オブジェクトに対応している. 表示パラメータを指定する場合も, grepl() の引数をそのまま使うことができる. つまり, ggs() と違い perl の正規表現も使える. また grep() 関数のように invert 引数を追加したので, 「マッチしない」条件の記述が簡単にできる.