Prostym sposobem jest zrasteryzowanie dziedziny integracji i obliczenie dyskretnego przybliżenia do całki.
Na kilka rzeczy należy uważać:
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).
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)
