Integrando o estimador de densidade do kernel em 2D

Estou saindo dessa pergunta , caso alguém queira seguir a trilha.

Basicamente, eu tenho um conjunto de dados $\Omega$ composto por $N$ objetos em que cada objeto tem um determinado número de valores medidos anexados (dois neste caso):

Eu preciso de uma maneira de determinar a probabilidade de um novo objeto $p[x_p, y_p]$ de pertencer a $\Omega$ então eu foi aconselhado nessa questão para obter uma densidade de probabilidade f através de um estimador de densidade kernel, que eu acredito que eu já tenho .$\hat{f}$

Desde que meu objetivo é obter a probabilidade de este novo objeto ( ) de pertencer a este conjunto de dados 2D Ω , disseram-me para integrar o pdf f sobre " valores do apoio para o qual a densidade é menor que o que você observou ". A densidade "observados" é f avaliada no novo objecto p , isto é: f ( x p , y p ) . Então, eu preciso resolver a equação:$p[x_p, y_p]$$\Omega$$\hat{f}$$\hat{f}$$p$$\hat{f}(x_p, y_p)$

O PDF do meu conjunto de dados 2D (obtido através do módulo stats.gaussian_kde do python ) fica assim:

onde o ponto vermelho representa o novo objeto plotado sobre o PDF do meu conjunto de dados.$p[x_p, y_p]$

Então a questão é: como eu posso calcular a integral acima para os limites quando o PDF se parece com isso?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

Adicionar

Fiz alguns testes para ver como o método de Monte Carlo mencionado em um dos comentários funcionou. Isto é o que eu tenho:

Os valores parecem variar um pouco mais para áreas de menor densidade, com ambas as larguras de banda mostrando mais ou menos a mesma variação. A maior variação na tabela ocorre para o ponto (x, y) = (2.4,1.5) comparando o valor da amostra 2500 vs 1000 do Silverman, que fornece uma diferença de 0.0126ou ~1.3%. No meu caso, isso seria amplamente aceitável.

Edit : Acabei de notar que, em 2 dimensões, a regra de Scott é equivalente à de Silverman, de acordo com a definição dada aqui .

Uma maneira simples é rasterizar o domínio de integração e calcular uma aproximação discreta à integral.

Há algumas coisas a serem observadas:

1. Certifique-se de cobrir mais do que a extensão dos pontos: você precisa incluir todos os locais onde a estimativa de densidade do kernel terá valores apreciáveis. Isso significa que você precisa expandir a extensão dos pontos em três a quatro vezes a largura de banda do kernel (para um kernel gaussiano).

2. O resultado varia um pouco com a resolução da varredura. A resolução precisa ser uma pequena fração da largura de banda. Como o tempo de cálculo é proporcional ao número de células na varredura, não leva quase tempo extra para executar uma série de cálculos usando resoluções mais grossas do que a pretendida: verifique se os resultados para os mais grossos estão convergindo no resultado para o melhor resolução. Caso contrário, uma resolução mais precisa pode ser necessária.

Aqui está uma ilustração para um conjunto de dados de 256 pontos:

Os pontos são mostrados como pontos pretos sobrepostos em duas estimativas de densidade do kernel. Os seis grandes pontos vermelhos são "sondas" nas quais o algoritmo é avaliado. Isso foi feito para quatro larguras de banda (um padrão entre 1,8 (verticalmente) e 3 (horizontalmente), 1/2, 1 e 5 unidades) com uma resolução de 1000 por 1000 células. A seguinte matriz de gráfico de dispersão mostra quão fortemente os resultados dependem da largura de banda para esses seis pontos de sonda, que cobrem uma ampla variedade de densidades:

A variação ocorre por dois motivos. Obviamente, as estimativas de densidade diferem, introduzindo uma forma de variação. Mais importante, as diferenças nas estimativas de densidade podem criar grandes diferenças em qualquer ponto único ("sonda"). A última variação é maior em torno das "franjas" de densidade média de aglomerados de pontos - exatamente nos locais em que esse cálculo provavelmente será mais utilizado.

Isso demonstra a necessidade de cautela substancial no uso e na interpretação dos resultados desses cálculos, porque eles podem ser muito sensíveis a uma decisão relativamente arbitrária (a largura de banda a ser usada).

Código R

O algoritmo está contido na meia dúzia de linhas da primeira função f,. Para ilustrar seu uso, o restante do código gera as 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)

Se você tiver um número decente de observações, talvez não seja necessário fazer nenhuma integração. Diga que seu novo ponto é${\bf{x}}_0$. Suponha que você tenha um estimador de densidade$\hat f$; resumir o número de observações${\bf{x}}$ para qual $\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$e divida pelo tamanho da amostra. Isso fornece uma aproximação à probabilidade necessária.

Isso pressupõe que $\hat f({\bf{x}}_0)$não é "muito pequeno" e o tamanho da amostra é grande o suficiente (e espalhado o suficiente) para fornecer uma estimativa decente nas regiões de baixa densidade. 20000 casos parecem grandes o suficiente, para bivariados${\bf{x}}$.