Γ分布函数

Γ分布函数算法新解及其应用

李世才吴戈堂林莺

(广西南宁水利电力设计院)

摘要从Γ函数与不完全Γ函数的恒等关系出发，导出了Γ(α)与lnΓ(α)的精确解析式，并在文[1]的基础上导出Γ分布函数新算法的精确解析式。把迄今Γ(α)、lnΓ(α)和Γ分布函数的计算只能应用各种逼近的近似公式现状，提高到精确解析式的计算水平，并归纳为收敛的级数展式和连分式展式的数值计算，使其算法统一成为现实。用新算法的通用数学模型设计的电算程序，对实际工程的计算和文[2～4]中的全部算例及文[5，6]中的有关数表进行了验证比较，结果表明新算法更为优越。

关键词Γ函数Γ分布函数算法新解精确解析式数学模型数值计算。

本文于1996年6月15日收到，广西自然科学基金资助项目，桂科［自］9912010.

统计学、分子结构论、特殊函数、工程水文分析与计算、水文学的汇流计算等应用和研究领域，经常会遇到Γ函数、Γ分布函数和其逆函数的数值计算问题。文[1]对几种常见的数值积分法进行了分析比较，并给出Γ分布函数通用算法的综合解析表达式

(1)

式中α为参变量(α>0)，x为自变量(x≥0)，T的表达式为

(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-x S(α,x) (5)

Γ(α,x)=xαe-x U(α,x) (6)

式中S(α,x)和U(α,x)分别为一无限幂级数展式和连分式展式，且有

(7)

(8)

式(8)中右边式子的分母是一无限连分式的一种简记法[1，8]。

2 Γ(α)和lnΓ(α)的算法新解

2.1 lnΓ(α)的双精度算法

计算lnΓ(α)的斯特林公式的渐近表达式为[1，8，

10]

(9)

为了使lnΓ(α)的计算得到双精度的数值结果，将式(9)和号部分取其前10项，并加整理得

(10)

式(9)和式(10)中的A为[1]

当α≥7时

当0<α<7时，取m=α+[7-α]，求出A后将m+1赋值于上述两式右边

表达式中的变量α，而左边lnΓ(α)中的α保持不变(仍为原值)，

这里的[]为取整符号。

B

2n

为第2n个伯努力数，C

n

与斯特林级数的项数和伯努力数B

2(n+1)

有关，

系数C

n

的取值见文[7]中的表1.应用恒等式Γ(α)=e lnΓ(α)，又可得到计算Γ(α)的双精度算法公式。

文[4]中计算Γ(α)和lnΓ(α)的双精度算法，其算法过程和有关系数的取值都显冗长。

2.2 Γ(α)和lnΓ(α)新算法的精确解析式

式(9)是一渐近级数，而不是收敛级数，式(10)也只能得到双精度数值结果的近似式，它们都不是精确解析式。要得到Γ(α)和lnΓ(α)新算法的精确解析式，可从恒等式(3)入手。

因为式(3)为一恒等式，所以对任意的x≥0都成立。通过大量的数值实验分析发现，当式(3)中的x取值α+1时，并用式(7)和式(8)计算其中的级数展式和连分式展式，则具有很好的收敛性和数值稳定性，还能保证计算结果的可靠性。将x=α+1代入式(3)，并由式(5)和式(6)，可得

Γ(α)=γ(α,α+1)+Γ(α,α+1)

(11)

=(α+1)αe-(α+1)[S(α,α+1)+U(α,α+1)]

=eαln(α+1)-(α+1)[S(α,α+1)+U(α,α+1)]

对式(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.