読了: Roy (2020) MCMC収束診断レビュー

Roy, V. (2020) Convergence Diagnostics for Markov Chain Monte Carlo. The Annual Review of Statistics and Its Application, 7(15), 1-26.

　仕事の都合でちょっと悩んだことがあって、現実逃避のために読んだ奴。
　題名通り、MCMCの収束診断についてのレビューである。正直、そんなに関心ある話題ではないのに、そしてすでにVats, et al.(2020), BDA3の該当部分を読んでいるのに、私にしては勉強熱心なことだ…

1. イントロダクション
　MCMCを実装するときに大事な問題がふたつある。

どこからはじめるか？つまり、マルコフ連鎖がここで定常性に収束したということをどうやって決めるのか？
いつやめるのか？つまり、モンテカルロ推定量が母集団の量に収束したとどうやって決めるのか？

　標準的な条件の下では、マルコフ連鎖$ \{X_n\} $の各点の分布は$n \rightarrow \infty$とともに定常分布$\pi$に収束する。しかし、初期分布が$\pi$から離れているときほど、$\pi$の近似には長い連鎖が必要になる。というわけで、初期の$n’$個の実現値を捨てるのがふつうである(burn-in)。ただしこの考え方には批判もある。 Geyer(2011, “Handbook of MCMC”)をみよ。
　関心があるのはたいてい、$\pi$のなんらかの関数$g$の平均 $E_\pi g \equiv \int_x g(x) \pi(x) dx $である。そのMCMC推定量は$\bar{g}_{n’,n} \equiv \frac{1}{n-n’}\sum^n_{n’+1} g(X_i)$である。$n$が足りないと推論が不正確になるおそれがあるし、長すぎるのも資源の無駄である。
　理論解析によって定常性への距離の分析的上界を求めることはできるだろう。またMCエラー$\bar{g}_{n’, n^*} – E_\pi g$の漸近分布に基づく標本サイズ計算で、誠実な停止値$n^*$を見つけることもできるだろう。そういう解析がない場面では、サンプラーと推定量の収束をチェックするツールが使われることが多い。もっとも、後述するように、収束を確実に決められるわけではない。
　90年代初頭以来、MCMCの収束診断ツールの開発に多大な努力が費やされてきた。本論文ではそれらをレビューする。

2. MCMC診断
2.1 誠実MCMC
　まずは、$n’, n^*$をrigorousにみつける方法について。
　$X_n$の密度を$f_n$とする。標準的な条件の下で、$ \frac{1}{2} \int_x | f_n(x) – \pi(x)| dx $は$n \rightarrow \infty$につれて0に近づく。つまり、$X_n$は全変動(TV)ノルムにおいて$\pi$に従う確率変数へと収束する。
　そこで、$ \frac{1}{2} \int_x | f_{n’}(x) – \pi(x) | dx$がたとえば0.01を下回るような$n’$をみつける。Jones & Hobert (2001)はこれをburn-inの誠実値と呼んだ。
　残念ながらTVノルムなんてふつうはわからない。TVノルムのboundを量的に構築するのは難しいんだけど、最近 drift and minorization テクニックというのがありましてですね。いろんなMCMCアルゴリズムを分析できるし、よくある要約指標のMCMC推定量の収束をチェックするのに使えます。さらに、TVノルムじゃなくて別の指標(たとえばWasserstein距離)を使うとか、$\log \pi$の導関数だけに依存する乖離指標を使う(つまり正規化定数がわからんような場合にも使える)とか、いろいろ進展があります。
　[よくわからん。最後の奴は、$\log \pi$の導関数がわかっている場面ってこと？それって$\pi$の分布形が既知だっていってんのと同じで、共役事前分布が用意できればもうMCMCなんていらないんじゃないかしらん…]

