需要帮助在 python 中优化三重嵌套二项式 CDF

Question

我必须计算三重求和（嵌套二项式 CDF），方程如下，

它是一个二项式求和，首先从 k = 0 到 C，然后 m = 0 到 k，f = 0 到 Ck，这里 s 是一个函数，它接受 0 和 1 之间的输入并给出 (0, 1) 之间的输出。

我需要帮助找到一种有效的方法来在 python 中执行此操作。现在我使用的是三重循环，它工作得很好，但对于大型 C 来说效率不高。我目前使用的三重循环如下，

这里的 's' 本质上是线性的，但它可以采取任何需要的形状。'r' 也被视为 0.1

import numpy as np
from math import comb

def s(d):
  return d * (0.99 - 0.01) + 0.01

S = 0
C = 100

for k in range(C):
  for m in range(k):
    for f in range(C-k):
      S += comb(C, k) * (((0.1)**(k))*((0.9)**(C-k))) * comb(k, m) * comb(C-k, f) * (2**(-C)) * s((f+m)/C)

我可以使用更有效的方法吗？

编辑： s 只是以下形式的函数，它接受 0 到 1 之间的输入，并根据其形状给出 0.01 到 0.99 之间的输出。在示例代码中它是线性的，但它可以是指数的或其他。

Edit2：函数“s”无法从第一次求和中得出，在我提供的等式中，由于拼写错误，它可能会出现。现在，随着输入 chrslg 的建议，该问题已得到修复。

Answer 1

所以这s只是一个仿射函数。这允许一些优化

\n

从...开始

\n

S=0\nfor k in range(C):\n  for m in range(k):\n    for f in range(C-k):\n      S += comb(C, k) * (((0.1)**(k))*((0.9)**(C-k))) * comb(k, m) * comb(C-k, f) * (2**(-C)) * s((f+m)/C)\n

\n

请注意，此代码中存在一个错误：数学符号中的 \xce\xa3 使用包含的边界（并且在这种公式中，通常有from 0 to n included，这意味着 n+1 次迭代。因为当你在一个中有 n 个袜子时抽屉里，并随机取出其中的数量，然后您可以取出 0 只袜子，1 只袜子，2 只袜子，...最多 n 只袜子，包括:D)。

\n

所以，代码应该是这样的

\n

S=0\nfor k in range(C+1):\n  for m in range(k+1):\n    for f in range(C-k+1):\n      S += comb(C, k) * (((0.1)**(k))*((0.9)**(C-k))) * comb(k, m) * comb(C-k, f) * (2**(-C)) * s((f+m)/C)\n

\n

正如我在评论中提到的，我们可以摆脱循环中不依赖于循环索引的所有内容

\n

S=0; r=0.1\nfor k in range(C+1):\n    Sm=0\n    for m in range(k+1):\n        Sf=0\n        for f in range(C-k+1):\n            Sf += comb(C-k, f) * s((f+m)/C)\n        Sm += Sf * comb(k,m) * 2**(-C)\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

在这里，我只是模仿数学方程而不是代码。一些操作已经以这种方式分解。

\n

当然，我不能像这样直接删除内部循环，因为现在你已经用s(m/C)which is not 替换了which 独立于f, s((m+f)/C)。

\n

但由于s是仿射，我们可以轻松地将其替换为

\n

S=0\nfor k in range(C+1):\n    Sm=0\n    for m in range(k+1):\n        Sf=0\n        for f in range(C-k+1):\n            Sf += comb(C-k, f) * (0.01 + 0.98*m/C + 0.98/C*f)\n        Sm += Sf * comb(k,m) * 2**(-C)\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

因此，一个常数项，一个与 f 项成正比

\n

由于 \xce\xa3comb(n,k) 是 2\xe1\xb5\x8f 且 \xce\xa3k.comb(n,k) 是 n\xc3\x972\xe1\xb5\x8f\xe2\x81\xbb\ xc2\xb9，我们可以将其替换为

\n

S=0\nfor k in range(C+1):\n    Sm=0\n    for m in range(k+1):\n        Sf= (0.01+0.98*m/C)*2**(C-k) + 0.98/C * (C-k)*2**(C-k-1)\n        Sm += Sf * comb(k,m) * 2**(-C)\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

所以这是一个循环！现在我们有 2 个循环而不是 3 个。但还没有结束！

\n

结合Sf=和Sm +=线

\n

S=0\nfor k in range(C+1):\n    Sm=0\n    for m in range(k+1):\n        Sm += comb(k,m) * ((0.01+0.98*m/C)*2**(-k) + 0.98/C * (C-k)*2**(-k-1))\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

请再次注意，我们乘以的是comb(k,m)常数项（0.01*2**(-k)或者0.98/C*(C-k)*2**(-k-1)就我们而言，是常数项m），并且与m项成比例。因此，使用相同的\xce\xa3comb(n,k) = 2\xe1\xb5\x8f和\xce\xa3k.comb(n,k) = k.2\xe1\xb5\x8f\xe2\x81\xbb\xc2\xb9公式，我们可以再次重写

