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:

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.

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("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");
}
catch (Exception ex)
{
Console.WriteLine("Fatal: " + ex.Message);
}
} // 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)
{

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

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
```