Metropolis Hastings 法の説明と Python による実装

MCMC（マルコフ連鎖モンテカルロ法）

ある分布 $\pi(x)$ からサンプルを得たいと思った時に、$\pi(x)$ が高次元分布の場合やそもそも解析的に表すのが難しい場合、直接的な手法で上手くサンプルを得ることは難しい。

そのようなケースの対処法として MCMC（マルコフ連鎖モンテカルロ法） がある。これは $\pi(x)$ を定常分布とするマルコフ連鎖をシミュレートするアルゴリズムである。このアルゴリズムを繰り返すことで、得られる値が $\pi(x)$ のサンプルに近いものになる。

マルコフ連鎖とは $P(X_i|X_1, ..., X_{i-1}) = P(X_i|X_{i-1})$ というように、ある状態が前の状態にのみ依存して決まる確率過程のこと。
定常分布 (stationary distribution) とは、確率過程によって次の状態に遷移しても変化しないような分布。もしマルコフ連鎖を繰り返してある分布に収束する場合（極限分布）、それは定常分布にもなることが知られている。

Metropolis Hastings 法

$\pi(x)$ を定常分布とするマルコフ連鎖をシミュレートする方法の一つとして、Metropolis Hastings 法 がある。
Metropolis Hastings 法においては、ユーザーが 提案分布 (proposal distribution) を指定してあげる必要がある。"提案" とは、サンプルが $x$ にあったときに、次のサンプル $y$ はどの位置にあるべきかを示すもので、その分布は $Q(y|x)$ と表記される。

単純な例として、$y$ は $x$ の位置の近くをランダムに選ぶとする。例えば正規分布に従う乱数を用いて

$$y | x \sim N(x, 1)$$

とおける。言い換えると、提案分布 $Q(y|x)$ は

$$Q(y|x) = \dfrac{1}{\sqrt{2 \pi}} \exp\left( - \dfrac{1}{2}(y-x)^2 \right)$$

という式になる。

この提案分布 $Q(y|x)$ を用いて、Metropolis Hastings 法 は以下のようなアルゴリズムになる。

1. 初期値 $X_1 = x_1$ を適当に定める
2. $t \geq 1$ に対して以下を十分な回数繰り返す
   • $Q(y|x_t)$ から $y$ という値をサンプルする（提案）
   • $A = \min \left(1, \dfrac{\pi(y)Q(x_t|y)}{\pi(x_t)Q(y|x_t)} \right)$ を計算する
   • 確率 $A$ で提案を受諾し、$x_{t+1} = y$ とする. 確率 $1-A$ で提案を拒否し $x_{t+1} = x_t$ に留まる

実装

例として、

$$\pi(x) = \begin{cases} \exp(-x) & (x \geq 0) \\ 0 & (x < 0) \end{cases}$$

という分布からサンプルを得たいとする。もちろん、このような単純な分布であれば MCMC をわざわざ使う必要はないが、あくまで簡単な例としてである。

提案分布 $Q(y|x)$ は上の例で挙げた、正規分布乱数による移動とする。正規分布乱数の場合、$Q(y|x)$ の式の形から $Q(x_t|y) = Q(y|x_t)$ が成り立つので、$A$ は

$$A = \min \left(1, \dfrac{\pi(y)}{\pi(x_t)}\right)$$

とよりシンプルに書ける。$A$ をこの式に置き換えた手法は Metropolis 法と呼ばれ、Metropolis Hastings 法の前身である（Hastings が後により一般的な $Q(y|x)$ に対して拡張したということである）。

この Metropolis 法を Python で実装してみると、以下のようになる。

import numpy as np

def pi(x):
    if x >= 0:
        return np.exp(-x)
    else:
        return 0
	
xs = []
xs.append(10) # 初期値を適当に 10 と定めた

for i in range(10000):
    y = np.random.normal(xs[i], 1) # 提案値
    a = min(1, pi(y) / pi(xs[i]))

    u = np.random.rand() # 0 ~ 1 の一様乱数を生成

    if u < a: # 受諾
        xs.append(y)
    else: # 拒否
        xs.append(xs[i])

さて、このように計算された xs は $\pi(x)$ のサンプルと考えられるので、このヒストグラムが $\pi(x)$ に近い形をしていると考えられる。これを確かめてみよう。

import matplotlib.pyplot as plt
%matplotlib inline

plt.hist(xs, bins=20)

確かに、右肩下がりの指数分布の形になっていることがわかる。

また、最初の100個のサンプルと次の100個のサンプルを比べてみる。

plt.hist(xs[0:100], bins=20)

plt.hist(xs[100:200], bins=20)

後者の方が、何となく右肩下がりがハッキリと確認できるように見える。このように試行回数を増やすことで、$\pi(x)$ のサンプルにより近づいていくことも分かる。

References

• マルコフ連鎖の定常分布と極限分布
• 第11章 サンプリング法
• MCMC チュートリアル 入門から多峰性分布の扱いとその応用まで
• The Metropolis Hastings Algorithm