メトロポリスヘイスティング法(MH法) - 機械学習基礎理論独習

本記事では以下で「メトロポリスヘイスティング法」を「MH法」と表記します。

詳細釣り合い条件

マルコフ連鎖の記事で取り上げたのは遷移核が明らかな状態で定常分布を求める問題でした。
しかし、やりたいことは、事後分布が明らかな状態で、事後分布に従う乱数を手に入れることです。
事後分布が定常分布になるように、遷移核を設計する必要があります。

サンプリングしたい分布に対して、それを「定常分布」とするマルコフ連鎖を構成する方法を
マルコフ連鎖モンテカルロ法(以下、MCMC法)といいます。

マルコフ連鎖が定常分布に収束させるには、詳細釣り合い条件を満たせばよいです。
詳細釣り合い条件は、任意の $\theta,\theta'$ に対して、

$\begin{eqnarray} f(\theta|\theta')f(\theta')=f(\theta'|\theta)f(\theta)\tag{1} \end{eqnarray}$

が成り立ちます。

$(1)$ の $f(|)$ が遷移核(未知)で、 $f(\theta)$ が目標分布(既知)です。
詳細釣り合い条件は定常分布へ収束するための十分条件です。

定常分布へ収束するためには、 $f(\theta)$ を $t$ で微分して $=0$ とおいたマスター方程式

$\begin{eqnarray} \frac{{\rm d}f(\theta)}{{\rm d}t}=\sum_{\theta'}(-f(\theta'|\theta)f(\theta)+f(\theta|\theta')f(\theta'))=0\tag{2} \end{eqnarray}$

を満たせばいいのですが、これを解くのは難しいようです。

なので、 $\sum$ の中身が0になるようにします。これが詳細釣り合い条件です。
(詳細釣り合い条件の詳細については参考リンクをご覧ください。もしくは他の文献をご覧ください。)

詳細釣り合い条件のイメージを下図に示します。

f:id:olj611:20210516022524p:plain:w400

MH法説明のための準備

サンプリングする分布を目標分布といい、 $f(\theta)$ とします。
遷移核をいきなり見つけるのは、難しいので、
$f(|)$ の代わりに適当な遷移核 $q(|)$ を提案分布として利用します。
ただし、提案分布は乱数発生が容易なものである必要があります。

ところで、提案分布を適当に決めるため、詳細釣り合い条件が成立しません。

$\begin{eqnarray} q(\theta|\theta')f(\theta') > q(\theta'|\theta)f(\theta)\tag{3} \end{eqnarray}$

$(3)$ の等式が成り立つように $c,c'$ で補正します。

$\begin{eqnarray} &&cq(\theta|\theta')f(\theta') = cq(\theta'|\theta)f(\theta)\\ &&\Leftrightarrow\frac{c}{c'}q(\theta|\theta')f(\theta') = \frac{c'}{c'}q(\theta'|\theta)f(\theta)\\ &&\Leftrightarrow rq(\theta|\theta')f(\theta') = r'q(\theta'|\theta)f(\theta)\tag{4} \end{eqnarray}$

$(4)$ で $r=c / c',r'=c' / c'$ としました。
$r$ は $(3)$ より、 $r\in[0,1]$ です。
よって、 $(4)$ は、 $rq(\theta|\theta')f(\theta')$ と $r'q(\theta'|\theta)f(\theta)$ を遷移核として採用すれば、詳細釣り合い条件を満たす遷移が実現できることを示しています。

MH法では $r$ を確率とみなして、確率的補正を行うことにより、 $(4)$ が成り立つようにします。

MH法アルゴリズム

1: 提案分布 $q(|\theta^{(t)})$ を利用し、乱数 $a$ を発生させる。

2: 以下の命題を判定します。

$\begin{eqnarray} q(a|\theta^{(t)})f(\theta^{(t)}) > q(\theta^{(t)}|a)f(a)\tag{5} \end{eqnarray}$

・ $(5)$ が真
$(3)$ の不等号条件より、 $\theta=a,\theta'=\theta^{(t)}$ の状態と判定されます。
この場合、提案分布は $q(\theta|\theta')$ なのですから、 $(4)$ 式左辺の係数 $r$ を使って確率的補正を行う必要があります。

$\begin{eqnarray} r=\frac{q(\theta^{(t)}|a)f(a)}{q(a|\theta^{(t)})f(\theta^{(t)})}\tag{6} \end{eqnarray}$

を計算し、
確率 $r$ で $a$ を受容し、 $\theta^{(t+1)}=a$ とします。
確率 $1-r$ で $a$ を破棄し、 $\theta^{(t+1)}=\theta^{(t)}$ とします。

・ $(5)$ が偽
$(3)$ の不等号条件より、 $\theta'=a,\theta=\theta^{(t)}$ の状態と判定されます。
この場合、提案分布は $q(\theta'|\theta)$ なのですから、 $(4)$ 式右辺の係数 $r'$ を使って確率的補正を行う必要があります。
確率 $1(=r')$ で $a$ を受容し、 $\theta^{(t+1)}=a$ とします。
事実上補正しません。

3: $t=t+1$ として、1へ戻ります。

MH法を噛み砕いて説明

私の自身の言葉で噛み砕いてMH法を説明してみたいと思います。(MH法を理解している方は読み飛ばしてください。)

※説明のために $q(\theta'|\theta)f(\theta),q(\theta|\theta')f(\theta')$ を流れと呼ぶことにします。

$q(|\theta)$ から乱数 $\theta'$ を発生させます。
$\theta'$ を提案分布から提案された遷移先の候補と考えます。

以下のようにして採用するか決めます。

・ $q(\theta'|\theta)f(\theta) > q(\theta|\theta')f(\theta')$ の場合
　　　↓
$\theta$ から $\theta'$ への流れが詳細釣り合い条件より強い。
　　　↓
提案分布が詳細釣り合い条件を満たしていれば、 $\theta'$ に遷移していなかった可能性がある。
　　　↓
確率的に補正し、採用するか決める。

・ $q(\theta'|\theta)f(\theta) < q(\theta|\theta')f(\theta')$ の場合
　　　↓
$\theta'$ から $\theta$ への流れが詳細釣り合い条件より弱い。
　　　↓
提案分布が詳細釣り合い条件を満たしていないにもかかわらず、候補である。
　　　↓
提案分布を満たしていれば、確実に遷移している。
　　　↓
無条件に採用する。

MH法の $r$ の計算について

目標分布が事後分布であることを考えて、 $r$ を式変形していきます。

$\begin{eqnarray} r&=&\frac{q(\theta^{(t)}|\theta_a)f(\theta_a|x)}{q(\theta_a|\theta^{(t)})f(\theta^{(t)}|x)}\\ &=&\frac{q(\theta^{(t)}|\theta_a)\frac{f(x|\theta_a)f(\theta_a)}{f(x)}}{q(\theta_a|\theta^{(t)})\frac{f(x|\theta^{(t)})f(\theta^{(t)})}{f(x)}}\\ &=&\frac{q(\theta^{(t)}|\theta_a)f(x|\theta_a)f(\theta_a)}{q(\theta_a|\theta^{(t)})f(x|\theta^{(t)})f(\theta^{(t)})}\tag{7} \end{eqnarray}$