PyMC3でMCMC入門(1) - old school magic

概要

先日、Tokyo.scipy というイベントがありました。

Python で科学技術計算を用いる方々の勉強会だそうです。
私は参加していないのですが、PyMC に関するセッションがあったそうです。

機械学習やデータマイニングで有名な神嶌先生の発表資料です。
非常に分かりやすいスライドでした。
このセッションで用いている PyMC は PyMC2 ですが、せっかくなので PyMC3 で書き換えてみました。PyMC3 の入門にちょうどいいのではないかと思います。

PyMC3 のインストール

こちらの記事をご参照ください。
Python3でPyMCのインストール - old school magic

サンプル1

サンプル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)で保存してます。

といった感じです。

結果はこんな感じです。
f:id:breakbee:20140804024407j:plain

真の平均は1.0ですが、この結果だと0.94くらいを中心に分布してるように見えますね...
でも元のサンプルも0.9を中心に分布してるので大丈夫...だと思います。

サンプル2

サンプル2(1入力，1出力のガウス単純ベイズ)

簡単なナイーブベイズを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 を用いて、ついでにコードを簡潔にしました。

(もし間違ってたらご指摘ください。)