Wie man Quantile (Isolinien?) Einer multivariaten Normalverteilung bestimmt

24

Bildbeschreibung hier eingeben

Mich interessiert, wie man ein Quantil einer multivariaten Verteilung berechnen kann. In den Abbildungen habe ich die 5% - und 95% -Quantile einer gegebenen univariaten Normalverteilung gezeichnet (links). Für die richtige multivariate Normalverteilung stelle ich mir vor, dass ein Analog eine Isolinie ist, die die Basis der Dichtefunktion umgibt. Unten ist ein Beispiel für meinen Versuch, dies mit dem Paket zu berechnen mvtnorm- aber ohne Erfolg. Ich nehme an, dies könnte durch Berechnen einer Kontur der Ergebnisse der multivariaten Dichtefunktion geschehen, aber ich habe mich gefragt, ob es eine andere Alternative gibt ( z. B. Analog von qnorm). Danke für Ihre Hilfe.

Beispiel:

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 in der Kiste
quelle
3
Eine Mathematica- Lösung (und eine Abbildung für den 3D-Fall) finden Sie unter mathematica.stackexchange.com/questions/21396/… . Es erkennt, dass die Konturebenen durch eine Chi-Quadrat-Verteilung gegeben sind.
Whuber
@whuber - würde es Ihnen etwas ausmachen, zu demonstrieren, was Sie mit "... das Vertrauensellipsoid ist eine Kontur der Inversen der Kovarianzmatrix" meinen? Prost.
Marc in der Box
2
Dies ist am einfachsten in einer Dimension zu sehen, in der die "Kovarianzmatrix" (für eine Stichprobenverteilung) eine Zahl , ihre Inverse also 1 / s 2 ist , was als quadratische Abbildung auf R 1 über x x 2 angesehen wird / s 2 . Eine Kontur auf der Ebene λ ist per Definition die Menge von x, für die x 2 / s 2 = λ ist ; das heißt x 2 = λ s 2 oder äquivalent x = ± s21/s2R1xx2/s2λxx2/s2=λx2=λs2. Wennλdas1-α-Quantil einerχ2(1)-Verteilung ist, istx=±λsλ1-αχ2(1) ist das1-α-Quantil einert(1)-Verteilung, aus der wir die üblichen Konfidenzgrenzen±t 1 - α gewinnen ; 1 s. λ1-αt(1)±t1-α;1s
Whuber
Sie können die erste Formel in dieser Antwort verwenden, indem Sie in ( 0 , 1 ) wählen , um die entsprechende Ellipse S α (die rote gestrichelte Linie in Ihren Plots) für xR 2α(0,1)SαxR2
user603 am

Antworten:

25

Die Konturlinie ist ein Ellipsoid. Der Grund dafür ist, dass Sie sich das Argument des Exponentials im PDF der multivariaten Normalverteilung ansehen müssen: Die Isolinien wären Linien mit demselben Argument. Dann bekommst du wobei Σ die Kovarianzmatrix ist. Das ist genau die Gleichung einer Ellipse; im einfachsten Fall ist μ = ( 0 , 0 ) und Σ ist diagonal, so dass Sie ( x erhalten

(x-μ)TΣ-1(x-μ)=c
Σμ=(0,0)Σ WennΣnicht diagonal ist, erhalten Sie beim Diagonalisieren das gleiche Ergebnis.
(xσx)2+(yσy)2=c
Σ

Nun müssten Sie die PDF-Datei der Multivariate innerhalb (oder außerhalb) der Ellipse integrieren und anfordern, dass diese dem gewünschten Quantil entspricht. Angenommen, Ihre Quantile sind nicht die üblichen, sondern im Prinzip elliptisch (dh Sie suchen nach dem HDR (Highest Density Region), wie die Antwort von Tim zeigt). Ich würde Variablen im pdf zu ändern , in den Winkel integrieren und dann für z von 0 nach z2=(x/σx)2+(y/σy)2z0 1-α=c Dazu Ersatz s = - z 2 / 2 :

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

Im Prinzip müssen Sie also nach der Ellipse suchen, die in zentriert ist , wobei die Achse über den Eigenvektoren von Σ und dem effektiven Radius - 2 liegtμΣ-2lnα

(x-μ)TΣ-1(x-μ)=-2lnα
chuse
quelle
4

Sie haben nach multivariaten Normalen gefragt, aber Ihre Frage begann mit der Frage nach dem "Quantil einer multivariaten Verteilung" im Allgemeinen. Aus dem Wortlaut Ihrer Frage und dem vorausgesetzten Beispiel geht hervor, dass Sie an Regionen mit der höchsten Dichte interessiert sind . Sie werden von Hyndman (1996) wie folgt definiert

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

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

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

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


Hyndman, RJ (1996). Berechnen und Zeichnen von Regionen mit der höchsten Dichte. The American Statistician, 50 (2), 120-126.

Tim
quelle
2

-2ln(α)

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

Sie könnten eine Ellipse zeichnen, die den Mahalanobis-Abständen entspricht.

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)

Oder mit Kreisen um 95%, 75% und 50% der Daten

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))
Gänseblümchen
quelle
4
Willkommen auf der Site @ user98114. Können Sie etwas Text bereitstellen, um zu erläutern, was dieser Code tut und wie das OP-Problem behoben wird?
gung - Wiedereinsetzung von Monica