001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math4.stat.inference;
019
020import java.math.RoundingMode;
021import java.util.Arrays;
022
023import org.apache.commons.rng.simple.RandomSource;
024import org.apache.commons.rng.UniformRandomProvider;
025import org.apache.commons.statistics.distribution.ContinuousDistribution;
026import org.apache.commons.numbers.combinatorics.BinomialCoefficientDouble;
027import org.apache.commons.numbers.fraction.BigFraction;
028import org.apache.commons.numbers.field.BigFractionField;
029import org.apache.commons.math4.distribution.EnumeratedRealDistribution;
030import org.apache.commons.math4.distribution.AbstractRealDistribution;
031import org.apache.commons.math4.exception.InsufficientDataException;
032import org.apache.commons.math4.exception.MathArithmeticException;
033import org.apache.commons.math4.exception.MaxCountExceededException;
034import org.apache.commons.math4.exception.NullArgumentException;
035import org.apache.commons.math4.exception.NumberIsTooLargeException;
036import org.apache.commons.math4.exception.OutOfRangeException;
037import org.apache.commons.math4.exception.TooManyIterationsException;
038import org.apache.commons.math4.exception.NotANumberException;
039import org.apache.commons.math4.exception.util.LocalizedFormats;
040import org.apache.commons.math4.linear.MatrixUtils;
041import org.apache.commons.math4.linear.RealMatrix;
042import org.apache.commons.math4.util.FastMath;
043import org.apache.commons.math4.util.MathArrays;
044import org.apache.commons.math4.util.MathUtils;
045import org.apache.commons.math4.field.linalg.FieldDenseMatrix;
046
047/**
048 * Implementation of the <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test">
049 * Kolmogorov-Smirnov (K-S) test</a> for equality of continuous distributions.
050 * <p>
051 * The K-S test uses a statistic based on the maximum deviation of the empirical distribution of
052 * sample data points from the distribution expected under the null hypothesis. For one-sample tests
053 * evaluating the null hypothesis that a set of sample data points follow a given distribution, the
054 * test statistic is \(D_n=\sup_x |F_n(x)-F(x)|\), where \(F\) is the expected distribution and
055 * \(F_n\) is the empirical distribution of the \(n\) sample data points. The distribution of
056 * \(D_n\) is estimated using a method based on [1] with certain quick decisions for extreme values
057 * given in [2].
058 * </p>
059 * <p>
060 * Two-sample tests are also supported, evaluating the null hypothesis that the two samples
061 * {@code x} and {@code y} come from the same underlying distribution. In this case, the test
062 * statistic is \(D_{n,m}=\sup_t | F_n(t)-F_m(t)|\) where \(n\) is the length of {@code x}, \(m\) is
063 * the length of {@code y}, \(F_n\) is the empirical distribution that puts mass \(1/n\) at each of
064 * the values in {@code x} and \(F_m\) is the empirical distribution of the {@code y} values. The
065 * default 2-sample test method, {@link #kolmogorovSmirnovTest(double[], double[])} works as
066 * follows:
067 * <ul>
068 * <li>When the product of the sample sizes is less than 10000, the method presented in [4]
069 * is used to compute the exact p-value for the 2-sample test.</li>
070 * <li>When the product of the sample sizes is larger, the asymptotic
071 * distribution of \(D_{n,m}\) is used. See {@link #approximateP(double, int, int)} for details on
072 * the approximation.</li>
073 * </ul><p>
074 * For small samples (former case), if the data contains ties, random jitter is added
075 * to the sample data to break ties before applying the algorithm above. Alternatively,
076 * the {@link #bootstrap(double[],double[],int,boolean,UniformRandomProvider)}
077 * method, modeled after <a href="http://sekhon.berkeley.edu/matching/ks.boot.html">ks.boot</a>
078 * in the R Matching package [3], can be used if ties are known to be present in the data.
079 * </p>
080 * <p>
081 * In the two-sample case, \(D_{n,m}\) has a discrete distribution. This makes the p-value
082 * associated with the null hypothesis \(H_0 : D_{n,m} \ge d \) differ from \(H_0 : D_{n,m} \ge d \)
083 * by the mass of the observed value \(d\). To distinguish these, the two-sample tests use a boolean
084 * {@code strict} parameter. This parameter is ignored for large samples.
085 * </p>
086 * <p>
087 * The methods used by the 2-sample default implementation are also exposed directly:
088 * <ul>
089 * <li>{@link #exactP(double, int, int, boolean)} computes exact 2-sample p-values</li>
090 * <li>{@link #approximateP(double, int, int)} uses the asymptotic distribution The {@code boolean}
091 * arguments in the first two methods allow the probability used to estimate the p-value to be
092 * expressed using strict or non-strict inequality. See
093 * {@link #kolmogorovSmirnovTest(double[], double[], boolean)}.</li>
094 * </ul>
095 * <p>
096 * References:
097 * <ul>
098 * <li>[1] <a href="http://www.jstatsoft.org/v08/i18/"> Evaluating Kolmogorov's Distribution</a> by
099 * George Marsaglia, Wai Wan Tsang, and Jingbo Wang</li>
100 * <li>[2] <a href="http://www.jstatsoft.org/v39/i11/"> Computing the Two-Sided Kolmogorov-Smirnov
101 * Distribution</a> by Richard Simard and Pierre L'Ecuyer</li>
102 * <li>[3] Jasjeet S. Sekhon. 2011. <a href="http://www.jstatsoft.org/article/view/v042i07">
103 * Multivariate and Propensity Score Matching Software with Automated Balance Optimization:
104 * The Matching package for R</a> Journal of Statistical Software, 42(7): 1-52.</li>
105 * <li>[4] Wilcox, Rand. 2012. Introduction to Robust Estimation and Hypothesis Testing,
106 * Chapter 5, 3rd Ed. Academic Press.</li>
107 * </ul>
108 * <br>
109 * Note that [1] contains an error in computing h, refer to <a
110 * href="https://issues.apache.org/jira/browse/MATH-437">MATH-437</a> for details.
111 *
112 * @since 3.3
113 */
114public class KolmogorovSmirnovTest {
115    /**
116     * Bound on the number of partial sums in {@link #ksSum(double, double, int)}
117     */
118    private static final int MAXIMUM_PARTIAL_SUM_COUNT = 100000;
119    /** Convergence criterion for {@link #ksSum(double, double, int)} */
120    private static final double KS_SUM_CAUCHY_CRITERION = 1e-20;
121    /** Convergence criterion for the sums in {@link #pelzGood(double, int)} */
122    private static final double PG_SUM_RELATIVE_ERROR = 1e-10;
123    /** 1/2 */
124    private static final BigFraction ONE_HALF = BigFraction.of(1, 2);
125
126    /**
127     * When product of sample sizes exceeds this value, 2-sample K-S test uses asymptotic
128     * distribution to compute the p-value.
129     */
130    private static final int LARGE_SAMPLE_PRODUCT = 10000;
131
132    /**
133     * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a one-sample <a
134     * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
135     * evaluating the null hypothesis that {@code data} conforms to {@code distribution}. If
136     * {@code exact} is true, the distribution used to compute the p-value is computed using
137     * extended precision. See {@link #cdfExact(double, int)}.
138     *
139     * @param distribution reference distribution
140     * @param data sample being being evaluated
141     * @param exact whether or not to force exact computation of the p-value
142     * @return the p-value associated with the null hypothesis that {@code data} is a sample from
143     *         {@code distribution}
144     * @throws InsufficientDataException if {@code data} does not have length at least 2
145     * @throws NullArgumentException if {@code data} is null
146     */
147    public double kolmogorovSmirnovTest(ContinuousDistribution distribution, double[] data, boolean exact) {
148        return 1d - cdf(kolmogorovSmirnovStatistic(distribution, data), data.length, exact);
149    }
150
151    /**
152     * Computes the one-sample Kolmogorov-Smirnov test statistic, \(D_n=\sup_x |F_n(x)-F(x)|\) where
153     * \(F\) is the distribution (cdf) function associated with {@code distribution}, \(n\) is the
154     * length of {@code data} and \(F_n\) is the empirical distribution that puts mass \(1/n\) at
155     * each of the values in {@code data}.
156     *
157     * @param distribution reference distribution
158     * @param data sample being evaluated
159     * @return Kolmogorov-Smirnov statistic \(D_n\)
160     * @throws InsufficientDataException if {@code data} does not have length at least 2
161     * @throws NullArgumentException if {@code data} is null
162     */
163    public double kolmogorovSmirnovStatistic(ContinuousDistribution distribution, double[] data) {
164        checkArray(data);
165        final int n = data.length;
166        final double nd = n;
167        final double[] dataCopy = new double[n];
168        System.arraycopy(data, 0, dataCopy, 0, n);
169        Arrays.sort(dataCopy);
170        double d = 0d;
171        for (int i = 1; i <= n; i++) {
172            final double yi = distribution.cumulativeProbability(dataCopy[i - 1]);
173            final double currD = FastMath.max(yi - (i - 1) / nd, i / nd - yi);
174            if (currD > d) {
175                d = currD;
176            }
177        }
178        return d;
179    }
180
181    /**
182     * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a two-sample <a
183     * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
184     * evaluating the null hypothesis that {@code x} and {@code y} are samples drawn from the same
185     * probability distribution. Specifically, what is returned is an estimate of the probability
186     * that the {@link #kolmogorovSmirnovStatistic(double[], double[])} associated with a randomly
187     * selected partition of the combined sample into subsamples of sizes {@code x.length} and
188     * {@code y.length} will strictly exceed (if {@code strict} is {@code true}) or be at least as
189     * large as (if {@code strict} is {@code false}) as {@code kolmogorovSmirnovStatistic(x, y)}.
190     *
191     * @param x first sample dataset.
192     * @param y second sample dataset.
193     * @param strict whether or not the probability to compute is expressed as
194     * a strict inequality (ignored for large samples).
195     * @return p-value associated with the null hypothesis that {@code x} and
196     * {@code y} represent samples from the same distribution.
197     * @throws InsufficientDataException if either {@code x} or {@code y} does
198     * not have length at least 2.
199     * @throws NullArgumentException if either {@code x} or {@code y} is null.
200     * @throws NotANumberException if the input arrays contain NaN values.
201     *
202     * @see #bootstrap(double[],double[],int,boolean,UniformRandomProvider)
203     */
204    public double kolmogorovSmirnovTest(double[] x, double[] y, boolean strict) {
205        final long lengthProduct = (long) x.length * y.length;
206        double[] xa = null;
207        double[] ya = null;
208        if (lengthProduct < LARGE_SAMPLE_PRODUCT && hasTies(x,y)) {
209            xa = Arrays.copyOf(x, x.length);
210            ya = Arrays.copyOf(y, y.length);
211            fixTies(xa, ya);
212        } else {
213            xa = x;
214            ya = y;
215        }
216        if (lengthProduct < LARGE_SAMPLE_PRODUCT) {
217            return exactP(kolmogorovSmirnovStatistic(xa, ya), x.length, y.length, strict);
218        }
219        return approximateP(kolmogorovSmirnovStatistic(x, y), x.length, y.length);
220    }
221
222    /**
223     * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a two-sample <a
224     * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
225     * evaluating the null hypothesis that {@code x} and {@code y} are samples drawn from the same
226     * probability distribution. Assumes the strict form of the inequality used to compute the
227     * p-value. See {@link #kolmogorovSmirnovTest(ContinuousDistribution, double[], boolean)}.
228     *
229     * @param x first sample dataset
230     * @param y second sample dataset
231     * @return p-value associated with the null hypothesis that {@code x} and {@code y} represent
232     *         samples from the same distribution
233     * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
234     *         least 2
235     * @throws NullArgumentException if either {@code x} or {@code y} is null
236     */
237    public double kolmogorovSmirnovTest(double[] x, double[] y) {
238        return kolmogorovSmirnovTest(x, y, true);
239    }
240
241    /**
242     * Computes the two-sample Kolmogorov-Smirnov test statistic, \(D_{n,m}=\sup_x |F_n(x)-F_m(x)|\)
243     * where \(n\) is the length of {@code x}, \(m\) is the length of {@code y}, \(F_n\) is the
244     * empirical distribution that puts mass \(1/n\) at each of the values in {@code x} and \(F_m\)
245     * is the empirical distribution of the {@code y} values.
246     *
247     * @param x first sample
248     * @param y second sample
249     * @return test statistic \(D_{n,m}\) used to evaluate the null hypothesis that {@code x} and
250     *         {@code y} represent samples from the same underlying distribution
251     * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
252     *         least 2
253     * @throws NullArgumentException if either {@code x} or {@code y} is null
254     */
255    public double kolmogorovSmirnovStatistic(double[] x, double[] y) {
256        return integralKolmogorovSmirnovStatistic(x, y)/((double)(x.length * (long)y.length));
257    }
258
259    /**
260     * Computes the two-sample Kolmogorov-Smirnov test statistic, \(D_{n,m}=\sup_x |F_n(x)-F_m(x)|\)
261     * where \(n\) is the length of {@code x}, \(m\) is the length of {@code y}, \(F_n\) is the
262     * empirical distribution that puts mass \(1/n\) at each of the values in {@code x} and \(F_m\)
263     * is the empirical distribution of the {@code y} values. Finally \(n m D_{n,m}\) is returned
264     * as long value.
265     *
266     * @param x first sample
267     * @param y second sample
268     * @return test statistic \(n m D_{n,m}\) used to evaluate the null hypothesis that {@code x} and
269     *         {@code y} represent samples from the same underlying distribution
270     * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
271     *         least 2
272     * @throws NullArgumentException if either {@code x} or {@code y} is null
273     */
274    private long integralKolmogorovSmirnovStatistic(double[] x, double[] y) {
275        checkArray(x);
276        checkArray(y);
277        // Copy and sort the sample arrays
278        final double[] sx = Arrays.copyOf(x, x.length);
279        final double[] sy = Arrays.copyOf(y, y.length);
280        Arrays.sort(sx);
281        Arrays.sort(sy);
282        final int n = sx.length;
283        final int m = sy.length;
284
285        int rankX = 0;
286        int rankY = 0;
287        long curD = 0l;
288
289        // Find the max difference between cdf_x and cdf_y
290        long supD = 0l;
291        do {
292            double z = Double.compare(sx[rankX], sy[rankY]) <= 0 ? sx[rankX] : sy[rankY];
293            while(rankX < n && Double.compare(sx[rankX], z) == 0) {
294                rankX += 1;
295                curD += m;
296            }
297            while(rankY < m && Double.compare(sy[rankY], z) == 0) {
298                rankY += 1;
299                curD -= n;
300            }
301            if (curD > supD) {
302                supD = curD;
303            }
304            else if (-curD > supD) {
305                supD = -curD;
306            }
307        } while(rankX < n && rankY < m);
308        return supD;
309    }
310
311    /**
312     * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a one-sample <a
313     * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
314     * evaluating the null hypothesis that {@code data} conforms to {@code distribution}.
315     *
316     * @param distribution reference distribution
317     * @param data sample being being evaluated
318     * @return the p-value associated with the null hypothesis that {@code data} is a sample from
319     *         {@code distribution}
320     * @throws InsufficientDataException if {@code data} does not have length at least 2
321     * @throws NullArgumentException if {@code data} is null
322     */
323    public double kolmogorovSmirnovTest(ContinuousDistribution distribution, double[] data) {
324        return kolmogorovSmirnovTest(distribution, data, false);
325    }
326
327    /**
328     * Performs a <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov
329     * test</a> evaluating the null hypothesis that {@code data} conforms to {@code distribution}.
330     *
331     * @param distribution reference distribution
332     * @param data sample being being evaluated
333     * @param alpha significance level of the test
334     * @return true iff the null hypothesis that {@code data} is a sample from {@code distribution}
335     *         can be rejected with confidence 1 - {@code alpha}
336     * @throws InsufficientDataException if {@code data} does not have length at least 2
337     * @throws NullArgumentException if {@code data} is null
338     */
339    public boolean kolmogorovSmirnovTest(ContinuousDistribution distribution, double[] data, double alpha) {
340        if ((alpha <= 0) || (alpha > 0.5)) {
341            throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL, alpha, 0, 0.5);
342        }
343        return kolmogorovSmirnovTest(distribution, data) < alpha;
344    }
345
346    /**
347     * Estimates the <i>p-value</i> of a two-sample
348     * <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test">Kolmogorov-Smirnov test</a>
349     * evaluating the null hypothesis that {@code x} and {@code y} are samples
350     * drawn from the same probability distribution.
351     * This method estimates the p-value by repeatedly sampling sets of size
352     * {@code x.length} and {@code y.length} from the empirical distribution
353     * of the combined sample.
354     * When {@code strict} is true, this is equivalent to the algorithm implemented
355     * in the R function {@code ks.boot}, described in <pre>
356     * Jasjeet S. Sekhon. 2011. 'Multivariate and Propensity Score Matching
357     * Software with Automated Balance Optimization: The Matching package for R.'
358     * Journal of Statistical Software, 42(7): 1-52.
359     * </pre>
360     *
361     * @param x First sample.
362     * @param y Second sample.
363     * @param iterations Number of bootstrap resampling iterations.
364     * @param strict Whether or not the null hypothesis is expressed as a strict inequality.
365     * @param rng RNG for creating the sampling sets.
366     * @return the estimated p-value.
367     */
368    public double bootstrap(double[] x,
369                            double[] y,
370                            int iterations,
371                            boolean strict,
372                            UniformRandomProvider rng) {
373        final int xLength = x.length;
374        final int yLength = y.length;
375        final double[] combined = new double[xLength + yLength];
376        System.arraycopy(x, 0, combined, 0, xLength);
377        System.arraycopy(y, 0, combined, xLength, yLength);
378        final ContinuousDistribution.Sampler sampler = new EnumeratedRealDistribution(combined).createSampler(rng);
379        final long d = integralKolmogorovSmirnovStatistic(x, y);
380        int greaterCount = 0;
381        int equalCount = 0;
382        double[] curX;
383        double[] curY;
384        long curD;
385        for (int i = 0; i < iterations; i++) {
386            curX = AbstractRealDistribution.sample(xLength, sampler);
387            curY = AbstractRealDistribution.sample(yLength, sampler);
388            curD = integralKolmogorovSmirnovStatistic(curX, curY);
389            if (curD > d) {
390                greaterCount++;
391            } else if (curD == d) {
392                equalCount++;
393            }
394        }
395        return strict ? greaterCount / (double) iterations :
396            (greaterCount + equalCount) / (double) iterations;
397    }
398
399    /**
400     * Calculates \(P(D_n &lt; d)\) using the method described in [1] with quick decisions for extreme
401     * values given in [2] (see above). The result is not exact as with
402     * {@link #cdfExact(double, int)} because calculations are based on
403     * {@code double} rather than {@link BigFraction}.
404     *
405     * @param d statistic
406     * @param n sample size
407     * @return \(P(D_n &lt; d)\)
408     * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
409     * {@link BigFraction} in expressing {@code d} as
410     * \((k - h) / m\) for integer {@code k, m} and \(0 \le h &lt; 1\)
411     */
412    public double cdf(double d, int n) {
413        return cdf(d, n, false);
414    }
415
416    /**
417     * Calculates {@code P(D_n < d)}. The result is exact in the sense that BigFraction/BigReal is
418     * used everywhere at the expense of very slow execution time. Almost never choose this in real
419     * applications unless you are very sure; this is almost solely for verification purposes.
420     * Normally, you would choose {@link #cdf(double, int)}. See the class
421     * javadoc for definitions and algorithm description.
422     *
423     * @param d statistic
424     * @param n sample size
425     * @return \(P(D_n &lt; d)\)
426     * @throws MathArithmeticException if the algorithm fails to convert {@code h} to a
427     * {@link BigFraction} in expressing {@code d} as
428     * \((k - h) / m\) for integer {@code k, m} and \(0 \le h &lt; 1\)
429     */
430    public double cdfExact(double d, int n) {
431        return cdf(d, n, true);
432    }
433
434    /**
435     * Calculates {@code P(D_n < d)} using method described in [1] with quick decisions for extreme
436     * values given in [2] (see above).
437     *
438     * @param d statistic
439     * @param n sample size
440     * @param exact whether the probability should be calculated exact using
441     *        {@link BigFraction} everywhere at the expense of
442     *        very slow execution time, or if {@code double} should be used convenient places to
443     *        gain speed. Almost never choose {@code true} in real applications unless you are very
444     *        sure; {@code true} is almost solely for verification purposes.
445     * @return \(P(D_n &lt; d)\)
446     * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
447     * {@link BigFraction} in expressing {@code d} as
448     * \((k - h) / m\) for integer {@code k, m} and \(0 \le h &lt; 1\).
449     */
450    public double cdf(double d, int n, boolean exact) {
451        final double ninv = 1 / ((double) n);
452        final double ninvhalf = 0.5 * ninv;
453
454        if (d <= ninvhalf) {
455            return 0;
456        } else if (ninvhalf < d && d <= ninv) {
457            double res = 1;
458            final double f = 2 * d - ninv;
459            // n! f^n = n*f * (n-1)*f * ... * 1*x
460            for (int i = 1; i <= n; ++i) {
461                res *= i * f;
462            }
463            return res;
464        } else if (1 - ninv <= d && d < 1) {
465            return 1 - 2 * Math.pow(1 - d, n);
466        } else if (1 <= d) {
467            return 1;
468        }
469        if (exact) {
470            return exactK(d, n);
471        }
472        if (n <= 140) {
473            return roundedK(d, n);
474        }
475        return pelzGood(d, n);
476    }
477
478    /**
479     * Calculates the exact value of {@code P(D_n < d)} using the method described in [1] (reference
480     * in class javadoc above) and {@link BigFraction} (see above).
481     *
482     * @param d statistic
483     * @param n sample size
484     * @return the two-sided probability of \(P(D_n &lt; d)\)
485     * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
486     * {@link BigFraction}.
487     */
488    private double exactK(double d, int n) {
489        final int k = (int) Math.ceil(n * d);
490
491        final FieldDenseMatrix<BigFraction> H;
492        try {
493            H = createExactH(d, n);
494        } catch (ArithmeticException e) {
495            throw new MathArithmeticException(LocalizedFormats.FRACTION);
496        }
497
498        final FieldDenseMatrix<BigFraction> Hpower = H.pow(n);
499
500        BigFraction pFrac = Hpower.get(k - 1, k - 1);
501
502        for (int i = 1; i <= n; ++i) {
503            pFrac = pFrac.multiply(i).divide(n);
504        }
505
506        /*
507         * BigFraction.doubleValue converts numerator to double and the denominator to double and
508         * divides afterwards. That gives NaN quite easy. This does not (scale is the number of
509         * digits):
510         */
511        return pFrac.bigDecimalValue(20, RoundingMode.HALF_UP).doubleValue();
512    }
513
514    /**
515     * Calculates {@code P(D_n < d)} using method described in [1] and doubles (see above).
516     *
517     * @param d statistic
518     * @param n sample size
519     * @return \(P(D_n &lt; d)\)
520     */
521    private double roundedK(double d, int n) {
522
523        final int k = (int) Math.ceil(n * d);
524        final RealMatrix H = this.createRoundedH(d, n);
525        final RealMatrix Hpower = H.power(n);
526
527        double pFrac = Hpower.getEntry(k - 1, k - 1);
528        for (int i = 1; i <= n; ++i) {
529            pFrac *= (double) i / (double) n;
530        }
531
532        return pFrac;
533    }
534
535    /**
536     * Computes the Pelz-Good approximation for \(P(D_n &lt; d)\) as described in [2] in the class javadoc.
537     *
538     * @param d value of d-statistic (x in [2])
539     * @param n sample size
540     * @return \(P(D_n &lt; d)\)
541     * @since 3.4
542     */
543    public double pelzGood(double d, int n) {
544        // Change the variable since approximation is for the distribution evaluated at d / sqrt(n)
545        final double sqrtN = FastMath.sqrt(n);
546        final double z = d * sqrtN;
547        final double z2 = d * d * n;
548        final double z4 = z2 * z2;
549        final double z6 = z4 * z2;
550        final double z8 = z4 * z4;
551
552        // Eventual return value
553        double ret = 0;
554
555        // Compute K_0(z)
556        double sum = 0;
557        double increment = 0;
558        double kTerm = 0;
559        double z2Term = MathUtils.PI_SQUARED / (8 * z2);
560        int k = 1;
561        for (; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
562            kTerm = 2 * k - 1;
563            increment = FastMath.exp(-z2Term * kTerm * kTerm);
564            sum += increment;
565            if (increment <= PG_SUM_RELATIVE_ERROR * sum) {
566                break;
567            }
568        }
569        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
570            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
571        }
572        ret = sum * FastMath.sqrt(2 * FastMath.PI) / z;
573
574        // K_1(z)
575        // Sum is -inf to inf, but k term is always (k + 1/2) ^ 2, so really have
576        // twice the sum from k = 0 to inf (k = -1 is same as 0, -2 same as 1, ...)
577        final double twoZ2 = 2 * z2;
578        sum = 0;
579        kTerm = 0;
580        double kTerm2 = 0;
581        for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
582            kTerm = k + 0.5;
583            kTerm2 = kTerm * kTerm;
584            increment = (MathUtils.PI_SQUARED * kTerm2 - z2) * FastMath.exp(-MathUtils.PI_SQUARED * kTerm2 / twoZ2);
585            sum += increment;
586            if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
587                break;
588            }
589        }
590        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
591            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
592        }
593        final double sqrtHalfPi = FastMath.sqrt(FastMath.PI / 2);
594        // Instead of doubling sum, divide by 3 instead of 6
595        ret += sum * sqrtHalfPi / (3 * z4 * sqrtN);
596
597        // K_2(z)
598        // Same drill as K_1, but with two doubly infinite sums, all k terms are even powers.
599        final double z4Term = 2 * z4;
600        final double z6Term = 6 * z6;
601        z2Term = 5 * z2;
602        final double pi4 = MathUtils.PI_SQUARED * MathUtils.PI_SQUARED;
603        sum = 0;
604        kTerm = 0;
605        kTerm2 = 0;
606        for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
607            kTerm = k + 0.5;
608            kTerm2 = kTerm * kTerm;
609            increment =  (z6Term + z4Term + MathUtils.PI_SQUARED * (z4Term - z2Term) * kTerm2 +
610                    pi4 * (1 - twoZ2) * kTerm2 * kTerm2) * FastMath.exp(-MathUtils.PI_SQUARED * kTerm2 / twoZ2);
611            sum += increment;
612            if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
613                break;
614            }
615        }
616        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
617            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
618        }
619        double sum2 = 0;
620        kTerm2 = 0;
621        for (k = 1; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
622            kTerm2 = k * k;
623            increment = MathUtils.PI_SQUARED * kTerm2 * FastMath.exp(-MathUtils.PI_SQUARED * kTerm2 / twoZ2);
624            sum2 += increment;
625            if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum2)) {
626                break;
627            }
628        }
629        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
630            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
631        }
632        // Again, adjust coefficients instead of doubling sum, sum2
633        ret += (sqrtHalfPi / n) * (sum / (36 * z2 * z2 * z2 * z) - sum2 / (18 * z2 * z));
634
635        // K_3(z) One more time with feeling - two doubly infinite sums, all k powers even.
636        // Multiply coefficient denominators by 2, so omit doubling sums.
637        final double pi6 = pi4 * MathUtils.PI_SQUARED;
638        sum = 0;
639        double kTerm4 = 0;
640        double kTerm6 = 0;
641        for (k = 0; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
642            kTerm = k + 0.5;
643            kTerm2 = kTerm * kTerm;
644            kTerm4 = kTerm2 * kTerm2;
645            kTerm6 = kTerm4 * kTerm2;
646            increment = (pi6 * kTerm6 * (5 - 30 * z2) + pi4 * kTerm4 * (-60 * z2 + 212 * z4) +
647                            MathUtils.PI_SQUARED * kTerm2 * (135 * z4 - 96 * z6) - 30 * z6 - 90 * z8) *
648                    FastMath.exp(-MathUtils.PI_SQUARED * kTerm2 / twoZ2);
649            sum += increment;
650            if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum)) {
651                break;
652            }
653        }
654        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
655            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
656        }
657        sum2 = 0;
658        for (k = 1; k < MAXIMUM_PARTIAL_SUM_COUNT; k++) {
659            kTerm2 = k * k;
660            kTerm4 = kTerm2 * kTerm2;
661            increment = (-pi4 * kTerm4 + 3 * MathUtils.PI_SQUARED * kTerm2 * z2) *
662                    FastMath.exp(-MathUtils.PI_SQUARED * kTerm2 / twoZ2);
663            sum2 += increment;
664            if (FastMath.abs(increment) < PG_SUM_RELATIVE_ERROR * FastMath.abs(sum2)) {
665                break;
666            }
667        }
668        if (k == MAXIMUM_PARTIAL_SUM_COUNT) {
669            throw new TooManyIterationsException(MAXIMUM_PARTIAL_SUM_COUNT);
670        }
671        return ret + (sqrtHalfPi / (sqrtN * n)) * (sum / (3240 * z6 * z4) +
672                + sum2 / (108 * z6));
673    }
674
675    /**
676     * Creates {@code H} of size {@code m x m} as described in [1] (see above).
677     *
678     * @param d statistic
679     * @param n sample size
680     * @return H matrix
681     * @throws NumberIsTooLargeException if fractional part is greater than 1.
682     * @throws ArithmeticException if algorithm fails to convert {@code h} to a
683     * {@link BigFraction}.
684     */
685    private FieldDenseMatrix<BigFraction> createExactH(double d,
686                                                       int n) {
687        final int k = (int) Math.ceil(n * d);
688        final int m = 2 * k - 1;
689        final double hDouble = k - n * d;
690        if (hDouble >= 1) {
691            throw new NumberIsTooLargeException(hDouble, 1.0, false);
692        }
693        BigFraction h = null;
694        try {
695            h = BigFraction.from(hDouble, 1e-20, 10000);
696        } catch (final ArithmeticException e1) {
697            try {
698                h = BigFraction.from(hDouble, 1e-10, 10000);
699            } catch (final ArithmeticException e2) {
700                h = BigFraction.from(hDouble, 1e-5, 10000);
701            }
702        }
703        final FieldDenseMatrix<BigFraction> Hdata = FieldDenseMatrix.create(BigFractionField.get(), m, m);
704
705        /*
706         * Start by filling everything with either 0 or 1.
707         */
708        for (int i = 0; i < m; ++i) {
709            for (int j = 0; j < m; ++j) {
710                if (i - j + 1 < 0) {
711                    Hdata.set(i, j, BigFraction.ZERO);
712                } else {
713                    Hdata.set(i, j, BigFraction.ONE);
714                }
715            }
716        }
717
718        /*
719         * Setting up power-array to avoid calculating the same value twice: hPowers[0] = h^1 ...
720         * hPowers[m-1] = h^m
721         */
722        final BigFraction[] hPowers = new BigFraction[m];
723        hPowers[0] = h;
724        for (int i = 1; i < m; i++) {
725            hPowers[i] = h.multiply(hPowers[i - 1]);
726        }
727
728        /*
729         * First column and last row has special values (each other reversed).
730         */
731        for (int i = 0; i < m; ++i) {
732            Hdata.set(i, 0,
733                      Hdata.get(i, 0).subtract(hPowers[i]));
734            Hdata.set(m - 1, i,
735                      Hdata.get(m - 1, i).subtract(hPowers[m - i - 1]));
736        }
737
738        /*
739         * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix should be (1 - 2*h^m +
740         * (2h - 1)^m )/m!" Since 0 <= h < 1, then if h > 1/2 is sufficient to check:
741         */
742        if (h.compareTo(ONE_HALF) > 0) {
743            Hdata.set(m - 1, 0,
744                      Hdata.get(m - 1, 0).add(h.multiply(2).subtract(1).pow(m)));
745        }
746
747        /*
748         * Aside from the first column and last row, the (i, j)-th element is 1/(i - j + 1)! if i -
749         * j + 1 >= 0, else 0. 1's and 0's are already put, so only division with (i - j + 1)! is
750         * needed in the elements that have 1's. There is no need to calculate (i - j + 1)! and then
751         * divide - small steps avoid overflows. Note that i - j + 1 > 0 <=> i + 1 > j instead of
752         * j'ing all the way to m. Also note that it is started at g = 2 because dividing by 1 isn't
753         * really necessary.
754         */
755        for (int i = 0; i < m; ++i) {
756            for (int j = 0; j < i + 1; ++j) {
757                if (i - j + 1 > 0) {
758                    for (int g = 2; g <= i - j + 1; ++g) {
759                        Hdata.set(i, j,
760                                  Hdata.get(i, j).divide(g));
761                    }
762                }
763            }
764        }
765        return Hdata;
766    }
767
768    /***
769     * Creates {@code H} of size {@code m x m} as described in [1] (see above)
770     * using double-precision.
771     *
772     * @param d statistic
773     * @param n sample size
774     * @return H matrix
775     * @throws NumberIsTooLargeException if fractional part is greater than 1
776     */
777    private RealMatrix createRoundedH(double d, int n)
778        throws NumberIsTooLargeException {
779
780        final int k = (int) Math.ceil(n * d);
781        final int m = 2 * k - 1;
782        final double h = k - n * d;
783        if (h >= 1) {
784            throw new NumberIsTooLargeException(h, 1.0, false);
785        }
786        final double[][] Hdata = new double[m][m];
787
788        /*
789         * Start by filling everything with either 0 or 1.
790         */
791        for (int i = 0; i < m; ++i) {
792            for (int j = 0; j < m; ++j) {
793                if (i - j + 1 < 0) {
794                    Hdata[i][j] = 0;
795                } else {
796                    Hdata[i][j] = 1;
797                }
798            }
799        }
800
801        /*
802         * Setting up power-array to avoid calculating the same value twice: hPowers[0] = h^1 ...
803         * hPowers[m-1] = h^m
804         */
805        final double[] hPowers = new double[m];
806        hPowers[0] = h;
807        for (int i = 1; i < m; ++i) {
808            hPowers[i] = h * hPowers[i - 1];
809        }
810
811        /*
812         * First column and last row has special values (each other reversed).
813         */
814        for (int i = 0; i < m; ++i) {
815            Hdata[i][0] = Hdata[i][0] - hPowers[i];
816            Hdata[m - 1][i] -= hPowers[m - i - 1];
817        }
818
819        /*
820         * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix should be (1 - 2*h^m +
821         * (2h - 1)^m )/m!" Since 0 <= h < 1, then if h > 1/2 is sufficient to check:
822         */
823        if (Double.compare(h, 0.5) > 0) {
824            Hdata[m - 1][0] += FastMath.pow(2 * h - 1, m);
825        }
826
827        /*
828         * Aside from the first column and last row, the (i, j)-th element is 1/(i - j + 1)! if i -
829         * j + 1 >= 0, else 0. 1's and 0's are already put, so only division with (i - j + 1)! is
830         * needed in the elements that have 1's. There is no need to calculate (i - j + 1)! and then
831         * divide - small steps avoid overflows. Note that i - j + 1 > 0 <=> i + 1 > j instead of
832         * j'ing all the way to m. Also note that it is started at g = 2 because dividing by 1 isn't
833         * really necessary.
834         */
835        for (int i = 0; i < m; ++i) {
836            for (int j = 0; j < i + 1; ++j) {
837                if (i - j + 1 > 0) {
838                    for (int g = 2; g <= i - j + 1; ++g) {
839                        Hdata[i][j] /= g;
840                    }
841                }
842            }
843        }
844        return MatrixUtils.createRealMatrix(Hdata);
845    }
846
847    /**
848     * Verifies that {@code array} has length at least 2.
849     *
850     * @param array array to test
851     * @throws NullArgumentException if array is null
852     * @throws InsufficientDataException if array is too short
853     */
854    private void checkArray(double[] array) {
855        if (array == null) {
856            throw new NullArgumentException(LocalizedFormats.NULL_NOT_ALLOWED);
857        }
858        if (array.length < 2) {
859            throw new InsufficientDataException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, array.length,
860                                                2);
861        }
862    }
863
864    /**
865     * Computes \( 1 + 2 \sum_{i=1}^\infty (-1)^i e^{-2 i^2 t^2} \) stopping when successive partial
866     * sums are within {@code tolerance} of one another, or when {@code maxIterations} partial sums
867     * have been computed. If the sum does not converge before {@code maxIterations} iterations a
868     * {@link TooManyIterationsException} is thrown.
869     *
870     * @param t argument
871     * @param tolerance Cauchy criterion for partial sums
872     * @param maxIterations maximum number of partial sums to compute
873     * @return Kolmogorov sum evaluated at t
874     * @throws TooManyIterationsException if the series does not converge
875     */
876    public double ksSum(double t, double tolerance, int maxIterations) {
877        if (t == 0.0) {
878            return 0.0;
879        }
880
881        // TODO: for small t (say less than 1), the alternative expansion in part 3 of [1]
882        // from class javadoc should be used.
883
884        final double x = -2 * t * t;
885        int sign = -1;
886        long i = 1;
887        double partialSum = 0.5d;
888        double delta = 1;
889        while (delta > tolerance && i < maxIterations) {
890            delta = FastMath.exp(x * i * i);
891            partialSum += sign * delta;
892            sign *= -1;
893            i++;
894        }
895        if (i == maxIterations) {
896            throw new TooManyIterationsException(maxIterations);
897        }
898        return partialSum * 2;
899    }
900
901    /**
902     * Given a d-statistic in the range [0, 1] and the two sample sizes n and m,
903     * an integral d-statistic in the range [0, n*m] is calculated, that can be used for
904     * comparison with other integral d-statistics. Depending whether {@code strict} is
905     * {@code true} or not, the returned value divided by (n*m) is greater than
906     * (resp greater than or equal to) the given d value (allowing some tolerance).
907     *
908     * @param d a d-statistic in the range [0, 1]
909     * @param n first sample size
910     * @param m second sample size
911     * @param strict whether the returned value divided by (n*m) is allowed to be equal to d
912     * @return the integral d-statistic in the range [0, n*m]
913     */
914    private static long calculateIntegralD(double d, int n, int m, boolean strict) {
915        final double tol = 1e-12;  // d-values within tol of one another are considered equal
916        long nm = n * (long)m;
917        long upperBound = (long)FastMath.ceil((d - tol) * nm);
918        long lowerBound = (long)FastMath.floor((d + tol) * nm);
919        if (strict && lowerBound == upperBound) {
920            return upperBound + 1l;
921        }
922        else {
923            return upperBound;
924        }
925    }
926
927    /**
928     * Computes \(P(D_{n,m} &gt; d)\) if {@code strict} is {@code true}; otherwise \(P(D_{n,m} \ge
929     * d)\), where \(D_{n,m}\) is the 2-sample Kolmogorov-Smirnov statistic. See
930     * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
931     * <p>
932     * The returned probability is exact, implemented by unwinding the recursive function
933     * definitions presented in [4] (class javadoc).
934     * </p>
935     *
936     * @param d D-statistic value
937     * @param n first sample size
938     * @param m second sample size
939     * @param strict whether or not the probability to compute is expressed as a strict inequality
940     * @return probability that a randomly selected m-n partition of m + n generates \(D_{n,m}\)
941     *         greater than (resp. greater than or equal to) {@code d}
942     */
943    public double exactP(double d, int n, int m, boolean strict) {
944       return 1 - n(m, n, m, n, calculateIntegralD(d, m, n, strict), strict) /
945           BinomialCoefficientDouble.value(n + m, m);
946    }
947
948    /**
949     * Uses the Kolmogorov-Smirnov distribution to approximate \(P(D_{n,m} &gt; d)\) where \(D_{n,m}\)
950     * is the 2-sample Kolmogorov-Smirnov statistic. See
951     * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
952     * <p>
953     * Specifically, what is returned is \(1 - k(d \sqrt{mn / (m + n)})\) where \(k(t) = 1 + 2
954     * \sum_{i=1}^\infty (-1)^i e^{-2 i^2 t^2}\). See {@link #ksSum(double, double, int)} for
955     * details on how convergence of the sum is determined.
956     * </p>
957     *
958     * @param d D-statistic value
959     * @param n first sample size
960     * @param m second sample size
961     * @return approximate probability that a randomly selected m-n partition of m + n generates
962     *         \(D_{n,m}\) greater than {@code d}
963     */
964    public double approximateP(double d, int n, int m) {
965        final double dm = m;
966        final double dn = n;
967        return 1 - ksSum(d * FastMath.sqrt((dm * dn) / (dm + dn)),
968                         KS_SUM_CAUCHY_CRITERION, MAXIMUM_PARTIAL_SUM_COUNT);
969    }
970
971    /**
972     * Fills a boolean array randomly with a fixed number of {@code true} values.
973     * The method uses a simplified version of the Fisher-Yates shuffle algorithm.
974     * By processing first the {@code true} values followed by the remaining {@code false} values
975     * less random numbers need to be generated. The method is optimized for the case
976     * that the number of {@code true} values is larger than or equal to the number of
977     * {@code false} values.
978     *
979     * @param b boolean array
980     * @param numberOfTrueValues number of {@code true} values the boolean array should finally have
981     * @param rng random data generator
982     */
983    private static void fillBooleanArrayRandomlyWithFixedNumberTrueValues(final boolean[] b, final int numberOfTrueValues, final UniformRandomProvider rng) {
984        Arrays.fill(b, true);
985        for (int k = numberOfTrueValues; k < b.length; k++) {
986            final int r = rng.nextInt(k + 1);
987            b[(b[r]) ? r : k] = false;
988        }
989    }
990
991    /**
992     * Uses Monte Carlo simulation to approximate \(P(D_{n,m} &gt; d)\) where \(D_{n,m}\) is the
993     * 2-sample Kolmogorov-Smirnov statistic. See
994     * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
995     * <p>
996     * The simulation generates {@code iterations} random partitions of {@code m + n} into an
997     * {@code n} set and an {@code m} set, computing \(D_{n,m}\) for each partition and returning
998     * the proportion of values that are greater than {@code d}, or greater than or equal to
999     * {@code d} if {@code strict} is {@code false}.
1000     * </p>
1001     *
1002     * @param d D-statistic value.
1003     * @param n First sample size.
1004     * @param m Second sample size.
1005     * @param iterations Number of random partitions to generate.
1006     * @param strict whether or not the probability to compute is expressed as a strict inequality
1007     * @param rng RNG used for generating the partitions.
1008     * @return proportion of randomly generated m-n partitions of m + n that result in \(D_{n,m}\)
1009     * greater than (resp. greater than or equal to) {@code d}.
1010     */
1011    public double monteCarloP(final double d,
1012                              final int n,
1013                              final int m,
1014                              final boolean strict,
1015                              final int iterations,
1016                              UniformRandomProvider rng) {
1017        return integralMonteCarloP(calculateIntegralD(d, n, m, strict), n, m, iterations, rng);
1018    }
1019
1020    /**
1021     * Uses Monte Carlo simulation to approximate \(P(D_{n,m} >= d / (n * m))\)
1022     * where \(D_{n,m}\) is the 2-sample Kolmogorov-Smirnov statistic.
1023     * <p>
1024     * Here {@code d} is the D-statistic represented as long value.
1025     * The real D-statistic is obtained by dividing {@code d} by {@code n * m}.
1026     * See also {@link #monteCarloP(double,int,int,boolean,int,UniformRandomProvider)}.
1027     *
1028     * @param d Integral D-statistic.
1029     * @param n First sample size.
1030     * @param m Second sample size.
1031     * @param iterations Number of random partitions to generate.
1032     * @param rng RNG used for generating the partitions.
1033     * @return proportion of randomly generated m-n partitions of m + n that result in \(D_{n,m}\)
1034     * greater than or equal to {@code d / (n * m))}.
1035     */
1036    private double integralMonteCarloP(final long d,
1037                                       final int n,
1038                                       final int m,
1039                                       final int iterations,
1040                                       UniformRandomProvider rng) {
1041        // ensure that nn is always the max of (n, m) to require fewer random numbers
1042        final int nn = FastMath.max(n, m);
1043        final int mm = FastMath.min(n, m);
1044        final int sum = nn + mm;
1045
1046        int tail = 0;
1047        final boolean b[] = new boolean[sum];
1048        for (int i = 0; i < iterations; i++) {
1049            fillBooleanArrayRandomlyWithFixedNumberTrueValues(b, nn, rng);
1050            long curD = 0l;
1051            for(int j = 0; j < b.length; ++j) {
1052                if (b[j]) {
1053                    curD += mm;
1054                    if (curD >= d) {
1055                        tail++;
1056                        break;
1057                    }
1058                } else {
1059                    curD -= nn;
1060                    if (curD <= -d) {
1061                        tail++;
1062                        break;
1063                    }
1064                }
1065            }
1066        }
1067        return (double) tail / iterations;
1068    }
1069
1070    /**
1071     * If there are no ties in the combined dataset formed from x and y,
1072     * this method is a no-op.
1073     * If there are ties, a uniform random deviate in
1074     * is added to each value in x and y, and this method overwrites the
1075     * data in x and y with the jittered values.
1076     *
1077     * @param x First sample.
1078     * @param y Second sample.
1079     * @throw NotANumberException if any of the input arrays contain
1080     * a NaN value.
1081     */
1082    private static void fixTies(double[] x, double[] y) {
1083        if (hasTies(x, y)) {
1084            // Add jitter using a fixed seed (so same arguments always give same results),
1085            // low-initialization-overhead generator.
1086            final UniformRandomProvider rng = RandomSource.create(RandomSource.SPLIT_MIX_64, 876543217L);
1087
1088            // It is theoretically possible that jitter does not break ties, so repeat
1089            // until all ties are gone.  Bound the loop and throw MIE if bound is exceeded.
1090            int fixedRangeTrials = 10;
1091            int jitterRange = 10;
1092            int ct = 0;
1093            boolean ties = true;
1094            do {
1095                // Nudge the data using a fixed range.
1096                jitter(x, rng, jitterRange);
1097                jitter(y, rng, jitterRange);
1098                ties = hasTies(x, y);
1099                ++ct;
1100            } while (ties && ct < fixedRangeTrials);
1101
1102            if (ties) {
1103                throw new MaxCountExceededException(fixedRangeTrials);
1104            }
1105        }
1106    }
1107
1108    /**
1109     * Returns true iff there are ties in the combined sample formed from
1110     * x and y.
1111     *
1112     * @param x First sample.
1113     * @param y Second sample.
1114     * @return true if x and y together contain ties.
1115     * @throw InsufficientDataException if the input arrays only contain infinite values.
1116     * @throw NotANumberException if any of the input arrays contain a NaN value.
1117     */
1118    private static boolean hasTies(double[] x, double[] y) {
1119       final double[] values = MathArrays.unique(MathArrays.concatenate(x, y));
1120
1121       // "unique" moves NaN to the head of the output array.
1122       if (Double.isNaN(values[0])) {
1123           throw new NotANumberException();
1124       }
1125
1126       int infCount = 0;
1127       for (double v : values) {
1128           if (Double.isInfinite(v)) {
1129               ++infCount;
1130           }
1131       }
1132       if (infCount == values.length) {
1133           throw new InsufficientDataException();
1134       }
1135
1136       if (values.length == x.length + y.length) {
1137           return false;  // There are no ties.
1138       }
1139
1140       return true;
1141    }
1142
1143    /**
1144     * Adds random jitter to {@code data} using deviates sampled from {@code dist}.
1145     * <p>
1146     * Note that jitter is applied in-place - i.e., the array
1147     * values are overwritten with the result of applying jitter.</p>
1148     *
1149     * @param data input/output data array - entries overwritten by the method
1150     * @param rng probability distribution to sample for jitter values
1151     * @param ulp ulp used when generating random numbers
1152     * @throws NullPointerException if either of the parameters is null
1153     */
1154    private static void jitter(double[] data,
1155                               UniformRandomProvider rng,
1156                               int ulp) {
1157        final int range = ulp * 2;
1158        for (int i = 0; i < data.length; i++) {
1159            final double orig = data[i];
1160            if (Double.isInfinite(orig)) {
1161                continue; // Avoid NaN creation.
1162            }
1163
1164            final int rand = rng.nextInt(range) - ulp;
1165            final double ulps = Math.ulp(orig);
1166            final double r = rand * ulps;
1167
1168            data[i] += r;
1169        }
1170    }
1171
1172    /**
1173     * The function C(i, j) defined in [4] (class javadoc), formula (5.5).
1174     * defined to return 1 if |i/n - j/m| <= c; 0 otherwise. Here c is scaled up
1175     * and recoded as a long to avoid rounding errors in comparison tests, so what
1176     * is actually tested is |im - jn| <= cmn.
1177     *
1178     * @param i first path parameter
1179     * @param j second path paramter
1180     * @param m first sample size
1181     * @param n second sample size
1182     * @param cmn integral D-statistic (see {@link #calculateIntegralD(double, int, int, boolean)})
1183     * @param strict whether or not the null hypothesis uses strict inequality
1184     * @return C(i,j) for given m, n, c
1185     */
1186    private static int c(int i, int j, int m, int n, long cmn, boolean strict) {
1187        if (strict) {
1188            return FastMath.abs(i*(long)n - j*(long)m) <= cmn ? 1 : 0;
1189        }
1190        return FastMath.abs(i*(long)n - j*(long)m) < cmn ? 1 : 0;
1191    }
1192
1193    /**
1194     * The function N(i, j) defined in [4] (class javadoc).
1195     * Returns the number of paths over the lattice {(i,j) : 0 <= i <= n, 0 <= j <= m}
1196     * from (0,0) to (i,j) satisfying C(h,k, m, n, c) = 1 for each (h,k) on the path.
1197     * The return value is integral, but subject to overflow, so it is maintained and
1198     * returned as a double.
1199     *
1200     * @param i first path parameter
1201     * @param j second path parameter
1202     * @param m first sample size
1203     * @param n second sample size
1204     * @param cnm integral D-statistic (see {@link #calculateIntegralD(double, int, int, boolean)})
1205     * @param strict whether or not the null hypothesis uses strict inequality
1206     * @return number or paths to (i, j) from (0,0) representing D-values as large as c for given m, n
1207     */
1208    private static double n(int i, int j, int m, int n, long cnm, boolean strict) {
1209        /*
1210         * Unwind the recursive definition given in [4].
1211         * Compute n(1,1), n(1,2)...n(2,1), n(2,2)... up to n(i,j), one row at a time.
1212         * When n(i,*) are being computed, lag[] holds the values of n(i - 1, *).
1213         */
1214        final double[] lag = new double[n];
1215        double last = 0;
1216        for (int k = 0; k < n; k++) {
1217            lag[k] = c(0, k + 1, m, n, cnm, strict);
1218        }
1219        for (int k = 1; k <= i; k++) {
1220            last = c(k, 0, m, n, cnm, strict);
1221            for (int l = 1; l <= j; l++) {
1222                lag[l - 1] = c(k, l, m, n, cnm, strict) * (last + lag[l - 1]);
1223                last = lag[l - 1];
1224            }
1225        }
1226        return last;
1227    }
1228}