Divergência de Kullback-Leibler para duas amostras

Tentei implementar uma estimativa numérica da divergência Kullback-Leibler para duas amostras. Para depurar a implementação, extrair as amostras de duas distribuições normais e .$\mathcal N (0,1)$$\mathcal N (1,2)$

Para uma estimativa simples, gerei dois histogramas e tentei aproximar numericamente a integral. Fiquei preso ao lidar com aquelas partes do histograma em que os compartimentos de um dos histogramas são zero, de modo que acabo dividindo por zero ou com o logaritmo de zero. Como eu manejo esse problema?

Uma pergunta relacionada surgiu em minha mente: como calcular exatamente a divergência KL entre duas distribuições uniformes diferentes? Preciso restringir a integral à união do suporte de ambas as distribuições?

A divergência Kullback-Leibler é definida como para calcular (estimar) isso a partir de dados empíricos, talvez precisássemos de algumas estimativas das funções de densidade . Portanto, um ponto de partida natural pode ser via estimativa de densidade (e depois disso, apenas integração numérica). Quão bom ou estável seria esse método, não sei.

Mas primeiro sua segunda pergunta, depois voltarei à primeira. Digamos que e sejam densidades uniformes em e respectivamente. Então enquanto é mais difícil de definir, mas o único valor razoável a ser dado é , tanto quanto posso ver, pois envolve integrando que podemos escolher para interpretar como . Estes resultados são razoáveis ​​a partir da interpretação que dou em Intuição sobre a divergência de Kullback-Leibler (KL)$p$$q$$[0,1]$$[0,10]$$\KL(p || q)=\log 10$$\KL(q || p)$$\infty$$\log(1/0)$$\log \infty$

Voltando à questão principal. É solicitado de uma maneira muito não paramétrica, e nenhuma suposição é declarada sobre as densidades. Provavelmente são necessárias algumas suposições. Mas, supondo que as duas densidades sejam propostas como modelos concorrentes para o mesmo fenômeno, provavelmente podemos assumir que elas tenham a mesma medida dominante: a divergência de KL entre uma distribuição de probabilidade contínua e uma discreta será sempre infinito, por exemplo. Um artigo abordando esta questão é o seguinte: https://pdfs.semanticscholar.org/1fbd/31b690e078ce938f73f14462fceadc2748bf.pdf Eles propõem um método que não precisa de estimativa preliminar de densidade e analisa suas propriedades.

(Existem muitos outros documentos). Voltarei e publicarei alguns detalhes desse trabalho, as idéias.

 EDIT               

Algumas idéias desse artigo, que trata da estimativa da divergência de KL com amostras de iid de distribuições absolutamente contínuas. Mostro sua proposta para distribuições unidimensionais, mas elas também fornecem uma solução para vetores (usando a estimativa de densidade de vizinhos mais próximos). Para provas, leia o artigo!

Eles propõem o uso de uma versão da função de distribuição empírica, mas interpolam linearmente entre os pontos da amostra para obter uma versão contínua. Eles definem onde é a função da etapa Heavyside, mas definida para que . Então essa função interpolada linearmente (e estendida horizontalmente além do intervalo) é ( para contínua). Então eles propõem estimar a divergência de Kullback-Leibler por onde e

O código R para a versão da função de distribuição empírica de que precisamos é

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
}          

note que rleé usado para cuidar do caso com duplicatas em x.

Então a estimativa da divergência KL é dada por

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              
}

Então eu mostro uma pequena simulação:

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

que fornece o seguinte histograma, mostrando (uma estimativa) a distribuição amostral desse estimador:

Para comparação, calculamos a divergência de KL neste exemplo por integração numérica:

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 ... a diferença é grande o suficiente para que haja muito aqui para investigar!

Expandindo um pouco a resposta de kjetil-b-halvorsen , e desculpe-me por não comentar, não tenho reputação:

1. Tenho a sensação de que a computação analítica deve ser (sem multiplicação por 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. Se eu estiver certo, o estimador não converge para a divergência KL, mas a convergência é declarada como: . A seta representa como convergência.$\hat D(P||Q)$$\hat D(P||Q)-1 \to D(P||Q)$

Depois que essas duas correções são feitas, os resultados parecem mais realistas.