hnorm <- function(p,mean=0,sd=1)
  return(c(qnorm((1-p)/2,mean,sd),qnorm((1+p)/2,mean,sd)))
x <- c(134,146,104,119,124,161,107,83,113,129,97,123)
y <- c(70,118,101,85,107,132,94)
m <- length(x)
n <- length(y)
xbar <- mean(x)
ybar <- mean(y)
phi <- var(x)
psi <- var(y)
cat("m =",m," n =",n,"xbar =",xbar,"ybar =",ybar,
    "phi =",phi,"psi =",psi,"\n")
meandiff <- xbar - ybar
vardiff <- phi/m + psi/n
s <- sqrt(vardiff)
cat("Mean diff",meandiff,"and variance",vardiff,"\n")
cat("Posterior 90% HDR",hnorm(0.9,meandiff,s),"\n")
cat("Posterior probability diff > 0 is",
    1-pnorm(0,meandiff,sqrt(vardiff)),"\n")
z <- (meandiff-0)/sqrt(vardiff)
omega <- phi + psi
B <- sqrt(1+omega/(phi/m+psi/n))*exp(-0.5*z^2/(1+(phi/m+psi/n)/omega))
cat("Bayes factor is",B,"\n")