package umontreal.iro.lecuyer.gof;
import cern.colt.list.*;
import umontreal.iro.lecuyer.util.*;
import umontreal.iro.lecuyer.probdist.*;
import java.util.Arrays;
/**
* This class provides methods to compute several types of EDF goodness-of-fit
* test statistics and to apply certain transformations to a set of
* observations. This includes the probability integral transformation
*
* Ui = F(Xi), as well as the power ratio and iterated spacings
* transformations. Here,
* U(0),..., U(N-1) stand
* for N observations
* U0,..., UN-1 sorted by increasing order, where
*
* 0 <= Ui <= 1.
*
*
* Note: This class uses the Colt library.
*
*/
public class GofStat {
private GofStat() {}
// Used in discontinuous distributions
private static double EPSILOND = 1.0E-15;
/**
* Applies the transformation
* Ui = F(Vi) for
* i = 0, 1,…, N - 1,
* where F is a continuous distribution function,
* and returns the result as an array of length N.
* V represents the N observations contained in data,
* and U, the returned transformed observations.
* If data contains random variables from the distribution function
* dist, then the result will contain uniform random variables
* over [0, 1].
*
* @param data array of observations to be transformed
*
* @param dist assumed distribution of the observations
*
* @return the array of transformed observations
*
*/
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);
}
/**
* Applies the transformation
* Ui = F(Vi) for
* i = 0, 1,…, N - 1,
* where F is a discrete distribution function,
* and returns the result as an array of length N.
* V represents the N observations contained in data,
* and U, the returned transformed observations.
*
*
* Note: If V are the values of random variables with
* distribution function dist, then the result will contain
* the values of discrete random variables distributed over the
* set of values taken by dist,
* not uniform random variables over [0, 1].
*
* @param data array of observations to be transformed
*
* @param dist assumed distribution of the observations
*
* @return the array of transformed observations
*
*/
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);
}
/**
* Assumes that the real-valued observations
* U0,..., UN-1
* contained in sortedData
* are already sorted in increasing order and computes the differences
* between the successive observations. Let D be the differences
* returned in spacings.
* The difference
* Ui - Ui-1 is put in Di for
* n1 < i <= n2, whereas
* Un1 - a is put into Dn1
* and
* b - Un2 is put into Dn2+1.
* The number of observations must be greater or equal than n2, we
* must have
* n1 < n2, and n1 and n2 are greater than 0.
* The size of spacings will be at least N + 1 after
* the call returns.
*
* @param sortedData array of sorted observations
*
* @param spacings pointer to an array object that will be filled with spacings
*
* @param n1 starting index, in sortedData, of the processed observations
*
* @param n2 ending index, in sortedData of the processed observations
*
* @param a minimum value of the observations
*
* @param b maximum value of the observations
*
*
*/
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];
}
/**
* Same as method
* {@link #diff diff}(IntArrayList,IntArrayList,int,int,int,int),
* but for the continuous case.
*
* @param sortedData array of sorted observations
*
* @param spacings pointer to an array object that will be filled with spacings
*
* @param n1 starting index of the processed observations in sortedData
*
* @param n2 ending index, in sortedData of the processed observations
*
* @param a minimum value of the observations
*
* @param b maximum value of the observations
*
*
*/
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];
}
/**
* Applies one iteration of the iterated spacings
* transformation.
* Let U be the N observations contained into data,
* and let S be the spacings contained into spacings,
* Assumes that S[0..N] contains the spacings
* between N real numbers
* U0,..., UN-1 in the interval [0, 1].
* These spacings are defined by
*
*
*
* Si = U(i) - U(i-1), 1 <= i < N,
*
* where U(0) = 0,
* U(N-1) = 1, and
*
* U(0),..., U(N-1), are the Ui sorted in increasing order.
* These spacings may have been obtained by calling
* {@link #diff(DoubleArrayList,DoubleArrayList,int,int,double,double) diff}.
* This method transforms the spacings into new
* spacings:
* it sorts
* S0,..., SN to obtain
*
* S(0) <= S(1) <= S(2) <= ... <= S(N),
* computes the weighted differences
*
*
*
*
* S0
* =
* (N + 1)S(0),
*
* S1
* =
* N(S(1) - S(0)),
*
* S2
* =
* (N - 1)(S(2) - S(1)),
*
*
* ...
*
*
* SN
* =
* S(N) - S(N-1),
*
*
*
*
* and computes
* Vi = S0 + S1 + ... + Si for
* 0 <= i < N.
* It then returns
* S0,..., SN in S[0..N] and
*
* V1,..., VN in V[1..N].
*
*
* Under the assumption that the Ui are i.i.d. U(0, 1), the new
* Si can be considered as a new set of spacings having the same
* distribution as the original spacings, and the Vi are a new sample
* of i.i.d. U(0, 1) random variables, sorted by increasing order.
*
*
* This transformation is useful to detect clustering in a data
* set: A pair of observations that are close to each other is transformed
* into an observation close to zero. A data set with unusually clustered
* observations is thus transformed to a data set with an
* accumulation of observations near zero, which is easily detected by
* the Anderson-Darling GOF test.
*
* @param data array of observations
*
* @param spacings spacings between the observations, will be filled with the new spacings
*
*
*/
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];
}
/**
* Applies the power ratios transformation W.
* Let U be the N observations contained into sortedData.
* Assumes that U contains N real numbers
*
* U(0),..., U(N-1) from the interval [0, 1],
* already sorted in increasing order, and computes the transformations:
*
*
*
* U'i = (U(i)/U(i+1))i+1, i = 0,..., N - 1,
*
* with
* U(N) = 1.
* These U'i are sorted in increasing order and put back in
* U[1...N].
* If the U(i) are i.i.d. U(0, 1) sorted by increasing order,
* then the U'i are also i.i.d. U(0, 1).
*
*
* This transformation is useful to detect clustering, as explained in
* {@link #iterateSpacings(DoubleArrayList,DoubleArrayList) iterateSpacings},
* except that here a pair of
* observations close to each other is transformed
* into an observation close to 1.
* An accumulation of observations near 1 is also easily detected by
* the Anderson-Darling GOF test.
*
* @param sortedData sorted array of real-valued observations in the interval [0, 1]
* that will be overwritten with the transformed observations
*
*/
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);
}
/**
* This class helps managing the partitions of possible outcomes
* into categories for applying chi-square tests.
* It permits one to automatically regroup categories to make sure that
* the expected number of observations in each category is large enough.
* To use this facility, one must first construct an
* OutcomeCategoriesChi2 object by passing to the constructor
* the expected number of observations for each original category.
* Then, calling the method {@link #regroupCategories regroupCategories} will regroup
* categories in a way that the expected number of observations in each
* category reaches a given threshold minExp.
* Experts in statistics recommend that minExp be always larger
* than or equal to 5 for the chi-square test to be valid. Thus,
* minExp = 10 is a safe value to use.
* After the call, nbExp gives the expected numbers in the new
* categories and loc[i] gives the relocation of category i,
* for each i. That is, loc[i] = j means that category i has
* been merged with category j because its original expected number was
* too small, and nbExp[i] has been added to nbExp[j]
* and then set to zero.
* In this case, all observations that previously belonged
* to category i are redirected to category j.
* The variable nbCategories gives the final number of categories,
* smin contains the new index of the lowest category,
* and smax the new index of the highest category.
*
*/
public static class OutcomeCategoriesChi2 {
/**
* Total number of categories.
*
*/
public int nbCategories;
/**
* Minimum index for valid expected numbers
* in the array nbExp.
*
*/
public int smin;
/**
* Maximum index for valid expected numbers
* in the array nbExp.
*
*/
public int smax;
/**
* Expected number of observations for each category.
*
*/
public double[] nbExp;
/**
* loc[i] gives the relocation of the category i in
* the nbExp array.
*
*/
public int[] loc;
/**
* Constructs an OutcomeCategoriesChi2 object
* using the array nbExp for the number of expected observations in
* each category. The smin and smax fields are set to 0 and
* (n - 1) respectively, where n is the length of array nbExp.
* The loc field is set such that loc[i]=i for each i.
* The field nbCategories is set to n.
*
* @param nbExp array of expected observations for each category
*
*
*/
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;
}
/**
* Constructs an OutcomeCategoriesChi2 object using the
* given nbExp expected observations array. Only the expected
* numbers from the smin to smax (inclusive) indices will be
* considered valid. The loc field is set such that loc[i]=i
* for each i in the interval [smin, smax]. All loc[i]
* for i <= smin are set to smin, and all loc[i] for
* i >= smax are set to smax.
* The field nbCategories is set to (smax - smin + 1).
*
* @param nbExp array of expected observations for each category
*
* @param smin Minimum index for valid expected number of observations
*
* @param smax Maximum index for valid expected number of observations
*
*
*/
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;
}
/**
* Constructs an OutcomeCategoriesChi2 object.
* The field nbCategories is set to nbCat.
*
* @param nbExp array of expected observations for each category
*
* @param smin Minimum index for valid expected number of observations
*
* @param smax Maximum index for valid expected number of observations
*
* @param loc array for which loc[i] gives the relocation of the category i
*
*
*/
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;
}
/**
* Regroup categories as explained earlier, so that the expected
* number of observations in each category is at least minExp.
* We usually choose minExp = 10.
*
* @param minExp mininum number of expected observations in each category
*
*
*/
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");
}
/**
* Provides a report on the categories.
*
* @return the categories represented as a string
*
*/
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();
}
}
/**
* Computes and returns the chi-square statistic for the
* observations oi in count[smin...smax], for which the
* corresponding expected values ei are in nbExp[smin...smax].
* Assuming that i goes from 1 to k, where k = smax-smin+1
* is the number of categories, the chi-square statistic is defined as
*
*
*
* X2 = ∑i=1k(oi - ei)2/ei.
*
* Under the hypothesis that the ei are the correct expectations and
* if these ei are large enough, X2 follows approximately the
* chi-square distribution with k - 1 degrees of freedom.
* If some of the ei are too small, one can use
* OutcomeCategoriesChi2 to regroup categories.
*
* @param nbExp number of expected observations in each category (or interval)
*
* @param count number of counted observations in each category
*
* @param smin index of the first valid data in count and nbExp
*
* @param smax index of the last valid data in count and nbExp
*
* @return the X2 statistic
*
*/
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;
}
/**
* Computes and returns the chi-square statistic for the
* observations stored in data, assuming that these observations follow
* the discrete distribution dist. For dist, we assume that
* there is one set
* S = {a, a + 1,..., b - 1, b}, where a < b and a >= 0,
* for which p(s) > 0 if s∈S and p(s) = 0 otherwise.
*
*
* Generally, it is not possible to divide the integers in intervals satisfying
*
* nP(a0 <= s < a1) = nP(a1 <= s < a2) = ... = nP(aj-1 <= s < aj)
* for a discrete distribution, where n is the sample size, i.e.,
* the number of
* observations stored into data.
* To perform a general chi-square test, the method starts
* from smin and finds the first non-negligible
* probability
* p(s) >= ε, where
* ε = {@link DiscreteDistributionInt#EPSILON DiscreteDistributionInt.EPSILON}.
* It uses smax to allocate an array storing the
* number of expected observations (np(s)) for each s >= smin.
* Starting from s = smin, the np(s) terms are computed and
* the allocated array grows if required until a negligible probability
* term is found.
* This gives the number of expected elements for
* each category, where an outcome category corresponds here to
* an interval in which sample observations could lie.
* The categories are regrouped to have at least
* minExp observations per category.
* The method then counts the number of samples in each categories and calls
* {@link #chi2(double[],int[],int,int) chi2} to get the chi-square test
* statistic. If numCat is not
* null, the number of categories after regrouping is returned
* in numCat[0]. The number of degrees of freedom is equal to
* numCat[0]-1. We usually choose minExp = 10.
*
* @param data observations, not necessarily sorted
*
* @param dist assumed probability distribution
*
* @param smin estimated minimum value of s for which p(s) > 0
*
* @param smax estimated maximum value of s for which p(s) > 0
*
* @param minExp minimum number of expected observations in each
* interval
*
* @param numCat one-element array that will be filled with the number of
* categories after regrouping
*
* @return the chi-square statistic for a discrete distribution
*
*/
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);
}
/**
* Similar to {@link #chi2(double[],int[],int,int) chi2},
* except that the expected
* number of observations per category is assumed to be the same for
* all categories, and equal to nbExp.
*
* @param nbExp number of expected observations in each category (or interval)
*
* @param count number of counted observations in each category
*
* @param smin index of the first valid data in count and nbExp
*
* @param smax index of the last valid data in count and nbExp
*
* @return the X2 statistic
*
*/
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;
}
/**
* Computes the chi-square statistic for a continuous distribution.
* Here, the equiprobable case can be used. Assuming that data contains
* observations coming from the uniform distribution, the [0, 1] interval
* is divided into 1/p subintervals, where p = minExp/n, n
* being the sample size, i.e., the number of observations stored in
* data. For each subinterval, the method counts the number of
* contained observations and the chi-square statistic is computed
* using {@link #chi2Equal(double,int[],int,int) chi2Equal}.
* We usually choose minExp = 10.
*
* @param data array of observations in [0, 1)
*
* @param minExp minimum number of expected observations in each subintervals
*
* @return the chi-square statistic for a continuous distribution
*
*/
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);
}
/**
* Equivalent to chi2Equal (data, 10).
*
* @param data array of observations in [0, 1)
*
* @return the chi-square statistic for a continuous distribution
*
*/
public static double chi2Equal (DoubleArrayList data) {
return chi2Equal (data, 10.0);
}
/**
* Computes and returns the scan statistic SN(d ),
* defined in {@link FBar#scan FBar.scan}.
* Let U be the N observations contained into sortedData.
* The N observations in U[0..N - 1] must be real numbers
* in the interval [0, 1], sorted in increasing order.
* (See {@link FBar#scan FBar.scan} for the distribution function of SN(d )).
*
* @param sortedData sorted array of real-valued observations in the interval [0, 1]
*
* @param d length of the test interval (∈(0, 1))
*
* @return the scan statistic
*
*/
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;
}
/**
* Computes and returns the Cramér-von Mises statistic WN2. It is
* defined by
*
*
*
* WN2 = 1/(12N) + ∑j=0N-1(U(j) - (j + 0.5)/N)2,
*
* assuming that sortedData contains
* U(0),..., U(N-1)
* sorted in increasing order.
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @return the Cramér-von Mises statistic
*
*/
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;
}
/**
* Computes and returns the Watson statistic GN. It is
* defined by
*
*
*
*
* GN
* =
* (N)1/2max0 <= j <= N-1{(j + 1)/N - U(j) + bar(U)N -1/2}
*
*
* =
* (N)1/2(DN+ + bar(U)N - 1/2),
*
*
*
*
* where
* bar(U)N is the average of the observations U(j),
* assuming that sortedData contains the sorted
* U(0),..., U(N-1).
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @return the Watson statistic GN
*
*/
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;
}
/**
* Computes and returns the Watson statistic UN2. It is
* defined by
*
*
*
*
* WN2
* =
* 1/(12N) + ∑j=0N-1{U(j) - (j + 0.5)/N}2,
*
* UN2
* =
* WN2 - N(bar(U)N -1/2)2.
*
*
*
*
* where
* bar(U)N is the average of the observations U(j),
* assuming that sortedData contains the sorted
*
* U(0),..., U(N-1).
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @return the Watson statistic UN2
*
*
*/
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.0;
/**
* Computes and returns the Anderson-Darling statistic AN2.
* It is
* defined by
*
*
*
*
* AN2
* =
* - N - 1/N ∑j=0N-1{(2j + 1)ln(U(j)) + (2N - 1 - 2j)ln(1 - U(j))},
*
*
*
*
* assuming that sortedData contains
* U(0),..., U(N-1).
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @return the Anderson-Darling statistic
*
*/
public static double andersonDarling (DoubleArrayList sortedData) {
double u1;
double u, a2;
double[] v = sortedData.elements();
int n = sortedData.size();
if (n <= 0) {
System.err.println ("andersonDarling: n <= 0");
return 0.0;
}
a2 = 0.0;
for (int i = 0; i < n; i++) {
u = v[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;
}
/**
* Computes the Kolmogorov-Smirnov (KS) test statistics
* DN+, DN-, and DN. It is
* defined by
*
*
*
*
* DN+
* =
* max0 <= j <= N-1((j + 1)/N - U(j)),
*
* DN-
* =
* max0 <= j <= N-1(U(j) - j/N),
*
* DN
* =
* max (DN+, DN-).
*
*
*
*
* and returns an array of length 3 that contains their values at positions
* 0, 1, and 2, respectively.
*
*
* These statistics compare the empirical distribution of
*
* U(1),..., U(N), which are assumed to be in sortedData,
* with the uniform distribution.
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @return an array containing DN+, DN- and DN, in that order
*
*/
public static double[] kolmogorovSmirnov (DoubleArrayList sortedData) {
double[] ret = new double[3];
double[] u = sortedData.elements();
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;
}
/**
* Compute the KS statistics DN+(a) and DN-(a) defined in
* the description of the method
* {@link FDist#kolmogorovSmirnovPlusJumpOne FDist.kolmogorovSmirnovPlusJumpOne}, assuming that F is the
* uniform distribution over [0, 1] and that
*
* U(1),..., U(N) are in sortedData.
* Returns an array of length 2 that contains their values at positions
* 0 and 1, respectively.
*
* @param sortedData array of sorted real-valued observations in the interval [0, 1]
*
* @param a size of the jump
*
* @return an array containing DN+ and DN-, in that order
*
*/
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;
}
/**
* Computes a variant of the p-value p whenever a test statistic
* has a discrete probability distribution.
* This p-value is defined as follows:
*
*
*
*
* pL
* =
* P[Y <= y]
*
* pR
* =
* P[Y >= y]
*
*
*
*
*
*
*
*
* p =
* pR,
* if pR < pL,
* p =
* 1 - pL,
* if pR >= pL and pL < 0.5,
* p =
* 0.5
* otherwise.
*
*
* The function takes pL and pR as input and returns p.
*
* @param pL left p-value
*
* @param pR right p-value
*
* @return the p-value for a test on a discrete distribution
*/
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;
}
}
Ads help maintain this website.