⛓️

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

2021/04/21に公開

機械学習

統計

tech

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

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

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

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 法は以下のようなアルゴリズムになる。

初期値 $X_1 = x_1$ を適当に定める
$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)$ のサンプルにより近づいていくことも分かる。

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

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

Metropolis Hastings 法

実装

References

Discussion