Kullback-Leibler Divergence pour deux échantillons

J'ai essayé d'implémenter une estimation numérique de la divergence de Kullback-Leibler pour deux échantillons. Pour déboguer l'implémentation, tirez les échantillons de deux distributions normales et .N ( 1 , 2 $\mathcal N (0,1)$$\mathcal N (1,2)$

Pour une estimation simple, j'ai généré deux histogrammes et essayé d'approximer numériquement l'intégrale. Je suis resté coincé avec la manipulation des parties de l'histogramme où les cases de l'un des histogrammes sont nulles, de sorte que je me retrouve avec une division par zéro ou le logarithme de zéro. Comment gérer ce problème?

Une question connexe m'est venue à l'esprit: comment calculer exactement la KL-Divergence entre deux distributions uniformes différentes? Dois-je restreindre l'intégrale à l'union du support des deux distributions?

La divergence Kullback-Leibler est définie comme donc pour calculer (estimer) ceci à partir de données empiriques, nous aurions besoin, peut-être, de quelques estimations des fonctions de densité . Ainsi, un point de départ naturel pourrait être via l'estimation de la densité (et après cela, juste l'intégration numérique). Je ne sais pas à quel point une telle méthode serait bonne ou stable.p ( x ) , q ( x

Mais d'abord votre deuxième question, puis je reviendrai à la première. Disons que et sont des densités uniformes sur et respectivement. Alors tandis que est plus difficile à définir, mais la seule valeur raisonnable à lui donner est , pour autant que je puisse voir, car elle implique intégrant que nous pouvons choisir d'interpréter comme . Ces résultats sont raisonnables d'après l'interprétation que je donne dans Intuition on the Kullback-Leibler (KL) Divergenceq [ 0 , 1 ] [ 0 , 10 ] KL ( p | | q ) = log 10 KL ( q | | p ) ∞ log ( 1 / 0 ) log $p$$q$$[0,1]$$[0,10]$$\KL(p || q)=\log 10$$\KL(q || p)$$\infty$$\log(1/0)$$\log \infty$

Revenons à la question principale. Elle est posée de manière très non paramétrique, et aucune hypothèse n'est formulée sur les densités. Certaines hypothèses sont probablement nécessaires. Mais en supposant que les deux densités sont proposées comme modèles concurrents pour le même phénomène, nous pouvons probablement supposer qu'elles ont la même mesure dominante: la divergence KL entre une distribution de probabilité continue et discrète serait toujours l'infini, par exemple. Un article traitant de cette question est le suivant: https://pdfs.semanticscholar.org/1fbd/31b690e078ce938f73f14462fceadc2748bf.pdf Ils proposent une méthode qui ne nécessite pas d'estimation préalable de la densité et analyse ses propriétés.

(Il existe de nombreux autres articles). Je reviendrai et publierai quelques détails de ce document, les idées.

 EDIT               

Quelques idées de cet article, qui concerne l'estimation de la divergence KL avec des échantillons iid à partir de distributions absolument continues. Je montre leur proposition de distributions unidimensionnelles, mais ils donnent également une solution pour les vecteurs (en utilisant l'estimation de la densité du plus proche voisin). Pour les épreuves, lisez le papier!

Ils proposent d'utiliser une version de la fonction de distribution empirique, mais interpolée linéairement entre les points d'échantillonnage pour obtenir une version continue. Ils définissent

Le code R pour la version de la fonction de distribution empirique dont nous avons besoin est

my.ecdf  <-  function(x)   {
    x   <-   sort(x)
    x.u <-   unique(x)
    n  <-  length(x) 
    x.rle  <-  rle(x)$lengths
    y  <-  (cumsum(x.rle)-0.5) / n
    FUN  <-  approxfun(x.u, y, method="linear", yleft=0, yright=1,
                           rule=2)
    FUN
}          

notez qu'il rleest utilisé pour prendre soin de l'affaire avec des doublons dans x.

L'estimation de la divergence KL est alors donnée par

KL_est  <-  function(x, y)   {
    dx  <-  diff(sort(unique(x)))
    dy  <-  diff(sort(unique(y)))
    ex  <-  min(dx) ; ey  <-  min(dy)
    e   <-  min(ex, ey)/2
    n   <-  length(x)    
    P  <-   my.ecdf(x) ; Q  <-  my.ecdf(y)
    KL  <-  sum( log( (P(x)-P(x-e))/(Q(x)-Q(x-e)))) / n
    KL              
}

Ensuite, je montre une petite simulation:

KL  <-  replicate(1000, {x  <-  rnorm(100)
                         y <- rt(100, df=5)
                         KL_est(x, y)})
hist(KL, prob=TRUE)

qui donne l'histogramme suivant, montrant (une estimation) de la distribution d'échantillonnage de cet estimateur:

Pour comparaison, nous calculons la divergence KL dans cet exemple par intégration numérique:

LR  <-  function(x) dnorm(x,log=TRUE)-dt(x,5,log=TRUE)
100*integrate(function(x) dnorm(x)*LR(x),lower=-Inf,upper=Inf)$value
[1] 3.337668

hmm ... la différence étant assez grande pour qu'il y ait beaucoup à étudier!

Développant un peu la réponse de kjetil-b-halvorsen , et désolé de ne pas avoir commenté, je n'ai pas la réputation:

1. J'ai le sentiment que le calcul analytique devrait être (sans multiplication par 100):

LR <- function(x) dnorm(x,log=TRUE)-dt(x,5,log=TRUE) integrate(function(x) dnorm(x)*LR(x),lower=-Inf,upper=Inf)$value

2. Si j'ai raison, l'estimateur ne converge pas vers la divergence KL, mais la convergence est indiquée comme suit: . La flèche représente la convergence. D (P||Q)-1→D(P||Q$\hat D(P||Q)$$\hat D(P||Q)-1 \to D(P||Q)$

Une fois ces deux corrections apportées, les résultats semblent plus réalistes.