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)
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 ellipse
fungsi, dan Anda siap.
Plot mudah dengan
ellipse()
fungsimixtools
paket untuk R:sumber
Pendekatan pertama
Anda dapat mencoba pendekatan ini di Mathematica.
Mari kita buat beberapa data bivariat:
Maka kita perlu memuat paket ini:
Dan sekarang:
memberikan output yang mendefinisikan elips kepercayaan 90%. Nilai yang Anda peroleh dari output ini adalah dalam format berikut:
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:
Bentuk parametrik umum elips adalah:
Dan Anda dapat memplotnya dengan cara ini:
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:
Kami mendapatkan distribusi kernel yang halus pada nilai data ini:
Kami memperoleh hasil numerik untuk setiap titik data:
Kami memperbaiki ambang batas dan kami memilih semua data yang lebih tinggi dari ambang ini:
Di sini kita mendapatkan data yang berada di luar wilayah:
Dan sekarang kita dapat memplot semua data:
Titik berwarna hijau adalah titik di atas ambang batas dan titik berwarna merah adalah titik di bawah ambang batas.
sumber
The
ellipse
fungsi dalamellipse
paket 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χ2
ellipse
internal fungsi menggunakan nilai untuk membuat elips, Anda bisa mulai dari sana untuk menemukan tinggi untuk digunakan.sumber
Saya menemukan jawabannya di: /programming/2397097/how-can-a-data-ellipse-be-superimposed-on-a-ggplot2-scatterplot
sumber