## Nov 5, 2007
## Let's calculate the ellipse from example 5.3
rm(list = ls())

## from the text
y.bar <- t(matrix(c(0.564,0.603),2,1))
S <- matrix(c(0.0144, 0.0117, 0.0117, 0.0146),2,2)
n <- 42
p <- 2

## calculate the eigen values and vectors
eig <- eigen(S)
eig.val.1 <- eig$values[1]
eig.val.2 <- eig$values[2]

eig.vec.1 <- eig$vectors[,1]
eig.vec.2 <- eig$vectors[,2]


## calculate the F-distribution modified by the appropriate constant
mod.F <- (((n-1)*p)/((n-p)))*qf(0.95, p,n-p) 

## determine the vectors
p.1 <- sqrt(eig.val.1)*sqrt(mod.F/n)*eig.vec.1
p.2 <- sqrt(eig.val.2)*sqrt(mod.F/n)*eig.vec.2

Ell <- rbind(y.bar, y.bar+p.1, y.bar+p.2, y.bar+(-p.1), y.bar+(-p.2))

# plot
plot(Ell[,1], Ell[,2], pch=16, cex=3, col="blue")

# now use the ellipse function
# b/c in how we do this!!
library(ellipse)
plot(ellipse(S*(1/n), centre=y.bar, t=sqrt(mod.F)), pch=16, xlab="mu1", ylab="mu2")
points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue")

## add the T^2 95% confidence intervals
mu1.T.sq.interval <- c(y.bar[1] - sqrt(mod.F*(S[1,1]/n)),  y.bar[1] + sqrt(mod.F*(S[1,1]/n)))
mu2.T.sq.interval <- c(y.bar[2] - sqrt(mod.F*(S[2,2]/n)),  y.bar[2] + sqrt(mod.F*(S[2,2]/n)))

mu1.T.sq.interval
mu2.T.sq.interval

abline(v=mu1.T.sq.interval, lwd=2)
abline(h=mu2.T.sq.interval, lwd=2)

## add the Bonferroni intervals
mu1.Bon.interval <- c(y.bar[1] - qt((1-(0.05/(2*2))), n-1)*sqrt(S[1,1]/n),  y.bar[1] + qt((1-(0.05/(2*2))), n-1)*sqrt(S[1,1]/n))
mu2.Bon.interval <- c(y.bar[2] - qt((1-(0.05/(2*2))), n-1)*sqrt(S[2,2]/n),  y.bar[2] + qt((1-(0.05/(2*2))), n-1)*sqrt(S[2,2]/n))

mu1.Bon.interval
mu2.Bon.interval

abline(v=mu1.Bon.interval, lwd=2, col="green")
abline(h=mu2.Bon.interval, lwd=2, col="green")


## asymptotically, S goes toward Sigma and we can use the chi-square distribution as before, 
## when we knew mu and Sigma!!  And the data do not need to be multivariate normal!!
points(ellipse(S*(1/n), centre=y.bar, level=0.95), pch=16, col="red")

## add the large sample asymptotic intervals
mu1.LS.interval <- c(y.bar[1] - sqrt(qchisq(0.95, p))*sqrt(S[1,1]/n),  y.bar[1] + sqrt(qchisq(0.95, p))*sqrt(S[1,1]/n))
mu2.LS.interval <- c(y.bar[2] - sqrt(qchisq(0.95, p))*sqrt(S[2,2]/n),  y.bar[2] + sqrt(qchisq(0.95, p))*sqrt(S[2,2]/n))

mu1.LS.interval
mu2.LS.interval

abline(v=mu1.LS.interval, lwd=2, col="red")
abline(h=mu2.LS.interval, lwd=2, col="red")