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.18 Synthetic data for generalized Poisson
====================================================

require(MASS)
require(R2jags)

â€‹

set.seed(160)
nobs <- 1000

â€‹

x1 <- runif(nobs)
xb <- 1 + 3.5*x1
exb <- exp(xb)

delta <- -0.3
gpy <- c()

â€‹

for (i in 1:nobs){
gpy[i] <- rgp(1, mu=(1-delta)*exb[i], delta = delta)
}

â€‹

gpdata <- data.frame(gpy, x1)

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

â€‹

â€‹

â€‹

Code 6.19 Bayesian generalized Poisson using JAGS

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

X <- model.matrix(~ x1, data = gpdata)
K <- ncol(X)

â€‹

model.data <- list( Y = gpdata\$gpy,                                            # response
X = X,                                                            # covariates
N = nrow(gpdata),                                         # sample size
K = K,                                                            # number of betas
Zeros = rep(0, nrow(gpdata)))

â€‹

sink("GP1reg.txt")

â€‹

cat("
model{
# Priors beta
for (i in 1:K) { beta[i] ~ dnorm(0, 0.0001)}

â€‹

# Prior for delta parameter of GP distribution
delta ~ dunif(-1, 1)

â€‹

C <- 10000
for (i in 1:N){
Zeros[i] ~ dpois(Zeros.mean[i])
Zeros.mean[i] <- -L[i] + C
l1[i] <- log(mu[i])
l2[i] <- (Y[i] - 1) * log(mu[i] + delta * Y[i])
l3[i] <- -mu[i] - delta * Y[i]
l4[i] <- -loggam(Y[i] + 1)
L[i] <- l1[i] + l2[i] + l3[i] + l4[i]
mu[i] <- (1 - delta)*exp(eta[i])
eta[i] <- inprod(beta[], X[i,])
}

â€‹

# Discrepancy measures: mean, variance, Pearson residuals
for (i in 1:N){
ExpY[i] <- mu[i] / (1 - delta)
VarY[i] <- mu[i] / ((1 - delta)^3)
Pres[i] <- (Y[i] - ExpY[i]) / sqrt(VarY[i])
} }
"
,fill = TRUE)

â€‹

sink()

â€‹

inits <- function () {
list(beta = rnorm(ncol(X), 0, 0.1),
delta = 0)}

â€‹

params <- c("beta", "delta")

â€‹

GP1 <- jags(data = model.data,
inits = inits,
parameters = params,
model = "GP1reg.txt",
n.thin = 1,
n.chains = 3,
n.burnin = 4000,
n.iter = 5000)

â€‹

print(GP1, intervals=c(0.025, 0.975), digits=3)

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

Anchor 1

Output on screen:

â€‹

Inference for Bugs model at "GP1reg.txt", fit using jags,

3 chains, each with 5000 iterations (first 4000 discarded)

n.sims = 3000 iterations saved

mu.vect    sd.vect                     2.5%                 97.5%       Rhat         n.eff

beta[1]                       1.029       0.017                    0.993                 1.062       1.017           130

beta[2]                       3.475       0.022                    3.432                  3.521      1.016           140

delta                         -0.287       0.028                   -0.344                 -0.233      1.002         1900

deviance      20005178.615       2.442      20005175.824     20005184.992      1.000               1

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 3.0 and DIC = 20005181.6

DIC is an estimate of expected predictive error (lower deviance is better).

bottom of page