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

3 min read読了の目安(約3100字

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) と表記される。

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

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