Bagaimana cara mendapatkan wilayah elips dari data terdistribusi normal bivariat?

11

Saya memiliki data yang terlihat seperti:

Angka

Saya mencoba menerapkan distribusi normal (estimasi kepadatan kernel berfungsi lebih baik, tetapi saya tidak membutuhkan ketelitian yang luar biasa) dan ini berfungsi dengan baik. Plot kepadatan membuat elips.

Saya perlu mendapatkan fungsi elips untuk memutuskan apakah suatu titik terletak di dalam wilayah elips atau tidak. Bagaimana cara melakukannya?

R atau kode Mathematica disambut.

matejuh
sumber

Jawaban:

18

Corsario memberikan solusi yang baik dalam komentar: gunakan fungsi kerapatan kernel untuk menguji inklusi dalam set level.

Interpretasi lain dari pertanyaan adalah bahwa ia meminta prosedur untuk menguji inklusi dalam elips yang dibuat oleh pendekatan normal bivariat terhadap data. Untuk memulai, mari buat beberapa data yang terlihat seperti ilustrasi dalam pertanyaan:

library(mvtnorm) # References rmvnorm()
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))

Elips ditentukan oleh momen pertama dan kedua data:

center <- apply(p, 2, mean)
sigma <- cov(p)

Rumus ini membutuhkan inversi dari matriks varians-kovarians:

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

Fungsi "tinggi" elips adalah negatif dari logaritma densitas normal bivariat :

ellipse <- function(s,t) {u<-c(s,t)-center; u %*% sigma.inv %*% u / 2}

(Saya telah mengabaikan konstanta aditif yang sama dengan .)log(2πdet(Σ))

Untuk menguji ini , mari kita menggambar beberapa konturnya. Itu membutuhkan menghasilkan kisi-kisi poin dalam arah x dan y:

n <- 50
x <- (0:(n-1)) * (500000/(n-1))
y <- (0:(n-1)) * (50000/(n-1))

Hitung fungsi ketinggian di kisi ini dan plot:

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
contour(x,y,matrix(z,n,n), levels=(0:10), col = terrain.colors(11), add=TRUE)

Plot kontur

Jelas itu berhasil. Oleh karena itu, tes untuk menentukan apakah titik terletak di dalam sebuah kontur elips di tingkat adalah(s,t)c

ellipse(s,t) <= c

Mathematica melakukan pekerjaan dengan cara yang sama: menghitung matriks varians-kovarians data, membalikkannya, membangun ellipsefungsi, dan Anda siap.

whuber
sumber
Terima kasih semuanya, terutama @whuber. Inilah yang saya butuhkan.
matejuh
Btw. apakah ada solusi sederhana untuk kontur estimasi kerapatan kernel? Karena jika saya ingin lebih ketat, data saya terlihat seperti: github.com/matejuh/doschecker_wiki_images/raw/master/… resp. github.com/matejuh/doschecker_wiki_images/raw/master/…
matejuh
Saya tidak dapat menemukan solusi sederhana di R. Pertimbangkan untuk menggunakan fungsi "SmoothKernelDistribution" Mathematica 8.
whuber
2
Apakah level tersebut sesuai dengan tingkat kepercayaan? Saya rasa tidak. Bagaimana saya bisa melakukan itu?
matejuh
Itu membutuhkan pertanyaan baru, karena Anda perlu menentukan apa yang Anda cari kepercayaannya dan - menilai dari plot Anda - ada kekhawatiran tentang apakah elips tersebut merupakan deskripsi data yang memadai di tempat pertama.
whuber
9

Plot mudah dengan ellipse()fungsi mixtoolspaket untuk R:

library(mixtools)
library(mvtnorm) 
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 

masukkan deskripsi gambar di sini

Stéphane Laurent
sumber
5

Pendekatan pertama

Anda dapat mencoba pendekatan ini di Mathematica.

Mari kita buat beberapa data bivariat:

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

Maka kita perlu memuat paket ini:

Needs["MultivariateStatistics`"]

Dan sekarang:

ellPar=EllipsoidQuantile[data, {0.9}]

memberikan output yang mendefinisikan elips kepercayaan 90%. Nilai yang Anda peroleh dari output ini adalah dalam format berikut:

{Ellipsoid[{x1, x2}, {r1, r2}, {{d1, d2}, {d3, d4}}]}

x1 dan x2 menentukan titik di mana elips di tengah, r1 dan r2 menentukan jari-jari semi-sumbu, dan d1, d2, d3 dan d4 menentukan arah penyelarasan.

Anda juga dapat merencanakan ini:

Show[{ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio -> 1],  Graphics[EllipsoidQuantile[data, 0.9]]}]

Bentuk parametrik umum elips adalah:

ell[t_, xc_, yc_, a_, b_, angle_] := {xc + a Cos[t] Cos[angle] - b Sin[t] Sin[angle],
    yc + a Cos[t] Sin[angle] + b Sin[t] Cos[angle]}

