マルコフ連鎖モンテカルロ法(2)Gibbsサンプラー

はじめに

引き続き、非常に多く利用されるMCMCアルゴリズムの1つであるGibbsサンプラーについて述べる。
Metropolis-Hastings(MH)アルゴリズムにおいて、多次元の提案分布からサンプリングすることが難しい場合に、サンプリングしたい変数をいくつかのまとまりに分割する多重ブロックMHアルゴリズムと呼ばれる方法があり、Gibbsサンプラーはその特別な場合として考えることができる。あとで見るように、Gibbsサンプラーでは常に提案されたサンプルが採択されることになる。サンプリングに使用する分布を特定することができるのであれば、MHアルゴリズムよりも実装が簡単だと個人的には感じる。

多重ブロックMHアルゴリズム

いま2次元のベクトル xをサンプリングしたいものとする。もし2次元確率分布を用いたMHアルゴリズムが困難な場合には、 x = (x_1, x_2)に分割してサンプリングを行う。そのために、 x^{(t)}の各要素の提案分布 q_1(x^{'}_1|x^{(t)}_1, x_2), q_2(x^{'}_2|x^{(t)}_2, x_1)を定義する。例えば前者の q_1は、「 x_2の現在の状態を固定した時に、 x_1の現在の値 x^{(t)}_1のもとで x^{'}_1を提案する確率」である。
さらに、MHアルゴリズムと同様、ブロックそれぞれの採択確率を次のように定める。

 \displaystyle
\begin{eqnarray}
\alpha(x^{'}_1|x^{(t)}_1, x_2) = min \left[ \frac{\pi(x^{'}_1|x^2)q_1(x^{'}_1|x^{(t)}_1, x_2)}{\pi(x^{(t)}_1|x^2)q_1(x^{(t)}_1|x^{'}_1, x_2)}, 1 \right] \\
\alpha(x^{'}_2|x^{(t)}_2, x_1) = min \left[ \frac{\pi(x^{'}_2|x^1)q_2(x^{'}_2|x^{(t)}_2, x_1)}{\pi(x^{(t)}_2|x^1)q_2(x^{(t)}_2|x^{'}_2, x_1)}, 1 \right]
\end{eqnarray}
一般に Kブロックに分割する場合も同様である(一般のKの場合には、必ずしも1変数ずつサンプリングを行う必要はなく、多変数のブロックを組み合わせてもよい)。サンプリングを行いブロックを更新する順番は、サンプリング全体を通して一貫していればどのような順番でもよい。

Gibbsサンプラー

Gibbsサンプラーは、多重ブロックMHアルゴリズムの特別な場合である。多重ブロックMHアルゴリズムでは、他のブロックで条件付けた提案分布を考え、そこからのサンプルが詳細釣り合い条件を満たすように決められた一定の確率で採択されることになるが、Gibbsサンプラーでは他のすべての変数で条件付けた目標分布そのものからサンプリングする(したがって各ブロックは1変数からなる)。つまり、MHアルゴリズムの条件付き提案分布について

 \displaystyle
q_k(x^{'}_k|x^{(t)}_k, x_{-k}) = \pi(x^{'}_k|x_{-k})  \\
q_k(x^{(t)}_k|x^{'}_k, x_{-k}) = \pi(x^{(t)}_k|x_{-k})
(ここで x_{-k}は、 k番目を除くすべての変数を意味する)であり、したがって提案分布からのサンプルは確率1で採択されるという場合である。
この方法は採択・棄却の手続きを踏むことなく、1次元の確率分布からのサンプリングを繰り返すことで目的の多次元分布をサンプリングできるため、実装は比較的容易と言えるだろう。ただし、「他のすべての変数で条件付けた目標分布」の具体的な形が分かっていないとサンプリングが困難という問題点がある(もしかしたら理論上は分布形を特定できるとしても、それを導出する数学力が足りなくてできないかもしれない)。そのような場合には、サンプリング困難なブロックだけMHアルゴリズムを用いればよい。このようにGibbsサンプラーとMHアルゴリズムを組み合わせる方法はMetropolis within Gibbsアルゴリズムと呼ばれる。

おわりに

Gibbsサンプラーの概要を紹介した。この手の話は抽象的な説明だけではどうしても理解が進みにくい。今後実例のシミュレーションも紹介していきたいと思う。

f:id:mstour:20200802190733p:plain
2変数の場合のGibbsサンプラー
参考文献