Statystyki: kombinacje w Pythonie


122

Muszę obliczyć combinatorials (NCR) w Pythonie, ale nie może znaleźć funkcji do zrobienia, że math, numpyczy stat bibliotek. Coś w rodzaju funkcji typu:

comb = calculate_combinations(n, r)

Potrzebuję liczby możliwych kombinacji, a nie rzeczywistych kombinacji, więc itertools.combinationsmnie to nie interesuje.

Na koniec chcę uniknąć silni, ponieważ liczby, dla których będę obliczać kombinacje, mogą być zbyt duże, a silnie będą potworne.

Wydaje się, że odpowiedź na to pytanie jest NAPRAWDĘ łatwa, jednak tonę w pytaniach o generowanie wszystkich rzeczywistych kombinacji, czego nie chcę.

Odpowiedzi:


122

Zobacz scipy.special.comb (scipy.misc.comb w starszych wersjach scipy). Gdy exactjest fałszywe, używa funkcji gammaln, aby uzyskać dobrą precyzję bez zajmowania dużo czasu. W dokładnym przypadku zwraca liczbę całkowitą o dowolnej precyzji, której obliczenie może zająć dużo czasu.


5
scipy.misc.combjest przestarzałe na korzyść scipy.special.combod wersji 0.10.0.
Dilawar

120

Dlaczego nie napisać tego samemu? To jeden wiersz lub coś takiego:

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 - drukowanie trójkąta Pascala:

>>> 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. edytowane w celu zastąpienia int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) przez, int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))więc nie będzie błądzić dla dużych N / K


26
+1 za sugestię napisania czegoś prostego, użycie
reduktora

6
-1, ponieważ ta odpowiedź jest błędna: wypisz silnię (54) / (silnię (54 - 27)) / silnię (27) == nCk (54, 27) daje Fałsz.
robert king

3
@robertking - Ok, byłeś drobny i technicznie poprawny. To, co zrobiłem, miało być ilustracją tego, jak napisać własną funkcję; Wiedziałem, że nie jest to dokładne dla wystarczająco dużych N i K ze względu na precyzję zmiennoprzecinkową. Ale możemy to naprawić - patrz powyżej, teraz nie powinno się mylić przy dużych liczbach
Nas Banov

9
Prawdopodobnie byłoby to szybkie w Haskell, ale niestety nie w Pythonie. W rzeczywistości jest to dość powolne w porównaniu z wieloma innymi odpowiedziami, np. @Alex Martelli, JF Sebastian i moją własną.
Todd Owen

9
W przypadku Pythona 3 też musiałem from functools import reduce.
Velizar Hristov

52

Szybkie wyszukiwanie w kodzie google daje (wykorzystuje formułę z odpowiedzi @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()jest 10 razy szybszy (testowany na wszystkich parach 0 <= (n, k) <1e3), niż scipy.misc.comb()gdybyś potrzebował dokładnej odpowiedzi.

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

Fajne rozwiązanie, które nie wymaga żadnego pakietu
Edward Newell

2
Do Twojej wiadomości: wspomniana formuła jest tutaj: en.wikipedia.org/wiki/ ...
jmiserez

Ta choosefunkcja powinna mieć znacznie więcej pozytywnych głosów! Python 3.8 ma math.comb, ale musiałem użyć Pythona 3.6 do wyzwania i żadna implementacja nie dała dokładnych wyników dla bardzo dużych liczb całkowitych. Ten robi i robi to szybko!
połącz ponownie

42

Jeśli chcesz uzyskać dokładne wyniki i szybkość, wypróbuj gmpy - gmpy.combpowinien robić dokładnie to, o co prosisz, i jest dość szybki (oczywiście jako gmpyautor jestem stronniczy ;-).


6
Rzeczywiście, gmpy2.comb()jest 10 razy szybszy niż choose()z mojej odpowiedzi dla kodu: for k, n in itertools.combinations(range(1000), 2): f(n,k)gdzie f()jest albo gmpy2.comb()albo choose()na Pythonie 3.
jfs

Ponieważ jesteś autorem pakietu, pozwolę Ci naprawić zepsuty link, aby wskazywał we właściwym miejscu ....
SeldomNeedy,

@SeldomNeedy, link do code.google.com to jedno właściwe miejsce (chociaż strona jest teraz w trybie archiwalnym). Oczywiście stamtąd łatwo jest znaleźć lokalizację github, github.com/aleaxit/gmpy i PyPI, pypi.python.org/pypi/gmpy2 , ponieważ prowadzi do obu! -)
Alex Martelli

@AlexMartelli Przepraszamy za zamieszanie. Strona wyświetla 404, jeśli javascript został (selektywnie) wyłączony. Myślę, że ma to zniechęcić nieuczciwe sztucznej inteligencji do tak łatwego włączania zarchiwizowanych źródeł projektu Google Code?
RzadkoNeedy

