CramerVonMisesDist.java

/*
 * Class:        CramerVonMisesDist
 * Description:  Cramér-von Mises 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 CramerVonMisesDist extends ContinuousDistribution {
   protected int n;
 
   private static class Function implements MathFunction {
      protected int n;
      protected double u;
 
      public Function(int n, double u) {
         this.n = n;
         this.u = u;
      }
 
      public double evaluate(double x) {
         return u - cdf(n, x);
      }
   }

   public CramerVonMisesDist(int n) {
      setN(n);
   }
 
   public double density(double x) {
      return density(n, x);
   }
 
   public double cdf(double x) {
      return cdf(n, x);
   }
 
   public double barF(double x) {
      return barF(n, x);
   }
 
   public double inverseF(double u) {
      return inverseF(n, u);
   }
 
   public double getMean() {
      return CramerVonMisesDist.getMean(n);
   }
 
   public double getVariance() {
      return CramerVonMisesDist.getVariance(n);
   }
 
   public double getStandardDeviation() {
      return CramerVonMisesDist.getStandardDeviation(n);
   }

   public static double density(int n, double x) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
 
      if (x <= 1.0 / (12.0 * n) || x >= n / 3.0)
         return 0.0;
 
      if (n == 1)
         return 1.0 / Math.sqrt(x - 1.0 / 12.0);
 
      if (x <= 0.002 || x > 3.95)
         return 0.0;
 
      throw new UnsupportedOperationException("density not implemented.");
   }

   public static double cdf(int n, double x) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
 
      if (n == 1) {
         if (x <= 1.0 / 12.0)
            return 0.0;
         if (x >= 1.0 / 3.0)
            return 1.0;
         return 2.0 * Math.sqrt(x - 1.0 / 12.0);
      }
 
      if (x <= 1.0 / (12.0 * n))
         return 0.0;
 
      if (x <= (n + 3.0) / (12.0 * n * n)) {
         double t = Num.lnFactorial(n) - Num.lnGamma(1.0 + 0.5 * n)
               + 0.5 * n * Math.log(Math.PI * (x - 1.0 / (12.0 * n)));
         return Math.exp(t);
      }
 
      if (x <= 0.002)
         return 0.0;
      if (x > 3.95 || x >= n / 3.0)
         return 1.0;
 
      final double EPSILON = Num.DBL_EPSILON;
      final int JMAX = 20;
      int j = 0;
      double Cor, Res, arg;
      double termX, termS, termJ;
 
      termX = 0.0625 / x; // 1 / (16x)
      Res = 0.0;
 
      final double A[] = { 1.0, 1.11803398875, 1.125, 1.12673477358, 1.1274116945, 1.12774323743, 1.1279296875,
            1.12804477649, 1.12812074678, 1.12817350091 };
 
      do {
         termJ = 4 * j + 1;
         arg = termJ * termJ * termX;
         termS = A[j] * Math.exp(-arg) * Num.besselK025(arg);
         Res += termS;
         ++j;
      } while (!(termS < EPSILON || j > JMAX));
 
      if (j > JMAX)
         System.err.println("cramerVonMises: iterations have not converged");
      Res /= Math.PI * Math.sqrt(x);
 
      // Empirical correction in 1/n
      if (x < 0.0092)
         Cor = 0.0;
      else if (x < 0.03)
         Cor = -0.0121763 + x * (2.56672 - 132.571 * x);
      else if (x < 0.06)
         Cor = 0.108688 + x * (-7.14677 + 58.0662 * x);
      else if (x < 0.19)
         Cor = -0.0539444 + x * (-2.22024 + x * (25.0407 - 64.9233 * x));
      else if (x < 0.5)
         Cor = -0.251455 + x * (2.46087 + x * (-8.92836 + x * (14.0988 - x * (5.5204 + 4.61784 * x))));
      else if (x <= 1.1)
         Cor = 0.0782122 + x * (-0.519924 + x * (1.75148 + x * (-2.72035 + x * (1.94487 - 0.524911 * x))));
      else
         Cor = Math.exp(-0.244889 - 4.26506 * x);
 
      Res += Cor / n;
      // This empirical correction is not very precise, so ...
      if (Res <= 1.0)
         return Res;
      else
         return 1.0;
   }

   public static double barF(int n, double x) {
      return 1.0 - cdf(n, x);
   }

   public static double inverseF(int n, double u) {
      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 n / 3.0;
      if (u <= 0.0)
         return 1.0 / (12.0 * n);
 
      if (n == 1)
         return 1.0 / 12.0 + 0.25 * u * u;
 
      Function f = new Function(n, u);
      return RootFinder.brentDekker(0.0, 10.0, f, 1e-6);
   }

   public static double getMean(int n) {
      return 1.0 / 6.0;
   }

   public static double getVariance(int n) {
      return (4.0 * n - 3.0) / (180.0 * n);
   }

   public static double getStandardDeviation(int n) {
      return Math.sqrt(getVariance(n));
   }

   public int getN() {
      return n;
   }

   public void setN(int n) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      this.n = n;
      supportA = 1.0 / (12.0 * n);
      supportB = n / 3.0;
   }

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

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