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

概要

前回は、PyMC2 向けのチュートリアルを PyMC3 に書き換えることでPyMC3 に入門してみました。
今回は、PyMC3 のチュートリアルを見て、実際にモデルを記述する時どういった流れになるか見てみようと思います。

チュートリアル

PyMC3 チュートリアル
http://nbviewer.ipython.org/github/pymc-devs/pymc/blob/master/pymc/examples/tutorial.ipynb

このチュートリアルを頭から見てみましょう。
まずはサンプルの生成です。

import pymc as pm
import numpy as np

ndims = 2
nobs = 20

xtrue = np.random.normal(scale=2., size=1)
ytrue = np.random.normal(loc=np.exp(xtrue), scale=1, size=(ndims, 1))
zdata = np.random.normal(loc=xtrue + ytrue, scale=.75, size=(ndims, nobs))

${x_{true}}$ は ${N(0, 2^2)}$ 、つまり平均0、標準偏差2の正規分布の1次元正規分布に従う1つの値です。
${y_{true}}$ は ${N(exp(x_{true}), 1)}$ 、つまり平均 ${exp(x_{true})}$ 、標準偏差1の正規分布の1次元正規分布に従う2つの値です。サンプリングされた2つの値を並べて2次元ベクトルとしています。
仮に ${y_{true} = (y_1, y_2)}$ とします。

${z_{data}}$ は ${N((x_{true}+y_1, x_{true}+y_2), 0.75)}$ つまり平均 ${(x_{true}+y_1, x_{true}+y_2)}$ 、共分散行列が単位行列の0.75倍である2次元正規分布に従う20個の値です。

今回の実験では、 ${x_{true} = -2.99736635, y_{true} = (-0.21840065, -0.53160934)}$ という値が得られました。

次にモデリングを見てみましょう。

with pm.Model() as model:
    x = pm.Normal('x', mu=0., sd=1)
    y = pm.Normal('y', mu=pm.exp(x), sd=2., shape=(ndims, 1))
    z = pm.Normal('z', mu=x + y, sd=.75, observed=zdata)

今回の目的は、この ${z_{data}}$ が得られていて、 ${N((x+y_1, x+y_2), 0.75)}$ というモデルに従っていると仮定した時、 ${x}$ と ${y}$ の事後分布を計算することです。
この事後分布が真の値を中心に広がっていれば成功、ということです。

${x}$ と ${y}$ がどんな分布に従っているかは分からないので、とりあえずそれぞれ

${x}$ : 平均0、標準偏差1の正規分布
${y}$ : 平均 ${exp(x)}$ 、標準偏差2の正規分布

を事前分布として仮定します。
モデリングする側からすれば知りようのないことですが、真の分布とは異なる分布に従っていますね。

このモデルを用いてサンプリングを行うのですが、このチュートリアルではその初期値をMAP推定法(最大事後確率推定法)を用いて設定しています。

with model:
    start = pm.find_MAP()

print("MAP found:")
print("x:", start['x'])
print("y:", start['y'])

print("Compare with true values:")
print("ytrue", ytrue)
print("xtrue", xtrue)

出力はこんな感じです。

MAP found:
x: -1.3848882353856375
y: [[-1.90266947]
 [-2.02702944]]
Compare with true values:
ytrue [[-0.21840065]
 [-0.53160934]]
xtrue [-2.99736635]

最大事後確率推定法は、簡単に言えば最尤推定法に罰則を付けたものです。
真の値からは少しずれていますが、初期値としては良い感じ、だと思います。

では、サンプリングを行いましょう。

with model:
    step = pm.NUTS()

PyMC3 ではサンプラーの設定を行います。
NUTS は No-U-Turn Sampler の略で、最近考案されたサンプリング手法です。

with model:
    trace = pm.sample(10000, step, start)

今回は1種類ですが、変数毎にサンプラを設定することも可能です。

このサンプラーと初期値を用いて、10000回のサンプリングを行います。(チュートリアルでは3000回ですが増やしました。)

結果を可視化するとこんな感じになります。

f:id:breakbee:20140809210636j:plain

今回の結果だと、 ${x}$ は少し真値からずれたところを中心に分布していますが、 ${y}$ は真値に近いところを中心にしていますね。サンプルが20個と少ないので、もう少し増やせば精度も上がるのではないかと思います。
こんな感じで PyMC3 ではMCMCサンプリングを行います。

PyMC3 の長所と短所

PyMC3 の長所ですが、

Pythonだけでサンプリングを行える
モデルの記述が簡潔

というところだと思います。実際、

with pm.Model() as model:
    x = pm.Normal('x', mu=0., sd=1)
    y = pm.Normal('y', mu=pm.exp(x), sd=2., shape=(ndims, 1))
    z = pm.Normal('z', mu=x + y, sd=.75, observed=zdata)
    start = pm.find_MAP()
    step = pm.NUTS()
    trace = pm.sample(10000, step, start)

の7行でモデルの記述と初期値をサンプラーの設定、サンプリングを行えています。
実際に自分でモデルを記述する時も、このチュートリアルを書き換える形で進めることができるんじゃないかと思います。

PyMC3 の短所ですが、

サンプラの設定が面倒
まだまだバグが多い

といったところだと思います。
BUGS や Stan といったソフトだと、サンプラの設定は自動でしてくれるものが多いので、これがちょっと面倒かと思います。実際、モデルに合わないサンプラを使うと容赦なくバグがでます。

また、PyMC3 はアルファ版なのでまだまだバグが多いです。いくつか自分でもモデリングしてみたのですが、何個かバグに引っかかりました。(サンプリング時に初期値から動かない等)

Github の issue と stackoverflow にバグレポートが集まっているので、モデリングした時挙動がおかしかったら調べてみるといいのではないかと思います。

感想

PyMC3 はモデルの記述が簡潔で便利なのですが、まだまだバグも多いので、これからの発展に期待したいと思います。
Github の README にもあるように、初めて PyMC に触る方は PyMC2 のほうが安定してていいかもしれません。

次は PyStan について調べたことをまとめようかと思います。