# Calculate the Kullback-Leibler Divergence for these 2 Gamma distributions

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?

• What are $a$, $d$ and $p$? – Carlos Campos Mar 26 '18 at 17:31
• @CarlosCampos I've edited the question to make it more clear. However, if you look at my second edit. Could you possibly tell me how he know the KL divergence is a difference in integrals of the form given. If I know how to show that then I dont need to undestand the above question. – user162934 Mar 26 '18 at 17:36
• I think the key to my question is that if I set $p=1$ then I get the gamma distribution that I want! – user162934 Mar 26 '18 at 17:36
• You will find an answer on page 8 of mast.queensu.ca/~linder/pdf/GiAlLi13.pdf – kjetil b halvorsen Mar 26 '18 at 17:38
• The online calculator is at wolframalpha.com/input/… (for example). – whuber Mar 27 '18 at 12:38

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).

• Per @whuber 's comment, you can do the calculation at the Wolfram link he provided, using LogGamma . – Mark L. Stone Mar 27 '18 at 12:49
• I am using R rather than MatLab or any other. So what your saying is to just take the log of all my gamma's? What about the Digamma function then? – user162934 Mar 27 '18 at 13:14
• 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. digamma is the logarithmic derivative of the gamma function, and not what you want for this. – Mark L. Stone Mar 27 '18 at 23:05
• For some reason the values are returning NaN, which I just dont get – user162934 Apr 2 '18 at 0:02
• I think there may be some confusion concerning the meanings of a,b,c,d. Using the formula I derived at stats.stackexchange.com/a/11668/919 and (with the notation of that answer) setting $a=1/114186.3$, $b=202$, $c=1/119237.3$, and $d=195$, I obtain a value of $0.60753126\ldots$ rather than $9.6433.$ – whuber Apr 2 '18 at 20:28