【ベイズ統計学再入門④】マルコフ連鎖モンテカルロ法（MCMC）

復習がてら、ベイズ統計の解説とPythonでの実装を行います。
第4回は「マルコフ連鎖モンテカルロ法（MCMC）」についてです。

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

マルコフ連鎖モンテカルロ法（Markov chain Monte Carlo methods、通称MCMC）とは、マルコフ連鎖を生成することで、任意の分布（正規分布でもなんでもOK）を定常分布とするような乱数を生成するモンテカルロ法のことを言います。

マルコフ連鎖とは確率過程（ステップ数などによってパラメータづけられた確率変数の族）の一種で、任意のデータ$x,x_1,\ldots, x_n$に対して

P(\theta_{n+1}=x|\theta_{n}=x_{n},\theta_{n-1}=x_{n-1},\cdots,\theta_1=x_1)=P(\theta_{n+1}=x|\theta_{n}=x_n)

を満たすようなものを言います。もう少し例えて表現すると、上記の左辺は

「$n$日分のデータ$x_1,\ldots, x_n$を観測したときに、次の日に$x$が起こる確率」

を表します（事後分布）。一方右辺は

「今日$x_n$を観測したときに、次の日に$x$が起こる確率」

と思うことができます。マルコフ連鎖の条件は、この2つの確率が等しいということを意味しているので、標語的に言うとマルコフ連鎖は「明日の確率は今日の確率で決まる」ような確率過程のことと言えるでしょう。漸化式のようなものと思ってもいいかもしれません。

マルコフ連鎖の例

マルコフ連鎖の例として、次のようなものを考えようと思います。

Aさんはじゃんけんの手の出し方に次のような癖があります：
①今グーを出しているとき
　次にグーを出す確率：$0.4$
　次にチョキを出す確率：$0.3$
　次にパーを出す確率：$0.3$
②今チョキを出しているとき
　次にグーを出す確率：$0.3$
　次にチョキを出す確率：$0.5$
　次にパーを出す確率：$0.2$
③今パーを出しているとき
　次にグーを出す確率：$0.6$
　次にチョキを出す確率：$0.1$
　次にパーを出す確率：$0.3$
このとき、一人でじゃんけんをしばらく続けたとき、Aさんはどの手を出している可能性が高いか？

問題の設定から、次に出す手の確率は今出している手で決まるので、これはマルコフ連鎖と言えます。この問題をプログラミングで解いてみます（厳密ではない）。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

np.random.seed(123)
p_trans = np.zeros([3,3])
# 0:グー, 1:チョキ, 2:パー
p_trans[0,0] = 0.4
p_trans[0,1] = 0.3
p_trans[0,2] = 0.3
p_trans[1,0] = 0.3
p_trans[1,1] = 0.5
p_trans[1,2] = 0.2
p_trans[2,0] = 0.6
p_trans[2,1] = 0.1
p_trans[2,2] = 0.3

# モンテカルロステップ
NMCS = 1000

# 初期値（最初はグー）
c_state = 0
c_arr = [c_state]
for i in range(NMCS):
    #0,1,2の中からランダムに1つだけチョイスする
    #「p=」とすることで指定の確率で数を選べる
    current = np.random.choice(3,1,p=p_trans[c_state, :])
    #current は array([選ばれた数])という形なので、
    #[0]をつけて選ばれた数を取り出す
    c_state = current[0]
    #モンテカルロステップを追加する
    c_arr.append(c_state)

# Aさんが出した手のデータ
df = pd.DataFrame(c_arr)

# ヒストグラムの作成
plt.hist([df[0][:10], df[0][:50], df[0][:100], df[0][:400], df[0][:700], df[0][:1000]],
        density=True, label=["10","50","100","400","700","1000"])
plt.xlabel("guu/choki/paa")
plt.ylabel("frequancy")
plt.legend()

じゃんけんの回数を増やしていくと、グー・チョキ・パーそれぞれの出る割合が収束していく様子が観察できました（この確率分布を定常分布といいます）。よって、じゃんけんの回数を増やしていくと、Aさんはグーを出している確率が高いようです。

詳細つり合い

先ほどの例では定常分布が存在することがわかりましたが、必ずしも定常分布が存在するかどうかはわかりません。そこで。定常分布が存在するための"十分条件"として「詳細つり合い」という条件があります。

P(\theta'|\theta)P(\theta)=P(\theta|\theta')P(\theta')

詳細つり合いの$P(\theta'|\theta)$の部分は遷移核といい、これの求め方は様々な方法が考えられていますが、この記事では「ランダムウォークM-Hアルゴリズム」を紹介します。

ランダムウォークM-Hアルゴリズムは次のようなアルゴリズムになっています。

• 初期値$\theta$を適当に決める
• ランダムウォークで新しい$\theta_{\text{new}}$を決める:

\theta_{\text{new}} = \theta + \varepsilon\text{Normal}(0, 1)

• $f(\theta_{\text{new}}) > f(\theta)$かどうかを判定する
• 真のとき:$\theta_{\text{new}}$を受け入れる
• 偽のとき:$r=\dfrac{f(\theta_{\text{new}})}{f(\theta)}$の確率で受け入れる
• この一連のステップを繰り返す

それでは試しに今回はベータ分布をサンプリングしてみようと思います。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
import pandas as pd
%matplotlib inline

# 今回はこのベータ分布からサンプリング
a, b = 1.5, 2.0
x = np.linspace(beta.ppf(0.001,a,b), beta.ppf(0.999,a,b), 100)
plt.plot(x, beta.pdf(x, a, b))

theta = 0.8
NMCS = 20000
# ランダムウォークの幅の調整は一般に難しい
ep = 0.5
theta_mcs = [theta]
for i in range(NMCS):
    theta_new = theta + ep * np.random.randn()
    if beta.pdf(theta_new, a, b) > beta.pdf(theta, a, b):
        theta = theta_new
    else:
        r = beta.pdf(theta_new, a, b) / beta.pdf(theta, a, b)
        if np.random.rand() < r:
            theta = theta_new
    theta_mcs.append(theta)

# シータの値がまんべんなく分布を動いていれば収束していると言える
plt.plot(df[0])
plt.xlabel("MCS")
plt.ylabel("$\Theta$")

今回は早い段階で安定しているので、1000より後のデータでヒストグラムを作ってみます。

plt.hist(df[0][1000:], density=True, bins=40)
x = np.linspace(beta.ppf(0.001,a,b), beta.ppf(0.999,a,b), 100)
plt.plot(x, beta.pdf(x, a, b))

いい感じにサンプリングできました。