library("mvtnorm")

Bayes.lm<-function(y,x,a=0.01,b=0.01,n.samples=50000){

   #This function returns n.samples samples from the model:
   #    y|alpha,beta,tau~N(x%*%beta,tau)
   #    beta~flat
   #    tau~Gamma(a,b)
   #
   # y is the vector of outcomes
   # x is the matrix predictors

   n<-length(y)
   p<-ncol(x)
   tXXinv<-solve(t(X)%*%X)
   betahat<-tXXinv%*%t(X)%*%y

   #intial values:
   tau<-rgamma(1,1,1)
   beta<-rnorm(p,0,10)

   #Initialize vectors to store the results:
   keep.tau<-rep(0,n.samples)
   keep.beta<-matrix(0,n.samples,p)

   #Start the MCMC sampler!
   for(i in 1:n.samples){
     #update tau:
     tau<-rgamma(1,n/2+a,rate=sum((y-x%*%beta)^2)/2+b)

     #update beta:
     beta<-as.vector(rmvnorm(1,betahat,tXXinv/tau))

     #store results:
     keep.tau[i]<-tau
     keep.beta[i,]<-beta

     #Display the chains every 10000 iterations:
     if(i%%10000==0){
       par(mfrow=c(2,1))
       plot(keep.tau[1:i],type="l",ylab=expression(tau),xlab="Iteration")
       plot(X%*%beta,y,ylab="Data",xlab="Predicted")
     }
   }

#return a list with the posterior samples:
list(tau=keep.tau,beta=keep.beta)}