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.1 Normal linear model in R using JAGS for accessing the relationship between central black hole mass and bulge velocity dispersion

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

require(R2jags)


# Data
path_to_data = "https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p1/M_sigma.csv"

# Read data
MS <- read.csv(path_to_data,header = T)

# Identify variables
N <- nrow(MS)                                                                  # number of data points
obsx <- MS$obsx                                                               # log black hole mass
errx <- MS$errx                                                                 # error on log black hole mass
obsy <- MS$obsy                                                               # log velocity dispersion
erry <- MS$erry                                                                 # error on log velocity dispersion

# Prepare data to JAGS
MS_data <- list(
  obsx = obsx,
  obsy = obsy,
  errx = errx,
  erry = erry,
  N = N
)

# Fit
NORM_errors <-"model{
    # Diffuse normal priors for predictors
    alpha ~ dnorm(0,1e-3)
    beta ~ dnorm(0,1e-3)

    # Gamma prior for scatter
    tau ~ dgamma(1e-3,1e-3) # precision
    epsilon <- 1/sqrt(tau) # intrinsic scatter

    # Diffuse normal priors for true x
    for (i in 1:N){ x[i] ~ dnorm(0,1e-3) }

    for (i in 1:N){
        obsx[i] ~ dnorm(x[i], pow(errx[i], -2))
        obsy[i] ~ dnorm(y[i], pow(erry[i], -2)) # likelihood function
        y[i] ~ dnorm(mu[i], tau)
        mu[i] <- alpha + beta*x[i] # linear predictor
    }
}"

# Define initial values
inits <- function () {
  list(alpha = runif(1,0,10),
       beta = runif(1,0,10))
}

# Identify parameters
params0 <- c("alpha","beta", "epsilon")

# Fit
NORM_fit <- jags(data = MS_data,
                               inits = inits,
                               parameters = params0,
                               model = textConnection(NORM_errors),
                               n.chains = 3,
                               n.iter = 50000,
                               n.thin = 10,
                               n.burnin = 30000
                               )


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

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

Output on screen:

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

   3 chains, each with 50000 iterations (first 30000 discarded), n.thin = 10

   n.sims = 6000 iterations saved

 

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

alpha               8.35         0.06          8.24        8.31           8.35         8.39          8.46       1.00       5300

beta                 4.44         0.75          3.78        4.23           4.46         4.68           5.10      1.04       3000

epsilon            0.27         0.06          0.16        0.23           0.27          0.31          0.39      1.00       2800

deviance    -225.79       15.33     -251.44   -235.79     -226.50     -216.91     -197.08      1.00      6000

 

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 = 117.6 and DIC = -108.2

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

© 2017 by Emille E. O. Ishida