2つの二項分布のパラメタの差の推定(PyMC3, PyStan)
二項分布から生成されると仮定される2つのデータがあったときにその差(2つの二項分布のパラメタの差)を推定することを考える。
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, Model, Deterministic
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
PyMC3でモデルを定義する。ここでは両方とも10回の施行で、それぞれ5回と7回が観察されたものとする。 パラメタの差を表現するにはDeterministicを使う。
In [3]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=10, p=theta1, observed=5)
k2 = Binomial("k2", n=10, p=theta2, observed=7)
MCMCする。
In [4]:
with model:
trace = sample(5000)
結果を表示する。
In [5]:
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
同じことを違うデータでやってみる。試行回数と観察回数をそれぞれ10倍する。
In [6]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=100, p=theta1, observed=80)
k2 = Binomial("k2", n=100, p=theta2, observed=70)
trace = sample(5000)
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
事後分布の左右の幅が狭くなったのがわかる。
次に試行回数が異なってそれぞれゼロ回の観察だった場合。
In [7]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=1, p=theta1, observed=0)
k2 = Binomial("k2", n=5, p=theta2, observed=0)
trace = sample(5000)
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
両方とも0なので、差が0のところが最もありそうなのだが、k1の試行回数が少ない分確信がないので右の裾が長い。
PyStan¶
同じことをPyStanでやってみる。
In [8]:
import pystan
In [9]:
code = """
data {
int<lower=0> k1;
int<lower=1> n1;
int<lower=0> k2;
int<lower=1> n2;
}
parameters {
real<lower=0,upper=1> theta1;
real<lower=0,upper=1> theta2;
}
transformed parameters {
real<lower=-1,upper=1> delta;
delta <- theta1 - theta2;
}
model {
theta1 ~ beta(1, 1);
theta2 ~ beta(1, 1);
k1 ~ binomial(n1, theta1);
k2 ~ binomial(n2, theta2);
}
"""
data = {'k1': 5, 'n1': 10, 'k2': 7, 'n2': 10}
In [10]:
fit = pystan.stan(model_code=code, data=data)
In [11]:
print(fit)
fit.plot()
Out[11]: