绘制具有大二项式系数的函数答案

【问题标题】：Plot function with large binomial coefficients绘制具有大二项式系数的函数
【发布时间】：2013-12-13 18:41:00
【问题描述】：

我想绘制一个涉及二项式系数的函数。我的代码是

#!/usr/bin/python
from __future__ import division
from scipy.special import binom
import matplotlib.pyplot as plt
import math

max = 500
ycoords = [sum([binom(n,w)*sum([binom(w,k)*(binom(w,k)/2**w)**(4*n/math.log(n)) for k in xrange(w+1)]) for w in xrange(1,n+1)]) for n in xrange(2,max)]

xcoords = range(2,max)

plt.plot(xcoords, ycoords)
plt.show()

不幸的是，这永远不会终止。如果您将最大值减少到 40，则表示它工作正常。有没有办法绘制这个函数？

我还担心 scipy.special.binom 可能无法给出准确的答案，因为它似乎在浮点中工作。

【问题讨论】：

您的代码似乎运行良好，我按原样运行它，虽然它确实花费了大量时间，但当您使用 500 个术语时，这是可以预料的。另外，剧情还不错。我不清楚问题是什么。
@TroyRockwood 我认为 binom(n,w) 只能精确到 53 位，所以你得到的图片不会是正确的吧？
这是一个不同于“它永远不会终止”或“你能画出来”的问题，显然结果将受到所使用的可变精度的准确性的限制。也许一点点说明问题的应用会有所帮助。

标签： python numpy scipy

【解决方案1】：

您可以通过使用 numpy 计算内循环来获得显着的加速。首先将max 更改为N（因为max 是内置函数）并将您的函数分解为更小、更易于管理的块：

N = 500
X = np.arange(2,N)

def k_loop(w,n):
    K = np.arange(0, w+1)
    return (binom(w,K)*(binom(w,K)/2**w)**(float(n)/np.log(n))).sum()

def w_loop(n):
    v = [binom(n,w)*k_loop(w,n) for w in range(1,n+1)]
    return sum(v)

Y = [w_loop(n) for n in X]

使用N=300 作为测试它需要3.932s 使用numpy 代码，但81.645s 使用你的旧代码。我什至没有计时 N=500 案例，因为你的旧代码花了这么长时间！

值得指出的是，您的函数基本上是指数增长的，可以这样近似。您可以在 semilogx 图中看到这一点：

【讨论】：

谢谢。我还稍微编辑了公式以使其更合理（请参阅新的 4 因子）。