Dan Anda dapat memplotnya dengan cara ini:

ParametricPlot[
    ell[t, ellPar[[1, 1, 1]], ellPar[[1, 1, 2]], ellPar[[1, 2, 1]], ellPar[[1, 2, 2]],
    ArcTan[ellPar[[1, 3, 1, 2]]/ellPar[[1, 3, 1, 1]]]], {t, 0, 2 \[Pi]},
    PlotRange -> {{0, 100}, {0, 100}}]

Anda dapat melakukan pemeriksaan berdasarkan informasi geometris murni: jika jarak Euclidean antara pusat elips (ellPar [[1,1]]) dan titik data Anda lebih besar daripada jarak antara pusat elips dan batas elips (jelas, dalam arah yang sama di mana titik Anda berada), maka titik data itu berada di luar elips.

Pendekatan kedua

Pendekatan ini didasarkan pada distribusi kernel yang halus.

Ini adalah beberapa data yang didistribusikan dengan cara yang mirip dengan data Anda:

data1 = RandomVariate[BinormalDistribution[{.3, .7}, {.2, .3}, .8], 500];
data2 = RandomVariate[BinormalDistribution[{.6, .3}, {.4, .15}, .8], 500];
data = Partition[Flatten[Join[{data1, data2}]], 2];

Kami mendapatkan distribusi kernel yang halus pada nilai data ini:

skd = SmoothKernelDistribution[data];

Kami memperoleh hasil numerik untuk setiap titik data:

eval = Table[{data[[i]], PDF[skd, data[[i]]]}, {i, Length[data]}];

Kami memperbaiki ambang batas dan kami memilih semua data yang lebih tinggi dari ambang ini:

threshold = 1.2;
dataIn = Select[eval, #1[[2]] > threshold &][[All, 1]];

Di sini kita mendapatkan data yang berada di luar wilayah:

dataOut = Complement[data, dataIn];

Dan sekarang kita dapat memplot semua data:

Show[ContourPlot[Evaluate@PDF[skd, {x, y}], {x, 0, 1}, {y, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, PlotPoints -> 50],
ListPlot[dataIn, PlotStyle -> Darker[Green]],
ListPlot[dataOut, PlotStyle -> Red]]

Titik berwarna hijau adalah titik di atas ambang batas dan titik berwarna merah adalah titik di bawah ambang batas.

masukkan deskripsi gambar di sini

VLC
sumber
Terima kasih, pendekatan kedua Anda banyak membantu saya dengan distribusi Kernel. Saya programmer, bukan statistik dan saya pemula di Mathmatica dan R jadi saya sangat menghargai bantuan Anda. Dalam pendekatan kedua Anda, jelas bagi saya bagaimana menguji satu titik di mana ia berada. Tetapi bagaimana melakukannya dalam pendekatan pertama? Saya kira saya harus membandingkan poin saya dengan definisi ellipsoid. Bisakah tou berikan bagaimana? Sekarang saya harus berharap bahwa ada definisi yang sama dalam R, karena saya perlu menggunakannya di RinRuby ...
matejuh
@matejuh Saya baru saja menambahkan beberapa baris lagi tentang pendekatan pertama yang mungkin mengarahkan Anda ke solusi.
VLC
2

The ellipsefungsi dalam ellipsepaket untuk R akan menghasilkan elips ini (sebenarnya poligon mendekati elips). Anda bisa menggunakan elips itu.

Yang mungkin sebenarnya lebih mudah adalah menghitung ketinggian kerapatan pada titik Anda dan melihat apakah itu lebih tinggi (di dalam elips) atau lebih rendah (di luar elips) daripada nilai kontur di elips. The ellipseinternal fungsi menggunakan nilai untuk membuat elips, Anda bisa mulai dari sana untuk menemukan tinggi untuk digunakan.χ2

Greg Snow
sumber
1

Saya menemukan jawabannya di: /programming/2397097/how-can-a-data-ellipse-be-superimposed-on-a-ggplot2-scatterplot

#bootstrap
set.seed(101)
n <- 1000
x <- rnorm(n, mean=2)
y <- 1.5 + 0.4*x + rnorm(n)
df <- data.frame(x=x, y=y, group="A")
x <- rnorm(n, mean=2)
y <- 1.5*x + 0.4 + rnorm(n)
df <- rbind(df, data.frame(x=x, y=y, group="B"))

#calculating ellipses
library(ellipse)
df_ell <- data.frame()
for(g in levels(df$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(df[df$group==g,], ellipse(cor(x, y), 
                                         scale=c(sd(x),sd(y)), 
                                         centre=c(mean(x),mean(y))))),group=g))
}
#drawing
library(ggplot2)
p <- ggplot(data=df, aes(x=x, y=y,colour=group)) + geom_point(size=1.5, alpha=.6) +
  geom_path(data=df_ell, aes(x=x, y=y,colour=group), size=1, linetype=2)

masukkan deskripsi gambar di sini

Guy L
sumber