概要
前回はMH法をまとめた。
今回はHMCで引っかかった点をまとめる。
MH法の問題点
MH法では遷移確率密度関数$f(a|b)$を解答者が任意に定めた。結果、設定の方法によって受容率が低く、初期値の影響を強く受けた。
受容率を高めるにはどうすればいいだろうか
定常条件が成立するように補正係数$r$をつけて
r f(a | b)f(b) = f(b|a)f(a)
を
(r f(a | b)) *f(b) = f(b|a)f(a)
として
f^\prime (a|b) = r f(a | b)
と考えるのがMH法だった。
rを受容率と考えることで定常条件を満たそうとした。
もっと原理的に受容率が高く、定常条件を満たすように$b$から$a$に遷移させる方法はないだろう。
HMC:ハミルトニアンモンテカルロ法 p.112
事後分布$f(\theta|x)$に標準正規分布に従う$p$を追加する。
$p,\theta$は互いに独立であるため、以下が言える
f(\theta,p|x) = f(\theta|x) f(p) \propto f(\theta|x) exp(-\frac{1}{2} p^2)
といえる。log,expで変形すると
f(\theta,p|x) \propto exp(log(f(\theta|x)) -\frac{1}{2} p^2)
ここで、
h(\theta) = -log(f(\theta|x))
と置きます。
f(\theta,p|x) \propto exp(-h(\theta) -\frac{1}{2} p^2)
となり、exp内をハミルトニアン$-H(\theta,p)$と置きます。
H(\theta,p) =h(\theta) + \frac{1}{2} p^2
とします。
アイデア
ここから話が飛躍していきます。
２つの条件（運動方程式と同じ、新たに時刻$t$が導入されていますが、マルコフ連鎖の遷移するのに時間がかかると見てください。）
p = \frac{d\theta}{dt}
-\frac{dh(\theta)}{d\theta} = \frac{dp}{dt}
の成立と
\frac{H(\theta,p)}{dt} = 0
は同値です。
さらに、前述の式
f(\theta,p|x) \propto exp(-h(\theta) -\frac{1}{2} p^2)
をもとにすれば、
\frac{df(\theta,p|x)}{dt} = 0
といえます。マルコフ連鎖では、遷移前に観測された母数$\theta$を元に遷移していきます。
マルコフ連鎖の定常条件は
\forall \theta',p'に対して
\frac{df(\theta,p|\theta',p',x)}{dt} = 0
の成立がマルコフ連鎖の定常条件になります。よって、まだ定常を言うには不十分です。
上の２つの条件（またはハミルトニアンの保存）はもう一つの性質を実現します。
\forall \theta,p,\theta',p'に対して
f(\theta,p|\theta ',p',x) = f(\theta ',p'|\theta,p,x)
ある母数$\theta$から遷移関数（上の２つの条件）に従って遷移した場合の分布関数の値（遷移確率）が逆方向の場合と一致するという点です。
(A点からB点に遷移する確率) = (B点からA点に遷移する確率)
であるということです。話がそれるので一度受け入れて証明は後にしましょう。
少しこんがらがっていると思います
まず、マルコフ連鎖で遷移前の分布$f(\theta|x)$があり、遷移則（前述２つの方程式）に従って遷移した結果の分布が$f(\theta'|\theta,x)$です。
\frac{df(\theta,p|x)}{dt} = 0
であったため
\forall \theta,p,\theta',p'に対して
f(\theta,p|x) = f(\theta',p'|x)
となります。
よって、
\forall \theta,p,\theta',p'に対して
f(\theta ',p'|\theta,p,x) f(\theta,p|x) = f(\theta,p|\theta ',p',x) f(\theta',p'|x)
となり、マルコフ連鎖の定常条件が満たされました。
マルコフ連鎖の定常条件
MH法で出てきました。
\forall a, \forall b に対して
f(a | b)f(b) = f(b|a)f(a)
が成立していると仮定すると、マルコフ連鎖の遷移によって発生する事後分布と事前分布が一致します。
ハミルトニアンの保存と可逆性についてp.114
ハミルトニアンの保存＝運動方程式
ma = F,
v = \frac{dx}{dt}
の成立。
終点での速度を正負逆にすれば、逆の経路をたどる。つまり、
f(a,b|a',b') = f(a',b'|a,b)
となる。
Pythonで実装
ベータ分布を実際に求めてみました。
受容率は99.85%でMH法より計算効率がよく、初期値の影響も少ないです。
import math
import numpy as np
import random
import matplotlib.pyplot as plt
import scipy.stats as st
def F(theta):
# if theta != 0:
# hoge = -13 + 10 / theta
# else:
# hoge = -100
if theta > 0:
hoge = -13 + 10 / theta
else:
hoge = 1000
return hoge
def h(theta):
if theta>0:
return 13*theta - 10 * math.log(theta)
else:
return 100
def H(theta,p):
return h(theta) + p*p/2
theta = 1
dt = 0.01
L=10
N = round(L/dt)
cul_rep_num = 0
rnd = np.empty([0]) #1次元から配列生成:保存用
theta_N = 2*10**4
for j in range(theta_N):
p = random.gauss(0,1)
theta_start = theta
p_start = p
for i in range(N):
p_05 = p + F(theta)*dt*0.5
theta_1 = theta + p_05*dt
p_1 = p_05 + F(theta_1)*dt*0.5
p = p_1
theta = theta_1
# 受容するか判定
r = math.exp( H(theta_start,p_start) - H(theta,p))
hoge = random.uniform(0,1)
if hoge < r:
rnd = np.append(rnd,theta)
cul_rep_num = cul_rep_num + 1
else:
theta = theta_start
p = p_start
if j%10**3 == 0:
print(j)
# ヒストグラム
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.hist(rnd, bins=100, histtype='barstacked', ec='black',density=True)
xx = np.linspace(0, 2.5,501)
ax.plot(xx, st.gamma(11, 0, 1/13.).pdf(xx))
plt.show()
#受容率
print(cul_rep_num)
print(theta_N)
print("---- accept num -----")
print(cul_rep_num/(theta_N))
最後に
書くの面倒になって糞記事になってしまった。