首页 文章

统计:Python中的组合

提问于
浏览
97

我需要在Python中计算组合(nCr),但无法在 mathnumpystat 库中找到执行此操作的函数 . 类似于类型函数的东西:

comb = calculate_combinations(n, r)

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

最后,我想避免使用阶乘,因为我将计算组合的数字会变得太大而且阶乘将变得非常可怕 .

这似乎是一个非常容易回答的问题,但是我被淹没在关于生成所有实际组合的问题中,这不是我想要的 . :)

非常感谢

14 回答

  • 45

    快速搜索谷歌代码给出(它使用@Mark Byers's answer中的公式):

    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
    
  • 19

    请参阅scipy.special.comb(旧版scipy中的scipy.misc.comb) . 当 exact 为False时,它使用gammaln函数获得良好的精度而不需要花费太多时间 . 在确切的情况下,它返回一个任意精度的整数,这可能需要很长时间才能计算出来 .

  • 1

    如果您的程序的上限为 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) 中返回 .

  • 23

    为什么不自己写呢?这是一个单行或类似的:

    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) )
    

    测试 - 打印Pascal的三角形:

    >>> 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
    >>>
    

    PS . 编辑后将 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错误

  • 1

    使用动态编程,时间复杂度为Θ(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]
    
  • 103

    用同情的话很容易 .

    import sympy
    
    comb = sympy.binomial(n, r)
    
  • 96

    这是另一种选择 . 这个最初是用C语言编写的,因此它可以被后向移植到C以获得有限精度的整数(例如__int64) . 优点是(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
    

    我们可以证明按此顺序操作的整数算术是精确的(即没有舍入误差) .

  • 4

    在很多情况下,数学定义的字面翻译是相当充分的(记住Python会自动使用大数字算术):

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

    对于我测试的一些输入(例如n = 1000 r = 500),这比另一个(当前最高投票)答案中提出的一个衬管 reduce 快10倍以上 . 另一方面,它由@ J.F提供的snippit表示 . 塞巴斯蒂安 .

  • 2

    仅使用随Python分发的标准库:

    import itertools
    
    def nCk(n, k):
        return len(list(itertools.combinations(range(n), k)))
    
  • 2

    对于相当大的输入,这可能与纯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
    
  • 9

    您可以编写2个简单的函数,实际上比使用scipy.special.comb快5-8倍 . 实际上,您不需要导入任何额外的包,并且该功能非常容易阅读 . 诀窍是使用memoization存储以前计算的值,并使用nCr的定义

    # create a memoization dict
    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*fact(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
    
  • 0

    当n大于20时,直接公式产生大整数 .

    那么,还有另一个回应:

    from math import factorial
    
    binomial = lambda n,r: reduce(long.__mul__, range(n-r, n+1), 1L) // factorial(r)
    

    简短,快速,高效 .

  • 0

    如果你想要确切的结果 and 速度,请尝试gmpy - gmpy.comb 应该完全按照你的要求做,而且速度非常快(当然,作为 gmpy 的原作者,我有偏见;-) .

  • 40

    如果您想要精确的结果,请使用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
    

相关问题