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