2010-06-11 4 views
88

Мне нужно вычислить комбинаторные (nCr) в Python, но не могу найти функцию для этого в math, numpy или stat библиотеках. Что-то вроде функции типа:Статистика: комбинации в Python

comb = calculate_combinations(n, r) 

мне нужно число возможных комбинаций, а не фактические комбинации, так itertools.combinations меня не интересует.

Наконец, я хочу избежать использования факториалов, так как числа, которые я буду вычислять для комбинаций, могут стать слишком большими, а факториалы будут чудовищными.

Это кажется ДЕЙСТВИТЕЛЬНО легким ответом на вопрос, однако я погружаюсь в вопросы о создании всех реальных комбинаций, чего я не хочу. :)

Большое спасибо

ответ

89

scipy.special.comb См (scipy.misc.comb в более старых версиях SciPy). Когда exact является ложным, он использует функцию gammaln для получения хорошей точности, не затрачивая много времени. В точном случае он возвращает целое число с произвольной точностью, что может занять много времени для вычисления.

+0

'scipy.misc.comb' устарел в пользу' scipy.special.comb' с версии '0.10.0'. – Dilawar

43

Быстрый поиск на Google Code дает (он использует формулу из @Mark Byers's answer):

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() в 10 раз быстрее (проверено на всех 0 < = (п, к) < 1E3 пар), чем scipy.misc.comb() если вам нужен точный ответ.

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

Хорошее решение, которое не требует каких-либо pkg –

+2

FYI: Указанная формула приведена здесь: https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula – jmiserez

40

Если вы хотите точные результаты и скорость, попробуйте gmpy - gmpy.comb должны делать то, что вы просите, и это довольно быстро (конечно, как gmpy «s оригинальный автор, я утра Пристрастный ;-).

+6

Действительно, 'gmpy2.comb()' в 10 раз быстрее, чем 'select()' из моего ответа для кода: 'для k, n в itertools.combinations (range (1000), 2): f (n, k)' где 'f()' либо 'gmpy2.comb () 'или' choose() 'на Python 3. – jfs

+0

Поскольку вы являетесь автором пакета, я разрешу * вы * исправить неработающую ссылку, чтобы она указывала на нужное место .... – SeldomNeedy

+0

@SeldomNeedy, ссылка на code.google.com является * одним * правильным местом (хотя сейчас сайт находится в режиме архивации). Конечно, оттуда легко найти местоположение github, https://github.com/aleaxit/gmpy и PyPI, https://pypi.python.org/pypi/gmpy2, поскольку он связывается с обоими! -) –

97

Почему бы не написать это самостоятельно? Это один вкладыш или такой:

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 - печать треугольник Паскаля:

>>> 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. отредактированный для замены int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1))) с номером int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1)), поэтому он не будет ошибочным для больших N/K

+22

+1 за предложение написать что-то простое, для использования сокращения, и для классной демонстрации с бонусными баллонами с паскалем –

+0

, если вы замените 'range (..)' на 'xrange (...)' в python 2.x; -) – hochl

+0

Я думаю, что плакат искал более эффективное решение, адаптированное к большим входам (поэтому использование цикла for, особенно с диапазоном, не будет работать) – Philip

8

Вот еще один вариант. Это было первоначально написано на C++, поэтому оно может быть обращено к C++ для цепочки конечной точности (например, __int64). Преимущество состоит в том, что (1) он включает только целые операции, и (2) он избегает раздувания целочисленного значения, выполняя последовательные пары умножения и деления.Я проверил результат с треугольником Паскаля Nas Банов, он получает правильный ответ:

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 

Обоснование: Для того, чтобы свести к минимуму # умножений и делений, перепишем выражение как

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

Чтобы избежать переполнение умножение как можно больше, мы будем оценивать в следующем ПРЯМУЮ порядке слева направо:

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

мы можем показать, что целое arithmatic работают в этом порядок точный (т. нет ошибки округления).

18

Буквальный перевод математического определения является вполне адекватной в большинстве случаев (напомним, что Python автоматически будет использовать большое число арифметических операций):

from math import factorial 

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

Для некоторых входов я тестировал (например, п = 1000 г = 500), это было более чем в 10 раз быстрее, чем один лайнер reduce, предложенный в другом (в настоящее время наиболее проголосовавшем) ответе. С другой стороны, он не выполняется сниппитом, предоставленным @ J.F. Себастьян.

22

Если вы хотите получить точный результат, используйте sympy.binomial. Кажется, это самый быстрый метод, руки вниз.

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 
4

Использование динамического программирования, время сложность Θ (п * м) и пространство сложность Θ (м):

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

Это, вероятно, так же быстро, как вы можете сделать это в чистом питоне для разумных большие входы:

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 
1

непосредственная формула производит большие целые числа, когда п больше, чем 20.

So, Еще один ответ:

from math import factorial 

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

короткий, быстрый и эффективный.

+0

Это неправильно! Если n == r, результат должен быть 1. Этот код возвращает 0. – reyammer

+0

Точнее, это должно быть 'range (n-r + 1, n + 1)' вместо 'range (n-r, n + 1)'. – reyammer

1

Это довольно легко с симпатичным.

import sympy 

comb = sympy.binomial(n, r) 
1

Использование только стандартной библиотеки распределены с Python:

import itertools 

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

Я не думаю, что его временная сложность (и использование памяти) приемлема. – xmcp

0

Если ваша программа имеет верхнюю границу n (скажем n <= N) и нужно повторно вычислить Ncr (предпочтительно для >>N раз), используя lru_cache, может дать вам огромное повышение производительности:

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) 

Построение кэша (что делается неявно) занимает до O(N^2) времени. Любые последующие звонки в nCr возвращаются в O(1).

Смежные вопросы