Gamma 分布，不完全 beta 函数

Question

我需要计算 Gamma 累积分布，这似乎与计算不完全 beta 函数相当。 Excel 确实有一个包含的计算器，但我没有发现使用过的算法的痕迹。 你们中有人知道计算此函数的准确方法吗？ 我尝试了以下，从网站翻译成 VB.NET，但它给出了愚蠢的结果：

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

谢谢！

Answer 1

Meta.Numerics 包含用于此任何其他特殊功能的经过良好测试和高性能的代码。其不完整的 Beta 功能记录在 here。底层代码可以研究here。它还有一个完整的 Gamma distribution 对象，除了计算 CDF 之外，它还会给出时刻、生成随机变量以及做其他与分布相关的事情。该软件包可通过NuGet 获得；只需在 VS NuGet 界面中搜索 Meta.Numerics 即可。