HypoExponentialDistQuick.java

/*
 * Class:        HypoExponentialDistQuick
 * Description:  Hypo-exponential distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Richard Simard
 * @since        January 2011
 *
 *
 * 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 java.util.Formatter;
import java.util.Locale;
import umontreal.ssj.util.*;
import umontreal.ssj.functions.MathFunction;

public class HypoExponentialDistQuick extends HypoExponentialDist {
   private double[] m_H;
 
   private static double[] computeH(double[] lambda) {
      int k = lambda.length;
      double[] H = new double[k];
      double tem;
      int j;
      for (int i = 0; i < k; i++) {
         tem = 1.0;
         for (j = 0; j < i; j++)
            tem *= lambda[j] / (lambda[j] - lambda[i]);
         for (j = i + 1; j < k; j++)
            tem *= lambda[j] / (lambda[j] - lambda[i]);
         H[i] = tem;
      }
      return H;
   }
 
   private static double m_density(double[] lambda, double[] H, double x) {
      if (x < 0)
         return 0;
 
      int k = lambda.length;
      double tem;
      double prob = 0;
      for (int j = 0; j < k; j++) {
         tem = Math.exp(-lambda[j] * x);
         if (tem > 0)
            prob += lambda[j] * H[j] * tem;
      }
 
      return prob;
   }
 
   private static double m_barF(double[] lambda, double[] H, double x) {
      if (x <= 0.)
         return 1.;
 
      int k = lambda.length;
      double tem;
      double prob = 0; // probability
      for (int j = 0; j < k; j++) {
         tem = Math.exp(-lambda[j] * x);
         if (tem > 0)
            prob += H[j] * tem;
      }
      return prob;
   }
 
   private static double m_cdf(double[] lambda, double[] H, double x) {
      if (x <= 0.)
         return 0.;
 
      int k = lambda.length;
      double tem = Math.exp(-lambda[0] * x);
      if (tem <= 0)
         return 1.0 - HypoExponentialDistQuick.m_barF(lambda, H, x);
 
      double prob = 0; // cumulative probability
      for (int j = 0; j < k; j++) {
         tem = Math.expm1(-lambda[j] * x);
         prob += H[j] * tem;
      }
      return -prob;
   }
 
   private static class myFunc implements MathFunction {
      // For inverseF
      private double[] m_lam;
      private double m_u;
 
      public myFunc(double[] lam, double u) {
         m_lam = lam;
         m_u = u;
      }
 
      public double evaluate(double x) {
         return m_u - HypoExponentialDistQuick.cdf(m_lam, x);
      }
   }

   public HypoExponentialDistQuick(double[] lambda) {
      super(lambda);
   }
 
   /*
    * These public methods are necessary so that methods cdf, barF and inverseF
    * used are those of the present class and not those of the mother class.
    */
   public double density(double x) {
      return m_density(m_lambda, m_H, x);
   }
 
   public double cdf(double x) {
      return m_cdf(m_lambda, m_H, x);
   }
 
   public double barF(double x) {
      return m_barF(m_lambda, m_H, x);
   }
 
   public double inverseF(double u) {
      return m_inverseF(m_lambda, m_H, u);
   }

   public static double density(double[] lambda, double x) {
      testLambda(lambda);
      double[] H = computeH(lambda);
      return m_density(lambda, H, x);
   }

   public static double cdf(double[] lambda, double x) {
      testLambda(lambda);
      double[] H = computeH(lambda);
      return m_cdf(lambda, H, x);
   }

   public static double barF(double[] lambda, double x) {
      testLambda(lambda);
      double[] H = computeH(lambda);
      return m_barF(lambda, H, x);
   }

   public static double inverseF(double[] lambda, double u) {
      testLambda(lambda);
      double[] H = computeH(lambda);
      return m_inverseF(lambda, H, u);
   }
 
   private static double m_inverseF(double[] lambda, double[] H, double u) {
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException("u not in [0,1]");
      if (u >= 1.0)
         return Double.POSITIVE_INFINITY;
      if (u <= 0.0)
         return 0.0;
 
      final double EPS = 1.0e-12;
      myFunc fonc = new myFunc(lambda, u);
      double x1 = getMean(lambda);
      double v = m_cdf(lambda, H, x1);
      if (u <= v)
         return RootFinder.brentDekker(0, x1, fonc, EPS);
 
      // u > v
      double x2 = 4.0 * x1 + 1.0;
      v = m_cdf(lambda, H, x2);
      while (v < u) {
         x1 = x2;
         x2 = 4.0 * x2;
         v = m_cdf(lambda, H, x2);
      }
      return RootFinder.brentDekker(x1, x2, fonc, EPS);
   }
 
   public void setLambda(double[] lambda) {
      if (lambda == null)
         return;
      super.setLambda(lambda);
      m_H = computeH(lambda);
   }
 
   public String toString() {
      StringBuilder sb = new StringBuilder();
      Formatter formatter = new Formatter(sb, Locale.US);
      formatter.format(getClass().getSimpleName() + " : lambda = {" + PrintfFormat.NEWLINE);
      int k = m_lambda.length;
      for (int i = 0; i < k; i++) {
         formatter.format("   %f%n", m_lambda[i]);
      }
      formatter.format("}%n");
      return sb.toString();
   }
}