ChiSquareDist.java

/*
 * Class:        ChiSquareDist
 * Description:  chi-square 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.*;
import umontreal.ssj.functions.MathFunction;

public class ChiSquareDist extends ContinuousDistribution {
   protected int n;
   protected double C1;
 
   private static class Function implements MathFunction {
      protected int n;
      protected double sumLog;
 
      public Function(double s, int n) {
         this.n = n;
         this.sumLog = s;
      }
 
      public double evaluate(double k) {
         if (k < 1.0)
            return 1.0e200;
         return (sumLog + n * (Num.lnGamma(k / 2.0) - 0.5 * Num.LN2 - Num.lnGamma((k + 1.0) / 2.0)));
      }
   }

   public ChiSquareDist(int n) {
      setN(n);
   }
 
   public double density(double x) {
      if (x <= 0)
         return 0.0;
      return Math.exp((n / 2.0 - 1) * Math.log(x) - x / 2.0 - C1);
   }
 
   public double cdf(double x) {
      return cdf(n, decPrec, x);
   }
 
   public double barF(double x) {
      return barF(n, decPrec, x);
   }
 
   public double inverseF(double u) {
      return inverseF(n, u);
   }
 
   public double getMean() {
      return ChiSquareDist.getMean(n);
   }
 
   public double getVariance() {
      return ChiSquareDist.getVariance(n);
   }
 
   public double getStandardDeviation() {
      return ChiSquareDist.getStandardDeviation(n);
   }

   public static double density(int n, double x) {
      if (x <= 0)
         return 0.0;
      return Math.exp((n / 2.0 - 1) * Math.log(x) - x / 2.0 - (n / 2.0) * Num.LN2 - Num.lnGamma(n / 2.0));
   }

   public static double cdf(int n, int d, double x) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (x <= 0.0)
         return 0.0;
      if (x >= XBIG * n)
         return 1.0;
      return GammaDist.cdf(n / 2.0, d, x / 2.0);
   }

   public static double barF(int n, int d, double x) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (x <= 0.0)
         return 1.0;
      return GammaDist.barF(n / 2.0, d, x / 2.0);
   }

   public static double inverseF(int n, double u) {
      /*
       * Returns an approximation of the inverse of Chi square cdf with n degrees of
       * freedom. As in Figure L.23 of P.Bratley, B.L.Fox, and L.E.Schrage. A Guide to
       * Simulation Springer-Verlag, New York, second edition, 1987.
       */
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException("u must be in [0,1]");
      if (u == 1.0)
         return Double.POSITIVE_INFINITY;
      if (u == 0.0)
         return 0.0;
 
      final double E = 0.5e-5; // Precision of this approximation
      final double AA = 0.6931471805;
      double A, XX, X, C, G, CH, Q, P1, P2, T, B, S1, S2, S3, S4, S5, S6;
 
      if (u < 0.00001 || u > 1.0 - 1.0e-5)
         return 2.0 * GammaDist.inverseF(n / 2.0, 7, u);
      if (u >= 1.0)
         return n * XBIG;
      if (u >= 0.999998)
         return (n + 4.0 * Math.sqrt(2.0 * n));
 
      G = Num.lnGamma(n / 2.0);
      XX = 0.5 * n;
      C = XX - 1.0;
 
      if (n >= -1.24 * Math.log(u)) {
         X = NormalDist.inverseF01(u);
         P1 = 0.222222 / n;
         Q = X * Math.sqrt(P1) + 1.0 - P1;
         CH = n * Q * Q * Q;
         if (CH > 2.2 * n + 6.0)
            CH = -2.0 * (Math.log1p(-u) - C * Math.log(0.5 * CH) + G);
 
      } else {
         CH = Math.pow(u * XX * Math.exp(G + XX * AA), 1.0 / XX);
         if (CH - E < 0)
            return CH;
      }
 
      Q = CH;
      P1 = 0.5 * CH;
      P2 = u - GammaDist.cdf(XX, 5, P1);
      if (GammaDist.cdf(XX, 5, P1) == -1.0)
         throw new IllegalArgumentException("RESULT = -1");
 
      T = P2 * Math.exp(XX * AA + G + P1 - C * Math.log(CH));
      B = T / CH;
      A = 0.5 * T - B * C;
      S1 = (210.0 + A * (140.0 + A * (105.0 + A * (84.0 + A * (70.0 + 60.0 * A))))) / 420.0;
      S2 = (420.0 + A * (735.0 + A * (966.0 + A * (1141.0 + 1278.0 * A)))) / 2520.0;
      S3 = (210.0 + A * (462.0 + A * (707.0 + 932.0 * A))) / 2520.0;
      S4 = (252.0 + A * (672.0 + 1182.0 * A) + C * (294.0 + A * (889.0 + 1740.0 * A))) / 5040.0;
      S5 = (84.0 + 264.0 * A + C * (175.0 + 606.0 * A)) / 2520.0;
      S6 = (120.0 + C * (346.0 + 127.0 * C)) / 5040.0;
      CH = CH + T * (1.0 + 0.5 * T * S1 - B * C * (S1 - B * (S2 - B * (S3 - B * (S4 - B * (S5 - B * S6))))));
 
      double temp;
      while (Math.abs(Q / CH - 1.0) > E) {
         Q = CH;
         P1 = 0.5 * CH;
         temp = GammaDist.cdf(XX, 6, P1);
         P2 = u - temp;
 
         if (temp == -1.0)
            return -1.0;
 
         T = P2 * Math.exp(XX * AA + G + P1 - C * Math.log(CH));
         B = T / CH;
         A = 0.5 * T - B * C;
         S1 = (210.0 + A * (140.0 + A * (105.0 + A * (84.0 + A * (70.0 + 60.0 * A))))) / 420.0;
         S2 = (420.0 + A * (735.0 + A * (966.0 + A * (1141.0 + 1278.0 * A)))) / 2520.0;
         S3 = (210.0 + A * (462.0 + A * (707.0 + 932.0 * A))) / 2520.0;
         S4 = (252.0 + A * (672.0 + 1182.0 * A) + C * (294.0 + A * (889.0 + 1740.0 * A))) / 5040.0;
         S5 = (84.0 + 264.0 * A + C * (175.0 + 606.0 * A)) / 2520.0;
         S6 = (120.0 + C * (346.0 + 127.0 * C)) / 5040.0;
         CH = CH + T * (1.0 + 0.5 * T * S1 - B * C * (S1 - B * (S2 - B * (S3 - B * (S4 - B * (S5 - B * S6))))));
      }
      return CH;
   }

   public static double[] getMLE(double[] x, int m) {
      if (m <= 0)
         throw new IllegalArgumentException("m <= 0");
 
      double[] parameters;
 
      parameters = getMomentsEstimate(x, m);
      double k = Math.round(parameters[0]) - 5.0;
      if (k < 1.0)
         k = 1.0;
 
      double sum = 0.0;
      for (int i = 0; i < m; i++) {
         if (x[i] > 0.0)
            sum += 0.5 * Math.log(x[i]);
         else
            sum -= 709.0;
      }
 
      Function f = new Function(sum, m);
      while (f.evaluate(k) > 0.0)
         k++;
      parameters[0] = k;
 
      return parameters;
   }

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

   public static double getMean(int n) {
      if (n <= 0)
         throw new IllegalArgumentException("degrees of freedom " + "must be non-null and positive.");
 
      return n;
   }

   public static double[] getMomentsEstimate(double[] x, int m) {
      double[] parameters = new double[1];
 
      double sum = 0.0;
      for (int i = 0; i < m; i++)
         sum += x[i];
      parameters[0] = sum / (double) m;
 
      return parameters;
   }

   public static double getVariance(int n) {
      if (n <= 0)
         throw new IllegalArgumentException("degrees of freedom " + "must be non-null and positive.");
 
      return (2 * n);
   }

   public static double getStandardDeviation(int n) {
      if (n <= 0)
         throw new IllegalArgumentException("degrees of freedom " + "must be non-null and positive.");
 
      return Math.sqrt(2 * n);
   }

   public int getN() {
      return n;
   }

   public void setN(int n) {
      if (n <= 0)
         throw new IllegalArgumentException("degrees of freedom " + "must be non-null and positive.");
 
      this.n = n;
      supportA = 0.0;
      C1 = 0.5 * n * Num.LN2 + Num.lnGamma(n / 2.0);
   }

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

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