#Bayesian analysis of the Behrens Fisher Problem
#Sampling from the posterior of difference of two means
#y1 ~ N(mu1, sigma1), y2 ~ N(mu2, sigma2) are independent.
#N=#simulations, enter summary statistics, ybar, sd and n.
#to complie type "source('ctnm.R')" followed by "ctnm()"
#Copyright(2000): Sujit K. Ghosh, NC State University.


ctnm <- function(N=5000,n=c(12, 7),ybar=c(120, 101), sd=c(21.38, 20.61)){
#Default data from p.146 of Peter Lee's book.

#Posterior simulation:
#Assumes flat prior for mu and sigma

       mu1 <- ybar[1] + (sd[1]/sqrt(n[1]))*rt(N,n[1]-1) 
       mu2 <- ybar[2] + (sd[2]/sqrt(n[2]))*rt(N,n[2]-1)
       
       diff <- as.vector(mu1-mu2)
       
#Plot the posterior histogram:
       par(mfrow=c(2,2))
       hist(mu1, xlab="mu1", cex=0.8)
       hist(mu2, xlab="mu2", cex=0.8)
       hist(diff, xlab="mu1-mu2", cex=0.8)
       plot(density(diff),xlab="mu1-mu2", cex=0.8)

#Output the summary values:
      cat("Posterior summary of diff=mu1-mu2",fill=T)
      print(round(summary(diff),3))
      cat("Extreme quantiles",fill=T)
      print(round(quantile(diff,c(0.025,0.05,0.95,0.975)),3))
      cat(c("sd=",round(sqrt(var(diff)),3)),fill=T)
      cat(c("MC error:",signif(sqrt(var(diff))/sqrt(N))),fill=T)
      cat("",fill=T)

      cat("Copyright(2000): Sujit K. Ghosh",fill=T)
      cat(c("Job finished at:",date()),fill=T)
                      }