## MVN Inference examples 
data(USArrests)
X = USArrests[,c(1,2)]
plot(X)

# Assume the dataset is from a mvn
# Estimates  
mu = colMeans(X)
S = cov(X)

# Confidence Intervals 
install.packages("car")
library(car)

n = dim(X)[1]
p = dim(X)[2]
radius = sqrt(p / (n-p) * qf(0.95, p, n - p) )
ellipse(mu, S, radius, log="", center.pch=19, center.cex=1.5,  draw=TRUE)

radius = sqrt(p / (n-p) * qf(0.99, p, n - p) )
ellipse(mu, S, radius, log="", center.pch=19, center.cex=1.5, draw=TRUE)

# Simulataneous CIs
a1 = c(1,0)
a2 = c(0,1)
a3 = c(1,0.01)
h = (n-1) * p / (n-p) * qf(0.95, p, n-p)
c(t(a1) %*% mu - sqrt(h * t(a1) %*% S %*% a1 / n ) , 
  t(a1) %*% mu + sqrt(h * t(a1) %*% S %*% a1 / n ) )

c(t(a2) %*% mu - sqrt(h * t(a2) %*% S %*% a2 / n ) , 
  t(a2) %*% mu + sqrt(h * t(a2) %*% S %*% a2 / n ) )

c(t(a3) %*% mu - sqrt(h * t(a3) %*% S %*% a3 / n ) ,
  t(a3) %*% mu + sqrt(h * t(a3) %*% S %*% a3 / n ) )

## Testing 
# Case I:  mu0 = (7, 140), Sigma = S known
mu0 = c(7,140)
W = n * t(mu - mu0) %*% solve(S) %*% (mu - mu0)
pvalue <- 1-pchisq(W,p) 

# Case II: mu0 = (7, 140), Sigmaunknown
W2 = (n - p) / p * n / (n-1) * t(mu - mu0) %*% solve(S) %*% (mu - mu0)
pvalue2 <- 1 - pf(W,p,n-p)

# Case III: 
Sigma0 = matrix(c(20,0,0,7000),2,2)
W = - n * log(det(solve(Sigma0)%*%S)) - n * p + n * sum(eigen( solve(Sigma0)%*%S )$values)
W
pvalue <- 1-pchisq(W,3)