Cómo determinar los cuantiles (¿aislamientos?) de una distribución normal multivariante

Question

Me interesa saber cómo se puede calcular un cuantil de una distribución multivariante. En las figuras, he dibujado los cuantiles del 5% y del 95% de una determinada distribución normal univariante (izquierda). Para la distribución normal multivariante de la derecha, estoy imaginando que un análogo sería una isolínea que rodea la base de la función de densidad. A continuación se muestra un ejemplo de mi intento de calcular esto utilizando el paquete mvtnorm - pero sin éxito. Supongo que esto podría hacerse calculando un contorno de los resultados de la función de densidad multivariante, pero me preguntaba si hay otra alternativa ( Por ejemplo análogo de qnorm ). Gracias por su ayuda.

Ejemplo:

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

La línea de contorno es un elipsoide. La razón es que hay que mirar el argumento de la exponencial, en la pdf de la distribución normal multivariante: las isolíneas serían líneas con el mismo argumento. Entonces se obtiene $$ ({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = c $$ donde $\Sigma$ es la matriz de covarianza. Esa es exactamente la ecuación de una elipse; en el caso más simple, $\mu=(0,0)$ y $\Sigma$ es diagonal, por lo que se obtiene $$ \left(\frac{x}{\sigma_x}\right)^2+\left(\frac{y}{\sigma_y}\right)^2=c $$ Si $\Sigma$ no es diagonal, diagonalizando se obtiene el mismo resultado.

Ahora, tendrías que integrar la pdf de la multivariable dentro (o fuera) de la elipse y pedir que ésta sea igual al cuantil que quieres. Digamos que tus cuantiles no son los habituales, sino que en principio son elípticos (es decir, buscas la región de mayor densidad, HDR, como señala Tim answer). Yo cambiaría las variables en el pdf por $z^2=(x/\sigma_x)^2+(y/\sigma_y)^2$ integrar en el ángulo y luego para $z$ de $0$ a $\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} $$ Entonces sustituye $s=-z^2/2$ : $$ \int_0^{\sqrt{c}}z\;e^{-z^2/2}=\int_{-c/2}^{0}e^sds=(1-e^{-c/2})$$

Así que, en principio, hay que buscar la elipse centrada en $\mu$ con eje sobre los vectores propios de $\Sigma$ y el radio efectivo $-2\ln\alpha$ : $$ ({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = -2\ln{\alpha} $$

Answer 2

Has preguntado por la normal multivariante, pero has empezado tu pregunta preguntando por el "cuantil de una distribución multivariante" en general. Por la redacción de su pregunta y el ejemplo proporcionado, parece que está interesado en regiones de mayor densidad . Hyndman (1996) los define de la siguiente manera

Dejemos que $f(z)$ sea la función de densidad de una variable aleatoria $X$ . Entonces el $100( 1 - \alpha )\%$ El HDR es el subconjunto $R(f_\alpha)$ del espacio muestral de $X$ tal que

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

donde $f_\alpha$ es la mayor constante tal que $\Pr(X \in R(f_\alpha)) \geq 1 - a$ .

Los HDR se pueden obtener por integración pero, como describe Hyndman, se puede hacer con un método numérico más sencillo. Si $Y = f(x)$ , entonces se puede obtener $f_\alpha$ tal que $\Pr(f(x) \geq f_\alpha) \geq 1 - \alpha$ simplemente tomando $\alpha$ cuantil de $Y$ . Se puede estimar utilizando Cuantiles de la muestra a partir de un conjunto de observaciones $y_1,...,y_m$ . El método se aplica incluso si no sabemos $f(x)$ pero sólo tienen un conjunto de observaciones i.i.d. Este método también funcionaría para distribuciones multimodales.

---

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

Answer 3

La respuesta correcta debería ser $-2*\ln(\alpha)$ . Había un error en el cálculo anterior. La versión corregida: $$ \int_0^\sqrt{c} z e^{-z^2/2} =\int_{-c/2}^0e^sds=(1-e^{-c/2}) $$

Answer 4

Se puede dibujar una elipse correspondiente a las distancias de Mahalanobis.

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)

O con círculos alrededor del 95%, 75% y 50% de los datos

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