マルコフ連鎖モンテカルロ法(MCMC)

 ある確率密度関数 f(z) が得られたけれど、その平均値や標準偏差、最大値などを解析的に計算できないことがよくあります。数値積分すれば良さそうですが、今度は数値積分の方法が問題になります。統計量や積分を直接計算する代わりに、もし、確率分布 f(z) に従うサンプル  z_1 , z_2 , \cdots , z_N が得られれば、ある量 g(z) の統計平均は  \langle g(z) \rangle = \int{g(z)f(z)dz}=\frac{1}{N} \sum_n {g(z_n)} と近似計算できます。

  最も単純なサンプリング法は、モンテカルロ法の一種である棄却サンプリングです。棄却サンプリングのうち、提案分布を一様分布としたものが、最も簡単なモンテカルロ法です。これは、確率変数  z を一様分布からサンプルし、得られた  z を確率  f(z)/M で採用し、確率  1-f(z)/M で棄却するものです(M は f(z) の最大値以上の値)。しかし、この方法では、棄却が非常に多くなってしまう問題や、適切な提案分布を選ぶのが難しいという問題があります。

マルコフ連鎖モンテカルロ

 以降、確率分布を離散分布とし、 p(z) と表します。 z 1〜k までの値をとる  k 状態変数とします。 \int{g(z)f(z)dz} = \sum^{k}_{z=1} {g(z)p(z)} という対応関係があります。


 ステップ  t に対して、変数が確率的に変化することを考えます。ステップ  t における変数の値が、それ以前の  m ステップ、すなわちステップ  t-m〜t-1 における変数の値に依存するとき、これを  mマルコフ過程と呼びます。


 特に1次マルコフ過程では、ステップ  t で変数が  z' であったときに、ステップ  t+1 z に遷移する確率  q(z|z') が定義されます。 q(z|z') は、マルコフ核とも呼ばれます。ステップ  t における確率分布を  q_t(z) と書くと、次の関係が成り立ちます。


 \begin{align} q_{t+1}(z)=\sum^{K}_{z'=1} q(z|z')q_t(z') \end{align} \tag{1}


 一般に、ステップ  t での分布  q_t(z) とステップ  t+1 での分布  q_{t+1}(z) とは異なる分布になりますが、ステップを進めることで特定の分布に収束し、それ以上変化しなくなることがあります(詳細は後述)。この分布を定常分布と呼び、ここでは  q_{\infty}(z) と書きます。すなわち、


 \begin{align} q_{\infty}(z)=\sum^{K}_{z'=1} q(z|z')q_{\infty}(z') \end{align} \tag{2}


となります。


 もし、当初の目的としていた確率分布  p(z) が定常分布  q_{\infty}(z) に一致するならば、適当な初期値  z_1 からはじめて、定常分布に収束するまで  z_1 , z_2 , \cdots とサンプリングを繰り返すことで、最終的に  p(z) に従う多数のサンプルを得ることができます。もちろん、定常分布に収束する前の前半のサンプルは捨てなければなりません。また、サンプル間の相関が小さくなるように、採用するサンプル間のステップ幅を離すことが必要です。


 これが、マルコフ連鎖モンテカルロ法(MCMC)です。問題は、どのようにして  p(z) を定常分布に持つマルコフ過程を準備するかです。


詳細釣り合い条件

 次の関係が成立しているとき、マルコフ核  q(z|z') は分布  p(z) に対して詳細釣り合い条件を満たしていると言います。


 \begin{align} q(z'|z)p(z) = q(z|z')p(z') \end{align} \tag{3}


 両辺それぞれで  z' に関して和をとると、


 \begin{align} p(z) = \sum^{K}_{z'=1}(z|z')p(z') \end{align} \tag{4}


となります。これは、(詳細でない)釣り合いの式です。式(2) と比較すると、 p(z) はこのマルコフ過程の定常分布になっていることがわかります。定常分布に対して、(詳細でない)釣り合いの式(2),(4) は必ず成り立ちますが、詳細釣り合いの式(3) は、可逆過程に対応する特殊な条件であり、常に成り立つとは限りません。


 詳細釣り合い条件が満たされていれば、確率分布  q_t(z) t \rightarrow \infty で目的の分布  p(z) に収束します(詳細は後述)。残る問題は、詳細釣り合い条件を満たすマルコフ核をどうやって求めるかです。


メトロポリスヘイスティングス法(MHアルゴリズム

 MHアルゴリズムは、適当な提案分布  \tilde{q}(z|z') からはじめて、詳細釣り合い条件を満たすマルコフ核  q(z|z') を作ることで、目的の分布  p(z) からのサンプリングを行うアルゴリズムです。


 アルゴリズムは次のようなものです。

  1. 適当な事前分布に従って、初期値  z_1 をサンプルする。
  2. ステップ  t において値  z' をサンプルした時、次のステップでのサンプル候補  z を、適当な提案分布  \tilde{q}(z|z') に従ってサンプルする。
  3. 次の値を計算する。 \rho(z, z') = \frac{\tilde{q}(z'|z) p(z)}{\tilde{q}(z|z') p(z')} \tag{5}
  4.  \rho(z,z') \geq 1 なら、サンプル候補  z を次のステップ  t+1 でのサンプル値とする。 \rho(z,z') < 1 なら、確率  \rho z を次のステップ  t+1 でのサンプル値とし、確率  1 - \rho(z,z') z' (遷移前と同じ値)を次のステップ  t+1 でのサンプル値とする。


 このアルゴリズムによって生成されたマルコフ連鎖は、目的の分布  p(x) に収束し、安定分布となります。提案分布  \tilde{q}(z|z') は何でも良いのですが、 z' を中心としたガウス分布等の単峰性の分布や矩形分布が用いられます。中心からの幅が狭すぎると、隣接ステップの相関が大きくなり、無駄なサンプルが増えます。一方で中心からの幅が広すぎると、マルコフ性の利点が失われてしまいます。提案分布の形状(単峰性か多峰性か等)が目的の分布に類似しているほど、マルコフ性の利点が増します。

MHアルゴリズムと詳細釣り合い条件との関係

 それでは、MHアルゴリズムが詳細釣り合い条件を満たすマルコフ核を作ることを見てみましょう。アルゴリズムマルコフ過程として書き下すと、 z \neq z' の場合は

 \begin{align} q_{t+1}(z)=\sum^K_{z'=1}  \min{(\rho(z, z'), 1)}\tilde{q}(z|z')q_{t}(z) \end{align} \tag{6}

となります。式(1) と比較することで、 z \neq z' の場合のマルコフ核は

 \begin{align} q(z|z') =  \min{(\rho(z, z'), 1)}\tilde{q}(z|z') \end{align} \tag{7}

であることが分かります。 q(z|z') は確率分布ですから、 z = z' の場合は

 \begin{align} q(z'|z') = 1 - \sum_{z \neq z'} q(z|z') \end{align} \tag{8}

となります。


 式(7) のマルコフ核の性質を調べます。式(5) から、 \rho(z, z')=\frac{1}{\rho(z', z)} です。従って、 \rho(z, z') \rho(z, z') の一方は 1 以上で、他方は 1 以下です。例えば  \rho(z, z') \leq 1 ならば

  \begin{align}  q(z|z') &=&  \rho(z, z') \tilde{q}(z|z') \\
 q(z'|z) &=&  \tilde{q}(z'|z) \end{align} \tag{9}

という関係が成り立っており、 \rho(z,z') に式(5) を代入すると、確かに詳細釣り合い条件  q(z'|z)p(z) = q(z|z')p(z') が成立しています。

ギブスサンプリング

 詳細釣り合いに依拠するMHアルゴリズムは汎用的で、目的の分布に確実に収束する方法です。しかし、変数  z の次元が大きければ、マルコフ核の底空間も巨大になり、サンプリングが追いつかなくなります。ギブスサンプリングはMHアルゴリズムほど汎用的ではありませんが、高次元空間においては非常に強力なサンプリング法です。


 ギブスサンプリングは、目的の分布から各変数ごとの条件付き分布を求め、各変数に対して逐次的にサンプリングを繰り返し、収束を狙うものです。まず、 z=(x_1, x_2, \cdots, x_d) として、目的の確率分布  p(x_1, x_2, \cdots, x_d) から、次のように変数ごとの条件付き分布を求めます。


 \begin{align}  p(x_1 | x_2, x_3, \cdots, x_d), \ \ \ 
 p(x_2 | x_1, x_3, \cdots, x_d), \ \ \ 
  \cdots ,  \ \ \ 
 p(x_d| x_1, x_2, \cdots, x_{d-1}) \end{align} \tag{10}


 そして、各変数を順番にサンプリングしながら更新し、収束するまで繰り返します。



もう少し詳しい議論


 これまでの説明の中には、いくつかごまかしがあります。下記の項目について、もう少し詳しく見てみましょう。

  • 安定分布は、遷移確率(マルコフ核)だけによって表現できるのか
  • 詳細釣り合い条件が成り立てば、なぜ定常分布に収束するのか
  • 詳細釣り合い条件が成り立たないのはどのような場合か
マルコフ過程と行列方程式

 マルコフ核は、 K_{zz'}=q(z|z') という  k\times k 行列  K で表せます。確率変数の分布  q_t(z) も同様に  k 次元ベクトル  q_t で表せます。式(1) は

 \begin{align} q_{t+1}=Kq_t \end{align} \tag{11}

と表せます。定常分布を  K 次元ベクトル  p と書くと、定常分布を表す行列方程式

 \begin{align} p=Kp \end{align} \tag{12}

が得られます。これが自明でない解を持つので、 I-K の各行(または列)は線形独立ではありません。実際、下記の拘束条件が含まれています。

 \begin{align} \sum_z K_{zz'} &=& 1 \\
\sum_z p_z &=& 1 \end{align} \tag{13}

 そこで、 p_k = 1 - p_1 - p_2 - \cdots - p_{k-1} を用いて、式(12) から  p_k を消去すると、次の行列方程式が得られます。

 \begin{align} p'=K'p'+b' \end{align} \tag{14}

 ただし、 K' (k-1)\times (k-1) 行列、 p', b' はいずれも  k-1 次元ベクトルで、

 \begin{align} p'_z &=& p_z \\
K'_{zz'} &=& K_{zz'} - K_{zk} \\
b'_z &=& K_{zk} \end{align} \tag{15}

を満たします。式(14) の各行は独立に与えられるので、 I-K'逆行列を持てば、

 \begin{align} p'=(I-K')^{-1}b' \end{align} \tag{16}

によって、遷移確率から安定分布を求めることができます。

詳細釣り合い条件と収束性

 式(1) で表されるマルコフ過程において、式(3) の詳細釣り合い条件が満たされていれば、 q_t(z) p(z) に収束することを示します。そのために、 q_t(z) p(z) の間の、 p(z) を基準としたカルバック・ライブラー・ダイバージェンス  D( p||q_t ) が、ステップ毎に単調に減少することを確認します。


 \begin{align} D(p||q_{t+1}) &=& -H(p) - \sum_z p(z) \log{q_{t+1}(z)}  \\
&=& -H(p) - \sum_z p(z) \log{\sum_{z'}{q(z|z')q_t(z')}} \\
&=& -H(p) - \sum_z p(z) \log{\sum_{z'}{\frac{q(z'|z)p(z)}{p(z')}q_t(z')}} \\
&=& - \sum_z p(z) \log{\sum_{z'}{\frac{q(z'|z)q_t(z')}{p(z')}}} \\
&\leq & - \sum_z p(z) \sum_{z'}{q(z'|z)\log{\frac{q_t(z')}{p(z')}}} \\
&=& - \sum_{z, z'} p(z')q(z|z')\log{\frac{q_t(z')}{p(z')}} \\
&=& - \sum_{z'} p(z')\log{\frac{q_t(z')}{p(z')}} \\
&=& D(p||q_t) \end{align} \tag{17}


となり、確かに  D(p||q_{t+1}) \leq D(p||q_t) となっています。 q_t(z) p(z) に一致するまで、 D(p||q_t) が減少し続けることが分かります。ただし、より厳密には、減少速度を評価する必要があります。 Hエントロピーを表します。なお、上記の式変形において、1行目から2行目へは式(1) を用い、2行目から3行目へは詳細釣り合いの式(3) を用いました。4行目と5行目の間の不等号は、イェンゼンの不等式によります。5行目から6行目へは再び詳細釣り合いの式(3) を用いました。

詳細釣り合い条件と自由度

 マルコフ核は、式(13) の拘束条件の範囲で自由な正値を取ることができます。式(13) の拘束条件は式(8) に反映されていますから、一般にマルコフ核の自由度は、 z=z' を除き、 k^2 - k 個です。


 一方で、詳細釣り合い条件が満たされている場合、式(3) より q(z,z')=\frac{q(z',z)p(z)}{p(z')} なので、 q(z',z) p(z) が決まれば  q(z', z) も決まります。従って、この場合の自由度はマルコフ核のうち  z < z' であるものと  p(z) のうち一つを除いたものの個数の和であり、 \frac{k(k+1)}{2}-1 個です。


 自由度が詳細釣り合い条件の有無によらないのは、 k=1, 2 の場合のみです。 k=2 の場合は、 q(2|1)p(1)=q(1|2)p(2) が必ず成り立ちます。 k \geq 3 では、詳細釣り合い条件を満たさないマルコフ過程が存在します。