0. 概要
マルコフ連鎖モンテカルロ法(MCMC)は確率分布のサンプリングを行うアルゴリズムの総称である。主に推定したい事象の分布（不確実性など： どれくらいの確度で推定できているのか）を知りたいときに用いる。今回はこのMCMC法の中でも最も有名なMH（メトロポリス・ヘイスティングス法）に入門する。
最終的には得られたデータ点$x$から、その$x$を発生させている確率密度関数（$\mu$）を推定する。
簡単に、「パラメータ（例えば、平均値や分散値）を正規分布等でサンプリング」→「そのパラメータを元にある関数（分布）を構成する」→「ある関数にデータ$x$を重ねて、その関数がデータ$x$を表現できているか評価する」→「表現できていそうなら、そのパラメータを採用して、最初に戻る」ということをしている。
以下を参考にした。
http://watanabe-www.math.dis.titech.ac.jp/~kohashi/document/bayes_51.pdf
https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%AB%E3%82%B3%E3%83%95%E9%80%A3%E9%8E%96%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95
https://www.igaku-shoin.co.jp/paper/archive/y2020/PA03354_02
https://ricrowl.hatenablog.com/entry/2020/08/22/121158
https://www.astr.tohoku.ac.jp/~chinone/Markov_Chain_Monte_Carlo/Markov_Chain_Monte_Carlo.pdf
https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%88%E3%83%AD%E3%83%9D%E3%83%AA%E3%82%B9%E3%83%BB%E3%83%98%E3%82%A4%E3%82%B9%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0%E3%82%B9%E6%B3%95
1. メトロポリスヘイスティングス法(M-H)
簡単な用語とアルゴリズム処理について説明する。
1.1. 用語と例題
まずは簡単な説明と例題を説明して勘所を掴む。
定常分布
定常分布（不変分布ともいう）$\pi$は、時間が経過しても確率過程の分布が変化しない分布のことである。
すなわち、常に同じ分布であることを指す（サイコロの確率等）。
次に離散時間マルコフ連鎖において、ある状態$i$からある状態$j$へ推移する確率を推移確率行列$P$という。
この時、
\pi P = \pi
を満たす分布$\pi$が存在すれば、$\pi$は定常分布となる。
詳細つり合い条件
次に詳細つり合い条件という考えがある。これは、A-Bという2つの状態が平衡状態(AとBの存在確率は変わらない)にある場合、A→Bの確率とB→Aの確率は等しくならないといけない定理である。
言い換えると、Aの箱にはx個の玉がないといけない、Bの箱にはy個の玉がないいけないという平衡状態の時、A→Bに1個の玉を移動したとする、その時、B→Aに1個の玉を移動しなければいけない。
この状態のことを詳細つり合い条件という。数式では以下になる。
\pi_A P(A \rightarrow B) = \pi_B P(B \rightarrow A)
ここで、$P(A \rightarrow B)$及び$P(B \rightarrow A)$のことを推移確率行列$P$という。
推移確率行列Pの求め方 (採択確率A)
次に、$P(A \rightarrow B)$及び$P(B \rightarrow A)$を以下の式によって求める。
A\left( \pi',\pi\right) = \min \left(1,\dfrac{P(\pi')}{P(\pi)}\dfrac{Q(\pi|\pi')}{Q(\pi'|\pi)}\right)
ここで、$A\left( \pi',\pi\right)=A\left( \pi,\pi'\right)$が一般性を失うことなく成り立つことが重要である。
例
定常分布$\pi$において晴れの確率$\dfrac{7}{10}$、雨の確率$\dfrac{3}{10}$とする。
次に推移確率行列$P$を晴れから雨を$P(A \rightarrow B)$及び雨から晴れを$P(B \rightarrow A)$とする。
以下の式から、
\pi P = \pi
詳細つり合い条件において次式が成り立つ
\begin{bmatrix}
1-P(A \rightarrow B) & P(B \rightarrow A) \\
P(A \rightarrow B) & 1-P(B \rightarrow A)
\end{bmatrix}
\begin{bmatrix}
\dfrac{7}{10} \\
\dfrac{3}{10}
\end{bmatrix}
=
\begin{bmatrix}
\dfrac{7}{10} \\
\dfrac{3}{10}
\end{bmatrix}
次に、次式から推移確率行列$P$の$P(A \rightarrow B)$及び$P(B \rightarrow A)$を求める。
A\left( \pi',\pi\right) = \min \left(1,\dfrac{P(\pi')}{P(\pi)}\dfrac{Q(\pi|\pi')}{Q(\pi'|\pi)}\right)
すなわち、
\min \left(1,\dfrac{7}{10} \div \dfrac{3}{10}\right) = 1 \\
\min \left(1,\dfrac{3}{10} \div \dfrac{7}{10}\right) = \dfrac{3}{7}
となり、$\pi P = \pi$は、
\begin{bmatrix}
1-\dfrac{3}{7} & 1 \ \\
\dfrac{3}{7} & 1-1
\end{bmatrix}
\begin{bmatrix}
\dfrac{7}{10} \\
\dfrac{3}{10}
\end{bmatrix}
=
\begin{bmatrix}
\dfrac{7}{10} \\
\dfrac{3}{10}
\end{bmatrix}
のように$\pi$が同じ結果になる。
1.2. アルゴリズムの処理
MHアルゴリズムの説明をする。
a) 初期化
$f(x)$を目標分布$P(x)$の確率密度関数に比例する関数とする。
まず初めに以下を決める。
- $x_{0}$: 初期値
- $Q(x|y)$: 任意の確率密度関数（提案分布という）
ここで、$Q(x|y)$は過去の状態$y$が与えられた時に、候補となる状態$x$を生成する関数である。
$Q$は対称性があり、$Q(x|y) = Q(y|x)$を満たす。典型的に$Q(x|y)$は、$y$を平均とした正規分布である。
b) 提案分布から遷移先の候補を見つける
確率分布$Q(\pi'|\pi)$から次の状態となる候補$\pi'$を生成する。
候補を見つけるため、採択率$A=\dfrac{f(\pi')}{f(\pi)}$を計算する。
この時、$f$は$P$の確率密度関数（推移確率行列）に比例している。
すなわち以下のような式になる。
A\left( \pi',\pi\right) = \min \left(1,\dfrac{P(\pi')}{P(\pi)}\dfrac{Q(\pi|\pi')}{Q(\pi'|\pi)}\right)
= \min \left(1,\dfrac{f(\pi')}{f(\pi)}\dfrac{Q(\pi|\pi')}{Q(\pi'|\pi)}\right)
因みに、$P$は事後確率$Pr(X|D)$である。事後確率は以下で求められる。
Pr(X|D) = \dfrac{Pr(D|X)Pr(x)}{\Sigma Pr(D|X)Pr(X)}
この時、$Pr(D|X)$は尤度、$Pr(X)$は事前確率である。
c) 遷移する
$A\geq1$の場合、$\pi'$を採択する。この場合$\pi \rightarrow \pi'$へと置き換わる。
ただし、そうでない場合は、$A$の確率で候補$\pi'$を採択する。
棄却された場合は、$\pi \rightarrow \pi$のままである。
例: 現在が晴れで、次の状態が「晴れの後、雨である確率」が100%以上であれば、雨とする。ただし、「晴れの後、雨である確率」が80%であれば、80%で雨の状態へとする。
更新処理
a) → b) → c)を繰り返す。
2. 実装
実際に数式をコードに落とし込んでいく。上記の数式を基本はそのまま使うが、プログラム化しにくい変数があるため置き換えている。適宜、読み替えて頂きたい。
最尤法や尤度関数との大きな違いは、分布の頂点（最もその事象が発生する確率）を知ることではなく、分布そのものを知ることである。
すなわち、データ点は得られているが、データ点元の分布が分からない時に使える。
今回は正規分布のデータ点$x$は得られているが、大元となる分布が不明のため、どのような分布となっているかをMCMCで推論したいと考える。ここで、$\sigma$は既知として、$\mu$をMCMCによって推論する。
全体の処理の流れは以下の通りである。
2.1. ②データ点の作成
正規分布からサンプル値$x$をランダムに取得する。
$mu=5.0$、$\sigma=3.0$とする。
# generate data x
mu_x = 5.0
sigma_x = 3.0
N = 200
x = np.random.normal(loc=mu_x, scale=sigma_x, size=N)
2.2. ③xの初期化
提案分布$Q(x'|x^{t})$から遷移先候補を取得する場合、$x^{t}$に従い、$x'$を生成する必要がある。
まずは、この$x^{t}$に初期値を与える。今回は$0$でよいだろう。
x_t=0
2.3. ④提案分布から遷移先を取得
$x^{t}$を中心とした正規分布から遷移先$x'$を取得する。分布はガウス分布である。
ここがMCMCの中でも重要で、遷移先がランダムなためモンテカルロなのである。
正規分布のため一定の確率で外れた場所をサンプリングする。
def random_sampling_q (_mu, _sigma=0.5):
return np.random.normal(_mu, _sigma)
ここでの出力は$x'$である。
2.4. ⑤遷移先を評価する
$A(x', x^{t})$に従い、$x'$を採択するか否かを決める。
A\left( x',x^{t}\right) = \min \left(1,\dfrac{P(x')}{P(x^{t})}\dfrac{Q(x^{t}|x')}{Q(x'|x^{t})}\right)
= \min \left(1,\dfrac{f(x')}{f(x^{t})}\dfrac{Q(x^{t}|x')}{Q(x'|x^{t})}\right)
ここで、
\dfrac{P(x')}{P(x^{t})} = \dfrac{f(x')}{f(x^{t})}
は、提案分布の遷移先である$x'$と、その1つ前のサンプル（データ点）$x^{t}$の尤度比である。
これにより、現在のサンプル点と、次の遷移先のサンプル点がどれくらいズレているものなのかを計算する。
次に、
\dfrac{Q(x^{t}|x')}{Q(x'|x^{t})}
は、2つの提案分布の比率である。対称の場合は$1$になる。
コードに起こすと以下のようになる。
def gpdf(_x, _mu=0, _sigma=1):
return (1/(np.sqrt(2*np.pi*_sigma**2)))*np.exp(-((_x-_mu)**2)/(2*_sigma**2))
def prior(_x, _mu=1.0, _sigma=10.0):
return gpdf(_x, _mu, _sigma)
def q(_x, _mu, _sigma=0.5):
return gpdf(_x, _mu, _sigma)
def f(_x, _mu, _sigma):
return np.prod(gpdf(_x, _mu, _sigma)) * prior(_mu)
a = (f(x, x_new, sigma_x) * q(x_t, x_new)) / (f(x, x_t, sigma_x) * q(x_new, x_t))
gpdf()は確率密度関数、prior()は事前分布である。
なお、提案分布の比率を求めると$1$になるはずである。
print(q(x_t, x_new) / q(x_new, x_t))
# 1
また、事前分布は正規分布を仮定している。
\frac{1}{\sqrt{2\pi\sigma^2}}\exp \left( -\frac{(x-\mu)^2}{2\sigma^2} \right)
2.5. ⑥遷移先候補の採択
$A(x', x^{t})$の結果に基づき、$x'$を採択するか否かが決まる。
採択のルールは以下である。
$(a\geq1)$の場合は採択となり、
x^{t} = x'
$(a<1)$の場合は、$a$の確率で採択となり、
x^{t} = x'
それ以外は、
x' = x^{t}
となる。
コードは以下である。
min_a = min(1, a)
accept_flag = False
if min_a > np.random.uniform():
accept_flag = True
# return accept_flag, x_new
後は繰り返し処理させていくと、以下のような結果を得られる。
赤が真値、分布がMCMCで推定した事後分布のサンプル結果である。
最終的なコード
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
from scipy.stats import norm
SIGMA_Q = 0.5
PRIOR_S = 10
PRIOR_M = 0
# gauss
def gpdf(_x, _mu=0, _sigma=1):
return (1/(np.sqrt(2*np.pi*_sigma**2)))*np.exp(-((_x-_mu)**2)/(2*_sigma**2))
# prior
def prior(_x, _mu=PRIOR_M, _sigma=PRIOR_S):
return gpdf(_x, _mu, _sigma)
# q
def q(_x, _mu, _sigma=SIGMA_Q):
return gpdf(_x, _mu, _sigma)
# get next point
def random_sampling_q (_mu, _sigma=SIGMA_Q):
return np.random.normal(_mu, _sigma)
# likelihood
def f(_x, _mu, _sigma):
return np.prod(gpdf(_x, _mu, _sigma)) * prior(_mu)
# mh 1step
def mh(_x, _x_t, _q, _random_sampling_q, _f, _sigma_x):
x_new = _random_sampling_q(_x_t) # get next point
a = (_f(_x, x_new, _sigma_x) * _q(_x_t, x_new)) / (_f(_x, _x_t, _sigma_x) * _q(x_new, _x_t)) # evaluation of x_new
min_a = min(1, a)
accept_flag = False
if min_a > np.random.uniform():
accept_flag = True
return accept_flag, x_new
if __name__ == "__main__":
### generate data
mu_x = 5.0
sigma_x = 3.0
N = 200
x = np.random.normal(loc=mu_x, scale=sigma_x, size=N)
print('Data x AVE.:', np.average(x))
print('Data x STD.:', np.std(x))
plt.hist(x, bins=50)
plt.savefig("input.png")
### prior
prior_x = np.linspace(-30, 30, 101)
prior_pdf = prior(_x=prior_x, _mu=1.0, _sigma=10)
plt.plot(prior_x, prior_pdf)
plt.savefig("prior.png")
### start sampling
sample_num = 2000
warmup = 500
x_t = 0.0
samples = []
while(len(samples) < sample_num + warmup):
accept_flag, x_new = mh(x, x_t, q, random_sampling_q, f, sigma_x)
if accept_flag:
x_t = x_new
samples.append(x_t) # x_t means mu of data x
mcmc_mu = samples[warmup:]
### visualize mcmc data
plt.gca().clear()
plt.hist(samples, bins=50)
plt.savefig("output.png")
print('M.H. (posterior distribution) AVE.:', np.average(mcmc_mu))
print('M.H. (posterior distribution) STD.:', np.std(mcmc_mu))
x_mean = np.average(x)
eq_mu_N = (sigma_x**2*PRIOR_M + N*PRIOR_S**2*x_mean)/(sigma_x**2 + N*PRIOR_S**2)
eq_sigma_N = np.sqrt((PRIOR_S**2*sigma_x**2)/(sigma_x**2 + N*PRIOR_S**2))
print('GroundTruth AVE.:', eq_mu_N)
print('GroundTruth STD.:', eq_sigma_N)
plt.gca().clear()
plt.hist(mcmc_mu, bins=50, density=True, alpha=0.5)
plt_arr = np.linspace(min(mcmc_mu), max(mcmc_mu), 1000)
plt.plot(plt_arr, norm.pdf(plt_arr, eq_mu_N, eq_sigma_N), c='r')
plt.savefig("comparison.png")
print("done")