Gamma Function

The gamma function has many applications in statistics and other "special functions" (e.g. the Beta function) can be derived from it.  It is a useful function to have as part of programming toolkit. 

It is often referred to as the non-integer version of the factorial function because of the relationship:

On this page we are only concerned with real numbers greater than zero.  As the value of the gamma function increases rapidly with increasing argument values, it is often convenient to work with an plot the log of the function as in the graph below:

The code is based on the Stirling series:

Where B_2i is the relevant Bernouli number, thus for N=4,

The Stirlng formula provides reasonable accuracy for values of x >= 7.  A property of the Gamma function is:

This relationship can be written as:

Thus for values less than 7 we can calculate log(Γ(x+7)), then derive log(Γ(x)).

Comment

As with all numerical methods, the user should ensure that the chosen algorithm provides adequate levels of accuracy and efficiency appropriate to the application.  There are many algorithms for the gamma function, the Stirling series was selected because there are many references to it in the literature and it provides an acceptable level of accuracy for many general applications.  The code has been written such that there is a reasonably clear relationship to underlying method rather than computational efficiency.

Test Data from Excel

The table below contains test data obtained by using Excel's gammaln function, the code below returns the same numbers to six d.p..

x	Excel gammaln(x)
0.914	0.055994
1.513	-0.120229
2.939	0.637598
4.375	2.282091
5.906	4.627922
9.077	10.769780
47.802	136.038708
48.444	138.518963
56.214	169.187366
82.891	281.993369

VB.Net code

Public Function GamLog(ByVal dblX As Double) As Double

'Declarations

Dim intI As Integer 'Loop Counter

Dim dblGL As Double 'Current Value of log(gamma(x))

'Bernoulli Numbers

Dim dblB(20) As Double

'Special Cases

If dblX <= 0 Then

Return Double.NaN

ElseIf dblX = 1 Then

Return 0

ElseIf dblX = 2 Then

Return 0

End If

'Shift

'Values less than the shift value will be incremented by

'the value of shift and the value obtained from the

'Stirling series modified.

Dim dblShift As Double = 7

If dblX < dblShift Then

Else

End If

'Stirling Series - elements 1, 2 & 3

'Stirling Series - elements 4......

For intI = 1 To 10

Next

'If x was less than shift, modify the result

'obtained from the Stirling series.

If dblShift > 0 Then

'change from log(gamma) to gamma

'Modify value from Stirling series

For intI = 0 To dblShift - 1

Next

'Back to log(gamma)

End If

'Return Value

Return dblGL

End Function

Page updated: 22-Jul-08

For more information: info@brighton-webs.co.uk