Intégration de l'estimateur de densité de noyau en 2D

Je viens de cette question au cas où quelqu'un voudrait suivre la piste.

Fondamentalement, j'ai un ensemble de données composé de objets où chaque objet a un nombre donné de valeurs mesurées qui lui sont attachées (deux dans ce cas):$\Omega$$N$

J'ai besoin d'un moyen de déterminer la probabilité qu'un nouvel objet appartienne à donc on m'a conseillé dans cette question d'obtenir une densité de probabilité grâce à un estimateur de densité de noyau, que je crois avoir J'ai déjà.Ω $p[x_p, y_p]$$\Omega$$\hat{f}$

Comme mon objectif est d'obtenir la probabilité d' appartenance de ce nouvel objet ( ) à cet ensemble de données 2D , on m'a dit d'intégrer le pdf sur "les valeurs du support pour lesquelles la densité est inférieure à celle que vous avez observée ". La densité "observée" est évaluée dans le nouvel objet , c'est-à-dire: . J'ai donc besoin de résoudre l'équation:Ω $p[x_p, y_p]$$\Omega$$\hat{f}$ p f (xp,yp$\hat{f}$$p$$\hat{f}(x_p, y_p)$

Le PDF de mon ensemble de données 2D (obtenu via le module stats.gaussian_kde de python ) ressemble à ceci:

où le point rouge représente le nouvel objet tracé sur le PDF de mon ensemble de données.$p[x_p, y_p]$

La question est donc: comment puis-je calculer l'intégrale ci-dessus pour les limites lorsque le pdf ressemble à ça?$x, y:\hat{f}(x, y) < \hat{f}(x_p, y_p)$

Ajouter

J'ai fait quelques tests pour voir si la méthode Monte Carlo que je mentionne dans l'un des commentaires fonctionnait bien. Voici ce que j'ai obtenu:

Les valeurs semblent varier un peu plus pour les zones de faible densité, les deux bandes passantes montrant plus ou moins la même variation. La plus grande variation dans le tableau se produit pour le point (x, y) = (2,4,1,5) comparant la valeur d'échantillon de Silverman 2500 vs 1000, ce qui donne une différence de 0.0126ou ~1.3%. Dans mon cas, cela serait largement acceptable.

Edit : je viens de remarquer qu'en 2 dimensions, la règle de Scott est équivalente à celle de Silverman selon la définition donnée ici .

Un moyen simple consiste à pixelliser le domaine de l'intégration et à calculer une approximation discrète de l'intégrale.

Il y a certaines choses à surveiller:

1. Assurez-vous de couvrir plus que l'étendue des points: vous devez inclure tous les emplacements où l'estimation de la densité du noyau aura des valeurs appréciables. Cela signifie que vous devez étendre l'étendue des points de trois à quatre fois la bande passante du noyau (pour un noyau gaussien).

2. Le résultat variera quelque peu avec la résolution du raster. La résolution doit être une petite fraction de la bande passante. Comme le temps de calcul est proportionnel au nombre de cellules dans le raster, il ne prend presque pas de temps supplémentaire pour effectuer une série de calculs en utilisant des résolutions plus grossières que celle prévue: vérifiez que les résultats pour les plus grossiers convergent sur le résultat pour le meilleure résolution. Si ce n'est pas le cas, une résolution plus fine peut être nécessaire.

Voici une illustration pour un ensemble de données de 256 points:

Les points sont représentés par des points noirs superposés à deux estimations de densité de noyau. Les six grands points rouges sont des "sondes" auxquelles l'algorithme est évalué. Cela a été fait pour quatre bandes passantes (une valeur par défaut entre 1,8 (verticalement) et 3 (horizontalement), 1/2, 1 et 5 unités) à une résolution de 1000 par 1000 cellules. La matrice de diagramme de dispersion suivante montre à quel point les résultats dépendent de la bande passante pour ces six points de sonde, qui couvrent une large gamme de densités:

La variation se produit pour deux raisons. De toute évidence, les estimations de densité diffèrent, introduisant une forme de variation. Plus important encore, les différences dans les estimations de densité peuvent créer de grandes différences en tout point ("sonde"). Cette dernière variation est plus importante autour des «franges» de densité moyenne de groupes de points - exactement les endroits où ce calcul est susceptible d'être le plus utilisé.

Cela démontre la nécessité d'une grande prudence dans l'utilisation et l'interprétation des résultats de ces calculs, car ils peuvent être si sensibles à une décision relativement arbitraire (la bande passante à utiliser).

Code R

L'algorithme est contenu dans la demi-douzaine de lignes de la première fonction f,. Pour illustrer son utilisation, le reste du code génère les figures précédentes.

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)

Si vous avez un nombre décent d'observations, vous n'aurez peut-être pas besoin de faire d'intégration du tout. Supposons que votre nouveau point soit . Supposons que vous ayez un estimateur de densité ; résumer le nombre d'observations pour lesquelles et diviser par la taille de l'échantillon. Cela vous donne une approximation de la probabilité requise.${\bf{x}}_0$$\hat f$f ( x )< f ( x 0${\bf{x}}$$\hat f({\bf{x}}) < \hat f({\bf{x}}_0)$

Cela suppose que n'est pas "trop ​​petit" et que la taille de votre échantillon est suffisamment grande (et suffisamment étalée) pour donner une estimation décente dans les régions à faible densité. Cependant, 20000 cas semblent assez volumineux, pour les .$\hat f({\bf{x}}_0)$${\bf{x}}$