Vats, D., Robertson, N., Flegal, J.M., Jones, G. (2020) Analyzing Markov Chain Monte Carlo output. WIREs Computational Statisitics. e1501.
都合で読んだ奴。MCMCの解説ではなくて、MCMCで得たサンプルをどうやって分析するか、という解説論文。
1. イントロダクション
サポート\( \mathcal{X} \supseteq \mathbb{R}^d (d \geq 1) \) を持つ確率分布 \(F\)について、その未知の特性を使って母集団について推論したい。たとえば、\( h:\mathcal{X} \rightarrow \mathbb{R} \)について、
$$ \mu_h = E_F[h(X)] = \int_\mathcal{X} h(x) F(dx) $$ に関心がある、とか。[ああ、数学得意な人ってなんでこういう難しい書き方をするんだろう… タンスの角に足の小指をぶつけますように…]
なんでもいいから、我々が求めている\( F \)の特性を\(p\)次元ベクトル\(\theta\)にまとめておく。
さて、多くの場合\( F \)は分析的に解けない。そこでモンテカルロ・サンプリングの登場である。
以下では、マルコフ連鎖\(\{X_1, X_2, \ldots, X_n\} \)を十分に大きな\( n \)だけ生成すれば \( X_n \dot{\sim} F \)となるようなMCMCについてのみ考える( \( \dot{\sim} \)は「漸近的に分布する」を表す)。つまり、\( n \)が十分に大きければ\( \hat{\theta}_n = \hat{\theta}(X_1, \ldots, X_n) \approx \theta \) となる。
本論文ではMCMC法そのものについては扱わず、MCMCサンプルを生成するなんらかの効率的な方法がもう手元にあるという場面について考える。
ふつうのMCMC実装では、初期状態\(X_1\)はなにか別の分布から抽出してくる。また、MCMCサンプルには系列相関がある。つまりマルコフ連鎖は独立でも同分布でもない。だから、MCMCサンプルの分析では次の2つが課題になる。
- 課題1. どこまで \(n\)を大きくすれば \( X_n \dot{\sim} F \) か。2章で扱う。
- 課題2. どこまで \(n\)を大きくすれば\( \hat{\theta}_n \approx \theta \) か。3章以降で扱う。
2. マルコフ連鎖理論と初期値
時間等質なマルコフ連鎖のダイナミクスは次のマルコフ遷移カーネルで記述できる。[またこういう難しい書き方を… 担々麺の汁がシャツに跳ねますように…]$$ P^n(x, A) := Pr(X_{n+j} \in A | X_j = x), n, j\geq 1$$ ただし \( P^1 \equiv P \)。\( X_1 \sim \nu \)ならば\(X_n\)の周辺分布は \( \nu P^n \) である。
MCMCでは目標分布 \( F \)が定常分布になるように、つまり\(FP^n = F\)となるようにマルコフ連鎖を生成する。だから\(X_1 \sim F\)としておけばそれぞれの\(X_n\)も\(F\)に従う… というのはもちろん夢物語ですけど、いま全分散距離を\( || \cdot ||\)と書くとして、標準的な条件のもとで、\(n \rightarrow \infty\)のときに\( || \nu P^n – F || \rightarrow 0\)となることが示せる(つまり\( X_n \)は\(F\)に分布収束する)。
さて、我々が突き止めたいのは、\(n \leq n^*\)ならば \(\epsilon > 0\)について \( || \nu P^n – F || < \epsilon \)となるような \(n^*(<\infty)\)である。それがわかれば胸を張って \( X_n \dot{\sim} F \) だといえるじゃないですか。
実は、こういう\(n^*\)を見つける方法もあるんだけど、多くの場合保守的すぎて実用性が乏しい[←へええ。やたら長く生成せよといわれちゃうってことでしょうね]。だから実務家は、トレースプロットを眺めたり、いわゆる収束診断をしたりするわけである。
厳密に言えば、有限の\(n\)で収束は起きない。収束してないという証拠がないからといって収束している証拠にはならない。
なお、多くの人が誤解してるんだけど、みなさんご存じGelman-Rubinの収束診断を含め、ほとんどすべての収束診断が取り組んでいるのは、実は課題1じゃなくて課題2である。
収束の速さは初期値によって変わる。極端にいえば、もし\(X_1 \sim F\)ならすべての\(n\)について\( ||FP^n -F|| = 0\)であり、マルコフ連鎖は(系列相関はあるけれど)\(F\)からの正確なドローである(夢物語ですけど)。というわけで、初期値は大事である。
定常な初期値からスタートするための方法として、perfect simulation, Bernoulli factories, (目標分布が低次元なら) 単純なaccept-reject sampler, などがある。[←ほひゃー。なにそれ全然わかんない]
最適化手法でもって確率が高い領域を見つけて初期値を得るという手もある。閉形式で書けるならMLEを使うのが早いし、他の手法もある(大域最適解を見つけるのが難しいけど、ともあれ確率の高い領域の値を初期値にするのは良いことである)。事前分布が注意深く選ばれているならばそこからドローしても良い。
Gibbsサンプラーのような方法だと、最初に更新するパラメータには初期値がいらないので、初期値が一番あてにならないパラメータから更新をはじめるという手もある。
3. 推定とサンプリング分布
我々はマルコフ連鎖の実現値 \( \{X_1, X_2, \ldots, X_n\} \)から、\(n \rightarrow \infty\)のときほぼ確実に\(\hat{\theta}_n \rightarrow \theta \)となるような推定量\(\hat{\theta}_n = \hat{\theta}(X_1, \ldots, X_n)\)をつくれるんだけど、いかに\(n\)が大きかろうが、未知のモンテカルロ誤差 \(\hat{\theta}_n – \theta\)は残る。ある種の中心極限定理に基づき、このモンテカルロ誤差の近似的な標本分布を求めることができる。以下の参考文献を参照なさい[メモ省略]。
いま関数 \( h:\mathcal{X} \rightarrow \mathbb{R} \) があり、\( E_F h(X) = \mu_h \)に関心があるとしよう。たとえば、\(h\)がidentity関数なら平均ベクトルに関心があることになる。
その推定は簡単で、$$ \bar{\mu}_n := \frac{1}{n} \sum_{t=1}^n h(X_t) $$ でよい。これは\(n \rightarrow \infty\)のとき\(\mu_h\)に概収束する。また中心極限定理の下で \( \sqrt{n} (\bar{\mu}_n – \mu_h) \)は\(N_p(0, \Sigma)\)に分布収束する。ここで$$\Sigma = \sum_{k=-\infty}^\infty Cov_F( h(X_1), h(X_{1+k}) ) $$ つまり\(\Sigma\)はは目標分布における\(h\)の共分散構造とマルコフ連鎖のせいで生じる系列ラグ共分散を表している。\(Cov_F\)と書いたのは、\(X_1 \sim F\)という仮定の下での期待値だという意味である。中心極限定理の成立のためには\( X_1 \sim F \)が必要だといっているわけではないことに注意。
[ここで分位点の推定と期待値関数の推定の話が1p弱。パス]
[なんだかそこはかとなく面倒くさくなってきたが、いやいや、論文の本題はここからであろう]
4. モンテカルロ誤差の推定
前節で述べたように、モンテカルロ誤差の近似的な標本分布を求めることができる。これを使ってシミュレーションの信頼性を評価できる。つまり課題2である。
\(\mu_h\)の信頼区間について考えよう。\(\Sigma_n\)を中心極限定理の下での\(\Sigma\)の推定量とする。次元\(p\), 自由度\(q\)のホテリングのT二乗分布の\(1-\alpha\)分位点を\(T^2_{1-\alpha, p, q}\)とする。漸近的な\(100(1-\alpha)%\)信頼区間は $$ C_\alpha(n) = \{ n(\bar{\mu}_n – \mu_h)^T \Sigma_n^{-1} (\bar{\mu}_n – \mu_h) < T^2_{1-\alpha, p, q} \} $$ 問題は \(\Sigma_n\)をどうするかである。大きく分けて、(1)スペクトラル分散推定量, (2)バッチ平均推定量, (3)初期系列推定量、がある。以下ではもっとも計算しやすいバッチ平均推定量について述べる。
多変量バッチ平均推定量は、マルコフ連鎖の重複しないバッチについて考え、それぞれのバッチの標本平均から、標本共分散行列を構築する。\(a\)をバッチ数、\(b\)をバッチサイズとしよう(\(n=ab\)である)。\(k=0, \ldots, a-1\)について$$ \hat{Y}_k := b^{-1} \sum_{t=1}^b h(X_{kb+t})$$と定義する。\(\Sigma\)のバッチ平均推定量は$$ \hat{\Sigma}_{n,b} := \frac{b}{a-1} \sum_{k=0}^{a-1} (\bar{Y}_k – \bar{\mu}_n) (\bar{Y}_k – \bar{\mu}_n)^T $$ である。[はいはい、バッチごとの平均の共分散を調べる訳ね]
バッチ平均推定量の漸近特性については研究が進んでいるけれど、相関が高いときの小標本のパフォーマンスは疑わしい。最近では有限標本でのパフォーマンスを改善する提案があって[…メモ省略]
バッチ平均推定量のパフォーマンスはバッチサイズの選択で決まる。最適なバッチサイズの決め方として[…メモ省略]
5. シミュレーションの停止
\(n\)の決め方にはいくつかのアプローチがある。もっとも単純なのは、とりあえずなにかに決めといて、モンテカルロ誤差を推定し、それが大きすぎたらチェーンを伸ばす、という方法である。
再び期待値ベクトルの推定に焦点を当てよう。\(\Lambda = Var_F(h(X_1))\)とし、その行列式を\(|\Lambda|\)とする。$$ ESS := n \left( \frac{|\Lambda|}{|\Sigma|} \right)^{1/p} $$ を有効標本サイズという。有効標本サイズが十分に大きいということは、目標分布の変動にくらべてモンテカルロ誤差が十分に小さいということである。というわけで、$$ \hat{ESS} := n \left( \frac{|\Lambda_n|}{|\Sigma_n|} \right)^{1/p} $$を求めればよい。なお、ESSは \(h\) に依存する、つまり関心ある量がなにかに依存するという点に注意。
ESSはどのくらい大きければよいのか。いま、\(\mu_h\)の\(100(1-\alpha)%\)信頼区間を構成して、そのvolumeが\(F\)の下での\(h\)のvariabilityの\(\epsilon\)番目のfractionであるようにしたい[←このくだり、全然理解できない…]。この場合、シミュレーションは次を満たしたときに停止できる。$$ \hat{ESS} \geq \frac{2^{2/p} \pi}{(p \Gamma(p/2)^{2/p}} \frac{\chi^2_{1-\alpha, p}}{\epsilon^2} := M_{\alpha, \epsilon, p}$$ このルールに従って停止すれば、\(\epsilon \rightarrow 0\)のとき、信頼区間がカバーする確率は\(1-\alpha\)に収束する。実務的には、なんらかの最小の\(n^{*}\)ステップに達するまで停止しないほうがよい(\(\Sigma, \Lambda \)の初期の推定はプアだから)。\(n^{*}\)は\(M_{\alpha, \epsilon, p}\)にするとよい。これは事前に計算できる。
よく、平均の信頼区間の幅をどうしたいか決めて、そのために必要な標本サイズを求めることがあるけれど、上の式も考え方は同じである。信頼区間の幅に相当するのが\(\epsilon\)である。典型的には0.05以下とするのがよい。
課題2の意味での収束評価にはGelman-Rubin-Brooksの診断がよく用いられている[これってPSRのことだよね…]。オリジナルでは、\(m\)本の独立なチェーンから\(\Sigma\)を推定する。あいにく、たいてい\(m\)は小さいので、この統計量は変動が大きいんだけど、前の節で述べた\(\Sigma\)を使えば改善できる。さらに、この改善された統計量を\( \hat{R}^p \)として、$$ \hat{R}^p \approx \sqrt{ 1+\frac{1}{\hat{ESS}} } $$ であることが示されている。これは一本のチェーンから求められるし、複数のチェーンから求める方法も提案されている。
\( \hat{R}^p \)がどのくらい小さくなったら停止するか。$$ \hat{R}^p \leq \sqrt{1+\frac{1}{M_{\alpha, \epsilon, p}}}$$で停止すればよいといわれている。しかしESSを使ったほうがよいだろう。解釈しやすいし、変動が小さいし、計算が楽だから。
\( \hat{R}^p \)も\( \hat{ESS} \)も、推定を目指しているのであって(つまり課題2)、標本の代表性を調べているのではない(課題1ではない)、という点に注意すること。
6. 拡張
thinningとは[…中略…]。推定量の分散は大きくなるが、\(h\)の評価に時間がかかるときとか、問題が高次元すぎてMCMC標本を全部メモリに置いておけないときには有用である。
複数の並列チェーンもよく使われる。でも複数チェーンだからといって短くてよいということにはならないので注意。
発展的な話題としては… 高次積率の推定 (精度がすごく下がる)。adaptive MCMC。trans-dimentional MCMCによる不確実性の定量化。
7. ワークフロー
最後に、MCMCの出力を視覚化しながら分析するワークフローを示そう。
液晶プロジェクターの信頼性をベイジアン信頼性モデルで評価するという事例を取り上げる。メーカーいわく、ランプの寿命は1500時間である。そこで同じランプを31種類のプロジェクターに取り付けて故障までの時間を測った。プロジェクタ \(i\)での故障時間を\(t_i\)とする。
モデル。\(t_i\)は $$ T_i \sim Weibull(\lambda, \beta)$$ の実現値とする。平均故障時間(MTTF)と信頼性関数\(R(t)\)は、モデルより、$$ MTTF = \lambda^{-1/\beta} \Gamma \left( 1+\frac{1}{\beta} \right) $$ $$ R(t) = \exp(-\lambda t^\beta)$$ となる。
事前分布を \(\lambda \sim Gamma(2,5, 2350), \beta \sim Gamma(1,1)\)とする。事後分布は[…面倒くさいのでメモ省略。以下、これこれのサンプラーとこれこれの初期値を使ったという話。えーい全部省略だ]
さて、ここで私たちは、$$ h(\lambda, \beta) = \left( \lambda^{-1/\beta} \Gamma \left( 1+\frac{1}{\beta} \right), \exp(-\lambda t^\beta) \right)^T $$ なる関数の事後平均に関心を持っているわけである。
最小ESSを求めてみると[…中略…] 7529である。そこで、まず7529ステップ走らせて、\(\lambda\)と\(\beta\)のトレースプロットを観察する。注目点はちゃんとミキシングしているかである。だいじょうぶだった。
ESSを推定すると\(\hat{ESS} = 1196\)、7529より随分小さいので、一気に100000ステップまで走らせてみる。\(\hat{ESS} = 11834\)、7529を越えたのでこれでよしとする。
次に、MTTFとR(1500)のACFプロットとCCFプロットをチェックする。さらに、MTTFとR(1500)を軸にとって信頼領域を描く。ここで特に大事なのはラグ0のCCFである(ラグつきのCCFは単変量分析なら無視して良い)。
最後に、MTTFとR(1500)の周辺密度をプロットする。結果、MTTFの95%信用区間は(434, 834)、メーカーの主張よりはるかに低い。
云々。
… 正直むずかしくて、途中で写経になっちゃった箇所もあるんだけど、自分なりに勉強になったのでよしとしよう。要はESSと\( \hat{R}^p \)についてちゃんと勉強すればいいのね。そしてそれを勉強したからといって、それはドローしたモンテカルロ標本から推定した分布特性が真の分布特性を近似しているかどうかという話であって、ドローが分布収束しているかどうかという話とはまた別なわけね。
バーンインの話が全然出てこなかったのはちょっと意外であった。そういうもんなんすか？
最後の節のチャートの描き方も勉強になった。いつまでもトレースプロットをぼけーっと眺めてないで、パラメータ間の交差相関をちゃんとみろ、という話であろう。へへー(平伏)。