Bagaimana cara menghitung tumpang tindih antara kepadatan probabilitas empiris?

Saya mencari metode untuk menghitung luas tumpang tindih antara dua perkiraan kepadatan kernel di R, sebagai ukuran kesamaan antara dua sampel. Untuk memperjelas, dalam contoh berikut, saya perlu mengukur luas wilayah keunguan yang tumpang tindih:

library(ggplot2)
set.seed(1234)
d <- data.frame(variable=c(rep("a", 50), rep("b", 30)), value=c(rnorm(50), runif(30, 0, 3)))
ggplot(d, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA)

Pertanyaan serupa didiskusikan di sini , perbedaannya adalah bahwa saya perlu melakukan ini untuk data empiris yang sewenang-wenang daripada distribusi normal yang telah ditentukan. The overlapalamat paket pertanyaan ini, tapi rupanya hanya untuk data timestamp, yang tidak bekerja untuk saya. Indeks Bray-Curtis (sebagaimana diterapkan dalam fungsi veganpaket vegdist(method="bray")) juga tampaknya relevan tetapi sekali lagi untuk data yang agak berbeda.

Saya tertarik pada pendekatan teoretis dan fungsi R yang mungkin saya terapkan untuk mengimplementasikannya.

Luas tumpang tindih dari dua estimasi kepadatan kernel dapat diperkirakan hingga tingkat akurasi yang diinginkan.

$\min(K_1(x),K_2(x))$

Jika keduanya berada di grid yang berbeda dan tidak dapat dengan mudah dihitung ulang di grid yang sama, interpolasi dapat digunakan.

$\frac{1}{h}K(\frac{x-x_i}{h})$

Namun , komentar whuber di atas harus diingat dengan jelas - ini belum tentu hal yang sangat berarti untuk dilakukan.

Demi kelengkapan, inilah cara saya akhirnya melakukan ini di R:

# simulate two samples
a <- rnorm(100)
b <- rnorm(100, 2)

# define limits of a common grid, adding a buffer so that tails aren't cut off
lower <- min(c(a, b)) - 1 
upper <- max(c(a, b)) + 1

# generate kernel densities
da <- density(a, from=lower, to=upper)
db <- density(b, from=lower, to=upper)
d <- data.frame(x=da$x, a=da$y, b=db$y)

# calculate intersection densities
d$w <- pmin(d$a, d$b)

# integrate areas under curves
library(sfsmisc)
total <- integrate.xy(d$x, d$a) + integrate.xy(d$x, d$b)
intersection <- integrate.xy(d$x, d$w)

# compute overlap coefficient
overlap <- 2 * intersection / total

Sebagaimana dicatat, ada ketidakpastian dan subjektivitas yang melekat yang terlibat dalam generasi KDE dan juga dalam integrasi.

Pertama, saya mungkin salah tetapi saya pikir solusi Anda tidak akan bekerja jika ada beberapa titik di mana Kernel Density Estimates (KDE) berpotongan. Kedua, meskipun overlappaket itu dibuat untuk digunakan dengan data timestamp, Anda masih dapat menggunakannya untuk memperkirakan area tumpang tindih dari dua KDE. Anda hanya perlu mengubah skala data Anda sehingga berkisar dari 0 hingga 2π.
Sebagai contoh :

# simulate two sample    
 a <- rnorm(100)
 b <- rnorm(100, 2)

# To use overplapTrue(){overlap} the scale must be in radian (i.e. 0 to 2pi)
# To keep the *relative* value of a and b the same, combine a and b in the
# same dataframe before rescaling. You'll need to load the ‘scales‘ library.
# But first add a "Source" column to be able to distinguish between a and b
# after they are combined.
 a = data.frame( value = a, Source = "a" )
 b = data.frame( value = b, Source = "b" )
 d = rbind(a, b)
 library(scales) 
 d$value <- rescale( d$value, to = c(0,2*pi) )

# Now you can created the rescaled a and b vectors
 a <- d[d$Source == "a", 1]
 b <- d[d$Source == "b", 1]

# You can then calculate the area of overlap as you did previously.
# It should give almost exactly the same answers.
# Or you can use either the overlapTrue() and overlapEst() function 
# provided with the overlap packages. 
# Note that with these function the KDE are fitted using von Mises kernel.
 library(overlap)
  # Using overlapTrue():
   # define limits of a common grid, adding a buffer so that tails aren't cut off
     lower <- min(d$value)-1 
     upper <- max(d$value)+1
   # generate kernel densities
     da <- density(a, from=lower, to=upper, adjust = 1)
     db <- density(b, from=lower, to=upper, adjust = 1)
   # Compute overlap coefficient
     overlapTrue(da$y,db$y)


  # Using overlapEst():            
    overlapEst(a, b, kmax = 3, adjust=c(0.8, 1, 4), n.grid = 500)

# You can also plot the two KDEs and the region of overlap using overlapPlot()
# but sadly I haven't found a way of changing the x scale so that the scale 
# range correspond to the initial x value and not the rescaled value.
# You can only change the maximum value of the scale using the xscale argument 
# (i.e. it always range from 0 to n, where n is set with xscale = n).
# So if some of your data take negative value, you're probably better off with
# a different plotting method. You can change the x label with the xlab
# argument.  
  overlapPlot(a, b, xscale = 10, xlab= "x metrics", rug=T)