Uploaded image for project: 'Commons Math'
  1. Commons Math
  2. MATH-92

[math][PATCH] inverseCumulativeProbability for PoissonDistributionImpl Fails for large lambda

    Details

    • Type: Bug
    • Status: Closed
    • Priority: Major
    • Resolution: Fixed
    • Affects Version/s: None
    • Fix Version/s: None
    • Labels:
      None
    • Environment:

      Operating System: other
      Platform: PC

    • Bugzilla Id:
      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)

      { 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; }

      }
      }

        Attachments

          Activity

            People

            • Assignee:
              Unassigned
              Reporter:
              mikael@amazon.com Mikael Weigelt
            • Votes:
              0 Vote for this issue
              Watchers:
              0 Start watching this issue

              Dates

              • Created:
                Updated:
                Resolved: