Γ分布函数.docx

Γ分布函数.docx

《Γ分布函数.docx》由会员分享，可在线阅读，更多相关《Γ分布函数.docx（10页珍藏版）》请在冰豆网上搜索。

Γ分布函数

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

李世才吴戈堂林莺

（广西南宁水利电力设计院）

摘要从Γ函数与不完全Γ函数的恒等关系出发，导出了Γ（α）与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-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Γ（α）

Γ（α）

lnΓ（α）

真值

1.772453850905516

0.5723649429247001

0.8862269254527580

-0.1207822376352452

文[4]双精度

1.772453850905516

0.5723649429247001

0.8862269254527580

-0.1207822376352453

本文双精度

1.772453850905517

0.5723649429247004

0.8862269254527582

-0.1207822376352450

新算法

1.772453850905516

0.5723649429247001

0.8862269254527580

-0.1207822376352452

注：

Γ（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

0.9932620530009146

0.9404266190477025

0.735970933014524

   注：

文献[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∶00∶38

00∶01∶13

00∶01∶14

   注：

ε为P（α,tp）的控制精度，ε′为tp迭代计算的控制精度，关于Φp值的计算方法和本表的一些计算条件，可参考文[5，7]，本表的计算时间是用奔腾Ⅱ（300MHz）计算机重新计算时而统计的。

4.2数值结果分析