library(mvtnorm)

set.seed(1234)

input.data <- read.table("http://www.biostat.umn.edu/~sudiptob/DukeSummerCourse/LinearModelExample.txt", header=T)

##Using classical least-squares or OLS

lm.out <- lm(Y ~ X1+X2+X3+X4+X5+X6, data=input.data)

summary(lm.out)

##Reading Y

Y <- data$Y

##Reading X

X <- as.matrix(data[,2:7])
ones <- rep(1,times=length(Y))
X <- cbind(ones,X)

N <- length(Y)

p <- ncol(X)

##Estimating quantities

tXX.inv <- solve(t(X)%*%X)
beta.hat <- tXX.inv%*%t(X)%*%Y
s.sq <- t(Y - X%*%beta.hat)%*%(Y - X%*%beta.hat) / (N - p)


##Simulation loops

NITER <- 500
beta <- matrix(0, nrow=NITER,ncol=p)


sigmasq <- 1/rgamma(NITER, (N-p)/2, (N-p)*s.sq/2)

for (i in 1:NITER) {
	Sigma <- sigmasq[i]*(tXX.inv)
	beta[i,] <- t(rmvnorm(1,beta.hat,Sigma))
} 

beta.qnt <- matrix(0, nrow=p, ncol=3)

for (i in 1:p) {
	beta.qnt[i,] <- quantile(beta[1:NITER,i], c(0.50, 0.025, 0.975))
}

sigmasq.qnt <- quantile(sigmasq, c(0.50, 0.025, 0.975))

## Prediction

row.extract <- c(30, 60, 90, 120, 150, 900, 930, 960, 990, 1020)

X.tilde <- X[row.extract,]

Y.tilde = matrix(0, nrow(X.tilde), NITER)


## Nested for loop segment to sample from posterior predictive distribution
## THIS IS BAD PROGRAMMING PRACTICE IN R
## Hence the segment is commented out.

##for (i in 1:nrow(X.tilde)) {
##
##  mu.tilde <- X.tilde%*%t(beta)
##  for (j in 1:NITER) {  
##    Y.tilde[i,j] <- rnorm(1, mu.tilde[i,j], sqrt(sigmasq[j]))
##  }

##}

## Sampling from the posterior predictive distribution in ONE LOOP:

for (i in 1:nrow(X.tilde)) {

  mu.tilde <- X.tilde%*%t(beta)
  Y.tilde[i,] <- rnorm(NITER, mu.tilde[i,], sqrt(sigmasq))

}



## The rows of Y.tilde contain the posterior predictive samples
## hist(Y.tilde[1,])

## Posterior predictive quantiles

Y.tilde.qnt <- matrix(0, nrow=nrow(Y.tilde), ncol=3)

for (i in 1:nrow(Y.tilde)) {
  Y.tilde.qnt[i,] <- quantile(Y.tilde[i, 1:NITER], c(0.50, 0.025, 0.975))
}


## Bayesian Model Comparisons with replicated data  

## Replication: "Prediction" from the original model

Y.rep <- matrix(0, length(Y), NITER)
Y.rep.post.mean <- rep(0, times=length(Y))
Y.rep.post.var <- rep(0, times=length(Y))

for (i in 1:nrow(Y.rep)) {

  mu.rep <- X%*%t(beta)
  Y.rep[i,] <- rnorm(NITER, mu.rep[i,], sqrt(sigmasq))
  Y.rep.post.mean[i] <- mean(Y.rep[i,])
  Y.rep.post.var[i] <- var(Y.rep[i,])

}

G <- sum((Y - Y.rep.post.mean)*(Y - Y.rep.post.mean))

P <- sum(Y.rep.post.var)

D <- G + P

##D is a posterior predictive measure of model comparisons

detach(package:mvtnorm)