Calculate the Kullback-Leibler Divergence for these 2 Gamma distributions

Question

I have 2 models

$P \sim Ga(115,1329.914) \\ Q \sim Ga(140,650.6775)$

and I'm looking to calculate the K-L divergence of these 2.

$D_{KL}(P||Q) = \int_\infty ^\infty p(x) log \frac{p(x)}{q(x)}\,dx$

where $p(x)$ and $q(x)$ are the probability densities.

So for my Gamma models I get

$D_{KL}(P||Q) = \int_0 ^\infty f_1(x|a_1,d_1,p_1) log \frac{f_1(x|a_1,d_1,p_1)}{f_2(x|a_2,d_2,p_2)} \,dx$

The 3 parameter gamma distribution is given by

$f(x | a,d,p) = \frac{p}{a^d} \frac{x^{d-1}}{\Gamma(d/p)} exp\Big\{-(\frac{x}{a})^p \Big\}$

I just dont know how to apply this to my problem, as the values I have are just a standard $Ga(\alpha, \beta)$

A further point is that

$\int_0 ^\infty f_1(x|a_1,d_1,p_1) log \frac{f_1(x|a_1,d_1,p_1)}{f_2(x|a_2,d_2,p_2)} \,dx \\ = log \frac{p_1d_2^{d_2}\Gamma(d_2/p_2)}{p_2d_1^{d_1}\Gamma(d_1/p_1)} + \Bigg\{\frac{\psi(d_1/p_1)}{p_1} + log a_1\Bigg\}(d_1-d_2)+\frac{\Gamma(\frac{d_1+p_2}{p_1})}{\Gamma(\frac{d_1}{p_1})} (\frac{a_1}{a_2})^{p_2}-\frac{d_1}{p_1}$

Where $\psi(\cdot)$ is the digamma function

But once again I don't know how to apply this to my original problem as I don't know how to parameterize correctly.

I was wondering if there is an easier way which I am overlooking as this seems very complex.

Edit: Just on a side note, I know that

$D_{KL}(P||Q) \neq D_{KL}(Q||P)$

But are they both still as valid as each other for KL divergence??

Second Edit: From here Kullback–Leibler divergence between two gamma distributions

One of the answers says that KL divergence is a difference in integrals of the form ..... but how do you know that?

Edit:

KL = function(a,b,c,d){
  return(((a-c)/c)*b + log((lgamma(d)*(c^d))/(lgamma(b)*(a^b)))+(b-d)*(log(a)+digamma(b)))
}

KL(202,114186.3,195,119237.3)

Gives answer NaN when it should be close to 1?

Answer 1

The comments by @kjetil b halvorse address the theoretical formula you need. But as you desxcribe, you are having numerical difficulties. That should be resolved as follows:

Use what I generically call "log gamma" function, whatever it may be called in your computing environment. For instance, gammaln in MATLAB https://www.mathworks.com/help/matlab/ref/gammaln.html . In Wolfram , it's called LogGamma http://reference.wolfram.com/language/ref/LogGamma.html .

log gamma function computes log(gamma(x)), but does it in a better way than first computing gamma function, then taking log. This function should avoid overflow (or underflow) from excessively large gamma before taking log. The algebraic rearrangement to use this should be straightforward.

Edit: Given that you want to do the calculations in R per your comment: In R, use lgamma(x) for log(gamma(x)). That should solve your numerical problems. Technically, lgamma is the natural log of the absolute value of the gamma function, but I don't think you need to worry about the abs in your usage.

Edit 2: Here is the complete R function (I also re-wrote the powers in the log to avoid problems).

KL = function(a,b,c,d){
((a-c)/c)*b + lgamma(d) - lgamma(b) + d*log(c) - b*log(a) + (b-d)*(log(a) + digamma(b)) 
}

KL(202,114186.3,195,119237.3) = 9.6433, as computed by the MATLAB analog of this formula (gammaln instead of lgamma, and psi instead of digamma).