任意の確率分布に従う標本生成：Markov連鎖　Metropolis Hasting法

棄却法では確率密度関数$${f({\bf x})}$$の上限を与える必要があったが、このMetropolis Hasting法では$${f({\bf x})}$$の値域に上限がない場合にも適用できる。
マルコフ連鎖過程は、$${{\bf x}_{t}}$$の状態が$${{\bf x}_{t-1}}$$で決まる過程であり、これを用いることで、定常状態が$${f(x)}$$の標本を得ることができるのである。
$${t}$において、$${t-1}$$での標本$${{\bf x}_{t-1}}$$に依存する確率密度$${g({\bf x}|{\bf x}_{t-1})}$$に従い、$${{\bf x}_t}$$を発生させる。
ここで、$${g({\bf x}|{\bf x}')=g({\bf x}'|{\bf x})}$$の平衡条件を与えることで、最終的に平衡状態$${f(x)}$$が成立する。
$${g({\bf x}|{\bf x}')}$$には、平衡条件に従い、かつ$${{\bf x}}$$は与えられている$${{\bf x}'}$$の近傍を分散$${\sigma^2}$$の正規分布に従って取るとして、
$${g({\bf x}_t|{\bf x}_{t-1})=\displaystyle{\frac{1}{(2\pi\sigma^2)^{d/2}} \exp\Big( -\frac{({\bf x}_t-{\bf x}_{t-1})^T({\bf x}_t-{\bf x}_{t-1})}{2\sigma^2} \Big)}}$$
と置く。ここで、$${d}$$は$${{\bf x}_t}$$の次元数である。
これで与えられた$${{\bf x}_t}$$が、$${(0,1]}$$の一様分布に従う乱数$${u}$$に対し、以下の条件
$${\displaystyle{u<\frac{f({\bf x}_t)g({\bf x}_{t-1}|{\bf x}_t)}{f({\bf x}_{t-1})g({\bf x}_t|{\bf x}_{t-1})}}}$$
である時、次の状態として更新される。
ここで、平衡条件$${g({\bf x}|{\bf x}')=g({\bf x}'|{\bf x})}$$を考慮に入れれば、条件式は、
$${\displaystyle{u<\frac{f({\bf x}_t)}{f({\bf x}_{t-1})}}}$$
と単純になる。

この方法を標準コーシー分布に従う標本生成に適用したのが以下のコードとなる。

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)
randnum=1_000_000

def f(x):
    return 1/(np.pi *(x**2 + 1))

x = []
x0=1
x.append(x0)
sig=1.0

for i in range(randnum):
    y = np.random.normal(x[i], sig) 
    a = min(1, f(y) / f(x[i]))
    u = rand.rand() 
    if u < a: 
        x.append(y)
    else: 
        x.append(x[i])
x0=[i for i in x if (i > -5) & (i<5)]
plt.hist(x0, density=True, bins='auto',label='Metropolis Hasting')
cau=rand.standard_cauchy(size=randnum)
plt.hist(cau[(cau>-5) & (cau<5)], density=True,bins='auto',histtype="step",label='standard caucy')
plt.ylabel('Probability')
plt.legend()
plt.xlabel('Data');