The KL-distance in R

I would like to compute the KL-distance from two gamma distribution using R.

Here is my开发者_运维知识库 R-code:

theta1 <- 0.2
theta2 <- 2

f <- function(u)
{
     dgamma(u, shape=1/theta1, scale=theta1) * 
      (dgamma(u, shape=1/theta1, scale=theta1, log=TRUE) -
       dgamma(u, shape=1/theta2, scale=theta2, log=TRUE)) 
}

f <- Vectorize(f)
integrate(f, lower=0, upper=Inf)

Do you have any comment on my R-code? Do you think it is the good way to compute the KL-distance?

Any suggestion will be appreciated,

Thx, Marco

I would define all arguments that are used in the function. What I mean is:

my.theta1 <- 0.2
my.theta2 <- 2

f <- function(u, theta1, theta2)
{
     dgamma(u, shape=1/theta1, scale=theta1) * 
      (dgamma(u, shape=1/theta1, scale=theta1, log=TRUE) -
       dgamma(u, shape=1/theta2, scale=theta2, log=TRUE)) 
}


f <- Vectorize(f)
integrate(f, lower=0, upper=Inf, theta1 = my.theta1, theta2 = my.theta2)

Being more explicit prevents "accidents", as your function searches theta1 and theta2 in higher (global) environments (which can get messy if you have this function buried deep inside a program).

Check out the explicit formula here:

http://en.wikipedia.org/wiki/Gamma_distribution#Kullback.E2.80.93Leibler_divergence