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

# Data from code 6.2

set.seed(18472)

nobs <- 750 

x1_2 <- rbinom(nobs,size=1,prob=0.7)
x2 <- rnorm(nobs,0,1)
xb <- 1 - 1.5*x1_2 - 3.5*x2

exb <- exp(xb)
py <- rpois(nobs, exb)
pois <- data.frame(py, x1_2, x2)


poi <- glm(py ~ x1_2 + x2, family=poisson, data=pois)

 

Code 6.3 Bayesian Poisson model using R

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

library(MCMCpack)

mypoisL <- MCMCpoisson(py ~ x1_2 + x2,
                                              burnin = 5000,
                                              mcmc = 10000,
                                              data = pois)
 
summary(mypoisL)

plot(mypoisL)

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

 

Output on screen:

Iterations = 5001:15000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000

 

1. Empirical mean and standard deviation for each variable,

    plus standard error of the mean:

 

                         Mean                  SD          Naive SE        Time-series SE

(Intercept)        1.003        0.010538          1.054e-04               0.0003470

x1_2               -1.507        0.004699           4.699e-05              0.0001528

x2                   -3.501        0.004214           4.214e-05              0.0001395

 

2. Quantiles for each variable:

                              2.5%           25%          50%         75%         97.5%

(Intercept)            0.982        0.9961        1.003        1.011          1.023

x1_2                   -1.516       -1.5102       -1.507      -1.504         -1.498

x2                       -3.509       -3.5036       -3.501      -3.498         -3.493