Как рассчитать перекрытие между эмпирическими плотностями вероятности?

Я ищу метод для расчета области перекрытия между двумя оценками плотности ядра в R, как мера сходства между двумя выборками. Чтобы уточнить, в следующем примере мне нужно было бы количественно определить площадь области пурпурного перекрытия:

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)

Подобный вопрос обсуждался здесь , с той разницей, что мне нужно сделать это для произвольных эмпирических данных, а не для предопределенных нормальных распределений. В overlapпакете рассматривает этот вопрос, но , видимо , только для данных временной метки, которая не работает для меня. Индекс Брей-Кертиса (как реализовано в функции veganпакета vegdist(method="bray")) также представляется актуальным, но опять же для несколько других данных.

Меня интересует как теоретический подход, так и функции R, которые я мог бы использовать для его реализации.

Область перекрытия двух оценок плотности ядра может быть аппроксимирована с любой желаемой степенью точности.

1) Поскольку исходные KDE, вероятно, были оценены по некоторой сетке, если сетка одинакова для обоих (или может быть легко сделана одинаковой), упражнение может быть так же просто, как просто взять$\min(K_1(x),K_2(x))$ в каждой точке, а затем с использованием правила трапеции или даже правила средней точки.

Если они находятся на разных сетках и не могут быть легко пересчитаны на одной сетке, можно использовать интерполяцию.

2) Вы можете найти точку (или точки) пересечения и интегрировать нижний из двух KDE в каждом интервале, где каждый ниже. На приведенной выше диаграмме вы бы интегрировали синюю кривую слева от пересечения и розовую справа, используя любые доступные вам средства. Это можно сделать по существу точно, рассматривая область под каждым компонентом ядра слева или справа от этой точки отсечения.$\frac{1}{h}K(\frac{x-x_i}{h})$

Тем не менее , вышеприведенные комментарии должны быть четко учтены - это не обязательно очень важно.

Для полноты, вот как я закончил делать это в 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

Как уже отмечалось, существует неопределенность и субъективность, связанные с генерацией KDE, а также с интеграцией.

Во-первых, я могу ошибаться, но я думаю, что ваше решение не сработает в том случае, если есть несколько точек, где пересекаются оценки плотности ядра (KDE). Во-вторых, хотя overlapпакет был создан для использования с данными временной метки, вы все равно можете использовать его для оценки области перекрытия любых двух KDE. Вам просто нужно изменить масштаб ваших данных, чтобы они варьировались от 0 до 2π.
Например :

# 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)