マルコフ連鎖モンテカルロ法(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}q(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 では、詳細釣り合い条件を満たさないマルコフ過程が存在します。



ハミルトニアンモンテカルロ法


 メトロポリスヘイスティングス法(MHアルゴリズム)では、遷移確率(マルコフ核)の提案分布  \tilde{q}(z|z') をどのように決めるかが問題です。 \tilde{q}(z|z') の決め方の一つが、ハミルトニアンモンテカルロ法です。


 そもそも、MHアルゴリズムの目的は、サンプリングを効率よく行うために、モンテカルロ法での棄却を抑制することでした。MHアルゴリズムにおいて、式(5) で表される  \rho(z, z') を1 に近づけることができれば、棄却率は低くなります。


 ハミルトニアンモンテカルロ法では、ハミルトン方程式に従って、時刻  t での位置座標  z' から時刻  t+1 での位置座標  z への変化をシミュレートします。通常、ハミルトン方程式による時間発展は確率的ではなく、確定的な変化です。そのため、運動量の初期値をサンプリングすることによって、確率的な遷移を表現します。このことから分かるように、マルコフ核  \tilde{q}(z|z') は明示的な関数形では与えられません。運動量の初期値の分布が決まっていて、そこから得られた運動量のサンプルと位置の初期値  z' を、ハミルトン方程式に従う時間発展という複雑な関数に通すことで、遷移後のサンプル  z が確率的に得られます。


 ハミルトニアンモンテカルロ法のポイントは、ハミルトニアンが次の二つの性質を持っていることです。一つめは、ハミルトニアンによる時間発展は可逆であるという性質です。二つめは、ハミルトニアンに従う時間発展において、エネルギーの値(ハミルトニアンに位置と運動量を代入した値)は不変であるという性質です。特に、二つめの性質によって、式(5) の  \rho (z,z') の値が 1 に近くなります。


 ハミルトニアンモンテカルロ法アルゴリズムを説明します。以下、 z を位置座標として扱い、運動量を  P とします(確率分布の記号ではないことに注意してください)。 まず、ハミルトニアン  H(k,z)=K(P)+U(z) は次のように定義します。

 \begin{align} K(P)=\frac{P^2}{2} \\
U(z)=-\log{p(z)} \tag{18} \end{align}

  p(z) は、目的の確率分布です。運動エネルギーの項は質量 1 の粒子に対応し、位置エネルギーの項は逆温度 1 に対応しています。


 運動量  P の初期値は、平均 0、標準偏差 1 のガウス分布からサンプリングします。なぜなら、 \exp{(-H(P,z))}=p(z)\exp{(-\frac{P^2}{2})} だからです。従って、運動エネルギーの項が質量  m の粒子を表すように決めることは、標準偏差  \sqrt{m^{-1}}ガウス分布を選ぶことに対応します。


 シミュレーションは、ハミルトン方程式を数値積分することで行います。ハミルトン方程式は以下の通りです。

 \begin{align}  \frac{dz}{dt} = \frac{\partial H}{\partial P} \\
 \frac{dP}{dt} = -\frac{\partial H}{\partial z} \end{align} \tag{18}

これから、エネルギーが時間不変であることは、次のようにして分かります。

 \begin{align} \frac{d H(P, z)}{dt}=\frac{dP}{dt} \frac{\partial H}{\partial P} + \frac{dz}{dt} \frac{\partial H}{\partial z}=0 \end{align} \tag{19}

 可逆性については、 t^*=-t, P^*=-P, z^*=z と置き換えても式(18) が成立することから分かります。


 式(11) を数値積分する方法としては、ルンゲ・クッタ法が最も汎用的ですが、ハミルトン力学系ではシンプレクティック法が良い性質を示します。ここでは、可逆性を持つ二次のシンプレクティック法を記載します。 \tau を時間ラベルとして、十分小さい値  \epsilon に対し、以下のように時間発展を行います。

 \begin{align} z(\tau + \frac{\epsilon}{2})=z(\tau)+\frac{\partial K (P(\tau))}{\partial P} \frac{\epsilon}{2}\\
P(\tau + \epsilon)=P(\tau)-\frac{\partial U(z(\tau + \frac{\epsilon}{2}))}{\partial z}\epsilon \\
z(\tau + \epsilon)=z(\tau + \frac{\epsilon}{2})+\frac{\partial K (P(\tau + \epsilon))}{\partial P} \frac{\epsilon}{2}
\end{align} \tag{20}

 この計算を  L 回繰り返すことで、時刻  \epsilon L だけ先の z が得られます。この値を、提案分布による遷移先のサンプルとします。 L=1 でも構いません。


 結局、ハミルトニアンモンテカルロ法アルゴリズムは次のようになります。

  1. サンプルの初期値  z_1 を決める。
  2. サンプル  z' が得られたとき、次の値をサンプリングするための運動量の初期値  P(0) を、平均 0、標準偏差 1 のガウス分布からサンプルする。
  3. 位置座標の初期値  z(0)=z' と運動量座標の初期値  P(0) から、式(20) の更新式を  L 回繰り返し適用する。
  4. 得られた  z(\epsilon L) を、サンプル  z の候補とし、 \rho(z,z')=\exp{[-\{H(P(\epsilon L), z(\epsilon L))-H(P(0), z(0))\}]} を計算する。
  5. 確率  \min{(\rho (z,z'), 1)} でサンプルを採択し、 z=z(\epsilon L) とする。確率  1-\min{(\rho (z,z'), 1)} でサンプルを棄却し、 z=z' とする。

 アルゴリズムの中で、 \rho(z,z')=\exp{[-\{H(P(\epsilon L), z(\epsilon L))-H(P(0), z(0))\}]} が非自明に見えるかもしれません。 \rho(z,z') の定義は、式(5) で与えられていました。そして、最終的な遷移確率が運動量の初期値の確率分布に依存することから、 \tilde{q}(z|z')=\frac{1}{\sqrt{2\pi}}\exp{(-\frac{P(0)^2}{2})} です。一方で、ハミルトン方程式の可逆性より、 \tilde{q}(z'|z)=\frac{1}{\sqrt{2\pi}}\exp{(-\frac{P(\epsilon L)^2}{2})} も成り立ちます。