#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