Integracja estymatora gęstości jądra w 2D

Pochodzę z tego pytania, na wypadek, gdyby ktoś chciał pójść szlakiem.

Zasadniczo mam zestaw danych złożony z obiektów, w których do każdego obiektu przypisana jest dana liczba zmierzonych wartości (w tym przypadku dwóch):$\Omega$$N$

Potrzebuję sposobu na określenie prawdopodobieństwa przynależności nowego obiektu do dlatego doradzono mi w tym pytaniu, aby uzyskać gęstość prawdopodobieństwa za pomocą estymatora gęstości jądra, który moim zdaniem już mam.Ω $p[x_p, y_p]$$\Omega$$\hat{f}$

Ponieważ moim celem jest uzyskanie prawdopodobieństwo tego nowego obiektu ( $p[x_p, y_p]$ ) przynależności do tego zestawu danych 2D $\Omega$ , powiedziano mi, aby zintegrować pdf f over " wartości wsparcia dla których gęstość jest mniejszy niż ten, który zaobserwowałeś . Przycisków "obserwuje" gęstość f oceniano w nowego obiektu , tj . Więc muszę rozwiązać równanie:$\hat{f}$$\hat{f}$f ( x p , y s$p$$\hat{f}(x_p, y_p)$

Plik PDF mojego zestawu danych 2D (uzyskany przez moduł stats.gaussian_kde Pythona ) wygląda następująco:

gdzie czerwona kropka reprezentuje nowy obiekt wykreślony nad PDF mojego zestawu danych.$p[x_p, y_p]$

Pytanie brzmi: jak obliczyć powyższą całkę dla granic kiedy pdf wygląda tak?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

Dodaj

Zrobiłem kilka testów, aby zobaczyć, jak dobrze działa metoda Monte Carlo, o której wspomniałem w jednym z komentarzy. Oto co mam:

Wydaje się, że wartości różnią się nieco bardziej dla obszarów o niższej gęstości, przy czym obie szerokości pasma pokazują mniej więcej taką samą zmienność. Największa zmiana w tabeli występuje dla punktu (x, y) = (2.4,1.5) porównując wartość próbki Silvermana 2500 vs 1000, co daje różnicę 0.0126lub ~1.3%. W moim przypadku byłoby to w dużej mierze do przyjęcia.

Edycja : Właśnie zauważyłem, że w 2 wymiarach reguła Scotta jest równoważna z regułą Silvermana zgodnie z podaną tutaj definicją.

Prostym sposobem jest zrasteryzowanie dziedziny integracji i obliczenie dyskretnego przybliżenia do całki.

Na kilka rzeczy należy uważać:

1. Pamiętaj, aby objąć więcej niż zakres punktów: musisz uwzględnić wszystkie lokalizacje, w których oszacowanie gęstości jądra będzie miało jakiekolwiek znaczące wartości. Oznacza to, że musisz zwiększyć zasięg punktów o trzy do czterech razy większą szerokość pasma jądra (dla jądra Gaussa).

2. Wynik różni się nieco w zależności od rozdzielczości rastra. Rozdzielczość musi stanowić niewielki ułamek szerokości pasma. Ponieważ czas obliczeń jest proporcjonalny do liczby komórek w rastrze, wykonanie serii obliczeń przy użyciu mniej dokładnych rozdzielczości nie zajmuje prawie więcej czasu: sprawdź, czy wyniki dla grubszych są zbieżne z wynikiem dla najlepsza rozdzielczość. Jeśli nie są, konieczne może być dokładniejsze rozwiązanie.

Oto ilustracja zestawu danych o 256 punktach:

Punkty są pokazane jako czarne kropki nałożone na dwa oszacowania gęstości jądra. Sześć dużych czerwonych punktów to „sondy”, przy których algorytm jest oceniany. Dokonano tego dla czterech szerokości pasma (domyślnie między 1,8 (pionowo) i 3 (poziomo), 1/2, 1 i 5 jednostek) przy rozdzielczości 1000 na 1000 komórek. Poniższa macierz wykresu rozrzutu pokazuje, jak silnie wyniki zależą od szerokości pasma dla tych sześciu punktów sondy, które pokrywają szeroki zakres gęstości:

Zmiana występuje z dwóch powodów. Oczywiście szacunki gęstości różnią się, wprowadzając jedną formę zmienności. Co ważniejsze, różnice w szacunkach gęstości mogą powodować duże różnice w dowolnym pojedynczym punkcie („sondzie”). Ta ostatnia odmiana jest największa wokół „obrzeży” skupisk punktów o średniej gęstości - dokładnie w tych lokalizacjach, w których to obliczenie prawdopodobnie będzie najczęściej używane.

Wskazuje to na konieczność zachowania znacznej ostrożności przy stosowaniu i interpretowaniu wyników tych obliczeń, ponieważ mogą one być tak wrażliwe na stosunkowo arbitralną decyzję (przepustowość do wykorzystania).

Kod R.

Algorytm ten jest zawarty w pół tuzina linii pierwszej funkcji f. Aby zilustrować jego użycie, reszta kodu generuje poprzednie liczby.

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)

Jeśli masz przyzwoitą liczbę obserwacji, być może nie będziesz musiał wcale dokonywać integracji. Powiedz, że twoim nowym punktem jest . Załóżmy, że masz estymator gęstości ; zsumuj liczbę obserwacji dla których i podziel przez wielkość próby. Daje to przybliżenie do wymaganego prawdopodobieństwa. f x f ( x )< f ( x 0${\bf{x}}_0$$\hat f$${\bf{x}}$$\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$

Zakłada się, że nie jest „zbyt mały”, a wielkość twojej próbki jest wystarczająco duża (i wystarczająco rozłożona), aby dać przyzwoite oszacowanie w regionach o niskiej gęstości. Wydaje się jednak, że 20000 przypadków jest wystarczająco dużych, dla dwóch zmiennych .$\hat f({\bf{x}}_0)$${\bf{x}}$