イントロダクション
MCMC（マルコフ連鎖モンテカルロ法）は、複雑な確率分布（多くの場合はベイズ統計の事後分布）から直接サンプルを生成するための手法である。回帰係数の解析解が求まらないようなモデルでも、MCMCを回せばそれらしいサンプル列が出力される。しかし、ここに大きな落とし穴がある。MCMCは「計算が終わった」ことと「正しい分布に収束した」ことを区別してくれない。途中で止めても、何かしらの数列は出力され続けるからである。
出力されたサンプルが、本当に目的の分布を代表しているのか。それとも、まだ特定の領域に偏ったまま探索が終わっていないのか。これを判断するために使われるのが、本記事で紹介するMCMCの評価指標（診断指標）である。
トレースプロット：最初に見るべきグラフ
最も基本的な診断は、サンプリングの反復回数を横軸に、パラメータの値を縦軸にとったトレースプロットを目で見ることである。よく混合しているチェーンは、特定の値の周りを高速に、毛虫のようにランダムに行き来する。逆に混合が悪いチェーンは、同じ値の周辺に長く留まったり、ゆっくりとしか値が変化しなかったりする。
図1: 左は値の周りを高速に往復しており、探索がうまくいっている。右は同じ値に長く留まってから別の値に飛んでおり、収束や混合に問題がある兆候。
トレースプロットは直感的だが、あくまで目視による確認にすぎない。そこで、数値として定量化する指標が必要になる。
R-hat（Gelman-Rubin統計量）
実務では、初期値を変えた複数のチェーン（例えば4本）を並列に走らせるのが一般的である。もし全てのチェーンが同じ分布に収束していれば、チェーン内のばらつきとチェーン間のばらつきはほぼ同じになるはずである。逆に、チェーンごとに違う場所をうろついているなら、チェーン間のばらつきの方が大きくなる。
この考え方を指標にしたのが$\widehat{R}$（R-hat、potential scale reduction factor）である。
$\widehat{R} = \sqrt{ {\widehat{V} \over W} }$
ここで $W$ はチェーン内分散の平均、$\widehat{V}$ はチェーン内分散とチェーン間分散を組み合わせた分散の推定値である。全チェーンが同じ分布から来ているなら $\widehat{R}$ は1に近づき、収束していなければ1より大きくなる。
実務上の目安として、Vehtari et al.（2021）は、少なくとも4本のチェーンを走らせた上で、$\widehat{R} \lt 1.01$ のときにだけそのサンプルを使うことを推奨している。これは Gelman and Rubin（1992） の推奨よりも大幅に厳しい基準で、従来は $\widehat{R} \lt 1.1$ が収束の目安として扱われていた。同論文は、1.01という数値はあくまで最初のチェックとしての目標値であり、用途に応じて調整すべきだとも述べている。
有効サンプルサイズ（ESS）
MCMCが生成するサンプルは、隣り合う反復同士が似た値になりやすい（自己相関を持つ）。そのため、たとえ1万個サンプルを取っても、独立なサンプル1万個分の情報量があるとは限らない。有効サンプルサイズ（Effective Sample Size, ESS）は、「自己相関を考慮すると、実質的に独立なサンプルが何個分に相当するか」を表す指標である。
$\mathrm{ESS} = \dfrac{n}{1 + 2\sum_{k=1}^{\infty}\rho_k}$
ここで $n$ は全サンプル数、$\rho_k$ はラグ $k$ の自己相関である。自己相関がすぐに0に近づくチェーンほど、ESSはサンプル数 $n$ に近づく。自己相関がなかなか減衰しないチェーンほど、ESSは $n$ よりずっと小さくなる。
図2: 緑は自己相関がすぐに小さくなり、少ないサンプルでも実質的な情報量が多い。赤は自己相関が長く残り、同じサンプル数でも実質的な情報量（ESS）が少ない。
なお近年では、事後分布の中心付近だけでなく、裾（分位点）の推定精度も別に確認するため、bulk-ESS（分布の中心付近の推定に関するESS）とtail-ESS（分布の端の推定に関するESS）を分けて報告することが推奨されている。目安として、Vehtari et al.（2021）は、順位正規化したESSが400を超えることを推奨している（安定したMCSEの推定に通常十分な値として、経験とシミュレーションから選んだ数字だと説明されている）。tail-ESSについても、同じ400が推奨水準として示されている。
モンテカルロ標準誤差（MCSE）
事後分布の平均を推定するとき、私たちが実際に使えるのは有限個のMCMCサンプルから計算した標本平均である。この標本平均には、有限サンプルであることに起因する誤差が乗る。この誤差の大きさを表すのがモンテカルロ標準誤差（Monte Carlo Standard Error, MCSE）である。
$\mathrm{MCSE} = {\mathrm{SD}(\theta\mid y) \over \sqrt{\mathrm{ESS}}}$
ESSが大きいほどMCSEは小さくなり、報告する推定値（事後平均など）の桁数をどこまで信じてよいかの目安になる。「小数点以下2桁まで報告しているが、MCSEを計算するとその桁は誤差に埋もれている」といったことは実務でもよく起きる。
ダイバージェンス（HMC・NUTS特有の診断）
StanやPyMCなどで使われるハミルトニアンモンテカルロ（HMC）やNUTSというアルゴリズムには、他にはない専用の診断指標がある。それがダイバージェンス（divergent transitions）である。事後分布の中に極端に曲率の高い領域（急に狭くなる谷など）があると、数値積分がその領域をうまく再現できず、警告として記録される。
$\widehat{R}$やESSが良好でも、ダイバージェンスが多数報告されている場合は要注意である。事後分布の特定の領域が正しく探索できていない可能性があり、モデルの再パラメータ化（例えば階層モデルの非中心化パラメータ化）などの対処が必要になることが多い。
実務でのチェックリスト
| 指標 | 見るポイント | 目安 |
|---|---|---|
| トレースプロット | 複数チェーンが同じ範囲を素早く往復しているか | 目視で「毛虫」状になっているか |
| $\widehat{R}$ | チェーン間・チェーン内のばらつきが同程度か | 1.01未満 |
| ESS（bulk / tail） | 自己相関を考慮した実質的なサンプル数 | それぞれ400以上 |
| MCSE | 推定値のどの桁まで信頼できるか | 報告する桁数より十分小さいか |
| ダイバージェンス（HMC/NUTS） | 数値積分が失敗した回数 | 0に近いことが望ましい |
まとめ
- MCMCは「止めた時点」と「収束した時点」が一致する保証がなく、専用の診断指標で確認する必要がある。
- トレースプロットは最も基本的な目視診断で、複数チェーンの混合具合を確認する。
- $\widehat{R}$は、チェーン間とチェーン内のばらつきを比較し、複数チェーンが同じ分布に収束しているかを数値化する。
- ESSは自己相関を考慮した実質的なサンプル数であり、MCSE（推定値の誤差）の計算にも使われる。
- HMC/NUTS系のサンプラーでは、ダイバージェンスという専用の警告も併せて確認する。
参考文献
- Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R̂ for assessing convergence of MCMC (with discussion). Bayesian Analysis, 16(2), 667–718. https://doi.org/10.1214/20-BA1221
- Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472. https://doi.org/10.1214/ss/1177011136
