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 で描画されるので実際きれいである.上記の参考ページにはもっと多くの画像がある.
ちなみに 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 引数を追加したので, 「マッチしない」条件の記述が簡単にできる.
2016/10/09: data.table の関数を使ってるのにパッケージを読み込んでなかったのを修正
2017/4/15: サンプルを thinning した場合を考慮していなかったので修正
gistcdd93fe3937d2c39a9f767845e6c7534
上記の ggs2() 関数を ggs() の代わりに使う.
また, デフォルトだと ggmcmc() 関数は全てのグラフを1つの pdf に出力してしまう. plot= 引数に対して出力するものを "ggs_histogram" などと指定できる (もしくは ggs_histogram() を直接呼び出す) ので, そちらも活用すると良い.
