複数の観察データからのパラメタ推定 (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()
上記では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()
ところで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()
この場合でもパラメタの推定値が同様であることがわかる。
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)
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)