Integrando el estimador de densidad del kernel en 2D

Vengo de esta pregunta en caso de que alguien quiera seguir el rastro.

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

Ω = o_{1} [x_{1}, y_{1}], o_{2} [x_{2}, y_{2}], . . ., o_{N} [x_{N}, y_{N}]

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

Necesito una manera de determinar la probabilidad de un nuevo objeto $p[x_p, y_p]$ de pertenencia a $\Omega$ por lo que estaba previsto en esa pregunta para obtener una densidad de probabilidad a través de un estimador de densidad de kernel, que creo que ya tengo . $\hat{f}$

Desde mi objetivo es obtener la probabilidad de que este nuevo objeto ( $p[x_p, y_p]$ ) de que pertenecen a este conjunto de datos en 2D $\Omega$ , me dijeron para integrar el pdf sobre " valores de la ayuda para el que la densidad es menor que la que observaste ". La densidad "observado" es evaluado en el nuevo objeto de , es . Entonces necesito resolver la ecuación: $\hat{f}$ $\hat{f}$ $p$ $\hat{f}(x_p, y_p)$

\iint_{x, y : \hat{f} (x, y) < \hat{f} (x_{p}, y_{p})} \hat{f} (x, y) d x d y

$\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 del módulo stats.gaussian_kde de Python ) se ve así:

ingrese la descripción de la imagen aquí

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

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

Añadir

Hice algunas pruebas para ver qué tan bien funcionaba el método de Monte Carlo que menciono en uno de los comentarios. Esto es lo que conseguí:

mesa

Los valores parecen variar un poco más para las áreas de menor densidad con ambos anchos de banda que muestran más o menos la misma variación. La mayor variación en la tabla se produce para el punto (x, y) = (2.4,1.5) que compara el valor de muestra de Silverman 2500 vs 1000, lo que da una diferencia de 0.0126o ~1.3%. En mi caso, esto sería en gran medida aceptable.

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

probability maximum-likelihood kernel-smoothing numerical-integration Gabriel
fuente

¿Ha notado que su estimador no es unimodal, pero que la recomendación que está siguiendo se aplica explícitamente solo a distribuciones "unimodales"? Eso no significa que esté haciendo algo mal, pero debería generar un pensamiento profundo sobre lo que podría significar la respuesta.

whuber

Hola @whuber, en realidad la respuesta en esa pregunta dice que se "comportó bien" para las distribuciones unimodales, así que pensé que tal vez podría funcionar en mi problema con alguna modificación. ¿"Comportarse bien" significa "solo funciona" en la jerga estadística (pregunta honesta)? Salud.

Gabriel

Mi principal preocupación es que el KDE puede ser sensible a la elección del ancho de banda y espero que su integral, especialmente para lugares marginales como el que se muestra en la ilustración, sea muy sensible a la elección. (El cálculo en sí mismo, por cierto, es fácil una vez que ha creado una imagen ráster como esta: es proporcional al valor promedio en la imagen entre puntos cuyo valor es menor que el del punto de "sonda".) Puede acercarse esto al calcular la respuesta para un rango completo de anchos de banda razonables y ver si cambia de alguna manera material dentro de ese rango. Si no, estás bien.

whuber

No comentaré sobre la solución, pero la integración se puede hacer con Monte Carlo simple: puntos de muestra de (eso es fácil, ya que el kde es una mezcla de densidades de las que es fácil tomar muestras), y contar la fracción de puntos que caen dentro de la región de integración (donde se mantiene la desigualdad).

\hat{f}

$\hat{f}$

Zen

¿Cuántas observaciones tiene en su conjunto de datos?

Hong Ooi

Respuestas:

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

Hay algunas cosas a tener en cuenta:

Asegúrese de cubrir más que la extensión de los puntos: debe incluir todas las ubicaciones donde la estimación de densidad del núcleo tendrá valores apreciables. Esto significa que necesita expandir la extensión de los puntos entre tres y cuatro veces el ancho de banda del núcleo (para un núcleo gaussiano).
El resultado variará un poco con la resolución del ráster. La resolución debe ser una pequeña fracción del ancho de banda. Debido a que el tiempo de cálculo es proporcional al número de celdas en el ráster, casi no se necesita tiempo adicional para realizar una serie de cálculos con resoluciones más gruesas que la prevista: verifique que los resultados para los más gruesos converjan en el resultado para La mejor resolución. Si no lo son, puede ser necesaria una resolución más fina.

Aquí hay una ilustración para un conjunto de datos de 256 puntos:

Figura 1

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 predeterminado 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 diagrama de dispersión muestra cuán fuertemente dependen los resultados del ancho de banda para estos seis puntos de sonda, que cubren una amplia gama de densidades:

Figura 2

La variación ocurre por dos razones. Obviamente, las estimaciones de densidad difieren, introduciendo una forma de variación. Más importante aún, las diferencias en las estimaciones de densidad pueden crear grandes diferencias en cualquier punto único ("sonda"). La última variación es mayor en torno a las "franjas" de densidad media de los grupos de puntos, exactamente aquellos lugares donde es probable que este cálculo sea el más utilizado.

Esto demuestra la necesidad de una precaución sustancial al usar 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 usar).

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)

whuber
fuente

Respuesta increíble, aunque no estoy seguro de entender el significado de las Defaulty Hp5los anchos de banda (supongo H1y H5media h=1y h=5) es Hp5el valor h=1/2? Si es así, ¿qué es Default?

Gabriel

Su comprensión es correcta. ("p5" significa ".5".) El valor predeterminado se calcula automáticamente mediante el kde2duso bandwidth.nrd. Para los datos de la muestra, es igual a en la dirección horizontal y en la dirección vertical, aproximadamente a la mitad entre los valores de y en la prueba. Tenga en cuenta que estos anchos de banda predeterminados son lo suficientemente grandes como para colocar una proporción apreciable de la densidad total mucho más allá de la extensión de los puntos en sí, por lo que esa extensión debe ampliarse independientemente del algoritmo de integración que elija utilizar.

3

$3$

1.85

$1.85$

1

$1$

5

$5$

whuber

Entonces, ¿entiendo correctamente si digo que a medida que aumente el ancho de banda utilizado, la extensión del resultado kdetambién aumenta (y por lo tanto necesito extender los límites de integración)? Dado que puedo vivir con un error <10%en el valor resultante de la integral, ¿qué piensa sobre el uso de la regla de Scott?

Gabriel

Creo que debido a que estas reglas se desarrollaron para objetivos completamente diferentes, debe sospechar que podrían no funcionar bien para sus propósitos, especialmente si se trata de implementar una sugerencia hecha en stats.stackexchange.com/questions/63263 . Es prematuro preocuparse por la regla general que podría usar para el KDE; en esta etapa, debe preocuparse seriamente si el enfoque completo incluso funcionará de manera confiable.

whuber

Rasca lo anterior. Me hacer tener una manera de saber si la aplicación está funcionando e incluso de la cuantificación de lo bien que está funcionando. Es un poco complicado y requiere mucho tiempo, pero puedo (debería poder) hacerlo.

Gabriel

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

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

Hong Ooi
fuente

Sería bienvenido algún análisis cuantitativo de esta recomendación, o al menos un ejemplo de una aplicación real. Sospecho que la precisión de su propuesta depende en gran medida de la forma del núcleo. Esto me haría reacio a confiar en tal cálculo sin un estudio sustancial de sus propiedades.

whuber