PyMC3とPyStanによる二項分布のパラメタ推定
前回の記事で二項分布のパラメタ推定をPyMC2で行った。 しかしPyMC2をいまから使っていくのも微妙な気がしてきたので新しいPyMC3で書き直す。 ついでにPyStanでも同じことをして比べてみる。
In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
%matplotlib inline
PyMC3¶
まずはPyMC3のモジュールをインポートする。 theano.config.warn.round=False
は警告を抑止するための設定。
In [2]:
import pymc3 as pm
from pymc3 import Beta, Binomial, Model
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
次にモデルを定義する。
In [3]:
model = Model()
with model:
theta = Beta('theta', alpha=1, beta=1)
k = Binomial("k", n=10, p=theta, observed=5)
PyMC3ではこのようにwith構文を使ってモデルを定義する。 何をしているのかわからなくて面食らうが、見た目上モデル定義がまとまりをなすようにするための工夫だろうか。 多分ダイナミックな文脈を作ってBetaやBinomialなどが現在の文脈と結びつくようにしているのかもしれない。
MCMCを行うにはsample関数を使う。
In [4]:
with model:
trace = sample(5000)
summary関数やtraceplot関数で結果を表示することができる。
In [5]:
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
PyStan¶
次に同じことをPyStanでやってみる。
In [6]:
import pystan
PyStanではモデルは完全にStanの文法で記述することになる。
In [7]:
code = """
data {
int<lower=0> k;
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1, 1);
k ~ binomial(10, theta);
}
"""
data = {'k': 5}
stan関数でMCMCが行われる。
In [8]:
fit = pystan.stan(model_code=code, data=data)
stan関数の戻り値をprintしたりplotメソッドを呼ぶと結果を表示できる。
In [9]:
print(fit)
fit.plot()
Out[9]:
なんで2つずつ出るのだろうか。これはよくわからない。