# WINBUGS in R for SLR
library("R2WinBUGS")

# comment out the lines below for Windows
WINE <- "/Applications/Darwine/Wine.bundle/Contents/bin/wine"
WINEPATH <- "/Applications/Darwine/Wine.bundle/Contents/bin/winepath"
BUGS.DIR = "/users/clyde/.wine/drive_c/Program Files/WinBUGS14/"


year =seq(1971, 1980, by=1)
y = c(41,52,18.7, 55, 40, 29.2,51,17.6, 46.6, 57)
x = c(23.9, 43.3, 36.3, 40.6, 57, 52.5, 46.1, 142, 112.6, 23.7)
x = x - mean(x)

rain = data.frame(year, y, x)
plot(rain)
summary(rain)
rain.lm = lm(y ~ x, data=rain)
summary(rain.lm)

data = list(n=length(y), y=y, x=x)

model = function() {
  df.5 <- 5*.5
  for ( i in 1:n) {
    lambda[i] ~ dgamma(df.5, .5)
    prec[i] <- lambda[i]*phi
    y[i] ~ dnorm(mu[i], prec[i])
    mu[i] <- alpha + beta*x[i]
  }
  prec.alpha <- .000001*phi
  alpha ~ dnorm(0, prec.alpha)
  prec.beta <- .000001*phi
  beta ~ dnorm(0, prec.beta)
  phi ~ dgamma(.000001, .000001)
  mu.pred <- alpha + beta*80
  lambda.pred ~ dgamma(df.5, .5)
  prec.pred <- lambda.pred*phi
  Y.pred ~ dnorm(mu.pred, prec.pred)
}


inits = function() {
  alpha = coef(rain.lm)[1]
  beta = coef(rain.lm)[2]
  phi = (length(y) - 2)/sum((rain.lm$residuals^2))
  lambda = rep(1, length(y))
  lambda.pred = 1

return(list(alpha=alpha, beta=beta, phi=phi, lambda=lambda,
            lambda.pred=lambda.pred))
}

parameters = c("alpha", "beta", "phi", "Y.pred", "lambda")
write.model(model, "model.txt")
model.file = "model.txt"
 
slr.sim = bugs(data,  inits=inits, parameters=parameters, model.file, n.chains=1, n.iter=5000, bugs.dir=BUGS.DIR, WINE=WINE, WINEPATH=WINEPATH, debug=T, DIC=F)

par(mfrow=c(2,2))

plot(density(slr.sim$sims.matrix[,"alpha"]), xlab=expression(alpha), ylab="Density", main=expression(paste("Posterior Distribution of ", alpha)))

plot(density(slr.sim$sims.matrix[,"beta"]), xlab=expression(beta), ylab="Density", main=expression(paste("Posterior Distribution of ", beta)))

plot(density(slr.sim$sims.matrix[,"phi"]), xlab=expression(phi), ylab="Density", main=expression(paste("Posterior Distribution of ", phi)))

plot(density(slr.sim$sims.matrix[,"Y.pred"]), xlab="Y", ylab="Density",
     main="Predictive Distribution of Y | X = 18")