はじめてのMCMC (メトロポリス・ヘイスティングス法)
こんにちはtatsyです。
最近はSRMの解説記事ばかりでしたが、たまにはもう少し実践的なことを書こうかと思います。というわけで特に理由はないのですがMCMCについて解説したいと思います。
MCMCとは?
マルコフ連鎖モンテカルロ法(MCMC)は多次元空間内の点を確率分布に基づいてサンプリングするための方法で、最もシンプルな応用例は適当な関数と確率密度関数のペアに対して期待値を計算するというものです。
積分を計算するときに、もし期待値を求めたい関数と確率密度関数の両方が分かっているのであれば、台形公式やシンプソン公式みたいな区分求積っぽいやり方をすれば良いわけです。ところが、あまり次元が高次元だったりする場合には、そもそも区分求積のようなやり方が使えない場合があります。そんなときに、ある確率密度関数に従うようなサンプル点をうまく得るための方法がMCMCです。
MCMCはその名前の通り、マルコフ連鎖を用いたモンテカルロ法です。マルコフ連鎖では任意の確率変数からスタートして順々に確率変数を得ていった時に、それらの確率変数の分布が不変になる(=一定の分布に収束する)ことを仮定しています。この収束先の分布を不変分布(invariant distribution)といいます。
今、サンプルを得たい確率密度関数が$p(x)$であるならば、これを不変分布に持つようなマルコフ連鎖を生成することで、実際のサンプリングを行おうというのがMCMCの考え方です。MCMCにおいて、サンプルされる確率変数$x$の分布がある不変分布に近づくための十分条件として詳細つり合い条件という条件があり、多くのMCMCでは詳細つり合い条件をうまく使って目的のマルコフ連鎖を得ます。
詳細つり合い条件とは、ある状態$x^t$から確率$q(x^{t+1}|x^t)$で別の状態$x^{t+1}$に移るとき、確率密度関数$p$と$q$の間に次の式が成り立つことを言います。
$$ p(x^t) q(x^{t+1} | x^t) = p(x^{t+1}) q(x^t | x^{t+1}) $$
特に、新しい状態に移る確率を表す確率密度関数$q$を提案分布と呼びます。
メトロポリス・ヘイスティングス法
メトロポリス・ヘイスティングス法(M-H法)では、必ずしも詳細つり合い条件を満たさない提案分布$q(x|y)$に対して、詳細釣り合い条件を成り立たせるための採択確率$a(x|y)$が次の式を満たすと考えます。
$$ p(x^t) q(x^{t+1} | x^t) a(x^{t+1} | x^t) = p(x^{t+1}) q(x^t | x^{t+1}) $$
ここでポイントなのは、サンプリングを試みたい不変分布から直接サンプリングは難しいけれども、提案分布からなら容易にサンプリングが可能という場合を考えているということです。ここで、提案分布$q(x | y)$により得られた$x$を$a(x | y)$の確率で採択することで、詳細つり合い条件を満たすように棄却サンプリング (rejection sampling)を行おうというのがM-H法のアイディアです。
具体的には採択確率$a(x | y)$は
$$ a(x^{t+1} | x^t) = \min \left[ 1, \frac{p(x^{t+1}) q(x^t | x^{t+1})}{p(x^t) q(x^{t+1} | x^t)} \right] $$
となります。もし提案分布にガウス分布のような対象の分布を用いる場合には、不変分布の比だけを考慮すればよいことになります。このように対称な提案分布を使う方法を特にメトロポリス法と呼びます。
以下のコードでは提案分布がsigma
を標準偏差とするガウス関数であるとして、M-H法を実行しています。
M-H法と焼きなまし法
またM-H法と焼きなまし法(simulated annealing)を組み合わせると、より確率の値が大きいところを重点的にサンプリングすることもできます。
焼きなまし法と組み合わせるときには、系の温度Tを徐々に1から0に近づけていきながら、次の採択確率でサンプリングを行います。
$$ a(x^{t+1} | x^t) = \min \left[ 1, \frac{p(x^{t+1})^{1/T} q(x^t | x^{t+1})}{p(x^t)^{1/T} q(x^{t+1} | x^t)} \right] $$
以下が焼きなまし法を用いたコードです。
まとめ
ちなみに、ここまで紹介しておいてなんですが、プログラムは正しく動きますが、論理的な考察は怪しいのでBishop本などを見ていただけると助かります。
またこの記事のつづきでギブス・サンプリングについても解説しておりますので、もしよろしければそちらもご覧ください。
最後までお読みいただき、ありがとうございました。