处理 Python 中大二项式的求和
Handle summation of big binomials in Python
我需要计算这个公式：
是这个积分的近似值
不过没关系，其实我只是想用PYTHON计算图1的值，这就是题目所关心的。
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开始结果不一样
在除法之前进行乘法会导致 sum1 和 sum2 的结果相差很大，例如 k=99。
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，因为它可以进行类似符号计算的操作，实际上它还 returns 分数形式的实际值。
其他k
近似问题（例如，Wolfram 计算的值与 python 中计算的值相差 10^0 或更多数量级）从 k~=60 开始出现。
此外，用 scipy.integrate 计算积分（图 2）会导致类似的近似误差。
问题：
你对处理这个计算有什么建议吗？增加小数精度似乎没有帮助。
我自己发现的问题：
执行scipy.special.binom(99,50)得到
5.044567227278209e+28
在 WolframAlpha 上计算二项式 (99,50) 时给出
5.0445672272782096667406248628e+28
有10^12数量级的绝对差异
这就是为什么 python 函数的结果对于高 k 值是不可靠的。所以我需要更改二项式的计算方式。
我不明白你为什么在这里涉及一个numpy函数，为什么你要转换成float对象。实际上，对于这个公式，如果您的输入始终是整数，那么只需坚持使用 int 和 fractions.Fraction，您的答案将始终是 精确 。实现您自己的 binom 功能很容易：
In [8]: def binom(n, k):
...: return (
...: factorial(n)
...: // (factorial(k)*factorial(n-k))
...: )
...:
注意，我使用整数除法：//。最后，您的总结：
In [9]: from fractions import Fraction
...: def F(k, a, s):
...: result = Fraction(0, 1)
...: for i in range(1, k+1):
...: b = binom(k, i)*pow(-1, i+1)
...: x = Fraction(s, (a-1)*i - 1)
...: result += b*x
...: return result
...:
结果：
In [10]: F(99, 3, 2)
Out[10]: Fraction(47372953498165579239913571130715220654368322523193013011418, 1421930192463933435386372127473055337225260516259631545875)
基于 wolfram-alpha 这似乎是正确的...
注意，如果说 alpha 可以是非整数，您可以使用 decimal.Decimal 进行任意精度浮点运算：
In [17]: from decimal import Decimal
...: def F(k, a, s):
...: result = Decimal('0')
...: for i in range(1, k+1):
...: b = binom(k, i)*pow(-1, i+1)
...: x = Decimal(s) / Decimal((a-1)*i - 1)
...: result += b*x
...: return result
...:
In [18]: F(99, 3, 2)
Out[18]: Decimal('33.72169506311642881389682714')
提高精度：
In [20]: import decimal
In [21]: decimal.getcontext().prec
Out[21]: 28
In [22]: decimal.getcontext().prec = 100
In [23]: F(99, 3, 2)
Out[23]: Decimal('33.31594880623309576443774363783112352607484321721989160481537847749994248174570647797323789728798446')
我需要计算这个公式：
是这个积分的近似值
不过没关系，其实我只是想用PYTHON计算图1的值，这就是题目所关心的。
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开始结果不一样
在除法之前进行乘法会导致 sum1 和 sum2 的结果相差很大，例如 k=99。
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，因为它可以进行类似符号计算的操作，实际上它还 returns 分数形式的实际值。
其他k
近似问题（例如，Wolfram 计算的值与 python 中计算的值相差 10^0 或更多数量级）从 k~=60 开始出现。
此外，用 scipy.integrate 计算积分（图 2）会导致类似的近似误差。
问题：
你对处理这个计算有什么建议吗？增加小数精度似乎没有帮助。
我自己发现的问题：
执行scipy.special.binom(99,50)得到
5.044567227278209e+28
在 WolframAlpha 上计算二项式 (99,50) 时给出
5.0445672272782096667406248628e+28
有10^12数量级的绝对差异
这就是为什么 python 函数的结果对于高 k 值是不可靠的。所以我需要更改二项式的计算方式。
我不明白你为什么在这里涉及一个numpy函数，为什么你要转换成float对象。实际上，对于这个公式，如果您的输入始终是整数，那么只需坚持使用 int 和 fractions.Fraction，您的答案将始终是 精确 。实现您自己的 binom 功能很容易：
In [8]: def binom(n, k):
...: return (
...: factorial(n)
...: // (factorial(k)*factorial(n-k))
...: )
...:
注意，我使用整数除法：//。最后，您的总结：
In [9]: from fractions import Fraction
...: def F(k, a, s):
...: result = Fraction(0, 1)
...: for i in range(1, k+1):
...: b = binom(k, i)*pow(-1, i+1)
...: x = Fraction(s, (a-1)*i - 1)
...: result += b*x
...: return result
...:
结果：
In [10]: F(99, 3, 2)
Out[10]: Fraction(47372953498165579239913571130715220654368322523193013011418, 1421930192463933435386372127473055337225260516259631545875)
基于 wolfram-alpha 这似乎是正确的...
注意，如果说 alpha 可以是非整数，您可以使用 decimal.Decimal 进行任意精度浮点运算：
In [17]: from decimal import Decimal
...: def F(k, a, s):
...: result = Decimal('0')
...: for i in range(1, k+1):
...: b = binom(k, i)*pow(-1, i+1)
...: x = Decimal(s) / Decimal((a-1)*i - 1)
...: result += b*x
...: return result
...:
In [18]: F(99, 3, 2)
Out[18]: Decimal('33.72169506311642881389682714')
提高精度：
In [20]: import decimal
In [21]: decimal.getcontext().prec
Out[21]: 28
In [22]: decimal.getcontext().prec = 100
In [23]: F(99, 3, 2)
Out[23]: Decimal('33.31594880623309576443774363783112352607484321721989160481537847749994248174570647797323789728798446')