Statistik: Kombinationen in Python

Ich brauche combinatorials (nCr) in Python zu berechnen , aber die Funktion nicht finden können , das zu tun in math, numpyoder stat Bibliotheken. So etwas wie eine Funktion des Typs:

comb = calculate_combinations(n, r)

Ich brauche die Anzahl der möglichen Kombinationen, nicht die tatsächlichen Kombinationen, also itertools.combinationsinteressiert mich das nicht.

Schließlich möchte ich die Verwendung von Fakultäten vermeiden, da die Zahlen, für die ich die Kombinationen berechnen werde, zu groß werden können und die Fakultäten monströs sein werden.

Dies scheint eine WIRKLICH einfach zu beantwortende Frage zu sein, aber ich bin in Fragen über das Generieren aller tatsächlichen Kombinationen ertrunken, was ich nicht will.

Siehe scipy.special.comb (scipy.misc.comb in älteren Versionen von scipy). Wenn exactFalse ist, wird die Gammaln-Funktion verwendet, um eine gute Präzision zu erzielen, ohne viel Zeit in Anspruch zu nehmen. Im genauen Fall wird eine Ganzzahl mit beliebiger Genauigkeit zurückgegeben, deren Berechnung möglicherweise lange dauert.

Warum schreibst du es nicht selbst? Es ist ein Einzeiler oder so:

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 - Drucken des Pascalschen Dreiecks:

>>> 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. bearbeitet, um int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) durch zu ersetzen, int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))damit es nicht für große N / K irrt

Eine schnelle Suche nach Google-Code ergibt (es wird eine Formel aus der Antwort von @Mark Byers verwendet ):

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()ist 10 mal schneller (getestet an allen 0 <= (n, k) <1e3 Paaren) als scipy.misc.comb()wenn Sie eine genaue Antwort benötigen.

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

Wenn Sie genaue Ergebnisse und Geschwindigkeit wünschen , versuchen Sie es mit gmpy - gmpy.combsollte genau das tun, wonach Sie fragen, und es ist ziemlich schnell (natürlich gmpybin ich als Originalautor voreingenommen ;-).

Wenn Sie ein genaues Ergebnis wünschen, verwenden Sie sympy.binomial. Es scheint zweifellos die schnellste Methode zu sein.

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

Eine wörtliche Übersetzung der mathematischen Definition ist in vielen Fällen völlig ausreichend (wobei zu beachten ist, dass Python automatisch eine Arithmetik mit großen Zahlen verwendet):

from math import factorial

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

Bei einigen von mir getesteten Eingaben (z. B. n = 1000 r = 500) war dies mehr als zehnmal schneller als der eine Liner reduce, der in einer anderen (derzeit am höchsten bewerteten) Antwort vorgeschlagen wurde. Auf der anderen Seite wird es von dem von @JF Sebastian bereitgestellten Snippit übertroffen.

Ab Python 3.8sofort enthält die Standardbibliothek jetzt die math.combFunktion zur Berechnung des Binomialkoeffizienten:

math.comb (n, k)

Dies ist die Anzahl der Möglichkeiten, k Elemente aus n Elementen ohne Wiederholung auszuwählen
n! / (k! (n - k)!):

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

Hier ist eine andere Alternative. Dieser wurde ursprünglich in C ++ geschrieben, sodass er für eine Ganzzahl mit endlicher Genauigkeit (z. B. __int64) nach C ++ zurückportiert werden kann. Der Vorteil besteht darin, dass (1) nur ganzzahlige Operationen erforderlich sind und (2) das Aufblähen des ganzzahligen Werts durch aufeinanderfolgende Multiplikations- und Divisionspaare vermieden wird. Ich habe das Ergebnis mit Nas Banovs Pascal-Dreieck getestet, es erhält die richtige Antwort:

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

Begründung: Um die Anzahl der Multiplikationen und Divisionen zu minimieren, schreiben wir den Ausdruck wie folgt um

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

Um einen Multiplikationsüberlauf so weit wie möglich zu vermeiden, werden wir in der folgenden STRICT-Reihenfolge von links nach rechts auswerten:

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

Wir können zeigen, dass die in dieser Reihenfolge betriebene Ganzzahlarithmetik genau ist (dh kein Rundungsfehler).

Bei dynamischer Programmierung beträgt die zeitliche Komplexität Θ (n * m) und die räumliche Komplexität Θ (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]

Wenn Ihr Programm eine Obergrenze für n(sagen wir n <= N) hat und nCr wiederholt berechnen muss (vorzugsweise für >> NZeiten), kann die Verwendung von lru_cache zu einer enormen Leistungssteigerung führen:

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)

Das Erstellen des Caches (was implizit erfolgt) nimmt O(N^2)Zeit in Anspruch. Alle nachfolgenden Anrufe an nCrwerden zurückgegeben O(1).

Sie können 2 einfache Funktionen schreiben, die sich tatsächlich als 5-8-mal schneller herausstellen als mit scipy.special.comb . Tatsächlich müssen Sie keine zusätzlichen Pakete importieren, und die Funktion ist recht gut lesbar. Der Trick besteht darin, die zuvor berechneten Werte mithilfe der Memoisierung zu speichern und die Definition von nCr zu verwenden

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

Wenn wir mal vergleichen

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

Mit Sympy ist das ziemlich einfach.

import sympy

comb = sympy.binomial(n, r)

Verwenden Sie nur die mit Python vertriebene Standardbibliothek :

import itertools

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

Die direkte Formel erzeugt große ganze Zahlen, wenn n größer als 20 ist.

Also noch eine Antwort:

from math import factorial

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

kurz, genau und effizient, da dies große Python-Ganzzahlen vermeidet, indem es bei Longs bleibt.

Es ist genauer und schneller im Vergleich zu 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

Dies ist der @ KillerT2333-Code, der den integrierten Memoization Decorator verwendet.

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

Hier ist ein effizienter Algorithmus für Sie

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

Zum Beispiel nCr (30,7) = Fakt (30) / (Fakt (7) * Fakt (23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24) / (1 * 2 * 3 * 4 * 5 * 6 * 7)

Führen Sie also einfach die Schleife von 1 nach r aus, um das Ergebnis zu erhalten.

Das ist wahrscheinlich so schnell, wie Sie es in reinem Python für relativ große Eingaben tun können:

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

Diese Funktion ist sehr optimiert.

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