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)
 
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -2.4075: 100%|██████████| 200000/200000 [00:24<00:00, 8220.59it/s]
Finished [100%]: Average ELBO = -2.4025
100%|██████████| 5000/5000 [00:04<00:00, 1246.86it/s]
 

summary関数やtraceplot関数で結果を表示することができる。

In [5]:
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
 
theta:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.494            0.129            0.003            [0.248, 0.756]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.234          0.407          0.497          0.581          0.745

 
 

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)
 
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_79f6d134e0766f90fac665338011aa8a NOW.
 

stan関数の戻り値をprintしたりplotメソッドを呼ぶと結果を表示できる。

In [9]:
print(fit)
fit.plot()
 
Inference for Stan model: anon_model_79f6d134e0766f90fac665338011aa8a.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
theta    0.5  2.5e-3   0.14   0.24    0.4    0.5    0.6   0.75   2907    1.0
lp__   -8.82    0.02   0.71 -10.93  -9.01  -8.55  -8.37  -8.32   2050    1.0

Samples were drawn using NUTS at Sat Jun  3 12:48:06 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Out[9]:
 
 

なんで2つずつ出るのだろうか。これはよくわからない。