Como determinar quantis (isolinhas?) De uma distribuição normal multivariada

24

insira a descrição da imagem aqui

Estou interessado em saber como calcular um quantil de uma distribuição multivariada. Nas figuras, desenhei os quantis de 5% e 95% de uma determinada distribuição normal univariada (esquerda). Para a distribuição normal multivariada correta, estou imaginando que um análogo seria uma isolina que circunda a base da função densidade. Abaixo está um exemplo da minha tentativa de calcular isso usando o pacote mvtnorm- mas sem sucesso. Suponho que isso possa ser feito calculando um contorno dos resultados da função de densidade multivariada, mas fiquei pensando se há outra alternativa ( por exemplo , analógica de qnorm). Obrigado pela ajuda.

Exemplo:

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()
Marc na caixa
fonte
3
Uma solução Mathematica é fornecida (e ilustrada para o caso 3D) em mathematica.stackexchange.com/questions/21396/… . Ele reconhece que os níveis de contorno são dados por uma distribuição qui-quadrado.
whuber
@ whuber - você se importaria de demonstrar o que quer dizer com "... o elipsóide de confiança é um contorno do inverso da matriz de covariância"? Felicidades.
Marc na caixa
2
Isso é mais fácil de ver em uma dimensão, onde a "matriz de covariância" (para uma distribuição de amostragem) é um número , então seu inverso é 1 / s 2 , pensado como um mapa quadrático em R 1 via x x 2 / s 2 . Um contorno no nível λ por definição é o conjunto de x para o qual x 2 / s 2 = λ ; ou seja, x 2 = λ s 2 ou equivalentemente x = ± s21/s2R1xx2/s2λxx2/s2=λx2=λs2. Quandoλé oquantil1-αde umadistribuiçãoχ2(1),x=±λsλ1-αχ2(1) é oquantil1-αde umadistribuiçãot(1), de onde recuperamos os limites de confiança usuais±t 1 - α ; 1 s. λ1-αt(1)±t1-α;1s
whuber
Você pode usar a primeira fórmula nesta resposta escolhendo em ( 0 , 1 ) para obter a elipse S α (a linha tracejada vermelha em suas plotagens) para qualquer xR 2α(0 0,1)SαxR2
user603

Respostas:

25

A linha de contorno é um elipsóide. O motivo é que você deve examinar o argumento da exponencial, no pdf da distribuição normal multivariada: as isolinhas seriam linhas com o mesmo argumento. Então você obtém onde Σ é a matriz de covariância. Essa é exatamente a equação de uma elipse; no caso mais simples, μ = ( 0 , 0 ) e Σ é diagonal, então você obtém ( x

(x-μ)TΣ-1(x-μ)=c
Σμ=(0 0,0 0)Σ SeΣnão for diagonal, diagonalizando você obtém o mesmo resultado.
(xσx)2+(yσy)2=c
Σ

Agora, você precisaria integrar o pdf do multivariado dentro (ou fora) da elipse e solicitar que isso seja igual ao quantil desejado. Digamos que seus quantis não sejam os habituais, mas elípticos em princípio (ou seja, você está procurando a Região de Maior Densidade, HDR, como a resposta de Tim aponta) Eu mudaria variáveis ​​no pdf para , integrei no ângulo e depois para z de 0 a z2=(x/σx)2+(y/σy)2z0 0 1-α=c Então você substitutos s = - z 2 / 2 :

1-α=0 0cdzze-z2/22π0 02πdθ=0 0cze-z2/2
s=-z2/2
0 0cze-z2/2=-c/20 0esds=(1-e-c/2)

μΣ-2emα

(x-μ)TΣ-1(x-μ)=-2emα
chuse
fonte
4

Você perguntou sobre o normal multivariado, mas iniciou sua pergunta perguntando sobre "quantil de uma distribuição multivariada" em geral. Pela redação da sua pergunta e pelo exemplo fornecido, parece que você está interessado em regiões de maior densidade . Eles são definidos por Hyndman (1996) da seguinte maneira

f(z)X100(1-α)%R(fα)X

R(fα)={x:f(x)fα}

fαPr(XR(fα))1-uma

Y=f(x)fαPr(f(x)fα)1-ααYy1,...,ymf(x)


Hyndman, RJ (1996). Computar e representar graficamente as regiões de maior densidade. The American Statistician, 50 (2), 120-126.

Tim
fonte
2

-2em(α)

0 0cze-z2/2=-c/20 0esds=(1-e-c/2)
chunjiw
fonte
1

Você pode desenhar uma elipse correspondente às distâncias 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)

Ou com círculos em torno de 95%, 75% e 50% dos dados

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))
margarida
fonte
4
Bem-vindo ao site @ user98114. Você pode fornecer algum texto para explicar o que esse código está fazendo e como ele resolve o problema do OP?
gung - Restabelece Monica