PyMCの初歩

Python
 

PyMCを試すことにする。 Pythonで体験するベイズ推論 を買った関係上、PyMC3ではなくてPyMC2を使っている。

この記事の内容は書籍と関係なく、書籍よりずっと初歩的なものである。

In [1]:
import pymc as pm
from pymc.Matplot import plot as mcplot
from pymc import Beta, Binomial, Uniform, Model, MCMC
from pymc.Matplot import plot as mcplot

from matplotlib import pyplot as plt
%matplotlib inline
 

データが二項分布によって生成されるモデルを考える。 二項分布では事象発生の確率を示すパラメタ$\theta$があって、これが推定の対象となる。

$\theta$は0から1までの連続値であり、特に事前の情報を持たないので、 0から1までのすべてが等確率（密度）となる事前分布を使う。 これは具体的にはパラメタ$\alpha$と$\beta$がそれぞれ1のベータ分布である。

学習に使うデータとして、10回中5回事象が起こったとしよう。

In [2]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=10, p=theta, value=5, observed=True)
model = Model([theta, k])
 

このモデルを使ってMCMCを行い、結果をプロットする。

In [3]:
mcmc = MCMC(model)
mcmc.sample(40000, 10000)

mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.6 secPlotting $\theta$
 
 

当然だが右側のプロットで示されるように、推定値は0.5の周りになる。 95%区間は点線で示されている。

10回中5回でなく、100回中50回だったら0.5に近いことにもっと確証を持つだろう。

In [4]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=100, p=theta, value=50, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 4.0 secPlotting $\theta$
 
 

95%区間が狭くなったことがわかる。

1回中1回、1回中0回、2回中1回も試してみる。

In [5]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=1, p=theta, value=1, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.9 secPlotting $\theta$
 
In [6]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=1, p=theta, value=0, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.8 secPlotting $\theta$
 
In [7]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=2, p=theta, value=1, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.9 secPlotting $\theta$
 
 

ところで二項分布のパラメタの事前分布としてベータ分布を使うというのは殆どの文献で慣例的にそうなっているのだが、 等確率分布ならばUniformを使ってもよいのだろうか。

In [8]:
theta = Uniform('$\\theta$', lower=0, upper=1)
k = Binomial("k", n=10, p=theta, value=5, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.7 secPlotting $\theta$
 
In [9]:
theta = Uniform('$\\theta$', lower=0, upper=1)
k = Binomial("k", n=1, p=theta, value=1, observed=True)
model = Model([theta, k])
mcmc = MCMC(model)
mcmc.sample(40000, 10000)
mcplot(mcmc.trace('$\\theta$'))
 
 [-----------------100%-----------------] 40000 of 40000 complete in 3.8 secPlotting $\theta$
 
 

結果に特に違いはないように見える。

In [ ]: