统计: Python 中的组合

我需要在 Python 中计算组合数(nCr) ，但是在 math、 numpy或 stat库中找不到这样做的函数。类似于这种类型的函数:

comb = calculate_combinations(n, r)

我需要的可能组合的数量，而不是实际的组合，所以 itertools.combinations不感兴趣我。

最后，我想避免使用阶乘，因为我将要计算的组合数可能会变得太大，而阶乘将会变得非常可怕。

这似乎是一个非常容易回答的问题，然而我被淹没在关于生成所有实际组合的问题中，这不是我想要的。

138042 次浏览
小开
最佳答案

请参见 特别的梳子(scypy.misc.comb 的旧版本)。当 exact为 False 时，采用 Gammaln 函数，不需要花费很多时间就可以获得很好的精度。在精确的情况下，它返回一个任意精度的整数，这可能需要很长的时间来计算。

小开

在谷歌代码上快速搜索一下(它使用的公式来自 @ Mark Byers 的回答) :

def choose(n, k):
"""
A fast way to calculate binomial coefficients by Andrew Dalke (contrib).
"""
if 0 <= k <= n:
ntok = 1
ktok = 1
for t in xrange(1, min(k, n - k) + 1):
ntok *= n
ktok *= t
n -= 1
return ntok // ktok
else:
return 0

如果需要精确的答案，choose()比 scipy.misc.comb()快10倍(在所有0 < = (n，k) < 1e3对上测试)。

def comb(N,k): # from scipy.comb(), but MODIFIED!
if (k > N) or (N < 0) or (k < 0):
return 0L
N,k = map(long,(N,k))
top = N
val = 1L
while (top > (N-k)):
val *= top
top -= 1
n = 1L
while (n < k+1L):
val /= n
n += 1
return val
小开

如果你想要精确的结果 还有的速度，尝试 Gmpy-gmpy.comb应该完全做你所要求的，还有它是相当快的(当然，作为 gmpy的原始作者，我 我的偏见;)。

小开

为什么不自己写呢? 这只是一句俏皮话:

from operator import mul    # or mul=lambda x,y:x*y
from fractions import Fraction


def nCk(n,k):
return int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )

测试打印帕斯卡三角形:

>>> for n in range(17):
...     print ' '.join('%5d'%nCk(n,k) for k in range(n+1)).center(100)
...
1
1     1
1     2     1
1     3     3     1
1     4     6     4     1
1     5    10    10     5     1
1     6    15    20    15     6     1
1     7    21    35    35    21     7     1
1     8    28    56    70    56    28     8     1
1     9    36    84   126   126    84    36     9     1
1    10    45   120   210   252   210   120    45    10     1
1    11    55   165   330   462   462   330   165    55    11     1
1    12    66   220   495   792   924   792   495   220    66    12     1
1    13    78   286   715  1287  1716  1716  1287   715   286    78    13     1
1    14    91   364  1001  2002  3003  3432  3003  2002  1001   364    91    14     1
1    15   105   455  1365  3003  5005  6435  6435  5005  3003  1365   455   105    15     1
1    16   120   560  1820  4368  8008 11440 12870 11440  8008  4368  1820   560   120    16     1
>>>

编辑以取代 int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) 与 int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))，所以它不会错误的大 N/K

小开

还有一个办法。这个函数最初是用 C + + 编写的，所以对于有限精度的整数(例如 _ _ int64) ，它可以向后移植到 C + + 。优点是(1)它只涉及整数运算，(2)它通过连续进行乘除运算避免了整数值的膨胀。我用 Nas Banov 的 Pascal 三角形测试了结果，它得到了正确的答案:

def choose(n,r):
"""Computes n! / (r! (n-r)!) exactly. Returns a python long int."""
assert n >= 0
assert 0 <= r <= n


c = 1L
denom = 1
for (num,denom) in zip(xrange(n,n-r,-1), xrange(1,r+1,1)):
c = (c * num) // denom
return c

基本原理: 为了最小化乘法和除法的 # ，我们将表达式重写为

    n!      n(n-1)...(n-r+1)
--------- = ----------------
r!(n-r)!          r!

为了尽可能避免乘法溢出，我们将按照以下 STRICT 顺序从左到右进行评估:

n / 1 * (n-1) / 2 * (n-2) / 3 * ... * (n-r+1) / r

我们可以证明按这个顺序运算的整数算术是精确的(即没有舍入误差)。

小开

在很多情况下，数学定义的直译是足够的(记住 Python 会自动使用大数算术) :

from math import factorial


def calculate_combinations(n, r):
return factorial(n) // factorial(r) // factorial(n-r)

对于我测试的一些输入(例如 n = 1000r = 500) ，这比另一个(目前投票最多的)答案中建议的一行 reduce快10倍以上。另一方面，它比@J 提供的片段执行得更好。F. Sebastian.

小开

如果你想要一个精确的结果，使用 sympy.binomial。这似乎是最快的方法，毫无疑问。

x = 1000000
y = 234050


%timeit scipy.misc.comb(x, y, exact=True)
1 loops, best of 3: 1min 27s per loop


%timeit gmpy.comb(x, y)
1 loops, best of 3: 1.97 s per loop


%timeit int(sympy.binomial(x, y))
100000 loops, best of 3: 5.06 µs per loop
小开

使用动态编程，时间复杂度为 Θ (n * m) ，空间复杂度为 Θ (m) :

def binomial(n, k):
""" (int, int) -> int


| c(n-1, k-1) + c(n-1, k), if 0 < k < n
c(n,k) = | 1                      , if n = k
| 1                      , if k = 0


Precondition: n > k


>>> binomial(9, 2)
36
"""


c = [0] * (n + 1)
c[0] = 1
for i in range(1, n + 1):
c[i] = 1
j = i - 1
while j > 0:
c[j] += c[j - 1]
j -= 1


return c[k]
小开

对于相当大的输入，这可能是在纯 python 中所能达到的最快速度:

def choose(n, k):
if k == n: return 1
if k > n: return 0
d, q = max(k, n-k), min(k, n-k)
num =  1
for n in xrange(d+1, n+1): num *= n
denom = 1
for d in xrange(1, q+1): denom *= d
return num / denom
小开

当 n 大于20时，直接公式产生大整数。

那么，另一个回答是:

from math import factorial


reduce(long.__mul__, range(n-r+1, n+1), 1L) // factorial(r)

短，准确和有效，因为这避免了巨蟒大整数坚持长。

相比之下，scipy.specal.comb 更准确、更快捷:

 >>> from scipy.special import comb
>>> nCr = lambda n,r: reduce(long.__mul__, range(n-r+1, n+1), 1L) // factorial(r)
>>> comb(128,20)
1.1965669823265365e+23
>>> nCr(128,20)
119656698232656998274400L  # accurate, no loss
>>> from timeit import timeit
>>> timeit(lambda: comb(n,r))
8.231969118118286
>>> timeit(lambda: nCr(128, 20))
3.885951042175293
小开

对付 Symy 很简单。

import sympy


comb = sympy.binomial(n, r)
小开

仅使用 使用 Python 发布的标准库:

import itertools


def nCk(n, k):
return len(list(itertools.combinations(range(n), k)))
小开

如果程序的上限是 n(比如说 n <= N) ，并且需要重复计算 nCr (最好是 > > N次) ，那么使用 Lru _ cache可以极大地提高性能:

from functools import lru_cache


@lru_cache(maxsize=None)
def nCr(n, r):
return 1 if r == 0 or r == n else nCr(n - 1, r - 1) + nCr(n - 1, r)

构建缓存(隐式地完成)需要花费大约 O(N^2)时间。对 nCr的任何后续调用将在 O(1)中返回。

小开

您可以编写两个简单的函数，它们实际上比使用 特别的梳子快5-8倍。实际上，您不需要导入任何额外的包，而且该函数非常容易阅读。诀窍是使用制表来存储以前计算的值，并使用 NCr的定义

# create a memoization dictionary
memo = {}
def factorial(n):
"""
Calculate the factorial of an input using memoization
:param n: int
:rtype value: int
"""
if n in [1,0]:
return 1
if n in memo:
return memo[n]
value = n*factorial(n-1)
memo[n] = value
return value


def ncr(n, k):
"""
Choose k elements from a set of n elements - n must be larger than or equal to k
:param n: int
:param k: int
:rtype: int
"""
return factorial(n)/(factorial(k)*factorial(n-k))

如果我们比较时间

from scipy.special import comb
%timeit comb(100,48)
>>> 100000 loops, best of 3: 6.78 µs per loop


%timeit ncr(100,48)
>>> 1000000 loops, best of 3: 1.39 µs per loop
小开

这是@kilerT2333代码，使用内置的制表修饰符。

from functools import lru_cache


@lru_cache()
def factorial(n):
"""
Calculate the factorial of an input using memoization
:param n: int
:rtype value: int
"""
return 1 if n in (1, 0) else n * factorial(n-1)


