source("hdr.txt")
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)
sx <- sqrt(phi)
sy <- sqrt(psi)
meandiff <- xbar - ybar
vardiff <- phi/m + psi/n
tantheta <- (sx/sqrt(m))/(sy/sqrt(n))
theta <- atan(tantheta)
thetadeg <- 180*theta/pi
hdrint <- hbehrens(0.9,m-1,n-1,thetadeg)
cat("Get tan(theta) =",tantheta,"so theta =",theta,
    "radians or",thetadeg,"degrees\n")
cat("Find 90% HDR for BF(",m-1,",",n-1,",",thetadeg,") is",
    hdrint,"\n")
cat("so posterior HDR is",meandiff+s*hdrint,"\n")