Statistik: kombinasi dengan Python

122

Saya perlu menghitung kombinatorial (nCr) dengan Python tetapi tidak dapat menemukan fungsi untuk melakukannya di math, numpyatau stat perpustakaan. Sesuatu seperti fungsi dari tipe:

comb = calculate_combinations(n, r)

Saya membutuhkan jumlah kombinasi yang mungkin, bukan kombinasi yang sebenarnya, jadi itertools.combinationssaya tidak tertarik.

Terakhir, saya ingin menghindari penggunaan faktorial, karena angka yang akan saya hitung kombinasinya bisa menjadi terlalu besar dan faktorialnya akan menjadi sangat besar.

Ini sepertinya pertanyaan yang BENAR-BENAR mudah dijawab, namun saya tenggelam dalam pertanyaan tentang menghasilkan semua kombinasi yang sebenarnya, yang bukan itu yang saya inginkan.

Morlock
sumber

Jawaban:

121

Lihat scipy.special.comb (scipy.misc.comb dalam versi scipy yang lebih lama). When exactis False, ia menggunakan fungsi gammaln untuk mendapatkan presisi yang baik tanpa menghabiskan banyak waktu. Dalam kasus yang tepat, ini mengembalikan bilangan bulat presisi arbitrer, yang mungkin membutuhkan waktu lama untuk dihitung.

Jouni K. Seppänen
sumber
5
scipy.misc.combtidak digunakan lagi scipy.special.combkarena versi 0.10.0.
Dilawar
120

Mengapa tidak menulisnya sendiri? Ini satu baris atau semacamnya:

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

Tes - mencetak segitiga 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. diedit untuk diganti int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) dengan int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))sehingga tidak akan salah untuk N / K besar

Nas Banov
sumber
26
+1 untuk menyarankan untuk menulis sesuatu yang sederhana, untuk menggunakan pengurangan, dan untuk demo keren dengan segitiga pascal
jon_darkstar
6
-1 karena jawaban ini salah: cetak faktorial (54) / (faktorial (54 - 27)) / faktorial (27) == nCk (54, 27) menghasilkan Salah.
robert king
3
@robertking - Oke, Anda berdua picik dan secara teknis benar. Apa yang saya lakukan itu dimaksudkan sebagai ilustrasi bagaimana menulis fungsinya sendiri; saya tahu itu tidak akurat untuk N dan K yang cukup besar karena presisi floating point. Tapi kita bisa memperbaikinya - lihat di atas, sekarang seharusnya tidak salah untuk angka besar
Nas Banov
9
Ini mungkin akan cepat di Haskell, tapi sayangnya tidak dengan Python. Sebenarnya ini cukup lambat dibandingkan dengan banyak jawaban lainnya, misalnya @Alex Martelli, JF Sebastian, dan saya sendiri.
Todd Owen
9
Untuk Python 3, saya juga harus from functools import reduce.
Velizar Hristov
52

