top of page

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.16 Negative binomial model in Python using pymc3

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

import numpy as np
import pandas
import pylab as plt
import pymc3 as pm
from scipy.stats import uniform, binom, nbinom
import statsmodels.api as sm

â€‹

# Data
np.random.seed(141)                                                    # set seed to replicate example
nobs= 2500                                                                   # number of obs in model

â€‹

x1 = binom.rvs(1, 0.6, size=nobs)                                # categorical explanatory variable
x2 = uniform.rvs(size=nobs)                                        # real explanatory variable

theta = 0.303

xb = 1 + 2 * x1 - 1.5 * x2                                             # linear predictor

â€‹

exb = np.exp(xb)
nby = nbinom.rvs(exb, theta)

â€‹

df = pandas.DataFrame({'x1': x1, 'x2':x2,'nby': nby})   # 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
linp = beta0 + beta1 * x1 + beta2 * x2
mu = np.exp(linp)
mu2 = mu * (1 - theta)/theta                                     # compensate for difference between
# parametrizations from pymc3 and scipy

y_obs = pm.NegativeBinomial('y_obs', mu2, theta, observed=nby)

â€‹

# inference
start = pm.find_MAP()                                             # Find starting value by optimization
step = pm.NUTS()
trace = pm.sample(niter, step, start, progressbar=True)

â€‹

# print summary to screen
pm.summary(trace)

â€‹

# show graphical output
pm.traceplot(trace)
plt.show()

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

Output on screen:

â€‹

beta0:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

1.020            0.089            0.002            [0.843, 1.189]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

0.849          0.960          1.020          1.077          1.197

beta1:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

1.988            0.077            0.001            [1.843, 2.143]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

1.840          1.936          1.987          2.038          2.141

beta2:

Mean             SD               MC Error         95% HPD interval
------------------------------------------------------------------------------

-1.516           0.129            0.002            [-1.765, -1.259]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

-1.774         -1.601         -1.516         -1.429         -1.265

bottom of page