import numpy as np
import scipy.misc
def binomial_prob(k, n, p):
"""Probability that X = k for X ~ B(n, p)."""
n_choose_k = scipy.misc.comb(n, k) #, exact=0)
return n_choose_k * p**k * (1.0 - p)**(n - k)
def binomial_prob_range(low, high, n, p):
prob = 0.0
for k in xrange(low, high + 1):
prob += binomial_prob(k, n, p)
return prob
assert np.isclose(binomial_prob(5, 10, 5.0 / 12.0), 0.2138, atol=.001)
binomial_prob_range(500010000, 500020000, 1000000000, 0.5)
nan