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

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 $\mathcal N (0,1)$ i $\mathcal 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?

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

$$\DeclareMathOperator{\KL}{KL} \KL(P || Q) = \int_{-\infty}^\infty p(x) \log \frac{p(x)}{q(x)} \; dx$$ $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.

$p$$q$$[0,1]$$[0,10]$$\KL(p || q)=\log 10$$\KL(q || p)$$\infty$$\log(1/0)$$\log \infty$

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_e(x) = \frac1{n}\sum_{i=1}^n U(x-x_i)$$ $U$$U(0)=0.5$$P_c$$c$ $$\hat{D}(P \| Q) = \frac1{n}\sum_{i=1}^n \log\left(\frac{\delta P_c(x_i)}{\delta Q_c(x_i)}\right)$$ $\delta P_c = P_c(x_i)-P_c(x_i-\epsilon)$$\epsilon$

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:

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!

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

2. Jeśli mam rację, estymator $\hat D(P||Q)$ nie jest zbieżny z rozbieżnością KL, ale zbieżność jest wyrażona jako: $\hat D(P||Q)-1 \to D(P||Q)$. Strzałka przedstawia jako zbieżność.

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