Kernel-Dichteschätzer in 2D integrieren

Ich komme von dieser Frage, falls jemand der Spur folgen möchte.

Grundsätzlich habe ich einen Datensatz $\Omega$ bestehend aus $N$ Objekten, an die an jedem Objekt eine bestimmte Anzahl von Messwerten angehängt ist (in diesem Fall zwei):

Ich brauche einen Weg , um die Wahrscheinlichkeit einen bestimmen neuen Objekt der Zugehörigkeit zu Ω so I in dieser Frage wurde geraten , eine Wahrscheinlichkeitsdichte zu erhalten $p[x_p, y_p]$$\Omega$$\hat{f}$ durch einen Kerndichteschätzer, die ich glaube , ich habe schon .

Da mein Ziel , die Wahrscheinlichkeit dieses neuen Objekts (zu erhalten , ist ) der Zugehörigkeit zu diesem 2D - Datensatz Ω , wurde ich gesagt , das PDF zu integrieren f über " Werte der Unterstützung der Dichte , für die ist kleiner als der, den Sie beobachtet haben ". Die "beobachtete" Dichte ist f in dem neuen Objekt ausgewertet p , das heißt: f ( x p , y p ) . Also muss ich die Gleichung lösen:$p[x_p, y_p]$$\Omega$$\hat{f}$$\hat{f}$$p$$\hat{f}(x_p, y_p)$

Das PDF meines 2D-Datensatzes (erhalten über das Modul stats.gaussian_kde von Python ) sieht folgendermaßen aus:

Dabei steht der rote Punkt für das neue Objekt , das über der PDF-Datei meines Datensatzes aufgezeichnet ist.$p[x_p, y_p]$

Die Frage ist also: Wie kann ich berechnen die obige Integral für die Grenzen , wenn die pdf sieht aus wie das?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

Hinzufügen

Ich habe einige Tests durchgeführt, um zu sehen, wie gut die in einem der Kommentare erwähnte Monte-Carlo-Methode funktioniert. Das habe ich bekommen:

Die Werte scheinen für Bereiche mit geringerer Dichte etwas stärker zu variieren, wobei beide Bandbreiten mehr oder weniger die gleiche Variation aufweisen. Die größte Variation in der Tabelle tritt für den Punkt (x, y) = (2.4, 1.5) auf, wenn der Wert von 2500 gegen 1000 von Silverman verglichen wird, was eine Differenz von 0.0126oder ergibt ~1.3%. In meinem Fall wäre dies weitgehend akzeptabel.

Edit : Mir ist gerade aufgefallen, dass Scotts Regel in 2 Dimensionen Silvermans gemäß der hier angegebenen Definition entspricht.

Ein einfacher Weg besteht darin, den Bereich der Integration zu rastern und eine diskrete Approximation an das Integral zu berechnen.

Es gibt einige Dinge zu beachten:

1. Stellen Sie sicher, dass Sie mehr als die Ausdehnung der Punkte abdecken: Sie müssen alle Stellen einschließen, an denen die Schätzung der Kerneldichte nennenswerte Werte aufweist. Dies bedeutet, dass Sie die Ausdehnung der Punkte um das Drei- bis Vierfache der Kernelbandbreite erweitern müssen (für einen Gaußschen Kernel).

2. Das Ergebnis hängt von der Auflösung des Rasters ab. Die Auflösung muss einen kleinen Bruchteil der Bandbreite ausmachen. Da die Berechnungszeit proportional zur Anzahl der Zellen im Raster ist, dauert es fast nicht länger, eine Reihe von Berechnungen mit einer gröberen Auflösung als der beabsichtigten durchzuführen: Überprüfen Sie, ob die Ergebnisse für die gröberen mit dem Ergebnis für das übereinstimmen feinste Auflösung. Ist dies nicht der Fall, ist möglicherweise eine feinere Auflösung erforderlich.

Hier ist eine Illustration für einen Datensatz von 256 Punkten:

Die Punkte werden als schwarze Punkte angezeigt, die zwei Kernel-Dichteschätzungen überlagert sind. Die sechs großen roten Punkte sind "Sonden", an denen der Algorithmus ausgewertet wird. Dies wurde für vier Bandbreiten (ein Standardwert zwischen 1,8 (vertikal) und 3 (horizontal), 1/2, 1 und 5 Einheiten) bei einer Auflösung von 1000 mal 1000 Zellen durchgeführt. Die folgende Streudiagramm-Matrix zeigt, wie stark die Ergebnisse von der Bandbreite dieser sechs Prüfpunkte abhängen, die einen weiten Bereich von Dichten abdecken:

