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

Morlok
źródło

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.

Jouni K. Seppänen
źródło
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

Nas Banov
źródło
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
jfs
źródło
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 ;-).

Alex Martelli
źródło
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
Jim Garrison
źródło
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.

Todd Owen
źródło
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
Xavier Guihot
źródło
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).

Wirawan Purwanto
źródło
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]
pantelis300
źródło
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).

yzn-pku
źródło
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
PyRsquared
źródło
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)
Konstabl
źródło
2

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

import itertools

def nCk(n, k):
    return len(list(itertools.combinations(range(n), k)))
MarianD
źródło
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
olivecoder
źródło
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))
obłąkany jeż
źródło
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.

kta
źródło
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
Rabih Kodeih
źródło
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
Santiago Coca Rojas
źródło