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

函数格式: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

内容版权声明:除非注明,否则皆为本站原创文章。

转载注明出处:https://www.heiqu.com/zzwyps.html