マルコフ連鎖モンテカルロ法(1)Metropolis-Hastingsアルゴリズム

はじめに

ちょっと趣向を変えて、ベイズ統計モデルでの推定に広く用いられているマルコフ連鎖モンテカルロ法（MCMC）について書いていこうと思う。統計解析をする時は、例えばSASという統計ソフトで言えばproc mixedなどのプロシージャ（組み込みのパッケージのようなもの）で一発で答えが出るようなシンプルで解釈のしやすい方法で済ませたいと思うのが人情なのだけれど、そうはいかず、いわばオーダーメイドの方法で問題を解決しないといけないこともある。さて、自分でベイズモデルを組むとなると、事後確率分布を簡単な形で書き下すなんてことはまずできないであろう。そのような場合にはMCMCがとても役に立つ。MCMCの計算も今の時代基本的に統計ソフトに組み込まれていると思うが、原理をきちんと知っておくことが大事だと思う。（なお、先ほど書いたSASにも「proc MCMC」という機能があったりする。僕は使ったことないけれど・・・）
まずは、代表的なアルゴリズムの1つであるMetropolis-Hastingsアルゴリズムを見ていきたい。

Metropolis-Hastings(MH)アルゴリズム

いま、サンプリングを行いたい確率分布 $\pi(x)$ があるのだが、この分布の具体的な形状が不明であるため直接サンプリングするのが困難だとする。ベイズモデルでデータを得た後の事後確率分布が典型的な例である。このような時に、サンプリングが簡単に行えるような代理の確率分布からサンプリングを行うこととし、両方の分布の違いを調整しながら本来得たい分布からのサンプルを得るというのがMHアルゴリズムの発想である。なお本来関心のある分布を目標分布、代理に用いる分布を提案分布と呼ぶ。
このアルゴリズムは、マルコフ連鎖が不変分布に収束する（つまりどの状態から出発してもある一つの確率分布に収束する）ための条件の一つである「詳細釣り合い条件」を満たすようにして構成されている。したがって、代理の提案分布を用いたとしても、きちんと目標分布をシミュレーションすることができるわけである。詳細釣り合い条件とは、全てのとりうる値 $x$ に対して

$\displaystyle \pi(x^{(t)}) p(x|x^{(t)}) = \pi(x) p(x^{(t)}|x)$

が満たされることをいう。ただし $\pi(.)$ は目標分布、 $p$ はマルコフ連鎖を構成する確率分布であり、 $p(x|x^{(t)})$ は $x^{(t)}$ から $x$ へ推移する確率を表す。 $x^{(t)}$ は現在の時点 $t$ における状態とする。詳細釣り合い条件とは、値 $x^{(t)}$ から $x$ への推移が可能であれば、逆に $x$ から $x^{(t)}$ への推移も同じ確率で可能なことを意味している。
さて、提案分布を $q$ とすると、提案分布は必ずしも目標分布に収束するようなマルコフ連鎖ではないので、詳細釣り合い条件を満たしているとは一般には言えない。そのような状況は

$\displaystyle \pi(x^{(t)}) q(x|x^{(t)}) > \pi(x) q(x^{(t)}|x)$

と表せる。この状況では、 $x^{(t)}$ は $x$ に推移することが多いが、逆に $x$ から $x^{(t)}$ に推移することは少ない。
詳細釣り合い条件を満たすようにするためには、提案分布 $q$ を調整する必要がある。そのために推移を調整する確率 $\alpha$ を用いて、以下のように定義した $p$ が詳細釣り合い条件を満たすようにする：

$\displaystyle p(x|x^{(t)}) = q(x|x^{(t)}) \times \alpha(x|x^{(t)})$

つまり、

$\displaystyle \begin{eqnarray} \pi(x^{(t)}) p(x|x^{(t)}) & = & \pi(x) p(x^{(t)}|x) \nonumber \\ \pi(x^{(t)}) q(x|x^{(t)}) \alpha(x|x^{(t)}) & = & \pi(x) q(x^{(t)}|x) \alpha(x^{(t)}|x) \nonumber \\ & = & \pi(x) q(x^{(t)}|x) \end{eqnarray}$

とする。ただし最後の等号は、 $x$ から $x^{(t)}$ への推移が少ないことから確率 $\alpha(x^{(t)}|x)$ を最大値の $1$ に基準化したことによる。よって、 $x^{(t)}$ から $x$ への推移確率を調整する $\alpha(x|x^{(t)})$ は

$\displaystyle \alpha(x|x^{(t)}) = \frac{ \pi(x) q(x^{(t)}|x) }{ \pi(x^{(t)}) q(x|x^{(t)}) }$

と設定すればよい。前述のように提案分布 $q$ における $x^{(t)}$ から $x$ への推移確率に $\alpha$ を掛けることで詳細釣り合い条件を満たすようにしているので、サンプリングの方式としては $q$ から得たサンプルを、さらに確率 $\alpha$ でのみ採用すればよいことになる。つまり、提案分布からの候補 $x$ をまずサンプリングし、以下の確率 $\alpha(x|x^{(t)})$ でその候補を採用する：

$\displaystyle \alpha(x|x^{(t)}) = min\left[\frac{ \pi(x) q(x^{(t)}|x) }{ \pi(x^{(t)}) q(x|x^{(t)}) }, 1\right]$

これを十分な回数繰り返せば、目標分布からのサンプルが得られることになる。アルゴリズムを整理すると、以下の手順となる。

初期値 $x^{(0)}$ を設定する。
$\tilde{x}$ を提案分布 $q(x|x^{(0)})$ からサンプリングする。
採択確率 $\alpha(x|x^{(0)})$ を計算する。
$[0, 1$ ]上の一様乱数 $u^{(1)}$ を発生させる。
$u^{(1)} \leq \alpha(x|x^{(0)})$ ならば $\tilde{x}$ を採択し $x^{(1)} = \tilde{x}$ とする。そうでなければ採択せず $x^{(1)} = x^{(0)}$ とする。
手順2に戻り、 $\tilde{x}$ を提案分布 $q(x|x^{(1)})$ からサンプリングする。以降手順2〜5を所定の反復回数繰り返す。