A simple way is to rasterize the domain of integration and compute a discrete approximation to the integral.
There are some things to watch out for:
Make sure to cover more than the extent of the points: you need to include all locations where the kernel density estimate will have any appreciable values. This means you need to expand the extent of the points by three to four times the kernel bandwidth (for a Gaussian kernel).
The result will vary somewhat with the resolution of the raster. The resolution needs to be a small fraction of the bandwidth. Because the calculation time is proportional to the number of cells in the raster, it takes almost no extra time to perform a series of calculations using coarser resolutions than the intended one: check that the results for the coarser ones are converging on the result for the finest resolution. If they are not, a finer resolution may be needed.
Here is an illustration for a dataset of 256 points:
The points are shown as black dots superimposed on two kernel density estimates. The six large red points are "probes" at which the algorithm is evaluated. This has been done for four bandwidths (a default between 1.8 (vertically) and 3 (horizontally), 1/2, 1, and 5 units) at a resolution of 1000 by 1000 cells. The following scatterplot matrix shows how strongly the results depend on bandwidth for these six probe points, which cover a wide range of densities:
The variation occurs for two reasons. Obviously the density estimates differ, introducing one form of variation. More importantly, the differences in density estimates can create large differences at any single ("probe") point. The latter variation is greatest around the medium-density "fringes" of clusters of points--exactly those locations where this calculation is likely to be used the most.
This demonstrates the need for substantial caution in using and interpreting the results of these calculations, because they can be so sensitive to a relatively arbitrary decision (the bandwidth to use).
R Code
The algorithm is contained in the half dozen lines of the first function, f. To illustrate its use, the rest of the code generates the preceding figures.
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)
Best Answer
The KDE is a mixture of Normal distributions. Let's look at a single one of them.
The definitions of $P(A)$ and $P(B)$ show their values are invariant under translations and rescalings in the plane, so it suffices to consider the standard Normal distribution with PDF $f$. The inequality
$$f(x,y) \le f(r,s)$$
is equivalent to
$$x^2 + y^2 \ge r^2 + s^2.$$
Introducing polar coordinates $\rho, \theta$ allows the integral to be rewritten
$$P(r,s) = \frac{1}{2\pi}\int_0^{2\pi}\int_\sqrt{r^2+s^2}^\infty \rho \exp(-\rho^2/2) d\rho d\theta= \exp(-(r^2+s^2)/2) = 2\pi f(r,s).$$
Now consider the mixture. Because it is linear,
$$\eqalign{ P(r,s) &= \frac{1}{n}\sum_i 2\pi f((r-x_i)/h, (s-y_i)/h) \\ &= 2\pi h^2\left(\frac{1}{n}\sum_i \frac{1}{h^2} f((r-x_i)/h, (s-y_i)/h)\right) \\ &=2\pi h^2 KDE(r,s). }$$
Indeed, $f$ and $P$ are proportional. The constant of proportionality is $2\pi h^2$.
That such a proportionality relationship between $P$ and $f$ is special can be appreciated by contemplating a simple counterexample. Let $f_1$ have a uniform distribution on a measurable set $A_1$ of unit area and $f_2$ have a uniform distribution on a measurable set $A_2$ which is disjoint from $A_1$ and has area $\mu\gt 1$. Then the mixture with PDF $f=f_1/2 + f_2/2$ has constant value $1/2$ on $A_1$, $1/(2\mu)$ on $A_2$, and is zero elsewhere. There are three cases to consider:
$(r,s)\in A_1$. Here $f(r,s)=1/2$ attains its maximum, whence $P(r,s)=1$. The ratio $f(r,s)/P(r,s) = 1/2$.
$(r,s)\in A_2$. Here $f(r,s)$ is strictly less than $1/2$ but greater than $0$. Thus the region of integration is the complement of $A_1$ and the resulting integral must equal $1/2$. The ratio $f(r,s)/P(r,s) = (1/(2\mu))/(1/2) = 1/\mu$.
Elsewhere, $f$ is zero and the integral $P$ is zero.
Evidently the ratio (where it is defined) is not constant and varies between $1$ and $1/\mu \ne 1$. Although this distribution is not continuous, it can be made so by adding a Normal$(0,\Sigma)$ distribution to it. By making both eigenvalues of $\Sigma$ small, this will change the distribution very little and produce qualitatively the same results--only now the values of the ratio $f/P$ will include all the numbers in the interval $[1,1/\mu]$.
This result also does not generalize to other dimensions. Essentially the same calculation that started this answer shows that $P$ is an incomplete Gamma function and that clearly is not the same as $f$. That two dimensions are special can be appreciated by noting that the integration in $P$ essentially concerns the distances and when those are Normally distributed, the distance function has a $\chi^2(2)$ distribution--which is the exponential distribution. The exponential function is unique in being proportional to its own derivative--whence the integrand $f$ and integral $P$ must be proportional.