## Setup
require(R2jags)
require(R2OpenBUGS)
bData <- with(swiss, list(x1=Agriculture, x2=Examination, x3=Education, x4=Catholic, x5=Infant.Mortality, y=Fertility, n=nrow(swiss)))
model <- function() {
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta[1] + beta[2]*x1[i] + beta[3]*x2[i] + beta[4]*x3[i] + beta[5]*x4[i] + beta[6]*x5[i]
  }
  for (j in 1:6) {
    beta[j] ~ dnorm(0, 0.0001)
  }
  tau ~ dgamma(0.001, 0.001)
}

## Slide 4-5
init <- function() list(beta=rnorm(6, sd=100), tau=runif(1,0,100)^(-2))
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=100, n.chain=3, n.burn=0, DIC=FALSE)
traceplot(fit, mfrow=c(4,2), ask=FALSE)
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=500, n.chain=3, n.burn=0, DIC=FALSE)
traceplot(fit, mfrow=c(4,2), ask=FALSE)
## Fully specified initial values
init <- list(list(beta=rep(100,6), tau=1), list(beta=rep(25,6), tau=1), list(beta=rep(-50,6), tau=1))

## Slide 9
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=1000, n.chain=3, n.burn=0, DIC=FALSE)
mcfit <- as.mcmc.list(fit)
gelman.diag(mcfit, autoburnin=FALSE)
## Similar to:
A <- fit$sims.array[,,1]
B <- sd(as.numeric(A))
W <- mean(apply(A, 2, sd))
B/W

## Slide 15
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=10000, n.chain=3, n.burn=5000, DIC=FALSE)
mcfit <- as.mcmc.list(fit)
acf(fit$sims.array[,2,1], ylim=c(0,1))
acfplot(mcfit, outer=TRUE, aspect="fill")

## Autocorrelation with thinning
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=10000, n.chain=3, n.burn=5000, n.thin=10, DIC=FALSE)
mcfit <- as.mcmc.list(fit)
acf(fit$sims.array[,2,1], ylim=c(0,1))
acfplot(mcfit, outer=TRUE, aspect="fill")

## Slide 22
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=10000, n.chain=3, n.burn=5000, n.thin=1, DIC=FALSE)
mcfit <- as.mcmc.list(fit)
effectiveSize(mcfit)
summary(mcfit)
apply(fit$sims.matrix, 2, sd)
batchSE(mcfit)
n.eff <- (apply(fit$sims.matrix, 2, sd)/batchSE(mcfit))^2

## Slide 24-25
fit <- bugs(model=model, init=init, param=c("beta", "tau"), data=bData, n.iter=10000, n.chain=3, n.burn=5000, n.thin=1, DIC=FALSE)
fit
plot(fit)