　目標密度のある特性を $E_\pi g$とする。マルコフ連鎖についての大数の強法則に基づき、$\{X_n\}_{n \geq 0}$がappropriately irreducibleならば[よくわからんけど、行けそうなところに実際に行けるならってことかな]、$ \bar{g}_{n’, n} \equiv \frac{1}{n-n’} \sum^n_{i=n’+1} g(X_i)$は$E_\pi g$の強一致推定量である。つまり、任意の固定された$n’$について、$n \rightarrow \infty$につれてほとんど確実に$\bar{g}_{n’, n} \rightarrow E_\pi g$である。
　停止ルールについて論じているときは、一般性を失うことなく$n’ =0$とできる。以下では$\bar{g}_{0, n}$を簡単に$\bar{g}_n$と書く。
　通常のモンテカルロ法と同じく、大数の強法則により、標本平均$\bar{g}_n$で$E_\pi g$を推定することが正当化できる。もし中心極限定理が成立するなら、$E_\pi g$の区間推定の幅に基づいて標本サイズを求めれば適切な$n^*$が決められるわけだ。実際、いくつかの条件のもとで$$ \sqrt{n} (\bar{g}_n – E_\pi g) \rightarrow N(0, \sigma^2_g) \ \ \mathrm{as} \ \ n \rightarrow \infty $$ が成り立つ。[最初の矢印の上に d が載っている。分布収束、つまり極限で左側の分布と右側の分布がぴったり一致するってことね]
　ここで $\sigma^2_g \equiv \mathrm{Var}_\pi (g (X_0)) + 2 \sum^{\infty}_{i=1} \mathrm{Cov}_\pi (g(X_0), g(X_i)) $である。マルコフ連鎖に自己相関があるので第2項がつく。
　もし$\sigma^2_g$の一致推定量$\hat{\sigma}_{g,n}$があれば、$\bar{g}_n$といっしょにSE$\frac{\hat{\sigma}_{g,n}}{\sqrt{n}}$も報告すればいい。ないし$E_\pi g$の$100(1-\alpha)$%信頼区間は標準正規分布の$(1-\alpha/2)$分位点を$z_{\alpha/2}$として$\bar{g}_n \pm z_{\alpha/2} \frac{\hat{\sigma}_{g,n}}{\sqrt{n}}$だから、それを報告してもいい。
　で、$100(1-\alpha)$%信頼区間の半分の幅が事前に決めた閾値$\epsilon$を下回ったらMCMCを止められばよい。Jones & Hobert (2001)はこれを誠実停止と呼んだ。

　実用的には、固定幅停止ルール(FWSR)というのが提案されている。最初は$\sigma^2_g$がわかんないから、とにかく$\tilde{n}$回は回す(ふつうは10000回くらい)。その後、適当な分位点を$t_*$として$$ t_* \frac{\hat{\sigma}_{g,n}}{\sqrt{n}} + \frac{1}{n} \leq \epsilon$$ が満たされ次第止める。[左辺の第1項はわかる、なんらかの信頼区間の幅だよね。第2項はなぜ必要なのだろうか？すごく小さな値になるから気にしなくていいんだろうけど、でも気になるなあ]
　誠実停止のためには、$\bar{g}_n$について中心極限定理が成り立つこと、そして一致推定量$\hat{\sigma}_{g,n}$があること、が必要になる。前者のためには、TVノルムがなんらかの速度でゼロに収束しないといけない。そのためにはマルコフ連鎖の幾何学的エルゴード性が必要で、その証明のためには…[以下、なんだかお経のようになってきたので中略…] というわけで、いろんなMCMCアルゴリズムについて中心極限定理の成立が証明されている。

