MCMCの収束診断に用いられるPSRとESSについて前からモヤモヤしていたので(←「不勉強でよく理解できていなかったので」の間接的表現)、このたびメモをとってみた。
以下、BDA3の10.5, 11.4, 11.5からのメモ。
ベイジアン・シミュレーションの目的は、事後分布から独立なドロー\(\theta^s \ \ ( s=1, \ldots, S ) \)の集合を手に入れ、関心ある量をリーズナブルなレベルの正確性で推定することだ。
たいていの場合\(S=100\)くらいあれば十分である。なぜかというと…
\(\theta\)はスカラーで、事後分布は\( N(\mu_\theta, \sigma^2_\theta) \)で近似できるとしよう。これを分析的に算出できないので、\(S\)回のシミュレーション・ドローから得た平均\(\bar{\theta}\)とSD \(s_\theta\)で推定するとしよう。事後平均の推定の正確性はおよそ \( s_\theta/\sqrt{S} \)である。したがって、計算されたパラメータ推定値のトータルのSD(モンテカルロ誤差、つまりシミュレーションのドローの数が有限でしかないことによって生じる不確実性を含む)は\(s_\theta \sqrt{1+1/S}\)である。\(S=100\)ならば\(\sqrt{1+1/S}\)は1.005となり、モンテカルロ誤差は実際の事後分散に由来する不確実性にほとんど何も付け加えていないとわかる。
[このくだり、いい加減に読み飛ばしてしまったんだけど、このメモを読み返して「なんだこれは???」とちょっと混乱した。えーと、\(\hat{Var}(\theta) + \hat{Var}(\bar{\theta}) = s^2_\theta + s^2_\theta/S = s^2_\theta (1+1/S) \)の平方根が\(s_\theta \sqrt{1+1/S}\)だということかしらん]
MCMCのような反復シミュレーションでも、考え方は上と基本的に同じなんだけど、新たな難題が2つある。その1、初期の反復は目標分布というより初期値を反映する。その2, 独立なドローではなく、系列相関がある。
対処は3つある。
- 初期値をパラメータ空間上で散らばらせた複数の系列をシミュレートする。少なくとも2系列は走らせる。
- 関心ある量のすべてについて収束を監視する。系列間分散と系列内分散を比べ、両者が大体同じになるまで走らせる。
- 関心ある量についての事後推測についてのおよその収束が得られるまで時間がかかりすぎる場合は、アルゴリズムを変える。
初期値の影響を取り除くため、初期の反復を捨てるのが良い。捨てる部分をwarm-upという。保守的な選択としては、各系列の前半を捨てる。
収束の監視にあたっては、関心あるパラメータとか事後密度の対数とかのそれぞれについて、複数系列が混合していてかつ定常になっているかどうかに注目する。
ひとつのやりかたはこうだ。チェーンが5本で、warm-upを捨てたあとの反復数が500だとしよう。それぞれを前半と後半に分けると、長さ\(n=250\)のチェーンが\(m=10\)本できる。
estimandはスカラーで\(\psi\), \(j\)本目のチェーンの\(i\)番目のシミュレーションを\(\psi_{ij}\)とする。あるチェーンの平均を\(\bar{\psi}_{.j} = \frac{1}{n} \sum_i \psi_{ij}\), 全平均を\( \bar{\psi}_{..} = \frac{1}{m} \sum_j \bar{\psi}_{.j}\), あるチェーンの分散を\( s^2_j = \frac{1}{n-1} (\psi_{ij} – \bar{\psi}_{.j})^2 \)とする。
系列間分散は $$ B = \frac{n}{m-1} \sum_j (\bar{\psi}_{.j} – \bar{\psi}_{..})^2$$ となる。系列内分散は $$ W = \frac{1}{m} \sum_j s^2_j $$ となる。
estimandの周辺事後分散\( var(\psi|y)\)は、以下のように推定できる。$$ \hat{var}^+ (\psi|y) = \frac{n-1}{n} W + \frac{1}{n} B $$ この量は、定常性の仮定(つまり初期値の分布が目標分布に等しいという仮定)の下で不偏である。ないし、\(n \rightarrow \infty\)の極限において不偏である。クラスタサンプリングの分散推定みたいなものである[←なるほど…]。しかし実際には、(初期値が適切に過分散を持たされていれば)大きめになる。
\(n \rightarrow \infty\)においては\(W\)も\( var(\psi|y)\)の不偏推定量である。実際には、目標分布を全部カバーできないので、小さめになる。
よって \( \frac{ \hat{var}^+ (\psi|y) }{ W }\)は、\(n \rightarrow \infty\)において1となる。そこで、$$ \hat{R} = \sqrt{\frac{ \hat{var}^+ (\psi|y) }{ W }} $$ に注目する(平方根をとったのはそのほうが表記上便利だから)。この量は、\(\psi\)の現在の分布のスケールが、もし\(n \rightarrow \infty\)まで続けたらどれだけ小さくなるかを表しているので、potential scale reductionという。
系列がうまく混合しているなら、任意のestimand \(\psi\)について近似的な「独立ドローの有効数」を計算できる。
もしドローが本当に独立だったら、\(B\)は事後分散\( var(\psi|y) \)の不偏推定である。実際には自己相関があるので、\(B\)は過大になっている。
\(E(\psi|y)\) の事後推定である\( \bar{\psi}_{..} \)に注目しよう(分布の要約の仕方はほかにもあるけれど)。その分散は漸近的に下式で与えられる。ラグ\(t\)の自己相関を\(\rho_t\)として、$$ lim_{n \rightarrow \infty} \ mn \ var(\bar{\psi}_{..}) = \left( 1+2\sum_{t=1}^{\infty} \rho_t \right) var(\psi|y) $$ ここから、$$ n_{eff} = \frac{mn}{1+2 \sum_{t=1}^{\infty} \rho_t } $$を有効標本サイズという。
漸近的な性質を使うのってどうなの、と思われるかもしれないけど、ここでは系列が十分に長く、目標分布へのおよその収束が起きているというのが前提だから、これでよいのである。
では、\( \rho_t \)をどう推定するか。
\( E(\psi_i – \psi_{i-t})^2 = 2(1-\rho_t) var(\psi) \)より、\( \rho_t = 1 – \frac{E(\psi_i – \psi_{i-t})^2}{2 var(\psi)} \) である。分母のほうはさきほど求めた \( \hat{var}^{+} \)を使い、分子のほうはバリオグラム $$ V_t = \frac{1}{m(n-t)} \sum_{j=1}^m \sum_{i=t+1}^n (\psi_{i,j} – \psi_{i-t, j})^2 $$ を使って $$ \hat{\rho}_t = 1 – \frac{V_t}{2 \hat{var}^{+}} $$ 残念ながら、\(t\)をあまり大きなところまで求めてしまうとノイズが大きくなってしまう。そこで、ラグ0からはじめて、隣り合う自己相関の和 \( \hat{\rho}_{2t} + \hat{\rho}_{2t+1} \) が最初に負になるまで足し続ける。こうして$$ \hat{n}_{eff} = \frac{mn}{1+2\sum_{t=1}^{T} \hat{\rho}_t } $$を求める。\(T\)は、\(\hat{\rho}_{T+1} + \hat{\rho}_{T+2} \)が最初に負になる奇数である。[←なぜ奇数でないといけないんだろう…?]
\( \hat{R} \) も \( \hat{n}_{eff} \)も平均と分散に基づいているので、事後分布が正規から離れているときにはうまくいかない。そういう場合は計算する前に変換するのが良い。すべての値が正である量なら対数をとるとか、すべての値が0-1に落ちる量ならロジットを採るとか、裾が長かったら順位をとるとか。
このように、関心あるスカラーそれぞれについて、PSRとESSを計算し、収束を監視するとよい。
関心あるすべてのスカラーのestimandについてPSRを求めることを勧める。すべての\( \hat{R} \)が1に近くなるまで(問題によるけれど、普通は1.1を閾値とする)、シミュレーションを続けるべし。
ESSはシミュレーションから得られる精度がどのくらいかを教えてくれる。\( \hat{n}_{eff} \)が\(5m\)になるまで(つまり系列あたり10個の独立ドローが得られたのと同じになるまで)シミュレーションを続けるのをお勧めする。
関心あるすべてのestimandについて、\(\hat{R}\)が1に近く、\(\hat{n}_{eff}\)がチェーンあたり10以上であっても、目標分布の重要な領域が初期分布によって捉えられておらず、シミュレーションのアルゴリズムでは容易にたどり着けないせいで、実はシミュレーションは収束にほど遠い、ということもありうる。
およその収束を宣言するということは、実際には、それぞれの系列は定常に見えますしお互いによく混合していますと結論することに過ぎない。このチェックは仮説検定ではない。p値も統計的有意性もない。あくまで実質的な有意性(ないしなんらかの慣習)の観点から、収束からの逸脱を評価しているに過ぎない。