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.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
==============================================================================