hnorm <- function(p,mean=0,sd=1)
  return(c(qnorm((1-p)/2,mean,sd),qnorm((1+p)/2,mean,sd)))
x <- c(22.5,20.1,23.3,22.9,23.1,22.0,22.3,23.6,24.7)
y <- c(17.0,14.9,16.0,17.4,17.4,16.5,17.2,17.2,18.0)
n <- length(x)
xbar <- mean(x)
ybar <- mean(y)
Sxx <- sum((x-xbar)^2)
Syy <- sum((y-ybar)^2)
Sxy <- sum((x-xbar)*(y-ybar))
cat("n =",n,"xbar =",xbar,"ybar =",ybar,"\n")
cat("Sxx =",Sxx,"Syy =",Syy,"Sxy =",Sxy,"\n")
r <- cor(x,y)
z <- atanh(r)
cat("Correlation",r,"with inverse hyperbolic tangent",z,"\n")
postint <- hnorm(0.95,z,sqrt(1/n))
cat("Interval for zeta",postint,"so for rho",tanh(postint),"\n")