FDist.java

/*
 * Class:        FDist
 * Description:  Empirical distributions
 * 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.gof;

public class FDist {
   private FDist() {
   }

   public static double kolmogorovSmirnovPlusJumpOne(int N, double a, double x) {
      final double EPSILONLR = 1.E-15;
      final double EPSILON = 1.0E-290;
      final double NXAPARAM = 6.5; // frontier: alternating series
      double LogCom;
      double q, p1, q1;
      double Sum = 0.0;
      double term;
      double Njreal;
      double jreal;
      int Sign;
      int j;
      int jmax;
 
      if (N < 1)
         throw new IllegalArgumentException("Calling kolmogorovSmirnovPlusJumpOne " + "with N < 1");
      if (a >= 1.0 || a <= 0.0)
         throw new IllegalArgumentException("Calling kolmogorovSmirnovPlusJumpOne " + "with a outside (0, 1)");
      if (x <= 0.0)
         return 0.0;
      if (x + a >= 1.0)
         return 1.0;
      LogCom = Math.log((double) N);
 
      // --------------------------------------------------------------------
      // the alternating series is stable and fast for N*(x + a) very small
      // --------------------------------------------------------------------
      if (N * (x + a) < NXAPARAM && a + x < 0.5) {
         jmax = (int) (N * (x + a));
         for (j = 1; j <= jmax; j++) {
            jreal = j;
            Njreal = N - j;
            q = jreal / N - x;
            if ((q < 0.0 && ((j & 1) != 0)) || ((q > 1.0) && (((N - j - 1) & 1) != 0)))
               Sign = -1;
            else
               Sign = 1;
 
            // we must avoid log (0.0)
            q1 = Math.abs(q);
            p1 = Math.abs(1.0 - q);
            if (q1 > EPSILON && p1 > EPSILON) {
               term = LogCom + jreal * Math.log(q1) + (Njreal - 1.0) * Math.log(p1);
               Sum += Sign * Math.exp(term);
            }
            LogCom += Math.log(Njreal / (jreal + 1.0));
         }
         // add the term j = 0
         Sum += Math.exp((N - 1) * Math.log(1.0 + x));
         return Sum * x;
      }
 
      // ---------------------------------------------
      // For N (x + a) >= NxaParam or (a + x) > 0.5,
      // use the non-alternating series.
      // ---------------------------------------------
 
      // EpsilonLR because the distribution has a jump
      jmax = (int) (N * (1.0 - a - x - EPSILONLR));
      for (j = 1; j <= jmax; j++) {
         jreal = j;
         Njreal = N - jreal;
         q = jreal / N + x;
         if (1.0 - q > EPSILON) {
            term = LogCom + (jreal - 1.0) * Math.log(q) + Njreal * Math.log(1.0 - q);
            Sum += Math.exp(term);
         }
         LogCom += Math.log(Njreal / (jreal + 1.0));
      }
      Sum *= x;
 
      // add the term j = 0
      if (1.0 - x > EPSILON)
         Sum += Math.exp(N * Math.log(1.0 - x));
      return 1.0 - Sum;
   }

   public static double scan(int N, double d, int m) {
      return 1.0 - FBar.scan(N, d, m);
   }
}