Gibbs Sampler and Metropolis Hastings Algorithm

先輩と二人でExample of Metropolis Hastings Algorithm - yasuhisa's blog付近の続き。

メトロポリスとかがどういうアルゴリズムかとかは分かった。実装もして欲しい乱数列を取得することもできた。が、なんでこれでうまく行くのかが分からない。よし、調べてみようというのが今回の流れ。ファイナンスのためのMCMC法によるベイズ分析の6章を参考にしてやりました。

マルコフ連鎖

応用確率論の付近で習ったけど大分忘れてしまっていた。以前の全ての状態で条件付けた確率が一つ前の期にしか依存しないよ、ということがとりあえずマルコフ連鎖の定義。これをうまく使ってあげると同時分布をchain ruleを使って
$\begin{align} f(x_0,\cdots,x_{t+1}) &= f_0(x_0) f(x_1|x_0) f(x_2|x_0,x_1) \times \cdots \times f(x_{t+1|x_0,\cdots,x_t}) \\ &= f_0(x_0) \prod^t_{s=0} f(x_{s+1}|x_s) \end{align}$
と書ける。ここで $f_0$ は $X_0$ の周辺密度。

推移核(transition kernel)

$f(x_{t+1}|x_t)$ をマルコフ連鎖の推移核(transition kernel)と呼ぶ。日本語で言えば、t期にこんな値が出た時のt+1期の密度がどうかってとですね。で、ちゃんとした感じで書くと $K(x_t,x_{t+1}) = f(x_{t+1}|x_t)$ みたいな感じになる。これを使うとさっきの同時分布は $f(x_0,\cdots,x_{t+1}) = f_0(x_0) \prod^t_{s=0} K(x_s,x_{s+1})$ で書き直される。

この推移核を使うと $X_{t+1}$ の周辺分布は $f_{t+1}(x_{t+1}) = \int_{\infty}^{\infty} f_0(x_0) K^t (x_0,x_{t+1}) d x_0$ と書ける(6.11式)。ここで重要なことはt期における周辺密度は $x_0$ と推移核さえ分かれば順番に求めていくことができる、ということである(後ろ向きに繰り返し適用していく)。

不変分布

これはとりあえず離散の状態を考えたほうが分かりやすいかもしれない。t期においてn個の状態(X1,X2,…Xn)があったとしよう。そこからt+1期にあるxという状態に行く確率は、X1からXに行く確率、X2からXに行く確率、…XnからXに行く確率を全部足してやればよい。そういうわけで $f(x) = \sum_{i=1}^n f(x_t) K(x_t,x) d x_t$ という感じで書き表わせる。で、これの状態が無限個あるような場合を考えれば(6.12)式が出てくるよね、という感じでごまかしておく。

エルゴート性

上に書いたようあ後ろ向きに $f_{t+1}(t+1)$ を書きくだしていく作業を無限回繰り返した時に不変分布に収束するならば、そのマルコフ連鎖はエルゴート性を持つというらしい。

Gibbs Sampler Algorithm

アルゴリズム6.4の説明が全てである。推移核 $K(\mathbf{\theta}^{(r-1)},\mathbf{\theta}^{(r)}) = f(\theta_2^{(r)}|\theta_1^{(r)}) f(\theta_1^{(r)}|\theta_2^{(r-1)})$ が不変分布であることが示されるので、Gibbs Sampler Algorithmは不変分布に収束する、という流れである。

Metropolis Hastings Algorithm

肝は「 $f(\mathbf{\theta}) > 0$ となる任意の $\mathbf{\theta} \in \mathbf{R}^k$ に対して $q(\mathbf{\phi},\mathbf{\theta}) > 0$ となる $\mathbf{\phi} \in \mathbf{R}^k$ が必ず存在すると仮定する。」という一文。この後は

Metropolis Hastings Algorithmの推移核を示す
その推移核の不変密度がfであることを示す

という流れになっていて、その流れでアルゴリズムも導出できるような流れになっている。前回のゼミでは「酔歩連鎖とか独立連鎖とかやってうまくいったけど、なんであれでうまくいったの?他の方法をやるためにはどういう性質を満たすものだったらいいの?」という疑問があったんだけど、さっきの一文が答えとなっている。

前回の独立連鎖のところでは正規乱数からのものを利用していた。 $q(\mathbf{\phi},\mathbf{\theta}) = g(\mathbf{\theta})$ としておいて、gを正規分布だとする。正規分布は全実数において0より大きい。そういうわけで $\mathbf{\phi}$ が何であれさっきの一文の条件を満たすものとなっている。条件を満たしているので、独立連鎖の方法は欲しい乱数列を生成できる、ということになる。じゃあ、t分布とかからも生成できるはずだよね、やってみようというのは次回までにやってみようかな。→やってみた。