set.seed(1234)

library(mcmc)
library(coda)

NITER <- 50000
start <- 20001
BatchLength <- 1

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

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

##Only to compare with standard linear model:

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

no.parameters = ncol(X)+1

##Defining the target function with IG priors on tau^2. Here theta = (beta, log(sigmasq))

NITER <- 150000
start <- 50001
BatchLength <- 1

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

## Jacobian:  logJ <- z
  
  out <- -(a+1+length(Y)/2.0) * z + z - (1/exp(z))*(b + 0.5*t(Y-X%*%beta)%*%(Y-X%*%beta))

  return(out)
}

lm.vcov = vcov(simple.lm)
myscale = diag(c(0.075*diag(lm.vcov), 2.5e-3))
##myscale <- diag(c(rep(0.0045,times=ncol(X)),0.0025))
myscale <- chol(myscale)
inits <- rep(10, times=no.parameters)

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


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)

p.samples <- as.mcmc(metrop.out$batch[1:nrow(metrop.out$batch),])
plot(p.samples, density=FALSE)

########Analysis with Uniform prior: Gelman (2006)#####

##Suppose sigma ~ U(a.b)

a <- 0
b <- 1000

## We transform z = log[(sigma-a)/(b-sigma)];
## z has support over entire real line;
## We can use the multivariate normal proposal
## This yields: sigma <- b - (b-a)/(1+exp(z))
## ATTENTION: We need to account for the jacobian of the transormation

NITER <- 150000
start <- 50001
BatchLength <- 1

target_uniform <- function (theta)
{
  beta <- theta[1:ncol(X)]
  z <- theta[length(theta)]
  sigma <- b - (b-a)/(1+exp(z))
  sigmasq <- sigma^2

  logJ <- z - 2*log(1+exp(z))
  
  out <- logJ - (length(Y)/2.0)*log(sigmasq) - 1/(2.0*sigmasq)*t(Y-X%*%beta)%*%(Y-X%*%beta)

  return(out)
}

lm.vcov = vcov(simple.lm)
myscale = diag(c(0.005*diag(lm.vcov), 0.015))
##myscale <- diag(c(rep(0.0045,times=ncol(X)),0.0025))
myscale <- chol(myscale)
inits <- rep(10, times=no.parameters)

metrop.out2 <- metrop(target_uniform, initial=inits, nbatch=NITER, blen=BatchLength, scale=myscale)

myqnt <- c(0.50, 0.025, 0.975)
beta2.qnt <- matrix(nrow=ncol(X),ncol=length(myqnt))

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

z <- metrop.out2$batch[start:nrow(metrop.out2$batch),no.parameters]
sigma <- b - (b-a)/(1+exp(z))
sigmasq <- sigma^2

sigmasq2.qnt <- quantile(sigmasq, myqnt)

p.samples2 <- as.mcmc(metrop.out2$batch[1:nrow(metrop.out2$batch),])
plot(p.samples2, density=FALSE)

########Another Analysis with Uniform prior: Gelman (2006)#####

## Suppose sigma ~ U(a.b)

a <- 0
b <- 1000

## We transform Z = -pi/2 + [pi/(b-a)]*(sigma-a);
## Now transform z = atan(Z), which has support over entire real line;
## This yields: sigma = a + [(b-a)/pi]*(atan(z) + pi/2) 
## ATTENTION: We need to account for the jacobian of the transormation

NITER <- 150000
start <- 50001
BatchLength <- 1

target3 <- function (theta)
{
  beta <- theta[1:ncol(X)]
  z <- theta[length(theta)]
  sigma <- a + (b-a)/pi*(atan(z) + pi/2)
  sigmasq <- sigma^2

  logJ <- -log(1+z^2)
  
  out <- logJ - (length(Y)/2.0)*log(sigmasq) - 1/(2.0*sigmasq)*t(Y-X%*%beta)%*%(Y-X%*%beta)

  return(out)
}

lm.vcov = vcov(simple.lm)
myscale = diag(c(0.085*diag(lm.vcov), 0.075))
##myscale <- diag(c(rep(0.05,times=ncol(X)),0.05))
myscale <- chol(myscale)
inits <- rep(10, times=no.parameters)

metrop.out3 <- metrop(target3, initial=inits, nbatch=NITER, blen=BatchLength, scale=myscale)

myqnt <- c(0.50, 0.025, 0.975)
beta3.qnt <- matrix(nrow=ncol(X),ncol=length(myqnt))

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

z <- metrop.out3$batch[start:nrow(metrop.out3$batch),no.parameters]
sigma <- a + (b-a)/pi*(atan(z) + pi/2)
sigmasq <- sigma^2

sigmasq3.qnt <- quantile(sigmasq, myqnt)

p.samples3 <- as.mcmc(metrop.out3$batch[1:nrow(metrop.out3$batch),])
plot(p.samples3, density=FALSE)

########Another Analysis with Uniform prior: Gelman (2006)#####

##Suppose sigma ~ U(a.b)
## We transform Z = -log[(sigma-a)/(b-a)], which ~ Exp(1);
## Now transform z = log(Z), which has support over entire real line;
## We can use the multivariate normal proposal
## ATTENTION: We need to account for the jacobian of the transormation
## This does not work that well

target_uniform2 <- function (theta)
{
  beta <- theta[1:ncol(X)]
  z <- theta[length(theta)]
  sigma <- a + (b-a)*exp(-exp(z))
  sigmasq <- sigma^2

  logJ <- z - exp(z)
  
  out <- logJ - (length(Y)/2.0)*log(sigmasq) - 1/(2.0*sigmasq)*t(Y-X%*%beta)%*%(Y-X%*%beta)

  return(out)
}

detach(package:mcmc)
detach(package:coda)