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)
df <- m+n-2
cat("Mean diff",meandiff,"and variance",vardiff,"\n")
cat("Posterior HDR",meandiff+s*ht(0.9,df),"\n")