\n

S=0\nfor k in range(C+1):\n    Sm=(0.01*2**(-k) + 0.98/C * (C-k)*2**(-k-1))*2**k + (0.98/C*2**(-k))*k*2**(k-1)\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

这是第二个循环！只剩下一张了。尚未结束！

\n

可以简化为

\n

S=0\nfor k in range(C+1):\n    Sm=0.01 + 0.49/C*(C-k) + (0.49/C)*k\n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

甚至

\n

S=0\nfor k in range(C+1):\n    Sm=0.01 + 0.49 - 0.49/C*k + (0.49/C)*k \n    S += Sm*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

所以，显然Sm=0.5

\n

S=0\nfor k in range(C+1):\n    S += 0.5*comb(C,k) * r**k * (1-r)**(C-k)\n

\n

但这\xce\xa3 comb(n,k) a\xe1\xb5\x8fb\xe2\x81\xbf\xe2\x81\xbb\xe1\xb5\x8f是众所周知的 (a+b)\xe2\x81\xbf 的展开。\n所以那就是 0.5\xc3\x97(r + (1-r))\xe1\xb6\x9c = 0.5\xc3\x971\ xe1\xb6\x9c = 0.5

\n

因此，您的代码得到了最终的简化（抱歉，没有比这更好的了！）

\n

S=0.5\n

\n

您可以轻松检查，一旦边界问题得到纠正，无论参数（C 和 r）是什么，您的代码也会返回 0.5。

\n

一些麻木的技巧

\n

由于这应该是一个 numpy 问题，而不是一个数学简化问题，因此以下是如何使用 numpy 加速此类计算的方法。

\n

假设您正在尝试计算内循环

\n

for f in range(C-k+1):\n     Sf += comb(C-k, f) * s((f+m)/C)\n

\n

（更有趣的一个。因此我想知道s，特别是它是否可矢量化。我最初的计划是给出我现在给出的答案，如果s可矢量化，即是否s可以使用两个标量和数组作为参数。但是当我看到它s甚至是仿射时，我不得不改变我的计划，这甚至更好。但是，好吧，让我们假设它s更复杂，例如包含sin，log甚至只是一些高阶多项式这使得我的正式计算变得不可能）

\n

如果我们可以使用C-k+1向量化函数（同样，如果需要的话，可以在整个数组上运行的函数）立即计算值向量，那么我们可以跳过循环for，并使用 numpy 的sum方法求和。现实中仍然是一个for循环。但有一个隐藏的，更重要的是，它发生在 numpy 的优化 C 代码中，即我们缓慢的纯 python 中的节点。

\n

f
我们可以从\n的值向量开始f=np.arange(C-k+1.0)。注意1.0，确保它是浮动的懒惰方法

\n

好消息是你的s（可能是s你能想到的所有，即使是sin我log提到的 或更高阶多项式）都被向量化了。所以我们可以很容易地计算 的向量s((f+m)/C)。\n这只是简单地按字面意思做：如果f是一个向量，那么也是(f+m)/C，所以也是，s((f+m)/C)因为s是向量化的。

\n

所以，最难的部分是梳子部分。

\n

使用 numpy 你可以用cumprod

\n

f=np.arange(C-k+1.0)\nnp.seterr(divide=\'ignore\') # Just because the 1st term I am about to compute is 1/0\ncombRatio=((C-k-f)/f) # That is also a vector, since f is\ncombRatio[0]=1 # neutral cumprod\nmycomb=combRatio.cumprod() #  All \xce\xa0combRatio aka comb(f,C-k)\n

\n

所以，现在我们有一个向量comb(f,C-k)和 一个向量s((f+m)/C)。要获得所有项的向量，我们只需将两者相乘

\n

f=np.arange(C-k+1.0)\nnp.seterr(divide=\'ignore\') # Just because the 1st term I am about to compute is 1/0\ncombRatio=((C-k-f)/f) # That is also a vector, since f is\ncombRatio[0]=1 # neutral cumprod\nmycomb=combRatio.cumprod() #  All \xce\xa0combRatio aka comb(f,C-k)\nterms = mycomb * s((f+m)/C)\n\n# And the \xce\xa3 result is simply\nterms.sum()\n

\n

请注意，如果您愿意使用 scipy，则更简单，因为 scipy 包含向量化comb函数

\n

f=np.arange(C-k+1.0)\nSf=(scipy.specials.comb(C-k,f) * s((f+m)/C)).sum()\n

\n

这一次，我没有数学技巧。我没有简化任何事情。这只是使循环消失的编码技巧。它仍然在那里（在 scipy\'s comb、 numpy\'s +、-和*内.sum()）。现在甚至有 5 个（非嵌套，并且比纯 python 循环快 5 倍以上）

\n