Merge pull request #4 from leonmkim/bug_binom_pmf

fix implementation of binom_pmf to handle large k/n

Author: HarukiNishimura-TRI

binomial_cis/binomial_helper.py

@@ -31,10 +31,26 @@ def binom_pmf(k, n, p):
     Returns
     pmf: binomial pmf evaluated at k, n, p
     """
-    pmf = binom_coeff(n, k) * p**k * (1-p)**(n-k)
-    
-    # safeguard against floating point errors
-    pmf = min(1, max(0, pmf)) # unable to use np.clip()
+    if p == 1:
+        if k == n:
+            pmf = 1
+        else:
+            pmf = 0
+    elif p == 0:
+        if k == 0:
+            pmf = 1
+        else:
+            pmf = 0
+    else:
+        # Use the multiplicative formula for computational efficiency:
+        # https://en.wikipedia.org/wiki/Binomial_coefficient#Multiplicative_formula
+        log_binom_coeff = 0
+        for i in range(1, int(min(k, n - k)) + 1):
+            log_binom_coeff += np.log(n - i + 1) - np.log(i)
+        log_pmf = log_binom_coeff + k * np.log(p) + (n-k) * np.log1p(-p)
+        pmf = np.exp(log_pmf)
+        # safeguard against floating point errors
+        pmf = min(1, max(0, pmf)) # unable to use np.clip()
     return pmf