Статистика: комбинации в Python
sample - неправильный инструмент, как вы видели. Вместо этого вы можете использовать choices :
choices([10, 20, 30, 40, 50], k=4)
5 ответов
См. scipy.special.comb (scipy.misc.comb в старых версиях scipy). Когда Exact имеет значение False, он использует функцию gammaln для получения хорошей точности, не занимая много времени. В точном случае он возвращает целое число произвольной точности, вычисление которого может занять много времени.
Если хотите точные результаты и скорость, попробуйте gmpy - gmpy.comb должен делать именно то, что вы просите, и довольно быстро (из Конечно, как первоначальный автор gmpy , я предвзято ;-).
Быстрый поиск в коде Google дает (он использует формулу из ответа @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 () работает в 10 раз быстрее (проверено на все 0 <= (n, k) <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
Почему бы не написать это самому? Это one-liner или типа того:
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) )
Тест - печать треугольника Паскаля:
>>> 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
Вот эффективный алгоритм для Вас
for i = 1.....r
p = p * ( n - i ) / i
print(p)
, Например, NCR (30,7) =, факт (30) / (факт (7) * факт (23)) = (30 * 29 * 28 * 27 * 26 * 25 * 24) / (1 * 2 * 3 * 4 * 5 * 6 * 7)
Поэтому просто работает, цикл от 1 до r может получить результат.