pymc3を始める。
ベイズ統計モデリングをpymc3を使って学ぶ
最近、ベイズ統計モデリングに興味があり、勉強をはじめた。学んだ結果の記録も兼ねて、ブログもやってみることに。
幾つかのWebや書籍を調べてみると、ツールとしてはstanが主流の模様。Rやpythonのラッパもあるようだが、一旦コンパイルが必要など少し面倒な印象。ここは、好きなpythonオンリーで書けるpymc3で学ぶことにする。ベイズ統計モデルとpymcの勉強だが、どちらかと言うとpymcに軸足を置き、stanコードのpymc化や、公式HPのexampleをトレースしつつ学んでいこう、という感じ。
pymc3のインストール
anacondaを使っていれば、簡単にインストールできる。
conda install -c conda-forge pymc3
現在のバージョンは3.5で、リリースノートによると幾つかの機能アップデートがあった模様。 個人的に大きいと感じた変更は以下。
簡単なサンプル（2枚のコインの差）
表が出る確率が異なる2枚のコインについて、コイン投げの結果からコイン2の方が表が出やすい確率をベイズ統計で求める。
# 2枚のコインのデータを生成する。 n = 50 # コイン投げ回数 trials1 = stats.bernoulli.rvs(0.5, size=n) trials2 = stats.bernoulli.rvs(0.55, size=n) #コイン2の方が表がちょっと出やすい。 h1 = sum(trials1) h2 = sum(trials2) # 最尤推定 print(h1/n, h2/n)
0.44 0.52
# pymc3によるモデル化 with pm.Model() as model: p1 = pm.Uniform('p1', lower=0.0, upper=1.0) y1 = pm.Binomial('y1', n=n, p=p1, observed=h1) p2 = pm.Uniform('p2', lower=0.0, upper=1.0) y2 = pm.Binomial('y2', n=n, p=p2, observed=h2) # 2枚目のコインと1枚目のコインの差 delta = pm.Deterministic('delta', p2-p1) # モデルのグラフィカル表示（3.5から追加された機能） pm.model_to_graphviz(model)
# MCMCの実行 with model: trace = pm.sample(2000, tune=500, step=pm.NUTS())
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [p2, p1]
Sampling 2 chains: 100%|██████████| 5000/5000 [00:03<00:00, 1275.23draws/s]
# 収束チェック
pm.traceplot(trace)
トレースプロットをみると収束していそう。
# 結果の出力
pm.summary(trace)
Rhatも問題なさそう。
| mean | sd | mc_error | hpd_2.5 | hpd_97.5 | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|
| p1 | 0.45 | 0.07 | 1.10e-03 | 0.31 | 0.57 | 4308.30 | 1.0 |
| p2 | 0.52 | 0.07 | 1.05e-03 | 0.38 | 0.64 | 4320.51 | 1.0 |
| delta | 0.07 | 0.09 | 1.36e-03 | -0.11 | 0.25 | 4576.03 | 1.0 |
# 事後分布の表示
pm.plot_posterior(trace)
# delta>0の確率は？ sum(trace['delta'] > 0)/len(trace['delta'])
0.78
50回のコイン投げ結果から、コイン2の方が表が出やすい確率は78%と推定された。
以上
