python - 在Python中处理大二项式的求和

Question

我需要计算这个公式:

这是该积分的近似值

不过没关系，其实我只是想计算图1的值与PYTHON，这就是主题所关心的。

K、alpha 和 sigma 是单次计算中的固定值，通常:

• 0 <= k <= 99;
• 阿尔法 = 3;
• 西格玛 = 2。

下面是我如何在 python 中计算这样的求和:

import decimal
from scipy.special import binom

def residual_time_mean(alpha, sigma=2, k=1):
    prev_prec = decimal.getcontext().prec
    D = decimal.Decimal
    decimal.getcontext().prec = 128

    a = float(alpha)
    s = float(sigma)
    sum1 = 0.0
    sum2 = 0.0
    sumD1 = D(0.0)
    sumD2 = D(0.0)

    for i in range(1, k + 1):
        sum1 += binom(k, i) * ((-1) ** (i + 1)) * (s / ((a - 1) * i - 1.0))
        sum2 += binom(k, i) * ((-1) ** (i + 1)) * s / ((a - 1) * i - 1.0)
        sumD1 += D(binom(k, i)) * (D(-1.0) ** (D(i) + D(1.0))) * (D(s) / ((D(a) - D(1.0)) * D(i) - D(1.0)))
        sumD2 += D(binom(k, i)) * (D(-1.0) ** (D(i) + D(1.0))) * D(s) / ((D(a) - D(1.0)) * D(i) - D(1.0))

    decimal.getcontext().prec = prev_prec

    return sum1, sum2, float(sumD1), float(sumD2)

正在运行

for k in [0, 1, 2, 4, 8, 20, 50, 99]:
    print("k={} -> {}".format(k, residual_time_mean(3, 2, k)))

结果是:

k=0 -> (0.0, 0.0, 0.0, 0.0)
k=1 -> (2.0, 2.0, 2.0, 2.0)
k=2 -> (3.3333333333333335, 3.3333333333333335, 3.3333333333333335, 3.3333333333333335)
k=4 -> (5.314285714285714, 5.314285714285714, 5.314285714285714, 5.314285714285714)
k=8 -> (8.184304584304588, 8.184304584304583, 8.184304584304584, 8.184304584304584)
k=20 -> (13.952692275798238, 13.952692275795965, 13.95269227579524, 13.95269227579524)
k=50 -> (23.134878809207617, 23.13390225415814, 23.134078892910786, 23.134078892910786)
k=99 -> (265412075330.96634, 179529505602.9507, 17667813427.20196, 17667813427.20196)

您可以看到，从 k=8 开始，结果有所不同。

例如，当 k=99 时，在除法之前进行乘法会导致 sum1 和 sum2 的结果出现很大差异。

sum1 += binom(k, i) * ((-1) ** (i + 1)) * (s / ((a - 1) * i - 1.0))
sum2 += binom(k, i) * ((-1) ** (i + 1)) * s / ((a - 1) * i - 1.0)

使用十进制不会出现此问题，但结果根本不正确。

在 WolframAlpha 上计算求和

对于 k = 99

(Here is the link for the computation on WolframAlpha)。它给出了33.3159488(...)，而对于我的python函数来说，它是17667813427.20196。我信任 WolframAlpha，因为它进行类似符号计算的操作，实际上它也以分数形式返回实际值。

对于其他k

近似问题(例如，Wolfram 计算的值与 Python 计算的值相差 10^0 或更多的数量级)从 k~=60 开始出现。

此外，使用 scipy.integrate 计算积分(图 2 )会导致类似的近似错误。

问题:

您对处理此计算有什么建议吗？增加小数精度似乎没有帮助。

Answer 1

我自己发现了这个问题:

执行scipy.special.binom(99,50)给出

5.044567227278209e+28

在 WolframAlpha 上计算二项式 (99,50) 时给出

5.0445672272782096667406248628e+28

绝对差异为10^12数量级。

这就是为什么 python 函数的结果对于高 k 值来说是不可靠的。所以我需要改变二项式的计算方式。