Edited at

数値例とコードから学ぶMH(メトロポリス・ヘイスティングス)法


はじめに

サンプリングについて学ぶためにMCMC法やMH法について調べ始めたものの正直よくわからん.

なので実際にプログラムを組みながら理解したい!という思いで記事を書きました.

高校の数学の知識があれば理解できます. python3で実装してます.


そもそもMH法とは

MH法(メトロポリス・ヘイスティングス法)とは, ある目標分布からのサンプリングが難しい場合に, 別の簡単な提案分布を用いてサンプリングを行う方法である(wikipedia).

とにかくサンプリングが難しい標本空間の中から頑張ってサンプリングする方法があるよということらしい.


数値例(分割表)

分割表-Wikipediaから

右利き
左利き

男性
43
9
52

女性
44
4
48


87
13
100


調べたいこと

男性と女性で左利きの割合に差があるか?

男性, 女性の左利きの割合を同じ$p$とすると, 男性の左利きが52人中9人となる確率は男性の左利きの人数を$M_L$とすると

P(M_L=9)={}_{52}\mathrm{C} _9 \, p^9(1-p)^{43}

同様に女性の左利きの数を$F_L$とすると女性の左利きが48人中4人となる確率は

P(F_L=4)={}_{48}\mathrm{C} _4 \, p^4(1-p)^{44}

よって, この表が得られる確率は

\begin{align}

P(M_L=9,F_L=4)&={}_{52}\mathrm{C} _9 \, p^9(1-p)^{43}{}_{48}\mathrm{C} _4 \, p^4(1-p)^{44}\\
&={}_{52}\mathrm{C} _9\,{}_{48}\mathrm{C} _4 \, p^{13}(1-p)^{87}
\end{align}

左利きの数が一定($M_L+F_L=13$)と仮定すると, $M_L=x\,$が得られる確率は

\begin{align}

[\pi]_x&=P(M_L=x|M_L+F_L=13)\\
&=\frac{{}_{52}\mathrm{C} _x\,{}_{48}\mathrm{C} _{13-x} \, p^{13}(1-p)^{87} }{\sum_{i=1}^{13}{}_{52}\mathrm{C} _i\,{}_{48}\mathrm{C} _{13-i} \, p^{13}(1-p)^{87}}\\
&=\frac{{}_{52}\mathrm{C} _x\,{}_{48}\mathrm{C} _{13-x}}{\sum_{i=1}^{13}{}_{52}\mathrm{C} _i\,{}_{48}\mathrm{C} _{13-i} }
\end{align}

この分布は分母に複雑な計算が含まれている. そこでMH法を使って$[\pi]_x$に従う乱数を生成し, ヒストグラムにすることで分布を近似する.


MH法

まずアルゴリズムについて説明する.

1. 入力として$x$を代入する.

2. 確率$Q(x^*|x)$で$x^*$を出力する.

3. 採択率$\alpha(x^*|x)$で$x^*$を出力, $1-\alpha(x^*|x)$で$x$を出力する.

$Q(x^*|x)$は提案分布と呼ばれている. 確率$Q(x^*|x)$は任意のもので構わないが, $Q(x^*|x)=Q(x|x^*)$である必要がある. 一般的にはランダムウォークが用いられている.

採択率$\alpha(x^*|x)$は

\alpha(x^*|x)=\min\left(1,\frac{[\pi]_xQ(x^*|x)}{[\pi]_{x^*}Q(x|x^*)}\right)

で表される.ここで$Q(x^*|x)=Q(x|x^*)$のため

\alpha(x^*|x)=\min\left(1,\frac{[\pi]_{x^*}}{[\pi]_{x}}\right)

となる.

結局のところ, 何が出力されるのかというと...


  • $\frac{[\pi]_{x^*}}{[\pi]_{x}}>1$のとき, $x^*$を出力


  • $\frac{[\pi]_{x^*}}{[\pi]_{x}}<1$のとき, $\frac{[\pi]_{x^*}}{[\pi]_x}$の確率で$x^*$, $1-\frac{[\pi]_{x^*}}{[\pi]_x}$の確率で$x$を出力


となる.

つまりは, $x^*$と$x$の起こりうる確率を調べて,


  • $x^*$の方が起こりそうなら$x^*$に移動


  • $x$の方が起こりそうならある確率で立ち止まる


ということになる. その確率というものがサンプリングしたい分布$[\pi]$になっているから, サンプル列$x$も$[\pi]$に従うという仕掛けになっているらしい. こうして見るとなかなかうまくできている(気がする).


数値例-解き方

ではこのアルゴリズムを数値例にあてはめてみよう.

まず提案分布$Q(x^*|x)$を実装する. 要はランダムウォーク. 入力$x$に対して$\frac{1}{2}$の確率で$x+1$または$x-1$を返す関数を作ればよい. ただし, 定義域の端点の場合はどちらかが$x$になる(端にいるなら外には出られない).


Q(x)

def Q(x):

if(x == 0):
if(np.random.uniform(0, 1) < 0.5):
y = 1
else:
y = x
elif(x == total):
if(np.random.uniform(0, 1) < 0.5):
y = total_lefty - 1
else:
y = x
else:
if(np.random.uniform(0, 1) < 0.5):
y = x - 1
else:
y = x + 1
return y

