PoissonDist.java

/*
 * Class:        PoissonDist
 * Description:  Poisson distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.probdist;
 
import umontreal.ssj.util.*;

public class PoissonDist extends DiscreteDistributionInt {
 
   private double lambda;


   public static double MAXLAMBDA = 100000;


   public PoissonDist(double lambda) {
      setLambda(lambda);
   }
 
   public double prob(int x) {
      if (x < 0)
         return 0.0;
      if (pdf == null)
         return prob(lambda, x);
      if (x > xmax || x < xmin)
         return prob(lambda, x);
      return pdf[x - xmin];
   }
 
   public double cdf(int x) {
      double Sum = 0.0;
      // int j;
 
      if (x < 0)
         return 0.0;
      if (lambda == 0.0)
         return 1.0;
 
      /*
       * For large lambda, we use the Chi2 distribution according to the exact
       * relation, with 2x + 2 degrees of freedom
       * 
       * cdf (lambda, x) = 1 - chiSquare (2x + 2, 2*lambda)
       * 
       * which equals also 1 - gamma (x + 1, lambda)
       */
      if (cdf == null)
         return GammaDist.barF(x + 1.0, 15, lambda);
 
      if (x >= xmax)
         return 1.0;
 
      if (x < xmin) {
         // Sum a few terms to get a few decimals far in the lower tail. One
         // could also call GammaDist.barF instead.
         final int RMAX = 20;
         int i;
         double term = prob(lambda, x);
         Sum = term;
         i = x;
         while (i > 0 && i >= x - RMAX) {
            term = term * i / lambda;
            i--;
            Sum += term;
         }
         return Sum;
      }
 
      if (x <= xmed)
         return cdf[x - xmin];
      else
         // We keep the complementary distribution in the upper part of cdf
         return 1.0 - cdf[x + 1 - xmin];
   }
 
   public double barF(int x) {
      /*
       * poisson (lambda, x) = 1 - cdf (lambda, x - 1)
       */
 
      if (x <= 0)
         return 1.0;
 
      /*
       * For large lambda, we use the Chi2 distribution according to the exact
       * relation, with 2x + 2 degrees of freedom
       * 
       * cdf (lambda, x) = 1 - ChiSquare.cdf (2x + 2, 2*lambda) cdf (lambda, x) = 1 -
       * GammaDist.cdf (x + 1, lambda)
       */
 
      if (cdf == null)
         return GammaDist.cdf((double) x, 15, lambda);
 
      if (x > xmax)
//         return GammaDist.cdf ((double)x, 15, lambda);
         return PoissonDist.barF(lambda, x);
      if (x <= xmin)
         return 1.0;
      if (x > xmed)
         // We keep the complementary distribution in the upper part of cdf
         return cdf[x - xmin];
      else
         return 1.0 - cdf[x - 1 - xmin];
   }
 
   public int inverseFInt(double u) {
      if ((cdf == null) || (u <= EPSILON))
         return inverseF(lambda, u);
      return super.inverseFInt(u);
   }
 
   public double getMean() {
      return PoissonDist.getMean(lambda);
   }
 
   public double getVariance() {
      return PoissonDist.getVariance(lambda);
   }
 
   public double getStandardDeviation() {
      return PoissonDist.getStandardDeviation(lambda);
   }

   public static double prob(double lambda, int x) {
      if (x < 0)
         return 0.0;
 
      if (lambda >= 100.0) {
         if ((double) x >= 10.0 * lambda)
            return 0.0;
      } else if (lambda >= 3.0) {
         if ((double) x >= 100.0 * lambda)
            return 0.0;
      } else {
         if ((double) x >= 200.0 * Math.max(1.0, lambda))
            return 0.0;
      }
 
      final double LAMBDALIM = 20.0;
      double Res;
      if (lambda < LAMBDALIM && x <= 100)
         Res = Math.exp(-lambda) * Math.pow(lambda, x) / Num.factorial(x);
      else {
         double y = x * Math.log(lambda) - Num.lnGamma(x + 1.0) - lambda;
         Res = Math.exp(y);
      }
      return Res;
   }

   public static double cdf(double lambda, int x) {
      /*
       * On our machine, computing a value using gamma is faster than the naive
       * computation for lambdalim > 200.0, slower for lambdalim < 200.0
       */
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
      if (lambda == 0.0)
         return 1.0;
      if (x < 0)
         return 0.0;
 
      if (lambda >= 100.0) {
         if ((double) x >= 10.0 * lambda)
            return 1.0;
      } else {
         if ((double) x >= 100.0 * Math.max(1.0, lambda))
            return 1.0;
      }
 
      /*
       * If lambda > LAMBDALIM, use the Chi2 distribution according to the exact
       * relation, with 2x + 2 degrees of freedom
       * 
       * poisson (lambda, x) = 1 - chiSquare (2x + 2, 2*lambda)
       * 
       * which also equals 1 - gamma (x + 1, lambda)
       */
 
      final double LAMBDALIM = 200.0;
      if (lambda > LAMBDALIM)
         return GammaDist.barF(x + 1.0, 15, lambda);
 
      if (x >= lambda)
         return 1 - PoissonDist.barF(lambda, x + 1);
 
      // Naive computation: sum all prob. from i = x
      double sum = 1;
      double term = 1;
      for (int j = 1; j <= x; j++) {
         term *= lambda / j;
         sum += term;
      }
      return sum * Math.exp(-lambda);
   }

   public static double barF(double lambda, int x) {
      if (lambda < 0)
         throw new IllegalArgumentException("lambda < 0");
      if (x <= 0)
         return 1.0;
 
      if (lambda >= 100.0) {
         if ((double) x >= 10.0 * lambda)
            return 0.0;
      } else {
         if ((double) x >= 100 + 100.0 * Math.max(1.0, lambda))
            return 0.0;
      }
 
      /*
       * If lambda > LAMBDALIM, we use the Chi2 distribution according to the exact
       * relation, with 2x + 2 degrees of freedom
       * 
       * cdf (lambda, x) = 1 - ChiSquare.cdf (2x + 2, 2*lambda)
       * 
       * which also equals 1 - GammaDist.cdf (x + 1, lambda)
       */
 
      final double LAMBDALIM = 200.0;
      if (lambda > LAMBDALIM)
         return GammaDist.cdf((double) x, 15, lambda);
 
      if (x <= lambda)
         return 1.0 - PoissonDist.cdf(lambda, x - 1);
 
      // Naive computation: sum all prob. from i = x to i = oo
      double term, sum;
      final int IMAX = 20;
 
      // Sum at least IMAX prob. terms from i = s to i = oo
      sum = term = PoissonDist.prob(lambda, x);
      int i = x + 1;
      while (term > EPSILON || i <= x + IMAX) {
         term *= lambda / i;
         sum += term;
         i++;
      }
      return sum;
   }

   public static int inverseF(double lambda, double u) {
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException("u is not in range [0,1]");
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
      if (u >= 1.0)
         return Integer.MAX_VALUE;
      if (u <= prob(lambda, 0))
         return 0;
      int i;
 
      final double LAMBDALIM = 700.0;
      if (lambda < LAMBDALIM) {
         double sumprev = -1.0;
         double term = Math.exp(-lambda);
         double sum = term;
         i = 0;
         while (sum < u && sum > sumprev) {
            i++;
            term *= lambda / i;
            sumprev = sum;
            sum += term;
         }
         return i;
 
      } else {
         i = (int) lambda;
         double term = PoissonDist.prob(lambda, i);
         while ((term >= u) && (term > Double.MIN_NORMAL)) {
            i /= 2;
            term = PoissonDist.prob(lambda, i);
         }
         if (term <= Double.MIN_NORMAL) {
            i *= 2;
            term = PoissonDist.prob(lambda, i);
            while (term >= u && (term > Double.MIN_NORMAL)) {
               term *= i / lambda;
               i--;
            }
         }
         int mid = i;
         double sum = term;
         double termid = term;
 
         while (term >= EPSILON * u && i > 0) {
            term *= i / lambda;
            sum += term;
            i--;
         }
 
         term = termid;
         i = mid;
         double prev = -1;
         if (sum < u) {
            while ((sum < u) && (sum > prev)) {
               i++;
               term *= lambda / i;
               prev = sum;
               sum += term;
            }
         } else {
            // The computed CDF is too big so we substract from it.
            sum -= term;
            while (sum >= u) {
               term *= i / lambda;
               i--;
               sum -= term;
            }
         }
      }
      return i;
   }

   public static double[] getMLE(int[] x, int n) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
 
      double parameters[];
      parameters = new double[1];
      double sum = 0.0;
      for (int i = 0; i < n; i++) {
         sum += x[i];
      }
 
      parameters[0] = (double) sum / (double) n;
      return parameters;
   }

   public static PoissonDist getInstanceFromMLE(int[] x, int n) {
      double parameters[] = getMLE(x, n);
      return new PoissonDist(parameters[0]);
   }

   public static double getMean(double lambda) {
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
 
      return lambda;
   }

   public static double getVariance(double lambda) {
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
 
      return lambda;
   }

   public static double getStandardDeviation(double lambda) {
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
 
      return Math.sqrt(lambda);
   }

   public double getLambda() {
      return lambda;
   }

   public void setLambda(double lambda) {
      supportA = 0;
      if (lambda < 0.0)
         throw new IllegalArgumentException("lambda < 0");
      this.lambda = lambda;
 
      // For lambda > MAXLAMBDAPOISSON, we do not use pre-computed arrays
      if (lambda > MAXLAMBDA) {
         pdf = null;
         cdf = null;
         return;
      }
 
      double epsilon;
      int i, mid, Nmax;
      int imin, imax;
      double sum;
      double[] P; // Poisson probability terms
      double[] F; // Poisson cumulative probabilities
 
      // In theory, the Poisson distribution has an infinite range. But
      // for i > Nmax, probabilities should be extremely small.
      Nmax = (int) (lambda + 16 * (2 + Math.sqrt(lambda)));
      P = new double[1 + Nmax];
 
      mid = (int) lambda;
      epsilon = EPSILON * EPS_EXTRA / prob(lambda, mid);
      // For large lambda, mass will lose a few digits of precision
      // We shall normalize by explicitly summing all terms >= epsilon
      sum = P[mid] = 1.0;
 
      // Start from the maximum and compute terms > epsilon on each side.
      i = mid;
      while (i > 0 && P[i] > epsilon) {
         P[i - 1] = P[i] * i / lambda;
         i--;
         sum += P[i];
      }
      xmin = imin = i;
 
      i = mid;
      while (P[i] > epsilon) {
         P[i + 1] = P[i] * lambda / (i + 1);
         i++;
         sum += P[i];
         if (i >= Nmax - 1) {
            Nmax *= 2;
            double[] nT = new double[1 + Nmax];
            System.arraycopy(P, 0, nT, 0, P.length);
            P = nT;
         }
      }
      xmax = imax = i;
      F = new double[1 + Nmax];
 
      // Renormalize the sum of probabilities to 1
      for (i = imin; i <= imax; i++)
         P[i] /= sum;
 
      // Compute the cumulative probabilities until F >= 0.5, and keep them in
      // the lower part of array, i.e. F[s] contains all P[i] for i <= s
      F[imin] = P[imin];
      i = imin;
      while (i < imax && F[i] < 0.5) {
         i++;
         F[i] = P[i] + F[i - 1];
      }
      // This is the boundary between F and 1 - F in the CDF
      xmed = i;
 
      // Compute the cumulative probabilities of the complementary distribution
      // and keep them in the upper part of the array. i.e. F[s] contains all
      // P[i] for i >= s
      F[imax] = P[imax];
      i = imax - 1;
      do {
         F[i] = P[i] + F[i + 1];
         i--;
      } while (i > xmed);
 
      /*
       * Reset imin because we lose too much precision for a few terms near imin when
       * we stop adding terms < epsilon.
       */
      i = imin;
      while (i < xmed && F[i] < EPSILON)
         i++;
      xmin = imin = i;
 
      /* Same thing with imax */
      i = imax;
      while (i > xmed && F[i] < EPSILON)
         i--;
      xmax = imax = i;
 
      pdf = new double[imax + 1 - imin];
      cdf = new double[imax + 1 - imin];
      System.arraycopy(P, imin, pdf, 0, imax - imin + 1);
      System.arraycopy(F, imin, cdf, 0, imax - imin + 1);
   }

   public double[] getParams() {
      double[] retour = { lambda };
      return retour;
   }

   public String toString() {
      return getClass().getSimpleName() + ": lambda = " + lambda;
   }
 
}