Estimador de densidad de núcleo integrado en 2D

Question

Vengo de esta pregunta por si alguien quiere seguir el rastro.

Básicamente tengo un conjunto de datos $\Omega$ compuesto por $N$ donde cada objeto tiene un número determinado de valores medidos asociados (dos en este caso):

$$\Omega = o_1[x_1, y_1], o_2[x_2, y_2], ..., o_N[x_N, y_N]$$

Necesito una forma de determinar la probabilidad de un nuevo objeto $p[x_p, y_p]$ de pertenecer a $\Omega$ por lo que se me aconsejó en esa pregunta obtener una densidad de probabilidad $\hat{f}$ a través de un estimador de densidad de núcleo, que creo que ya tengo.

Como mi objetivo es obtener la probabilidad de este nuevo objeto ( $p[x_p, y_p]$ ) de pertenecer a este conjunto de datos 2D $\Omega$ Me dijeron que integrara el pdf $\hat{f}$ sobre " valores del soporte para los que la densidad es menor que la observada ". La densidad "observada" es $\hat{f}$ evaluado en el nuevo objeto $p$ , es decir: $\hat{f}(x_p, y_p)$ . Así que tengo que resolver la ecuación:

$$\iint_{x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)} \hat{f}(x,y)\,dx\,dy$$

El PDF de mi conjunto de datos 2D (obtenido a través de la función de python stats.gaussian_kde ) tiene el siguiente aspecto:

donde el punto rojo representa el nuevo objeto $p[x_p, y_p]$ trazado sobre el PDF de mi conjunto de datos.

Entonces la pregunta es: ¿cómo puedo calcular la integral anterior para los límites $x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$ cuando el pdf se ve así?

---

Añadir

Hice algunas pruebas para ver qué tal funcionaba el método de Montecarlo que menciono en uno de los comentarios. Esto es lo que obtuve:

Los valores parecen variar un poco más para las zonas de menor densidad, y ambos anchos de banda muestran más o menos la misma variación. La mayor variación de la tabla se produce en el punto (x,y)=(2,4,1,5), comparando el valor de muestra de Silverman de 2.500 frente al de 1.000, lo que da una diferencia de 0.0126 o ~1.3% . En mi caso, esto sería ampliamente aceptable.

Editar : Acabo de comprobar que en 2 dimensiones la regla de Scott es equivalente a la de Silverman según la definición dada aquí .

Answer 1

Una forma sencilla es rasterizar el dominio de integración y calcular una aproximación discreta a la integral.

Hay que tener en cuenta algunas cosas:

1. Asegúrese de cubrir más que la extensión de los puntos: hay que incluir todos los lugares en los que la estimación de la densidad del núcleo tendrá algún valor apreciable. Esto significa que debe ampliar la extensión de los puntos en tres o cuatro veces el ancho de banda del kernel (para un kernel gaussiano).

2. El resultado variará un poco con la resolución de la trama. La resolución debe ser una pequeña fracción del ancho de banda. Dado que el tiempo de cálculo es proporcional al número de celdas del ráster, casi no se necesita tiempo adicional para realizar una serie de cálculos utilizando resoluciones más gruesas que la prevista: compruebe que los resultados de las más gruesas convergen en el resultado de la resolución más fina. Si no es así, es posible que se necesite una resolución más fina.

A continuación se muestra una ilustración para un conjunto de datos de 256 puntos:

Los puntos se muestran como puntos negros superpuestos en dos estimaciones de densidad del núcleo. Los seis puntos rojos grandes son "sondas" en las que se evalúa el algoritmo. Esto se ha hecho para cuatro anchos de banda (un valor por defecto entre 1,8 (verticalmente) y 3 (horizontalmente), 1/2, 1 y 5 unidades) con una resolución de 1000 por 1000 celdas. La siguiente matriz de dispersión muestra la fuerte dependencia de los resultados del ancho de banda para estos seis puntos de sondeo, que cubren una amplia gama de densidades:

La variación se produce por dos razones. Obviamente, las estimaciones de densidad difieren, lo que introduce una forma de variación. Y lo que es más importante, las diferencias en las estimaciones de densidad pueden crear grande diferencias en cualquier punto ("sonda"). Esta última variación es mayor en torno a las "franjas" de densidad media de los grupos de puntos, precisamente los lugares en los que es más probable que se utilice este cálculo.

Esto demuestra la necesidad de tener mucha precaución al utilizar e interpretar los resultados de estos cálculos, ya que pueden ser muy sensibles a una decisión relativamente arbitraria (el ancho de banda a utilizar).

---

Código R

El algoritmo está contenido en la media docena de líneas de la primera función, f . Para ilustrar su uso, el resto del código genera las figuras anteriores.

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)

Answer 2

Si tiene un número decente de observaciones, es posible que no necesite hacer ninguna integración. Digamos que su nuevo punto es ${\bf{x}}_0$ . Suponga que tiene un estimador de densidad $\hat f$ ; sumar el número de observaciones ${\bf{x}}$ para lo cual $\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$ y dividir por el tamaño de la muestra. Esto le da una aproximación a la probabilidad requerida.

Esto supone que $\hat f({\bf{x}}_0)$ no es "demasiado pequeño" y el tamaño de la muestra es lo suficientemente grande (y está lo suficientemente repartido) como para ofrecer una estimación decente en las regiones de baja densidad. Sin embargo, 20000 casos parecen ser lo suficientemente grandes, para la bivariante ${\bf{x}}$ .