Why does the following code that calls comb produce the wrong output?
from scipy.special import comb
from sympy import binomial
def coup(n=100, m=10):
expectation = 0
for j in range(1, n + 1):
coeff = ((-1)**(j-1)) * comb(n, j)
denominator = 1 - (comb(n - j, m) / comb(n, m))
expectation += coeff / denominator
print(expectation)
# output: 1764921496446.9807
# correct output when using binomial: 49.9445660566641
You need to include the exact parameter as True for the comb function if you want it to have full precision and not just floating point precision. (It defaults to False) for example -
comb(100,15) - binomial(100,15)
would give an output of -32, but
comb(100,15, exact=True) - binomial(100,15)
gives an output of 0
This is as per SciPy docs - https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.comb.html