　ここまでの話は単変量の平均の推定であった。ベクトルの場合にはどうするか。
　Vats, Flegal, & Jones (2019 Biometrika)はこう提案している。$$ \sqrt{n} (\bar{g}_n – E_\pi g) \rightarrow N(0, \Sigma_g) \ \ \mathrm{as} \ \ n \rightarrow \infty $$ が成り立つならば、$\Sigma_g$の一致推定量$\hat{\Sigma}_{g,n}$を使って、$E_\pi g$の$100(1-\alpha)$%漸近信頼領域$C_\alpha(n)$を構成できる。で、$\tilde{n}$回回した後で、$$ (\mathrm{Vol} \{ C_{\alpha}(n) \})^{\frac{1}{p}}
+ \frac{1}{n} \leq \varepsilon $$ が満たされたら止めればよい。
　単変量のときの$$ t_* \frac{\hat{\sigma}_{g,n}}{\sqrt{n}} + \frac{1}{n} \leq \epsilon$$と比べると、$\frac{1}{n}$の項を別にすれば、$\varepsilon = 2\epsilon$となっていることに注意。[しばらく固まったあとでようやく飲み込めた。$t_* \frac{\hat{\sigma}_{g,n}}{\sqrt{n}}$は信頼区間の幅の半分という意味なのだ。いっぽう$\mathrm{Vol}$は単変量でいえば信頼区間の幅そのものだもんね]

　なお、マルコフ連鎖の中心極限定理が理論的には証明されていない場合でも、バッチ平均と分散のスペクトル推定量[?] を使ったFWSRが実装されている。Rだとmcmcseパッケージがそれだ。[←へえー？ってことは、難しいことはよくわかんないけど、ときどきMCMCを中断してmcmcseパッケージでパラメータ推定量の信頼区間を求めて、それが十分短くなったら終了、っていうのがアリなわけね。そんなら俺にもできちゃうなあ]

2.2 相対的固定幅停止ルール
　上述のFWSRの固定幅を相対的に決めるという考え方もある。$$ t_* \frac{\hat{\sigma}_{g,n}}{\sqrt{n}} + \frac{1}{n} \leq \epsilon \bar{g}_n$$ が満たされたら止めるとか。母標準偏差 $\lambda_g$の強一致推定量を$\hat{\lambda}_{g,n}$として$$ t_* \frac{\hat{\sigma}_{g,n}}{\sqrt{n}} + \frac{1}{n} \leq \epsilon \hat{\lambda}_{g,n}$$ が満たされたら止めるとか(SDFWSRという)。
　
　多変量の場合はどうするか。まず、周辺チェーンについてSDFWSRを使うという手がある。いっぽうVats, et al. (2019)は交差相関を無視してはあかんよと批判し、こういうのを提案した。標本共分散行列を$\hat{\Lambda}_{g,n}$として、$$ (\mathrm{Vol} \{ C_{\alpha}(n) \})^{\frac{1}{p}}
+ \frac{1}{n} \leq \varepsilon (|\hat{\Lambda}_{g,n}|)^{\frac{1}{2p}}$$が満たされたら止める。

2.3 有効標本サイズ
　マルコフ連鎖標本と同じSEを持つ独立な標本の個数を有効標本サイズ(ESS)という。もっとも一般的な定義は$$ ESS = \frac{n}{1 + 2\sum_{i=1}^\infty \mathrm{Corr}_\pi (g(X_0), g(X_i))} $$というもの。別の定義として$$ ESS = n \frac{\lambda^2_g}{\sigma^2_g} $$ というのもある。まあとにかく、ESSの一致推定量が事前に指定した値を超えたら止めよう、という考え方もできる。

　多変量の場合、Vats et al.(2019)は多変量ESSというのを定義している。$g$が長さ$p$のベクトルを返し、その母共分散行列を$\Lambda_g$, 標本共分散行列を$\Sigma_g$として、$$ mESS = n \left( \frac{|\Lambda_g|}{|\Sigma_g|} \right)^{1/p} $$ どこで止めるかというと、$\chi^2_p$の$(1-\alpha)$分位点を$\chi-2_{1-\alpha,p}$とし、$100(1-\alpha)\%$漸近信頼領域のボリュームに求められる正確性レベルを$\epsilon$として、$$ \hat{mESS} = n \left( \frac{|\hat{\Lambda}_{g,n}|}{|\hat{\Sigma}_{g,n}|} \right)^{1/p} \geq \frac{2^{2/p} \pi}{(p \Gamma(p/2))^{2/p}} \frac{\chi^2_{1-\alpha, p}}{\varepsilon^2}$$を満たしたら止める。実はこれ、上述の多変量版SDFWRSとほぼ等価だということがわかっている。[おおおお。なんだかよくわからんが、これ、俺が欲しいものに近い…]
　単変量ESSはRではcodaパッケージとmcmcseパッケージに、多変量ESSはmcmcseパッケージに入っている。

