pythonだけで、ベイズ推定したい(pymcを使う_入門)

ベイジアンモデリングといえば、pystanでstanをキックすることが多いですが、Pythonで書かれたpymcはPythonだけで完結できます。

pymcの使い方を学習していくことにします。

基本的な書き方

withの中でモデルを定義していきます。
各種分布がpymcの中には用意されており、それを順に定義していきます。
データに関してはobservedパラメータから渡すことができます。

そのあと、そのモデルからサンプリングすることができます。

import pymc3 as pm

#モデルの定義
with pm.Model() as mdl:
  theta = pm.Beta('theta',alpha=1,beta=1)
  y = pm.Bernoulli('y',p=theta,observed=data)

#サンプリング
with mdl:
  start = pm.find_MAP()
  step = pm.Metropolis()
  trace = pm.sample(1000,start=start,step=step,cores=4)

やってみる

データの発生

ベルヌーイ分布から100ほどデータを発生させます。

import numpy as np
import pymc3 as pm
import scipy.stats as stats

theta_real = 0.3
data = stats.bernoulli.rvs(p=theta_real,size=100)

モデリングと結果の確認

上記と同じようにモデルを定義します。
thetaはベータ分布にしておき、それをベルヌーイ分布に渡してサンプリングします。

with pm.Model() as mdl:
  theta = pm.Beta('theta',alpha=1,beta=1)
  y = pm.Bernoulli('y',p=theta,observed=data)

with mdl:
  start = pm.find_MAP()
  step = pm.Metropolis()
  trace = pm.sample(1000,start=start,step=step,cores=4)

バーンイン期間はスライスで落とします。
様々な可視化が簡単にできます。

burnin = 100
chain = trace[burnin:]

#トレースプロット
pm.traceplot(chain)
#フォレストプロット
pm.forestplot(chain,varnames=['theta'])
#自己相関プロット
pm.autocorrplot(chain)

推定結果もみてみます。

pm.summary(chain)

 mean    sd  hdi_3%  hdi_97% mcse_mean   mcse_sd ess_mean    ess_sd  ess_bulk    ess_tail    r_hat
theta   0.325   0.237   0.0 0.751   0.009   0.006   694.0   694.0   641.0   668.0   1.0

それっぽいのが推定できました。