Jak obliczyć nakładanie się między gęstościami prawdopodobieństwa empirycznego?

14

Szukam metody do obliczenia obszaru nakładania się dwóch oszacowań gęstości jądra w R, jako miary podobieństwa między dwiema próbkami. Aby to wyjaśnić, w poniższym przykładzie musiałbym określić ilościowo obszar pokrywającego się regionu fioletowo:

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)

wprowadź opis zdjęcia tutaj

Omówiono tutaj podobne pytanie , z tą różnicą, że muszę to zrobić dla arbitralnych danych empirycznych, a nie dla predefiniowanych rozkładów normalnych. Te overlapadresy pakietów to pytanie, ale najwyraźniej tylko dla danych datownika, który nie działa dla mnie. Indeks Bray-Curtis (zaimplementowany w funkcji veganpakietu vegdist(method="bray")) również wydaje się istotny, ale znowu dla nieco innych danych.

Interesuje mnie zarówno podejście teoretyczne, jak i funkcje R, które mogę wykorzystać do jego wdrożenia.

mmk
źródło
2
„kwantyfikacja fioletowego obszaru” jest problemem w szacowaniu, a nie w testowaniu hipotez, więc nie możesz mieć nadziei, że „osiągniesz to przy użyciu standardowego cytowanego testu statystycznego ”. Zaprzeczasz sobie. Wyjaśnij, czego naprawdę chcesz. Jeśli wszystko, czego potrzebujesz, to oszacowanie obszaru nakładania się dwóch KDE, to proste obliczenie.
Glen_b
@Glen_b dzięki za komentarz, pomógł wyjaśnić moje niestatystyczne myślenie. Wierzę, że obszar nakładania się KDE jest rzeczywiście tym, czego szukam - zredagowałem to pytanie, aby to odzwierciedlić.
mmk
2
(0,1)
To samo pytanie pojawiło się kilka miesięcy później, ale dotyczyło punktów przecięcia, jednak były pewne ważne notatki, które można wziąć pod uwagę. W zadanym pytaniu chodzi o dwa rozkłady empiryczne. Dodaję link, ponieważ ten post odpowiada na to tylko poprzez oszacowanie gęstości jądra i dla normalnych dystrybucji. Myślę, że poniższy link rozciąga się na pytanie dotyczące par rozkładów empirycznych. stats.stackexchange.com/questions/122857/… - Barnaby 7 godzin temu
Barnaby

Odpowiedzi:

9

Obszar nakładania się dwóch oszacowań gęstości jądra może być przybliżony do dowolnego pożądanego stopnia dokładności.

min(K1(x),K2(x))

Jeśli oba są na różnych siatkach i nie można ich łatwo przeliczyć na tej samej siatce, można zastosować interpolację.

2) Możesz znaleźć punkt (punkty) przecięcia i zintegrować dolne z dwóch KDE w każdym przedziale, gdzie każdy jest niższy. Na powyższym diagramie zintegrowałbyś niebieską krzywą po lewej stronie skrzyżowania i różową po prawej w dowolny sposób, jaki chcesz / masz do dyspozycji. Można to zrobić zasadniczo dokładnie, biorąc pod uwagę obszar pod każdym z nich1hK(xxih)

jednak pamiętać o powyższych uwagach Whubera - niekoniecznie jest to bardzo znacząca rzecz.

Glen_b - Przywróć Monikę
źródło
Jak obliczyć błąd związany z metodą pierwszą i metodą 2?
olliepower
W normalnych okolicznościach oba będą maleńkie w porównaniu z błędem w szacunkach gęstości jądra, więc nie martwiłbym się zbytnio. Granice błędów można obliczyć metodami trapezoidalnymi i innymi liczbowymi całkami oczywiście - takie obliczenia są dość standardowe - ale nie ma sensu martwić, biorąc pod uwagę, że KDE mają duże niepewności. Metoda 2 będzie dokładna do skumulowanego błędu zaokrąglenia obliczeń.
Glen_b
1
Te sugestie metodologiczne mają sens, dziękuję bardzo za odpowiedź. Będę pracował nad implementacją tego w języku R, ale jako nowicjusz byłbym zainteresowany sugestiami, jak kodować to w czysty sposób.
mmk
10

Dla kompletności, oto jak skończyłem robić to w 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

Jak wspomniano, generowanie KDE, a także integracja, wiąże się z nieodłączną niepewnością i podmiotowością.

mmk
źródło
2
Dostępny jest teraz pakiet o nazwie CRAN, overlappingktóry szacuje obszar nakładania się 2 (lub więcej) rozkładów empirycznych. Sprawdź dokumentację tutaj: rdocumentation.org/packages/overlapping/versions/1.5.0/topics/…
Stefan Avey
x,dx,dx,re
@mmk czy możesz to zrobić dla gęstości 2D?
OverFlow Police
4

Po pierwsze, mogę się mylić, ale myślę, że twoje rozwiązanie nie zadziałałoby w przypadku, gdy istnieje wiele punktów, w których przecinają się szacunki gęstości jądra (KDE). Po drugie, chociaż overlappakiet został stworzony do użytku z danymi znaczników czasu, nadal możesz go użyć do oszacowania obszaru nakładania się dowolnych dwóch KDE. Musisz po prostu przeskalować dane, aby zawierały się w przedziale od 0 do 2π.
Na przykład :

# 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)
S. Venne
źródło