#Highland Statistics Ltd.
#File to test whether JAGS runs properly on
#your computer.
#1. Install the JAGS binary file
#Go to: http://sourceforge.net/projects/mcmc-jags/files/JAGS/
#Use the executable for windows...or the .dmg file for Mac
#Mac users: One of the instructors (Alain Zuur) is using a Mac with Yosemite...
# and he managed to install JAGS. He had to ensure that both JAGS and
# R were created with the same compiler (and updating quartz may be a good thing too).
# He also had to go for the second latest version.
# In case of trouble (and you will know that within a minute) just google the error message.
#
#2. Install and load the package R2jags
install.packages("R2jags", dependencies = TRUE)
library(R2jags)
#3. RUN ALL CODE BELOW. IF YOU DO NOT GET A GRAPH WITH SOME
#RESULTS, FIGURE OUT WHAT YOU HAVE DONE WRONG. IT IS
#UNLIKELY THAT WE CAN SORT IT OUT DURING THE COURSE
set.seed(1)
X <- runif(1000)
alpha <- 1
beta <- 2
eps <- rnorm(1000, mean = 0, sd = 0.5)
Y <- alpha + beta * X + eps
M1 <- lm(Y ~ X)
summary(M1)
#And now the JAGS bit:
JAGS.data <- list(Y = Y,
X = X,
N = 1000)
sink("JAGStest.txt")
cat("
model{
#Diffuse Prior
b[1] ~ dnorm(0, 0.001)
b[2] ~ dnorm(0, 0.001)
sigma ~ dunif(0, 100)
tau <- 1 / (sigma * sigma)
for (i in 1:N){
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- b[1] + b[2] * X[i]
}
}
",fill = TRUE)
sink()
#Set the initial values for the betas and sigma
inits <- function () {
list(
b = rnorm(2, 0, 0.1),
sigma = rlnorm(1)) }
#Parameters to estimate
params <- c("b", "sigma")
#MCMC settings
nc <- 3 #Number of chains
ni <- 500 #Number of draws from posterior (for each chain)
nb <- 50 #Number of draws to discard as burn-in
nt <- 5 #Thinning rate
Test <- jags(data = JAGS.data,
inits = inits,
parameters = params,
model = "JAGStest.txt",
n.thin = nt,
n.chains = nc,
n.burnin = nb)
print(Test, digits = 3)
# You should get output like this (though your results will differ slightly):
# Inference for Bugs model at "JAGStest.txt", fit using jags,
# 3 chains, each with 2000 iterations (first 50 discarded), n.thin = 5
# n.sims = 1170 iterations saved
# mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
# b[1] 0.983 0.033 0.919 0.961 0.983 1.005 1.043 1.000 1200
# b[2] 2.010 0.057 1.903 1.970 2.010 2.047 2.121 1.001 1200
# sigma 0.518 0.011 0.495 0.510 0.517 0.525 0.540 1.000 1200
# deviance 1519.358 2.302 1516.559 1517.641 1518.766 1520.494 1525.519 1.003 1200
#
# 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.7 and DIC = 1522.0
# DIC is an estimate of expected predictive error (lower deviance is better).
plot(Test)
# This should make a graph.
# It doesn't matter how it looks like. A graph is good, no graph is bad