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

##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
sigmasq <- rep(0,times=NITER)
beta <- matrix(0, nrow=NITER,ncol=p)

library(MASS)

for (i in 1:NITER) {

	sigmasq[i] <- 1/rgamma(1, (N-p)/2, (N-p)*s.sq/2)
	Sigma <- sigmasq[i]*(tXX.inv)
	beta[i,] <- mvrnorm(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)

for (i in 1:NITER) {

  mu.tilde <- X.tilde%*%beta[i,]
  Sigma.tilde <- sigmasq[i]*diag(nrow(X.tilde))
  Y.tilde[,i] <- mvrnorm(1, mu.tilde, Sigma.tilde)
}

## 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))
}

detach(package:MASS)