読者です 読者をやめる 読者になる 読者になる

ill-identified diary

所属組織の業務や見解とは一切無関係なアフィリエイト付きメモ帳。所属とは関係ないけどここを見て所属先にも興味を持っていただけると喜びます。

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

概要

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

参考: Using the ggmcmc package

内容

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

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

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

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

ちなみに ggmcmcggplot2 のほか, dplyrtidyr に依存しているので呼び出した時点でこれらも読み込まれる.

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

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

これらの問題を解決するために, 自分なりに改変した関数 ggs2 を書いておく. rstan パッケージの stanfit, rjagsmcmc.list オブジェクト, runjagsrunjags オブジェクトに対応している. 表示パラメータを指定する場合も, 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() を直接呼び出す) ので, そちらも活用すると良い.


*1:定番のモデルならいいのだろうが, パラメータが膨大にあるモデルを作る必要があって, そういう場合は1つの画像に大量のグラフが凝集されるので大変見づらかった.

*2:そして自己相関の計算の際に「標準偏差がゼロになりました」と大量の警告を発する……

*3:100MB未満程度なら早いが, 500MB を超えたものだと自分で dplyr で処理した場合と比べて倍近く遅くなった.