[MATH-92] [math][PATCH] inverseCumulativeProbability for PoissonDistributionImpl Fails for large lambda

• Type: Bug
• Status: Closed
• Priority: Major
• Resolution: Fixed

Environment

Operating System: other
Platform: PC

• Bugzilla Id: 36105
• Assignee: Unassigned
• Reporter: Mikael Weigelt
• Created: 10/Aug/05 02:44
• Updated: 23/Apr/09 02:26
• Resolved: 23/Apr/09 02:26

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)

{ if (y == 0 || this.cumulativeProbability(y - 1) < prob) return y; y--; }

} else { // missed it to the left, search to the right
while(true)

{ y++; if (this.cumulativeProbability(y) >= prob) return y; }

}
}