Work in the log domain to compute combination and exponentiation functions and then raise them to exponent.
Something like this:
combination_num = range(k+1, n+1)
combination_den = range(1, n-k+1)
combination_log = np.log(combination_num).sum() - np.log(combination_den).sum()
p_k_log = k * np.log(p)
neg_p_K_log = (n - k) * np.log(1 - p)
p_log = combination_log + p_k_log + neg_p_K_log
probability = np.exp(p_log)
Gets rid of numeric underflow/overflow because of large numbers. On your example with n=450000
and p = 0.5, k = 17
, it returns p_log = -311728.4
, i. e., the log of final probability is pretty small and hence underflow occurs while taking np.exp
. However, you can still work with log probability.