2.4 Gelman-Rubin診断
　いまもっとも人気の方法はGR診断であろう。これは初期値の異なる複数のマルコフ連鎖があるときに使える方法である。連鎖が$m$本あるとして、$X (\sim \pi)$の分散の推定量として次のふたつを考える。
　ひとつめ、連鎖内推定量$W$。連鎖のなかの平均を$\bar{X}_{i.}$として、$$ W = \frac{\sum_{i=1}^m \sum_{j=0}^{n-1} (X_{ij} – \bar{x}_{i.})^2}{m(n-1)}$$ ふたつめ、プールした分散推定量$\hat{v}$。すべての連鎖を通した平均を$\bar{X}_{..}$として、連鎖間の分散は $\frac{1}{m-1} \sum_{i=1}^m (\bar{X}_{i.} – \bar{X}_{..})^2$なので、プールすると $$ \hat{V} = \frac{n-1}{n} W + \frac{1}{m-1} \sum_{i=1}^m (\bar{X}_{i.} – \bar{X}_{..})^2$$ となる。この2つの比$$ \hat{R} = \frac{\hat{V}}{W} $$ をpotential scale reduction factor (PSRF)と呼ぶ。初期値は分散が大きすぎるはずだから、最初$\hat{V}$は大きすぎ、いっぽう$W$は小さすぎるはずである。そこで、$\hat{R}$が1に近づくまで回しましょう、というわけである。閾値として1.1が良く使われている。
　なお、上述の$\hat{R}$の定義が広く用いられているが、実はこれ、オリジナルのGelman & Rubin(1992)による定義とはちょっぴり異なっている。オリジナルの定義では、最初に連鎖の前半を捨てることになっていた。さすがにそれは標本の無駄遣いなのでお勧めできない。[へー]
　Vats & Knudson (2018 arXiv) は、連鎖間の分散のところを、$\bar{X}_n$の漸近分散のバッチ平均推定量に差し替えることを提案している。するとESSと関連付けることが可能になるし、連鎖が一本しかなくても算出可能になる。
　[ああ、シングル・チェーンのMCMCでも$\hat{R}$を出力するソフトがあるのはそういうわけか。恥ずかしながらはじめて知った]

　多変量だったらどうするか。Brooks & Gelman (1998 J.Comp.Graph.Stat.)は多変量PSRFというのを提案している。プールした共分散行列を$\hat{V}^*$, 連鎖内共分散行列を$W^*$, 連鎖間共分散行列を$B^*$として$$ \hat{R}_p = \max_a \frac{a^\top \hat{V}^* a}{a^\top W^* a} = \frac{n-1}{n} + \left( 1+\frac{1}{m} \right) \lambda_1 $$を求める。ただし、$\lambda_1$とは行列 $ \frac{1}{n} (W^*)^{-1} B^* $の最大の固有値である。[へええええ？かんたんに計算できちゃうじゃん。なんだかキツネにつままれたような気分だ]

　GR診断は、codaをはじめとしてたくさんのRパッケージに入っている。[いや、$\hat{R}_p$を出すパッケージはなくない？俺が気づいてないだけ？]

2.5 2スペクトル密度ベースの手法
　ここでは、漸近的分散推定値に基づく診断ツールをふたつ紹介しよう。

