CODEWITHSUNDEEP
r
Marco
Marco

Reputation: 9614

Kullback-Leibler divergence

I have written a function that computes the Kullback-Leibler divergence from N(mu2, sigma2) to N(0, 1). 

mu1 <- 0
sigma1 <- 1
f <- function(mu2, sigma2)
{
      g <- function(x)
      {
            (dnorm(x, mean=mu1, sd=sigma1, log=TRUE) -
             dnorm(x, mean=mu2, sd=sigma2, log=TRUE)) *
             dnorm(x, mean=mu1, sd=sigma1)
      }
      return(integrate(g, -Inf, Inf)$value)   
}

For example, the KL divergence from N(5, 1) to N(0, 1) is

> f(5, 1)
[1] 12.5

I am sure that this result is correct because I computed at hand a closed form expression that gives the KL divergence from N(mu2, sigma2) to N(mu1, sigma1).

My question is about the KLdiv function from the flexmix package. Why doesn't it yield the same result ? What does it actually compute ?

> library(flexmix)
> x <- seq(-4, 12, length=200)
> y <- cbind(norm1=dnorm(x, mean=0, sd=1), norm2=dnorm(x, mean=5, sd=1))
> KLdiv(cbind(y))
         norm1    norm2
norm1 0.000000 7.438505
norm2 7.438375 0.000000

Instead of using KLdiv, what do you think of the following procedure :

> x <- rnorm(1000)
> dist <- mean(dnorm(x, mean=0, sd=1, log=TRUE)) - 
+ mean(dnorm(x, mean=5, sd=1, log=TRUE))
> print(dist)
[1] 12.40528

???

Thank you in advance !

Upvotes: 3

Views: 8826

Answers (2)

user2561915
user2561915

Reputation: 21

Check the eps parameter in the manual page ?KLdiv,matrix-method:

> KLdiv(cbind(y),eps=1e-16)
         norm1    norm2
norm1  0.00000 12.49908
norm2 12.49941  0.00000

Upvotes: 1

Seth
Seth

Reputation: 4795

In the last part you write

 x <- rnorm(1000)
 dist <- mean(dnorm(x, mean=0, sd=1, log=TRUE)) - 

   mean(dnorm(x, mean=5, sd=1, log=TRUE))

   print(dist)

[1] 12.40528

This is the divergence for a random sample of size 1000. The closed form expression is the limiting value as sample size goes to infinity. If you change your sample size you will get closer. or if you do the same calculation repeatedly you can see that the mean of the estimates is 12.5 like you want.

Upvotes: 6

Related Questions