library(MASS)
library(car)
data(stackloss)


### R interface to JAGS:
library(R2jags)

# Create a data list with inputs for JAGS

n = nrow(stackloss)
## scale X such that X^TX has ones on the diagonal;
## scale divides by the standard deviation so we need
## to divide by the sqrt(n-1)
scaled.X = scale(as.matrix(stackloss[, -4]))/sqrt(n-1)
# are diagonal elements 1?
#
t(scaled.X) %*% scaled.X
data = list(Y = stackloss$stack.loss, X=scaled.X, p=ncol(scaled.X))
data$n = n   #check

data$scales = attr(scaled.X, "scaled:scale")*sqrt(n-1) # fix scale
data$Xbar = attr(scaled.X, "scaled:center")

# define a function that returns the Model

rr.model = function() {
  df <- 9
  shape <- df/2

  for (i in 1:n) {
    mu[i] <- alpha0 + inprod(X[i,], alpha)
    lambda[i] ~ dgamma(shape, shape)
    prec[i] <- phi*lambda[i]
    Y[i] ~ dnorm(mu[i], prec[i])
  }
  phi ~ dgamma(1.0E-6, 1.0E-6)
  alpha0 ~ dnorm(0, 1.0E-6)

  for (j in 1:p) {
    prec.beta[j] <- lambda.beta[j]*phi
    alpha[j] ~ dnorm(0, prec.beta[j])
    beta[j] <- alpha[j]/scales[j]
    lambda.beta[j] ~ dgamma(.5, .5)
  }


  beta0 <- alpha0 - inprod(beta[1:p], Xbar)

  sigma <- pow(phi, -.5)
}


parameters = c("beta0", "beta", "sigma","lambda.beta", "lambda")


bf.sim = jags(data, inits=NULL, par=parameters, model=rr.model,  n.iter=5000)

bf.bugs = as.mcmc(bf.sim$BUGSoutput$sims.matrix)  # create an MCMC object

plot(bf.sim)
summary(bf.sim)  # names of objects in bf.sim
bf.sim  # print gives summary
par(mfrow=c(1,1))
quantile(bf.bugs[,"beta0"], c(.025, .5, .975))
HPDinterval(bf.bugs)



par(mfrow=c(2,2))
hist(bf.bugs[,"beta0"], prob=T, xlab=expression(beta[0]),
     main="Posterior Distribution")
lines(density(bf.bugs[,"beta0"]))
#densplot(bf.bugs[,"beta0"])

hist(bf.bugs[,"beta[1]"], prob=T, xlab=expression(beta[1]),
     main="Posterior Distribution")
lines(density(bf.bugs[,"beta[1]"]))
#densplot(bf.bugs[,"beta[1]"])

hist(bf.bugs[,"beta[2]"], prob=T, xlab=expression(beta[2]),
     main="Posterior Distribution")
lines(density(bf.bugs[,"beta[2]"]))
#densplot(bf.bugs[,"beta[2]"])
hist(bf.bugs[,"beta[3]"], prob=T, xlab=expression(beta[3]),
     main="Posterior Distribution")
lines(density(bf.bugs[,"beta[3]"]))
#densplot(bf.bugs[,"beta[3]"])

pdf("lambda21.pdf")
hist(bf.bugs[,"lambda[21]"], prob=T, xlab=expression(lambda[39]),
     main="Posterior Distribution")
lines(density(bf.bugs[,"lambda[21]"]))
# add prior density
lines(seq(0.01, 2.5, length=100), dgamma(seq(0.01, 2.5, length=100), 9/2,rate=9/2))
dev.off()