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

14

Я ищу метод для расчета области перекрытия между двумя оценками плотности ядра в 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, которые я мог бы использовать для его реализации.

Магнитка
источник
2
«Количественная оценка фиолетовой области» является проблемой при оценке, а не при проверке гипотез, поэтому вы не можете надеяться «выполнить это с помощью стандартного статистического критерия цитируемости ». Вы противоречите себе. Пожалуйста, уточните, что вы на самом деле хотите. Если все, что вам нужно, это оценка области перекрытия двух KDE, это простой расчет.
Glen_b
@Glen_b спасибо за комментарий, помогли прояснить мое нестатистическое мышление. Я считаю, что область совпадений между KDE - это то, что я ищу - я отредактировал вопрос, чтобы отразить это.
ммк
2
Я был бы очень обеспокоен риском произвола в этом методе. В зависимости от пропускной способности ядра вычисленное перекрытие между любыми двумя наборами данных может быть сделано равным любому выбранному значению в интервале . Пропускная способность по умолчанию не оптимизирована для этой цели и, следовательно, может дать неожиданные, произвольные или противоречивые результаты. Наборы данных с естественными границами (такими как неотрицательные данные или пропорции и т. Д.) Могут привести к нежелательным краевым эффектам. Что делать вместо этого? Начните с причины такого расчета: что означает это «сходство»? (0,1)
whuber
Тот же вопрос появился несколько месяцев спустя, но касался точек пересечения, однако были некоторые важные замечания, которые можно было бы принять во внимание. В упомянутом вопросе речь идет о двух эмпирических распределениях. Я добавляю ссылку, так как этот пост отвечает на это только через оценку плотности ядра и для нормальных распределений. Ссылка ниже, я думаю, распространяется на вопрос о парах эмпирических распределений. stats.stackexchange.com/questions/122857/… - Барнаби 7 часов назад
Барнаби

Ответы:

9

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

1) Поскольку исходные KDE, вероятно, были оценены по некоторой сетке, если сетка одинакова для обоих (или может быть легко сделана одинаковой), упражнение может быть так же просто, как просто взятьмин(К1(Икс),К2(Икс)) в каждой точке, а затем с использованием правила трапеции или даже правила средней точки.

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

2) Вы можете найти точку (или точки) пересечения и интегрировать нижний из двух KDE в каждом интервале, где каждый ниже. На приведенной выше диаграмме вы бы интегрировали синюю кривую слева от пересечения и розовую справа, используя любые доступные вам средства. Это можно сделать по существу точно, рассматривая область под каждым компонентом ядра слева или справа от этой точки отсечения.1часК(Икс-Иксячас)

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

Glen_b - Восстановить Монику
источник
Как рассчитать ошибку, связанную с методом 1 и 2?
olliepower
В обычных обстоятельствах, оба будут незначительными по сравнению с ошибкой в ​​оценках плотности ядра, поэтому я бы не стал сильно беспокоиться. Границы ошибок могут быть рассчитаны на основе трапециевидных методов и других числовых интегрирований, конечно - такие вычисления довольно стандартны - но это бесполезно, учитывая, что KDE имеют большую неопределенность. Метод 2 будет с точностью до накопленной ошибки округления расчетов.
Glen_b
1
Эти методические предложения имеют смысл, большое спасибо за ваш ответ. Я буду работать над реализацией этого в R, но, как новичок, меня заинтересуют предложения о том, как правильно это кодировать.
ММК
10

Для полноты, вот как я закончил делать это в 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, а также с интеграцией.

Магнитка
источник
2
В настоящее время в CRAN существует пакет, overlappingкоторый оценивает область перекрытия 2 (или более) эмпирических распределений. Ознакомьтесь с документацией здесь: rdocumentation.org/packages/overlapping/versions/1.5.0/topics/…
Стефан Авей,
Икс,dИкс,dИкс,d
@mmk вы можете сделать это для 2D плотности?
Нет лжи
4

Во-первых, я могу ошибаться, но я думаю, что ваше решение не сработает в том случае, если есть несколько точек, где пересекаются оценки плотности ядра (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)
С. Венне
источник