Given a $5 \times 2$ dataset $\mathbf{X} =\left( \begin{array}{rr}-0.9&0.2\\2.4&0.7\\-1.4&1.0\\2.9&-0.5\\2.0&-1.0 \end{array} \right)$. Assume that $X\sim N_2(\mu, \Sigma)$.

Using $\hat{\mu}_{MLE}$ and $\hat{\Sigma}_{MLE}$, I need to calculate the eigenvalues and eigenvectors and draw an 80% probability ellipse. In the statistical books I've seen, it is usually stated that to draw such an ellipse, you use $\mathbf{S} = \frac{1}{n-1}\sum(\mathbf{x_j} - \mathbf{\bar{x}})(\mathbf{x_j} - \mathbf{\bar{x}})'$ and calculate the eigenvalues and eigenvectors for this matrix. Then you can draw the ellipse, using these eigenvectors and the calculated $\mathbf{\bar{x}}$.

My Question Now I now that $\hat{\Sigma}_{MLE} = \frac{n-1}{n}\mathbf{S}$. So what do I need to do here? Is it correct to use $\hat{\Sigma}_{MLE}$ to calculate the eigenvalues and -vectors and use them to find the 80% confidence interval, or do you always need to use $\mathbf{S}$?

  • $\begingroup$ If no one forces me to use $\hat{\Sigma}_{MLE}$, I will use $S$, because $S$ is unbiased estimate of $\Sigma$. $\endgroup$ – user158565 Oct 17 '18 at 15:22
  • $\begingroup$ I see that posting the question here worked better to generate some response :-) $\endgroup$ – StubbornAtom Oct 17 '18 at 15:56
  • $\begingroup$ Is what you want an ellipse that contains $X$ with 80% probability under the assumption that the true parameter values equals the MLEs? Or do you perhaps want an 80% en.wikipedia.org/wiki/Confidence_region for the unknown mean vector $\mu$? $\endgroup$ – Jarle Tufto Oct 18 '18 at 9:03
  • $\begingroup$ @JarleTufto An ellipse that contains $X$ with 80% probability. I believe that for an 80% confidence region for the unknown mean vector $\mu$, you do need to use $S$, and also the axes are divided by $sqrt(n)$? $\endgroup$ – Whizkid95 Oct 18 '18 at 9:36
  • 1
    $\begingroup$ If you want a prediction region containing $X$ with a certain probability taking into account uncertainty in estimated parameters, you can probably develop something analogous to a prediction interval for $\mu$ in the univariate normal case, based on Hotelling's T-squared distribution $\endgroup$ – Jarle Tufto Oct 26 '18 at 15:02

Assume first the the parameters $\boldsymbol\mu$ and $\boldsymbol\Sigma$ are known. Just as $\frac{x-\mu}\sigma$ is standard normal and $\frac{(x-\mu)^2}{\sigma^2}$ chi-square with 1 degree of freedom in the univariate case, the quadratic form $(\mathbf{x}-\boldsymbol\mu)^T\boldsymbol\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)$ is chi-square with $p$ degrees of freedom when $\mathbf{x}$ is multivariate normal. Hence, this pivot $$ (\mathbf{x}-\boldsymbol\mu)^T\boldsymbol\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)\le \chi_{p,\alpha}^2 \tag{1} $$ with probability $(1-\alpha)$. A probability region for $\mathbf{x}$ is found by inverting (1) with respect to $\mathbf{x}$. For points at the boundary of this set, ${\mathbf{L}^{-1}}(\mathbf{x}-\boldsymbol{\mu})$ lies on a circle with radius $\sqrt{\chi^2_{p,\alpha}}$ where $\mathbf L$ is the cholesky factor of $\boldsymbol\Sigma$ (or some other square root) such that $$ \mathbf{L}^{-1}(\mathbf{x}-\boldsymbol{\mu})=\sqrt{\chi^2_{p,\alpha}} \left[ \begin{matrix} \cos(\theta)\\ \sin(\theta) \end{matrix} \right]. $$ Hence, the boundary of the set (an ellipse) is described by the parametric curve $$ \mathbf{x}(\theta)= \boldsymbol{\mu} + \sqrt{\chi^2_{p,\alpha}}\mathbf{L} \left[ \begin{matrix} \cos(\theta)\\ \sin(\theta) \end{matrix} \right], $$ for $0<\theta <2\pi$.

If the parameters are unknown and we we use $\bar{\mathbf{x}}$ to estimate $\boldsymbol\mu$, $\mathbf{x}-\bar{\mathbf{x}} \sim N_p(0,(1+1/n))\boldsymbol{\Sigma})$. Hence, $(1+1/n)^{-1}(\mathbf{x}-\bar{\mathbf{x}})^T\boldsymbol\Sigma^{-1}(\mathbf{x}-\bar{\mathbf{x}})$ is chi-square with $p$ degrees of freedom. Substituting $\boldsymbol\Sigma$ by its estimate $\hat{\boldsymbol\Sigma}=\frac1{n-1}\mathbf{X}^T \mathbf{X}$ the resulting pivot is instead Hotelling $T$-squared distributed with $p$ and $n-p$ degrees of freedom (analogous to the $F_{1,n-1}$ distributed squared $t$-statistic in the univariate case) such that $$ \Big(1+\frac1n\Big)^{-1}(\mathbf{x}-\bar{\mathbf{x}})^T\hat{\boldsymbol\Sigma}^{-1}(\mathbf{x}-\bar{\mathbf{x}}) \le T^2_{p,n-p,\alpha} \tag{2} $$ with probability $(1-\alpha)$. Because the Hotelling $T$-squared is just a rescaled $F$-distribution, the above quantile equals $\frac{p(n-1)}{n-p}F_{p,n-p,\alpha}$.

Inverting (2) with respect to $\mathbf{x}$ leads to a prediction region with boundary described by the parametric curve $$ \mathbf{x}(\theta)= \bar{\mathbf x} + \sqrt{\Big(1+\frac1n\Big)\frac{p(n-1)}{n-p}F_{p,n-p,\alpha}}\hat{\mathbf{L}} \left[ \begin{matrix} \cos(\theta)\\ \sin(\theta) \end{matrix} \right] $$ where $\hat{\mathbf L}$ is the cholesky factor of the sample variance matrix $\hat{\boldsymbol\Sigma}$.

Code computing this for the data in the original question:

pred.int.mvnorm <- function(x, alpha=.05) {
  p <- ncol(x)
  n <- nrow(x)
  Sigmahat <- var(x)
  xbar <- apply(x,2,mean)
  theta <- seq(0, 2*pi, length=100)
  polygon <- xbar + 
    sqrt(p*(n-1)/(n-p)*(1 + 1/n)*qf(alpha, p, n - p, lower.tail = FALSE))*
    t(chol(Sigmahat)) %*% 
    rbind(cos(theta), sin(theta))
x <- matrix(c(-0.9,2.4,-1.4,2.9,2.0,0.2,0.7,1.0,-0.5,-1.0),ncol=2)
plot(pred.int.mvnorm(x), type="l",xlab=expression(x[1]),ylab=expression(x[2]))

enter image description here

More code testing the coverage

hits <- 0
for (i in 1:1e+5) {
  x <- rmvnorm(6, sigma = diag(2))
  pred.int <- pred.int.mvnorm(x[-1,])
  x <- x[1,]
  if (point.in.polygon(x[1], x[2], pred.int[,1], pred.int[,2])==1)
    hits <- hits + 1
[1] 94955
| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.