Gamma distribution, incomplete beta function

Question

I need to calculate Gamma cumulative distribution, and it seems this is fairly equivalent to calculating the incomplete beta function. Excel does have an inculded calculator, but I found no trace of the used algorithm. Do any of you know an accurate way to calculate this function? I tried the following, translated into VB.NET from a website, but it gives stupid results:

Function IncompleteBetaFunc(x As Double, a As Double, b As Double) As Double
    If x <= 0 Or x >= 1 Then Return 0
    Dim bt As Double
    bt = Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + a * Math.Log(x) + b * Math.Log(1.0 - x))
    If x < (a + 1.0) / (a + b + 2.0) Then

        Return bt * betacf(a, b, x) / a
    Else
        Return 1.0 - bt * betacf(b, a, 1.0 - x) / b
    End If
End Function

Function betacf(x As Double, a As Double, b As Double) As Double
    Const MAXIT As Integer = 100
    Const EPS As Double = 0.0000003
    Const FPMIN As Double = 1.0E-30
    Dim aa, c, d, del, h, qab, qam, qap As Double
    Dim m, m2 As Integer
    qab = a + b
    qap = a + 1.0
    qam = a - 1.0
    c = 1.0
    d = 1.0 - qab * x / qap
    If (Math.Abs(d) < FPMIN) Then d = FPMIN
    d = 1.0 / d
    h = d
    For m = 1 To MAXIT
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d
        If (Math.Abs(d) < FPMIN) Then d = FPMIN
        c = 1.0 + aa / c
        If (Math.Abs(c) < FPMIN) Then c = FPMIN
        d = 1.0 / d
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d
        If (Math.Abs(d) < FPMIN) Then d = FPMIN
        c = 1.0 + aa / c
        If (Math.Abs(c) < FPMIN) Then c = FPMIN
        d = 1.0 / d
        del = d * c
        h *= del
        If (Math.Abs(del - 1.0) < EPS) Then Exit For
    Next
    Return h
End Function

Thanks!

"this is fairly equivalent to calculating the incomplete beta function" -- that seems incorrect. What you need is the incomplete gamma function, is it not? A web search finds [Incomplete gamma function - ALGLIB](http://www.alglib.net/specialfunctions/incompletegamma.php) -- maybe that's useful, I don't know. — Robert Dodier, Apr 20 '17 at 20:15
@RobertDodier : you're right, it's me getting confused... I meant Student distribution, which is related to Gamma function and Beta function... — Pierre, Apr 21 '17 at 08:31

score 2 · Accepted Answer · answered Apr 20 '17 at 21:19

2

Meta.Numerics includes well-tested and performant code for this any many other special functions. Its incomplete Beta function is documented here. The underlying code can be studied here. It also has a full-on Gamma distribution object, which will give moments, generate random variates, and do other distribution-related stuff in addition to computing the CDF. The package available via NuGet; just search for Meta.Numerics in the VS NuGet interface.

answered Apr 20 '17 at 21:19

David Wright

765
6
15

Thanks. I made a stupid typo in header, I meant Student , not Gamma (got confused, I coded gamma with Lanczos method just before), but I bet I'll find that too on that website. – Pierre Apr 21 '17 at 08:41
PS: I could translate all the needed code into the right language, works fine! :) – Pierre Apr 21 '17 at 14:08
@Pierre The implementation of numerical functions is pretty sensitive to details, so although it might be easy to get pretty close to right, it is often quite difficult to get it exactly right. My advice is to use code that someone else has already written, and, we hope, tested also. – Robert Dodier Apr 21 '17 at 17:06
@RobertDodier : I agree, but VB.NET and C# (the library given by David is in C#) are totally equivalent. That's why I could translate it into VB.NET with no risk of damaging it. ;-) I also tested the precision of the results, and they are more precise than those I used when using VBA built-in functions. – Pierre Apr 23 '17 at 09:04

Gamma distribution, incomplete beta function

1 Answers1