How to determine quantiles (isolines?) of a multivariate normal distribution

Question

I'm interested in how one can calculate a quantile of a multivariate distribution. In the figures, I have drawn the 5% and 95% quantiles of a given univariate normal distribution (left). For the right multivariate normal distribution, I am imagining that an analog would be an isoline that encircles the base of the density function. Below is an example of my attempt to calculate this using the package mvtnorm - but to no success. I suppose this could be done by calculating a contour of the results of the multivariate density function, but I was wondering if there is another alternative (e.g., analog of qnorm). Thanks for your help.

Example:

mu <- 5
sigma <- 2 
vals <- seq(-2,12,,100)
ds <- dnorm(vals, mean=mu, sd=sigma)

plot(vals, ds, t="l")
qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)
abline(v=qs, col=2, lty=2)


#install.packages("mvtnorm")
require(mvtnorm)
n <- 2
mmu <- rep(mu, n)
msigma <- rep(sigma, n)
mcov <- diag(msigma^2)
mvals <- expand.grid(seq(-2,12,,100), seq(-2,12,,100))
mvds <- dmvnorm(x=mvals, mean=mmu, sigma=mcov)

persp(matrix(mvds,100,100), axes=FALSE)
mvqs <- qmvnorm(0.95, mean=mmu, sigma=mcov, tail = "both") #?

#ex. plot   
png("tmp.png", width=8, height=4, units="in", res=400)
par(mfcol=c(1,2))

#univariate
plot(vals, ds, t="l")
qs <- qnorm(c(0.05, 0.95), mean=mu, sd=sigma)
abline(v=qs, col=2, lty=2)

#multivariate
pmat <- persp(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), axes=FALSE, shade=TRUE, lty=0)
cont <- contourLines(seq(-2,12,,100), seq(-2,12,,100), matrix(mvds,100,100), levels=0.05^2)
lines(trans3d(cont[[1]]$x, cont[[1]]$y, cont[[1]]$level, pmat), col=2, lty=2)

dev.off()

Answer 1

The contour line is an ellipsoid. The reason is because you have to look at the argument of the exponential, in the pdf of the multivariate normal distribution: the isolines would be lines with the same argument. Then you get $$ ({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = c $$ where $\Sigma$ is the covariance matrix. That is exactly the equation of an ellipse; in the simplest case, $\mu=(0,0)$ and $\Sigma$ is diagonal, so you get $$ \left(\frac{x}{\sigma_x}\right)^2+\left(\frac{y}{\sigma_y}\right)^2=c $$ If $\Sigma$ is not diagonal, diagonalizing you get the same result.

Now, you would have to integrate the pdf of the multivariate inside (or outside) the ellipse and request that this is equal to the quantile you want. Let's say that your quantiles are not the usual ones, but elliptical in principle (i.e. you are looking for the Highest Density Region, HDR, as Tim answer points out). I would change variables in the pdf to $z^2=(x/\sigma_x)^2+(y/\sigma_y)^2$, integrate in the angle and then for $z$ from $0$ to $\sqrt{c}$ $$ 1-\alpha=\int_0^{\sqrt{c}}dz\frac{z\;e^{-z^2/2}}{2\pi}\int_0^{2\pi}d\theta=\int_0^{\sqrt{c}}z\;e^{-z^2/2} $$ Then you substitute $s=-z^2/2$: $$ \int_0^{\sqrt{c}}z\;e^{-z^2/2}=\int_{-c/2}^{0}e^sds=(1-e^{-c/2})$$

So in principle, you have to look for the ellipse centered in $\mu$, with axis over the eigenvectors of $\Sigma$ and effective radius $-2\ln\alpha$: $$ ({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = -2\ln{\alpha} $$

Answer 2

You asked about multivariate normal, but started your question with asking about "quantile of a multivariate distribution" in general. From wording of your question and the example provided it seems that you are interested in highest density regions. They are defined by Hyndman (1996) as following

Let $f(z)$ be the density function of a random variable $X$ . Then the $100( 1 - \alpha )\%$ HDR is the subset $R(f_\alpha)$ of the sample space of $X$ such that

$$ R(f_\alpha) = \{ x : f(x) \geq f_\alpha\}$$

where $f_\alpha$ is the largest constant such that $\Pr(X \in R(f_\alpha)) \geq 1 - a$.

HDR's can be obtained by integration but, as described by Hyndman, you can do it using a simpler, numerical method. If $Y = f(x)$, then you can obtain $f_\alpha$ such that $\Pr(f(x) \geq f_\alpha) \geq 1 - \alpha$ simply by taking $\alpha$ quantile of $Y$. It can be estimated using sample quantiles from a set of observations $y_1,...,y_m$. The method applies even if we do not know $f(x)$, but have only a set of i.i.d. observations. This method would work also for multimodal distributions.

---

Hyndman, R.J. (1996). Computing and graphing highest density regions. The American Statistician, 50(2), 120-126.

Answer 3

The correct answer should be $-2*\ln(\alpha)$. There was a mistake in the calculation above. The corrected version: $$ \int_0^\sqrt{c} z e^{-z^2/2} =\int_{-c/2}^0e^sds=(1-e^{-c/2}) $$

Answer 4

You could draw an ellipses corresponding to Mahalanobis distances.

library(chemometrics)
data(glass)
data(glass.grp)
x=glass[,c(2,7)]
require(robustbase)
x.mcd=covMcd(x)
drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=0.90)

Or with circles around 95%, 75%, and 50% of data

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))