source("hdr.txt")
x <- c(2.31035,2.31026,2.31024,2.31012,2.31027,2.31017,2.30986,2.31010,
       2.31001,2.31024,2.31010,2.31028)
y <- c(2.30143,2.29890,2.29816,2.30182,2.29869,2.29940,2.29849,2.29889)
m <- length(x)
n <- length(y)
nux <- m-1
nuy <- n-1
xbar <- mean(x)
ybar <- mean(y)
phi <- var(x)
psi <- var(y)
cat("m =",m," n =",n,"xbar =",xbar,"ybar =",ybar,"\n")
cat("phi =",phi,"psi =",psi,"\n")
k <- phi/psi
hdrint <- hf(0.9,nux,nuy)
postmed <- k*(nux/(nux+1))*(nuy-2)/nuy
cat("Find interval for F corresponding to 90% HDR for log F(",
    m-1,",",n-1,") is\n")
cat(hdrint,"so posterior HDR is",sort(k/hdrint),"\n")
cat("Posterior mode is",postmed,"\n")