Python Binomial Coefficient

Python Binomial Coefficient

This question is old but as it comes up high on search results I will point out that scipy has two functions for computing the binomial coefficients:

1. scipy.special.binom()

2. scipy.special.comb()

   import scipy.special
   
   # the two give the same results 
   scipy.special.binom(10, 5)
   # 252.0
   scipy.special.comb(10, 5)
   # 252.0
   
   scipy.special.binom(300, 150)
   # 9.375970277281882e+88
   scipy.special.comb(300, 150)
   # 9.375970277281882e+88
   
   # ...but with `exact == True`
   scipy.special.comb(10, 5, exact=True)
   # 252
   scipy.special.comb(300, 150, exact=True)
   # 393759702772827452793193754439064084879232655700081358920472352712975170021839591675861424
   

Note that scipy.special.comb(exact=True) uses Python integers, and therefore it can handle arbitrarily large results!

Speed-wise, the three versions give somewhat different results:

num = 300

%timeit [[scipy.special.binom(n, k) for k in range(n + 1)] for n in range(num)]
# 52.9 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit [[scipy.special.comb(n, k) for k in range(n + 1)] for n in range(num)]
# 183 ms ± 814 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)each)

%timeit [[scipy.special.comb(n, k, exact=True) for k in range(n + 1)] for n in range(num)]
# 180 ms ± 649 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

(and for n = 300, the binomial coefficients are too large to be represented correctly using float64 numbers, as shown above).

printing binomial coefficient using numpy

You may want to use: scipy.special.binom()

or, since Python 3.8: math.comb()

---

EDIT

I am not quite sure why you would not want to use SciPy but you are OK with NumPy, as SciPy is a well-established library from essentially the same folks developing NumPy.

Anyway, here a couple of other methods:

• using math.factorial:

import math

def binom(n, k):
    return math.factorial(n) // math.factorial(k) // math.factorial(n - k)

• using prod() and math.factorial() (theoretically more efficient, but not in practice):

def prod(items, start=1):
    for item in items:
        start *= item
    return start

def binom_simplified(n, k):
    if k > n - k:
        return prod(range(k + 1, n + 1)) // math.factorial(n - k)
    else:
        return prod(range(n - k + 1, n + 1)) // math.factorial(k)

• using numpy.prod():

import numpy as np

def binom_np(n, k):
    return 1 if k == 0 or k == n else np.prod([(n + 1 - i) / i for i in range(1, k + 1)])

---

Speed-wise, scipy.special.binom() is the fastest by far and large, but if you need the exact value also for very large numbers, you may prefer binom() (somewhat surprisingly even over math.comb()).

%timeit scipy.special.binom(600, 298)
# 1000000 loops, best of 3: 1.56 µs per loop
print(scipy.special.binom(600, 298))
# 1.3332140543730587e+179

%timeit math.comb(600, 298)
# 10000 loops, best of 3: 75.6 µs per loop
print(math.binom(600, 298))
# 133321405437268991724586879878020905773601074858558174180536459530557427686938822154484588609548964189291743543415057988154692680263088796451884071926401665548516571367537285901600

%timeit binom(600, 298)
# 10000 loops, best of 3: 36.5 µs per loop
print(binom(600, 298))
# 133321405437268991724586879878020905773601074858558174180536459530557427686938822154484588609548964189291743543415057988154692680263088796451884071926401665548516571367537285901600

%timeit binom_np(600, 298)
# 10000 loops, best of 3: 45.8 µs per loop
print(binom_np(600, 298))
# 1.3332140543726893e+179

%timeit binom_simplified(600, 298)
# 10000 loops, best of 3: 41.9 µs per loop
print(binom_simplified(600, 298))
# 133321405437268991724586879878020905773601074858558174180536459530557427686938822154484588609548964189291743543415057988154692680263088796451884071926401665548516571367537285901600

How to perform binomial-coefficient and factorial calculation with more precision?

The problem is that you have exact=False in the factorial.

>>> import numpy as np
>>> from decimal import *
>>> import scipy.special
>>> from scipy.special import factorial
>>> getcontext().prec = 30
>>> 
>>> i = 500
>>> sum(np.array([scipy.special.comb(Decimal(i), (r), exact=True)*pow(-1, r)/Decimal(factorial(r, exact=False)) for r in range(i+1)]))
Decimal('-7.13859062099388393889008217957')
>>> sum(np.array([scipy.special.comb(Decimal(i), (r), exact=True)*pow(-1, r)/Decimal(factorial(r, exact=True)) for r in range(i+1)]))
Decimal('0.196589352363439561009074161963')

Python: print out the combinations of heads and tails based on binomial coefficient

Here's a straightforward solution using itertools.combinations:

from itertools import combinations

def coin_flips(n, k):
    for c in combinations(range(n), k):
        out = ['H'] * n
        for i in c: out[i] = 'T'
        yield tuple(out)

Example:

>>> for c in coin_flips(3, 2):
...     print(c)
... 
('T', 'T', 'H')
('T', 'H', 'T')
('H', 'T', 'T')

---

Related Topics