data <- read.table("~/public_html/ph5470/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))

detach(package:MASS)