Mengintegrasikan penduga kepadatan kernel dalam 2D

Saya datang dari pertanyaan ini kalau-kalau ada yang ingin mengikuti jejak.

Pada dasarnya saya memiliki kumpulan data $\Omega$ terdiri dari objek $N$ mana setiap objek memiliki sejumlah nilai terukur yang melekat padanya (dua dalam kasus ini):

Saya perlu cara untuk menentukan probabilitas baru objek $p[x_p, y_p]$ dari milik $\Omega$ jadi saya disarankan dalam pertanyaan itu untuk mendapatkan kepadatan probabilitas f melalui estimator densitas kernel, yang saya percaya saya sudah punya .$\hat{f}$

Karena tujuan saya adalah untuk mendapatkan probabilitas objek baru ini ( $p[x_p, y_p]$ ) dari milik ini 2D kumpulan data $\Omega$ , saya diberitahu untuk mengintegrasikan pdf f lebih " nilai-nilai dari dukungan yang densitas kurang dari yang Anda amati ". The "mengamati" kepadatan f dievaluasi dalam objek baru p , yaitu: f ( x p , y p ) . Jadi saya perlu memecahkan persamaan:$\hat{f}$$\hat{f}$$p$$\hat{f}(x_p, y_p)$

PDF dari kumpulan data 2D saya (diperoleh melalui modul stats.gaussian_kde python ) terlihat seperti ini:

di mana titik merah mewakili objek baru diplot pada PDF dari kumpulan data saya.$p[x_p, y_p]$

Jadi pertanyaannya adalah: bagaimana saya bisa menghitung integral di atas untuk batas ketika pdf terlihat seperti itu?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

Menambahkan

Saya melakukan beberapa tes untuk melihat seberapa baik metode Monte Carlo yang saya sebutkan di salah satu komentar bekerja. Inilah yang saya dapat:

Nilai-nilai tampak sedikit lebih bervariasi untuk daerah kepadatan lebih rendah dengan kedua bandwidth menunjukkan variasi yang kurang lebih sama. Variasi terbesar dalam tabel terjadi untuk titik (x, y) = (2,4,1.5) membandingkan nilai sampel 2500,000 1000 Silverman, yang memberikan perbedaan 0.0126atau ~1.3%. Dalam kasus saya ini sebagian besar dapat diterima.

Sunting : Saya baru saja memperhatikan bahwa dalam 2 dimensi aturan Scott sama dengan Silverman sesuai dengan definisi yang diberikan di sini .

Cara sederhana adalah merasterisasi domain integrasi dan menghitung pendekatan diskrit ke integral.

Ada beberapa hal yang harus diperhatikan:

1. Pastikan untuk mencakup lebih dari batas poin: Anda harus menyertakan semua lokasi di mana estimasi kepadatan kernel akan memiliki nilai yang cukup besar. Ini berarti Anda perlu memperluas titik dengan tiga hingga empat kali bandwidth kernel (untuk kernel Gaussian).

2. Hasilnya akan sedikit berbeda dengan resolusi raster. Resolusi perlu sebagian kecil dari bandwidth. Karena waktu perhitungan proporsional dengan jumlah sel dalam raster, hampir tidak ada waktu ekstra untuk melakukan serangkaian perhitungan menggunakan resolusi yang lebih kasar daripada yang dimaksud: periksa bahwa hasil untuk yang lebih kasar konvergen pada hasil untuk resolusi terbaik. Jika tidak, resolusi yang lebih baik mungkin diperlukan.

Berikut ini adalah ilustrasi untuk dataset 256 poin:

Poin ditampilkan sebagai titik-titik hitam yang ditumpangkan pada dua perkiraan kepadatan kernel. Enam poin merah besar adalah "probe" di mana algoritma dievaluasi. Ini telah dilakukan untuk empat bandwidth (default antara 1,8 (vertikal) dan 3 (horizontal), 1/2, 1, dan 5 unit) pada resolusi 1000 oleh 1000 sel. Matriks scatterplot berikut menunjukkan seberapa kuat hasil bergantung pada bandwidth untuk enam titik penyelidikan ini, yang mencakup berbagai kepadatan:

