概率統計Γ(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