bodyfat = read.table("bodyfat.txt", header=T)

### WinBugs:
library(R2WinBUGS)

# Create a data list with inputs for WinBugs

bf.data = list(Y = bodyfat$Bodyfat, X=bodyfat$Abdomen)
bf.data$n = length(bf.data$Y)
bf.data$Xbar = mean(bf.data$X)

# define a function that returns the Model for Winbugs

rr.model = function() {
  df <- 9
  shape <- df/2
  
  for (i in 1:n) {
    mu[i] <- alpha0 + alpha1*(X[i] - Xbar)
    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)
  alpha1 ~ dnorm(0,1.0E-6)
  beta0 <- alpha0 - alpha1*Xbar
  beta1 <- alpha1
  sigma <- pow(phi, -.5)
  mu34 <- beta0 + beta1*2.54*34
  y34 ~ dt(mu34,phi, df) 
}

# create a function that provides intial values for WinBUGS
rr.inits = function() {
    bf.lm <- lm(bf.data$Y ~ I(bf.data$X - bf.data$Xbar))
    coefs = coef(bf.lm)
    alpha1=coefs[2]
    alpha0 = coefs[1]
    phi = (1/summary(bf.lm)$sigma)^2
 lambda = rep(1, bf.data$n)
return(list(alpha0=alpha0, alpha1 = alpha1, phi=phi, lambda=lambda))
}

# a list of all the parameters to save

parameters = c("beta0", "beta1", "sigma", "mu34", "y34", "lambda[39]")

#uncomment for MAC OSX  
WINE <- "/Applications/Darwine/Wine.bundle/Contents/bin/wine"
WINEPATH <- "/Applications/Darwine/Wine.bundle/Contents/bin/winepath"

#uncomment for Linux or MAC OSX using WINE
BUGS.DIR <- "/users/clyde/.wine/drive_c/Program Files/WinBUGS14/"

write.model(rr.model, "rr-model.txt")
model.file = "rr-model.txt"

# For Windows
#bf.sim = bugs(bf.data, rr.inits, parameters, model.file, n.chains=2, n.iter=5000, debug=T, DIC=F)

# For LINUX
#bf.sim = bugs(bf.data, rr.inits, parameters, model.file, n.chains=2, n.iter=5000,  bugs.dir=BUGS.DIR, debug=T, DIC=F)

# For MAC OSX
bf.sim = bugs(bf.data, rr.inits, parameters, model.file, n.chains=2, n.iter=5000,  bugs.dir=BUGS.DIR, WINE=WINE, WINEPATH=WINEPATH, debug=T, DIC=F)


plot(bf.sim)
summary(bf.sim)  # names of objects in bf.sim
bf.sim$summary
par(mfrow=c(1,1))
quantile(bf.sim$sims.matrix[,"beta1"], c(.025, .5, .975))
HPDinterval(as.mcmc(bf.sim$sims.matrix[,"mu34"]))
HPDinterval(as.mcmc(bf.sim$sims.matrix[,"y34"]))

hist(bf.sim$sims.matrix[,"beta1"], prob=T, xlab=expression(beta[1]),
     main="Posterior Distribution")
lines(density(bf.sim$sims.matrix[,"beta1"]))

par(mfrow=c(1,2))
hist(bf.sim$sims.matrix[,"mu34"], prob=T, xlab=expression(mu),
     main="Posterior Distribution of Expected Bodyfat\n for Men with 34 inch Circumference ")
lines(density(bf.sim$sims.matrix[,"mu34"]))

hist(bf.sim$sims.matrix[,"y34"], prob=T, xlab="Abdominal Circumference",
     main="Predictive Distribution of Bodyfat\n for Men with 34 inch Circumference ")
lines(density(bf.sim$sims.matrix[,"y34"]))

postscript("lambda39.ps")
hist(bf.sim$sims.matrix[,"lambda[39]"], prob=T, xlab=expression(lambda[39]),
     main="Posterior Distribution")
lines(density(bf.sim$sims.matrix[,"lambda[39]"]))
dev.off()

# compare to Normal results
bodyfat.lm = lm(Bodyfat ~ I(Abdomen - 2.54*34), data=bodyfat)
coef = summary(bodyfat.lm)$coef

c(coef[2,1] - coef[2,2]*qt(.975,250), coef[2,1] + coef[2,2]*qt(.975,250))
# interval for mu
c(coef[1,1] - coef[1,2]*qt(.975,250), coef[1,1] + coef[1,2]*qt(.975,250))
# interval for y
std.pred = sqrt(coef[1,2]^2 + (summary(bodyfat.lm)$sigma)^2)
c(coef[1,1] - std.pred*qt(.975,250), coef[1,1] + std.pred*qt(.975,250))


# compare to Normal results w/out case 39
bodyfat.lm = lm(Bodyfat ~ I(Abdomen - 2.54*34), subset=c(-39), data=bodyfat)
coef = summary(bodyfat.lm)$coef

c(coef[2,1] - coef[2,2]*qt(.975,250), coef[2,1] + coef[2,2]*qt(.975,250))
# interval for mu
c(coef[1,1] - coef[1,2]*qt(.975,250), coef[1,1] + coef[1,2]*qt(.975,250))
# interval for y
std.pred = sqrt(coef[1,2]^2 + (summary(bodyfat.lm)$sigma)^2)
c(coef[1,1] - std.pred*qt(.975,250), coef[1,1] + std.pred*qt(.975,250))


#############Gibbs Code
bf.data$nu = 9
Msims = 500
theta = rr.inits()
thetasim = matrix(NA, nrow = Msims, ncol= length(unlist(theta)))

for (m in 1:Msims) {
  theta = gen.alpha0(theta, bf.data)
  theta = gen.alpha1(theta,bf.data)
  theta = gen.phi(theta,bf.data)
  theta = gen.lambda(theta,bf.data)
  thetasim[m,] = unlist(theta)
# print( thetasim[m,1:3])
}

gen.alpha0 = function(theta, data) {
  attach(theta)
  attach(data)
  theta$alpha0 = rnorm(1, sum(lambda*(Y - alpha1*(X -Xbar)))/sum(lambda),
                          sqrt(1/(sum(lambda)*phi)))
  detach(theta)
  detach(data)

  return(theta)

}

gen.alpha1 = function(theta, data) {
  attach(theta)
  attach(data)
  theta$alpha1 = rnorm(1, sum(lambda*(Y-alpha0)*(X - Xbar))/
                          sum(lambda*(X-Xbar)^2),
                          sqrt(1/(sum(lambda*(X- Xbar)^2)*phi)))
  detach(theta)
  detach(data)
  return(theta)

}

gen.lambda = function(theta, data) {
attach(theta)
attach(data)
theta$lambda = rgamma(n, .5*(nu + 1),
                      rate=.5*(phi*(Y - alpha0 -alpha1*(X - Xbar))^2 + nu))
  detach(theta)
  detach(data)
return(theta)
  
}

gen.phi = function(theta, data) {
  attach(theta)
  attach(data)
  theta$phi = rgamma(1, .5*n,
                    rate=.5*(sum(lambda*(Y - alpha0 -alpha1*(X-Xbar))^2)))
  detach(theta)
  detach(data)
  return(theta)
}