概率统计Γ(a,b)函数伽马函数

关于伽马函数

1、

函数格式:GammaP(a,b,value)

a,b:服从Γ(a,b)分布的函数.

value:代表服从Γ(a,b)分布的函数在x≥value时的概率值.

函数功能:概率统计的gamma函数在x≥value时的概率值.

2、

函数格式:Gammln(a)

a:Γ(a)当中的值

当a≤1的时候,采用多项式逼近的方式来求得Gammln(a)

函数功能:函数成功返回ln(Γ(a))的值.

3、

函数格式:Gamma(a)

a:Γ(a)当中的值

函数原理:

Γ(a+1)=(a+1)*Γ(a)

Γ(a)*Γ(1-a)=π/sin(π*a)

Γ(a)=Exp(Gammln(a))

函数功能:函数成功返回Γ(a)的值.

4、

函数格式:Gamma(a,value)

a:Γ(a,x)当中的a,a>0

value:Γ(a,value)当中的x,value>0

注意:上面a、b参数的值必须都大于0,且都必须小于171,否则程序会崩溃

函数功能:返回不完全伽马函数Γ(a,value)的值.

例子:

问:已知某元件寿命服从X~Γ(2,0.5)分布[单位小时],随便取一个元件,求该元件寿命大于4小时的概率.

答:执行本函数GammaP(2,0.5,4)后返回0.40600584970982即此值即为其概率值.

源代码:

Public Function Gamma(ByVal a As Double, ByVal value As Double) As Double '注意a必须大于0
        '不完全伽马函数求值P(a,value)
        If value = 0 Then
            Return 0
        ElseIf value > 1.0E+35 Then
            Return 1
        ElseIf value < 0 Then
            value = -value
            Return -Gamma(a, value)
        End If
        Dim i As Integer
        Dim d, p, p0, p1, q, q0, q1, s, s1, qq As Double
        q = Math.Log(value)
        q *= a
        qq = Math.Exp(q)
        If value < a + 1 Then
            p = a
            d = 1 / a
            s = d
            For i = 1 To 100
                p += 1
                d = d * value / p
                s += d
                If Math.Abs(d) < Math.Abs(s) * 0.0000001 Then
                    s = s * Math.Exp(-value) * qq / Gamma(a)
                    Return s
                End If
            Next
        Else
            s = 1 / value
            p0 = 0
            p1 = 1
            q0 = 1
            q1 = value
            For i = 1 To 100
                p0 = p1 + (i - a) * p0
                q0 = q1 + (i - a) * q0
                p = p0 * value + i * p1
                q = q0 * value + i * q1
                If q <> 0 Then
                    s1 = p / q
                    p1 = p
                    q1 = q
                    If Math.Abs((s1 - s) / s1) < 0.0000001 Then
                        s = s1 * Math.Exp(-value) * qq / Gamma(a)
                        Return 1 - s
                    End If
                    s = s1
                End If
                p1 = p
                q1 = q
            Next
        End If
        s = 1 - s * Math.Exp(-value) * qq / Gamma(a)
        Return s
    End Function

Public Function Gammln(ByVal value As Double) As Double
        Dim temp(2) As Double
        Dim i As Integer
        Dim data() As Double = {76.180091729471457, -86.505320329416776, 24.014098240830911, -1.231739572450155, 0.001208650973866179, -0.000005395239384953}
        temp(0) = value
        temp(1) = value + 5.5
        temp(1) = temp(1) - (value + 0.5) * Math.Log(temp(1))
        temp(2) = 1.0000000001900149
        For i = 0 To 5
            temp(0) += 1
            temp(2) += data(i) / temp(0)
        Next
        temp(0) = Math.Log(2.5066282746310007 * temp(2) / value) - temp(1)
        Return temp(0)
    End Function

    Public Function Gamma_Call(ByVal value As Double) As Double
        If value < 1 Then
            Return Math.Exp(Gammln(value))
        ElseIf value = 1 Then
            Return 1
        End If
        value -= 1
        Return value * Gamma_Call(value)
    End Function

    Public Function Gamma(ByVal value As Double) As Double '概率统计里的gamma函数.注意value的值不要大过171且不能且当value<=0时《sin(value*pi)<>0
        Dim ret As Double = 0
        If value > 171 Then
            IsTrue = False
        ElseIf value > 0 Then
            ret = Gamma_Call(value)
        Else
            ret = Math.Sin(Math.PI * value)
            If ret = 0 Then
            Else
                ret = Math.PI / ret
                ret = ret / Gamma_Call(1 - value)
            End If
        End If
        Return ret
    End Function

    Public Function GammaP(ByVal a As Double, ByVal b As Double, ByVal value As Double) As Double
        If a <= 0 Or a > 17 Or b < 0 Or value < 0 Then
            Return -1
        End If
        Dim temp As Double = Gamma(a)
        temp = Math.Pow(b, a) / temp
        '*****************采用复化后的2点高斯积分方法
        Dim x As Double = 0
        Dim b_a_ As Double = value
        Dim h_ As Double = b_a_ / 1000
        Dim h2_ As Double = h_ / 2
        Dim h23_ As Double = h2_ / Math.Sqrt(3)
        Dim n_ As Integer = 1000
        Dim ret_ As Double = 0
        Dim x_ As Double = 0
        a -= 1
        While n_ > 0
            x = x_ + h2_ + h23_
            ret_ += Math.Exp(-b * x) * x ^ a
            x = x_ + h2_ - h23_
            ret_ += Math.Exp(-b * x) * x ^ a
            x_ += h_
            n_ -= 1
        End While
        ret_ *= h_
        ret_ /= 2
        '*******************
        value = 1 - temp * ret_
        Return value
    End Function