事後分布の期待値を推定値とする
$$ \hat{\theta_{eap}} = E[\theta|D]\\ = \int d\theta \cdot \theta \cdot p(\theta|D) \\ = \int d\theta \cdot \theta \cdot \frac{p(D|\theta) p(\theta)}{p(D)} $$事後分布の最大値を期待値の推定値とする
$$ \hat{\theta_{map}} = max_\theta p(\theta|D) $$事後分布の中央値を期待値の推定値とする -> 分布関数 $F$が0.5になるメディアン値を推定値とする.
$$ F(\hat{\theta_{med}} | D) = \int^{\theta_{med}} d\theta \cdot f(\theta | D ) = \frac{1}{2} $$幅をもたせてパラメータの推定をおこなう方法で,確率分布の極端な両端の$\alpha$%をのぞいた$(1 - \alpha)$%の部分を推定幅として用いる.伝統的な統計学では,$\alpha = 5$%を使うことが多い.
なぜモンテカルロサンプリングが必要なのか.それは,事後分布のランダムサンプリングが可能ならば,簡単にモデルのパラメータを推定できるから,である.
次の事後期待値を考える.
$$ \int d\theta \cdot f(\theta) p(\theta) $$結局,以下のように,事後分布のランダムサンプリングが可能ならば,簡単にモデルのパラメータを推定できることがわかった.
$$ \int d\theta \cdot f(\theta) p(\theta) = \frac{1}{N} \sum_{\theta_j \sim p(\theta_j)} f(\theta_j) $$MCMCとは,マルコフチェーンを利用して,欲しい確率分布(目標分布)の乱数を生成するアルゴリズムのことをさす.
マルコフチェーンとは,以下のような,一期前の変数にだけ依存するような条件付き確率で表現できる確率過程をいいます.
$$ p(X^{(t)}|X^{(t-1)}, X^{(t-2)}, \dots, X^{(1)}) = p(X^{(t)}|X^{(t-1)}) $$マルコフチェーンを規定する条件付き確率事態を遷移確率の行列を遷移核といいいます.
マルコフチェーンを利用する理由は,目標分布である定常状態に収束することが確率論的に保証されているからです.収束するまえの状態にあるチェーンは「バーン・イン, burn-in」といって,通常バーン・インにあるチェーンは取り除いて利用する.
MCMC法のアルゴリズムには以下のものがある.
- メトロポリス・ヘイスティング法
- ランダムメトロポリス・ヘイスティング法
- ギブスサンプリング法
- スライスサンプリング
- ハミルトニアンモンテカルロ法
- NUTS(No-U-Turn Sampler)
Initialize $x^{(1)}$
For i=2 to N
$x^{(i+1)}=\xs$
注意:もし遷移確率が以下のように対称なら,
$$q(\xs|x^{(i)}) = q(x^{(i)}|\xs)$$受け入れ確率を以下のようにすること.
$$ A(x,\xs) = \min\left[1, \frac{p(\xs)}{p(x) } \right] $$%matplotlib inline
import numpy as np
import matplotlib.pylab as plt
import pandas as pd
pd.set_option('display.width', 500)
pd.set_option('display.max_columns', 100)
import seaborn as sns
## FUNCTIONS
# target distribution p(x)
p = lambda x: 6*x*(1-x)
# number of samples
n = 100000
sig =0.10
#intitialize the sampling. Start somewhere from 0..1
x0 = np.random.uniform()
x_prev = x0
x=[]
k=1
i=0
while i<n:
x_star = np.random.normal(x_prev, sig)
while (x_star <0) | (x_star > 1): # MAKE SURE WE STAY WITHIN BOUNDS
x_star = np.random.normal(x_prev, sig)
P_star = 6*x_star*(1-x_star) #p(x_star);
P_prev = 6*x_prev*(1-x_prev) #p(x_prev);
U = np.random.uniform()
A = P_star/P_prev
if U < A:
x.append(x_star)
i = i + 1
x_prev = x_star
else :
x.append(x_prev)
x_prev = x[i]
i = i + 1
k=k+1
e,q,h=plt.hist(x,10, alpha=0.4, label=u'MCMC distribution')
xx= np.linspace(0,1,100)
plt.plot(xx, 0.67*np.max(e)*p(xx), 'r', label=u'True dsitribution')
plt.legend()
<matplotlib.legend.Legend at 0x7f1554274860>
from IPython.display import YouTubeVideo
# Bayesian Inference and MCMC (Bob Carpenter) 2:03:10 (t=7390)
YouTubeVideo('qQFF4tPgeWI')
# Hamiltonian Monte Carlo and Stan (Michael Betancourt) 58:25 (t=3507)
YouTubeVideo('pHsuIaPbNbY')
# HT: Thomas Wiecki
from IPython.display import IFrame
IFrame('http://twiecki.github.io/blog/2014/01/02/visualizing-mcmc/', width='100%', height=1000)
# John Salvatier: "Bayesian inference with PyMC 3" in PyData Seattle 2015
YouTubeVideo('VVbJ4jEoOfU')
# Thomas Wiecki: "Probabilistic Programming in Python" in EuroPython 2014
YouTubeVideo('KqTUNJ1smYM')
# Thomas Wiecki: "Bayesian Data Analysis with PyMC3" in PyData NYC 2013
IFrame('http://twiecki.github.io/blog/2013/12/12/bayesian-data-analysis-pymc3/', width='100%', height=400)
# Ehsan Karim: "Stan for the beginners [Bayesian inference] in 6 mins "
YouTubeVideo('tLprFqSWS1w')
# Bob Carpenter: "Bayesian Inference and MCMC" in MLSS Sydney 2015
YouTubeVideo('6NXRCtWQNMg')
# Michael Betancourt: "Hamiltonian Monte Carlo and Stan" in MLSS Iceland 2014
YouTubeVideo('xWQpEAyI5s8')