複数の観察データからのパラメタ推定 (PyMC3, PyStan)

 

今度は二項分布のパラメタが1つだけあって、同じ二項分布から複数のセッション(n回中k回)がデータとして取られる場合に、 その1つだけのパラメタを推定するということを考える。

In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
%matplotlib inline
 

PyMC3

In [2]:
import pymc3 as pm
from pymc3 import Beta, Binomial, Bernoulli, Model, Deterministic
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
 

PyMC3でモデルを定義する。 観察されるデータはk1をk2の2つあり、同じパラメタthetaの二項分布から生み出されているとする。2つの観察データは10回中5回と10回中7回であった。

In [3]:
model = Model()
with model:
    # Prior
    theta = Beta('theta', alpha=1, beta=1)

    # Observed
    k1 = Binomial("k1", n=10, p=theta, observed=5)
    k2 = Binomial("k2", n=10, p=theta, observed=7)

    trace = sample(5000)

summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
 
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -4.4786: 100%|██████████| 200000/200000 [00:23<00:00, 8437.65it/s]
Finished [100%]: Average ELBO = -4.4647
100%|██████████| 5000/5000 [00:03<00:00, 1276.28it/s]
 
theta:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.591            0.099            0.003            [0.386, 0.778]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.381          0.529          0.596          0.658          0.776

 
 

上記ではk1とk2という2つのBinomialを記述したが、下記のように1つのBinomialでobservedを配列にする書き方もできるようだ。

In [4]:
model = Model()
with model:
    # Prior
    theta = Beta('theta', alpha=1, beta=1)

    # Observed
    k = Binomial("k", n=10, p=theta, observed=[5, 7])
    trace = sample(5000)
    
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
 
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -4.4718: 100%|██████████| 200000/200000 [00:26<00:00, 7629.06it/s]
Finished [100%]: Average ELBO = -4.4784
100%|██████████| 5000/5000 [00:04<00:00, 1172.66it/s]
 
theta:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.590            0.101            0.002            [0.394, 0.788]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.380          0.525          0.594          0.660          0.780

 
 

ところで10回中5回、7回の2つの二項分布のセッションがあるというのは、20個のベルヌーイ分布で12回が1であるという見方をしてもいいのだろうか。

In [5]:
model = Model()
with model:
    # Prior
    theta = Beta('theta', alpha=1, beta=1)

    # Observed
    k = Bernoulli("k", p=theta, observed=[1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0])
    trace = sample(5000)
    
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
 
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -14.793: 100%|██████████| 200000/200000 [00:24<00:00, 8070.62it/s]
Finished [100%]: Average ELBO = -14.791
100%|██████████| 5000/5000 [00:04<00:00, 1211.14it/s]
 
theta:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.595            0.100            0.002            [0.396, 0.778]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.393          0.528          0.597          0.666          0.777

 
 

この場合でもパラメタの推定値が同様であることがわかる。

 

PyStan

同じことをPyStanでやってみる。まずはk1とk2を別々に定義する書き方から。

In [6]:
import pystan
In [7]:
code = """
data { 
  int<lower=0> k1;
  int<lower=1> n1;
  int<lower=0> k2;
  int<lower=1> n2;
} 
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1, 1);
  k1 ~ binomial(n1, theta);
  k2 ~ binomial(n2, theta);
}
"""

data = {'k1': 5, 'n1': 10, 'k2': 7, 'n2': 10}
In [8]:
fit = pystan.stan(model_code=code, data=data)
print(fit)
 
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_8bbd120c2f3e5793a7ec522da7d4265d NOW.
 
Inference for Stan model: anon_model_8bbd120c2f3e5793a7ec522da7d4265d.
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.59  1.9e-3    0.1   0.38   0.53    0.6   0.66   0.78   2781    1.0
lp__  -15.39    0.02   0.71 -17.42 -15.55 -15.11 -14.93 -14.88   2032    1.0

Samples were drawn using NUTS at Sun Jun  4 01:18: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).
 

Stanでも配列を使うことができる。forを使ってインデクシングするようだ。

In [9]:
code = """
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
} 
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1, 1);
  for (i in 1:N)
    k[i] ~ binomial(n[i], theta);
}
"""

data = {'N': 2, 'k': [5, 7], 'n': [10, 10]}

fit = pystan.stan(model_code=code, data=data)
print(fit)
 
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_8ee3da0050cf86e0c9c1a53e4765a6c7 NOW.
 
Inference for Stan model: anon_model_8ee3da0050cf86e0c9c1a53e4765a6c7.
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.59  2.1e-3    0.1   0.39   0.52   0.59   0.66   0.77   2295    1.0
lp__  -15.38    0.02   0.72 -17.34 -15.52  -15.1 -14.93 -14.88   1978    1.0

Samples were drawn using NUTS at Sun Jun  4 01:19:16 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).