jstat/jstat — Issue #34: binomial cdf returns values >1
Author: richarddmorey
Created: 2012-02-03 17:02:33+00:00
State: open

There's some imprecision here creeping in from somewhere.

> jStat.binomial.cdf(21,22,.3)
1.000000000005807

R returns:

> print(pbinom(21,22,.3),digits=22)
[1] 0.9999999999968618435986

The problem appears to be combinationln, because using no logs:

> var d = [];
> for(var i=0;i<22;i++) d[i] = jStat.combination(22,i) * Math.pow(.3,i) * Math.pow(1-.3,22-i)
> jStat.sum(d);
0.9999999999968605

agrees with R, and using logs for only the likelihood term:

> for(var i=0;i<22;i++) d[i] = jStat.combination(22,i) * Math.exp(i * Math.log(.3) + (22-i)*Math.log(1-.3))
1.6108943932619987e-10
> jStat.sum(d);
0.99999999999686

also agrees with R. jStat.combination() fails, however, for moderately sized N and k:

> jStat.combination(200,150)
Infinity

So the logs are necessary for large N or small/large p situations. We should probably add a check to makes sure that N or p are "appropriately" sized before taking logs. And maybe if the cdf ever goes above 1 or below 0, we make it 1 or 0, respectively?