DiscreteDistributionInt.java

/*
 * Class:        DiscreteDistributionInt
 * Description:  discrete distributions over the integers
 * 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;

 
public abstract class DiscreteDistributionInt implements Distribution {

   public static double EPSILON = 1.0e-16;
 
   /*
    * For better precision in the tails, we keep the cumulative probabilities (F)
    * in cdf[x] for x <= xmed (i.e. cdf[x] is the sum off all the probabi- lities
    * pdf[i] for i <= x), and the complementary cumulative probabilities (1 - F) in
    * cdf[x] for x > xmed (i.e. cdf[x] is the sum off all the probabilities pdf[i]
    * for i >= x).
    */
   protected final static double EPS_EXTRA = 1.0e-6;
   protected double cdf[] = null; // cumulative probabilities
   protected double pdf[] = null; // probability terms or mass distribution
   protected int xmin = 0; // pdf[x] < EPSILON for x < xmin
   protected int xmax = 0; // pdf[x] < EPSILON for x > xmax
 
   // xmed is such that cdf[xmed] >= 0.5 and cdf[xmed - 1] < 0.5.
   protected int xmed = 0; /*
                            * cdf[x] = F(x) for x <= xmed, and cdf[x] = bar_F(x) for x > xmed
                            */
   protected int supportA = Integer.MIN_VALUE;
   protected int supportB = Integer.MAX_VALUE;

   public abstract double prob(int x);

   public double cdf(double x) {
      return cdf((int) x);
   }

   public abstract double cdf(int x);

   public double barF(double x) {
      return barF((int) x);
   }

   public double barF(int x) {
      return 1.0 - cdf(x - 1);
   }

   public int getXinf() {
      return supportA;
   }

   public int getXsup() {
      return supportB;
   }

   public double inverseF(double u) {
      return inverseFInt(u);
   }

   public int inverseFInt(double u) {
      int i, j, k;
 
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException("u is not in [0,1]");
      if (u <= 0.0)
         return supportA;
      if (u >= 1.0)
         return supportB;
 
      // Remember: the upper part of cdf contains the complementary distribu-
      // tion for xmed < s <= xmax, and the lower part of cdf the
      // distribution for xmin <= x <= xmed
 
      if (u <= cdf[xmed - xmin]) {
         // In the lower part of cdf
         if (u <= cdf[0])
            return xmin;
         i = 0;
         j = xmed - xmin;
         while (i < j) {
            k = (i + j) / 2;
            if (u > cdf[k])
               i = k + 1;
            else
               j = k;
         }
      } else {
         // In the upper part of cdf
         u = 1 - u;
         if (u < cdf[xmax - xmin])
            return xmax;
 
         i = xmed - xmin + 1;
         j = xmax - xmin;
         while (i < j) {
            k = (i + j) / 2;
            if (u < cdf[k])
               i = k + 1;
            else
               j = k;
         }
         i--;
      }
 
      return i + xmin;
   }
 
}