Calculating the Gamma Function using Estimated Integration

The gamma function is a strange math function that pops up in many areas of probability and combinatorics. The gamma function extends the idea of the factorial function.

The factorial(4) = 4*3*2*1 = 24. And factorial(3) = 3*2*1 = 6. And so on. But what if you wanted factorial(3.5)?

The math behind the gamma function is very deep. It turns out that gamma(n) = factorial(n-1) if n is an integer (whole number). So gamma(4) = factorial(3) = 6, and gamma(3) = factorial(2) = 2. Therefore gamma(3.5) should be something between 2.0 and 6.0.

The math definition of gamma(t) is:

GammaFunctionEquation

So, gamma(3.5) is the area under the curve of x^(2.5) * e^(-x), from 0 to infinity. Just for fun I thought I’d write a short computer program to estimate the value of gamma, by adding up the areas of many small rectangles. This technique is known as the Riemann sum technique.

I coded up a demo using the C# language (below). To test my demo I found an online calculator that could do the gamma function.

GammaFunctionEstimation

My estimated values agree pretty closely with the online calculator’s values. Well, I’m not sure what the point is. I don’t see any practical value in estimating gamma using a Riemann sum, but the exercise was interesting for me and gave me some new insights into the gamma function.

using System;
namespace GammaRiemann
{
  class GammaRiemannProgram
  {
    static void Main(string[] args)
    {
      try
      {
        Console.WriteLine("Begin Gamma function demo\n");
        Console.WriteLine("Estimate the Gamma function using");
        Console.WriteLine("approx integration adding up areas");
        Console.WriteLine("of many small rectangles \n");

        double[] tValues = new double[] { 2.0, 2.5, 3.0, 3.1,
          3.2, 3.3, 3.4, 3.5, 4.0, 4.5, 5.0, 5.5 };
        double[] gValues = new double[] { 1.0, 1.329340388179,
          2.0,2.1976202783925, 2.4239654799354, 2.683437381956,
          2.98120642681, 3.323350970448, 6.0, 11.631728396567,
          24.0, 52.34277778455 };

        for (int i = 0; i < tValues.Length; ++i)
        {
          double t = tValues[i];
          double gExternal = gValues[i];
          double gCalc = Gamma(t);
          Console.WriteLine("t = " + t.ToString("F2"));
          Console.WriteLine("Gamma(t) from reference = " +
            gExternal.ToString("F14"));
          Console.WriteLine("Gamma(t) my calculation = " +
            gCalc.ToString("F14"));
          Console.WriteLine("");
        }

        Console.WriteLine("\nEnd demo\n");
        Console.ReadLine();
      }
      catch (Exception ex)
      {
        Console.WriteLine("Fatal: " + ex.Message);
        Console.ReadLine();
      }
    } // Main

    public static double Gamma(double t)
    {
      double width = 0.002;
      double mid = width / 2;
      double area = 0.0;

      for (int i = 0; i < 100000; ++i)
      {
        double height = GammaDensityFunc(t, mid);

        if (height < double.MinValue + 1.0E-30)
        {
          Console.WriteLine("ht = " + height +
            " at i = " + i);
          Console.ReadLine();
        }

        double a = width * height;
        area += a;
        mid = mid + width;

      }
      return area;
    }

    public static double GammaDensityFunc(double t, double x)
    {
      return Math.Pow(x, t - 1) * Math.Exp(-x);
    }

  } // Program 

} // ns
This entry was posted in Machine Learning. Bookmark the permalink.