You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ps...@apache.org on 2014/02/27 00:14:08 UTC
svn commit: r1572335 - in /commons/proper/math/trunk/src: changes/
main/java/org/apache/commons/math3/distribution/
main/java/org/apache/commons/math3/stat/inference/ test/R/
test/java/org/apache/commons/math3/stat/inference/
Author: psteitz
Date: Wed Feb 26 23:14:08 2014
New Revision: 1572335
URL: http://svn.apache.org/r1572335
Log:
Added KolmogorovSmirnovTest, deprecated KolmogorovSmirnovDistributio. JIRA: MATH-437
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java
- copied, changed from r1533924, commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java
commons/proper/math/trunk/src/test/R/KolmogorovSmirnovTestCases.R
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java
- copied, changed from r1488409, commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistributionTest.java
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java
commons/proper/math/trunk/src/test/R/testAll
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1572335&r1=1572334&r2=1572335&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Wed Feb 26 23:14:08 2014
@@ -51,6 +51,9 @@ If the output is not quite correct, chec
</properties>
<body>
<release version="3.3" date="TBD" description="TBD">
+ <action dev="psteitz" type="update" issue="MATH-437">
+ Added KolmogorovSmirnovTest class, deprecating KolmogorovSmirnovDistribution.
+ </action>
<action dev="luc" type="add" issue="MATH-1101">
Improved documentation of QR decomposition handling of singular matrices.
</action>
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java?rev=1572335&r1=1572334&r2=1572335&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java Wed Feb 26 23:14:08 2014
@@ -69,6 +69,8 @@ import org.apache.commons.math3.util.Fas
* @see <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test">
* Kolmogorov-Smirnov test (Wikipedia)</a>
* @version $Id$
+ * @deprecated to be removed in version 4.0 -
+ * use {@link org.apache.commmons.math3.stat.inference.KolmogorovSmirnovTest}
*/
public class KolmogorovSmirnovDistribution implements Serializable {
Copied: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java (from r1533924, commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java?p2=commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java&p1=commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java&r1=1533924&r2=1572335&rev=1572335&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java Wed Feb 26 23:14:08 2014
@@ -15,14 +15,19 @@
* limitations under the License.
*/
-package org.apache.commons.math3.distribution;
+package org.apache.commons.math3.stat.inference;
-import java.io.Serializable;
import java.math.BigDecimal;
+import java.util.Arrays;
+import java.util.Iterator;
+import org.apache.commons.math3.distribution.RealDistribution;
+import org.apache.commons.math3.exception.InsufficientDataException;
import org.apache.commons.math3.exception.MathArithmeticException;
-import org.apache.commons.math3.exception.NotStrictlyPositiveException;
+import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
+import org.apache.commons.math3.exception.OutOfRangeException;
+import org.apache.commons.math3.exception.TooManyIterationsException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.fraction.BigFraction;
import org.apache.commons.math3.fraction.BigFractionField;
@@ -31,171 +36,413 @@ import org.apache.commons.math3.linear.A
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.FieldMatrix;
import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.random.RandomGenerator;
+import org.apache.commons.math3.random.Well19937c;
+import org.apache.commons.math3.util.CombinatoricsUtils;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
/**
- * Implementation of the Kolmogorov-Smirnov distribution.
- *
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test">
+ * Kolmogorov-Smirnov (K-S) test</a> for equality of continuous distributions.
+ * <p>
+ * The K-S test uses a statistic based on the maximum deviation of the empirical distribution of
+ * sample data points from the distribution expected under the null hypothesis. For one-sample tests
+ * evaluating the null hypothesis that a set of sample data points follow a given distribution, the
+ * test statistic is \(D_n=\sup_x |F_n(x)-F(x)|\), where \(F\) is the expected distribution and
+ * \(F_n\) is the empirical distribution of the \(n\) sample data points. The distribution of
+ * \(D_n\) is estimated using a method based on [1] with certain quick decisions for extreme values
+ * given in [2].
+ * </p>
* <p>
- * Treats the distribution of the two-sided {@code P(D_n < d)} where
- * {@code D_n = sup_x |G(x) - G_n (x)|} for the theoretical cdf {@code G} and
- * the empirical cdf {@code G_n}.
+ * Two-sample tests are also supported, evaluating the null hypothesis that the two samples
+ * {@code x} and {@code y} come from the same underlying distribution. In this case, the test
+ * statistic is \(D_{n,m}=\sup_t | F_n(t)-F_m(t)|\) where \(n\) is the length of {@code x}, \(m\) is
+ * the length of {@code y}, \(F_n\) is the empirical distribution that puts mass \(1/n\) at each of
+ * the values in {@code x} and \(F_m\) is the empirical distribution of the {@code y} values. The
+ * default 2-sample test method, {@link #kolmogorovSmirnovTest(double[], double[])} works as
+ * follows:
+ * <ul>
+ * <li>For very small samples (where the product of the sample sizes is less than
+ * {@value #SMALL_SAMPLE_PRODUCT}), the exact distribution is used to compute the p-value for the
+ * 2-sample test.</li>
+ * <li>For mid-size samples (product of sample sizes greater than or equal to
+ * {@value #SMALL_SAMPLE_PRODUCT} but less than {@value #LARGE_SAMPLE_PRODUCT}), Monte Carlo
+ * simulation is used to compute the p-value. The simulation randomly generates partitions of \(m +
+ * n\) into an \(m\)-set and an \(n\)-set and reports the proportion that give \(D\) values
+ * exceeding the observed value.</li>
+ * <li>When the product of the sample sizes exceeds {@value #LARGE_SAMPLE_PRODUCT}, the asymptotic
+ * distribution of \(D_{n,m}\) is used. See {@link #approximateP(double, int, int)} for details on
+ * the approximation.</li>
+ * </ul>
* </p>
* <p>
- * This implementation is based on [1] with certain quick decisions for extreme
- * values given in [2].
+ * In the two-sample case, \(D_{n,m}\) has a discrete distribution. This makes the p-value
+ * associated with the null hypothesis \(H_0 : D_{n,m} \ge d \) differ from \(H_0 : D_{n,m} > d \)
+ * by the mass of the observed value \(d\). To distinguish these, the two-sample tests use a boolean
+ * {@code strict} parameter. This parameter is ignored for large samples.
* </p>
* <p>
- * In short, when wanting to evaluate {@code P(D_n < d)}, the method in [1] is
- * to write {@code d = (k - h) / n} for positive integer {@code k} and
- * {@code 0 <= h < 1}. Then {@code P(D_n < d) = (n! / n^n) * t_kk}, where
- * {@code t_kk} is the {@code (k, k)}'th entry in the special matrix
- * {@code H^n}, i.e. {@code H} to the {@code n}'th power.
+ * The methods used by the 2-sample default implementation are also exposed directly:
+ * <ul>
+ * <li>{@link #exactP(double, int, int, boolean)} computes exact 2-sample p-values</li>
+ * <li>{@link #monteCarloP(double, int, int, boolean, int)} computes 2-sample p-values by Monte
+ * Carlo simulation</li>
+ * <li>{@link #approximateP(double, int, int)} uses the asymptotic distribution The {@code boolean}
+ * arguments in the first two methods allow the probability used to estimate the p-value to be
+ * expressed using strict or non-strict inequality. See
+ * {@link #kolmogorovSmirnovTest(double[], double[], boolean)}.</li>
+ * </ul>
* </p>
* <p>
* References:
* <ul>
- * <li>[1] <a href="http://www.jstatsoft.org/v08/i18/">
- * Evaluating Kolmogorov's Distribution</a> by George Marsaglia, Wai
- * Wan Tsang, and Jingbo Wang</li>
- * <li>[2] <a href="http://www.jstatsoft.org/v39/i11/">
- * Computing the Two-Sided Kolmogorov-Smirnov Distribution</a> by Richard Simard
- * and Pierre L'Ecuyer</li>
+ * <li>[1] <a href="http://www.jstatsoft.org/v08/i18/"> Evaluating Kolmogorov's Distribution</a> by
+ * George Marsaglia, Wai Wan Tsang, and Jingbo Wang</li>
+ * <li>[2] <a href="http://www.jstatsoft.org/v39/i11/"> Computing the Two-Sided Kolmogorov-Smirnov
+ * Distribution</a> by Richard Simard and Pierre L'Ecuyer</li>
* </ul>
- * Note that [1] contains an error in computing h, refer to
- * <a href="https://issues.apache.org/jira/browse/MATH-437">MATH-437</a> for details.
+ * <br/>
+ * Note that [1] contains an error in computing h, refer to <a
+ * href="https://issues.apache.org/jira/browse/MATH-437">MATH-437</a> for details.
* </p>
*
- * @see <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test">
- * Kolmogorov-Smirnov test (Wikipedia)</a>
+ * @since 3.3
* @version $Id$
*/
-public class KolmogorovSmirnovDistribution implements Serializable {
+public class KolmogorovSmirnovTest {
+
+ /**
+ * Bound on the number of partial sums in {@link #ksSum(double, double, int)}
+ */
+ protected static final int MAXIMUM_PARTIAL_SUM_COUNT = 100000;
+
+ /** Convergence criterion for {@link #ksSum(double, double, int)} */
+ protected static final double KS_SUM_CAUCHY_CRITERION = 1E-20;
+
+ /** When product of sample sizes is less than this value, 2-sample K-S test is exact */
+ protected static final int SMALL_SAMPLE_PRODUCT = 200;
+
+ /**
+ * When product of sample sizes exceeds this value, 2-sample K-S test uses asymptotic
+ * distribution for strict inequality p-value.
+ */
+ protected static final int LARGE_SAMPLE_PRODUCT = 10000;
+
+ /** Default number of iterations used by {@link #monteCarloP(double, int, int, boolean, int)} */
+ protected static final int MONTE_CARLO_ITERATIONS = 1000000;
+
+ /** Random data generator used by {@link #monteCarloP(double, int, int, boolean, int)} */
+ private final RandomGenerator rng;
+
+ /**
+ * Construct a KolmogorovSmirnovTest instance with a default random data generator.
+ */
+ public KolmogorovSmirnovTest() {
+ rng = new Well19937c();
+ }
+
+ /**
+ * Construct a KolmogorovSmirnovTest with the provided random data generator.
+ *
+ * @param rng random data generator used by {@link #monteCarloP(double, int, int, boolean, int)}
+ */
+ public KolmogorovSmirnovTest(RandomGenerator rng) {
+ this.rng = rng;
+ }
+
+ /**
+ * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a one-sample <a
+ * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
+ * evaluating the null hypothesis that {@code data} conforms to {@code distribution}. If
+ * {@code exact} is true, the distribution used to compute the p-value is computed using
+ * extended precision. See {@link #cdfExact(double, int)}.
+ *
+ * @param distribution reference distribution
+ * @param data sample being being evaluated
+ * @param exact whether or not to force exact computation of the p-value
+ * @return the p-value associated with the null hypothesis that {@code data} is a sample from
+ * {@code distribution}
+ * @throws InsufficientDataException if {@code data} does not have length at least 2
+ * @throws NullArgumentException if {@code data} is null
+ */
+ public double kolmogorovSmirnovTest(RealDistribution distribution, double[] data, boolean exact) {
+ return 1d - cdf(kolmogorovSmirnovStatistic(distribution, data), data.length, exact);
+ }
+
+ /**
+ * Computes the one-sample Kolmogorov-Smirnov test statistic, \(D_n=\sup_x |F_n(x)-F(x)|\) where
+ * \(F\) is the distribution (cdf) function associated with {@code distribution}, \(n\) is the
+ * length of {@code data} and \(F_n\) is the empirical distribution that puts mass \(1/n\) at
+ * each of the values in {@code data}.
+ *
+ * @param distribution reference distribution
+ * @param data sample being evaluated
+ * @return Kolmogorov-Smirnov statistic \(D_n\)
+ * @throws InsufficientDataException if {@code data} does not have length at least 2
+ * @throws NullArgumentException if {@code data} is null
+ */
+ public double kolmogorovSmirnovStatistic(RealDistribution distribution, double[] data) {
+ checkArray(data);
+ final int n = data.length;
+ final double nd = n;
+ final double[] dataCopy = new double[n];
+ System.arraycopy(data, 0, dataCopy, 0, n);
+ Arrays.sort(dataCopy);
+ double d = 0d;
+ for (int i = 1; i <= n; i++) {
+ final double yi = distribution.cumulativeProbability(dataCopy[i - 1]);
+ final double currD = FastMath.max(yi - (i - 1) / nd, i / nd - yi);
+ if (currD > d) {
+ d = currD;
+ }
+ }
+ return d;
+ }
- /** Serializable version identifier. */
- private static final long serialVersionUID = -4670676796862967187L;
+ /**
+ * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a two-sample <a
+ * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
+ * evaluating the null hypothesis that {@code x} and {@code y} are samples drawn from the same
+ * probability distribution. Specifically, what is returned is an estimate of the probability
+ * that the {@link #kolmogorovSmirnovStatistic(double[], double[])} associated with a randomly
+ * selected partition of the combined sample into subsamples of sizes {@code x.length} and
+ * {@code y.length} will strictly exceed (if {@code strict} is {@code true}) or be at least as
+ * large as {@code strict = false}) as {@code kolmogorovSmirnovStatistic(x, y)}.
+ * <ul>
+ * <li>For very small samples (where the product of the sample sizes is less than
+ * {@value #SMALL_SAMPLE_PRODUCT}), the exact distribution is used to compute the p-value. This
+ * is accomplished by enumerating all partitions of the combined sample into two subsamples of
+ * the respective sample sizes, computing \(D_{n,m}\) for each partition and returning the
+ * proportion of partitions that give \(D\) values exceeding the observed value.</li>
+ * <li>For mid-size samples (product of sample sizes greater than or equal to
+ * {@value #SMALL_SAMPLE_PRODUCT} but less than {@value #LARGE_SAMPLE_PRODUCT}), Monte Carlo
+ * simulation is used to compute the p-value. The simulation randomly generates partitions and
+ * reports the proportion that give \(D\) values exceeding the observed value.</li>
+ * <li>When the product of the sample sizes exceeds {@value #LARGE_SAMPLE_PRODUCT}, the
+ * asymptotic distribution of \(D_{n,m}\) is used. See {@link #approximateP(double, int, int)}
+ * for details on the approximation.</li>
+ * </ul>
+ *
+ * @param x first sample dataset
+ * @param y second sample dataset
+ * @param strict whether or not the probability to compute is expressed as a strict inequality
+ * (ignored for large samples)
+ * @return p-value associated with the null hypothesis that {@code x} and {@code y} represent
+ * samples from the same distribution
+ * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
+ * least 2
+ * @throws NullArgumentException if either {@code x} or {@code y} is null
+ */
+ public double kolmogorovSmirnovTest(double[] x, double[] y, boolean strict) {
+ if (x.length * y.length < SMALL_SAMPLE_PRODUCT) {
+ return exactP(kolmogorovSmirnovStatistic(x, y), x.length, y.length, strict);
+ }
+ if (x.length * y.length < LARGE_SAMPLE_PRODUCT) {
+ return monteCarloP(kolmogorovSmirnovStatistic(x, y), x.length, y.length, strict, MONTE_CARLO_ITERATIONS);
+ }
+ return approximateP(kolmogorovSmirnovStatistic(x, y), x.length, y.length);
+ }
- /** Number of observations. */
- private int n;
+ /**
+ * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a two-sample <a
+ * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
+ * evaluating the null hypothesis that {@code x} and {@code y} are samples drawn from the same
+ * probability distribution. Assumes the strict form of the inequality used to compute the
+ * p-value. See {@link #kolmogorovSmirnovTest(RealDistribution, double[], boolean)}.
+ *
+ * @param x first sample dataset
+ * @param y second sample dataset
+ * @return p-value associated with the null hypothesis that {@code x} and {@code y} represent
+ * samples from the same distribution
+ * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
+ * least 2
+ * @throws NullArgumentException if either {@code x} or {@code y} is null
+ */
+ public double kolmogorovSmirnovTest(double[] x, double[] y) {
+ return kolmogorovSmirnovTest(x, y, true);
+ }
/**
- * @param n Number of observations
- * @throws NotStrictlyPositiveException if {@code n <= 0}
+ * Computes the two-sample Kolmogorov-Smirnov test statistic, \(D_{n,m}=\sup_x |F_n(x)-F_m(x)|\)
+ * where \(n\) is the length of {@code x}, \(m\) is the length of {@code y}, \(F_n\) is the
+ * empirical distribution that puts mass \(1/n\) at each of the values in {@code x} and \(F_m\)
+ * is the empirical distribution of the {@code y} values.
+ *
+ * @param x first sample
+ * @param y second sample
+ * @return test statistic \(D_{n,m}\) used to evaluate the null hypothesis that {@code x} and
+ * {@code y} represent samples from the same underlying distribution
+ * @throws InsufficientDataException if either {@code x} or {@code y} does not have length at
+ * least 2
+ * @throws NullArgumentException if either {@code x} or {@code y} is null
*/
- public KolmogorovSmirnovDistribution(int n)
- throws NotStrictlyPositiveException {
- if (n <= 0) {
- throw new NotStrictlyPositiveException(LocalizedFormats.NOT_POSITIVE_NUMBER_OF_SAMPLES, n);
+ public double kolmogorovSmirnovStatistic(double[] x, double[] y) {
+ checkArray(x);
+ checkArray(y);
+ // Copy and sort the sample arrays
+ final double[] sx = MathArrays.copyOf(x);
+ final double[] sy = MathArrays.copyOf(y);
+ Arrays.sort(sx);
+ Arrays.sort(sy);
+ final int n = sx.length;
+ final int m = sy.length;
+
+ // Find the max difference between cdf_x and cdf_y
+ double supD = 0d;
+ // First walk x points
+ for (int i = 0; i < n; i++) {
+ final double cdf_x = (i + 1d) / n;
+ final int yIndex = Arrays.binarySearch(sy, sx[i]);
+ final double cdf_y = yIndex >= 0 ? (yIndex + 1d) / m : (-yIndex - 1d) / m;
+ final double curD = FastMath.abs(cdf_x - cdf_y);
+ if (curD > supD) {
+ supD = curD;
+ }
+ }
+ // Now look at y
+ for (int i = 0; i < m; i++) {
+ final double cdf_y = (i + 1d) / m;
+ final int xIndex = Arrays.binarySearch(sx, sy[i]);
+ final double cdf_x = xIndex >= 0 ? (xIndex + 1d) / n : (-xIndex - 1d) / n;
+ final double curD = FastMath.abs(cdf_x - cdf_y);
+ if (curD > supD) {
+ supD = curD;
+ }
}
+ return supD;
+ }
+
+ /**
+ * Computes the <i>p-value</i>, or <i>observed significance level</i>, of a one-sample <a
+ * href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov test</a>
+ * evaluating the null hypothesis that {@code data} conforms to {@code distribution}.
+ *
+ * @param distribution reference distribution
+ * @param data sample being being evaluated
+ * @return the p-value associated with the null hypothesis that {@code data} is a sample from
+ * {@code distribution}
+ * @throws InsufficientDataException if {@code data} does not have length at least 2
+ * @throws NullArgumentException if {@code data} is null
+ */
+ public double kolmogorovSmirnovTest(RealDistribution distribution, double[] data) {
+ return kolmogorovSmirnovTest(distribution, data, false);
+ }
- this.n = n;
+ /**
+ * Performs a <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test"> Kolmogorov-Smirnov
+ * test</a> evaluating the null hypothesis that {@code data} conforms to {@code distribution}.
+ *
+ * @param distribution reference distribution
+ * @param data sample being being evaluated
+ * @param alpha significance level of the test
+ * @return true iff the null hypothesis that {@code data} is a sample from {@code distribution}
+ * can be rejected with confidence 1 - {@code alpha}
+ * @throws InsufficientDataException if {@code data} does not have length at least 2
+ * @throws NullArgumentException if {@code data} is null
+ */
+ public boolean kolmogorovSmirnovTest(RealDistribution distribution, double[] data, double alpha) {
+ if ((alpha <= 0) || (alpha > 0.5)) {
+ throw new OutOfRangeException(LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL, alpha, 0, 0.5);
+ }
+ return kolmogorovSmirnovTest(distribution, data) < alpha;
}
/**
- * Calculates {@code P(D_n < d)} using method described in [1] with quick
- * decisions for extreme values given in [2] (see above). The result is not
- * exact as with
- * {@link KolmogorovSmirnovDistribution#cdfExact(double)} because
- * calculations are based on {@code double} rather than
- * {@link org.apache.commons.math3.fraction.BigFraction}.
+ * Calculates \(P(D_n < d)\) using the method described in [1] with quick decisions for extreme
+ * values given in [2] (see above). The result is not exact as with
+ * {@link KolmogorovSmirnovDistribution#cdfExact(double)} because calculations are based on
+ * {@code double} rather than {@link org.apache.commons.math3.fraction.BigFraction}.
*
* @param d statistic
- * @return the two-sided probability of {@code P(D_n < d)}
- * @throws MathArithmeticException if algorithm fails to convert {@code h}
- * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
- * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
- */
- public double cdf(double d) throws MathArithmeticException {
- return this.cdf(d, false);
+ * @param n sample size
+ * @return \(P(D_n < d)\)
+ * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\) for integer {@code k, m} and \(0 \le h < 1\)
+ */
+ public double cdf(double d, int n)
+ throws MathArithmeticException {
+ return cdf(d, n, false);
}
/**
- * Calculates {@code P(D_n < d)} using method described in [1] with quick
- * decisions for extreme values given in [2] (see above). The result is
- * exact in the sense that BigFraction/BigReal is used everywhere at the
- * expense of very slow execution time. Almost never choose this in real
- * applications unless you are very sure; this is almost solely for
- * verification purposes. Normally, you would choose
- * {@link KolmogorovSmirnovDistribution#cdf(double)}
+ * Calculates {@code P(D_n < d)}. The result is exact in the sense that BigFraction/BigReal is
+ * used everywhere at the expense of very slow execution time. Almost never choose this in real
+ * applications unless you are very sure; this is almost solely for verification purposes.
+ * Normally, you would choose {@link KolmogorovSmirnovDistribution#cdf(double)}. See the class
+ * javadoc for definitions and algorithm description.
*
* @param d statistic
- * @return the two-sided probability of {@code P(D_n < d)}
- * @throws MathArithmeticException if algorithm fails to convert {@code h}
- * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
- * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
+ * @param n sample size
+ * @return \(P(D_n < d)\)
+ * @throws MathArithmeticException if the algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\) for integer {@code k, m} and \(0 \le h < 1\)
*/
- public double cdfExact(double d) throws MathArithmeticException {
- return this.cdf(d, true);
+ public double cdfExact(double d, int n)
+ throws MathArithmeticException {
+ return cdf(d, n, true);
}
/**
- * Calculates {@code P(D_n < d)} using method described in [1] with quick
- * decisions for extreme values given in [2] (see above).
+ * Calculates {@code P(D_n < d)} using method described in [1] with quick decisions for extreme
+ * values given in [2] (see above).
*
* @param d statistic
+ * @param n sample size
* @param exact whether the probability should be calculated exact using
- * {@link org.apache.commons.math3.fraction.BigFraction} everywhere at the
- * expense of very slow execution time, or if {@code double} should be used
- * convenient places to gain speed. Almost never choose {@code true} in real
- * applications unless you are very sure; {@code true} is almost solely for
- * verification purposes.
- * @return the two-sided probability of {@code P(D_n < d)}
- * @throws MathArithmeticException if algorithm fails to convert {@code h}
- * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
- * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
+ * {@link org.apache.commons.math3.fraction.BigFraction} everywhere at the expense of
+ * very slow execution time, or if {@code double} should be used convenient places to
+ * gain speed. Almost never choose {@code true} in real applications unless you are very
+ * sure; {@code true} is almost solely for verification purposes.
+ * @return \(P(D_n < d)\)
+ * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\) for integer {@code k, m} and \(0 \le h < 1\).
*/
- public double cdf(double d, boolean exact) throws MathArithmeticException {
+ public double cdf(double d, int n, boolean exact)
+ throws MathArithmeticException {
final double ninv = 1 / ((double) n);
final double ninvhalf = 0.5 * ninv;
if (d <= ninvhalf) {
-
return 0;
-
} else if (ninvhalf < d && d <= ninv) {
-
double res = 1;
- double f = 2 * d - ninv;
-
+ final double f = 2 * d - ninv;
// n! f^n = n*f * (n-1)*f * ... * 1*x
for (int i = 1; i <= n; ++i) {
res *= i * f;
}
-
return res;
-
} else if (1 - ninv <= d && d < 1) {
-
return 1 - 2 * Math.pow(1 - d, n);
-
} else if (1 <= d) {
-
return 1;
}
-
- return exact ? exactK(d) : roundedK(d);
+ return exact ? exactK(d, n) : roundedK(d, n);
}
/**
- * Calculates the exact value of {@code P(D_n < d)} using method described
- * in [1] and {@link org.apache.commons.math3.fraction.BigFraction} (see
+ * Calculates the exact value of {@code P(D_n < d)} using the method described in [1] (reference
+ * in class javadoc above) and {@link org.apache.commons.math3.fraction.BigFraction} (see
* above).
*
* @param d statistic
- * @return the two-sided probability of {@code P(D_n < d)}
- * @throws MathArithmeticException if algorithm fails to convert {@code h}
- * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
- * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
+ * @param n sample size
+ * @return the two-sided probability of \(P(D_n < d)\)
+ * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\) for integer {@code k, m} and \(0 \le h < 1\).
*/
- private double exactK(double d) throws MathArithmeticException {
+ private double exactK(double d, int n)
+ throws MathArithmeticException {
final int k = (int) Math.ceil(n * d);
- final FieldMatrix<BigFraction> H = this.createH(d);
+ final FieldMatrix<BigFraction> H = this.createH(d, n);
final FieldMatrix<BigFraction> Hpower = H.power(n);
BigFraction pFrac = Hpower.getEntry(k - 1, k - 1);
@@ -205,50 +452,45 @@ public class KolmogorovSmirnovDistributi
}
/*
- * BigFraction.doubleValue converts numerator to double and the
- * denominator to double and divides afterwards. That gives NaN quite
- * easy. This does not (scale is the number of digits):
+ * BigFraction.doubleValue converts numerator to double and the denominator to double and
+ * divides afterwards. That gives NaN quite easy. This does not (scale is the number of
+ * digits):
*/
return pFrac.bigDecimalValue(20, BigDecimal.ROUND_HALF_UP).doubleValue();
}
/**
- * Calculates {@code P(D_n < d)} using method described in [1] and doubles
- * (see above).
+ * Calculates {@code P(D_n < d)} using method described in [1] and doubles (see above).
*
* @param d statistic
- * @return the two-sided probability of {@code P(D_n < d)}
- * @throws MathArithmeticException if algorithm fails to convert {@code h}
- * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing
- * {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
+ * @param n sample size
+ * @return the two-sided probability of \(P(D_n < d)\)
+ * @throws MathArithmeticException if algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\ for integer {@code k, m} and \(0 <= h < 1\).
*/
- private double roundedK(double d) throws MathArithmeticException {
+ private double roundedK(double d, int n)
+ throws MathArithmeticException {
final int k = (int) Math.ceil(n * d);
- final FieldMatrix<BigFraction> HBigFraction = this.createH(d);
+ final FieldMatrix<BigFraction> HBigFraction = this.createH(d, n);
final int m = HBigFraction.getRowDimension();
/*
- * Here the rounding part comes into play: use
- * RealMatrix instead of FieldMatrix<BigFraction>
+ * Here the rounding part comes into play: use RealMatrix instead of
+ * FieldMatrix<BigFraction>
*/
final RealMatrix H = new Array2DRowRealMatrix(m, m);
-
for (int i = 0; i < m; ++i) {
for (int j = 0; j < m; ++j) {
H.setEntry(i, j, HBigFraction.getEntry(i, j).doubleValue());
}
}
-
final RealMatrix Hpower = H.power(n);
-
double pFrac = Hpower.getEntry(k - 1, k - 1);
-
for (int i = 1; i <= n; ++i) {
pFrac *= (double) i / (double) n;
}
-
return pFrac;
}
@@ -256,37 +498,32 @@ public class KolmogorovSmirnovDistributi
* Creates {@code H} of size {@code m x m} as described in [1] (see above).
*
* @param d statistic
+ * @param n sample size
* @return H matrix
* @throws NumberIsTooLargeException if fractional part is greater than 1
- * @throws FractionConversionException if algorithm fails to convert
- * {@code h} to a {@link org.apache.commons.math3.fraction.BigFraction} in
- * expressing {@code d} as {@code (k - h) / m} for integer {@code k, m} and
- * {@code 0 <= h < 1}.
+ * @throws FractionConversionException if algorithm fails to convert {@code h} to a
+ * {@link org.apache.commons.math3.fraction.BigFraction} in expressing {@code d} as \((k
+ * - h) / m\) for integer {@code k, m} and \(0 <= h < 1\).
*/
- private FieldMatrix<BigFraction> createH(double d)
- throws NumberIsTooLargeException, FractionConversionException {
-
- int k = (int) Math.ceil(n * d);
-
- int m = 2 * k - 1;
- double hDouble = k - n * d;
+ private FieldMatrix<BigFraction> createH(double d, int n)
+ throws NumberIsTooLargeException, FractionConversionException {
+ final int k = (int) Math.ceil(n * d);
+ final int m = 2 * k - 1;
+ final double hDouble = k - n * d;
if (hDouble >= 1) {
throw new NumberIsTooLargeException(hDouble, 1.0, false);
}
-
BigFraction h = null;
-
try {
h = new BigFraction(hDouble, 1.0e-20, 10000);
- } catch (FractionConversionException e1) {
+ } catch (final FractionConversionException e1) {
try {
h = new BigFraction(hDouble, 1.0e-10, 10000);
- } catch (FractionConversionException e2) {
+ } catch (final FractionConversionException e2) {
h = new BigFraction(hDouble, 1.0e-5, 10000);
}
}
-
final BigFraction[][] Hdata = new BigFraction[m][m];
/*
@@ -303,8 +540,8 @@ public class KolmogorovSmirnovDistributi
}
/*
- * Setting up power-array to avoid calculating the same value twice:
- * hPowers[0] = h^1 ... hPowers[m-1] = h^m
+ * Setting up power-array to avoid calculating the same value twice: hPowers[0] = h^1 ...
+ * hPowers[m-1] = h^m
*/
final BigFraction[] hPowers = new BigFraction[m];
hPowers[0] = h;
@@ -321,23 +558,19 @@ public class KolmogorovSmirnovDistributi
}
/*
- * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix
- * should be (1 - 2*h^m + (2h - 1)^m )/m!" Since 0 <= h < 1, then if h >
- * 1/2 is sufficient to check:
+ * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix should be (1 - 2*h^m +
+ * (2h - 1)^m )/m!" Since 0 <= h < 1, then if h > 1/2 is sufficient to check:
*/
if (h.compareTo(BigFraction.ONE_HALF) == 1) {
Hdata[m - 1][0] = Hdata[m - 1][0].add(h.multiply(2).subtract(1).pow(m));
}
/*
- * Aside from the first column and last row, the (i, j)-th element is
- * 1/(i - j + 1)! if i - j + 1 >= 0, else 0. 1's and 0's are already
- * put, so only division with (i - j + 1)! is needed in the elements
- * that have 1's. There is no need to calculate (i - j + 1)! and then
- * divide - small steps avoid overflows.
- *
- * Note that i - j + 1 > 0 <=> i + 1 > j instead of j'ing all the way to
- * m. Also note that it is started at g = 2 because dividing by 1 isn't
+ * Aside from the first column and last row, the (i, j)-th element is 1/(i - j + 1)! if i -
+ * j + 1 >= 0, else 0. 1's and 0's are already put, so only division with (i - j + 1)! is
+ * needed in the elements that have 1's. There is no need to calculate (i - j + 1)! and then
+ * divide - small steps avoid overflows. Note that i - j + 1 > 0 <=> i + 1 > j instead of
+ * j'ing all the way to m. Also note that it is started at g = 2 because dividing by 1 isn't
* really necessary.
*/
for (int i = 0; i < m; ++i) {
@@ -349,7 +582,198 @@ public class KolmogorovSmirnovDistributi
}
}
}
-
return new Array2DRowFieldMatrix<BigFraction>(BigFractionField.getInstance(), Hdata);
}
+
+ /**
+ * Verifies that {@code array} has length at least 2.
+ *
+ * @param array array to test
+ * @throws NullArgumentException if array is null
+ * @throws InsufficientDataException if array is too short
+ */
+ private void checkArray(double[] array) {
+ if (array == null) {
+ throw new NullArgumentException(LocalizedFormats.NULL_NOT_ALLOWED);
+ }
+ if (array.length < 2) {
+ throw new InsufficientDataException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, array.length,
+ 2);
+ }
+ }
+
+ /**
+ * Computes \( 1 + 2 \sum_{i=1}^\infty (-1)^i e^{-2 i^2 t^2} \) stopping when successive partial
+ * sums are within {@code tolerance} of one another, or when {@code maxIterations} partial sums
+ * have been computed. If the sum does not converge before {@code maxIterations} iterations a
+ * {@link TooManyIterationsException} is thrown.
+ *
+ * @param t argument
+ * @param tolerance Cauchy criterion for partial sums
+ * @param maxIterations maximum number of partial sums to compute
+ * @return Kolmogorov sum evaluated at t
+ * @throws TooManyIterationsException if the series does not converge
+ */
+ public double ksSum(double t, double tolerance, int maxIterations) {
+ // TODO: for small t (say less than 1), the alternative expansion in part 3 of [1]
+ // from class javadoc should be used.
+ final double x = -2 * t * t;
+ int sign = -1;
+ long i = 1;
+ double partialSum = 0.5d;
+ double delta = 1;
+ while (delta > tolerance && i < maxIterations) {
+ delta = FastMath.exp(x * i * i);
+ partialSum += sign * delta;
+ sign *= -1;
+ i++;
+ }
+ if (i == maxIterations) {
+ throw new TooManyIterationsException(maxIterations);
+ }
+ return partialSum * 2;
+ }
+
+ /**
+ * Computes \(P(D_{n,m} > d)\) if {@code strict} is {@code true}; otherwise \(P(D_{n,m} \ge
+ * d)\), where \(D_{n,m}\) is the 2-sample Kolmogorov-Smirnov statistic. See
+ * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
+ * <p>
+ * The returned probability is exact, obtained by enumerating all partitions of {@code m + n}
+ * into {@code m} and {@code n} sets, computing \(D_{n,m}\) for each partition and counting the
+ * number of partitions that yield \(D_{n,m}\) values exceeding (resp. greater than or equal to)
+ * {@code d}.
+ * </p>
+ * <p>
+ * <strong>USAGE NOTE</strong>: Since this method enumerates all combinations in \({m+n} \choose
+ * {n}\), it is very slow if called for large {@code m, n}. For this reason,
+ * {@link #kolmogorovSmirnovTest(double[], double[])} uses this only for {@code m * n < }
+ * {@value #SMALL_SAMPLE_PRODUCT}.
+ * </p>
+ *
+ * @param d D-statistic value
+ * @param n first sample size
+ * @param m second sample size
+ * @param strict whether or not the probability to compute is expressed as a strict inequality
+ * @return probability that a randomly selected m-n partition of m + n generates \(D_{n,m}\)
+ * greater than (resp. greater than or equal to) {@code d}
+ */
+ public double exactP(double d, int n, int m, boolean strict) {
+ Iterator<int[]> combinationsIterator = CombinatoricsUtils.combinationsIterator(n + m, n);
+ long tail = 0;
+ final double[] nSet = new double[n];
+ final double[] mSet = new double[m];
+ while (combinationsIterator.hasNext()) {
+ // Generate an n-set
+ final int[] nSetI = combinationsIterator.next();
+ // Copy the n-set to nSet and its complement to mSet
+ int j = 0;
+ int k = 0;
+ for (int i = 0; i < n + m; i++) {
+ if (j < n && nSetI[j] == i) {
+ nSet[j++] = i;
+ } else {
+ mSet[k++] = i;
+ }
+ }
+ final double curD = kolmogorovSmirnovStatistic(nSet, mSet);
+ if (curD > d) {
+ tail++;
+ } else if (curD == d && !strict) {
+ tail++;
+ }
+ }
+ return (double) tail / (double) CombinatoricsUtils.binomialCoefficient(n + m, n);
+ }
+
+ /**
+ * Uses the Kolmogorov-Smirnov distribution to approximate \(P(D_{n,m} > d)\) where \(D_{n,m}\)
+ * is the 2-sample Kolmogorov-Smirnov statistic. See
+ * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
+ * <p>
+ * Specifically, what is returned is \(1 - k(d \sqrt{mn / (m + n)})\) where \(k(t) = 1 + 2
+ * \sum_{i=1}^\infty (-1)^i e^{-2 i^2 t^2}\). See {@link #ksSum(double, double, int)} for
+ * details on how convergence of the sum is determined. This implementation passes {@code ksSum}
+ * {@value #KS_SUM_CAUCHY_CRITERION} as {@code tolerance} and
+ * {@value #MAXIMUM_PARTIAL_SUM_COUNT} as {@code maxIterations}.
+ * </p>
+ *
+ * @param d D-statistic value
+ * @param n first sample size
+ * @param m second sample size
+ * @return approximate probability that a randomly selected m-n partition of m + n generates
+ * \(D_{n,m}\) greater than {@code d}
+ */
+ public double approximateP(double d, int n, int m) {
+ final double dm = m;
+ final double dn = n;
+ return 1 - ksSum(d * FastMath.sqrt((dm * dn) / (dm + dn)), KS_SUM_CAUCHY_CRITERION, MAXIMUM_PARTIAL_SUM_COUNT);
+ }
+
+ /**
+ * Uses Monte Carlo simulation to approximate \(P(D_{n,m} > d)\) where \(D_{n,m}\) is the
+ * 2-sample Kolmogorov-Smirnov statistic. See
+ * {@link #kolmogorovSmirnovStatistic(double[], double[])} for the definition of \(D_{n,m}\).
+ * <p>
+ * The simulation generates {@code iterations} random partitions of {@code m + n} into an
+ * {@code n} set and an {@code m} set, computing \(D_{n,m}\) for each partition and returning
+ * the proportion of values that are greater than {@code d}, or greater than or equal to
+ * {@code d} if {@code strict} is {@code false}.
+ * </p>
+ *
+ * @param d D-statistic value
+ * @param n first sample size
+ * @param m second sample size
+ * @param iterations number of random partitions to generate
+ * @param strict whether or not the probability to compute is expressed as a strict inequality
+ * @return proportion of randomly generated m-n partitions of m + n that result in \(D_{n,m}\)
+ * greater than (resp. greater than or equal to) {@code d}
+ */
+ public double monteCarloP(double d, int n, int m, boolean strict, int iterations) {
+ final int[] nPlusMSet = MathArrays.natural(m + n);
+ final double[] nSet = new double[n];
+ final double[] mSet = new double[m];
+ int tail = 0;
+ for (int i = 0; i < iterations; i++) {
+ copyPartition(nSet, mSet, nPlusMSet, n, m);
+ final double curD = kolmogorovSmirnovStatistic(nSet, mSet);
+ if (curD > d) {
+ tail++;
+ } else if (curD == d && !strict) {
+ tail++;
+ }
+ MathArrays.shuffle(nPlusMSet, rng);
+ Arrays.sort(nPlusMSet, 0, n);
+ }
+ return (double) tail / iterations;
+ }
+
+ /**
+ * Copies the first {@code n} elements of {@code nSetI} into {@code nSet} and its complement
+ * relative to {@code m + n} into {@code mSet}. For example, if {@code m = 3}, {@code n = 3} and
+ * {@code nSetI = [1,4,5,2,3,0]} then after this method returns, we will have
+ * {@code nSet = [1,4,5], mSet = [0,2,3]}.
+ * <p>
+ * <strong>Precondition:</strong> The first {@code n} elements of {@code nSetI} must be sorted
+ * in ascending order.
+ * </p>
+ *
+ * @param nSet array to fill with the first {@code n} elements of {@code nSetI}
+ * @param mSet array to fill with the {@code m} complementary elements of {@code nSet} relative
+ * to {@code m + n}
+ * @param nSetI array whose first {@code n} elements specify the members of {@code nSet}
+ * @param n number of elements in the first output array
+ * @param m number of elements in the second output array
+ */
+ private void copyPartition(double[] nSet, double[] mSet, int[] nSetI, int n, int m) {
+ int j = 0;
+ int k = 0;
+ for (int i = 0; i < n + m; i++) {
+ if (j < n && nSetI[j] == i) {
+ nSet[j++] = i;
+ } else {
+ mSet[k++] = i;
+ }
+ }
+ }
}
Added: commons/proper/math/trunk/src/test/R/KolmogorovSmirnovTestCases.R
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/R/KolmogorovSmirnovTestCases.R?rev=1572335&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/R/KolmogorovSmirnovTestCases.R (added)
+++ commons/proper/math/trunk/src/test/R/KolmogorovSmirnovTestCases.R Wed Feb 26 23:14:08 2014
@@ -0,0 +1,172 @@
+# Licensed to the Apache Software Foundation (ASF) under one or more
+# contributor license agreements. See the NOTICE file distributed with
+# this work for additional information regarding copyright ownership.
+# The ASF licenses this file to You 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.
+#
+#------------------------------------------------------------------------------
+# R source file to validate ChiSquare tests in
+# org.apache.commons.math.stat.inference.ChiSquareTestTest
+#
+# To run the test, install R, put this file and testFunctions
+# into the same directory, launch R from this directory and then enter
+# source("<name-of-this-file>")
+#
+#------------------------------------------------------------------------------
+tol <- 1E-14 # error tolerance for tests
+#------------------------------------------------------------------------------
+# Function definitions
+
+source("testFunctions") # utility test functions
+
+verifyOneSampleGaussian <- function(data, expectedP, expectedD, mean, sigma, exact, tol, desc) {
+ results <- ks.test(data, "pnorm", mean, sigma, exact = exact)
+ if (assertEquals(expectedP, results$p.value, tol, "p-value")) {
+ displayPadded(c(desc," p-value test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " p-value test"), FAILED, WIDTH)
+ }
+ if (assertEquals(expectedD, results$statistic, tol, "D statistic value")) {
+ displayPadded(c(desc," D statistic test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " D statistic test"), FAILED, WIDTH)
+ }
+}
+
+verifyOneSampleUniform <- function(data, expectedP, expectedD, min, max, exact, tol, desc) {
+ results <- ks.test(data, "punif", min, max, exact = exact)
+ if (assertEquals(expectedP, results$p.value, tol, "p-value")) {
+ displayPadded(c(desc," p-value test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " p-value test"), FAILED, WIDTH)
+ }
+ if (assertEquals(expectedD, results$statistic, tol, "D statistic value")) {
+ displayPadded(c(desc," D statistic test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " D statistic test"), FAILED, WIDTH)
+ }
+}
+
+verifyTwoSampleLargeSamples <- function(sample1, sample2, expectedP, expectedD, tol, desc) {
+ results <- ks.test(sample1, sample2)
+ if (assertEquals(expectedP, results$p.value, tol, "p-value")) {
+ displayPadded(c(desc," p-value test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " p-value test"), FAILED, WIDTH)
+ }
+ if (assertEquals(expectedD, results$statistic, tol, "D statistic value")) {
+ displayPadded(c(desc," D statistic test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " D statistic test"), FAILED, WIDTH)
+ }
+}
+
+verifyTwoSampleSmallSamplesExact <- function(sample1, sample2, expectedP, expectedD, tol, desc) {
+ results <- ks.test(sample1, sample2, exact = TRUE)
+ if (assertEquals(expectedP, results$p.value, tol, "p-value")) {
+ displayPadded(c(desc," p-value test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " p-value test"), FAILED, WIDTH)
+ }
+ if (assertEquals(expectedD, results$statistic, tol, "D statistic value")) {
+ displayPadded(c(desc," D statistic test"), SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(c(desc, " D statistic test"), FAILED, WIDTH)
+ }
+}
+
+
+cat("KolmogorovSmirnovTest test cases\n")
+
+gaussian <- c(0.26055895, -0.63665233, 1.51221323, 0.61246988, -0.03013003, -1.73025682,
+ -0.51435805, 0.70494168, 0.18242945, 0.94734336, -0.04286604, -0.37931719, -1.07026403, -2.05861425,
+ 0.11201862, 0.71400136, -0.52122185, -0.02478725, -1.86811649, -1.79907688, 0.15046279, 1.32390193,
+ 1.55889719, 1.83149171, -0.03948003, -0.98579207, -0.76790540, 0.89080682, 0.19532153, 0.40692841,
+ 0.15047336, -0.58546562, -0.39865469, 0.77604271, -0.65188221, -1.80368554, 0.65273365, -0.75283102,
+ -1.91022150, -0.07640869, -1.08681188, -0.89270600, 2.09017508, 0.43907981, 0.10744033, -0.70961218,
+ 1.15707300, 0.44560525, -2.04593349, 0.53816843, -0.08366640, 0.24652218, 1.80549401, -0.99220707,
+ -1.14589408, -0.27170290, -0.49696855, 0.00968353, -1.87113545, -1.91116529, 0.97151891, -0.73576115,
+ -0.59437029, 0.72148436, 0.01747695, -0.62601157, -1.00971538, -1.42691397, 1.03250131, -0.30672627,
+ -0.15353992, -1.19976069, -0.68364218, 0.37525652, -0.46592881, -0.52116168, -0.17162202, 1.04679215,
+ 0.25165971, -0.04125231, -0.23756244, -0.93389975, 0.75551407, 0.08347445, -0.27482228, -0.4717632,
+ -0.1867746, -0.1166976, 0.5763333, 0.1307952, 0.7630584, -0.3616248, 2.1383790,-0.7946630,
+ 0.0231885, 0.7919195, 1.6057144, -0.3802508, 0.1229078, 1.5252901, -0.8543149, 0.3025040)
+
+shortGaussian <- gaussian[1:50]
+
+gaussian2 <- c(2.88041498038308, -0.632349445671017, 0.402121295225571, 0.692626364613243, 1.30693446815426,
+ -0.714176317131286, -0.233169206599583, 1.09113298322107, -1.53149079994305, 1.23259966205809,
+ 1.01389927412503, 0.0143898711497477, -0.512813545447559, 2.79364360835469, 0.662008875538092,
+ 1.04861546834788, -0.321280099931466, 0.250296656278743, 1.75820367603736, -2.31433523590905,
+ -0.462694696086403, 0.187725700950191, -2.24410950019152, 2.83473751105445, 0.252460174391016,
+ 1.39051945380281, -1.56270144203134, 0.998522814471644, -1.50147469080896, 0.145307533554146,
+ 0.469089457043406, -0.0914780723809334, -0.123446939266548, -0.610513388160565, -3.71548343891957,
+ -0.329577317349478, -0.312973794075871, 2.02051909758923, 2.85214308266271, 0.0193222002327237,
+ -0.0322422268266562, 0.514736012106768, 0.231484953375887, -2.22468798953629, 1.42197716075595,
+ 2.69988043856357, 0.0443757119128293, 0.721536984407798, -0.0445688839903234, -0.294372724550705,
+ 0.234041580912698, -0.868973119365727, 1.3524893453845, -0.931054600134503, -0.263514296006792,
+ 0.540949457402918, -0.882544288773685, -0.34148675747989, 1.56664494810034, 2.19850536566584,
+ -0.667972122928022, -0.70889669526203, -0.00251758193079668, 2.39527162977682, -2.7559594317269,
+ -0.547393502656671, -2.62144031572617, 2.81504147017922, -1.02036850201042, -1.00713927602786,
+ -0.520197775122254, 1.00625480138649, 2.46756916531313, 1.64364743727799, 0.704545210648595,
+ -0.425885789416992, -1.78387854908546, -0.286783886710481, 0.404183648369076, -0.369324280845769,
+ -0.0391185138840443, 2.41257787857293, 2.49744281317859, -0.826964496939021, -0.792555379958975,
+ 1.81097685787403, -0.475014580016638, 1.23387615291805, 0.646615294802053, 1.88496377454523, 1.20390698380814,
+ -0.27812153371728, 2.50149494533101, 0.406964323253817, -1.72253451309982, 1.98432494184332, 2.2223658560333,
+ 0.393086362404685, -0.504073151377089, -0.0484610869883821)
+
+uniform <- c(0.7930305, 0.6424382, 0.8747699, 0.7156518, 0.1845909, 0.2022326,
+ 0.4877206, 0.8928752, 0.2293062, 0.4222006, 0.1610459, 0.2830535, 0.9946345, 0.7329499,
+ 0.26411126, 0.87958133, 0.29827437, 0.39185988, 0.38351185, 0.36359611, 0.48646472, 0.05577866,
+ 0.56152250, 0.52672013, 0.13171783, 0.95864085, 0.03060207, 0.33514887, 0.72508148, 0.38901437,
+ 0.9978665, 0.5981300, 0.1065388, 0.7036991, 0.1071584, 0.4423963, 0.1107071, 0.6437221,
+ 0.58523872, 0.05044634, 0.65999539, 0.37367260, 0.73270024, 0.47473755, 0.74661163, 0.50765549,
+ 0.05377347, 0.40998009, 0.55235182, 0.21361998, 0.63117971, 0.18109222, 0.89153510, 0.23203248,
+ 0.6177106, 0.6856418, 0.2158557, 0.9870501, 0.2036914, 0.2100311, 0.9065020, 0.7459159,
+ 0.56631790, 0.06753629, 0.39684629, 0.52504615, 0.14199103, 0.78551120, 0.90503321, 0.80452362,
+ 0.9960115, 0.8172592, 0.5831134, 0.8794187, 0.2021501, 0.2923505, 0.9561824, 0.8792248,
+ 0.85201008, 0.02945562, 0.26200374, 0.11382818, 0.17238856, 0.36449473, 0.69688273, 0.96216330,
+ 0.4859432,0.4503438, 0.1917656, 0.8357845, 0.9957812, 0.4633570, 0.8654599, 0.4597996,
+ 0.68190289, 0.58887855, 0.09359396, 0.98081979, 0.73659533, 0.89344777, 0.18903099, 0.97660425)
+
+smallSample1 <- c(6, 7, 9, 13, 19, 21, 22, 23, 24)
+smallSample2 <- c(10, 11, 12, 16, 20, 27, 28, 32, 44, 54)
+
+shortUniform <- uniform[1:20]
+
+verifyOneSampleGaussian(gaussian, 0.3172069207622391, 0.0932947561266756, 0, 1,
+TRUE, tol, "One sample gaussian - gaussian values")
+
+verifyOneSampleGaussian(shortGaussian, 0.683736463728347, 0.09820779969463278, 0, 1,
+TRUE, tol, "One sample gaussian - gaussian values - small sample")
+
+verifyOneSampleGaussian(uniform, 8.881784197001252E-16, 0.5117493931609258, 0, 1,
+TRUE, tol, "One sample gaussian - uniform values")
+
+verifyOneSampleUniform(uniform, 8.881784197001252E-16, 0.5400666982352942, -0.5, 0.5,
+TRUE, tol, "One sample uniform - uniform values")
+
+verifyOneSampleUniform(shortUniform, 4.117594598618268E-9, 0.6610459, -0.5, 0.5,
+TRUE, tol, "One sample uniform - uniform values - small sample")
+
+verifyOneSampleUniform(gaussian, 4.9405812774239166E-11, 0.3401058049019608, -0.5, 0.5,
+TRUE, tol, "One sample uniform - unit normal values")
+
+verifyTwoSampleLargeSamples(gaussian, gaussian2, 0.0319983962391632, 0.202352941176471, tol,
+"Two sample N(0, 1) vs N(0, 1.6)")
+
+verifyTwoSampleSmallSamplesExact(smallSample1, smallSample2, 0.105577085453247, .5, tol,
+"Two sample small samples exact")
+
+displayDashes(WIDTH)
+
+
Modified: commons/proper/math/trunk/src/test/R/testAll
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/R/testAll?rev=1572335&r1=1572334&r2=1572335&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/R/testAll (original)
+++ commons/proper/math/trunk/src/test/R/testAll Wed Feb 26 23:14:08 2014
@@ -48,6 +48,7 @@ source("regressionTestCases")
# inference
source("chiSquareTestCases")
source("anovaTestCases")
+source("KolmogorovSmirnovTestCases")
# descriptive
source("descriptiveTestCases")
Copied: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java (from r1488409, commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistributionTest.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java?p2=commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java&p1=commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistributionTest.java&r1=1488409&r2=1572335&rev=1572335&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistributionTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java Wed Feb 26 23:14:08 2014
@@ -15,104 +15,368 @@
* limitations under the License.
*/
-package org.apache.commons.math3.distribution;
+package org.apache.commons.math3.stat.inference;
+import org.apache.commons.math3.distribution.NormalDistribution;
+import org.apache.commons.math3.distribution.UniformRealDistribution;
+import org.apache.commons.math3.random.Well19937c;
+import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
- * Test cases for {@link KolmogorovSmirnovDistribution}.
- *
+ * Test cases for {@link KolmogorovSmirnovTest}.
+ *
* @version $Id$
+ * @since 3.3
*/
-public class KolmogorovSmirnovDistributionTest {
-
+public class KolmogorovSmirnovTestTest {
+
private static final double TOLERANCE = 10e-10;
+ // Random N(0,1) values generated using R rnorm
+ private final double[] gaussian = {
+ 0.26055895, -0.63665233, 1.51221323, 0.61246988, -0.03013003, -1.73025682, -0.51435805, 0.70494168, 0.18242945,
+ 0.94734336, -0.04286604, -0.37931719, -1.07026403, -2.05861425, 0.11201862, 0.71400136, -0.52122185,
+ -0.02478725, -1.86811649, -1.79907688, 0.15046279, 1.32390193, 1.55889719, 1.83149171, -0.03948003,
+ -0.98579207, -0.76790540, 0.89080682, 0.19532153, 0.40692841, 0.15047336, -0.58546562, -0.39865469, 0.77604271,
+ -0.65188221, -1.80368554, 0.65273365, -0.75283102, -1.91022150, -0.07640869, -1.08681188, -0.89270600,
+ 2.09017508, 0.43907981, 0.10744033, -0.70961218, 1.15707300, 0.44560525, -2.04593349, 0.53816843, -0.08366640,
+ 0.24652218, 1.80549401, -0.99220707, -1.14589408, -0.27170290, -0.49696855, 0.00968353, -1.87113545,
+ -1.91116529, 0.97151891, -0.73576115, -0.59437029, 0.72148436, 0.01747695, -0.62601157, -1.00971538,
+ -1.42691397, 1.03250131, -0.30672627, -0.15353992, -1.19976069, -0.68364218, 0.37525652, -0.46592881,
+ -0.52116168, -0.17162202, 1.04679215, 0.25165971, -0.04125231, -0.23756244, -0.93389975, 0.75551407,
+ 0.08347445, -0.27482228, -0.4717632, -0.1867746, -0.1166976, 0.5763333, 0.1307952, 0.7630584, -0.3616248,
+ 2.1383790, -0.7946630, 0.0231885, 0.7919195, 1.6057144, -0.3802508, 0.1229078, 1.5252901, -0.8543149, 0.3025040
+ };
+
+ // Random N(0, 1.6) values generated using R rnorm
+ private final double[] gaussian2 = {
+ 2.88041498038308, -0.632349445671017, 0.402121295225571, 0.692626364613243, 1.30693446815426,
+ -0.714176317131286, -0.233169206599583, 1.09113298322107, -1.53149079994305, 1.23259966205809,
+ 1.01389927412503, 0.0143898711497477, -0.512813545447559, 2.79364360835469, 0.662008875538092,
+ 1.04861546834788, -0.321280099931466, 0.250296656278743, 1.75820367603736, -2.31433523590905,
+ -0.462694696086403, 0.187725700950191, -2.24410950019152, 2.83473751105445, 0.252460174391016,
+ 1.39051945380281, -1.56270144203134, 0.998522814471644, -1.50147469080896, 0.145307533554146,
+ 0.469089457043406, -0.0914780723809334, -0.123446939266548, -0.610513388160565, -3.71548343891957,
+ -0.329577317349478, -0.312973794075871, 2.02051909758923, 2.85214308266271, 0.0193222002327237,
+ -0.0322422268266562, 0.514736012106768, 0.231484953375887, -2.22468798953629, 1.42197716075595,
+ 2.69988043856357, 0.0443757119128293, 0.721536984407798, -0.0445688839903234, -0.294372724550705,
+ 0.234041580912698, -0.868973119365727, 1.3524893453845, -0.931054600134503, -0.263514296006792,
+ 0.540949457402918, -0.882544288773685, -0.34148675747989, 1.56664494810034, 2.19850536566584,
+ -0.667972122928022, -0.70889669526203, -0.00251758193079668, 2.39527162977682, -2.7559594317269,
+ -0.547393502656671, -2.62144031572617, 2.81504147017922, -1.02036850201042, -1.00713927602786,
+ -0.520197775122254, 1.00625480138649, 2.46756916531313, 1.64364743727799, 0.704545210648595,
+ -0.425885789416992, -1.78387854908546, -0.286783886710481, 0.404183648369076, -0.369324280845769,
+ -0.0391185138840443, 2.41257787857293, 2.49744281317859, -0.826964496939021, -0.792555379958975,
+ 1.81097685787403, -0.475014580016638, 1.23387615291805, 0.646615294802053, 1.88496377454523, 1.20390698380814,
+ -0.27812153371728, 2.50149494533101, 0.406964323253817, -1.72253451309982, 1.98432494184332, 2.2223658560333,
+ 0.393086362404685, -0.504073151377089, -0.0484610869883821
+ };
+
+ // Random uniform (0, 1) generated using R runif
+ private final double[] uniform = {
+ 0.7930305, 0.6424382, 0.8747699, 0.7156518, 0.1845909, 0.2022326, 0.4877206, 0.8928752, 0.2293062, 0.4222006,
+ 0.1610459, 0.2830535, 0.9946345, 0.7329499, 0.26411126, 0.87958133, 0.29827437, 0.39185988, 0.38351185,
+ 0.36359611, 0.48646472, 0.05577866, 0.56152250, 0.52672013, 0.13171783, 0.95864085, 0.03060207, 0.33514887,
+ 0.72508148, 0.38901437, 0.9978665, 0.5981300, 0.1065388, 0.7036991, 0.1071584, 0.4423963, 0.1107071, 0.6437221,
+ 0.58523872, 0.05044634, 0.65999539, 0.37367260, 0.73270024, 0.47473755, 0.74661163, 0.50765549, 0.05377347,
+ 0.40998009, 0.55235182, 0.21361998, 0.63117971, 0.18109222, 0.89153510, 0.23203248, 0.6177106, 0.6856418,
+ 0.2158557, 0.9870501, 0.2036914, 0.2100311, 0.9065020, 0.7459159, 0.56631790, 0.06753629, 0.39684629,
+ 0.52504615, 0.14199103, 0.78551120, 0.90503321, 0.80452362, 0.9960115, 0.8172592, 0.5831134, 0.8794187,
+ 0.2021501, 0.2923505, 0.9561824, 0.8792248, 0.85201008, 0.02945562, 0.26200374, 0.11382818, 0.17238856,
+ 0.36449473, 0.69688273, 0.96216330, 0.4859432, 0.4503438, 0.1917656, 0.8357845, 0.9957812, 0.4633570,
+ 0.8654599, 0.4597996, 0.68190289, 0.58887855, 0.09359396, 0.98081979, 0.73659533, 0.89344777, 0.18903099,
+ 0.97660425
+ };
+
@Test
public void testCumulativeDensityFunction() {
-
- KolmogorovSmirnovDistribution dist;
-
- /* The code below is generated using the R-script located in
- * /src/test/R/KolmogorovSmirnovDistributionTestCases.R
- */
-
- /* R version 2.11.1 (2010-05-31) */
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+
+ /*
+ * The code below is generated using the R-script located in
+ * /src/test/R/KolmogorovSmirnovDistributionTestCases.R R version 2.11.1 (2010-05-31)
+ */
- /* formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(200), PACKAGE = "stats")$p, 40) gives
- * 4.907829957616471622388047046469198862537e-86
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(200), PACKAGE =
+ * "stats")$p, 40) gives 4.907829957616471622388047046469198862537e-86
*/
- dist = new KolmogorovSmirnovDistribution(200);
- Assert.assertEquals(4.907829957616471622388047046469198862537e-86, dist.cdf(0.005, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(200), PACKAGE = "stats")$p, 40) gives
- * 5.151982014280041957199687829849210629618e-06
+ Assert.assertEquals(4.907829957616471622388047046469198862537e-86, test.cdf(0.005, 200, false), TOLERANCE);
+
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(200), PACKAGE =
+ * "stats")$p, 40) gives 5.151982014280041957199687829849210629618e-06
*/
- dist = new KolmogorovSmirnovDistribution(200);
- Assert.assertEquals(5.151982014280041957199687829849210629618e-06, dist.cdf(0.02, false), TOLERANCE);
+ Assert.assertEquals(5.151982014280041957199687829849210629618e-06, test.cdf(0.02, 200, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(200), PACKAGE = "stats")$p, 40) gives
- * 0.01291614648162886340443389343590752105229
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(200), PACKAGE =
+ * "stats")$p, 40) gives 0.01291614648162886340443389343590752105229
*/
- dist = new KolmogorovSmirnovDistribution(200);
- Assert.assertEquals(0.01291614648162886340443389343590752105229, dist.cdf(0.031111, false), TOLERANCE);
+ Assert.assertEquals(0.01291614648162886340443389343590752105229, test.cdf(0.031111, 200, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(200), PACKAGE = "stats")$p, 40) gives
- * 0.1067137011362679355208626930107129737735
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(200), PACKAGE =
+ * "stats")$p, 40) gives 0.1067137011362679355208626930107129737735
*/
- dist = new KolmogorovSmirnovDistribution(200);
- Assert.assertEquals(0.1067137011362679355208626930107129737735, dist.cdf(0.04, false), TOLERANCE);
+ Assert.assertEquals(0.1067137011362679355208626930107129737735, test.cdf(0.04, 200, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(341), PACKAGE = "stats")$p, 40) gives
- * 1.914734701559404553985102395145063418825e-53
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(341), PACKAGE =
+ * "stats")$p, 40) gives 1.914734701559404553985102395145063418825e-53
*/
- dist = new KolmogorovSmirnovDistribution(341);
- Assert.assertEquals(1.914734701559404553985102395145063418825e-53, dist.cdf(0.005, false), TOLERANCE);
+ Assert.assertEquals(1.914734701559404553985102395145063418825e-53, test.cdf(0.005, 341, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(341), PACKAGE = "stats")$p, 40) gives
- * 0.001171328985781981343872182321774744195864
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(341), PACKAGE =
+ * "stats")$p, 40) gives 0.001171328985781981343872182321774744195864
*/
- dist = new KolmogorovSmirnovDistribution(341);
- Assert.assertEquals(0.001171328985781981343872182321774744195864, dist.cdf(0.02, false), TOLERANCE);
+ Assert.assertEquals(0.001171328985781981343872182321774744195864, test.cdf(0.02, 341, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(341), PACKAGE = "stats")$p, 40) gives
- * 0.1142955196267499418105728636874118819833
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(341), PACKAGE =
+ * "stats")$p, 40) gives 0.1142955196267499418105728636874118819833
*/
- dist = new KolmogorovSmirnovDistribution(341);
- Assert.assertEquals(0.1142955196267499418105728636874118819833, dist.cdf(0.031111, false), TOLERANCE);
+ Assert.assertEquals(0.1142955196267499418105728636874118819833, test.cdf(0.031111, 341, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(341), PACKAGE = "stats")$p, 40) gives
- * 0.3685529520496805266915885113121476024389
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(341), PACKAGE =
+ * "stats")$p, 40) gives 0.3685529520496805266915885113121476024389
*/
- dist = new KolmogorovSmirnovDistribution(341);
- Assert.assertEquals(0.3685529520496805266915885113121476024389, dist.cdf(0.04, false), TOLERANCE);
+ Assert.assertEquals(0.3685529520496805266915885113121476024389, test.cdf(0.04, 341, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(389), PACKAGE = "stats")$p, 40) gives
- * 1.810657144595055888918455512707637574637e-47
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.005), n = as.integer(389), PACKAGE =
+ * "stats")$p, 40) gives 1.810657144595055888918455512707637574637e-47
*/
- dist = new KolmogorovSmirnovDistribution(389);
- Assert.assertEquals(1.810657144595055888918455512707637574637e-47, dist.cdf(0.005, false), TOLERANCE);
+ Assert.assertEquals(1.810657144595055888918455512707637574637e-47, test.cdf(0.005, 389, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(389), PACKAGE = "stats")$p, 40) gives
- * 0.003068542559702356568168690742481885536108
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.02), n = as.integer(389), PACKAGE =
+ * "stats")$p, 40) gives 0.003068542559702356568168690742481885536108
*/
- dist = new KolmogorovSmirnovDistribution(389);
- Assert.assertEquals(0.003068542559702356568168690742481885536108, dist.cdf(0.02, false), TOLERANCE);
+ Assert.assertEquals(0.003068542559702356568168690742481885536108, test.cdf(0.02, 389, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(389), PACKAGE = "stats")$p, 40) gives
- * 0.1658291700122746237244797384846606291831
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.031111), n = as.integer(389), PACKAGE =
+ * "stats")$p, 40) gives 0.1658291700122746237244797384846606291831
*/
- dist = new KolmogorovSmirnovDistribution(389);
- Assert.assertEquals(0.1658291700122746237244797384846606291831, dist.cdf(0.031111, false), TOLERANCE);
+ Assert.assertEquals(0.1658291700122746237244797384846606291831, test.cdf(0.031111, 389, false), TOLERANCE);
- /* formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(389), PACKAGE = "stats")$p, 40) gives
- * 0.4513143712128902529379104180407011881471
+ /*
+ * formatC(.C("pkolmogorov2x", p = as.double(0.04), n = as.integer(389), PACKAGE =
+ * "stats")$p, 40) gives 0.4513143712128902529379104180407011881471
*/
- dist = new KolmogorovSmirnovDistribution(389);
- Assert.assertEquals(0.4513143712128902529379104180407011881471, dist.cdf(0.04, false), TOLERANCE);
+ Assert.assertEquals(0.4513143712128902529379104180407011881471, test.cdf(0.04, 389, false), TOLERANCE);
+ }
+
+ /** Unit normal distribution, unit normal data */
+ @Test
+ public void testOneSampleGaussianGaussian() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final NormalDistribution unitNormal = new NormalDistribution(0d, 1d);
+ // Uncomment to run exact test - takes about a minute. Same value is used in R tests and for
+ // approx.
+ // Assert.assertEquals(0.3172069207622391, test.kolmogorovSmirnovTest(unitNormal, gaussian,
+ // true), TOLERANCE);
+ Assert.assertEquals(0.3172069207622391, test.kolmogorovSmirnovTest(unitNormal, gaussian, false), TOLERANCE);
+ Assert.assertFalse(test.kolmogorovSmirnovTest(unitNormal, gaussian, 0.05));
+ Assert.assertEquals(0.0932947561266756, test.kolmogorovSmirnovStatistic(unitNormal, gaussian), TOLERANCE);
+ }
+ /** Unit normal distribution, unit normal data, small dataset */
+ @Test
+ public void testOneSampleGaussianGaussianSmallSample() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final NormalDistribution unitNormal = new NormalDistribution(0d, 1d);
+ final double[] shortGaussian = new double[50];
+ System.arraycopy(gaussian, 0, shortGaussian, 0, 50);
+ Assert.assertEquals(0.683736463728347, test.kolmogorovSmirnovTest(unitNormal, shortGaussian, false), TOLERANCE);
+ Assert.assertFalse(test.kolmogorovSmirnovTest(unitNormal, gaussian, 0.05));
+ Assert.assertEquals(0.09820779969463278, test.kolmogorovSmirnovStatistic(unitNormal, shortGaussian), TOLERANCE);
+ }
+
+ /** Unit normal distribution, uniform data */
+ @Test
+ public void testOneSampleGaussianUniform() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final NormalDistribution unitNormal = new NormalDistribution(0d, 1d);
+ // Uncomment to run exact test - takes a long time. Same value is used in R tests and for
+ // approx.
+ // Assert.assertEquals(0.3172069207622391, test.kolmogorovSmirnovTest(unitNormal, uniform,
+ // true), TOLERANCE);
+ Assert.assertEquals(8.881784197001252E-16, test.kolmogorovSmirnovTest(unitNormal, uniform, false), TOLERANCE);
+ Assert.assertFalse(test.kolmogorovSmirnovTest(unitNormal, gaussian, 0.05));
+ Assert.assertEquals(0.5117493931609258, test.kolmogorovSmirnovStatistic(unitNormal, uniform), TOLERANCE);
+ }
+
+ /** Uniform distribution, uniform data */
+ // @Test - takes about 6 seconds, uncomment for
+ public void testOneSampleUniformUniform() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final UniformRealDistribution unif = new UniformRealDistribution(-0.5, 0.5);
+ Assert.assertEquals(8.881784197001252E-16, test.kolmogorovSmirnovTest(unif, uniform, false), TOLERANCE);
+ Assert.assertTrue(test.kolmogorovSmirnovTest(unif, uniform, 0.05));
+ Assert.assertEquals(0.5400666982352942, test.kolmogorovSmirnovStatistic(unif, uniform), TOLERANCE);
}
+ /** Uniform distribution, uniform data, small dataset */
+ @Test
+ public void testOneSampleUniformUniformSmallSample() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final UniformRealDistribution unif = new UniformRealDistribution(-0.5, 0.5);
+ final double[] shortUniform = new double[20];
+ System.arraycopy(uniform, 0, shortUniform, 0, 20);
+ Assert.assertEquals(4.117594598618268E-9, test.kolmogorovSmirnovTest(unif, shortUniform, false), TOLERANCE);
+ Assert.assertTrue(test.kolmogorovSmirnovTest(unif, shortUniform, 0.05));
+ Assert.assertEquals(0.6610459, test.kolmogorovSmirnovStatistic(unif, shortUniform), TOLERANCE);
+ }
+
+ /** Uniform distribution, unit normal dataset */
+ @Test
+ public void testOneSampleUniformGaussian() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final UniformRealDistribution unif = new UniformRealDistribution(-0.5, 0.5);
+ // Value was obtained via exact test, validated against R. Running exact test takes a long
+ // time.
+ Assert.assertEquals(4.9405812774239166E-11, test.kolmogorovSmirnovTest(unif, gaussian, false), TOLERANCE);
+ Assert.assertTrue(test.kolmogorovSmirnovTest(unif, gaussian, 0.05));
+ Assert.assertEquals(0.3401058049019608, test.kolmogorovSmirnovStatistic(unif, gaussian), TOLERANCE);
+ }
+
+ /** Small samples - exact p-value, checked against R */
+ @Test
+ public void testTwoSampleSmallSampleExact() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final double[] smallSample1 = {
+ 6, 7, 9, 13, 19, 21, 22, 23, 24
+ };
+ final double[] smallSample2 = {
+ 10, 11, 12, 16, 20, 27, 28, 32, 44, 54
+ };
+ // Reference values from R, version 2.15.3 - R uses non-strict inequality in null hypothesis
+ Assert
+ .assertEquals(0.105577085453247, test.kolmogorovSmirnovTest(smallSample1, smallSample2, false), TOLERANCE);
+ Assert.assertEquals(0.5, test.kolmogorovSmirnovStatistic(smallSample1, smallSample2), TOLERANCE);
+ }
+
+ /**
+ * Checks exact p-value computations using critical values from Table 9 in V.K Rohatgi, An
+ * Introduction to Probability and Mathematical Statistics, Wiley, 1976, ISBN 0-471-73135-8.
+ */
+ @Test
+ public void testTwoSampleExactP() {
+ checkExactTable(4, 6, 5d / 6d, 0.01d);
+ checkExactTable(4, 7, 17d / 28d, 0.2d);
+ checkExactTable(6, 7, 29d / 42d, 0.05d);
+ checkExactTable(4, 10, 7d / 10d, 0.05d);
+ checkExactTable(5, 15, 11d / 15d, 0.02d);
+ checkExactTable(9, 10, 31d / 45d, 0.01d);
+ checkExactTable(7, 10, 43d / 70d, 0.05d);
+ }
+
+ @Test
+ public void testTwoSampleApproximateCritialValues() {
+ final double tol = .01;
+ final double[] alpha = {
+ 0.10, 0.05, 0.025, 0.01, 0.005, 0.001
+ };
+ // From Wikipedia KS article - TODO: get (and test) more precise values
+ final double[] c = {
+ 1.22, 1.36, 1.48, 1.63, 1.73, 1.95
+ };
+ final int k[] = {
+ 60, 100, 500
+ };
+ double n;
+ double m;
+ for (int i = 0; i < k.length; i++) {
+ for (int j = 0; j < i; j++) {
+ n = k[i];
+ m = k[j];
+ for (int l = 0; l < alpha.length; l++) {
+ final double dCrit = c[l] * FastMath.sqrt((n + m) / (n * m));
+ checkApproximateTable(k[i], k[j], dCrit, alpha[l], tol);
+ }
+ }
+ }
+ }
+
+ /** Verifies large sample approximate p values against R */
+ @Test
+ public void testTwoSampleApproximateP() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ // Reference values from R, version 2.15.3
+ Assert.assertEquals(0.0319983962391632, test.kolmogorovSmirnovTest(gaussian, gaussian2), TOLERANCE);
+ Assert.assertEquals(0.202352941176471, test.kolmogorovSmirnovStatistic(gaussian, gaussian2), TOLERANCE);
+ }
+
+ /**
+ * Verifies that Monte Carlo simulation gives results close to exact p values. This test is a
+ * little long-running (more than two minutes on a fast machine), so is disabled by default.
+ */
+ // @Test
+ public void testTwoSampleMonteCarlo() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest(new Well19937c(1000));
+ final int sampleSize = 14;
+ final double tol = .001;
+ final double[] shortUniform = new double[sampleSize];
+ System.arraycopy(uniform, 0, shortUniform, 0, sampleSize);
+ final double[] shortGaussian = new double[sampleSize];
+ final double[] shortGaussian2 = new double[sampleSize];
+ System.arraycopy(gaussian, 0, shortGaussian, 0, sampleSize);
+ System.arraycopy(gaussian, 10, shortGaussian2, 0, sampleSize);
+ final double[] d = {
+ test.kolmogorovSmirnovStatistic(shortGaussian, shortUniform),
+ test.kolmogorovSmirnovStatistic(shortGaussian2, shortGaussian)
+ };
+ for (double dv : d) {
+ double exactPStrict = test.exactP(dv, sampleSize, sampleSize, true);
+ double exactPNonStrict = test.exactP(dv, sampleSize, sampleSize, false);
+ double montePStrict = test.monteCarloP(dv, sampleSize, sampleSize, true,
+ KolmogorovSmirnovTest.MONTE_CARLO_ITERATIONS);
+ double montePNonStrict = test.monteCarloP(dv, sampleSize, sampleSize, false,
+ KolmogorovSmirnovTest.MONTE_CARLO_ITERATIONS);
+ Assert.assertEquals(exactPStrict, montePStrict, tol);
+ Assert.assertEquals(exactPNonStrict, montePNonStrict, tol);
+ }
+ }
+
+ /**
+ * Verifies the inequality exactP(criticalValue, n, m, true) < alpha < exactP(criticalValue, n,
+ * m, false).
+ *
+ * Note that the validity of this check depends on the fact that alpha lies strictly between two
+ * attained values of the distribution and that criticalValue is one of the attained values. The
+ * critical value table (reference below) uses attained values. This test therefore also
+ * verifies that criticalValue is attained.
+ *
+ * @param n first sample size
+ * @param m second sample size
+ * @param criticalValue critical value
+ * @param alpha significance level
+ */
+ private void checkExactTable(int n, int m, double criticalValue, double alpha) {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ Assert.assertTrue(test.exactP(criticalValue, n, m, true) < alpha);
+ Assert.assertTrue(test.exactP(criticalValue, n, m, false) > alpha);
+ }
+
+ /**
+ * Verifies that approximateP(criticalValue, n, m) is within epsilon of alpha.
+ *
+ * @param n first sample size
+ * @param m second sample size
+ * @param criticalValue critical value (from table)
+ * @param alpha significance level
+ * @param epsilon tolerance
+ */
+ private void checkApproximateTable(int n, int m, double criticalValue, double alpha, double epsilon) {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ Assert.assertEquals(alpha, test.approximateP(criticalValue, n, m), epsilon);
+ }
}