Bagaimana menentukan kuantil (isoline?) Dari distribusi normal multivariat

24

masukkan deskripsi gambar di sini

Saya tertarik pada bagaimana seseorang dapat menghitung kuantil dari distribusi multivariat. Dalam gambar, saya telah menggambar 5% dan 95% kuantil dari distribusi normal univariat yang diberikan (kiri). Untuk distribusi normal multivariat yang tepat, saya membayangkan analog akan menjadi isoline yang mengelilingi basis fungsi kerapatan. Di bawah ini adalah contoh dari upaya saya untuk menghitung ini menggunakan paket mvtnorm- tetapi tidak berhasil. Saya kira ini bisa dilakukan dengan menghitung kontur hasil dari fungsi kepadatan multivarian, tetapi saya bertanya-tanya apakah ada alternatif lain ( misalnya , analog qnorm). Terima kasih atas bantuan Anda.

Contoh:

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 di dalam kotak
sumber
3
Sebuah Mathematica solusi diberikan (dan diilustrasikan untuk kasus 3D) di mathematica.stackexchange.com/questions/21396/... . Ia mengakui bahwa level kontur diberikan oleh distribusi chi-squared.
whuber
@whuber - bisakah Anda menunjukkan apa yang Anda maksud dengan "... ellipsoid kepercayaan adalah kontur dari kebalikan dari matriks kovarians"? Tepuk tangan.
Marc di dalam kotak
2
Ini paling mudah dilihat dalam satu dimensi, di mana "matriks kovarians" (untuk distribusi sampel) adalah angka , jadi kebalikannya adalah 1 / s 2 , dianggap sebagai peta kuadrat pada R 1 melalui x x 2 / s 2 . Sebuah kontur pada tingkat λ menurut definisi adalah himpunan x yang x 2 / s 2 = λ ; yaitu, x 2 = λ s 2 atau ekuivalen x = ± s21/s2R1xx2/s2λxx2/s2=λx2=λs2. Ketikaλadalah1-αkuantil daridistribusiχ2(1),x=±λsλ1αχ2(1) adalah1-αkuantil daridistribusit(1), di mana kita memulihkan batas kepercayaan yang biasa±t 1 - α ; 1 s. λ1αt(1)±t1α;1s
Whuber
Anda dapat menggunakan rumus pertama dalam jawaban ini dengan memilih dalam ( 0 , 1 ) untuk mendapatkan ellipse S α (garis putus-putus merah pada plot Anda) untuk xR 2α(0,1)SαxR2
pengguna603

Jawaban:

25

Garis kontur adalah ellipsoid. Alasannya adalah karena Anda harus melihat argumen eksponensial, dalam pdf distribusi normal multivarian: isolin akan menjadi garis dengan argumen yang sama. Maka Anda mendapatkan mana Σ adalah matriks kovarians. Persisnya persamaan elips; dalam kasus paling sederhana, μ = ( 0 , 0 ) dan Σ adalah diagonal, sehingga Anda mendapatkan ( x

(x-μ)TΣ-1(x-μ)=c
Σμ=(0,0)Σ JikaΣtidak diagonal, mendiagonalisasi Anda mendapatkan hasil yang sama.
(xσx)2+(yσy)2=c
Σ

Sekarang, Anda harus mengintegrasikan pdf dari multivariat di dalam (atau di luar) elips dan meminta ini sama dengan jumlah yang Anda inginkan. Katakanlah bahwa kuantil Anda bukan yang biasa, tetapi pada prinsipnya berbentuk elips (yaitu Anda mencari Wilayah Kepadatan Tertinggi, HDR, seperti yang ditunjukkan oleh jawaban Tim). Saya akan mengubah variabel dalam pdf ke , berintegrasi dalam sudut dan kemudian untuk z dari 0 ke z2=(x/σx)2+(y/σy)2z0 1-α=c Kemudian Anda pengganti 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)

Jadi pada prinsipnya, Anda harus mencari elips yang berpusat di , dengan sumbu di atas vektor eigen Σ dan radius efektif - 2 ln α : ( x - μ ) T Σ - 1 ( x - μ ) = - 2 ln αμΣ-2dalamα

(x-μ)TΣ-1(x-μ)=-2dalamα
chuse
sumber
4

Anda bertanya tentang multivarian normal, tetapi memulai pertanyaan Anda dengan bertanya tentang "kuantil distribusi multivarian" secara umum. Dari kata-kata pertanyaan Anda dan contoh yang diberikan tampaknya Anda tertarik pada daerah dengan kepadatan tertinggi . Mereka didefinisikan oleh Hyndman (1996) sebagai berikut

Biarkan menjadi fungsi kepadatan dari variabel acak X . Maka 100 ( 1 - α ) % HDR adalah subset R ( f α ) dari ruang sampel X sedemikian rupa sehinggaf(z)X100(1-α)%R(fα)X

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

di mana adalah konstanta terbesar sehingga Pr ( X R ( f α ) ) 1 - a .fαPr(XR(fα))1-Sebuah

HDR dapat diperoleh dengan integrasi tetapi, seperti dijelaskan oleh Hyndman, Anda dapat melakukannya menggunakan metode numerik yang lebih sederhana. Jika Y=f(x) , maka Anda dapat memperoleh sehingga Pr ( f ( x ) f α ) 1 - α hanya dengan mengambil α kuantil dari Y . Hal ini dapat diperkirakan dengan menggunakan quantiles sampel dari serangkaian pengamatan y 1 , . . . , Y mfαPr(f(x)fα)1-ααYy1,...,ym. Metode ini berlaku bahkan jika kita tidak tahu , tetapi hanya memiliki satu set pengamatan iid. Metode ini akan bekerja juga untuk distribusi multimoda.f(x)


Hyndman, RJ (1996). Komputasi dan grafik daerah kepadatan tertinggi. The American Statistician, 50 (2), 120-126.

Tim
sumber
2

-2dalam(α)

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

Anda bisa menggambar elips yang sesuai dengan jarak 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)

Atau dengan lingkaran sekitar 95%, 75%, dan 50% data

drawMahal(x,center=x.mcd$center,covariance=x.mcd$cov,quantile=c(0.95,.75,.5))
bunga aster
sumber
4
Selamat datang di situs @ user98114. Bisakah Anda memberikan beberapa teks untuk menjelaskan apa yang dilakukan kode ini & bagaimana menyelesaikan masalah OP?
gung - Reinstate Monica