Integrando o estimador de densidade do kernel em 2D

12

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

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

Ω=o1[x1,y1],o2[x2,y2],...,oN[xN,yN]

Eu preciso de uma maneira de determinar a probabilidade de um novo objeto p[xp,yp] de pertencer a Ω 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 .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[xp,yp]Ωf^f^pf^(xp,yp)

x,y:f^(x,y)<f^(xp,yp)f^(x,y)dxdy

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

insira a descrição da imagem aqui

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

Então a questão é: como eu posso calcular a integral acima para os limites quando o PDF se parece com isso?x,y:f^(x,y)<f^(xp,yp)


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:

mesa

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 .

Gabriel
fonte
2
Você percebeu que seu estimador não é unimodal, mas que a recomendação que você está seguindo se aplica explicitamente apenas a distribuições "unimodais"? Isso não significa que você esteja fazendo algo errado, mas deve gerar algumas reflexões sobre o que a resposta pode significar.
whuber
Olá @whuber, na verdade a resposta nessa pergunta diz que é "bem comportada" para distribuições unimodais, então pensei que talvez pudesse funcionar no meu problema com algumas modificações. "Bem comportado" significa "só funciona" no jargão estatístico (pergunta honesta)? Felicidades.
Gabriel
Minha principal preocupação é que o KDE possa ser sensível à escolha da largura de banda e espero que sua integral, especialmente em locais marginais como o mostrado na ilustração, seja muito sensível à escolha. (O cálculo propriamente dito, a propósito, é fácil depois que você cria uma imagem rasterizada como esta: é proporcional ao valor médio na imagem entre pontos cujo valor é menor que o do ponto "sonda".) isso calculando a resposta para uma gama completa de larguras de banda razoáveis ​​e ver se isso muda de alguma maneira material dentro desse intervalo. Se não, você está bem.
whuber
Não vou comentar sobre a solução, mas a integração pode ser feito por simples Monte Carlo: pontos de amostragem de ff^ é fácil, uma vez que o KDE é uma mistura de densidades que são fáceis de amostra), e contar a fração de pontos que se enquadram na região de integração (onde a desigualdade se mantém).
Zen
Quantas observações você tem no seu conjunto de dados?
Hong Ooi 27/07/2013

Respostas:

11

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:

figura 1

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:

Figura 2

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)
whuber
fonte
Resposta surpreendente, embora eu não tenha certeza de que entendi o significado das larguras de banda Defaulte Hp5(presumo H1e H5quero dizer h=1e h=5). É Hp5o valor h=1/2? Se sim, o que é Default?
Gabriel
1
Sua compreensão está correta. ("p5" significa ".5".) O padrão é calculado automaticamente kde2dusando bandwidth.nrd. Para os dados de amostra, é igual a3 na direção horizontal e 1,85 na direção vertical, aproximadamente a meio caminho entre os valores de 1 e 5no teste. Observe que essas larguras de banda padrão são grandes o suficiente para colocar uma proporção apreciável da densidade total muito além da extensão dos pontos em si, e é por isso que essa extensão precisa ser expandida independentemente do algoritmo de integração que você escolher usar.
whuber
Então, estou entendendo corretamente se digo que à medida que aumenta a largura de banda usada, a extensão do resultado kdetambém aumenta (e, portanto, preciso estender os limites de integração)? Dado que eu posso viver com um erro <10%no valor resultante da integral, o que você acha de usar a regra de Scott?
Gabriel
Acho que, como essas regras foram desenvolvidas para objetivos totalmente diferentes, você deve suspeitar que elas podem não ter um bom desempenho para seus propósitos, especialmente se for para implementar uma sugestão feita em stats.stackexchange.com/questions/63263 . É prematuro se preocupar com a regra de ouro que você pode usar para o KDE; Nesse estágio, você deve estar seriamente preocupado se a abordagem inteira funcionará de maneira confiável.
whuber
1
Risque o acima. Eu faço tem uma maneira de saber se a implementação está a trabalhar e até mesmo de quantificar o quão bem ele está trabalhando. É um pouco complicado e demorado, mas eu posso (deveria) fazer isso.
Gabriel
1

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 éx0 0. Suponha que você tenha um estimador de densidadef^; resumir o número de observaçõesx para qual f^(x)<f^(x0 0)e divida pelo tamanho da amostra. Isso fornece uma aproximação à probabilidade necessária.

Isso pressupõe que f^(x0 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 bivariadosx.

Hong Ooi
fonte
Alguma análise quantitativa desta recomendação, ou pelo menos um exemplo de aplicação real, seria bem-vinda. Suspeito que a precisão da sua proposta dependa fortemente da forma do kernel. Isso me deixaria relutante em confiar em tal cálculo sem um estudo substancial de suas propriedades.
whuber