Brighton Webs Ltd.
Statistics for Energy and the Environment

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 B2i 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 Excelgammaln(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

dblB(0) = 1

dblB(1) = -1 / 2

dblB(2) = 1 / 6

dblB(4) = -1 / 30

dblB(6) = 1 / 42

dblB(8) = -1 / 30

dblB(10) = 5 / 66

dblB(12) = -691 / 2730

dblB(14) = 7 / 6

dblB(16) = -3617 / 510

dblB(18) = 43867 / 798

dblB(20) = -174611 / 330

'Special Cases

If dblX <= 0 Then

Return Double.NaN

ElseIf dblX = 1

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

dblShift = 7

Else

dblShift = 0

End If

'Stirling Series - elements 1, 2 & 3

dblGL = ((dblX + dblShift) - 0.5) * Math.Log((dblX + dblShift))

dblGL = dblGL - (dblX + dblShift)

dblGL = dblGL + Math.Log(2 * Math.PI) / 2

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

For intI = 1 To 10

dblGL = dblGL + dblB(2 * intI) / 2 / intI / (2 * intI - 1) / Math.Pow((dblX + dblShift), 2 * intI - 1)

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

dblGL = Math.Exp(dblGL)

'Modify value from Stirling series

For intI = 0 To dblShift - 1

dblGL = dblGL / (dblX + intI)

Next

'Back to log(gamma)

dblGL = Math.Log(dblGL)

End If

'Return Value

Return dblGL

End Function

Page update: 30-Aug-2011