Brighton Webs Ltd.
Statistics for Energy and the Environment
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)).
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..
Public Function GamLog(ByVal dblX As Double) As Double
Dim intI As Integer 'Loop Counter
Dim dblGL As Double 'Current Value of log(gamma(x))
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
If dblX <= 0 Then
ElseIf dblX = 1
ElseIf dblX = 2 Then
'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
'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
'If x was less than shift, modify the result
'obtained from the Stirling series.
If dblShift > 0 Then
|Page update: 30-Aug-2011|