　まずはGewekeの手法。これは、時系列にはスペクトル密度というものが存在する、という仮定に基づいている。[周波数領域の話ってほんとに苦手なので、ここからしばらくは虚心に写経…]
　$E_\pi g$の推定において、$\bar{g}_n$の漸近分散は、$ \{ g(X_n), n \geq 0\}$を時系列として見たときの、時点ゼロにおけるスペクトル密度$S_g(0)$である。マルコフ連鎖を$n$回回し、最初の$n_A$個の平均を$\bar{g}_{n_A}$, 最後の$n_B$個の平均を$\bar{g}_{n_B}$としよう。Gewekeの統計量は $$ Z_n = \frac{\bar{g}_{n_A} – \bar{g}_{n_B}}{\sqrt{ \frac{\hat{S}_g (0)}{n_A} + \frac{\hat{S}_g (0)}{n_B}}} $$ である。つまり、平均の差を、$S_g(0)$のノンパラメトリック推定値を使って求めた標準誤差で基準化したものである。[学力不足で$\hat{S}_g (0)$の求め方がわかんないんだけど、うん、まあ主旨はわかるからいいや]
　Gewekeさんは最初の1割と後半(5割)を使うことを勧めている。
　このZ得点は2つの部分の独立性を仮定して求めている。つまり、Gewekeの収束診断は平均の等値性のZ検定で、標本の自己相関は、SEを求めるときに考慮しているわけだ。

　いっぽう、Heidelberger & Welch (1983)はスペクトル密度推定に基づきこういう方法を提案している。[数式が出てくるんだけど、まっっっったく理解できない。なにやってんのこれ。潔くあきらめます]

　どちらの方法もcodaパッケージにはいっている。
　多変量の場合についてはCowles & Carlin (1996 JASA) の提案がある。$g$を-2log(目標密度)としてGewekeの方法を使うというものである。[わからん。なにいってんの??? パラメータ推定値のチェーンの定常性じゃなくて、モデルの対数尤度のチェーンの定常性を診断するということだろうか。それマルコフ連鎖じゃなくないですか]
　まあとにかくどちらの方法も、ESSやGRと同様に、マルコフ連鎖について中心極限定理が成り立つことを前提としていることに注意してください。

2.6 Raftery-Lewis診断
　やりたいことが$g(X)$の分位点の推定だとしよう。$ W_n \equiv I (g(X_n) \leq u) $ という二値過程について考える。これ自体はマルコフ連鎖ではないけれど、Raftery & Lewisいわく、十分に大きな$k$をとれば、$ \{W_{nk} \}_{n \geq 0} $は近似的にマルコフ連鎖だ。sその遷移確率 $P(W_{nk} = j | W_{n-1, k} = i)$は単純に$$ \frac{ \sum_{l=1}^{n} I(W_{lk} = j, W_{l-1,k} = i) }{\sum_{l=1}^n I(W_{lk} = i)}$$ で推定できる。モデル選択で$k$を選び、2状態の遷移行列を推定すれば、burn-inの推定値もわかるし、中心極限定理を仮定すれば、どこで止めたらいいかもわかる。[…中略…]
　この診断もcodaパッケージに載ってます。
　[正直、この節は話に全然ついていけなかった。$k$ってのは適当なthinningパラメータってこと？ thinningしたマルコフ連鎖標本が適当に決めた閾値$u$を上回るかどうかの2状態マルコフ連鎖を考え、その遷移行列を推定する、ってこと？それがわかったら、その定常状態の推定値がわかり、さらにその信頼区間もわかり、その幅が求められた幅より小さくなるようなイテレーション数がわかる、ってこと？じゃあburn-inのサイズがわかるのはなんで？
　どうやらこの論文の簡潔な説明では理解できそうにないので、別の機会を待とう。というわけで、かなり中略してしまった。まああれだ、codaパッケージでRaftery-Lewis診断をやろうとして「おまえの連鎖は短すぎて無理だ帰れ」といわれたことがあるんだけど、その理由がなんとなくわかったので良しとしよう。なんか大変そうだもんね、これ]

