GammaDist.java

/*
 * Class:        GammaDist
 * Description:  gamma 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.Num;
import umontreal.ssj.functions.MathFunction;
import umontreal.ssj.util.RootFinder;

public class GammaDist extends ContinuousDistribution {
   private double alpha; // shape parameter
   private double lambda; // rate parameter
   private double logFactor; // Log (lambda^alpha / Gamma (alpha))
   private static final double ALIM = 1.0E5;
 
   private static class Function implements MathFunction {
      // For MLE
      private int n;
      private double empiricalMean;
      private double sumLn;
 
      public Function(int n, double empiricalMean, double sumLn) {
         this.n = n;
         this.empiricalMean = empiricalMean;
         this.sumLn = sumLn;
      }
 
      public double evaluate(double x) {
         if (x <= 0.0)
            return 1.0e200;
         return (n * Math.log(empiricalMean / x) + n * Num.digamma(x) - sumLn);
      }
   }
 
   private static class myFunc implements MathFunction {
      // For inverseF
      protected int d;
      protected double alp, u;
 
      public myFunc(double alp, int d, double u) {
         this.alp = alp;
         this.d = d;
         this.u = u;
      }
 
      public double evaluate(double x) {
         return u - GammaDist.cdf(alp, d, x);
      }
   }

   public GammaDist(double alpha) {
      setParams(alpha, 1.0, decPrec);
   }

   public GammaDist(double alpha, double lambda) {
      setParams(alpha, lambda, decPrec);
   }

   public GammaDist(double alpha, double lambda, int d) {
      setParams(alpha, lambda, d);
   }
 
   static double mybelog(double x)
   /*
    * This is the function 1 + (1 - x*x + 2*x*log(x)) / ((1 - x)*(1 - x))
    */
   {
      if (x < 1.0e-30)
         return 0.0;
      if (x > 1.0e30)
         return 2.0 * (Math.log(x) - 1.0) / x;
      if (x == 1.0)
         return 1.0;
 
      double t = 1.0 - x;
      if (x < 0.9 || x > 1.1) {
         double w = (t + x * Math.log(x)) / (t * t);
         return 2.0 * w;
      }
 
      // For x near 1, use a series expansion to avoid loss of precision.
      double term;
      final double EPS = 1.0e-12;
      double tpow = 1.0;
      double sum = 0.5;
      int j = 3;
      do {
         tpow *= t;
         term = tpow / (j * (j - 1));
         sum += term;
         j++;
      } while (Math.abs(term / sum) > EPS);
      return 2.0 * sum;
   }
 
   public double density(double x) {
      if (x <= 0)
         return 0.0;
      double z = logFactor + (alpha - 1.0) * Math.log(x) - lambda * x;
      if (z > -XBIGM)
         return Math.exp(z);
      else
         return 0.0;
   }
 
   public double cdf(double x) {
      return cdf(alpha, lambda, decPrec, x);
   }
 
   public double barF(double x) {
      return barF(alpha, lambda, decPrec, x);
   }
 
   public double inverseF(double u) {
      return inverseF(alpha, decPrec, u) / lambda;
   }
 
   public double getMean() {
      return GammaDist.getMean(alpha, lambda);
   }
 
   public double getVariance() {
      return GammaDist.getVariance(alpha, lambda);
   }
 
   public double getStandardDeviation() {
      return GammaDist.getStandardDeviation(alpha, lambda);
   }

   public static double density(double alpha, double lambda, double x) {
      if (alpha <= 0)
         throw new IllegalArgumentException("alpha <= 0");
      if (lambda <= 0)
         throw new IllegalArgumentException("lambda <= 0");
 
      if (x <= 0)
         return 0.0;
      double z = alpha * Math.log(lambda * x) - lambda * x - Num.lnGamma(alpha);
      if (z > -XBIGM)
         return Math.exp(z) / x;
      else
         return 0.0;
   }

   public static double cdf(double alpha, double lambda, int d, double x) {
      return cdf(alpha, d, lambda * x);
   }

   public static double cdf(double alpha, int d, double x) {
      if (alpha <= 0.0)
         throw new IllegalArgumentException("alpha <= 0");
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
      if (x <= 0.0)
         return 0.0;
      if (1.0 == alpha)
         return ExponentialDist.cdf(1.0, x);
 
      if (alpha > 10.0) {
         if (x > alpha * 10.0)
            return 1.0;
      } else {
         if (x > XBIG)
            return 1.0;
      }
 
      if (alpha >= ALIM) {
         double d2 = x + 1.0 / 3.0 - alpha - 0.02 / alpha;
         double S = alpha - 1.0 / 2.0;
         double w = mybelog(S / x);
         double y = d2 * Math.sqrt(w / x);
         return NormalDist.cdf01(y);
      }
 
      if (x <= 1.0 || x < alpha) {
         double factor, z, rn, term;
         factor = Math.exp(alpha * Math.log(x) - x - Num.lnGamma(alpha));
         final double EPS = EPSARRAY[d];
         z = 1.0;
         term = 1.0;
         rn = alpha;
         do {
            rn += 1.0;
            term *= x / rn;
            z += term;
         } while (term >= EPS * z);
         return z * factor / alpha;
 
      } else
         return 1.0 - barF(alpha, d, x);
   }

   public static double barF(double alpha, double lambda, int d, double x) {
      return barF(alpha, d, lambda * x);
   }

   public static double barF(double alpha, int d, double x) {
      if (alpha <= 0.0)
         throw new IllegalArgumentException("alpha <= 0");
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
      if (x <= 0.0)
         return 1.0;
      if (1.0 == alpha)
         return ExponentialDist.barF(1.0, x);
 
      if (alpha >= 70.0) {
         if (x >= alpha * XBIG)
            return 0.0;
      } else {
         if (x >= XBIGM)
            return 0.0;
      }
 
      if (alpha >= ALIM) {
         double d2 = x + 1.0 / 3.0 - alpha - 0.02 / alpha;
         double S = alpha - 1.0 / 2.0;
         double w = mybelog(S / x);
         double y = d2 * Math.sqrt(w / x);
         return NormalDist.barF01(y);
      }
 
      if (x <= 1.0 || x < alpha)
         return 1.0 - cdf(alpha, d, x);
 
      double[] V = new double[6];
      final double EPS = EPSARRAY[d];
      final double RENORM = 1.0E100;
      double R, dif;
      int i;
      double factor = Math.exp(alpha * Math.log(x) - x - Num.lnGamma(alpha));
 
      double A = 1.0 - alpha;
      double B = A + x + 1.0;
      double term = 0.0;
      V[0] = 1.0;
      V[1] = x;
      V[2] = x + 1.0;
      V[3] = x * B;
      double res = V[2] / V[3];
 
      do {
         A += 1.0;
         B += 2.0;
         term += 1.0;
         V[4] = B * V[2] - A * term * V[0];
         V[5] = B * V[3] - A * term * V[1];
         if (V[5] != 0.0) {
            R = V[4] / V[5];
            dif = Math.abs(res - R);
            if (dif <= EPS * R)
               return factor * res;
            res = R;
         }
         for (i = 0; i < 4; i++)
            V[i] = V[i + 2];
         if (Math.abs(V[4]) >= RENORM) {
            for (i = 0; i < 4; i++)
               V[i] /= RENORM;
         }
      } while (true);
   }

   public static double inverseF(double alpha, double lambda, int d, double u) {
      return inverseF(alpha, d, u) / lambda;
   }

   public static double inverseF(double alpha, int d, double u) {
      if (alpha <= 0.0)
         throw new IllegalArgumentException("alpha <= 0");
      if (u > 1.0 || u < 0.0)
         throw new IllegalArgumentException("u not in [0,1]");
      if (u <= 0.0)
         return 0;
      if (u >= 1.0)
         return Double.POSITIVE_INFINITY;
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
      if (d > 15)
         d = 15;
      final double EPS = Math.pow(10.0, -d);
 
      double sigma = GammaDist.getStandardDeviation(alpha, 1.0);
      double x = NormalDist.inverseF(alpha, sigma, u);
      if (x < 0.)
         x = 0.;
      double v = GammaDist.cdf(alpha, d, x);
      double xmax;
      if (alpha < 1.0)
         xmax = 100.0;
      else
         xmax = alpha + 40.0 * sigma;
      myFunc f = new myFunc(alpha, d, u);
 
      if (u <= 1.0e-8 || alpha <= 1.5) {
         if (v < u)
            return RootFinder.bisection(x, xmax, f, EPS);
         else
            return RootFinder.bisection(0, x, f, EPS);
      } else {
         if (v < u)
            return RootFinder.brentDekker(x, xmax, f, EPS);
         else
            return RootFinder.brentDekker(0, x, f, EPS);
      }
   }

   public static double[] getMLE(double[] x, int n) {
      double parameters[];
      double sum = 0.0;
      double sumLn = 0.0;
      double empiricalMean;
      double alphaMME;
      double a;
      final double LN_EPS = Num.LN_DBL_MIN - Num.LN2;
 
      parameters = new double[2];
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      for (int i = 0; i < n; i++) {
         sum += x[i];
         if (x[i] <= 0.0)
            sumLn += LN_EPS;
         else
            sumLn += Math.log(x[i]);
      }
      empiricalMean = sum / (double) n;
 
      sum = 0.0;
      for (int i = 0; i < n; i++) {
         sum += (x[i] - empiricalMean) * (x[i] - empiricalMean);
      }
 
      alphaMME = (empiricalMean * empiricalMean * (double) n) / sum;
      if ((a = alphaMME - 10.0) <= 0) {
         a = 1.0e-5;
      }
 
      Function f = new Function(n, empiricalMean, sumLn);
      parameters[0] = RootFinder.brentDekker(a, alphaMME + 10.0, f, 1e-7);
      parameters[1] = parameters[0] / empiricalMean;
 
      return parameters;
   }

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

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

   public static double getVariance(double alpha, double lambda) {
      if (alpha <= 0.0)
         throw new IllegalArgumentException("alpha <= 0");
      if (lambda <= 0.0)
         throw new IllegalArgumentException("lambda <= 0");
 
      return (alpha / (lambda * lambda));
   }

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

   public double getAlpha() {
      return alpha;
   }

   public double getLambda() {
      return lambda;
   }
 
   public void setParams(double alpha, double lambda, int d) {
      if (alpha <= 0)
         throw new IllegalArgumentException("alpha <= 0");
      if (lambda <= 0)
         throw new IllegalArgumentException("lambda <= 0");
 
      this.alpha = alpha;
      this.lambda = lambda;
      this.decPrec = d;
      logFactor = alpha * Math.log(lambda) - Num.lnGamma(alpha);
      supportA = 0.0;
   }

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

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