Details
-
Bug
-
Status: Closed
-
Major
-
Resolution: Fixed
-
None
-
None
-
None
-
None
-
Operating System: other
Platform: PC
-
36105
Description
The inverseCumulativeProbability for the Poisson distribution is calculated
through the standard implementation in AbstractIntegerDistribution. This
implementation fails with a continue fraction stack overflow for large lambda.
Instead, I suggest the following override in PoissonDistributionImpl
/**
- Calculate the inverse cumulative probability for the Poisson distribution
- based on a Cornish-Fisher expansion approximation followed up by a line
- search.
- For the Cornish-Fisher expansion see Abramawitz and Stegun
- 'Handbook of Mathmatical Functions' pages 935 and 928
- @param prob the target probability
- @return the desired quantile
- @throws MathException when things go wrong
*/
private int inverseCumulativeProbability(double prob) throws MathException{
if (prob < 0.0 || prob >= 1.0)
throw new MathException("Probability must be in right-side open interval
[0.0, 1.0)");
if (prob == 0) return 0; // there is nothing to calculate
// Use the Cornish-Fisher Expansion with two terms to get a very good
approximation
// see Abramawitz and Stegun 'Handbook of Mathmatical Functions'
// pages 935 and 928
double mu = this.getMean(); // mean
double sigma = Math.sqrt(mu); // standard deviation
double gamma = 1.0/Math.sqrt(mu); // skewness
double z = new NormalDistributionImpl(0.0,
1.0).inverseCumulativeProbability(prob);
// this is the actual expansion
// the floor(... + 0.5) operation is the continuity correction
int y = (int) Math.floor(mu + sigma*(z + gamma*(z*z - 1.0)/6.0) + 0.5);
// Given this starting point, line search to the right or left.
// Bisection search is not necessary, since the approximation is rarely
// off by more than 1 or 2
z = this.cumulativeProbability
;
if ( z > prob) { // missed it to the right, search to the left
while(true)
} else { // missed it to the left, search to the right
while(true)
}
}