HypoExponentialDist.java

/*
 * Class:        HypoExponentialDist
 * 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 cern.colt.matrix.*;
import cern.colt.matrix.impl.*;
import umontreal.ssj.util.*;
import umontreal.ssj.functions.MathFunction;

public class HypoExponentialDist extends ContinuousDistribution {
   protected double[] m_lambda;
 
   protected static void testLambda(double[] lambda) {
      int m = lambda.length;
      for (int j = 0; j < m; ++j) {
         if (lambda[j] <= 0)
            throw new IllegalArgumentException("lambda_j <= 0");
      }
   }
 
   // Builds the bidiagonal matrix A out of the lambda
   private static DoubleMatrix2D buildMatrix(double[] lambda, double x) {
      int m = lambda.length;
      testLambda(lambda);
      DoubleFactory2D F2 = DoubleFactory2D.dense;
      DoubleMatrix2D A = F2.make(m, m);
      for (int j = 0; j < m - 1; j++) {
         A.setQuick(j, j, -lambda[j] * x);
         A.setQuick(j, j + 1, lambda[j] * x);
      }
      A.setQuick(m - 1, m - 1, -lambda[m - 1] * x);
      return A;
   }
 
   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 - HypoExponentialDist.cdf(m_lam, x);
      }
   }

   public HypoExponentialDist(double[] lambda) {
      supportA = 0.0;
      setLambda(lambda);
   }
 
   public double density(double x) {
      return density(m_lambda, x);
   }
 
   public double cdf(double x) {
      return cdf(m_lambda, x);
   }
 
   public double barF(double x) {
      return barF(m_lambda, x);
   }
 
   public double inverseF(double u) {
      return inverseF(m_lambda, u);
   }
 
   public double getMean() {
      return getMean(m_lambda);
   }
 
   public double getVariance() {
      return getVariance(m_lambda);
   }
 
   public double getStandardDeviation() {
      return getStandardDeviation(m_lambda);
   }

   public static double density(double[] lambda, double x) {
      if (x < 0)
         return 0;
      DoubleMatrix2D Ax = buildMatrix(lambda, x);
      DoubleMatrix2D T = DMatrix.expBidiagonal(Ax);
      int m = lambda.length;
      return lambda[m - 1] * T.getQuick(0, m - 1);
   }

   public static double cdf(double[] lambda, double x) {
      if (x <= 0.0)
         return 0.0;
      if (x >= Double.MAX_VALUE)
         return 1.0;
      double mean = HypoExponentialDist.getMean(lambda);
      double std = HypoExponentialDist.getStandardDeviation(lambda);
      double LOW = mean - 1.5 * std;
      if (x > LOW) {
         double p = 1.0 - HypoExponentialDist.barF(lambda, x);
         double LIMIT = 1.0e-3;
         if (p > LIMIT)
            return p;
      }
 
      DoubleMatrix2D T = buildMatrix(lambda, x);
      DoubleFactory1D fac1 = DoubleFactory1D.dense;
      int m = lambda.length;
      DoubleMatrix1D C = fac1.make(m, 1.0);
      DoubleMatrix1D B = DMatrix.expmiBidiagonal(T, C);
      return Math.abs(B.getQuick(0));
   }

   public static double cdf2(double[] lambda, double x) {
      if (x <= 0.0)
         return 0.0;
      if (x >= Double.MAX_VALUE)
         return 1.0;
      return (1.0 - HypoExponentialDist.barF(lambda, x));
   }

   public static double barF(double[] lambda, double x) {
      if (x <= 0.0)
         return 1.0;
      if (x >= Double.MAX_VALUE)
         return 0.0;
      DoubleMatrix2D T = buildMatrix(lambda, x);
      DoubleFactory1D fac1 = DoubleFactory1D.dense;
      int m = lambda.length;
      DoubleMatrix1D C = fac1.make(m, 1.0);
      DoubleMatrix1D B = DMatrix.expBidiagonal(T, C);
      return B.getQuick(0);
   }

   public static double inverseF(double[] lambda, 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 = cdf(lambda, x1);
      if (u <= v)
         return RootFinder.brentDekker(0, x1, fonc, EPS);
 
      // u > v
      double x2 = 4.0 * x1 + 1.0;
      v = cdf(lambda, x2);
      while (v < u) {
         x1 = x2;
         x2 = 4.0 * x2;
         v = cdf(lambda, x2);
      }
      return RootFinder.brentDekker(x1, x2, fonc, EPS);
   }

   public static double getMean(double[] lambda) {
      testLambda(lambda);
      int k = lambda.length;
      double sum = 0;
      for (int j = 0; j < k; j++)
         sum += 1.0 / lambda[j];
      return sum;
   }

   public static double getVariance(double[] lambda) {
      testLambda(lambda);
      int k = lambda.length;
      double sum = 0;
      for (int j = 0; j < k; j++)
         sum += 1.0 / (lambda[j] * lambda[j]);
      return sum;
   }

   public static double getStandardDeviation(double[] lambda) {
      double s = getVariance(lambda);
      return Math.sqrt(s);
   }

   public double[] getLambda() {
      return m_lambda;
   }

   public void setLambda(double[] lambda) {
      if (lambda == null)
         return;
      int k = lambda.length;
      m_lambda = new double[k];
      testLambda(lambda);
      System.arraycopy(lambda, 0, m_lambda, 0, k);
   }

   public double[] getParams() {
      return m_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("   %g%n", m_lambda[i]);
      }
      formatter.format("}%n");
      return sb.toString();
   }
 
}