统计学：Python中的组合

Question

我需要在Python中计算组合数（nCr），但是在math、numpy或stat库中找不到相关的函数。我想要的函数大概是这样的：

comb = calculate_combinations(n, r)

我只需要可能的组合数量，而不是实际的组合，所以itertools.combinations对我来说没什么用。

最后，我想避免使用阶乘，因为我计算的组合数可能会很大，阶乘的数值会变得非常庞大。

这个问题看起来应该很简单，但我却被关于生成所有实际组合的问题淹没了，而这并不是我想要的。

Answer 1

在谷歌代码上快速搜索一下，得到了一个结果（这个结果使用了@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

Answer 2

为什么不自己写呢？
这其实只需要一行代码或者类似的东西：

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时就不会出错了。

Answer 3

2023年的更新回答：可以使用 math.comb 这个函数，它从Python 3.8开始就有了，并且在3.11版本中变得 更快了很多。

---

旧的回答：可以查看 scipy.special.comb（在旧版本的scipy中是scipy.misc.comb）。当 exact 设置为False时，它会使用gammaln函数来获得较好的精度，同时不会花太多时间。如果设置为精确计算，它会返回一个任意精度的整数，这可能需要较长的计算时间。