In [284]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm2
import pymc3 as pm
import time
import math
import numpy.random as rd
import pandas as pd
from pymc3 import summary
In [271]:
from pymc3.backends.base import merge_traces
import theano.tensor as T
メトロポリス法
（１）パラメーターqの初期値を選ぶ
（２）qを増やすか減らすかをランダムに決める
（３）q（新）において尤度が大きくなるならqの値をq（新）に変更する
（４）q（新）で尤度が小さくなる場合であっても、確率rでqの値をq（新）に変更する
In [ ]:
In [116]:
def comb(n, r):
if n == 0 or r == 0: return 1
return comb(n, r - 1) * (n - r + 1) / r
def prob(n, y, q):
p = comb(n, y) * q ** y * (1 - q) ** (n - y)
return p
def likelighood(n, y, q):
p = 1.0
for i in y:
p = p*prob(n, i, q)
return p
def metropolis(n, y, q, b, num):
qlist = np.array([q])
for i in range(num):
old_q = q
q = q+np.random.choice([b, -b])
old_l = likelighood(n, y, old_q)
new_l = likelighood(n, y, q)
if new_l > old_l:
old_q = q
else:
r = new_l/old_l
q = np.random.choice([q, old_q], p=[r, 1.0-r])
q = round(q, 5)
qlist = np.append(qlist, q)
return q, qlist
In [122]:
y = [4, 3, 4, 5, 5, 2, 3, 1, 4, 0, 1, 5, 5, 6, 5, 4, 4, 5, 3, 4]
q, qlist = metropolis(8, y, 0.3, 0.01, 10000)
In [129]:
plt.plot(qlist)
Out[129]:
In [128]:
plt.hist(qlist)
Out[128]:
In [245]:
qlist.mean()
Out[245]:
In [348]:
N = 40
X = np.random.uniform(10, size=N)
Y = X*30 + 4 + np.random.normal(0, 16, size=N)
plt.plot(X, Y, "o")
Out[348]:
In [349]:
multicore = False
saveimage = False
itenum = 1000
t0 = time.clock()
chainnum = 3
with pm.Model() as model:
alpha = pm.Normal('alpha', mu=0, sd =20)
beta = pm.Normal('beta', mu=0, sd=20)
sigma = pm.Uniform('sigma', lower=0)
y = pm.Normal('y', mu=beta*X + alpha, sd=sigma, observed=Y)
start = pm.find_MAP()
step = pm.NUTS(state=start)
In [352]:
with model:
if(multicore):
trace = pm.sample(itenum, step, start=start,
njobs=chainnum, random_seed=range(chainnum),
progress_bar=False)
else:
ts = [pm.sample(itenum, step, chain=i, progressbar=False) for i in range(chainnum)]
trace = merge_traces(ts)
if(saveimage):
pm.tracepot(trace).savefig("simple_linear_trace.png")
print "Rhat="+str(pm.gelman_rubin(trace))
t1=time.clock()
print "elapsed time=" + str(t1-t0)
In [293]:
if(not multicore):
trace = ts[0]
with model:
pm.traceplot(trace, model.vars)
In [294]:
pm.forestplot(trace)
Out[294]:
In [302]:
summary(trace)
In [353]:
multicore = True
t0 = time.clock()
with model:
if(multicore):
trace = pm.sample(itenum, step, start=start,
njobs=chainnum, random_seed=range(chainnum),
progress_bar=False)
else:
ts = [pm.sample(itenum, step, chain=i, progressbar=False) for i in range(chainnum)]
trace = merge_traces(ts)
if(saveimage):
pm.tracepot(trace).savefig("simple_linear_trace.png")
print "Rhat="+str(pm.gelman_rubin(trace))
t1=time.clock()
print "elapsed time=" + str(t1-t0)
if(not multicore):
trace = ts[0]
with model:
pm.traceplot(trace, model.vars)
例題：個体差と生存種子数
ある植物を考える。i番目の個体の生存種子数をyiとする。yiは０以上8以下である。以下はヒストグラムである。
In [297]:
data = pd.read_csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/data7a.csv")
plt.bar(range(9), data.groupby('y').sum().id)
data.groupby('y').sum().T
Out[297]:
種子生存確率が9通の生存確率qの二項分布で説明できるとする。
実際のデータを説明できていない。
・・・個体差を考慮できるGLMMを用いる。logit(qi) = β + ri
切片βは全個体に共通するパラメーター、riは個体差で平均0、標準偏差sの正規分布に従う。事後分布∝p(Y|β, {ri})×事前分布
βの事前分布には無情報事前分布を指定する。
p(β)=1/√2π×100^2 × exp(-β^2/2×100^2)
In [324]:
plt.hist(np.random.normal(0, 100, 1000))
Out[324]:
riの事前分布には平均0、標準偏差ｓの正規分布を仮定する。
p(ri|s)=1/√2π×s^2 × exp(-ri^2/2×s^2)sの事前分布には無情報事前分布を指定する。
p(s)=(０から10^4までの連続一様分布)
In [298]:
Y = np.array(data.y)[:6]
yの数が大きく全て使うと時間がかかりすぎるので6体だけ選び出す。
In [299]:
def invlogit(v):
return T.exp(v)/(T.exp(v) + 1)
with pm.Model() as model_hier:
s = pm.Uniform('s', 0, 1.0E+2)
beta = pm.Normal('beta', 0, 1.0E+2)
r = pm.Normal('r', 0, s, shape=len(Y))
q = invlogit(beta+r)
y = pm.Binomial('y', 8, q, observed=Y) #p(q|Y)
step = pm.Slice([s, beta, r])
trace_hier = pm.sample(1000, step)
with model_hier:
pm.traceplot(trace_hier, model_hier.vars)
In [322]:
summary(trace_hier)
In [323]:
trace_hier
Out[323]:
In [278]:
x_sample = np.random.normal(loc=1.0, scale=1.0, size=1000)
with pm.Model() as model:
mu = pm.Normal('mu', mu=0., sd=0.1)
x = pm.Normal('x', mu=mu, sd=1., observed=x_sample)
with model:
start = pm.find_MAP()
step = pm.NUTS()
trace = pm.sample(10000, step, start)
pm.traceplot(trace)
plt.savefig("result1.jpg")
In [282]:
ndims = 2
nobs = 20
n = 1000
y_sample = np.random.binomial(1, 0.5, size=(n,))
x_sample=np.empty(n)
x_sample[y_sample==0] = np.random.normal(-1, 1, size=(n, ))[y_sample==0]
x_sample[y_sample==1] = np.random.normal(1, 1, size=(n, ))[y_sample==1]
with pm.Model() as model:
p = pm.Beta('p', alpha=1.0, beta=1.0)
y = pm.Bernoulli('y', p=p, observed=y_sample)
mu0 = pm.Normal('mu0', mu=0., sd=1.)
mu1 = pm.Normal('mu1', mu=0., sd=1.)
mu = pm.Deterministic('mu', mu0 * (1-y_sample) + mu1 * y_sample)
x = pm.Normal('x', mu=mu, sd=1., observed=x_sample)
with model:
start = pm.find_MAP()
step = pm.NUTS()
trace = pm.sample(10000, step, start)
pm.traceplot(trace)
plt.savefig("result2.jpg")