Estadísticas: combinaciones en Python

Necesito calcular combinatorials (nCr) en Python pero no puedo encontrar la función para hacer eso en las bibliotecas de math , numpy o stat . Algo así como una función del tipo:

 comb = calculate_combinations(n, r) 

Necesito el número de combinaciones posibles, no las combinaciones reales, por lo que itertools.combinations no me interesa.

Finalmente, quiero evitar el uso de factoriales, ya que los números para los que calcularé las combinaciones pueden ser demasiado grandes y los factoriales serán monstruosos.

Esto parece una pregunta REALMENTE fácil de responder, sin embargo, me están ahogando las preguntas sobre cómo generar todas las combinaciones reales, que no es lo que quiero. 🙂

Muchas gracias

Consulte scipy.special.comb (scipy.misc.comb en versiones anteriores de scipy). Cuando exact es Falso, utiliza la función de juego para obtener una buena precisión sin tomar mucho tiempo. En el caso exacto, devuelve un entero de precisión arbitraria, que puede tardar mucho tiempo en calcularse.

¿Por qué no lo escribes tú mismo? Es una sola línea o tal:

 from operator import mul # or mul=lambda x,y:x*y from fractions import Fraction def nCk(n,k): return int( reduce(mul, (Fraction(ni, i+1) for i in range(k)), 1) ) 

Prueba – imprimiendo el 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 >>> 

PD. editado para reemplazar int(round(reduce(mul, (float(ni)/(i+1) for i in range(k)), 1))) con int(reduce(mul, (Fraction(ni, i+1) for i in range(k)), 1)) por lo que no errará para N / K grande

Una búsqueda rápida en el código de Google da (utiliza la fórmula de la respuesta 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() es 10 veces más rápido (probado en todos los 0 <= (n, k) <1e3 pares) que scipy.misc.comb() si necesita una respuesta exacta.

 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 > (Nk)): val *= top top -= 1 n = 1L while (n < k+1L): val /= n n += 1 return val 

Si quieres resultados y velocidad exactos, prueba gmpy – gmpy.comb debería hacer exactamente lo que pides, y es bastante rápido (por supuesto, como autor original de gmpy , tengo prejuicios ;-).

Si quieres un resultado exacto, usa sympy.binomial . Parece ser el método más rápido, sin duda.

 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 

Una traducción literal de la definición matemática es bastante adecuada en muchos casos (recordando que Python usará automáticamente la aritmética de números grandes):

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

Para algunas entradas que probé (por ejemplo, n = 1000 r = 500) esto fue más de 10 veces más rápido que la reduce un forro sugerida en otra respuesta (actualmente votada con el mayor número de votos). Por otro lado, es superado por el fragmento proporcionado por @JF Sebastian.

Aquí hay otra alternativa. Este se escribió originalmente en C ++, por lo que se puede volver a C ++ para un entero de precisión finita (por ejemplo, __int64). La ventaja es (1) que involucra solo operaciones de enteros, y (2) evita la hinchazón del valor entero al hacer pares sucesivos de multiplicación y división. He probado el resultado con el triángulo de Pascal de Nas Banov, se obtiene la respuesta correcta:

 def choose(n,r): """Computes n! / (r! (nr)!) exactly. Returns a python long int.""" assert n >= 0 assert 0 <= r <= n c = 1L denom = 1 for (num,denom) in zip(xrange(n,nr,-1), xrange(1,r+1,1)): c = (c * num) // denom return c 

Justificación: para minimizar el # de multiplicaciones y divisiones, reescribimos la expresión como

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

Para evitar el desbordamiento de la multiplicación tanto como sea posible, evaluaremos en el siguiente orden ESTRICTO, de izquierda a derecha:

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

Podemos mostrar que la aritmática de enteros operada en este orden es exacta (es decir, sin error de redondeo).

Usando la progtwigción dinámica, la complejidad del tiempo es Θ (n * m) y la complejidad del espacio Θ (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] 

Si su progtwig tiene un límite superior a n (digamos n <= N ) y necesita computar repetidamente nCr (preferiblemente por >> N veces), usar lru_cache puede darle un gran impulso de rendimiento:

 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) 

La construcción de la memoria caché (que se realiza de manera implícita) lleva hasta O(N^2) tiempo. Cualquier llamada subsiguiente a nCr regresará en O(1) .

Puede escribir 2 funciones simples que en realidad resultan ser 5 a 8 veces más rápidas que usar scipy.special.comb . De hecho, no es necesario importar ningún paquete adicional, y la función es bastante fácil de leer. El truco es usar la memoria para almacenar valores previamente calculados y usar la definición de 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(nk)) 

Si comparamos tiempos

 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 

La fórmula directa produce grandes enteros cuando n es mayor que 20.

Entonces, otra respuesta más:

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

Corto, rápido y eficiente.

Es bastante fácil con sympy.

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

Este es el código @ killerT2333 que utiliza el decorador de memorización 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)) 

Probablemente sea lo más rápido que pueda hacerlo en python puro para entradas razonablemente grandes:

 def choose(n, k): if k == n: return 1 if k > n: return 0 d, q = max(k, nk), min(k, nk) 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 

Usando solo la biblioteca estándar distribuida con Python :

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