Bayesian Poisson model using R

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

# 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)

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

Code 6.3

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

HSI