HSI
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 10.14
require(R2jags)
require(jagstools)
# Data
path_to_data = "https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p7/GCs.csv"
​
# Read data
GC_dat = read.csv(file=path_to_data,header = T,dec=".")
​
# 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)
​
​
​
​
​
​
Code 10.16 NB-P model in R using JAGS, for modeling the relationship between globular cluster population and host galaxy visual magnitude
=================================================================================
# Fit
model.NBP <- "model{
# Priors for regression coefficients
# Diffuse normal priors betas
for (i in 1:K) { beta[i] ~ dnorm(0, 1e-5)}
​
# Prior for size
theta ~ dgamma(1e-3,1e-3)
# Uniform prior for Q
Q ~ dunif(0,3)
​
# Likelihood
for (i in 1:N){
eta[i]<-inprod(beta[], X[i,])
mu[i] <- exp(eta[i])
theta_eff[i]<- theta*(mu[i]^Q)
p[i]<-theta_eff[i]/(theta_eff[i]+mu[i])
Y[i]~dnegbin(p[i],theta_eff[i])
​
# Discrepancy
expY[i] <- mu[i] # mean
varY[i] <- mu[i] + pow(mu[i],2-Q)/theta #variance
PRes[i] <- ((Y[i] - expY[i])/sqrt(varY[i]))^2
}
Dispersion <- sum(PRes)/(N-4)#
}"
​
# Set initial values
inits <- function () {
list(
beta = rnorm(K, 0, 0.1),
theta = runif(1,0.1,5),
Q = runif(1,0,1)
)
}
​
# Identify parameters
params <- c("Q","beta","theta","Dispersion")
​
# Start JAGS
NBP_fit <- jags(data = JAGS_data,
inits = inits,
parameters = params,
model = textConnection(model.NBP),
n.thin = 1,
n.chains = 3,
n.burnin = 5000,
n.iter = 20000)
​
# Output
# Plot posteriors
MyBUGSHist(out,c("Dispersion",uNames("beta",K),"theta"))
​
# Screen output
print(NBP_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 20000 iterations (first 5000 discarded)
n.sims = 45000 iterations saved
mu.vect sd.vect 2.5% 97.5% Rhat n.eff
Dispersion 1.925 0.208 1.547 2.355 1.002 1700
Q 0.018 0.015 0.001 0.056 1.001 5500
beta[1] -11.797 0.328 -12.454 -11.172 1.012 180
beta[2] -0.883 0.016 -0.916 -0.852 1.012 180
theta 0.996 0.104 0.775 1.184 1.001 32000
deviance 5193.038 2.920 5189.225 5200.284 1.001 45000
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 = 4.3 and DIC = 5197.3
DIC is an estimate of expected predictive error (lower deviance is better).