Python 不完全伽马函数|极客教程

Python 计算不完全伽马函数，不完全伽马函数具有级数展开式，这意味着需要计算一系列的值然后对它们进行求和。更多相关信息，请访问http://dlmf.nist.gov/8。

Python 不完全伽马函数

这个级数有无穷多个项，这些值最终会变得很小以至于相关性不大。可以建立一个极小值 $\epsilon$ ，并在下一项小于该值时停止计算。
${(-1)}^k$ 删除的计算会产生交替的符号：

Python 不完全伽马函数

当 s=1 且 z=2 时，序列中的项如下所示：

Python 不完全伽马函数

在某一时刻，额外增加的每一项都不会显著影响结果。

当回顾累积分布函数 F(x,k) 时，可以考虑使用fractions.Fraction值。自由度 k 除以2后是一个整数。值 x 既可以是Fraction，也可以是float值，很少会是整型值。

在计算 $\gamma$ (s,z)中的项时，值 $\frac{{(-1)}^k}{k!}$ 会包含整数，并且可以用一个合适的Fraction值来表示。然而，整个表达式并不总是一个Fraction对象。值 $z^{s+k}$ 可以是Fraction或float值。如果 s+k 不是整数，则会导致结果为无理数。对于无理数，可以用形式较为复杂的Fraction对象来近似表示它们的值。

在伽马函数中使用Fraction值是可行的，但似乎帮助不大。然而对于完全伽马函数，Fraction对象具有潜在优势。鉴于此，即使Fraction对象能得到无理数值的近似，我们还是会在此实现中使用该对象。

对上面阐述的级数展开的实现如下所示：

from typing import Iterator, Iterable, Callable, cast
def gamma(s: Fraction, z: Fraction) -> Fraction:

    def terms(s: Fraction, z: Fraction) -> Iterator[Fraction]:
        """Terms for computing partial gamma"""
        for k in range(100):
            t2 = Fraction(z**(s+k))/(s+k)
            term = Fraction((-1)**k, fact(k))*t2
            yield term
        warnings.warn("More than 100 terms")

    def take_until(
            function: Callable[..., bool], source: Iterable
        ) -> Iterator:
        """Take from source until function is false."""
        for v in source:
            if test(v):
                return
        yield v

ε = 1E-8
g = sum(take_until(lambda t: abs(t) < ε, terms(s, z)))
# sum() from Union[Fraction, int] to Fraction
return cast(Fraction, g)

这样就定义了一个term()函数用于生成一系列项，其中限制了for语句最多生成100项。也可以使用itertools.count()函数来生成无限项序列，但这里使用带上限的循环似乎更简单一些。

上面还计算了一个可能是无理数的 $z^{s+k}$ ，并基于这个值创建了一个Fraction值。除法 $\frac{z^{s+k}}{s+k}$ 会涉及两个Fraction相除，并生成一个新的Fraction值赋给变量t2。之后变量term的值是这两个Fraction对象的乘积。

此外定义了一个take_until()的函数，它从一个可迭代对象中取值，直到给定的函数为真时结束。一旦该函数为真，可迭代对象便不再生成新值。上面还定义了一个较小的阈值 $\epsilon$ ： $10^{-8}$ 。我们会从term()函数中取值，直到取到的值小于 $\epsilon$ 。这些值的和是对伽马函数部分的近似。

请注意，阈值变量是希腊字母 $\epsilon$ 。Python 3允许变量名使用任何Unicode字符。

可以用如下测试用例验证计算结果是否正常：

Python 不完全伽马函数