次に採択確率$\alpha(x^*|x)$を実装する.

採択率$\alpha(x^*|x)$は

\alpha(x^*|x)=\min\left(1,\frac{[\pi]_{x^*}}{[\pi]_{x}}\right)

で表される. また, $[\pi]_x$は

\begin{align}

[\pi]_x&=\frac{{}_{52}\mathrm{C} _x\,{}_{48}\mathrm{C} _{13-x}}{\sum_{i=1}^{13}{}_{52}\mathrm{C} _i\,{}_{48}\mathrm{C} _{13-i} }
\end{align}

であった. この式の分母は定数であるため,

\begin{align}

\frac{[\pi]_{x^*}}{[\pi]_x}&=\frac{{}_{52}\mathrm{C} _{x^*}\,{}_{48}\mathrm{C} _{13-{x^*}}}{{}_{52}\mathrm{C} _x\,{}_{48}\mathrm{C} _{13-x}}
\end{align}

となり厄介だった分母が消えるというのがミソである.

しかも$x^*$は結局のところ, $x-1, x, x+1$のどれかなので, もっと簡単にすることができる. 気合で計算すると

  \frac{[\pi]_{x^*}}{[\pi]_x} = \begin{cases}

\frac{(13-x)(52-x)}{(x+1)(36+x)} & (x^*=x+1) \\
1 & (x^*=x)\\
\frac{x(35+x)}{(14-x)(53-x)} &(x^*=x-1)
\end{cases}

となり, かなり簡単に計算できる. 実装すると,


alpha(x,y)

x = current

y = candidate
if(y == x + 1):
a = ((total_lefty - x) * (male - x)) / (y * (female - total_lefty + y))
elif(y == x - 1):
a = (x * (female - total_lefty + x)) / ((total_lefty - y) * (male - y))
else:
a = 1

こんな感じかな?

これらをもとにMH法を指定された回数$N$回繰り返す関数を実装する.


metropolis(N)

def metropolis(N):

current = 0
sample = []
sample.append(current)
accept_ratio = 0

for i in range(N):
candidate = Q(current)
x = current
y = candidate
if(y == x + 1):
a = ((total_lefty - x) * (male - x)) / (y * (female - total_lefty + y))
elif(y == x - 1):
a = (x * (female - total_lefty + x)) / ((total_lefty - y) * (male - y))
else:
a = 1
#確率aによって新しいcandidateを出力
if(a > np.random.uniform(0, 1)):
# Update state
current = candidate
sample.append(current)
accept_ratio+=1


accept_ratioで受け入れ比率を求める. なくても問題はない.

これらをまとめると, こんな感じ


metropolis_hastings.py

import copy

import numpy as np
import matplotlib.pyplot as plt
import math

total_lefty = 13
male = 52
female = 48

#真の値の計算
def combinations_count(n, r):
return math.factorial(n) / (math.factorial(n - r) * math.factorial(r))

def pi(x):
sum = 0
for i in range(total_lefty + 1):
sum+=combinations_count(male, i) * combinations_count(female, total_lefty - i)
return combinations_count(male, x) * combinations_count(female, total_lefty - x) / sum

# Q(x,y) : ranfom walk
# ランダムにxから1ずれるようなyを出力
def Q(x):
if(x == 0):
if(np.random.uniform(0, 1) < 0.5):
y = 1
else:
y = x
elif(x == total_lefty):
if(np.random.uniform(0, 1) < 0.5):
y = total_lefty - 1
else:
y = x
else:
if(np.random.uniform(0, 1) < 0.5):
y = x - 1
else:
y = x + 1
return y

def metropolis(N):
current = 0
sample = []
sample.append(current)
accept_ratio = 0

for i in range(N):
candidate = Q(current)
x = current
y = candidate
if(y == x + 1):
a = ((total_lefty - x) * (male - x)) / (y * (female - total_lefty + y))
elif(y == x - 1):
a = (x * (female - total_lefty + x)) / ((total_lefty - y) * (male - y))
else:
a = 1
#確率aによって新しいcandidateを出力
if(a > np.random.uniform(0, 1)):
# Update state
current = candidate
sample.append(current)
accept_ratio+=1

print('Accept ratio:', float(accept_ratio / N))
return np.array(sample)

def main():
#MH法
N = 100000
sample = metropolis(N)
target_dist = plt.hist(sample, bins=total_lefty + 1,range=(0,total_lefty),normed=True)
plt.title('x(MH method)')
plt.show()
#真の値を表示
ax_x = range(0,total_lefty + 1)
answer = []
for i in ax_x:
answer.append(pi(i))
plt.title('x(target distribution)')
plt.bar(ax_x, answer)
plt.show()

if __name__ == '__main__':
main()


ついでに真の分布との比較も行った.出てきた表は以下のようになった.

Figure_1.png

Figure_1-1.png

...合っているような気がする.

male,female,total_leftyをそれぞれ10倍にしてみる.

Figure_2.png

Figure_2-2.png

ほぼ一致した.


終わりに

MH法はアルゴリズムだけ聞かされるとややこしく感じるけど, 実際にやっていることや考え方は単純で理にかなっていると思います.

ちなみに$x$(男性の左利き)は130人中60~70人をとる確率が高いので多分男女で差はないです(知ってた).