Gayer, C.J. (2011) Introduction to Markov Chain Monte Carlo. in Brooks, S., et al. (eds) “Handbook of Markov Chain Monte Carlo“, CRC Press.
MCMC解説のメモ、第二弾。2回で終わるかと思ったのだが、全然終わんないぞ。今回は12-16節。
12. MCMCの基本理論
我々は指定した分布を保存してくれるような更新メカニズムに関心がある。そういうメカニズムについて考えよう。
12.1 Metropolis-Hastings更新
ほしい定常分布の、基準化されていない密度を\(h\)とする。Metropolis-Hastings更新は以下である。
- 現在状態が\(x\)のとき、\(x\)のもとでの条件付き確率密度\(q(x, \cdot)\)から提案状態\(y\)をドローする。[前から疑問に思っているのだが、ここ、なぜ\( q(\cdot|x) \)と書かないのだろうか？]
- Hastings比を求める: $$ r(x, y) = \frac{h(y)q(y,x)}{h(x)q(x, y)}$$
- 次の確率で提案状態\(y\)への移動を受理する(Metropolis棄却): $$ a(x, y) = \min(1, r(x, y)) $$
Metropolis棄却と、OMCでいう棄却サンプリングの棄却をごっちゃにしないように。棄却サンプリングの棄却は、受理が起きるまで繰り返される(つまり状態は常に新しい)。Metropolis棄却は、状態は\(y\)になるかもしれないし、\(x\)のままかもしれない。
Hasting比の式をみるとわかるように、\(y\)があり得ないような値であっても大丈夫。
[Rコード例。略]
12.2 Metropolis-Hastings定理
では、M-H更新が\(h\)に関して可逆であることを証明しよう。もし棄却されたら、\(X_n = X_{n+1}\) だから、当然可逆である。もし受理されたらどうか。
[ここから逐語訳]
我々が示す必要があるのは次の式である: $$ E[f(X_n, Y_n) a(X_n, Y_n)] = E[f(Y_n, X_n) a(X_n, Y_n)] $$ ただし、\(f\)は(\(X_n\)が目指している定常分布を持つとして)期待値を持つ任意の関数である。つまり、我々は$$ \int \int f(x, y)h(x)a(x,y)q(x,y) dx dy \ \ \ \ \mathrm{(1.22)}$$の\(f\)の引数が交換可能であることを示す必要がある(状態が離散ならインテグラルを総和で置き換えること)。
このことは、もし$$ h(x)a(x,y)q(x,y) $$ の\(x\)と\(y\)が交換可能ならば成り立つ。なぜなら、\(x, y\)がダミー変数ならば(1.22)式の\(x\), \(y\)が交換可能になるからである。
あきらかに、\(x\)と\(y\)のうち積分ないし総和に貢献するのは\(h(x) \gt 0, q(x, y) \gt 0, a(x,y) \gt 0\)の集合のみである。さらに、不等式であるから\(h(y) \gt 0, q(x, y) \gt 0\)である。こうして、これらの不等式が成り立つ\(x, y\)については次の式が成り立つ:$$ r(y, x) = \frac{1}{r(x,y)} $$ \( r(x,y) \leq 1\)のとき、\(r(x, y) = a(x, y), a(y, x) = 1\)であるから、$$ h(x) a(x,y) q(x, y) = h(x) r(x, y) q(x, y) = h(y) q(y, x) = h(y) q(y, x) a(y,x) $$ \( r(x,y) \geq 1\)のとき、\(a(x, y) = 1, a(y, x) = r(y, x)\)であるから、$$ h(x) a(x,y) q(x, y) = h(x) q(x, y) = h(y) r(y, x) q(y, x) = h(y) q(y, x) a(y,x) $$ いずれの場合も\(x, y\)が交換可能であるので、証明できた。
[思いあまって逐語訳してみたが、全然！わからないーーー！
\(X_n\)と\(Y_n\)の同時分布上で\(X_n\)と\(Y_n\)をいれかえても同時密度が変わらないということが示せればいいんですよね？ 同時密度は、定常分布上の密度と提案密度と受理密度の積ですよね？ だから\(h(x)a(x,y)q(x,y)\)の\(x,y\)が交換可能であるということが示せればいいんじゃないの？ なんでこんな難しい話になるのか、全然理解できない！
これはもちろん私がアホなせいであって、状態が不可算連続である場合を視野に入れると、きっとこういうややこしい前置きが必要になってしまうのであろう。なぜ必要なのか全然わかんないけど。悲しいなあ。家庭教師とかを雇った方がよいのかもしれない…
ここの証明がもっとも知りたい場所だったのに、「わからなかった」では悲しすぎる。花田・松浦(2020)「ゼロからできるMCMC」p.128を見ながら再チャレンジ。
この本では、
- 提案確率\(q(x, y)\)のことを\(f(\{x\} \rightarrow \{x’\})\)
- 定常分布上の密度\(h(x)\)のことを\( e^{-S(x)} \)
- 遷移確率\(q(x, y)a(x, y)\)のことを\(T(x \rightarrow x’)\)
- 可逆性のことを詳細釣り合い条件
と呼んでいるので、いちいち読み替えながらメモする。
可逆性とは \(h(x) q(x, y) a(x, y) = h(y) q(y, x) a(y, x)\)のことです。
\(q(x, y) = q(y, x) = 0\)であれば遷移確率は\(q(x, y)a(x, y) = q(y, x)a(y, x) = 0\)なので、可逆性は自明です。そこで以下では、\(q(x, y) \gt, q(y, x) \gt 0\)の場合を考えることにしましょう。遷移確率は$$ q(x, y)a(x, y) = q(x, y) \times \min \left(1, \frac{h(y) q(y, x)}{h(x) q(x, y)} \right) $$ $$ q(y, x)a(y, x) = q(y, x) \times \min \left(1, \frac{h(x) q(x, y)}{h(y) q(y, x)} \right) $$ です。
\(h(y) q(y, x) \geq h(x) q(x, y)\) のときは、[\(x\)から\(y\)への遷移のHastings比は1を越え, 逆の遷移のHastings比は1を下回るから] $$ h(x) q(x, y) a(x, y) = h(x) q(x, y) \times 1$$ $$ h(y) q(y, x) a(y, x) = h(y) q(y, x) \frac{h(x) q(x, y)}{h(y) q(y, x)} $$ なので、$$ h(x)q(x,y)a(x,y) = h(y)q(y,x)a(y,x) $$ が成り立っており、可逆性が成り立っています。\(h(y) q(y, x) \le h(x) q(x, y)\) の場合も\(x, y\)を入れ替えて同じ計算を繰り返せば可逆性が示せます
なるほど、これなら理解できる…
練習のために、受理関数をBarker(1965)の受理関数 $$ a_B(x, y) = \frac{h(y)q(y,x)}{h(y)q(y,x) + h(x)q(x,y)}$$ に差し替えた場合について考えてみよう。$$ h(x)q(x,y)a_B(x,y) = \frac{h(x)q(x,y) \times h(y)q(y,x)}{h(y)q(y,x) + h(x)q(x,y)}$$ $$ h(y)q(y,x)a_B(y,x) = \frac{h(y)q(y,x) \times h(x)q(x,y)}{h(y)q(y,x) + h(x)q(x,y)}$$ だから、やっぱり可逆性がなりたっている。なるほど。
12.3 Metropolis更新
すべての\(x, y\)について\(q(x, y) = q(y, x)\)である場合、つまり$$ r(x, y) = \frac{h(y)}{h(y)} $$ になっている場合をMetroplis更新という。ちょっぴり計算が速くなるけど、ほかにメリットはない。
\(q(x, y) = q(y, x)\)となる例のひとつは、\(e\)を\(x\)とは独立で平均0の確率変数にして\(y = x+e\)とするやり方である。\(e\)を正規分布とか一様分布とかにして使うことが多い。
12.4 Gibbs更新
ほしい均衡分布の条件付き分布から提案をドローし、必ず受理するのがGibbs更新である。
これは可逆である。なぜなら、ほしい定常分布が\(X_n\)だとして、\(f(X_n)\)の下での\(X_{n+1}\)の条件付き分布が\(f(X_n\)の下での\(X_{n}\)の条件付き分布と同じなら、\((X_n, X_{n+1})\)は\(f(X_n\)\)のもとで条件付きに交換可能だからである。従って\((X_n, X_{n+1})\)は無条件に交換可能である。
一般的な用語としては、状態ベクトルのある要素の、他の要素の下での条件付き分布を使う場合をGibbs更新という。この場合は、その要素のみ除外して\(f(X_n)\)を使っていることになる。この形の条件付き分布を”full conditional”という。こういう条件つき分布が好まれているは、伝統としかいいようがない。[←へええええ！]
他の方法も検討されている。いくつかの要素を除外して\(f(X_n)\)を使うのを「ブロックGibbs」という。これもまた、好む理由は伝統のほかにない。
もしfull conditionalのみを指してGibbs更新と呼ぶべきだという人がいるのなら、その人はここで述べているGibbs更新を「一般化Gibbs更新」とでも呼ぶのかもしれない。全然一般化じゃないんだけど。
Gibbs更新は、M-H更新にはない面白い特徴を持っている。複数の更新の効果が、一回の更新の効果と同じだという特徴である。これをidempointという。なぜかというと、更新によって\(f(X_n\)が変わるわけではないので、何回やっても、\(f(X_n)\)の下での\(X_{n+1}\)の条件付き分布は1回やったときと変わらない。Gibbs更新はそれ自体は役に立たないもので、なにか別の更新と組み合わせる必要がある。
12.5 variable-at-a-time Metroplis-Hastings
M-H更新を修正して、一回の更新では状態空間の一部だけを変えるようにすることもできる。これをvariable-at-a-time M-H更新と呼ぼう。
[…中略…]
この呼び方は実はちょっと誤っている。Metroplis et al.(1953)が提案したのはvariable-a-time サンプラーだったからだ。歴史的には、Metropolisのアルゴリズムは、本書でいうところのM-H更新、そしてvariable-at-a-time M-H更新を含んでいた。
12.6 GibbsはM-Hの特殊ケースである
[略]
12.7 更新の合成
更新メカニズム(というかコード)が\(k\)種類あるとする。それぞれを\(P_1, \ldots, P_k\)とし、順番に適応していくとしよう。それぞれの更新メカニズムが分布を保持するなら、その組み合わせ\(P_1 P_2 \ldots P_k\)も分布を保持する。
各メカニズムが、ほしい均衡分布の”full conditionals”を使ったGibbs更新だったら、組み合わせのことをfixed can Gibbs サンプラーという。
[…中略…]
もし\(P_1 P_2 \ldots P_k\)の遷移確率と\(P_k P_{k-1} \ldots P_1\)の遷移確率がちがっていたら可逆性が失われる。しかし\(P_i = P_{k-1}\)としておけば可逆性は保たれる。こういう組み合わせのことをpalindromicだという。
12.8 状態独立な混合
更新メカニズム(コード) \(P_y\)があって、それらのメカニズムからランダムにひとつを選んで更新することを\(E(P_Y)\)としよう。これを混合と呼ぶことにする。
選ばれるメカニズム\(Y\)が現在状態と独立に決まり、かつそれぞれの\(P_y\)が同じ分布を保存するならば、\(E(P_Y)\)も同じ分布を保存する。\(X_n\)の分布がほしい均衡分布だったら\(X_{n+1}\)の分布も無条件にそうだということになる。さらに、それぞれの\(P_y\)が可逆なら\(E(P_Y)\)も可逆である。
たとえば、状態ベクトルの中から選ぶ要素を変えた”full conditional” Gibbs更新たちを用意し、一様分布で選ぶ、ということがある。これをrandom scan Gibbsサンプラーという。
合成と混合は組み合わせてもよい。[…中略…]
我々はこういう”scan”技術があまり好きでない。そういう技術は合成と混合の組み合わせのうちいくつかの特殊ケースだけに注目しているが、注目しているケースに特に良い点があるわけではない。それらが注目されている理由はただの伝統でしかない。[…]
12.9 サブサンプリング
[略]
12.10 Gibbs と Metropolis を再訪する
本書の用語法はちょっと世間と違うので、世間の用語法を整理しておこう。
- Gibbsサンプラー: すべての基本更新がGibbsで、なんらかの合成とか(fixed scan Gibbs), 混合とか(random scan Gibbs), どの両方とか(random sequence scan Gibbs)を伴っているMCMCサンプラー。
- Metropolisアルゴリズム: すべての基本更新がMetropolisで、合成や混合やその両方を伴っているもの。
- Metropolis-Hastingsアルゴリズム: すべての基本更新がM-Hで、合成や混合やその両方を伴っているもの。
- Metropolis-within-Gibbsサンプラー: 上と同じ。なぜならGibbsはM-Hの特殊ケースだから。
- 独立Metropolis-Hastingsアルゴリズム: M-Hアルゴリズムで、提案状態が現在状態に依存しないケース。
- ランダムウォークMetropolis-Hastingsアルゴリズム: M-Hアルゴリズムで、提案が\(x – e\)の形になっていて、\(e\)が\(x\)と確率的に独立なもの。
Gibbsサンプラーはナイーブなユーザにいまでも人気がある。もし私が、MCMCの収束が遅いんで助けてくれと頼まれてGibbsを使うのをやめろとアドバイスするたびに5セントもらっていたら、今頃大金持ちになっていた。
Gibbsが使われている理由の一つは自動的であることだろう。M-Hだと提案分布を選ばないといけないし、正規ランダムウォークMHだと決めたら今度は分散行列を選ばないといけない。しかし自動的であるというのは幻想にすぎない。もしscan技術について知っていたらfixed scanにするかrandom scanにするか選ばないといけないし、block Gibbsとか、もっと一般的なGibbsもある。
にもかかわらず、GibbsはM-Hより自動的だと思っているユーザが多い。問題は、選択肢がないことが良いことなのか悪いことなのか、だ。それがいいことだというならGibbsでいいけど、もしそうでないならGibbsはよくない。
13. Metropolisの事例
[mcmcパッケージのmetrop()を使ったrandom-walk Metropolisの説明。親切なコード例がついている。メモは省略するけど、なるほどな、忘れないようにしようと思った箇所についてメモしておく]
二値の反応ベクトル y、モデル行列 x、ロジスティック回帰係数 betaがあって、対数尤度を得たいとしよう。
eta <- as.numeric(x %*% beta) として、$$ p = \frac{\exp(\eta)}{1+\exp(\eta)} = \frac{1}{1+\exp(-\eta)}$$ $$q = \frac{1}{1+\exp(\eta)} = \frac{-\exp(-\eta)}{1+\exp(-\eta)}$$ の対数が入ったベクトル logpとlogqを求めたい。そうすれば、sum(logp[y == 1]) + sum(logq[y == 0])を求めればいいでしょ。
こういうとき、そのまま logp <- eta - log(1+exp(eta)) という風にコードを書くのは良くない。\(\eta\)が正の大きな値の時にexp(eta)がオーバーフローするからだ。こういうときは、\(\eta\)の符号で場合分けするとよい。負ならば、$$\log(p) = \eta - \log(1+\exp(\eta))$$ $$\log(q) = -\log(1+\exp(\eta))$$を使う。正ならば$$\log(p) = -\log(1+\exp(-\eta))$$ $$\log(q) = -\eta - \log(1+\exp(-\eta)) $$を使う。
それに、log(1+x)じゃなくてlog1p(x)を使ったほうがよい。xが0に近いときも正確になる。
というわけで、こう書きましょう。
logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
14. チェックポイントを設ける
[連鎖の途中で出力するRコード例。パス]
15. MCMCコードのデザイン
MCMCアルゴリズムのデザインほど簡単なものはない。トイ・プロブレムを別にすれば、すべての問題はMetropolis-Hastings-Greenアルゴリズムでよい。
新しいサンプラーを開発したら、それがM-H-Gアルゴリズムの特殊ケースであることを証明すれば良い。難しくはないが頭を使わないといけない(証明というのはそういうものだ)。よくある誤りは、マルコフチェーンの状態がどうなるかをきちんと考えないことだ。たとえば、ある提案分布がある変数に依存しているのに、Hastings比を計算する際にそれを考慮しないことである。つまり、ある部分について考えているときの状態空間と別の部分について考えているときの状態空間を変えてしまうことだ。[...中略...]
[コードを書くときの心かまえ。restart出来るようにしろ、バッチ平均を出せるようにしろ、云々]
16. MCMCコードの妥当性確認とデバッギング
[出力だけチェックしても無意味だ。内部変数の動きをちゃんと調べろ。云々]
----
やれやれ、先は長いな...