source("hdr.txt")
m <- 7
n <- 7
xbar <- 448.6
ybar <- 408.8
sx <- sqrt(308.6)
sy <- sqrt(1251.3)
meandiff <- xbar - ybar
vardiff <- sx^2/m + sy^2/n
s <- sqrt(vardiff)
tantheta <- (sx/sqrt(m))/(sy/sqrt(n))
theta <- atan(tantheta)
thetadeg <- 180*theta/pi
hdrint <- hbehrens(0.5,m-1,n-1,thetadeg)
cat("Get tan(theta) =",tantheta,"so theta =",theta,
    "radians or",thetadeg,"degrees\n")
cat("Find 50% HDR for BF(",m-1,",",n-1,",",thetadeg,") is",
    hdrint,"\n")
cat("so posterior HDR is",meandiff+s*hdrint,"\n")