jstat/jstat — Issue #35: gammap approximation breaks down for small k
Author: richarddmorey
Created: 2012-02-04 16:04:49+00:00
State: open

The gammap function has issues with moderately small values of a, where the function doesn't increase monotonically as x increases. For instance,

jStat( 0.5, 3, 10 ).gamma( .15, 2 ).cdf().toArray();

yields

0.8253665805712278,
0.84370429375926,
0.8380937268336934,
0.8179624512257149,
0.7883836589089586,
0.75262694142181,
0.7129903285130516,
0.9738102454941774,
0.9785160308224301,
0.9823072094250467

If you look at a plot, it appears to have a discontinuity. The problem gets worse for progressively smaller values of a. The problem appears to be in lines 128-132 of special.js:

} else if ( x < a + 1 ) {
for ( ; i <= ITMAX; i++ ) {
sum += del *= x / ++ap;
}
endval = sum * Math.exp( -x + a * Math.log( x ) - ( aln ));
}

because it appears to happen when x < a + 1. I had a brief look at the Numerical Recipes algorithm but didn't immediately see what might be wrong. Could be not enough terms (ITMAX) in the series?