Multivariate 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 6.6 Multivariate Poisson model in Python
===================================================

import numpy as np
from scipy.stats import norm, binom, poisson
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 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

my_data = {}
my_data['x1_2'] = x1_2
my_data['x2'] = x2
my_data['py'] = py

# build model
myp = smf.glm('py ~ x1_2 + x2', data=my_data, family=sm.families.Poisson())

# find parameter values
res = myp.fit()

# Output
print(res.summary())

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

Output on screen:

Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: py No. Observations: 750
Model: GLM Df Residuals: 747
Model Family: Poisson Df Model: 2
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -1171.0
Date:    Sat, 24 Dec 2016 Deviance: 556.26
Time: 01:50:37 Pearson chi2: 643.
No. Iterations: 12
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------------------------------------------------------------
Intercept    1.0021 0.012 83.174 0.000 0.979 1.026
x1_2 -1.5003 0.006   -252.682 0.000 -1.512 -1.489
x2 -3.5004 0.004 -796.736 0.000    -3.509 -3.492
==============================================================================

GET SOURCE

HSI

HSI