PyMC3でMCMC入門(1)
概要
先日、Tokyo.scipy というイベントがありました。
Python で科学技術計算を用いる方々の勉強会だそうです。
私は参加していないのですが、PyMC に関するセッションがあったそうです。
機械学習やデータマイニングで有名な神嶌先生の発表資料です。
非常に分かりやすいスライドでした。
このセッションで用いている PyMC は PyMC2 ですが、せっかくなので PyMC3 で書き換えてみました。PyMC3 の入門にちょうどいいのではないかと思います。
PyMC3 のインストール
こちらの記事をご参照ください。
Python3でPyMCのインストール - old school magic
サンプル1
単純なガウス分布の平均パラメータをMCMCで計算しています。
PyMC3 のコードはこちら。(コードが表示されなかったら更新してみてください)
This is PyMC3 version of http://nbviewer.ipython.o ...
PyMC2 のコードとの主な違いですが、
- PyMC3 では with 構文の中でモデルを構築し、サンプリングを行っています。
- モデルを構築する時、PyMC2 では観測変数は value に値を、 observed に観測変数かどうかのフラグを代入しますが、 PyMC3 では observed に直接変数を代入してます。
- このコードでは、サンプリングの初期値として、MAP推定の結果を用いています。
MAP推定とは「最大事後確率推定」の略で、簡単に言うと最尤推定に罰則を付けたものです。(ベイズ推定とは異なります。)
詳しくはPyMC3のチュートリアルをご参照ください。
http://nbviewer.ipython.org/github/pymc-devs/pymc/blob/master/pymc/examples/tutorial.ipynb
- 結果を画像(result1.jpg)で保存してます。
といった感じです。
結果はこんな感じです。
真の平均は1.0ですが、この結果だと0.94くらいを中心に分布してるように見えますね...
でも元のサンプルも0.9を中心に分布してるので大丈夫...だと思います。
サンプル2
簡単なナイーブベイズをMCMCで計算しています。
PyMC3 のコードはこちら。
This is PyMC3 version of http://nbviewer.ipython.o ...
PyMC2 との主な違いですが、
- PyMC3 では deterministic デコレータがありません。
PyMC2 のコードでは、
@deterministic(plot=False) def mu(y=y, mu0=mu0, mu1=mu1): out = np.empty_like(y, dtype=np.float) out[y == 0] = mu0 out[y == 1] = mu1 return out
として確定的なオブジェクトを扱っていましたが、 PyMC3では、pymc.Deterministic を使用します。
mu = pm.Deterministic('mu', mu0 * (1-y_sample) + mu1 * y_sample)
参考: What's @pm.deterministic in pyMC3?
https://github.com/pymc-devs/pymc/issues/434
一応等価なコードになっている...はずです。
PyMC2 のコードでは、y をそのまま計算に使えているのですが、PyMC3 ではこの y はそのまま使えません。
というのも、y は ObservedRV というデータ型で、ようするに観測変数なので計算に使う用途で設計されてないようです。
なので y に束縛してる y_sample を用いて、ついでにコードを簡潔にしました。
(もし間違ってたらご指摘ください。)
- 結果を画像(result2.jpg)で保存してます。
といった感じです。
結果はこんな感じです。
各パラメータともに真値のまわりで分布しているように見えるので多分あってるはず...
感想
PyMC3 はまだ alpha 版ですし、ドキュメントもあまり整備されてないのでちょっと大変ですが、これからも頑張ってみようと思います。
あと、一応続きます。次は公式のチュートリアルの和訳と解説でもやろうかな...
最後に、神嶌先生、素晴らしい資料をありがとうございます。
もし問題があればお知らせください。