Merge pull request probmods#455 from longouyang/dont-approx-binomial-…

…score

Dont approx binomial score

src/dists.ad.js

@@ -574,6 +574,8 @@ function binomialG(x) {
   return (1 - (x * x) + (2 * x * Math.log(x))) / (d * d);
 }
 
+// see lemma 6.1 from Ahrens & Dieter's
+// Computer Methods for Sampling from Gamma, Beta, Poisson and Binomial Distributions
 function binomialSample(p, n) {
   var k = 0;
   var N = 10;
@@ -614,32 +616,9 @@ var Binomial = makeDistributionType({
     'use ad';
     var p = this.params.p;
     var n = this.params.n;
-    if (n > 20 && n * p > 5 && n * (1 - p) > 5) {
-      // large n, reasonable p approximation
-      var s = val;
-      var inv2 = 1 / 2;
-      var inv3 = 1 / 3;
-      var inv6 = 1 / 6;
-      if (s >= n) {
-        return -Infinity;
-      }
-      var q = 1 - p;
-      var S = s + inv2;
-      var T = n - s - inv2;
-      var d1 = s + inv6 - (n + inv3) * p;
-      var d2 = q / (s + inv2) - p / (T + inv2) + (q - inv2) / (n + 1);
-      d2 = d1 + 0.02 * d2;
-      var num = 1 + q * binomialG(S / (n * p)) + p * binomialG(T / (n * q));
-      var den = (n + inv6) * p * q;
-      var z = num / den;
-      var invsd = Math.sqrt(z);
-      z = d2 * invsd;
-      return gaussianScore(0, 1, z) + Math.log(invsd);
-    } else {
-      // exact formula
-      return (lnfact(n) - lnfact(n - val) - lnfact(val) +
-          val * Math.log(p) + (n - val) * Math.log(1 - p));
-    }
+    // exact formula
+    return (lnfact(n) - lnfact(n - val) - lnfact(val) +
+            val * Math.log(p) + (n - val) * Math.log(1 - p));
   },
   support: function() {
     return _.range(this.params.n).concat(this.params.n);