Die Variation tritt aus zwei Gründen auf. Offensichtlich unterscheiden sich die Dichteschätzungen, was eine Form der Variation einführt. Noch wichtiger ist, dass die Unterschiede bei den Dichteschätzungen zu großen Unterschieden an jedem einzelnen Punkt ("Sondenpunkt") führen können. Die letztere Variation ist um die "Ränder" von Punkthaufen mittlerer Dichte am größten - genau an den Stellen, an denen diese Berechnung wahrscheinlich am häufigsten verwendet wird.

Dies zeigt, dass bei der Verwendung und Interpretation der Ergebnisse dieser Berechnungen große Vorsicht geboten ist, da sie für eine relativ willkürliche Entscheidung (die zu verwendende Bandbreite) so empfindlich sein können.

R-Code

Der Algorithmus ist in einem halben Dutzend Zeilen der ersten Funktion enthalten f. Zur Veranschaulichung der Verwendung werden im Rest des Codes die vorhergehenden Abbildungen generiert.

library(MASS)     # kde2d
library(spatstat) # im class
f <- function(xy, n, x, y, ...) {
  #
  # Estimate the total where the density does not exceed that at (x,y).
  #
  # `xy` is a 2 by ... array of points.
  # `n`  specifies the numbers of rows and columns to use.
  # `x` and `y` are coordinates of "probe" points.
  # `...` is passed on to `kde2d`.
  #
  # Returns a list:
  #   image:    a raster of the kernel density
  #   integral: the estimates at the probe points.
  #   density:  the estimated densities at the probe points.
  #
  xy.kde <- kde2d(xy[1,], xy[2,], n=n, ...)
  xy.im <- im(t(xy.kde$z), xcol=xy.kde$x, yrow=xy.kde$y) # Allows interpolation $
  z <- interp.im(xy.im, x, y)                            # Densities at the probe points
  c.0 <- sum(xy.kde$z)                                   # Normalization factor $
  i <- sapply(z, function(a) sum(xy.kde$z[xy.kde$z < a])) / c.0
  return(list(image=xy.im, integral=i, density=z))
}
#
# Generate data.
#
n <- 256
set.seed(17)
xy <- matrix(c(rnorm(k <- ceiling(2*n * 0.8), mean=c(6,3), sd=c(3/2, 1)), 
               rnorm(2*n-k, mean=c(2,6), sd=1/2)), nrow=2)
#
# Example of using `f`.
#
y.probe <- 1:6
x.probe <- rep(6, length(y.probe))
lims <- c(min(xy[1,])-15, max(xy[1,])+15, min(xy[2,])-15, max(xy[2,]+15))
ex <- f(xy, 200, x.probe, y.probe, lim=lims)
ex$density; ex$integral
#
# Compare the effects of raster resolution and bandwidth.
#
res <- c(8, 40, 200, 1000)
system.time(
  est.0 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, lims=lims)$integral))
est.0
system.time(
  est.1 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1, lims=lims)$integral))
est.1
system.time(
  est.2 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=1/2, lims=lims)$integral))
est.2
system.time(
  est.3 <- sapply(res, 
           function(i) f(xy, i, x.probe, y.probe, h=5, lims=lims)$integral))
est.3
results <- data.frame(Default=est.0[,4], Hp5=est.2[,4], 
                      H1=est.1[,4], H5=est.3[,4])
#
# Compare the integrals at the highest resolution.
#
par(mfrow=c(1,1))
panel <- function(x, y, ...) {
  points(x, y)
  abline(c(0,1), col="Red")
}
pairs(results, lower.panel=panel)
#
# Display two of the density estimates, the data, and the probe points.
#
par(mfrow=c(1,2))
xy.im <- f(xy, 200, x.probe, y.probe, h=0.5)$image
plot(xy.im, main="Bandwidth=1/2", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

xy.im <- f(xy, 200, x.probe, y.probe, h=5)$image
plot(xy.im, main="Bandwidth=5", col=terrain.colors(256))
points(t(xy), pch=".", col="Black")
points(x.probe, y.probe, pch=19, col="Red", cex=.5)

Wenn Sie eine anständige Anzahl von Beobachtungen haben, müssen Sie möglicherweise überhaupt keine Integration durchführen. Sagen Sie, Ihr neuer Punkt ist${\bf{x}}_0$$\hat f$${\bf{x}}$$\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$

$\hat f({\bf{x}}_0)$${\bf{x}}$