commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From bugzi...@apache.org
Subject DO NOT REPLY [Bug 36105] New: - [PATCH] inverseCumulativeProbability for PoissonDistributionImpl Fails for large lambda
Date Tue, 09 Aug 2005 17:44:53 GMT
DO NOT REPLY TO THIS EMAIL, BUT PLEASE POST YOUR BUG·
RELATED COMMENTS THROUGH THE WEB INTERFACE AVAILABLE AT
<http://issues.apache.org/bugzilla/show_bug.cgi?id=36105>.
ANY REPLY MADE TO THIS MESSAGE WILL NOT BE COLLECTED AND·
INSERTED IN THE BUG DATABASE.

http://issues.apache.org/bugzilla/show_bug.cgi?id=36105

           Summary: [PATCH] inverseCumulativeProbability for
                    PoissonDistributionImpl Fails for large lambda
           Product: Commons
           Version: unspecified
          Platform: PC
        OS/Version: other
            Status: NEW
          Severity: normal
          Priority: P2
         Component: Math
        AssignedTo: commons-dev@jakarta.apache.org
        ReportedBy: mikael@amazon.com


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(y);
		
    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;
        }
    }
}

-- 
Configure bugmail: http://issues.apache.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are the assignee for the bug, or are watching the assignee.

---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message