set.seed(1234)

library(mcmc)

NITER <- 70000
start <- 50001
BatchLength <- 1

Y <- read.table("http://www.biostat.umn.edu/~sudiptob/pubh5485/Ytoughnut.txt", header=F)
Y <- Y[,1]

X <- as.matrix(read.table("http://www.biostat.umn.edu/~sudiptob/pubh5485/Xtoughnut.txt", header=F))

##Only to compare with standard linear model:

simple.lm <- lm(Y ~ X[,2]+X[,3]+X[,4]+X[,5]+X[,6]+X[,7])
simple.lm.summary <- summary(simple.lm)

no.parameters = ncol(X)+1

##Defining the target function. Here theta = (beta, log(sigmasq))

target <- function (theta, a, b)
{
  beta <- theta[1:ncol(X)]
  sigmasq <- exp(theta[length(theta)])

  out <- -(a+1+length(Y)/2.0) * log(sigmasq) - (1/sigmasq)*(b + 0.5*t(Y-X%*%beta)%*%(Y-X%*%beta))

  return(out)
}

myscale = vcov(simple.lm)
myscale = diag(c(diag(myscale),4.75e-6))
##myscale <- diag(c(rep(1.5e-6,times=ncol(X)),2.5e-5))
myscale <- chol(myscale)
inits <- rep(0, times=no.parameters)

metrop.out <- metrop(target, initial=inits, nbatch=NITER, blen=BatchLength, scale=myscale, a=2.0, b=0.000001)

myqnt <- c(0.50, 0.025, 0.975)
beta.qnt <- matrix(nrow=ncol(X),ncol=length(myqnt))
sigmasq.qnt <- rep(0,times=length(myqnt))

for (i in 1:nrow(beta.qnt)) {
  beta.qnt[i,] <- quantile(metrop.out$batch[start:nrow(metrop.out$batch),i],myqnt)
}

sigmasq.qnt <- quantile(exp(metrop.out$batch[start:nrow(metrop.out$batch),no.parameters]),myqnt)

detach(package:mcmc)