@lru_cache()
def ncr(n, k):
"""
Choose k elements from a set of n elements,
n must be greater than or equal to k.
:param n: int
:param k: int
:rtype: int
"""
return factorial(n) / (factorial(k) * factorial(n - k))


print(ncr(6, 3))
小开

这个函数是非常优化的。

def nCk(n,k):
m=0
if k==0:
m=1
if k==1:
m=n
if k>=2:
num,dem,op1,op2=1,1,k,n
while(op1>=1):
num*=op2
dem*=op1
op1-=1
op2-=1
m=num//dem
return m
小开

从 Python 3.8开始，标准库现在包括了计算二项式系数的 math.comb函数:

Comb (n，k)

这是从 n 个项目中选择 k 个项目而不重复的方法的数量 < br > n! / (k! (n - k)!):

import math
math.comb(10, 5) # 252
小开

这里有一个有效的算法

for i = 1.....r


p = p * ( n - i ) / i


print(p)

例如 nCr (30,7) = 事实(30)/(事实(7) * 事实(23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24)/(1 * 2 * 3 * 4 * 5 * 6 * 7)

所以只要运行从1到 r 的循环就可以得到结果。

在巨蟒中:

n,r=5,2
p=n
for i in range(1,r):
p = p*(n - i)/i
else:
p = p/(i+1)
print(p)
小开

我从这个线程和链接到这里的库中计时了17个不同的函数。

因为我觉得在这里转储有点多，所以我把函数的代码放在这里的粘贴文件中。

我做的第一个测试是把帕斯卡三角形建到第100行。我用时间做了100次。下面的数字是构建一次三角形所花费的平均时间(秒)。

gmpy2.gmpy2.comb 0.0012259269999998423
math.comb 0.007063110999999935
__main__.stdfactorial2 0.011469491
__main__.scipybinom 0.0120114319999999
__main__.stdfactorial 0.012105122
__main__.scipycombexact 0.012569045999999844
__main__.andrewdalke 0.01825201100000015
__main__.rabih 0.018472497000000202
__main__.kta 0.019374668000000383
__main__.wirawan 0.029312811000000067
scipy.special._basic.comb 0.03221609299999954
__main__.jfsmodifiedscipy 0.04332894699999997
__main__.rojas 0.04395155400000021
sympy.functions.combinatorial.factorials.binomial 0.3233529779999998
__main__.nasbanov 0.593365528
__main__.pantelis300 1.7780402499999999

您可能注意到这里只有16个函数。这是因为 recursive()函数甚至不能在合理的时间内完成一次，所以我不得不从它的测试时间中排除它。说真的，已经好几个小时了。

我还计算了并非上述所有函数都支持的各种其他类型的输入的时间。请记住，我只对每个测试运行了10次，因为 nCr 的计算开销很大，而且我没有耐心

N 的分数值

__main__.scipybinom 0.011481370000000001
__main__.kta 0.01869513999999999
sympy.functions.combinatorial.factorials.binomial 6.33897291

R 的分数值

__main__.scipybinom 0.010960040000000504
scipy.special._basic.comb 0.03681254999999908
sympy.functions.combinatorial.factorials.binomial 3.2962564499999987

N 和 r 的分数值

__main__.scipybinom 0.008623409999998444
sympy.functions.combinatorial.factorials.binomial 3.690936439999999

N 的负值

gmpy2.gmpy2.comb 0.010770989999997482
__main__.kta 0.02187850000000253
__main__.rojas 0.05104292999999984
__main__.nasbanov 0.6153183200000001
sympy.functions.combinatorial.factorials.binomial 3.0460310799999943

N 的负分数值，r 的分数值

sympy.functions.combinatorial.factorials.binomial 3.7689941699999965

目前实现最大速度和通用性的最佳解决方案是一个混合函数，根据输入在不同算法之间进行选择

def hybrid(n: typing.Union[int, float], k: typing.Union[int, float]) -> typing.Union[int, float]:
# my own custom hybrid solution
def is_integer(n):
return isinstance(n, int) or n.is_integer()
if k < 0:
raise ValueError("k cannot be negative.")
elif n == 0:
return 0
elif k == 0 or k == n:
return 1
elif is_integer(n) and is_integer(k):
return int(gmpy2.comb(int(n), int(k)))
elif n > 0:
return scipy.special.binom(n, k)
else:
return float(sympy.binomial(n, k))

由于 sympy.binomial()是如此之慢，真正理想的解决方案将是结合的代码的 scipy.special.binom()表现良好的分数和 gmpy2.comb()表现良好的整数。Scipy 的功能和 健身房2的功能都是用我不太熟悉的 C 语言写的。


提交答案