Bayesian random intercept Poisson model in Python

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 8.13 Bayesian random intercept Poisson model in Python

===============================================

import numpy as np
import pymc3 as pm

from scipy.stats import norm, uniform, poisson

# Data
np.random.seed(1656) # set seed to replicate example
N = 2000 # number of obs in model
NGroups = 10

x1 = uniform.rvs(size=N)
x2 = uniform.rvs(size=N)

Groups = np.array([200 * [i] for i in range(NGroups)]).flatten()

a = norm.rvs(loc=0, scale=0.5, size=NGroups)
eta = 1 + 0.2 * x1 - 0.75 * x2 + a[list(Groups)]

mu = np.exp(eta)

y = poisson.rvs(mu, size=N)

with pm.Model() as model:
# Define priors
sigma_a = pm.Uniform('sigma_a', 0, 100)
beta1 = pm.Normal('beta1', 0, sd=100)
beta2 = pm.Normal('beta2', 0, sd=100)
beta3 = pm.Normal('beta3', 0, sd=100)

# priors for random intercept (RI) parameters
a_param = pm.Normal('a_param',
np.repeat(0, NGroups), # mean
sd=np.repeat(sigma_a, NGroups), # standard deviation
shape=NGroups) # number of RI parameters

eta = beta1 + beta2*x1 + beta3*x2 + a_param[Groups]

# Define likelihood
y = pm.Poisson('y', mu=np.exp(eta), observed=y)

# Fit
start = pm.find_MAP() # Find starting value by optimization
step = pm.NUTS(scaling=start)    # Initiate sampling

trace = pm.sample(20000, step, start=start, progressbar=False)

# Print summary to screen
pm.summary(trace)

===============================================

GET SOURCE

Output on screen:

Applied interval-transform to sigma_a and added transformed sigma_a_interval to model.

sigma_a_interval:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-5.202 0.268 0.003 [-5.716, -4.660]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-5.684 -5.388 -5.223 -5.035 -4.612

beta1:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.966 0.187 0.006 [0.592, 1.331]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.579 0.852 0.969 1.088 1.322

beta2:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.214 0.052 0.001 [0.115, 0.315]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.113 0.179 0.214 0.249 0.315

beta3:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
-0.742 0.053 0.000 [-0.846, -0.638]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
-0.847 -0.779 -0.742 -0.707 -0.639

a_param:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.373 0.187 0.006 [0.008, 0.755]
0.346 0.187 0.006 [-0.015, 0.732]
0.147 0.188 0.006 [-0.219, 0.531]
-1.117 0.199 0.006 [-1.518, -0.726]
0.152 0.187 0.006 [-0.206, 0.542]
-0.268 0.189 0.006 [-0.634, 0.125]
0.192 0.188 0.006 [-0.173, 0.580]
0.471 0.186 0.006 [0.106, 0.846]
0.122 0.188 0.006 [-0.241, 0.509]
-0.446 0.192 0.006 [-0.827, -0.063]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.014 0.252 0.371 0.489 0.761
-0.012 0.225 0.344 0.462 0.736
-0.216 0.025 0.144 0.262 0.537
-1.514 -1.242 -1.118 -0.988 -0.721
-0.208 0.031 0.149 0.268 0.540
-0.634 -0.389 -0.270 -0.150 0.125
-0.173 0.070 0.190 0.307 0.581
0.116 0.352 0.468 0.585 0.860
-0.238 0.000 0.121 0.238 0.514
-0.817 -0.572 -0.447 -0.326 -0.052

sigma_a:

Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.568 0.164 0.002 [0.310, 0.900]

Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.339 0.455 0.536 0.646 0.983

HSI

HSI