Estatísticas: combinações em Python

Eu preciso calcular combinatorials (nCr) em Python, mas não consegue encontrar a função de fazer isso em math, numpyou stat bibliotecas. Algo como uma função do tipo:

comb = calculate_combinations(n, r)

Preciso do número de combinações possíveis, não das combinações reais, por itertools.combinationsisso não me interessa.

Por fim, quero evitar o uso de fatoriais, pois os números para os quais vou calcular as combinações podem ficar muito grandes e os fatoriais serão monstruosos.

Parece uma pergunta REALMENTE fácil de responder, no entanto, estou sendo afogado em perguntas sobre como gerar todas as combinações reais, o que não é o que eu quero.

Consulte scipy.special.comb (scipy.misc.comb em versões mais antigas do scipy). Quando exactFalse, ele usa a função gammaln para obter boa precisão sem levar muito tempo. No caso exato, ele retorna um número inteiro de precisão arbitrária, que pode levar muito tempo para ser computado.

Por que não escrever você mesmo? É uma linha ou algo assim:

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

Teste - impressão do triângulo de 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. editado para substituir int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) com int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))isso, não será err para grande N / K

Uma rápida pesquisa no código do google fornece (ele usa a fórmula da resposta de @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()é 10 vezes mais rápido (testado em todos os pares 0 <= (n, k) <1e3) do que scipy.misc.comb()se você precisar de uma resposta exata.

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

Se você deseja resultados e velocidade exatos , tente o gmpy - gmpy.combfaça exatamente o que você pede e é muito rápido (é claro, como gmpyautor original do site, sou tendencioso ;-).

Se você deseja um resultado exato, use sympy.binomial. Parece ser o método mais rápido, sem dúvida.

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

Uma tradução literal da definição matemática é bastante adequada em muitos casos (lembrando que o Python usará automaticamente a aritmética de grandes números):

from math import factorial

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

Para algumas entradas que testei (por exemplo, n = 1000 r = 500), isso foi mais de 10 vezes mais rápido do que o liner reducesugerido em outra resposta (atualmente com o maior voto). Por outro lado, é superado pelo snippit fornecido por @JF Sebastian.

Começando Python 3.8, a biblioteca padrão agora inclui a math.combfunção para calcular o coeficiente binomial:

math.comb (n, k)

qual é o número de maneiras de escolher k itens de n itens sem repetição
n! / (k! (n - k)!):

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

Aqui está outra alternativa. Este foi originalmente escrito em C ++, para que possa ser portado em C ++ para um número inteiro de precisão finita (por exemplo, __int64). A vantagem é (1) envolver apenas operações com números inteiros e (2) evitar inchar o valor inteiro, fazendo pares sucessivos de multiplicação e divisão. Testei o resultado com o triângulo Pascal de Nas Banov, ele obtém a resposta correta:

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

Fundamentação da petição: Para minimizar o número de multiplicações e divisões, reescrevemos a expressão como

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

Para evitar o excesso de multiplicação, tanto quanto possível, avaliaremos na seguinte ordem STRICT, da esquerda para a direita:

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

Podemos mostrar que a aritmética inteira operada nesta ordem é exata (ou seja, nenhum erro de arredondamento).

Usando programação dinâmica, a complexidade do tempo é Θ (n * m) e a complexidade do espaço Θ (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]

Se o seu programa tiver um limite superior para n(digamos n <= N) e precisar calcular repetidamente a nCr (de preferência por >> Nvezes), o uso do lru_cache poderá oferecer um enorme aumento de desempenho:

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)

Construir o cache (que é feito implicitamente) leva O(N^2)tempo. Quaisquer chamadas subseqüentes nCrretornarão O(1).

Você pode escrever duas funções simples que, na verdade, são cerca de 5 a 8 vezes mais rápidas do que usar scipy.special.comb . De fato, você não precisa importar nenhum pacote extra, e a função é facilmente legível. O truque é usar a memorização para armazenar valores previamente calculados e usar a definição de 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))

Se compararmos os tempos

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

É bem fácil com o sympy.

import sympy

comb = sympy.binomial(n, r)

Usando apenas biblioteca padrão distribuída com Python :

import itertools

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

A fórmula direta produz grandes números inteiros quando n é maior que 20.

Então, mais uma resposta:

from math import factorial

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

curto, preciso e eficiente, porque isso evita inteiros grandes em python, permanecendo com longs.

É mais preciso e mais rápido quando comparado ao scipy.special.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

Este é o código @ killerT2333 usando o decorador de memorização incorporado.

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

Aqui está um algoritmo eficiente para você

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

Por exemplo, nCr (30,7) = fato (30) / (fato (7) * fato (23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24) / (1 * 2 * 3 * 4 * 5 * 6 * 7)

Portanto, basta executar o loop de 1 a r para obter o resultado.

Provavelmente é o mais rápido que você pode fazer em python puro para entradas razoavelmente grandes:

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

Esta função é muito otimizada.

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