Gamma 分布,不完整的 beta 函数
Gamma distribution, incomplete beta function
我需要计算 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
谢谢!
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;只需在 VS NuGet 界面中搜索 Meta.Numerics。
我需要计算 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
谢谢!
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;只需在 VS NuGet 界面中搜索 Meta.Numerics。