### Wichern & Johnson - Microwave Data

dclosed <- scan("http://www.stat.ufl.edu/~winner/sta4702/data/wichern/T4-1.DAT")
dopen <- scan("http://www.stat.ufl.edu/~winner/sta4702/data/wichern/T4-5.DAT")
       
microwv <- data.frame(dclosed,dopen)
attach(microwv)

microwv


#### Points a constant distance from the origin: drive distance/fairway pct
x1 <- dclosed^0.25; x2 <- dopen^0.25
n <- length(x1); p <- 2
alpha <- 0.05
crit.dist <- sqrt((p*(n-1)/(n*(n-p)))*qf(1-alpha,p,n-p))
A <- rbind(cbind(var(x1),cov(x1,x2)),cbind(cov(x1,x2),var(x2)))
mu.test <- c(0.562,0.589)

ctr <- c(mean(x1),mean(x2))
angles <- seq(0, 2*pi, length.out=200)


eigVal  <- eigen(A)$values
eigVec  <- eigen(A)$vectors
eigScl  <- eigVec %*% diag(sqrt(eigVal))  # scale eigenvectors to length = square-root
xMat    <- rbind(ctr[1] + eigScl[1, ]*crit.dist, ctr[1] - eigScl[1, ]*crit.dist)
yMat    <- rbind(ctr[2] + eigScl[2, ]*crit.dist, ctr[2] - eigScl[2, ]*crit.dist)
ellBase <- cbind(sqrt(eigVal[1])*crit.dist*cos(angles), sqrt(eigVal[2])*crit.dist*sin(angles)) # normal ellipse
ellRot  <- eigVec %*% t(ellBase)                                          # rotated ellipse

plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2,
main="100(1-a)% Confidence Ellipsoid",
xlab="x1", ylab="x2")
matlines(xMat, yMat, lty=1, lwd=2, col="blue")
points(ctr[1], ctr[2], pch=4, col="orange", lwd=3)
points(mu.test[1],mu.test[2], pch=8, col="green")
bon11 <- ctr[1] - qt(1-alpha/(2*p),n-1)*sqrt(A[1,1]/n)
bon12 <- ctr[1] + qt(1-alpha/(2*p),n-1)*sqrt(A[1,1]/n)
bon21 <- ctr[2] - qt(1-alpha/(2*p),n-1)*sqrt(A[2,2]/n)
bon22 <- ctr[2] + qt(1-alpha/(2*p),n-1)*sqrt(A[2,2]/n)
rect(bon11,bon21,bon12,bon22,border="violetred3")