title: Kullback-Leibler divergence
body: <p>I have written a function that computes the Kullback-Leibler divergence from N(mu2, sigma2) to N(0, 1). </p>
  
  <pre><code>mu1 &lt;- 0
  sigma1 &lt;- 1
  f &lt;- function(mu2, sigma2)
  {
        g &lt;- 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)   
  } 
  </code></pre>
  
  <p>For example, the KL divergence from N(5, 1) to N(0, 1) is</p>
  
  <pre><code>&gt; f(5, 1)
  [1] 12.5
  </code></pre>
  
  <p>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). </p>
  
  <p>My question is about the <strong>KLdiv</strong> function from the flexmix package. Why doesn't it yield the same result ? What does it actually compute ?</p>
  
  <pre><code>&gt; library(flexmix)
  &gt; x &lt;- seq(-4, 12, length=200)
  &gt; y &lt;- cbind(norm1=dnorm(x, mean=0, sd=1), norm2=dnorm(x, mean=5, sd=1))
  &gt; KLdiv(cbind(y))
           norm1    norm2
  norm1 0.000000 7.438505
  norm2 7.438375 0.000000
  </code></pre>
  
  <p>Instead of using KLdiv, what do you think of the following procedure :</p>
  
  <pre><code>&gt; x &lt;- rnorm(1000)
  &gt; dist &lt;- mean(dnorm(x, mean=0, sd=1, log=TRUE)) - 
  + mean(dnorm(x, mean=5, sd=1, log=TRUE))
  &gt; print(dist)
  [1] 12.40528
  </code></pre>
  
  <p>???</p>
  
  <p>Thank you in advance !</p>