Bayes.t.test<-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(alpha+x*beta,tau)
   #    alpha,beta~flat
   #    tau~Gamma(a,b)
   #
   # y is the vector of outcomes
   # x is the vector of binary predictors

   n<-length(y)
   nX1<-sum(x)

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

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

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

     #update alpha:
     alpha<-rnorm(1,mean(y-x*beta),1/sqrt(n*tau))

     #update beta:
     beta<-rnorm(1,mean(y[x==1]-alpha),1/sqrt(nX1*tau))

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

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

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

#Generate data:
n<-100                     #samples size of the simulated data
true.alpha<-0              #true intercept
true.beta<-1               #true slope
true.tau<-2                #true precision
n.samples<-50000           #number of MCMC samples to generate

set.seed(2008)
x<-rbinom(n,1,0.5)
y<-rnorm(n,true.alpha+x*true.beta,1/sqrt(true.tau))

#fit the model
fit<-Bayes.t.test(y,x,n.samples=50000)


#summarize the posterior distributions:
par(mfrow=c(2,2))
hist(1/sqrt(fit$tau[10000:50000]),xlab="1/sqrt(tau)",main="Posterior of the error sd")
hist(fit$alpha[10000:50000],xlab="alpha",main="Posterior of the intercept")
hist(fit$beta[10000:50000],xlab="beta",main="Posterior of the slope")
plot(fit$alpha[10000:50000],fit$beta[10000:50000],xlab="alpha",ylab="beta",
   main="Joint Post of the int and slope")

print("Posterior mean of beta, sd of beta, cor(alpha,beta)")
print(mean(fit$beta[10000:50000]))
print(sd(fit$beta[10000:50000]))
print(cor(fit$alpha[10000:50000],fit$beta[10000:50000]))

#Compare with least squares:
print(summary(lm(y~x)))