Purpose
The objective is to give out bands for the mean vector in multivariate hypothesis testing framework

Simulate a trivariate independent normal distribution
Use Hotelling T square to test out the hypotheses

> library(mnormt)
> runs <- 10000
> data <- cbind(rnorm(runs, 3, 1), rnorm(runs, 3, 1), rnorm(runs,
+     3, 1))
> k <- ncol(data)
> n <- nrow(data)
> xbar <- apply(data, 2, mean)
> mubar <- rep(3, k)
> dbar <- xbar - mubar
> v <- var(data)
> t2 <- n * dbar %*% solve(v) %*% dbar
> F <- (n - k) * t2/((n - 1) * k)
> P <- 1 - pf(F, k, n - k)
> print(P)
          [,1]
[1,] 0.6345915

Use mnormt and check the hypothesis using Hotelling T square

> runs <- 10000
> data <- rmnorm(runs, mean = c(1, 2, 3), varcov = diag(3))
> k <- ncol(data)
> n <- nrow(data)
> xbar <- apply(data, 2, mean)
> mubar <- c(1, 2, 3)
> dbar <- xbar - mubar
> v <- var(data)
> t2 <- n * dbar %*% solve(v) %*% dbar
> F <- (n - k) * t2/((n - 1) * k)
> P <- 1 - pf(F, k, n - k)
> print(P)
          [,1]
[1,] 0.5551741

So the function for Hoteling T square is

> HotellingTsquare <- function(xbar, mubar, covariance.mat, k,
+     n) {
+     dbar <- xbar - mubar
+     t2 <- n * dbar %*% solve(covariance.mat) %*% dbar
+     F <- (n - k) * t2/((n - 1) * k)
+     P <- 1 - pf(F, k, n - k)
+     return(P)
+ }
> HotellingTsquare(xbar, mubar, v, k, n)
          [,1]
[1,] 0.5551741