『Bayesian Methods For Hackers』Ch. 3 まとめ
Ch. 3です。
ベイズ推定と変数空間
N変数のベイズ推定をするとき、我々はN次元空間を考えている。最初、このN次元空間には事前分布により定義されたsurfaceがある(2次元の場合を考えるとイメージしやすい)。ベイズ推定とは、このsurfaceをデータにより変形していくプロセスのことである。すべてのデータを使われた後、パラメータがどこに住んでいるかを反映した事後分布ができる。
事後分布を"押し上げる"レートは事前分布によって決まる。事前確率が低いということは、データに対して抵抗が強い(上がりにくい)ことになる。特に重要な点として、事前確率が0ならその座標での事後確率も必ず0になる(ベイズの定理を考えれば明らか)。
MCMCでは、"賢く"N次元空間を探索することができる。MCMCは事後分布を直接求めるのではなく、事後分布からのサンプルを得る。これは一見非効率なやり方のような気がするが、
- 事後分布の関数形を求める(めちゃくちゃ大変)
- 事後分布のピーク位置を求める(どのくらいの不確定性があるかわからない)
などと比べると良いアプローチである。
MCMCのアルゴリズム
MCMCでは以下の作業によりサンプリングを行う。
- 現在位置から出発
- 次の位置を提案する
- 次の位置への移動をaccept/rejectして1に戻る
- 上記を何回も繰り返して、すべての座標を返す
PyMCによるMCMC
MCMC
の引数
mcmc.sample(10000, burn = 5000, thinning = 5)
- burn: 最初の方のサンプルは初期位置の影響を受けているので捨てたい(burn-in)。全サンプルの半分くらいが良いらしい(サンプル数が多い時は90%にしてもOK)
- thinning: MCMCによるサンプリングは自己相関が強いので、これを弱める(分布を"混ぜ"てサンプルを独立に近づける)ために何個かおきにサンプルする。10以上にしても意味が薄いのでそれ以下の適当な値がいいらしい
MAP
MCMCの結果は初期値に強く依存し、悪い初期値からサンプルを開始すると定常分布に収束する(初期値を"忘れる")までに時間がかかる。これを防ぐために、事後分布のピーク(maximum a posterior, MAP)からサンプルをスタートしたい。これを行うのがpm.MAP
である。
map_ = pm.MAP(model)
map_.fit()
mcmc = pm.MCMC(model)
mcmc.sample(10000)
MAPを求めるアルゴリズムはいくつかあるがデフォルトはscipy
のfmin
。基本はこれで良いが、うまくいかなそうだったらPowell’s methodをfit
の引数に指定して試してみると良い。
pymc.Matplot.plot()
Matplot
を使うとパラメータのヒストグラム、信用区間。自己相関などを一気にプロットしてくれて便利。引数はtrace
。
from pymc.Matplot import plot as mcplot mcplot(mcmc.trace('centers', 2))