I was sitting at home one Saturday morning, waiting for the rain to stop so I could walk my dogs. For reasons that I don’t remember now, I decided to look at some old ANOVA (analysis of variance test) code.

In ANOVA, you must compute the area under the curve of the F-distribution. To do that, you must compute an incomplete regularized beta function. To do that, you must compute a LogBeta function. To do that, you must compute a LogGamma function. Whew!

The LogGamma function cannot be computed exactly, it must be estimated. There are many algorithms to estimate the LogGamma function. I implemented several approximations using Python. Here’s the output of a demo run:

Approximations to the LogGamma(z) function

Setting z = 0.5

Using scipy loggamma()       : 0.572364942925
Using my old Lanczos(5,7 )   : 0.572364942923
Using Wikipedia Lanczos(5,7) : 0.572364942923
Using Nemes approx           : 0.563828887759
Using A&S Stirling(6)        : 0.512589326855

End demo

The most popular algorithm is called the Lanczos approximation. There are different versions of Lanczos. I used the g=5, n=7 version. While I was reviewing my code, I looked up the Wikipedia article on the Lanczos approximation, and noticed that the version in the article used a slightly different version from my old implementation that I’ve been using for many years.

The output shows that the Nemes approximation is significantly less accurate than Lanczos, and the Stirling approximation with 6 terms is quite poor (especially for small z values).

Good fun.

The LogBeta function is old but still interesting (to me anyway). Here is a fascinating old automata made by Henry Phalibois in approximately 1920. The magician points to the door of the closet and it opens to reveal his assistant. He waves and the door closes. He hits the gong and the door opens to reveal the assistant has disappeared. He gongs again and the door to the small red cube opens to reveal the assistant sitting cross-legged.

Demo program. Replace “lt” (less than) with the Boolean operator symbol.

# log_gamma_approximations.py

import numpy as np
import scipy as sp
import math

# -----------------------------------------------------------

def logGammaWikiLanczos(z):
  # Lanczos approximation g=5, n=7, Wikipedia
  coefs = [ 1.0000000001900148240,
    76.180091729471463483, -86.505320329416767652,
    24.014098240830910490, -1.2317395724501553875,
    0.0012086509738661785061, -5.3952393849531283785e-6 ]

  if z "lt" 0.5:
    y = np.pi / (np.sin(np.pi * z) * logGammaWikiLanczos(1-z))
  else:
    z -= 1.0
    x = coefs[0]
    for i in range(1, len(coefs)):
      x += coefs[i] / (z + i)
      t = z + 5 + 0.5
      y = np.sqrt(2 * np.pi) * (t**(z + 0.5)) * np.exp(-t) * x
  return np.log(y)

# -----------------------------------------------------------

def logGammaMyOldLanczos(z):
  # Lanczos approximation g=5, n=7, my old version
  coefs = [ 1.000000000190015,
    76.18009172947146, -86.50532032941677,
    24.01409824083091, -1.231739572450155,
    0.1208650973866179e-2, -0.5395239384953e-5 ]

  logSqrtTwoPi = 0.91893853320467274178
  if z "lt" 0.5:
    return np.log(np.pi / np.sin(np.pi * z)) * \
      logGammaMyLanczos(1.0 - z);
  else:
    zz = z - 1.0
    b = zz + 5.5
    sum = coefs[0]

    for i in range(1,len(coefs)):
      sum += coefs[i] / (zz + i)

    return (logSqrtTwoPi + np.log(sum) - b) + \
      (np.log(b) * (zz + 0.5))

# -----------------------------------------------------------

def logGammaNemes(z):
  a = 0.5 * (np.log(2 * np.pi) - np.log(z))
  b = (12 * z) - (1.0 / 10 * z)
  c = np.log(z + (1 / b))
  d = z * (c - 1)
  return a + d

# -----------------------------------------------------------

def logGammaStirling(z):
  # A&S 6.1.41 (Stirling's approximation)
  x1 = (z - 0.5) * np.log(z)
  x3 = 0.5 * np.log(2 * np.pi)
  x4 = 1 / (12 * z)
  x5 = 1 / (360 * z * z * z)
  x6 = 1 / (1260 * z * z * z * z * z)
  x7 = 1 / (1680 * z * z * z * z * z * z * z)
  # more terms possible
  return x1 - z + x3 + x4 - x5 + x6 - x7

# -----------------------------------------------------------

print("\nApproximations to the LogGamma(z) function ")

z = 0.5
print("\nSetting z = %0.1f \n" % z)

result = sp.special.loggamma(z)
print("Using scipy loggamma()       : %0.12f " % result)

result = logGammaMyOldLanczos(z)
print("Using my old Lanczos(5,7 )   : %0.12f " % result)

result = logGammaWikiLanczos(z)
print("Using Wikipedia Lanczos(5,7) : %0.12f " % result)

result = logGammaNemes(z)
print("Using Nemes approx           : %0.12f " % result)

result = logGammaStirling(z)
print("Using A&S Stirling(6)        : %0.12f " % result)

print("\nEnd demo ")