¿Cómo determinar los cuantiles (isolinas?) De una distribución normal multivariada

24

ingrese la descripción de la imagen aquí

Estoy interesado en cómo se puede calcular un cuantil de una distribución multivariada. En las figuras, he dibujado los cuantiles de 5% y 95% de una distribución normal univariada dada (izquierda). Para la distribución normal multivariada correcta, estoy imaginando que un análogo sería una isolina que rodea la base de la función de densidad. A continuación se muestra un ejemplo de mi intento de calcular esto usando el paquete mvtnorm, pero sin éxito. Supongo que esto podría hacerse calculando un contorno de los resultados de la función de densidad multivariada, pero me preguntaba si hay otra alternativa ( por ejemplo , análogo de qnorm). Gracias por tu 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()

r pdf quantiles multivariate-normal multivariate-distribution Marc en la caja
fuente

3

Se proporciona una solución de Mathematica (y se ilustra para el caso 3D) en Mathica.stackexchange.com/questions/21396/… . Reconoce que los niveles de contorno están dados por una distribución de chi-cuadrado.

whuber

@whuber: ¿le importaría demostrar lo que quiere decir con "... el elipsoide de confianza es un contorno del inverso de la matriz de covarianza"? Aclamaciones.

Marc en la caja el

2

Esto es más fácil de ver en una dimensión, donde la "matriz de covarianza" (para una distribución de muestreo) es un número

, por lo que su inverso es

, considerado como un mapa cuadrático en

través de

. Un contorno en el nivel

por definición es el conjunto de

para el cual

; es decir,

o equivalente

s^{2}

$s^2$

1 / s^{2}

$1/s^2$

R^{1}

$\mathbb{R}^1$

x \to x^{2} / s^{2}

$x\to x^2/s^2$

λ

$\lambda$

x

$x$

x^{2} / s^{2} = λ

$x^2/s^2=\lambda$

x^{2} = λ s^{2}

$x^2=\lambda s^2$

. Cuando

es elcuantil

de unadistribución

,

x = \pm \sqrt{λ} s

$x=\pm\sqrt{\lambda}s$

λ

$\lambda$

1 - α

$1-\alpha$

χ^{2} (1)

$\chi^2(1)$

es elcuantil

de unadistribución

, de donde recuperamos los límites de confianza habituales

.

\sqrt{λ}

$\sqrt{\lambda}$

1 - α

$1-\alpha$

t (1)

$t(1)$

\pm t_{1 - α; 1} s

$\pm t_{1-\alpha; 1}s$

Whuber

Puede usar la primera fórmula en esta respuesta eligiendo en

para obtener la elipse

(la línea punteada roja en sus gráficos) para cualquier

α

$\alpha$

(0, 1)

$(0,1)$

S_{α}

$S_\alpha$

x \in R^{2}

$\mathbf{x}\in \mathbb{R}^2$

usuario603

25

La línea de contorno es un elipsoide. La razón es porque hay que mirar el argumento de la exponencial, en el pdf de la distribución normal multivariante: las isolinas serían líneas con el mismo argumento. Entonces obtienes donde es la matriz de covarianza. Esa es exactamente la ecuación de una elipse; en el caso más simple, y es diagonal, por lo que obtienes

(x - μ)^{T} Σ^{- 1} (x - μ) = c

$({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = c$

Σ

$\Sigma$

μ = (0, 0)

$\mu=(0,0)$

Σ

$\Sigma$

Si

no es diagonal, diagonalizando se obtiene el mismo resultado.

{(\frac{x}{σ_{x}})}^{2} + {(\frac{y}{σ_{y}})}^{2} = c

$\left(\frac{x}{\sigma_x}\right)^2+\left(\frac{y}{\sigma_y}\right)^2=c$

Σ

$\Sigma$

Ahora, tendría que integrar el pdf del multivariado dentro (o fuera) de la elipse y solicitar que sea igual al cuantil que desea. Digamos que sus cuantiles no son los habituales, sino elípticos en principio (es decir, está buscando la Región de mayor densidad, HDR, como señala la respuesta de Tim). Cambiaría las variables en el pdf a , me integraría en el ángulo y luego para de a $z^2=(x/\sigma_x)^2+(y/\sigma_y)^2$ $z$ $0$ $\sqrt{c}$ A continuación, sustitutos :

1 - α = \int_{0 0}^{\sqrt{do}} re z \frac{z {mi}^{- z^{2} / / 2}}{2 π} \int_{0 0}^{2 π} re θ = \int_{0 0}^{\sqrt{do}} z {mi}^{- z^{2} / / 2}

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

s = - z^{2} / 2

$s=-z^2/2$

\int_{0 0}^{\sqrt{do}} z {mi}^{- z^{2} / / 2} = \int_{- do / / 2}^{0 0} {mi}^{s} re s = (1 - {mi}^{- do / / 2})

$\int_0^{\sqrt{c}}z\;e^{-z^2/2}=\int_{-c/2}^{0}e^sds=(1-e^{-c/2})$

$\mu$ $\Sigma$ $-2\ln\alpha$

(X - μ)^{T} Σ^{- 1} (X - μ) = - 2 En α

$({\bf x}-\mu)^T\Sigma^{-1}({\bf x}-\mu) = -2\ln{\alpha}$

chuse
fuente

4

Usted preguntó sobre multivariante normal, pero comenzó su pregunta preguntando sobre "cuantil de una distribución multivariada" en general. Según la redacción de su pregunta y el ejemplo proporcionado, parece que está interesado en las regiones de mayor densidad . Son definidos por Hyndman (1996) de la siguiente manera

$f(z)$ $X$ $100( 1 - \alpha )\%$ $R(f_\alpha)$ $X$

$R (F_{α}) = {X : F (X) \geq F_{α}}$ $R(f_\alpha) = \{ x : f(x) \geq f_\alpha\}$
$f_\alpha$ $\Pr(X \in R(f_\alpha)) \geq 1 - a$

Los HDR se pueden obtener por integración, pero, como lo describe Hyndman, puede hacerlo utilizando un método numérico más simple. Si , entonces se puede obtener de tal manera que simplemente tomando cuantil de . Se puede estimar utilizando cuantiles de muestra de un conjunto de observaciones . El método se aplica incluso si no conocemos , pero solo tenemos un conjunto de observaciones iid. Este método funcionaría también para distribuciones multimodales. $Y = f(x)$ $f_\alpha$ $\Pr(f(x) \geq f_\alpha) \geq 1 - \alpha$ $\alpha$ $Y$ $y_1,...,y_m$ $f(x)$

Hyndman, RJ (1996). Calcular y graficar regiones de mayor densidad. El estadístico estadounidense, 50 (2), 120-126.

Tim
fuente

2

La respuesta correcta debería ser . Hubo un error en el cálculo anterior. La versión corregida: $-2*\ln(\alpha)$

\int_{0 0}^{\sqrt{do}} z {mi}^{- z^{2} / / 2} = \int_{- do / / 2}^{0 0} {mi}^{s} re s = (1 - {mi}^{- do / / 2})

$\int_0^\sqrt{c} z e^{-z^2/2} =\int_{-c/2}^0e^sds=(1-e^{-c/2})$

chunjiw
fuente

1

Podrías dibujar una elipsis 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))

margarita
fuente

44

Bienvenido al sitio @ user98114. ¿Puede proporcionar algún texto para explicar qué está haciendo este código y cómo resuelve el problema del OP?

gung - Restablece a Monica

¿Cómo determinar los cuantiles (isolinas?) De una distribución normal multivariada

Ejemplo:

Respuestas: