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 10.14 Poisson model, in R using JAGS, for modeling the relation between globular  clusters population and host galaxy visual magnitude

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

require(R2jags)
require(jagstools)

# Data
path_to_data =

# Prepare data to JAGS
N <- nrow(GC_dat)
x <- GC_dat\$MV_T
y <- GC_dat\$N_GC
X <- model.matrix(~ x, data=GC_dat)
K = ncol(X)

JAGS_data <- list(
Y = y,
X = X,
N = N,
K = K)

# Fit
model.pois <- "model{
# Diffuse normal priors betas
for (i in 1:K) { beta[i] ~ dnorm(0, 1e-5)}

for (i in 1:N){
# Likelihood
eta[i]<-inprod(beta[], X[i,])
mu[i] <- exp(eta[i])
Y[i]~dpois(mu[i])

# Discrepancy
expY[i] <- mu[i] # mean
varY[i] <- mu[i] # variance
PRes[i] <- ((Y[i] - expY[i])/sqrt(varY[i]))^2
}

Dispersion <- sum(PRes)/(N-2)
}"

# Define initial values
inits <- function () {
list(beta = rnorm(K, 0, 0.1))
}

# Identify parameters
params <- c("beta","Dispersion")

# Start JAGS
pois_fit <- jags(data = JAGS_data ,
inits = inits,
parameters = params,
model = textConnection(model.pois),
n.thin = 1,
n.chains = 3,
n.burnin = 3500,
n.iter = 7000)

# Output
print(pois_fit , intervals=c(0.025, 0.975), digits=3)

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

â€‹

Output on screen:

â€‹

Inference for Bugs model at "3", fit using jags,

3 chains, each with 7000 iterations (first 3500 discarded)

n.sims = 10500 iterations saved

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

Dispersion        1080.232        0.205          1079.881         1080.603     1.018      10000

beta[1]                 -11.910        0.038            -11.966            -11.855      1.016      10000

beta[2]                   -0.918        0.002              -0.920              -0.915      1.016      10000

deviance       497171.835     935.473     495628.137     498733.429      1.005      10000

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 = 437637.2 and DIC = 934809.1

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

bottom of page