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.5 Multivariate normal model in R using JAGS for accessing the relationship between period, luminosity, and color in early-type contact binaries

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

library(R2jags)

â€‹

# Data

â€‹

# Prepare data for JAGS
nobs = nrow(PLC)                                                                  # number of data points
x1 <- PLC\$logP                                                                      # log period
x2 <- PLC\$V_I                                                                       # V-I color
y <- PLC\$M_V                                                                       # V magnitude
type <- as.numeric(PLC\$type)                                               # type NC/GC
X <- model.matrix(~ 1 + x1+x2)                                            # covariate matrix
K <- ncol(X)                                                                           # number of covariates per type

â€‹

jags_data <- list(Y = y,
X = X,
K = K,
type = type,
N = nobs)

â€‹

# Fit
NORM <-"model{
# Shared hyperprior
tau0 ~ dgamma(0.001,0.001)
mu0 ~ dnorm(0,1e-3)

â€‹

# Diffuse normal priors for predictors
for(j in 1:2){
for (i in 1:K) {
beta[i,j] ~ dnorm(mu0, tau0)
}
}

â€‹

# Uniform prior for standard deviation
for(i in 1:2) {
tau[i] <- pow(sigma[i], -2)                         # precision
sigma[i] ~ dgamma(1e-3, 1e-3)                 # standard deviation
}

â€‹

# Likelihood function
for (i in 1:N){
Y[i]~dnorm(mu[i],tau[type[i]])
mu[i] <- eta[i]
eta[i] <- beta[1, type[i]] * X[i, 1] + beta[2, type[i]] * X[i, 2] +
beta[3, type[i]] * X[i, 3]
}
}"

â€‹

# Determine initial values
inits <- function () {
list(beta = matrix(rnorm(6,0, 0.01),ncol=2))
}

â€‹

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

â€‹

# Fit
jagsfit <- jags(data = jags_data,
inits = inits,
parameters = params,
model = textConnection(NORM),
n.chains = 3,
n.iter = 5000,
n.thin = 1,
n.burnin = 2500)

â€‹

# Output
print(jagsfit,justify = "left", digits=2)

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

Output screen:

â€‹

â€‹

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

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

n.sims = 7500 iterations saved

mu.vect     sd.vect     2.5%      25%      50%       75%     97.5%      Rhat       n.eff

beta[1,1]    -1.01         0.27    -1.55      -1.18     -1.00      -0.83      -0.47             1       7500

beta[2,1]    -3.31         0.97    -5.21      -3.94     -3.30      -2.67      -1.41             1       7500

beta[3,1]     7.24         1.29     4.68        6.39      7.26       8.10        9.73             1       7500

beta[1,2]     -0.41        0.16    -0.73      -0.52     -0.42      -0.31      -0.10             1       7500

beta[2,2]     -3.19        0.58    -4.31      -3.58     -3.19      -2.81       -2.01            1       7500

beta[3,2]      8.48        0.82     6.85       7.94       8.48       9.03       10.08            1       6900

sigma[1]      0.62        0.09     0.47      0.55        0.61       0.67         0.82            1       1500

sigma[2]      0.43        0.05     0.34      0.39        0.42       0.46         0.55            1       4300

deviance    91.96       4.34      85.5     88.78     91.27     94.39     102.15            1       7300

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 = 9.4 and DIC = 101.4

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

bottom of page