22
$\begingroup$

enter image description here

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()
$\endgroup$
  • 3
    $\begingroup$ A Mathematica solution is given (and illustrated for the 3D case) at mathematica.stackexchange.com/questions/21396/…. It recognizes that the contour levels are given by a chi-squared distribution. $\endgroup$ – whuber Jul 18 '13 at 13:49
  • $\begingroup$ @whuber - would you mind demonstrating what you mean by "... the confidence ellipsoid is a contour of the inverse of the covariance matrix"? Cheers. $\endgroup$ – Marc in the box Aug 19 '13 at 8:53
  • 2
    $\begingroup$ This is easiest to see in one dimension, where the "covariance matrix" (for a sampling distribution) is a number $s^2$, so its inverse is $1/s^2$, thought of as a quadratic map on $\mathbb{R}^1$ via $x\to x^2/s^2$. A contour at level $\lambda$ by definition is the set of $x$ for which $x^2/s^2=\lambda$; that is, $x^2=\lambda s^2$ or equivalently $x=\pm\sqrt{\lambda}s$. When $\lambda$ is the $1-\alpha$ quantile of a $\chi^2(1)$ distribution, $\sqrt{\lambda}$ is the $1-\alpha$ quantile of a $t(1)$ distribution, whence we recover the usual confidence limits $\pm t_{1-\alpha; 1}s$. $\endgroup$ – whuber Aug 19 '13 at 14:11
  • $\begingroup$ You could use the first formula in this answer by choosing $\alpha$ in $(0,1)$ to obtain the corresponding ellipse $S_\alpha$ (the red dashed line in your plots) for any $\mathbf{x}\in \mathbb{R}^2$ $\endgroup$ – user603 Jan 18 '17 at 13:26
25
$\begingroup$

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} $$

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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}) $$

$\endgroup$
1
$\begingroup$

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))
$\endgroup$
  • 4
    $\begingroup$ Welcome to the site @user98114. Can you provide some text to explicate what this code is doing & how it resolves the OP's issue? $\endgroup$ – gung Dec 13 '15 at 23:04

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.