Γ分布函数Word下载.docx
《Γ分布函数Word下载.docx》由会员分享，可在线阅读，更多相关《Γ分布函数Word下载.docx（10页珍藏版）》请在冰豆网上搜索。
（2）
电算实践表明：
一些特殊问题的计算中，需要程序参与计算的各个变量（含常数）都按双精度（16位有效数字）或高精度（任意指定精度）运行，而计算Γ（α）和lnΓ（α）的各种逼近算法公式[1～5]最多只能求得10至12位有效数字，这样即使应用双精度计算P（α,x），最多也只能达到与Γ（α）或lnΓ（α）同样的精度，并且有时会加速计算过程的误差传播和积累，从而导致死循环、迭代过程不收敛、计算结果失真等不良的现象。
要解决这些问题，可以将Γ（α）和lnΓ（α）算法公式中的系数扩充[4，7]，但文[4]中的系数较为冗长。
由于应用式
（1）解决实际问题的计算时，通常只能应用各种逼近算法的近似公式[1～14]计算其中的lnΓ（α），所以上述式
（1）还不是精确的解析表达式，不能使求解的精度和算法程序实现达到统一。
通过对文[1]中有关公式进行恒等变换，找到了满足上述要求的精确解析式，把它归纳为收敛的幂级数展式和连分式展式的数值计算，并设计成一标准的子程序，已大量地应用于实际工程的数值计算中，取得了满意的效果。
综上可知，提高式
（1）的计算精度，简化编程手续，出路在于建立具有统一算法基础的Γ（α）和lnΓ（α）精确解析式，寻求Γ分布函数算法新解，以取代现行Γ（α）、lnΓ（α）和Γ分布函数算法公式的精度和算法程序的不统一性，这无疑会有理论价值和实用意义。
1数学原理
1.1特殊函数定义
根据文[1，8，9]，下列含参变量α的定积分式
分别称为Γ函数、不完全Γ函数、普利姆函数、Γ分布函数、余Γ分布函数，它们是数学里的特殊函数，都为高等超越函数。
前两个是最基本的特殊函数，第一个也被称为完全Γ函数，第三个是另一不完全Γ函数，并有如下恒等式
Γ（α）=γ（α,x）+Γ（α,x）
（3）
P（α,x）+Q（α,x）=1
（4）
1.2不完全Γ函数的级数展式和连分式展式
根据文[1]中的式（14）与式（22），可得
γ（α,x）=xαe-xS（α,x）
（5）
Γ（α,x）=xαe-xU（α,x）
（6）
式中S（α,x）和U（α,x）分别为一无限幂级数展式和连分式展式，且有
（7）
（8）
式（8）中右边式子的分母是一无限连分式的一种简记法[1，8]。
2Γ（α）和lnΓ（α）的算法新解
2.1lnΓ（α）的双精度算法
计算lnΓ（α）的斯特林公式的渐近表达式为[1，8，10]
（9）
为了使lnΓ（α）的计算得到双精度的数值结果，将式（9）和号部分取其前10项，并加整理得
（10）
式（9）和式（10）中的A为[1]
当α≥7时
当0<
α<
7时，取m=α+[7-α]，求出A后将m+1赋值于上述两式右边表达式中的变量α，而左边lnΓ（α）中的α保持不变（仍为原值），这里的[]为取整符号。
B2n为第2n个伯努力数，Cn与斯特林级数的项数和伯努力数B2（n+1）有关，系数Cn的取值见文[7]中的表1.应用恒等式Γ（α）=elnΓ（α），又可得到计算Γ（α）的双精度算法公式。
文[4]中计算Γ（α）和lnΓ（α）的双精度算法，其算法过程和有关系数的取值都显冗长。
2.2Γ（α）和lnΓ（α）新算法的精确解析式
式（9）是一渐近级数，而不是收敛级数，式（10）也只能得到双精度数值结果的近似式，它们都不是精确解析式。
要得到Γ（α）和lnΓ（α）新算法的精确解析式，可从恒等式（3）入手。
因为式（3）为一恒等式，所以对任意的x≥0都成立。
通过大量的数值实验分析发现，当式（3）中的x取值α+1时，并用式（7）和式（8）计算其中的级数展式和连分式展式，则具有很好的收敛性和数值稳定性，还能保证计算结果的可靠性。
将x=α+1代入式（3），并由式（5）和式（6），可得
Γ（α）=γ（α,α+1）+Γ（α,α+1）
=（α+1）αe-（α+1）[S（α,α+1）+U（α,α+1）]
=eαln（α+1）-（α+1）[S（α,α+1）+U（α,α+1）]
（11）
对式（11）两边取自然对数，可得
lnΓ（α）=αln（α+1）-（α+1）+ln[S（α,α+1）+U（α,α+1）]
（12）
上述式（11）和（12）就是分别求解Γ（α）与lnΓ（α）的一种新算法的精确解析式，其中S（α,α+1）和U（α，α+1）的值是分别将x=α+1代入式（7）和式（8）中的x计算而得，具体的计算表达式这里从略。
3Γ分布函数算法新解
以上已导出求解Γ（α）和lnΓ（α）新算法的精确解析式，可通过式（11）和式（12）以及式（7）和式（8）的计算而得，将所求的lnΓ（α）代入引言部分中的式
（2），再代入式
（1）后，就可得到计算Γ分布函数P（α,x）的一种新算法的精确解析式。
当α较大时，则Γ分布P（α,x）渐近于标准化正态分布N（0，1）=Φ（u）；
当α→+∞时，则Γ分布P（α,x）就转化为标准化正态分布φ（u）.标准化正态分布函数的表达式为
（13）
令x=u2/2，通过变量代换并加整理得
（14）
从式（14）可看出：
标准化正态分布函数Φ（u）也可用Γ分布函数的表达式表示，其式中P（0.5，0.5u2）就是取α=0.5，x=0.5u2，代入Γ分布函数P（α,x）的表达式而得，sgn（u）为u的符号函数。
因此，当α较大时，Γ分布函数也可用式（14）作近似计算；
对于α→+∞时的P（α,x）和标准化正态分布函数Φ（u）完全可以用式（14）作精确的计算；
当α为正整数时，P（α,x）具有初等函数解析表达式[1]。
有关Γ分布函数截断误差的估计和收敛速度的论述可见文[1]中的4.5与4.6.
当α很大时，则计算P（α,x）的级数和连分式展式的收敛速度都变慢，且α越大，收敛速度就越慢，特别是在x→α+1时，收敛速度也就更慢，为了加快计算速度，可根据P（α,x）≈Φ（u）的关系，用式（14）计算P（α,x）的近似值，其中u用以下两式计算[15]
（15）
或
（16）
式（15）和式（16）根据文[11～13]中的χ2分布与N（0，1）分布的近似关系，再转化成Γ分布与N（0，1）的近似关系[7，12]，然后反推而得。
在多数情况下，应用式（16）比用式（15）代入式（14）求P（α,x）近似值的效果更佳，且α越大，P（α,x）就越接近于Φ（u），在此种情况下求解P（α,x）的计算表达式，也见文[15]。
综上可知，本文所介绍的5个特殊函数都可用新算法的解析式来表示，而计算P（α,x）的关键是需要计算出当α固定时，不同x的式（7）和式（8）的级数展式和连分式展式的数值。
因此，可以根据所用程序设计语言的特点，将式（7）和式（8）的算法设计成标准子程序（或过程），供解题时方便地调用。
对于式（7）和式（8）两式的详细推导过程和算法设计可分别参考文[1，2～5，8，13，14]。
4数值实验分析与算法应用
4.1数值实验
应用所介绍的新算法对Γ（α）、lnΓ（α）、P（α,x）等特殊函数进行数值计算，计算程序全按双精度运行，并将计算结果与文[2～6]中的全部算例和有关数表进行比较，部分结果见表1和表2.
表1Γ（α）与lnΓ（α）数值计算比较
计算方法
α=0.5
α=1.5
Γ（α）
lnΓ（α）
真值
1.772453850905516
0.5723649429247001
0.8862269254527580
-0.1207822376352452
文[4]双精度
-0.1207822376352453
本文双精度
1.772453850905517
0.5723649429247004
0.8862269254527582
-0.1207822376352450
新算法
注：
Γ（0.5）和Γ（1.5）的精确解析表达式是
和
，用新算法计算时，程序中的控制精度取绝对误差限ε<
10-18。
表2Γ分布函数P（α,x）数值计算比较
α=0.1，x=0.031623
α=1,x=5
α=11,x=16.58312
α=41,x=44.82187
文献[2]
0.742026
0.993262
0.940427
0.735970
双精度算法
0.7420268385459224
0.9932620530009146
0.9404266190477025
0.735970933014524
0.7420268385459223
注：
文献[2]的控制精度取ε<
3×
10-7，新算法的控制精度取ε<
10-17。
lnΓ（α）双精度算法公式为有限项的近似解析式，lnΓ（α）新算法公式为无限项的精确解析式。
因此，可以肯定在计算精度相同的条件下，一般后者比前者的计算速度稍慢一些，Γ分布函数P（α,x）的计算也是如此。
以下表3仅列出应用两种算法计算文[6]中P—Ⅲ型分布曲线Φp（或Kp）值表计算时间的比较结果。
表3Φp（或Kp）值表两种算法的计算时间比较
计算项目
控制精度：
ε<
10-9,ε′<
10-6
10-16,ε′<
10-13
Φp值表
00∶00∶19
00∶00∶20
00∶00∶38
00∶00∶39
Kp值表
00∶00∶37
00∶01∶13
00∶01∶14
ε为P（α,tp）的控制精度，ε′为tp迭代计算的控制精度，关于Φp值的计算方法和本表的一些计算条件，可参考文[5，7]，本表的计算时间是用奔腾Ⅱ（300MHz）计算机重新计算时而统计的。
4.2数值结果分析
升级会员 