Thursday, March 31, 2016

Log Gamma approximation

By the standard I set out yesterday, C#'s math library is less than decent. I was a bit surprised that even the stats library didn't have a function for the log of gamma. Fortunately, there's no secret to the approximation. Most implementations use a variant of the Stirling Approximation, which can be super accurate if you want it to be. I don't really need much more than a few decimal places, so I coded up the Nemes version which only uses a few terms and, according to my tests, gives 6-digits when z is greater than 1.5. The CISS confidence intervals are pretty meaningless without at least 1 block sampled, so the lower bound on the range I care about is 2.0. You can read the wiki page (or the source article, if you know German) for the mathematical details. I'll simply post the code for my C# brethren:

        // log gamma function using simplified Stirling's formula
        // GergÅ‘ Nemes: New asymptotic expansion for the Gamma function
        // Archiv der Mathematik, August 2010, Volume 95, Issue 2, pp 161-169
        public static double LogGamma(double z)
        {
            return (Math.Log(2.0 * Math.PI) - Math.Log(z)) / 2.0 +
                z * (Math.Log(z + 1.0 / (12.0 * z - 1.0 / (10.0 * z))) - 1.0);
        }


1 comment: