Jak określić kwantyle (izoliny?) O wielowymiarowym rozkładzie normalnym

24

wprowadź opis zdjęcia tutaj

Interesuje mnie, jak można obliczyć kwantyl rozkładu wielowymiarowego. Na rysunkach narysowałem 5% i 95% kwantyli danego rozkładu jednowymiarowego normalnego (po lewej). Dla właściwego wielowymiarowego rozkładu normalnego wyobrażam sobie, że analog byłby izoliną otaczającą podstawę funkcji gęstości. Poniżej znajduje się przykład mojej próby obliczenia tego za pomocą pakietu mvtnorm- ale bez powodzenia. Przypuszczam, że można tego dokonać, obliczając kontur wyników funkcji wielowymiarowej gęstości, ale zastanawiałem się, czy istnieje inna alternatywa ( np. Analog qnorm). Dzięki za pomoc.

Przykład:

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 w pudełku
źródło
3
Mathematica roztwór podano (zilustrowanych dla przypadku, 3d) w mathematica.stackexchange.com/questions/21396/... . Rozpoznaje, że poziomy konturu są podawane przez rozkład chi-kwadrat.
whuber
@ whuber - czy mógłbyś zademonstrować, co masz na myśli mówiąc „... elipsoida zaufania jest konturem odwrotności macierzy kowariancji”? Twoje zdrowie.
Marc w pudełku,
2
Najłatwiej jest to zobaczyć w jednym wymiarze, gdzie „macierz kowariancji” (dla rozkładu próbkowania) jest liczbą , więc jej odwrotność wynosi 1 / s 2 , uważana za mapę kwadratową na R 1 poprzez x x 2 / s 2 . Kontur na poziomie λ z definicji jest zbiorem x, dla którego x 2 / s 2 = λ ; to znaczy x 2 = λ s 2 lub równoważnie x = ± s21/s2R1xx2/s2λxx2/s2=λx2=λs2. Gdyλjestkwantylem1-αrozkładuχ2(1),x=±λsλ1αχ2(1) jestkwantylem1-αrozkładut(1), skąd odzyskujemy zwykłe granice ufności±t 1 - α ; 1 s. λ1αt(1)±t1α;1s
whuber
Możesz użyć pierwszej formuły w tej odpowiedzi, wybierając w ( 0 , 1 ), aby uzyskać odpowiednią elipsę S α (czerwoną przerywaną linią na wykresach) dla dowolnego xR 2α(0,1)SαxR2
użytkownik603

Odpowiedzi:

25

Linia konturu jest elipsoidą. Powodem jest to, że musisz spojrzeć na argument wykładniczy w pliku pdf wielowymiarowego rozkładu normalnego: izolinie byłyby liniami z tym samym argumentem. Otrzymujesz gdzie Σ jest macierzą kowariancji. To jest dokładnie równanie elipsy; w najprostszym przypadku μ = ( 0 , 0 ) i Σ jest przekątna, więc otrzymujesz ( x

(xμ)TΣ1(xμ)=c
Σμ=(0,0)Σ JeśliΣnie jest przekątna, przekątna daje ten sam wynik.
(xσx)2+(yσy)2=c
Σ

Teraz musiałbyś zintegrować pdf wielowymiarowy w elipsie (lub poza nią) i zażądać, aby była ona równa żądanemu kwantylowi. Powiedzmy, że twoje kwantyle nie są zwykłymi, ale w zasadzie eliptycznymi (tzn. Szukasz regionu o największej gęstości, HDR, jak wskazuje Tim). Zmieniłbym zmienne w pdf na , całkuj w kącie, a następnie dla z od 0 do z2=(x/σx)2+(y/σy)2z0 1-α=c Następnie zastąpić s = - Z 2 / 2 :

1α=0cdzzez2/22π02πdθ=0czez2/2
s=z2/2
0czez2/2=c/20esds=(1ec/2)

μΣ2lnα

(xμ)TΣ1(xμ)=2lnα
chuse
źródło
4

Pytałeś o normalną wielowymiarową normę, ale zacząłeś od pytania o „kwantyl rozkładu wielowymiarowego” w ogóle. Z treści pytania i podanego przykładu wynika, że ​​interesują Cię regiony o największej gęstości . Są one zdefiniowane przez Hyndman (1996) w następujący sposób

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

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

fαPr(XR(fα))1a

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


Hyndman, RJ (1996). Obliczanie i tworzenie wykresów regionów o największej gęstości. The American Statistician, 50 (2), 120-126.

Tim
źródło
2

2ln(α)

0czez2/2=c/20esds=(1ec/2)
Czunjiw
źródło
1

Możesz narysować elipsy odpowiadające odległościom Mahalanobisa.

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)

Lub z kręgami około 95%, 75% i 50% danych

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))
stokrotka
źródło
4
Witamy na stronie @ user98114. Czy możesz podać tekst wyjaśniający, co robi ten kod i jak rozwiązuje problem PO?
gung - Przywróć Monikę