第3章はMCMCの原理について。
MCMCとは
詳しくはマルコフ連鎖モンテカルロ法 - Wikipediaを見るか、ベイズ統計と統計物理を参照。
MCMCはマルコフ連鎖モンテカルロ法(Markov chain Monte Carlo methods)の略称で、求めたい確率分布に従って動くマルコフ過程によってその分布をサンプリングするアルゴリズム。色々細かく違いがあるアルゴリズムがあるがそれらを総称している(厳密にはマルコフ過程でないものも含む)。具体的なアルゴリズムは簡単で、
- どこか適当にスタート地点を決める。
- 次の移動地点を決める。
- 分布から計算した確率でその点に移動するか决める。
- 繰り返す
だけである。
MCMC以外の分布近似アルゴリズム
があるらしい
例: 混合分布の教師なし分類
モデリングとパラメータの推定
上のような分布があるとき、この分布は2つのピークがあるように見える。これが2つの正規分布を混合したものと考えてそれぞれを推定することを考える。この分布が
- 確率pで分布1から1-pで分布2からサンプリングする。
- それぞれの分布は平均、標準偏差とする。
- pの事前分布は[0, 1]の一様分布。
- はそれぞれ120、190あたりがもっともらしいのでそこを中心とする分散の逆数が0.01の正規分布()。
- はが[0,100]の一様分布とする。
と考えると、
p = pm.Uniform("p", 0, 1) # 0か1を指定した確率に従ってランダムにサンプリングする assignment = pm.Categorical("assignment", [p, 1 - p], size=data.shape[0]) centers = pm.Normal("centers", [120, 190], [0.01, 0.01], size=2) taus = 1 / pm.Uniform("stds", 0, 1, size=2) ** 2 @pm.deterministic def center_i(assignment=assignment, centers=centers): return centers[assignment] @pm.deterministic def tau_i(assignment=assignment, taus=taus): return taus[assignment] observation = pm.Normal("obs", center_i, tau_i, value=data, observed=True) model = pm.Model([p, assignment, taus, centers]) mcmc = pm.MCMC(model) mcmc.sample(50000) # 結果を見て追加のサンプリング mcmc.sample(100000) # 2回めのサンプル centers_trace = mcmc.trace("centers", chain=1)[:] # 1回めのサンプル prev_center_trace = mcmc.trace("centers", chain=0)[:]
MCMC.traceはchainを引数に取る。chainはデフォルトが-1でこのとき最後にMCMC.sampleを読んだ時の結果を返す。 やの最尤推定はtraceしたものの平均を取ればいい。
収束判定
- 自己相関:
- MCMCは原理的にサンプリングするときに自己相関があるサンプルを取ってくる。この自己相関が高すぎる場合に分布が行ったり来たりを繰り返す状態になって収束しないことがある。
- nサンプルに1つしかサンプルとして採用しないことで自己相関を抑えることができる。この場合より多くのサンプリングが必要になりトレードオフが生じる。
- pymc.Matplot.plotを使うことで自己相関をプロットして視覚的に見ることができる。
from pymc.Matplot import plot as mcplot mcplot(mcmc.trace("param"), commmnon_scale=False)
左上がtrace、左下が自己相関、右が分布をプロットしたもの。
MAP(maximum a prior)を使う
MCMC.sampleだけを利用した場合、初期値の位置がうまく選択できていない場合burn in期間を長く取ることになる。初期値を最尤値に近いところで自動的にうまくとってくれるのがMAPである。
# pm.MAP(model).fit()でも良さそう
map_ = pm.MAP(model)
map_.fit()
補足