Rozbieżność Kullbacka-Leiblera dla dwóch próbek

10

Próbowałem zaimplementować oszacowanie liczbowe dywergencji Kullbacka-Leiblera dla dwóch próbek. Aby debugować implementację, narysuj próbki z dwóch rozkładów normalnych N(0,1) i N(1,2) .

Dla prostego oszacowania wygenerowałem dwa histogramy i próbowałem liczbowo aproksymować całkę. Utknąłem z obsługą tych części histogramu, w których przedziały jednego z histogramów mają zero, tak że albo skończę z dzieleniem przez zero, albo logarytmem zera. Jak poradzić sobie z tym problemem?

Przyszło mi do głowy podobne pytanie: jak dokładnie obliczyć dywergencję KL między dwoma różnymi rozkładami jednolitych? Czy muszę ograniczyć całkę do unii wsparcia obu dystrybucji?

Jimbob
źródło
Cóż, wsparcie rozkładu normalnego jest zbiorem liczb rzeczywistych. W czystej matematyce nie ma problemu, ale tak, dla twojego przybliżenia liczbowego musisz upewnić się, że wielkość próbki jest wystarczająco duża w stosunku do regionu, który chcesz zintegrować. Nie będziesz w stanie zintegrować się z (-inf, + inf) tak jak w czystej matematyce ... Poszukać czegoś rozsądnego? Jeśli dzieli Cię więcej niż 3 standardowe odchylenia od średniej, będzie dość cienko ...
Matthew Gunn
1
W odniesieniu do drugiego pytania rozbieżność KL między dwoma różnymi jednolitymi rozkładami jest niezdefiniowana ( jest niezdefiniowany). Podobnie, rozbieżność KL dla dwóch rozkładów empirycznych jest niezdefiniowana, chyba że każda próbka ma co najmniej jedną obserwację o tej samej wartości co każda obserwacja w drugiej próbce. log(0)
jbowman
@jbowman Mała notatka. Chociaż masz rację, że jest niezdefiniowany (lub - ), w teorii informacji zwyczajowo traktuje się log ( 0 ) 0 jako 0 . log(0)log(0)00
Luca Citi,
Podobne pytanie: mathoverflow.net/questions/119752/…
kjetil b halvorsen

Odpowiedzi:

9

Rozbieżność Kullbacka-Leiblera jest zdefiniowana jako więc aby obliczyć (oszacować) to na podstawie danych empirycznych, potrzebowalibyśmy być może niektórych szacunków funkcji gęstości p ( x ) , q ( x )

KL(P.||Q)=-p(x)logp(x)q(x)rex
p(x),q(x) . Tak więc naturalnym punktem wyjścia może być oszacowanie gęstości (a następnie po prostu całkowanie numeryczne). Jak dobra lub stabilna byłaby taka metoda, nie wiem.

pq[0,1][0,10]KL(p||q)=log10KL(q||p)log(1/0)log

Wracając do głównego pytania. Pytanie jest zadawane w bardzo nieparametryczny sposób i nie podano żadnych założeń dotyczących gęstości. Prawdopodobnie potrzebne są pewne założenia. Ale zakładając, że dwie gęstości są proponowane jako konkurujące modele dla tego samego zjawiska, możemy prawdopodobnie założyć, że mają one tę samą dominującą miarę: rozbieżność KL między ciągłym a dyskretnym rozkładem prawdopodobieństwa zawsze byłaby na przykład nieskończonością. Jeden artykuł dotyczący tego pytania jest następujący: https://pdfs.semanticscholar.org/1fbd/31b690e078ce938f73f14462fceadc2748bf.pdf Proponują metodę, która nie wymaga wstępnego oszacowania gęstości, i analizuje jej właściwości.

(Jest wiele innych artykułów). Wrócę i opublikuję kilka szczegółów z tego artykułu, pomysłów.

 EDIT               

Kilka pomysłów z tego artykułu, który dotyczy oszacowania rozbieżności KL z próbkami z absolutnie ciągłych rozkładów. Pokazuję ich propozycję rozkładów jednowymiarowych, ale dają one również rozwiązanie dla wektorów (wykorzystując oszacowanie gęstości najbliższego sąsiada). Aby uzyskać dowody, przeczytaj artykuł!

P.mi(x)=1nja=1nU(x-xja)
UU(0)=0,5P.dodo
re^(P.Q)=1nja=1nlog(δP.do(xja)δQdo(xja))
δP.do=P.do(xja)-P.do(xja-ϵ)ϵ

Kod R dla wersji funkcji rozkładu empirycznego, której potrzebujemy, to

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
}          

Uwaga: rlesłuży do załatwiania sprawy z duplikatami w x.

Następnie oszacowanie dywergencji KL podaje

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              
}

Następnie pokazuję małą symulację:

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

co daje następujący histogram, pokazujący (oszacowanie) rozkład próbkowania tego estymatora:

Rozkład próbkowania estymatora KL

Dla porównania obliczamy dywergencję KL w tym przykładzie przez całkowanie numeryczne:

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 ... różnica jest na tyle duża, że ​​jest tu wiele do zbadania!

kjetil b halvorsen
źródło
5

Rozwijając trochę odpowiedź kjetil-b-halvorsena i przepraszam, że nie komentuję, nie mam reputacji:

  1. Mam wrażenie, że obliczenia analityczne powinny być (bez pomnożenia przez 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. Jeśli mam rację, estymator re^(P.||Q) nie jest zbieżny z rozbieżnością KL, ale zbieżność jest wyrażona jako: re^(P.||Q)-1re(P.||Q). Strzałka przedstawia jako zbieżność.

Po dokonaniu tych dwóch poprawek wyniki wydają się bardziej realistyczne.

ColibriIO
źródło
Dzięki, przyjrzę się temu i zaktualizuję swoją odpowiedź.
kjetil b halvorsen