2.7 カーネル密度ベースの手法
　ここでは、2本の連鎖(単一の連鎖ならその2つの部分)についてそれぞれカーネル密度推定して、その間の距離が0に近かったら収束と判断する、という方法を紹介する。GR診断の場合は連鎖のなんらかの積率を比較しているのに対して、この方法では分布全体を比べる。
　いろいろあるけど、Dixit & Roy (2017 J.Stat.Comp.Sim.)の提案を紹介しよう。

　2本のマルコフ連鎖$ \{X_{1j}\}, \{X_{2j}\} $があるとする。それぞれについてのカーネル密度推定を$p_{1n}, p_{2n}$とする。KLダイバージェンスを$$ KL(p_{in}|p_{jn}) = \int_x p_{in}(x) \log \frac{p_{in}(x)}{p_{jn}(x)} dx $$ として、対称的KLダイバージェンス $$ \frac{1}{2} (KL(p_{1n}|p_{2n}) +KL(p_{2n}|p_{1n}) )$$ を求める。これを使って検定する[小さかったら収束ってことね]。タイプIエラーは、マルコフ連鎖が収束してないときに「収束した」と判断しちゃうこととなる。連鎖が3本以上あるときは、すべてのペアについてのKLダイバージェンスを出す。これをツール1と呼ぼう。
　多変量の場合はどうするか。同時分布のカーネル密度推定はあてにならないので、それぞれの変数の周辺分布について調べ、検定の際の有意水準をボンフェローニ調整するという手がある。[うわー、検定力落ちそう… それって過度に楽観的な収束診断になりませんかね…]

　悲しいシナリオとして、目標分布が多峰で、すべての連鎖が同じ峰にスタックしちゃってて多峰であることに気づかない、という場合が考えられる。そこでツール2というのが提案されている。MCMC標本のカーネル密度推定と、目標密度とのKLダイバージェンスを求めるのである。$\pi(x) = f(x)/c$で正規化定数$c$が未知だとしますね。連鎖のカーネル密度推定と$\pi$のKLダイバージェンスに基づいて求めた$c$の推定値を$\hat{c}$とし、数値積分で求めた$c$の推定値を$c^*$として $$ T^*_2 = \frac{|\hat{c} – c^*|}{c^*} $$を求める。これは目標分布のうちマルコフ連鎖で捉えきれていない割合を表す。これが0.05を下回ったら止める、というアイデアである。これは単変量のときのみ使える。[ええええ？これって$f(x)$は既知ってことですよね？？？ううう、適用場面が想像できないよう]

　Dixit & Roy はKLダイバージョンスの視覚化を推奨している。たとえば連鎖が4本なら3×3の三角行列を描き、各セルに色を塗る。KLダイバージェンスが閾値を超えたら濃い色、でなかったら薄い色、というふうに。[なるほど、すべてのタイルが薄くなるのを確認しろってわけね]

2.8 グラフィカルな手法
　もっとも一般的なのはトレース・プロット。横軸がイテレーション、縦軸がマルコフ連鎖の実現値である。[…見方について。中略…]
　自己相関プロットも大事。横軸がラグ、縦軸がACFである。マルコフ連鎖の長さは、(ACFがゼロに近づくまでのラグ)x(大きな値)だけ必要である。
　running meanプロットというのも使われている。横軸がイテレーション、縦軸がその時点でのモンテカルロ推定値である。マルコフ連鎖が定常性に収束したかどうかはわからないが、いつ止めるかを判断する助けになる。もっともこれには批判もある。Geyer(2011 in “Handbook of MCMC”)をみよ。

　多変量の場合、各変数の周辺分布についてトレース・自己相関・running meanをみると同時に、変数間の相関も視覚化するとよい。

