Chib, S., Greenberg, E. (1995) Understanding the Metropolis-Hastings Algorithm. American Statistician, 327-335.
MCMCの古典的手法であるMetropolis-Hastingsアルゴリズムについての解説。Hitchcock(2003) によれば、この手法の普及させる立役者となった論文なのだそうだ。
現在では日本語でも優れた解説が山のようにあるのに、いまになってこれを読んでいるのは無駄なような気もするんだけど… なんとなく目を通してしまった。
1. イントロダクション
[略]
2. 受容-棄却サンプリング
[略]
3. MCMCシミュレーション
[この節のみ細かくメモする]
連続状態空間におけるマルコフ連鎖理論に対する通常のアプローチは、遷移カーネル\(P(x, A)\)が出発点となる。ただし\(x \in \mathcal{R}, A \in \mathcal{B}\), \(\mathcal{B}\)は\(\mathcal{R}^d\)上のボレル\(\sigma\)加法族である。[いきなりこういう難しい言い方すんの、ほんとにやめてほしい]
遷移カーネルは、\(x\)から集合\(A\)上のある点へと移動する確率を表現する条件つき確率関数である。\(P(x, \mathcal{R}^d) = 1\)である。\(P(x, \{x\})\)は必ずしも0でない。
マルコフ連鎖理論の主たる関心は、不変分布\(\pi^*\)が存在する条件、そして遷移カーネルの反復によってその不変分布に収束する条件、を決めることである。
不変分布は$$ \pi^*(dy) = \int_{\mathcal{R}^d} P(x, dy) \pi(x) dx $$ を満たす。ただし、\(\pi\)は\(\pi^*\)のルベーグ測度に関する密度であり、従って\(\pi^*(dy) = \pi(y)dy\)である。
[うぐぐぐ… よくわかんないけど、左辺の\(\pi^*(\cdot)\)は確率密度関数ではなくて確率関数で、だから引数として微小な幅を与えているんだけど、なんなら確率密度関数\(\pi(y)\)に微小な幅を掛けたものだと思ってくれてもいいよ、ということなのかな。右辺は遷移カーネル\(P(x, A)\)の\(x\)側を潰した周辺分布になっている。つまり不変分布は遷移カーネルを移動先の周辺分布に潰したものですよ、ということではないかと。あっているかな？？]
\(n\)番目の反復は、\(P^{(1)}(x, dy) = P(x, dy)\)として$$ P^{(n)}(x, A) = \int_{\mathcal{R}^d} P^{(n-1)} (x, dy) P(y, A) $$となる。以下に述べる条件の下で、\(n\)番目の反復が\(n \rightarrow \infty\)のときに不変分布に収束することを示せる。
MCMCはこの理論を裏側からみたものである。不変分布は既知で(目標密度\(\pi(\cdot)\))、遷移カーネルがわからない。大きな\(n\)について\(n\)番目の反復が\(\pi(\cdot)\)に収束するような遷移カーネル\(P(x, dy)\)を見つけたい。
適切な遷移カーネルを探すのは至難の業だが、以下のように考えると道が開ける。
遷移カーネルが以下のように書けるとしよう。$$ P(x, dy) = p(x, y)dy + r(x) \delta_x(dy)$$ \(p(x, y)\)はなんらかの関数で、\(p(x, x)=0\)である。\(\delta_x(dy)\)は、\(x \in dy\)のときに1, そうでないときに0となる関数である。\(r(x) = 1 – \int_{\mathcal{R}^d} p(x, y)dy \)は連鎖が\(x\)に留まる確率である。\(r(x) \neq 0\)となる確率は0でないので、\(p(x, y)\)の\(y\)を通じた積分は必ずしも1にならない。
[ここで訳が分からなくなってしまったんだけど… おさらいすると、左辺の\(P(x, dy)\)は移動元の点\(x\)と移動先の点\(y\)の同時確率密度関数ではなくて、点\(x\)がきまったら(任意の\(y\)の確率密度ではなくて)\(y\)が任意の区間に落ちる確率が決まるような関数である。右辺は、第一項が「動いたとき」、第二項が「動かなかったとき」を表していて、それぞれの頭にそうなる確率がかかっている(第一項に\(1-r(x)\)がかかっているとわかりやすいんだけど、\(p(x,y)\)にまとめてしまっているのだと思う)。さて、第一項の\(p(x, y)\)は\(x, y\)の同時確率密度関数もどきの関数で、\(p(x, y)dy\)が\(x\)から\(y\)周りの微小な区間へと移動する確率を表すんだけど、動いたときのことしか表さないので、\(p(x, x)\)は0にしてあるし、\((x, y)\)平面を水平に切って積分しても1より小さくなる。第2項は、動かないんだから、\(dy\)が\(x\)を含む区間ならば1、そうでなければ0であるような関数を用意すればよい。ということですよね]
さて、この式の\(p(x, y)\)が$$ \pi(x)p(x,y) = \pi(y)p(y,x)$$を満たすとき、\(\pi(\cdot)\)は\(P(x, \cdot)\)の不変密度になる。この条件のことを、可逆性reversibilityとか、詳細釣り合いdetailed balanceとか、微視的可逆性 microscopic reversibilityとか、時間可逆性time reversibilityなどという。
証明しよう。不変分布の右辺からスタートする。$$ \int P(x, A) \pi(x) dx = \int \left[ \int_A p(x,y) dy \right] \pi(x) dx + \int r(x) \delta_x(A) \pi(x) dx$$ 第1項だけみよう。さきに\(x\)で積分する。$$ \int \left[ \int_A p(x,y) dy \right] \pi(x) dx = \int_A \left[ \int p(x,y) \pi(x) dx \right] dy$$ 可逆性が成り立っていれば $$ = \int_A \left[ \int p(y,x) \pi(y) dx \right] dy $$ \(r(x) = 1 – \int p(x, y)dy \) より $$ = \int_A (1-r(y)) \pi(y) dy$$ 第2項をみよう。$$ \int r(x) \delta_x(A) \pi(x) dx = \int_A r(x) \pi(x) dx$$ 足すと $$ \int_A (1-r(y)) \pi(y) dy + \int_A r(x) \pi(x) dx = \int_A \pi(y) dy$$ ね？ つまり可逆性が満たされれば\(\pi^*(dy) = \pi(y)dy\)は不変なのだ。
[うううう… 納得できるような、そうでもないような… 連続状態空間を考えているせいで、話が初心者向けでなくなっているように思う。Haggstromの本の該当箇所を読み直したほうがよさそう]
4. Metropolis-Hastingsアルゴリズム
[M-Hの説明。メモ省略]
5. 実装上の諸問題: \(q(x, y)\)の選択
候補を生成する密度関数\(q(x,y)\)としてはいろんな路線がある。
- Metropolis(1953)はなんらかの多変量密度\(q_1(\cdot)\)を使って\(q(x,y) = q_1(y-x)\)を考えた。いいかえると、\(q_1\)に従うランダム変数を\(z\)として\(y = x + z\)とするわけである。ランダム・ウォーク連鎖といってもよい。\(q_1\)としては、MVNとか多変量\(t\)分布とかが考えられる。
- \(q(x, y) = q_2(y)\)という路線。現在位置\(x\)とは独立に決まるわけで、独立連鎖と呼んでよい。
- もし\(\pi(\cdot)\)の形式が既知なら、それを使うという手もある。たとえば、一様分布\(\psi(t)\), なんらかサンプリングできる分布を\(h(t)\)として\(\pi(t) = \psi(t)h(t)\)と書けるとしよう。このときは\(q(x, y) = h(y)\)とすればよい。
- 受容-棄却サンプリングを使い、擬似支配密度を使うという手もある。[\(\pi(x)\)を覆うような\(f(x)\)を用意して…という話ですね]
- 自己回帰連鎖。分布\(q\)に従う確率変数を\(z\)として\(y = a + B(x -z)+ z\)とする。\(q(x,y) = q(y – a -B(x-a))\)となる。\(B = -I\)とすれば負の自己相関を持つ連鎖になる。[へえええ？ そんな考え方があるのか]
肝心なのは、\(q(x,y)\)の広がりの決め方である。広すぎると変なところに飛んじゃって受容率が下がるし狭すぎるとサポートをカバーできなくなる。
ランダムウオーク連鎖\(q_1\)の場合、目標密度も提案密度も正規なら、受容率は1次元なら0.45くらい、6次元くらいなら0.25くらいがよいといわれている。0.5くらいがいいという話もある。
独立連鎖\(q_2\)の場合、提案分布の裾が目標分布より高いことが大事である。
6. M-Hアルゴリズムの応用
6.1 M-H受容棄却アルゴリズム
[パス]
6.2 Block-at-a-Timeアルゴリズム
たとえば、\(x = (x_1, x_2)\)というように分けられるとしよう(\(x_i \in \mathbb{R}^{d_i}\))。さらに、条件付き遷移カーネル\(P_1(x_1, dy_1|x_2)\)があって、\(x_2\)を固定したら\(\pi^*_{1|2}(dy_1 | x_2)\)がちゃんと不変分布になる、つまり$$ \pi^*_{1|2}(dy_1 | x_2) = \int P_1(x_1, dy_1|x_2) \pi_{1|2}(x_1|x_2)dx_1$$ となるし、\(P_2(x_2, dt_2 | x_1)\)についてもそうだとしよう。このとき、驚いたことに、遷移カーネルの積は不変密度\(\pi(x_1, x_2)\)を持つ。これをカーネル積の原理という。
この原理は実践的にみてものすごく重要である。それぞれのカーネルからドローしていけばいいことになるからだ。それに、同時密度に収束するようなカーネルを見つけるより、それぞれの条件付き密度に収束するようないくつかの条件付きカーネルを見つける方が容易である。
カーネル積の原理を証明するためには、\(x\)の要素を通じて「スキャン」するということの性質を特定する必要がある。遷移カーネル\(P_1(\cdot, \cdot|x_2)\)が\(x_1, x_2\)の下で\(y_1\)を生み、遷移カーネル\(P_2(\cdot, \cdot|y_1)\)が\(x_2, y_1\)の下で\(y_2\)を生むとしよう。その積を\(\pi^*(\cdot, \cdot)\)として、$$ \int \int P(x_1, dy_1|x_2) P(x, dy_2|y_1) \pi(x_1, x_2) dx_1 dx_2$$ $$ = \int P_2(x_2, dy_2|y_1) \left[ P_1(x_1, dy_1 | x_2) \pi_{1|2}(x_1 | x_2) dx_1 \right] \times \pi_2(x_2) dx_2 $$ カッコのなかは \( \pi^*_{1|2}(dy_1 | x_2)\)になっているので $$ = \int P_2(x_2, dy_2|y_1) \pi^*_{1|2}(dy_1 | x_2) \pi_2(x_2) dx_2$$ ベイズの定理より$$ = \int P_2(x_2, dy_2|y_1) \frac{\pi_{2|1}(x_2|y_1) \pi^*_1(dy_1)}{\pi_2(x_2)} \pi_2(x_2) dx_2$$ \(\pi^*_1(dy_1)\)を前に出して $$ = \pi^*_1(dy_1) \int P_2(x_2, dy_2|y_2) \pi_{2|1}(x_2|y_1) dx_2$$ 積分しているところは\( \pi^*_{2|1}(dy_2 | y_1)\) になっているので$$ = \pi^*_1(dy_1) \pi^*_{2|1}(dy_2 | y_1)$$ $$ =\pi^*(dy_1, dy_1)$$ というわけで、このとき\(\pi^*(dy_1, dy_2)\)は不変分布である。[ほ、ほんとだーーー]
この結果を念頭において、M-Hアルゴリズムのいくつかの重要な特殊ケースを紹介しておこう。
- いわゆるGibbsサンプラー。遷移カーネルを\(P_1(x_1, dy_1 | x_2) = \pi^*_{1|2}(dy_1 | x_2), P_2(x_1, dy_1|y_1) = \pi^*_{2|1}(dy_2 | y_1)\)とする。つまり、サンプルを完全条件付き分布から直接に生成するわけだ。
- “M-H within Gibbs”アルゴリズム。たとえば\(\pi_{1|2}(y_1|x_2)\)がintractableなとき、これはM-Hアルゴリズムからサンプリングし、他は完全条件付き分布から直接に生成する。
もともとHastingsは一回に一部だけをスキャンするという方法について議論していたわけで、だからGibbsサンプラーはM-Hの特殊ケースに過ぎないし、”M-H within Gibbs”なんていういい方はおかしい。
7. 事例
[ニ変数正規分布のシミュレーション、AR(2)時系列モデルの事後分布のシミュレーション。パス]
8. 結論
[略]