GofStat.java

/*
 * Class:        GofStat
 * Description:  Goodness-of-fit test statistics
 * 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;
 
import cern.colt.list.*;
import umontreal.ssj.util.*;
import umontreal.ssj.probdist.*;
import java.util.Arrays;

public class GofStat {
   private GofStat() {
   }

 
   // Used in discontinuous distributions
   private static double EPSILOND = 1.0E-15;

   public static DoubleArrayList unifTransform(DoubleArrayList data, ContinuousDistribution dist) {
      double[] v = data.elements();
      int n = data.size();
 
      double[] u = new double[n];
      for (int i = 0; i < n; i++)
         u[i] = dist.cdf(v[i]);
      return new DoubleArrayList(u);
   }

   public static DoubleArrayList unifTransform(DoubleArrayList data, DiscreteDistribution dist) {
      double[] v = data.elements();
      int n = data.size();
 
      double[] u = new double[n];
      for (int i = 0; i < n; i++)
         u[i] = dist.cdf((int) v[i]);
      return new DoubleArrayList(u);
   }

   public static void diff(IntArrayList sortedData, IntArrayList spacings, int n1, int n2, int a, int b) {
      if (n1 < 0 || n2 < 0 || n1 >= n2 || n2 >= sortedData.size())
         throw new IllegalArgumentException("n1 and n2 not valid.");
      int[] u = sortedData.elements();
      int n = sortedData.size();
      if (spacings.size() <= (n2 + 2))
         spacings.setSize(n2 + 2);
      int[] d = spacings.elements();
 
      d[n1] = u[n1] - a;
      for (int i = n1 + 1; i <= n2; i++)
         d[i] = u[i] - u[i - 1];
      d[n2 + 1] = b - u[n2];
   }

   public static void diff(DoubleArrayList sortedData, DoubleArrayList spacings, int n1, int n2, double a, double b) {
 
      if (n1 < 0 || n2 < 0 || n1 >= n2 || n2 >= sortedData.size())
         throw new IllegalArgumentException("n1 and n2 not valid.");
      double[] u = sortedData.elements();
      int n = sortedData.size();
      if (spacings.size() <= (n2 + 2))
         spacings.setSize(n2 + 2);
      double[] d = spacings.elements();
 
      d[n1] = u[n1] - a;
      for (int i = n1 + 1; i <= n2; i++)
         d[i] = u[i] - u[i - 1];
      d[n2 + 1] = b - u[n2];
   }

   public static void iterateSpacings(DoubleArrayList data, DoubleArrayList spacings) {
      if (spacings.size() < (data.size() + 1))
         throw new IllegalArgumentException("Invalid array sizes.");
      double[] v = data.elements();
      spacings.quickSortFromTo(0, data.size());
      double[] s = spacings.elements();
      int n = data.size();
 
      for (int i = 0; i < n; i++)
         s[n - i] = (i + 1) * (s[n - i] - s[n - i - 1]);
      s[0] = (n + 1) * s[0];
      v[0] = s[0];
      for (int i = 1; i < n; i++)
         v[i] = v[i - 1] + s[i];
   }

   public static void powerRatios(DoubleArrayList sortedData) {
 
      double[] u = sortedData.elements();
      int n = sortedData.size();
 
      for (int i = 0; i < (n - 1); i++) {
         if (u[i + 1] == 0.0 || u[i + 1] == -0.0)
            u[i] = 1.0;
         else
            u[i] = Math.pow(u[i] / u[i + 1], (double) i + 1);
      }
 
      u[n - 1] = Math.pow(u[n - 1], (double) n);
      sortedData.quickSortFromTo(0, sortedData.size() - 1);
   }



   public static class OutcomeCategoriesChi2 {

      public int nbCategories;

      public int smin;

      public int smax;

      public double[] nbExp;

      public int[] loc;

      public OutcomeCategoriesChi2(double[] nbExp) {
         this.nbExp = nbExp;
         smin = 0;
         smax = nbExp.length - 1;
         nbCategories = nbExp.length;
         loc = new int[nbExp.length];
         for (int i = 0; i < nbExp.length; i++)
            loc[i] = i;
      }

      public OutcomeCategoriesChi2(double[] nbExp, int smin, int smax) {
         this.nbExp = nbExp;
         this.smin = smin;
         this.smax = smax;
         nbCategories = smax - smin + 1;
         loc = new int[nbExp.length];
         for (int i = 0; i < smin; i++)
            loc[i] = smin;
         for (int i = smin; i < smax; i++)
            loc[i] = i;
         for (int i = smax; i < nbExp.length; i++)
            loc[i] = smax;
      }

      public OutcomeCategoriesChi2(double[] nbExp, int[] loc, int smin, int smax, int nbCat) {
         this.nbExp = nbExp;
         this.smin = smin;
         this.smax = smax;
         this.nbCategories = nbCat;
         this.loc = loc;
      }

      public void regroupCategories(double minExp) {
         int s0 = 0, j;
         double somme;
 
         nbCategories = 0;
         int s = smin;
         while (s <= smax) {
            /*
             * Merge categories to ensure that the number expected in each category is >=
             * minExp.
             */
            if (nbExp[s] < minExp) {
               s0 = s;
               somme = nbExp[s];
               while (somme < minExp && s < smax) {
                  nbExp[s] = 0.0;
                  ++s;
                  somme += nbExp[s];
               }
               nbExp[s] = somme;
               for (j = s0; j <= s; j++)
                  loc[j] = s;
 
            } else
               loc[s] = s;
 
            ++nbCategories;
            ++s;
         }
         smin = loc[smin];
 
         // Special case: the last category, if nbExp < minExp
         if (nbExp[smax] < minExp) {
            if (s0 > smin)
               --s0;
            nbExp[s0] += nbExp[smax];
            nbExp[smax] = 0.0;
            --nbCategories;
            for (j = s0 + 1; j <= smax; j++)
               loc[j] = s0;
            smax = s0;
         }
         if (nbCategories <= 1)
            throw new IllegalStateException("nbCategories < 2");
      }

      public String toString() {
         int s, s0;
         double somme;
         final double EPSILON = 5.0E-16;
         StringBuffer sb = new StringBuffer();
         sb.append("-----------------------------------------------" + PrintfFormat.NEWLINE);
         if (nbExp[smin] < EPSILON)
            sb.append("Only expected numbers larger than " + PrintfFormat.g(6, 1, EPSILON) + "  are printed"
                  + PrintfFormat.NEWLINE);
         sb.append("Number of categories: " + PrintfFormat.d(4, nbCategories) + PrintfFormat.NEWLINE
               + "Expected numbers per category:" + PrintfFormat.NEWLINE + PrintfFormat.NEWLINE
               + "Category s      nbExp[s]" + PrintfFormat.NEWLINE);
 
         // Do not print values < EPSILON
         s = smin;
         while (nbExp[s] < EPSILON)
            s++;
         int s1 = s;
 
         s = smax;
         while (nbExp[s] < EPSILON)
            s--;
         int s2 = s;
 
         somme = 0.0;
         for (s = s1; s <= s2; s++)
            if (loc[s] == s) {
               somme += nbExp[s];
               sb.append(PrintfFormat.d(4, s) + " " + PrintfFormat.f(18, 4, nbExp[s]) + PrintfFormat.NEWLINE);
            }
         sb.append(PrintfFormat.NEWLINE + "Total expected number = " + PrintfFormat.f(18, 2, somme)
               + PrintfFormat.NEWLINE + PrintfFormat.NEWLINE + "The groupings:" + PrintfFormat.NEWLINE
               + " Category s      loc[s]" + PrintfFormat.NEWLINE);
         for (s = smin; s <= smax; s++) {
            if ((s == smin) && (s > 0))
               sb.append("<= ");
            else if ((s == smax) && (s < loc.length - 1))
               sb.append(">= ");
            else
               sb.append("   ");
            sb.append(PrintfFormat.d(4, s) + " " + PrintfFormat.d(12, loc[s]) + PrintfFormat.NEWLINE);
         }
 
         sb.append(PrintfFormat.NEWLINE + PrintfFormat.NEWLINE);
         return sb.toString();
      }
   }



   public static double chi2(double[] nbExp, int[] count, int smin, int smax) {
      double diff, khi = 0.0;
 
      for (int s = smin; s <= smax; s++) {
         if (nbExp[s] <= 0.0) {
            if (count[s] != 0)
               throw new IllegalArgumentException("nbExp[s] = 0 and count[s] > 0");
         } else {
            diff = count[s] - nbExp[s];
            khi += diff * diff / nbExp[s];
         }
      }
      return khi;
   }

   public static double chi2(OutcomeCategoriesChi2 cat, int[] count) {
      int[] newcount = new int[1 + cat.smax];
      for (int s = cat.smin; s <= cat.smax; s++) {
         newcount[cat.loc[s]] += count[s];
      }
 
      double diff, khi = 0.0;
 
      for (int s = cat.smin; s <= cat.smax; s++) {
         if (cat.nbExp[s] > 0.0) {
            diff = newcount[s] - cat.nbExp[s];
            khi += diff * diff / cat.nbExp[s];
         }
      }
      newcount = null;
      return khi;
   }

   public static double chi2(IntArrayList data, DiscreteDistributionInt dist, int smin, int smax, double minExp,
         int[] numCat) {
      int i;
      int n = data.size();
 
      // Find the first non-negligible probability term and fix
      // the real smin. The linear search starts from the given smin.
      i = smin;
      while (dist.prob(i) * n <= DiscreteDistributionInt.EPSILON)
         i++;
      smin = i--;
 
      // smax > smin is required
      while (smax <= smin)
         smax = 2 * smax + 1;
 
      // Allocate and fill the array of expected observations
      // Each category s corresponds to a value s for which p(s)>0.
      double[] nbExp = new double[smax + 1];
      do {
         i++;
         if (i > smax) {
            smax *= 2;
            double[] newNbExp = new double[smax + 1];
            System.arraycopy(nbExp, smin, newNbExp, smin, nbExp.length - smin);
            nbExp = newNbExp;
         }
         nbExp[i] = dist.prob(i) * n;
      } while (nbExp[i] > DiscreteDistributionInt.EPSILON);
      smax = i - 1;
 
      // Regroup the expected observations intervals
      // satisfying np(s)>=minExp
      OutcomeCategoriesChi2 cat = new OutcomeCategoriesChi2(nbExp, smin, smax);
      cat.regroupCategories(minExp);
      if (numCat != null)
         numCat[0] = cat.nbCategories;
 
      // Count the number of observations in each categories.
      int[] count = new int[cat.smax + 1];
      for (i = 0; i < count.length; i++)
         count[i] = 0;
      for (i = 0; i < n; i++) {
         int s = data.get(i);
         while (cat.loc[s] != s)
            s = cat.loc[s];
         count[s]++;
      }
 
      // Perform the chi-square test
      return chi2(cat.nbExp, count, cat.smin, cat.smax);
   }

   public static double chi2Equal(double nbExp, int[] count, int smin, int smax) {
 
      double diff, khi = 0.0;
      for (int s = smin; s <= smax; s++) {
         diff = count[s] - nbExp;
         khi += diff * diff;
      }
      return khi / nbExp;
   }

   public static double chi2Equal(DoubleArrayList data, double minExp) {
      int n = data.size();
      if (n < (int) Math.ceil(minExp))
         throw new IllegalArgumentException("Not enough observations");
      double p = minExp / n;
      int m = (int) Math.ceil(1.0 / p);
      // to avoid an exception when data[i] = 1/p, reserve one element more
      int[] count = new int[m + 1];
      for (int i = 0; i < n; i++) {
         int j = (int) Math.floor(data.get(i) / p);
         count[j]++;
      }
      // put the elements in count[m] where they belong: in count[m-1]
      count[m - 1] += count[m];
      return chi2Equal(minExp, count, 0, m - 1);
   }

   public static double chi2Equal(DoubleArrayList data) {
      return chi2Equal(data, 10.0);
   }

   public static int scan(DoubleArrayList sortedData, double d) {
 
      double[] u = sortedData.elements();
      int n = sortedData.size();
 
      int m = 1, j = 0, i = -1;
      double High = 0.0;
 
      while (j < (n - 1) && High < 1.0) {
         ++i;
 
         High = u[i] + d;
         while (j < n && u[j] < High)
            ++j;
         // j is now the index of the first obs. to the right of High.
         if (j - i > m)
            m = j - i;
      }
      return m;
   }

   public static double cramerVonMises(DoubleArrayList sortedData) {
      double w, w2;
      double[] u = sortedData.elements();
      int n = sortedData.size();
 
      if (n <= 0) {
         System.err.println("cramerVonMises:  n <= 0");
         return 0.0;
      }
 
      w2 = 1.0 / (12 * n);
      for (int i = 0; i < n; i++) {
         w = u[i] - (i + 0.5) / n;
         w2 += w * w;
      }
      return w2;
   }

   public static double watsonG(DoubleArrayList sortedData) {
      double[] u = sortedData.elements();
      int n = sortedData.size();
      double sumZ;
      double d2;
      double dp, g;
      double unSurN = 1.0 / n;
 
      if (n <= 0) {
         System.err.println("watsonG: n <= 0");
         return 0.0;
      }
 
      // degenerate case n = 1
      if (n == 1)
         return 0.0;
 
      // We assume that u is already sorted.
      dp = sumZ = 0.0;
      for (int i = 0; i < n; i++) {
         d2 = (i + 1) * unSurN - u[i];
         if (d2 > dp)
            dp = d2;
         sumZ += u[i];
      }
      sumZ = sumZ * unSurN - 0.5;
      g = Math.sqrt((double) n) * (dp + sumZ);
      return g;
   }

   public static double watsonU(DoubleArrayList sortedData) {
      double sumZ, w, w2, u2;
      double[] u = sortedData.elements();
      int n = sortedData.size();
 
      if (n <= 0) {
         System.err.println("watsonU: n <= 0");
         return 0.0;
      }
 
      // degenerate case n = 1
      if (n == 1)
         return 1.0 / 12.0;
 
      sumZ = 0.0;
      w2 = 1.0 / (12 * n);
      for (int i = 0; i < n; i++) {
         sumZ += u[i];
         w = u[i] - (i + 0.5) / n;
         w2 += w * w;
      }
      sumZ = sumZ / n - 0.5;
      u2 = w2 - sumZ * sumZ * n;
      return u2;
   }

   public static double EPSILONAD = Num.DBL_EPSILON / 2;

   public static double andersonDarling(DoubleArrayList sortedData) {
      double[] v = sortedData.elements();
      return andersonDarling(v);
   }

   public static double andersonDarling(double[] sortedData) {
      double u1;
      double u, a2;
      int n = sortedData.length;
 
      if (n <= 0) {
         System.err.println("andersonDarling: n <= 0");
         return 0.0;
      }
 
      a2 = 0.0;
      for (int i = 0; i < n; i++) {
         u = sortedData[i];
         u1 = 1.0 - u;
         if (u < EPSILONAD)
            u = EPSILONAD;
         else if (u1 < EPSILONAD)
            u1 = EPSILONAD;
         a2 += (2 * i + 1) * Math.log(u) + (1 + 2 * (n - i - 1)) * Math.log(u1);
      }
      a2 = -n - a2 / n;
      return a2;
   }

   public static double[] andersonDarling(double[] data, ContinuousDistribution dist) {
      int n = data.length;
      double[] U = new double[n];
      for (int i = 0; i < n; i++) {
         U[i] = dist.cdf(data[i]);
      }
 
      Arrays.sort(U);
      double x = GofStat.andersonDarling(U);
      double v = AndersonDarlingDistQuick.barF(n, x);
      double[] res = { x, v };
      return res;
   }

   public static double[] kolmogorovSmirnov(double[] sortedData) {
      DoubleArrayList v = new DoubleArrayList(sortedData);
      return kolmogorovSmirnov(v);
   }

   public static double[] kolmogorovSmirnov(DoubleArrayList sortedData) {
      double[] ret = new double[3];
      int n = sortedData.size();
 
      if (n <= 0) {
         ret[0] = ret[1] = ret[2] = 0.0;
         System.err.println("kolmogorovSmirnov:   n <= 0");
         return ret;
      }
 
      double[] retjo = kolmogorovSmirnovJumpOne(sortedData, 0.0);
      ret[0] = retjo[0];
      ret[1] = retjo[1];
      if (ret[1] > ret[0])
         ret[2] = ret[1];
      else
         ret[2] = ret[0];
 
      return ret;
   }

   public static void kolmogorovSmirnov(double[] data, ContinuousDistribution dist, double[] sval, double[] pval) {
      int n = data.length;
      double[] T = new double[n];
      for (int i = 0; i < n; i++) {
         T[i] = dist.cdf(data[i]);
      }
 
      Arrays.sort(T);
      double[] statks = GofStat.kolmogorovSmirnov(T);
      for (int i = 0; i < 3; i++) {
         sval[i] = statks[i];
      }
      pval[2] = KolmogorovSmirnovDistQuick.barF(n, sval[2]);
      pval[1] = KolmogorovSmirnovPlusDist.barF(n, sval[1]);
      pval[0] = KolmogorovSmirnovPlusDist.barF(n, sval[0]);
   }

   public static double[] kolmogorovSmirnovJumpOne(DoubleArrayList sortedData, double a) {
      /*
       * Statistics KS+ and KS-. Case with 1 jump at a, near the lower tail of the
       * distribution.
       */
 
      double[] u = sortedData.elements();
      int n = sortedData.size();
      int j, i;
      double d2, d1, unSurN;
      double[] ret = new double[2];
 
      if (n <= 0) {
         ret[0] = ret[1] = 0.0;
         System.err.println("kolmogorovSmirnovJumpOne: n <= 0");
         return ret;
      }
 
      ret[0] = 0.0;
      ret[1] = 0.0;
      unSurN = 1.0 / n;
      j = 0;
 
      while (j < n && u[j] <= a + EPSILOND)
         ++j;
 
      for (i = j - 1; i < n; i++) {
         if (i >= 0) {
            d1 = (i + 1) * unSurN - u[i];
            if (d1 > ret[0])
               ret[0] = d1;
         }
         if (i >= j) {
            d2 = u[i] - i * unSurN;
            if (d2 > ret[1])
               ret[1] = d2;
         }
      }
      return ret;
   }

   public static double pDisc(double pL, double pR) {
      double p;
 
      if (pR < pL)
         p = pR;
      else if (pL > 0.5)
         p = 0.5;
      else
         p = 1.0 - pL;
      // Note: si p est tres proche de 1, on perd toute la precision ici!
      // Note2: je ne pense pas que cela puisse se produire a cause des if (RS)
      return p;
   }
}