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

, uważana za mapę kwadratową na

poprzez

. Kontur na poziomie

z definicji jest zbiorem

dla którego

; to znaczy

lub równoważnie

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$

. Gdy

jestkwantylem

rozkładu

,

x = \pm \sqrt{λ} s

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

λ

$\lambda$

1 - α

$1-\alpha$

χ^{2} (1)

$\chi^2(1)$

jestkwantylem

rozkładu

, skąd odzyskujemy zwykłe granice ufności

.

\sqrt{λ}

$\sqrt{\lambda}$

1 - α

$1-\alpha$

t (1)

$t(1)$

\pm t_{1 - α; 1} s

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

— whuber

Możesz użyć pierwszej formuły w tej odpowiedzi, wybierając

w

aby uzyskać odpowiednią elipsę

(czerwoną przerywaną linią na wykresach) dla dowolnego

α

$\alpha$

(0, 1)

$(0,1)$

S_{α}

$S_\alpha$

x \in R^{2}

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

— użytkownik603

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 i jest przekątna, więc otrzymujesz

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

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

Σ

$\Sigma$

μ = (0, 0)

$\mu=(0,0)$

Σ

$\Sigma$

Jeśli

nie jest przekątna, przekątna daje ten sam wynik.

{(\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$

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 od do $z^2=(x/\sigma_x)^2+(y/\sigma_y)^2$ $z$ $0$ $\sqrt{c}$ Następnie zastąpić :

1 - α = \int_{0}^{\sqrt{c}} d z \frac{z e^{- z^{2} / 2}}{2 π} \int_{0}^{2 π} d θ = \int_{0}^{\sqrt{c}} z e^{- 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}^{\sqrt{c}} z e^{- z^{2} / 2} = \int_{- c / 2}^{0} e^{s} d s = (1 - e^{- c / 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 \ln α

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

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

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

$-2*\ln(\alpha)$

\int_{0}^{\sqrt{c}} z e^{- z^{2} / 2} = \int_{- c / 2}^{0} e^{s} d s = (1 - e^{- c / 2})

$\int_0^\sqrt{c} z e^{-z^2/2} =\int_{-c/2}^0e^sds=(1-e^{-c/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ę