3. 事例
3.1 目標分布が指数分布のとき
　$\pi(x) = \exp(-x) \ (x > 0)$として、提案密度$q(x) = \theta \exp(-\theta x)$のメトロポリス・サンプラーを考えよう。$X_0 = 0.1$として、$\theta = 0.5$のチェーンと$\theta = 5$のチェーンをそれぞれ1本走らせた。1000イテレーション。
　トレースプロットと自己相関プロットをみてみよう。$\theta = 0.5$のチェーンでは、どうやらburn-inはいらなさそうで、ACFもラグ4あたりから無視できる。いっぽう$\theta = 5$のチェーンでは、なんだかburn-inが必要な感じがするし、ACFはラグ50でもまだ高めである。[←ここまで、ものすごく直観的な話をしていて、心が和む]
　では、codaパッケージで、$g(x) = x$として診断してみよう。GewekeのZ (とるのは頭1割と尻5割)では、平均の等値性はどちらも棄却。Heidelberger & Welch の定常性検定はどちらも通過。
　burn-inについて。Raftery-Lewisは38415イテレーションまで伸ばさないと計算できず、伸ばしたうえでやってみると、$\theta = 5$のほうが必要なburn-inが長い。誠実burn-inは…[めんどくさそうな話になってきたので中略]
　停止について。[めんどくさくなってきたのでメモは省略するけど、平均の真値1に対して、チェーンのrunning meanは$\theta = 0.5$ならさっさと近づくけど、$\theta = 5$だと全然ダメで、10万イテレーションあたりで0.8あたりから動けなくなってしまう。つまり、真の分布が平均1の指数分布のとき、提案分布の平均が2なら辿りつけるけど、0.2ならどんだけ回しても辿りつかないってことだよね、これ。そりゃそういうもんなんだろうけどさ、実例をみちゃうとびびるな…]

3.2 目標分布が6峰のとき
[力尽きたのでパス]

3.3 ベイジアン・ロジスティック・モデル
[力尽きたのでパスするけど、面白そうね。MCMCの診断で悩んだときにお手本にするといいかもしれない]

4. 結論と考察
　事例でご覧いただいたように、標本平均の実証的な収束診断は往々にして失敗し、真値からかけ離れたところでシミュレーションを止めてしまう。固定幅停止ルールかESS停止ルールのほうがよい。中心極限定理の成立が示せてないといけないんでしょう？その証明って大変なんでしょ？という心配はごもっともだが、幸い多くのモデルで証明されている。
　どの事例でもthinningしなかった。thinningは結局標本の無駄遣いなので、マルコフ連鎖のサンプリングよりも関心ある関数$g$の推定のほうが大変だ、とか、メモリだかディスクだかが足りない、といった場合を別にして、お勧めできない。
　目標分布が多峰の場合、並列チェーンに与えた初期値が高密度領域でなかったり、チェーンが短すぎて峰の間で動けなかったりすると、収束の検出はできない。よって、最終的には単一の長ーいチェーンが望ましい。いっぽう、マルコフチェーンのカップリング[?]を用いて並列MCMCをするという提案もある。

　実務家のみなさん、実証的な収束診断ツールに頼るときには気を付けて。多峰性が疑われるときにはとくに。
　過去20年間、よく使われる統計モデルの、いろんなMCMCアルゴリズムに関し、誠実モンテカルロ標本サイズを求める研究が進められてきた。しかし、MCMCアルゴリズムの理論解析はいまだ進展中の研究領域である。今後の魅力的な研究分野としては、Dixit & RoyのKLダイバージェンス・ベースの統計量の収束を理論的に検証することとか、超・高次元な場面でのMCMC収束診断とかが挙げられる。云々。
—————–
　やれやれ…
　３つの事例のうち最初の奴しか読んでないけど、Gelman-Rubinの収束診断が全然うまくいってないのでがっくりした。$\hat{R}$が1.1を切ったんで大丈夫っす！なんて、なかなか胸を張って言えないものなのね。参っちゃうなあ。モデルがブラックボックスなのはよくないけど、せめてパラメータ推定くらいは、ユーザからみてブラックボックスにしてほしい。統計学者のみなさん、そこんとこ、お願いしますよ… (愚痴)

読書日記

読んだ本を淡々と記録します

読了: Roy (2020) MCMC収束診断レビュー