!date
Mon Jun 16 08:14:43 PDT 2014
http://stackoverflow.com/questions/24242660/pymc3-multiple-observed-values
import matplotlib.pylab as plt
%matplotlib inline
import pymc3 as pm # fork of pymc to allow importing pymc2 as well, when necessary
# I prefer to use pm.*, but this keeps consistency with the SO question
from pymc3 import Model, Uniform, Normal, Poisson, Metropolis, traceplot
from pymc3 import sample
# Simulate data
import scipy.stats
totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates)
successRate = 0.1*totalRates
successCounts = scipy.stats.poisson.rvs(mu=successRate)
with Model() as success_model:
total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
total = Poisson('total', mu=total_lambda, observed=totalCounts)
loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
error = Poisson('error', mu=total_lambda*loss_lambda_factor, observed=successCounts)
with success_model:
step = Metropolis()
success_samples = sample(20000, step) #, start)
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples)
[-----------------100%-----------------] 20000 of 20000 complete in 41.3 sec
<matplotlib.figure.Figure at 0x7f95c40c5150>
Problem: including the burn-in in the plot makes it hard to see what is happening.
Solution: don't include.
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples[10000:])
<matplotlib.figure.Figure at 0x7f95a894ced0>
Now we can see problem 2: non-convergence
Solution: change priors, start from MAP(?), use step method that fits space better
with Model() as success_model:
total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100)
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000)
total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
total = Poisson('total', mu=total_lambda, observed=totalCounts)
loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
error = Poisson('error', mu=total_lambda*loss_lambda_factor, observed=successCounts)
with success_model:
step = Metropolis()
success_samples = sample(20000, step) #, start)
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples[10000:])
[-----------------100%-----------------] 20000 of 20000 complete in 38.4 sec
<matplotlib.figure.Figure at 0x7f95a3b38950>
with success_model:
start = pm.find_MAP()
step = pm.CompoundStep([pm.Metropolis([total_lambda_mu]),
pm.Metropolis([total_lambda_tau]),
pm.Metropolis([total_lambda]),
pm.Metropolis([loss_lambda_factor]),
])
success_samples = sample(20000, step, start)
[-----------------100%-----------------] 20000 of 20000 complete in 156.9 sec
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples[10000:])
<matplotlib.figure.Figure at 0x7f95a81e9dd0>
with success_model:
step = pm.NUTS()
success_samples = sample(20000, step)
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples[10000:])
[-----------------100%-----------------] 20000 of 20000 complete in 1788.7 sec
/homes/abie/anaconda/lib/python2.7/site-packages/Theano-0.6.0-py2.7.egg/theano/scan_module/scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility from scan_perform.scan_perform import *
<matplotlib.figure.Figure at 0x7f95a89afc90>
Problem 3: The model
import theano.tensor as T
with Model() as success_model:
loss_lambda_rate = pm.Flat('loss_lambda_rate')
error = Poisson('error', mu=totalCounts*T.exp(loss_lambda_rate), observed=successCounts)
with success_model:
step = Metropolis()
success_samples = sample(20000, step) #, start)
plt.figure(figsize=(10, 10))
_ = traceplot(success_samples[10000:])
[-----------------100%-----------------] 20000 of 20000 complete in 34.3 sec
<matplotlib.figure.Figure at 0x7f959a1c2090>
import numpy as np
plt.hist(np.exp(success_samples['loss_lambda_rate'][10000:]));