HSI
From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press
(c) 2017, Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida
you are kindly asked to include the complete citation if you used this material in a publication
​
​
Code 6.7: Bayesian Poisson model using pymc3
====================================================
import numpy as np
import pandas
import pylab as plt
import pymc3 as pm
from scipy.stats import norm, binom, poisson
​
# Data
np.random.seed(18472) # set seed to replicate example
nobs= 750 # number of obs in model
​
x1_2 = binom.rvs(1, 0.7, size=nobs)
x2 = norm.rvs(loc=0, scale=1.0, size=nobs)
xb = 1 - 1.5 * x1_2 - 3.5 * x2 # linear predictor
exb = np.exp(xb)
py = poisson.rvs(exb) # create y as adjusted
​
df = pandas.DataFrame({'x1_2': x1_2, 'x2':x2, 'py': py}) # re-write data
​
# Fit
niter = 10000 # parameters for MCMC
​
with pm.Model() as model_glm:
# define priors
beta0 = pm.Flat('beta0')
beta1 = pm.Flat('beta1')
beta2 = pm.Flat('beta2')
​
# define likelihood
mu = np.exp(beta0 + beta1*x1_2 + beta2 * x2)
y_obs = pm.Poisson('y_obs', mu, observed=py)
​
# inference
start = pm.find_MAP() # Find starting value by optimization
step = pm.NUTS()
trace = pm.sample(niter, step, start, progressbar=True)
​
# Output
pm.summary(trace)
​
# show graphical output
pm.traceplot(trace)
plt.show()
====================================================
Output on screen:
​
beta0:
Mean SD MC Error 95% HPD interval
------------------------------------------------------------------------------
1.002 0.012 0.000 [0.978, 1.024]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.978 0.994 1.002 1.010 1.024
beta1:
Mean SD MC Error 95% HPD interval
------------------------------------------------------------------------------
-1.500 0.006 0.000 [-1.512, -1.489]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-1.512 -1.504 -1.500 -1.496 -1.489
beta2:
Mean SD MC Error 95% HPD interval
------------------------------------------------------------------------------
-3.501 0.004 0.000 [-3.509, -3.492]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-3.509 -3.503 -3.501 -3.498 -3.492