28

Jeśli chcesz uzyskać dokładny wynik, użyj sympy.binomial. Wydaje się, że jest to najszybsza metoda.

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

22

Dosłowne tłumaczenie definicji matematycznej jest wystarczające w wielu przypadkach (pamiętając, że Python automatycznie użyje arytmetyki dużych liczb):

from math import factorial

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

Dla niektórych testowanych przeze mnie danych wejściowych (np. N = 1000 r = 500) było to ponad 10 razy szybsze niż jedna linijka reducesugerowana w innej odpowiedzi (aktualnie najwyżej głosowanej). Z drugiej strony wyprzedza go fragment dostarczony przez @JF Sebastian.


11

Zaczynając Python 3.8, biblioteka standardowa zawiera teraz math.combfunkcję obliczania współczynnika dwumianu:

math.comb (n, k)

czyli liczba sposobów wyboru k elementów z n elementów bez powtórzeń
n! / (k! (n - k)!):

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

10

Oto inna alternatywa. Ten został pierwotnie napisany w C ++, więc można go przenieść do C ++ w celu uzyskania liczby całkowitej o skończonej precyzji (np. __Int64). Zaletą jest to, że (1) obejmuje tylko operacje na liczbach całkowitych, a (2) pozwala uniknąć powiększania wartości całkowitej poprzez wykonywanie kolejnych par mnożenia i dzielenia. Przetestowałem wynik za pomocą trójkąta Pascala Nas Banova, otrzymuje poprawną odpowiedź:

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

Uzasadnienie: aby zminimalizować liczbę mnożeń i dzieleń, przepisujemy wyrażenie jako

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

Aby uniknąć przepełnienia mnożenia w jak największym stopniu, będziemy oceniać w następującej kolejności STRICT, od lewej do prawej:

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

Możemy pokazać, że arytmatyka liczb całkowitych wykonywana w tej kolejności jest dokładna (tj. Nie ma błędu zaokrąglenia).


5

Używając programowania dynamicznego, złożoność czasowa wynosi Θ (n * m), a złożoność przestrzeni Θ (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]

4

Jeśli twój program ma górną granicę n(powiedzmy n <= N) i musi wielokrotnie obliczać nCr (najlepiej >> Nrazy), użycie lru_cache może dać ogromny wzrost wydajności:

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)

Konstruowanie pamięci podręcznej (co jest wykonywane niejawnie) zajmuje trochę O(N^2)czasu. Wszelkie kolejne połączenia z numerem nCrpowrócą za O(1).


4

Możesz napisać 2 proste funkcje, które w rzeczywistości okazują się być około 5-8 razy szybsze niż przy użyciu scipy.special.comb . W rzeczywistości nie musisz importować żadnych dodatkowych pakietów, a funkcja jest dość łatwa do odczytania. Sztuczka polega na tym, aby użyć zapamiętywania do przechowywania wcześniej obliczonych wartości i użyć definicji 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))

Jeśli porównamy czasy

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

W dzisiejszych czasach istnieje dekorator memoize w functools o nazwie lru_cache, który może uprościć twój kod?
obłąkany jeż

2

Z Sympy jest to całkiem proste.

import sympy

comb = sympy.binomial(n, r)

2

Używając tylko standardowej biblioteki dystrybuowanej z Pythonem :

import itertools

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

3
Nie sądzę, aby jego złożoność czasowa (i zużycie pamięci) była akceptowalna.
xmcp

2

Formuła bezpośrednia daje duże liczby całkowite, gdy n jest większe niż 20.

A więc kolejna odpowiedź:

from math import factorial

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

krótkie, dokładne i wydajne, ponieważ pozwala to uniknąć dużych liczb całkowitych w Pythonie poprzez trzymanie się długich liczb.

Jest dokładniejszy i szybszy w porównaniu do 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

To jest źle! Jeśli n == r, wynik powinien wynosić 1. Ten kod zwraca 0.
reyammer

Dokładniej, powinno być range(n-r+1, n+1)zamiast range(n-r,n+1).
reyammer

1

To jest kod @ killerT2333 wykorzystujący wbudowany dekorator zapamiętywania.

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

1

Oto skuteczny algorytm dla Ciebie

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

Na przykład nCr (30,7) = fakt (30) / (fakt (7) * fakt (23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24) / (1 * 2 * 3 * 4 * 5 * 6 * 7)

Więc wystarczy uruchomić pętlę od 1 do r, aby uzyskać wynik.


0

To prawdopodobnie tak szybko, jak możesz to zrobić w czystym Pythonie dla dość dużych danych wejściowych:

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

0

Ta funkcja jest bardzo zoptymalizowana.

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
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.