Statistiques: combinaisons en Python

Je dois calculer combinatoires (nCr) en Python mais ne peut pas trouver la fonction de le faire dans math, numpyou les stat bibliothèques. Quelque chose comme une fonction du type:

comb = calculate_combinations(n, r)

J'ai besoin du nombre de combinaisons possibles, pas des combinaisons réelles, donc itertools.combinationscela ne m'intéresse pas.

Enfin, je veux éviter d'utiliser des factorielles, car les nombres pour lesquels je vais calculer les combinaisons peuvent devenir trop gros et les factorielles vont être monstrueuses.

Cela semble être une question VRAIMENT facile à répondre, mais je suis noyé dans des questions sur la génération de toutes les combinaisons réelles, ce qui n'est pas ce que je veux.

Voir scipy.special.comb (scipy.misc.comb dans les anciennes versions de scipy). Quand exactest False, il utilise la fonction gammaln pour obtenir une bonne précision sans prendre beaucoup de temps. Dans le cas exact, il renvoie un entier de précision arbitraire, dont le calcul peut prendre beaucoup de temps.

Pourquoi ne pas l'écrire vous-même? C'est un one-liner ou tel:

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

Test - impression du triangle 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. modifié pour remplacer int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) par int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))afin qu'il ne se trompe pas pour le grand N / K

Une recherche rapide sur le code google donne (elle utilise la formule de la réponse 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()est 10 fois plus rapide (testé sur toutes les paires 0 <= (n, k) <1e3) que scipy.misc.comb()si vous avez besoin d'une réponse exacte.

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

Si vous voulez des résultats et une vitesse exacts , essayez gmpy - gmpy.combdevrait faire exactement ce que vous demandez, et c'est assez rapide (bien sûr, en tant gmpyqu'auteur original, je suis partial ;-).

Si vous voulez un résultat exact, utilisez sympy.binomial. Cela semble être la méthode la plus rapide, haut la main.

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

Une traduction littérale de la définition mathématique est tout à fait adéquate dans de nombreux cas (en se rappelant que Python utilisera automatiquement l'arithmétique des grands nombres):

from math import factorial

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

Pour certaines entrées que j'ai testées (par exemple n = 1000 r = 500), c'était plus de 10 fois plus rapide que la ligne reducesuggérée dans une autre réponse (actuellement la plus votée). En revanche, il est surpassé par l'extrait de code fourni par @JF Sebastian.

Au départ Python 3.8, la bibliothèque standard comprend désormais la math.combfonction de calcul du coefficient binomial:

math.comb (n, k)

qui est le nombre de façons de choisir k éléments parmi n éléments sans répétition
n! / (k! (n - k)!):

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

Voici une autre alternative. Celui-ci a été écrit à l'origine en C ++, il peut donc être rétroporté en C ++ pour un entier de précision finie (par exemple __int64). L'avantage est (1) qu'il n'implique que des opérations entières, et (2) il évite de gonfler la valeur entière en faisant des paires successives de multiplication et de division. J'ai testé le résultat avec le triangle Pascal de Nas Banov, il obtient la bonne réponse:

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

Justification: pour minimiser le nombre de multiplications et de divisions, nous réécrivons l'expression comme

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

Pour éviter autant que possible le débordement de multiplication, nous évaluerons dans l'ordre STRICT suivant, de gauche à droite:

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

Nous pouvons montrer que l'arithmatique entière opérée dans cet ordre est exacte (c'est-à-dire pas d'erreur d'arrondi).

En utilisant la programmation dynamique, la complexité temporelle est Θ (n * m) et la complexité spatiale Θ (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 votre programme a une limite supérieure à n(par exemple n <= N) et doit calculer à plusieurs reprises nCr (de préférence pour >> Nfois), l'utilisation de lru_cache peut vous donner un énorme gain de performances:

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 construction du cache (ce qui est fait implicitement) prend du O(N^2)temps. Tous les appels ultérieurs à nCrretourneront O(1).

Vous pouvez écrire 2 fonctions simples qui s'avèrent être environ 5 à 8 fois plus rapides qu'en utilisant scipy.special.comb . En fait, vous n'avez pas besoin d'importer de packages supplémentaires et la fonction est assez facilement lisible. L'astuce consiste à utiliser la mémorisation pour stocker les valeurs précédemment calculées et à utiliser la définition 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))

Si nous comparons les temps

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

C'est assez facile avec sympy.

import sympy

comb = sympy.binomial(n, r)

En utilisant uniquement la bibliothèque standard distribuée avec Python :

import itertools

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

La formule directe produit de grands entiers lorsque n est supérieur à 20.

Alors, encore une autre réponse:

from math import factorial

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

court, précis et efficace car cela évite les grands entiers python en s'en tenant aux longs.

Il est plus précis et plus rapide par rapport à 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

Il s'agit du code @ killerT2333 utilisant le décorateur de mémorisation intégré.

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

Voici un algorithme efficace pour vous

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

Par exemple nCr (30,7) = fact (30) / (fact (7) * fact (23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24) / (1 * 2 * 3 * 4 * 5 * 6 * 7)

Donc, exécutez simplement la boucle de 1 à r pour obtenir le résultat.

C'est probablement aussi rapide que vous pouvez le faire en python pur pour des entrées raisonnablement volumineuses:

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

Cette fonction est très optimisée.

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