#
# Exact Bayes inference for normal data with a conjugate prior.
#
bayesnorm <- function(n,xbar,s,m,c,a,b){
  cpost <- c+n
  apost <- ifelse(c>0,a+n,a+n-1)
  mpost <- (c*m+n*xbar)/cpost
  bpost <- b+(n-1)*(s^2)+c*n*((xbar-m)^2)/(c+n)
  scalepar <- sqrt(bpost/(cpost*apost))
  cpostx <- 1/(1+1/cpost)
  scaleparx <- sqrt(bpost/(cpostx*apost))
  return(list('mupars'=c(mpost,scalepar,apost),'xpars'=c(mpost,scaleparx,apost),'taupars'=c(apost*0.5,bpost*0.5)))
}  
#
# Inference via Gibbs sampling with a semi-conjugate prior.
#
bayesncnorm <- function(n,xbar,s,m,c,a,b,iters=1000){
  mu <- rep(0,iters)
  tau <- rep(0,iters)
  mu[1]<-xbar
  tau[1] <- rgamma(1,0.5*(a+n),0.5*(b+(n-1)*(s^2)))
  for (i in 2:iters){
    mu[i] <- rnorm(1,(c*m+n*tau[i-1]*xbar)/(c+n*tau[i-1]),sqrt(1/(c+n*tau[i-1])))
    tau[i] <- rgamma(1,0.5*(a+n),0.5*(b+(n-1)*(s^2)+n*((xbar-mu[i])^2)))
  }
  sigma <- 1/sqrt(tau)
  xpred <- rnorm(iters,mu,sigma)
  return(list("mu"=mu,"sigma"=sigma,"tau"=tau,"xpred"=xpred))
}
#
# The Bayesian solution to the Behrens Fisher problem via simulation.
#
bayesbfsim <- function(n,xbar,s,simnum=1000){
  simt1 <- s[1]*rt(simnum,n[1]-1)/sqrt(n[1])
  simt2 <- s[2]*rt(simnum,n[2]-1)/sqrt(n[2])
  delta <- xbar[1]-xbar[2]+simt1-simt2
  return(delta)
}



#
# Example
# -------
# Copy the data in the text file and use the following commands to read it in.
#
normtemp <-read.table(stdin())
names(normtemp)
degreesF <- normtemp$V1
#
# Calculate a classical 95% c.i.
#
t.test(degreesF)
#
# Calculate summary statistics for a Bayesian analysis.
#
n <- length(degreesF)
xbar <- mean(degreesF)
s <- sd(degreesF)
m <- 98.6
c <- 1
a <- 1
b <- 1
v <- bayesnorm(n,xbar,s,m,c,a,b)
mupars <- v$mupars
bayesci <- mupars[1]+mupars[2]*qt(c(0.025,0.975),mupars[3])
#
# Also compare predictive cis
#
xclass <- qnorm(c(0.025,0.975),xbar,s*(1+1/sqrt(n)))
xpars <- v$xpars
xci <- xpars[1]+xpars[2]*qt(c(0.025,0.975),xpars[3])
#
# Now do a semi-conjugate analysis.
#
w <- bayesncnorm(n,xbar,s,m,c,a,b,10000)
mu <- w$mu
quantile(mu,c(0.025,0.975))
quantile(w$xpred,c(0.025,0.975))
#
# Two sample problem.
#
sex <- normtemp$V2
men <- degreesF[sex==1]
women <- degreesF[sex==2]
t.test(men,women)
n <- c(length(men),length(women))
xbar <- c(mean(men),mean(women))
s <- c(sd(men),sd(women))
simnum <- 10000
delta <- bayesbfsim(n,xbar,s,simnum)
hist(delta,xlab='delta',ylab='f',probability=TRUE,main='')
quantile(delta,c(0.025,0.975))