Pencarian cepat di kode google memberikan (menggunakan rumus dari jawaban @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()adalah 10 kali lebih cepat (diuji pada semua 0 <= (n, k) <1e3 pasangan) dibandingkan scipy.misc.comb()jika Anda membutuhkan jawaban yang tepat.

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
sumber
Solusi bagus yang tidak memerlukan pkg apa pun
Edward Newell
2
FYI: Rumus yang disebutkan ada di sini: en.wikipedia.org/wiki/…
jmiserez
chooseFungsi ini seharusnya mendapatkan lebih banyak suara! Python 3.8 memiliki math.comb, tetapi saya harus menggunakan Python 3.6 untuk tantangan dan tidak ada implementasi yang memberikan hasil yang tepat untuk integer yang sangat besar. Yang ini melakukannya dan melakukannya dengan cepat!
sambungkan kembali
42

Jika Anda ingin hasil yang tepat dan kecepatan, coba gmpy - gmpy.combharus melakukan apa yang Anda minta, dan itu cukup cepat (tentu saja, sebagai gmpy's penulis asli, saya sedang bias ;-).

Alex Martelli
sumber
6
Memang, gmpy2.comb()10 kali lebih cepat daripada choose()dari jawaban saya untuk kode: di for k, n in itertools.combinations(range(1000), 2): f(n,k)mana f()ada gmpy2.comb()atau choose()di Python 3.
jfs
Karena Anda adalah pembuat paketnya, saya akan membiarkan Anda memperbaiki tautan yang rusak sehingga mengarah ke tempat yang tepat ....
SeldomNeedy
@SeldomNeedy, link ke code.google.com adalah salah satu tempat yang tepat (meskipun situs ini dalam mode arsip sekarang). Tentu saja dari sana mudah untuk menemukan lokasi github, github.com/aleaxit/gmpy , dan yang PyPI, pypi.python.org/pypi/gmpy2 , karena terhubung ke keduanya! -)
Alex Martelli
@AlexMartelli Maaf atas kebingungannya. Halaman ini menampilkan 404 jika javascript telah (secara selektif) dinonaktifkan. Saya kira itu untuk mencegah AI nakal memasukkan sumber Proyek Google Code yang diarsipkan dengan mudah?
JarangNeedy
28

Jika Anda menginginkan hasil yang tepat, gunakan sympy.binomial. Tampaknya ini metode tercepat, tangan ke bawah.

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
sumber
22

Terjemahan literal dari definisi matematika cukup memadai dalam banyak kasus (mengingat bahwa Python akan secara otomatis menggunakan aritmatika angka besar):

from math import factorial

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

Untuk beberapa masukan yang saya uji (misalnya n = 1000 r = 500) ini lebih dari 10 kali lebih cepat daripada satu baris yang reducedisarankan dalam jawaban lain (saat ini dengan suara terbanyak). Di sisi lain, itu dilakukan oleh cuplikan yang disediakan oleh @JF Sebastian.

Todd Owen
sumber
11

Memulai Python 3.8, pustaka standar sekarang menyertakan math.combfungsi untuk menghitung koefisien binomial:

math.comb (n, k)

yaitu banyaknya cara untuk memilih k item dari n item tanpa pengulangan
n! / (k! (n - k)!):

import math
math.comb(10, 5) # 252
Xavier Guihot
sumber
10

Berikut alternatif lain. Yang ini aslinya ditulis dalam C ++, sehingga dapat di-backport ke C ++ untuk integer presisi hingga (mis. __Int64). Keuntungannya adalah (1) hanya melibatkan operasi integer, dan (2) menghindari pembengkakan nilai integer dengan melakukan pasangan perkalian dan pembagian yang berurutan. Saya telah menguji hasilnya dengan segitiga Pascal Nas Banov, hasilnya benar:

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

Rasional: Untuk meminimalkan # perkalian dan pembagian, kami menulis ulang ekspresi sebagai

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

Untuk menghindari kelimpahan perkalian sebanyak mungkin, kami akan mengevaluasi dalam urutan STRICT berikut, dari kiri ke kanan:

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

Kita dapat menunjukkan bahwa aritmatika integer yang dioperasikan dalam urutan ini adalah tepat (yaitu tidak ada kesalahan pembulatan).

Wirawan Purwanto
sumber
5

Menggunakan pemrograman dinamis, kompleksitas waktu adalah Θ (n * m) dan kompleksitas ruang Θ (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
sumber
4

Jika program Anda memiliki batas atas n(katakanlah n <= N) dan perlu menghitung nCr berulang kali (sebaiknya untuk >> Nkali), menggunakan lru_cache dapat memberi Anda peningkatan kinerja yang besar:

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)

Membangun cache (yang dilakukan secara implisit) membutuhkan O(N^2)waktu. Setiap panggilan berikutnya ke nCrakan kembali masuk O(1).

yzn-pku
sumber
4

Anda dapat menulis 2 fungsi sederhana yang ternyata sekitar 5-8 kali lebih cepat daripada menggunakan scipy.special.comb . Faktanya, Anda tidak perlu mengimpor paket tambahan apa pun, dan fungsinya cukup mudah dibaca. Triknya adalah dengan menggunakan memoization untuk menyimpan nilai yang dihitung sebelumnya, dan menggunakan definisi 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))

Jika kita membandingkan waktu

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
sumber
Hari-hari ini ada penghias memoize di functools bernama lru_cache yang dapat menyederhanakan kode Anda?
landak gila
2

Sangat mudah dengan sympy.

import sympy

comb = sympy.binomial(n, r)
Polisi
sumber
2

Hanya menggunakan pustaka standar yang didistribusikan dengan Python :

import itertools

def nCk(n, k):
    return len(list(itertools.combinations(range(n), k)))
MarianD
sumber
3
Saya tidak berpikir kompleksitas waktunya (dan penggunaan memori) dapat diterima.
xmcp
2

Rumus langsung menghasilkan bilangan bulat besar jika n lebih besar dari 20.

Jadi, tanggapan lain:

from math import factorial

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

pendek, akurat dan efisien karena ini menghindari bilangan bulat besar python dengan bertahan dengan rindu.

Ini lebih akurat dan lebih cepat jika dibandingkan dengan 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.dll
sumber
Ini salah! Jika n == r, hasilnya harus 1. Kode ini mengembalikan 0.
reyammer
Lebih tepatnya, itu seharusnya range(n-r+1, n+1)bukan range(n-r,n+1).
reyammer
1

Ini adalah kode @ killerT2333 menggunakan dekorator memoisasi bawaan.

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))
landak gila
sumber
1

Berikut adalah algoritma yang efisien untuk Anda

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

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

Jadi jalankan saja loop dari 1 ke r bisa mendapatkan hasilnya.

kta
sumber
0

Itu mungkin secepat Anda bisa melakukannya dengan python murni untuk input yang cukup besar:

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
sumber
0

Fungsi ini sangat dioptimalkan.

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
sumber