Variasi terjadi karena dua alasan. Jelas perkiraan kepadatan berbeda, memperkenalkan satu bentuk variasi. Lebih penting lagi, perbedaan dalam estimasi kepadatan dapat menciptakan perbedaan besar pada setiap titik ("penyelidikan"). Variasi yang terakhir adalah yang terbesar di sekitar "pinggiran" kepadatan menengah dari kelompok titik - tepatnya lokasi tempat perhitungan ini paling sering digunakan.

Ini menunjukkan perlunya kehati-hatian yang substansial dalam menggunakan dan menafsirkan hasil perhitungan ini, karena mereka bisa sangat sensitif terhadap keputusan yang relatif sewenang-wenang (bandwidth untuk digunakan).

Kode R

Algoritma ini terkandung dalam setengah lusin baris fungsi pertama f,. Untuk menggambarkan penggunaannya, sisa kode menghasilkan angka-angka sebelumnya.

library(MASS)     # kde2d
library(spatstat) # im class
f <- function(xy, n, x, y, ...) {
  #
  # Estimate the total where the density does not exceed that at (x,y).
  #
  # `xy` is a 2 by ... array of points.
  # `n`  specifies the numbers of rows and columns to use.
  # `x` and `y` are coordinates of "probe" points.
  # `...` is passed on to `kde2d`.
  #
  # Returns a list:
  #   image:    a raster of the kernel density
  #   integral: the estimates at the probe points.
  #   density:  the estimated densities at the probe points.
  #
  xy.kde <- kde2d(xy[1,], xy[2,], n=n, ...)
  xy.im <- im(t(xy.kde$z), xcol=xy.kde$x, yrow=xy.kde$y) # Allows interpolation $
  z <- interp.im(xy.im, x, y)                            # Densities at the probe points
  c.0 <- sum(xy.kde$z)                                   # Normalization factor $
  i <- sapply(z, function(a) sum(xy.kde$z[xy.kde$z < a])) / c.0
  return(list(image=xy.im, integral=i, density=z))
}
#
# Generate data.
#
n <- 256
set.seed(17)
xy <- matrix(c(rnorm(k <- ceiling(2*n * 0.8), mean=c(6,3), sd=c(3/2, 1)), 
               rnorm(2*n-k, mean=c(2,6), sd=1/2)), nrow=2)
#
# Example of using `f`.
#
y.probe <- 1:6
x.probe <- rep(6, length(y.probe))
lims <- c(min(xy[1,])-15, max(xy[1,])+15, min(xy[2,])-15, max(xy[2,]+15))
ex <- f(xy, 200, x.probe, y.probe, lim=lims)
ex$density; ex$integral
#
# Compare the effects of raster resolution and bandwidth.
#
res <- c(8, 40, 200, 1000)
system.time(
  est.0 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, lims=lims)$integral))
est.0
system.time(
  est.1 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1, lims=lims)$integral))
est.1
system.time(
  est.2 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1/2, lims=lims)$integral))
est.2
system.time(
  est.3 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=5, lims=lims)$integral))
est.3
results <- data.frame(Default=est.0[,4], Hp5=est.2[,4], 
                      H1=est.1[,4], H5=est.3[,4])
#
# Compare the integrals at the highest resolution.
#
par(mfrow=c(1,1))
panel <- function(x, y, ...) {
  points(x, y)
  abline(c(0,1), col="Red")
}
pairs(results, lower.panel=panel)
#
# Display two of the density estimates, the data, and the probe points.
#
par(mfrow=c(1,2))
xy.im <- f(xy, 200, x.probe, y.probe, h=0.5)$image
plot(xy.im, main="Bandwidth=1/2", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

xy.im <- f(xy, 200, x.probe, y.probe, h=5)$image
plot(xy.im, main="Bandwidth=5", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

${\bf{x}}_0$$\hat f$${\bf{x}}$$\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$

$\hat f({\bf{x}}_0)$${\bf{x}}$