Skip to content
Permalink
Browse files

Allow binomial.cdf to handle large numbers

Also add some checks to make sure NaN is returned when it should be.

Fixes: #228
  • Loading branch information...
trevnorris committed Jun 4, 2019
1 parent 28c5207 commit eeb696a4390c31d85525595a5d5bce1f1e5574ef
Showing with 59 additions and 11 deletions.
  1. +51 −11 src/distribution.js
  2. +8 −0 test/distribution/binomial-test.js
@@ -724,6 +724,34 @@ jStat.extend(jStat.uniform, {
});


// Got this from http://www.math.ucla.edu/~tom/distributions/binomial.html
function betinc(x, a, b, eps) {
var a0 = 0;
var b0 = 1;
var a1 = 1;
var b1 = 1;
var m9 = 0;
var a2 = 0;
var c9;

while (Math.abs((a1 - a2) / a1) > eps) {
a2 = a1;
c9 = -(a + m9) * (a + b + m9) * x / (a + 2 * m9) / (a + 2 * m9 + 1);
a0 = a1 + c9 * a0;
b0 = b1 + c9 * b0;
m9 = m9 + 1;
c9 = m9 * (b - m9) * x / (a + 2 * m9 - 1) / (a + 2 * m9);
a1 = a0 + c9 * a1;
b1 = b0 + c9 * b1;
a0 = a0 / b1;
b0 = b0 / b1;
a1 = a1 / b1;
b1 = 1;
}

return a1 / a;
}


// extend uniform function with static methods
jStat.extend(jStat.binomial, {
@@ -734,18 +762,30 @@ jStat.extend(jStat.binomial, {
},

cdf: function cdf(x, n, p) {
var binomarr = [],
k = 0;
if (x < 0) {
var betacdf;
var eps = 1e-10;

if (x < 0)
return 0;
}
if (x < n) {
for (; k <= x; k++) {
binomarr[ k ] = jStat.binomial.pdf(k, n, p);
}
return jStat.sum(binomarr);
}
return 1;
if (x >= n)
return 1;
if (p < 0 || p > 1 || n <= 0)
return NaN;

x = Math.floor(x);
var z = p;
var a = x + 1;
var b = n - x;
var s = a + b;
var bt = Math.exp(jStat.gammaln(s) - jStat.gammaln(b) -
jStat.gammaln(a) + a * Math.log(z) + b * Math.log(1 - z));

if (z < (a + 1) / (s + 2))
betacdf = bt * betinc(z, a, b, eps);
else
betacdf = 1 - bt * betinc(1 - z, b, a, eps);

return Math.round((1 - betacdf) * (1 / eps)) / (1 / eps);
}
});

@@ -20,6 +20,14 @@ suite.addBatch({
var tol = 0.0000001;
assert.epsilon(tol, jStat.binomial.cdf(10, 25, 0.5), 0.2121781);
assert.epsilon(tol, jStat.binomial.cdf(50, 1000, 0.05), 0.537529);

assert.ok(Number.isNaN(jStat.binomial.cdf(10, 12, -1)));
assert.ok(Number.isNaN(jStat.binomial.cdf(10, 12, 2)));
assert.equal(jStat.binomial.cdf(12, 10, 0.5), 1);
assert.equal(jStat.binomial.cdf(12, -1, 0.5), 1);
assert.epsilon(tol,
jStat.binomial.cdf(101073, 101184, 0.9988219676207195),
0.7857313651);
}
}
});

0 comments on commit eeb696a

Please sign in to comment.
You can’t perform that action at this time.