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 4.8 Synthetic normal data in R with errors in measurements

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

require(R2jags)


# Data
set.seed(1056)                                                        # set seed to replicate example
nobs = 1000                                                           # number of obs in model

sdobsx <- 1.25
truex <- rnorm(nobs, 02.5)                                  # normal variable
errx <- rnorm(nobs, 0, sdobsx)
obsx <- truex + errx

 

beta1 <- -4
beta2 <- 7
sdy <- 1.25
sdobsy <- 2.5

erry <- rnorm(nobs, 0, sdobsy)
truey <- rnorm(nobs, beta1 + beta2 * truex, sdy)
obsy <- truey + erry

Code 4.9 Normal linear model in R using JAGS and ignoring errors in measurements

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

K <-  2
model.data <- list(obsy = obsy,
                              obsx = obsx,
                              K = K,
                              N = nobs)

NORM <- " 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 function
    for (i in 1:N){
        obsy[i] ~ dnorm(mu[i], tau)
        mu[i] <- eta[i]
        eta[i] <- beta[1] + beta[2] * obsx[i]
    }
}"

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

 

# Parameters to display and save
params <- c("beta", "sigma")

 

normefit <- jags(data = model.data,
                            inits = inits,
                            parameters = params,
                            model = textConnection(NORM),
                            n.chains = 3,
                            n.iter = 10000,
                            n.thin = 1,
                            n.burnin = 5000)

 

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

 

Output on screen:

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

    3 chains, each with 10000 iterations (first 5000 discarded) .

    n.sims = 15000 iterations saved

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

beta[1]             -4.077          0.268          -4.607          -3.552       1.001      15000

beta[2]              5.637          0.094            5.453           5.822       1.001      15000

sigma                8.555          0.194            8.185           8.942       1.001      15000

deviance     7128.687          2.467      7125.886      7134.988      1.002        3100

 

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 = 3.0 and DIC = 7131.7

 

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