Wie berechnet man die Überlappung zwischen den empirischen Wahrscheinlichkeitsdichten?

14

Ich suche nach einer Methode zur Berechnung der Überlappungsfläche zwischen zwei Kerndichteschätzungen in R als Maß für die Ähnlichkeit zwischen zwei Stichproben. Um dies zu verdeutlichen, müsste ich im folgenden Beispiel die Fläche des violett überlappenden Bereichs quantifizieren:

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)

Bildbeschreibung hier eingeben

Eine ähnliche Frage wurde hier diskutiert , mit dem Unterschied, dass ich dies eher für willkürliche empirische Daten als für vordefinierte Normalverteilungen tun muss. Das overlapPaket behandelt diese Frage, aber anscheinend nur für Zeitstempeldaten, was bei mir nicht funktioniert. Der Bray-Curtis-Index (wie in vegander vegdist(method="bray")Funktion des Pakets implementiert ) scheint ebenfalls relevant zu sein, jedoch wiederum für etwas andere Daten.

Ich interessiere mich sowohl für den theoretischen Ansatz als auch für die R-Funktionen, mit denen ich ihn implementieren könnte.

mmk
quelle
2
"Quantifizierung der violetten Fläche" ist ein Problem bei der Schätzung und nicht beim Testen von Hypothesen. Sie können also nicht hoffen, "dies mithilfe eines standardmäßigen zitierfähigen statistischen Tests zu erreichen ". Sie widersprechen sich. Bitte klären Sie, was Sie tatsächlich wollen. Wenn Sie lediglich eine Schätzung des Überlappungsbereichs von zwei KDEs wünschen, ist dies eine einfache Berechnung.
Glen_b
@ Glen_b danke für den Kommentar, der mir geholfen hat, mein nicht-statistisches Denken zu verdeutlichen. Ich glaube, der Bereich der Überlappung zwischen KDEs ist in der Tat das, wonach ich suche - ich habe die Frage bearbeitet, um dies widerzuspiegeln.
MMK
2
Ich wäre sehr besorgt über das Risiko der Willkür bei dieser Methode. Abhängig von der Kernelbandbreite kann die berechnete Überlappung zwischen zwei beliebigen Datensätzen so eingestellt werden, dass sie einem beliebigen Wert im Intervall . Die Standardbandbreiten sind nicht für diesen Zweck optimiert und können daher möglicherweise überraschende, willkürliche oder inkonsistente Ergebnisse liefern. Datensätze mit natürlichen Grenzen (wie nicht-negative Daten oder Proportionen usw.) würden außerdem unerwünschte Kanteneffekte hervorrufen. Was ist stattdessen zu tun? Beginnen Sie mit dem Grund für diese Berechnung: Was soll diese "Ähnlichkeit" bedeuten? (0,1)
whuber
Die gleiche Frage tauchte einige Monate später auf, bezog sich jedoch auf Kreuzungspunkte, es gab jedoch einige gültige Hinweise, die berücksichtigt werden konnten. In der genannten Frage geht es um zwei empirische Verteilungen. Ich füge den Link hinzu, da dieser Beitrag dies nur über die Schätzung der Kerneldichte und für Normalverteilungen beantwortet. Der Link unten, denke ich, erstreckt sich auf die Frage nach empirischen Verteilungspaaren. stats.stackexchange.com/questions/122857/… - Barnaby vor 7 Stunden
Barnaby

Antworten:

9

Der Überlappungsbereich von zwei Kerndichteschätzungen kann auf jeden gewünschten Genauigkeitsgrad angenähert werden.

1) Da die ursprünglichen KDEs wahrscheinlich über ein Raster ausgewertet wurden, kann die Übung so einfach sein, als ob man einfach an jedem Punkt und dann unter Verwendung der Trapezregel oder sogar einer Mittelpunktsregel.min(K1(x),K2(x))

Wenn sich die beiden in unterschiedlichen Gittern befinden und nicht einfach im selben Gitter neu berechnet werden können, kann Interpolation verwendet werden.

1hK(xxih)

Die obigen Ausführungen von whuber sollten jedoch klar beachtet werden - dies ist nicht unbedingt eine sehr bedeutsame Sache.

Glen_b - Setzen Sie Monica wieder ein
quelle
Wie berechnen Sie den Fehler in Verbindung mit Methode 1 und Methode 2?
olliepower
Unter normalen Umständen sind beide im Vergleich zu den Fehlern bei der Schätzung der Kerneldichte winzig, sodass ich mir keine Sorgen machen würde. Fehlergrenzen können natürlich mit trapezförmigen Methoden und anderen numerischen Integrationen berechnet werden - solche Berechnungen sind ziemlich üblich -, aber es ist sinnlos, sich Sorgen zu machen, da KDEs große Unsicherheiten aufweisen. Methode 2 ist auf den akkumulierten Rundungsfehler der Berechnungen genau.
Glen_b
1
Diese Methodenvorschläge sind sinnvoll, vielen Dank für Ihre Antwort. Ich werde daran arbeiten, dies in R zu implementieren, aber als Neuling würde ich mich für Vorschläge interessieren, wie dies sauber codiert werden kann.
MMK
10

Der Vollständigkeit halber habe ich Folgendes in R getan:

# 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

Wie bereits erwähnt, ist die KDE-Generation und auch die Integration mit Unsicherheit und Subjektivität verbunden.

mmk
quelle
2
Es gibt jetzt ein Paket auf CRAN, overlappingdas den Bereich der Überlappung von 2 (oder mehr) empirischen Verteilungen schätzt. Schauen Sie sich die Dokumentation hier an: rdocumentation.org/packages/overlapping/versions/1.5.0/topics/…
Stefan Avey
x,dx,dx,d
@mmk können Sie dies für 2D-Dichten tun?
Keine Lüge
4

Erstens könnte ich mich irren, aber ich denke, Ihre Lösung würde nicht funktionieren, wenn es mehrere Punkte gibt, an denen sich die Kernel Density Estimates (KDE) überschneiden. Zweitens, obwohl das overlapPaket für die Verwendung mit Zeitstempeldaten erstellt wurde, können Sie es dennoch zum Schätzen des Überlappungsbereichs von zwei beliebigen KDEs verwenden. Sie müssen Ihre Daten nur so skalieren, dass sie zwischen 0 und 2π liegen.
Zum Beispiel :

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