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 5.3

set.seed(1056)                                                                # set seed to replicate example

nobs = 5000                                                                   # number of observations in model

x1 <- runif(nobs)                                                           # random uniform variable

xb <- 2 + 3*x1                                                               # linear predictor, xb

y <- rlnorm(nobs, xb, sdlog=1)                                     # create y as random lognormal variate

Code 5.4 Lognormal model in R using JAGS

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

require(R2jags)

X <- model.matrix(~ 1 + x1)
K <- ncol(X)
model_data <- list(Y = y, X = X, K = K, N = nobs,
                               Zeros = rep(0, nobs))


LNORM <-"
model{
    # Diffuse normal priors for predictors
    for (i in 1:K) { beta[i] ~ dnorm(0, 0.0001) }

    # Uniform prior for standard deviation
    tau <- pow(sigma, -2)                                              # precision
    sigma ~ dunif(0, 100)                                              # standard deviation

    # Likelihood
    for (i in 1:N){
        Y[i] ~ dlnorm(mu[i],tau)
        mu[i] <- eta[i]
        eta[i] <- inprod(beta[], X[i,])
    }
}"

inits <- function () { list(beta = rnorm(K, 0, 0.01)) }

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

LN <- jags(data = model_data,
                   inits = inits,
                   parameters = params,
                   model = textConnection(LNORM),
                   n.chains = 3,
                   n.iter = 5000,
                   n.thin = 1,
                   n.burnin = 2500)

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

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

# plot
source("https://raw.githubusercontent.com/astrobayes/BMAD/master/auxiliar_functions/CH-Figures.R")

out <- LN$BUGSoutput
MyBUGSHist(out,c(uNames("beta",K),"sigma"))

MyBUGSChains(out,c(uNames("beta",K),"sigma"))

Output on 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%             97.5%        Rhat       n.eff

beta[1]                 1.994       0.028          1.940               2.048      1.002       2200

beta[2]                 3.005      0.048           2.911               3.099      1.001       4700

sigma                   0.999      0.010          0.980                1.019      1.001       7500

deviance      49162.565      2.419  49159.819        49168.725      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 = 2.9 and DIC = 49165.5

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