Divergência de Kullback-Leibler para duas amostras

10

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 .N(0,1)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?

Jimbob
fonte
Bem, o suporte da distribuição normal é o conjunto de números reais. Não há problema em matemática pura, mas sim, para sua aproximação numérica, é necessário garantir que o tamanho da amostra seja grande o suficiente em relação à região na qual você deseja integrar. Você não será capaz de integrar sobre (-inf, + inf) como é possível em matemática pura ... Escolha algo razoável? Se você é mais de 3 desvios padrão de distância da média, que vai ser muito fina ...
Matthew Gunn
11
Com relação à sua segunda pergunta, a divergência KL entre duas distribuições uniformes diferentes é indefinida ( é indefinido). Da mesma forma, a divergência KL para duas distribuições empíricas é indefinida, a menos que cada amostra tenha pelo menos uma observação com o mesmo valor que todas as observações na outra amostra. log(0)
jbowman
@jbowman Nota pequena. Embora você esteja certo de que é indefinido (ou ), é habitual na teoria da informação tratar como . log(0)log(0)00
Luca Citi
Uma pergunta semelhante: mathoverflow.net/questions/119752/…
kjetil b halvorsen

Respostas:

9

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.

KL(P||Q)=p(x)logp(x)q(x)dx
p(x),q(x)

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)pq[0,1][0,10]KL(p||q)=log10KL(q||p)log(1/0)log

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

Pe(x)=1ni=1nU(xxi)
UU(0)=0.5Pcc
D^(PQ)=1ni=1nlog(δPc(xi)δQc(xi))
δPc=Pc(xi)Pc(xiϵ)ϵ é um número menor que o menor espaçamento das amostras.

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:

Distribuição amostral do estimador KL

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!

kjetil b halvorsen
fonte
5

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

  1. 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.D^(P||Q)D^(P||Q)1D(P||Q)

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

ColibriIO
fonte
Obrigado, vou analisar isso e atualizar minha resposta.
Kjetil b halvorsen