Computing log gamma differences

Statistical computing often involves working with ratios of factorials. These factorials are often too big to fit in a floating point number, and so we work with logarithms. So if we need to compute log(a! / b!), we call software that calculates log(a!) and log(b!) directly without computing a! or  b! first. More on that here. But sometimes this doesn’t work.

Suppose a = 10^40 and b = a – 10^10. Our first problem is that we cannot even compute b directly. Since ab is 30 orders of magnitude smaller than a, we’d need about 100 bits of precision to begin to tell a and b apart, about twice as much as a standard floating point number. (Why 100? That’s the log base 2 of 10^30. And how much precision does a floating point number have? See here.)

Our next problem is that even if we could accurately compute b, the log gamma function is going to be approximately the same for a and b, and so the difference will be highly inaccurate.

So what do we do? We could use some sort of extended precision package, but that is not necessary. There’s an elegant solution using ordinary precision.

Let f(x) = log(Γ(x)). The mean value theorem says that

f(a+1) – f(b+1) = (a-b) f‘(c+1)

for some c between a+1 and b+1. We don’t need to compute b, only ab, which we know is 10^10. The derivative of the log of the gamma function is called the digamma function, written ψ(x). So

log(a! / b!) = 10^10 ψ(c+1).

But what is c? The mean value theorem only says that the equation above is true for some c. What c do we use? It doesn’t matter. Since a and b are so relatively close, the digamma function will take on essentially the same value for any value between a+1 and b+1. Therefore

log(a! / b!) ≈ 10^10 ψ(a).

We could compute this in two lines of Python:

from scipy.special import digamma
print 1e10*digamma(1e40)

Related posts:

How to compute log factorial
Math library functions that seem unnecessary

Tagged with: ,
Posted in Computing, Python, Statistics
4 comments on “Computing log gamma differences
  1. Paul A. Clayton says:

    “the difference between a and b is 30 orders of magnitude”

    Shouldn’t that be something more like: “the difference between a and b is 30 orders of magnitude smaller than a (b is the same order of magnitude–40–as a)”?

  2. John says:

    Thanks. I reworded that part to say “Since a – b is 30 orders of magnitude smaller than a, …”

  3. Lee McCuller says:

    It’s also worth noting that this is the same as using a first-order taylor expansion. Your comment that c does not matter is equivalent to the second-derivative of log-gamma at a or b being much smaller than (a-b)^2 (bounding the error term).

  4. Luis Carvalho says:

    Using the digamma function in this case can be “expensive”. In your example, since b = a – x, with x = 10^10, we can approximate the factorial ratio by

    a!/(a-x)! ~ (a-x)^x = a^x * (1-x/a)^x ~ a^x * (1 – x^2 / a)

    since x/a is small. In fact, 1 – x^2 /a ~ 1 since x^2 / a is still very small, so

    log(a!/(a-x)!) ~ x * log(a) = 10^10 * log(a)

    which is somehow obvious in hindsight since the digamma function can be roughly approximated by the (natural) log function… :)

3 Pings/Trackbacks for "Computing log gamma differences"
  1. [...] Computing log gamma differences John D. Cook, al solito ::: The Endeavour [...]

  2. [...] The Endeavour, John Cook explains how to computationally deal with logarithms of large [...]

  3. [...] When I first heard of software that could do extended precision calculation, I thought it would be very useful. Years later, I haven’t had much use for it. When I’ve though I needed extended precision, I’ve usually found a better way to solve my problem using ordinary precision and a little pencil-and-paper math. Avoiding extended precision calculation has caused me to understand my problems better. (Here’s a recent example.) [...]