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 2015/12/31 23:02:20 UTC
[math] Improved performance and accuracy of 2-sample KS tests. JIRA:
MATH-1310.
Repository: commons-math
Updated Branches:
refs/heads/MATH_3_X 56a628776 -> 2983ff81b
Improved performance and accuracy of 2-sample KS tests. JIRA: MATH-1310.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/2983ff81
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/2983ff81
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/2983ff81
Branch: refs/heads/MATH_3_X
Commit: 2983ff81b17638ea041d73676049d002b9096781
Parents: 56a6287
Author: Phil Steitz <ph...@gmail.com>
Authored: Thu Dec 31 15:02:11 2015 -0700
Committer: Phil Steitz <ph...@gmail.com>
Committed: Thu Dec 31 15:02:11 2015 -0700
----------------------------------------------------------------------
src/changes/changes.xml | 3 +
.../stat/inference/KolmogorovSmirnovTest.java | 161 +++++++++----------
src/test/R/KolmogorovSmirnovTestCases.R | 14 +-
.../inference/KolmogorovSmirnovTestTest.java | 44 ++++-
4 files changed, 134 insertions(+), 88 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/2983ff81/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 2f94e04..8b7a409 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -51,6 +51,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</properties>
<body>
<release version="3.6" date="XXXX-XX-XX" description="">
+ <action dev="psteitz" type="fix" issue="MATH-1310">
+ Improved performance and accuracy of 2-sample KolmogorovSmirnov tests.
+ </action>
<action dev="luc" type="fix" issue="MATH-1297">
Detect start failures with multi-step ODE integrators.
</action>
http://git-wip-us.apache.org/repos/asf/commons-math/blob/2983ff81/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java b/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java
index b9f2ec0..6b70e9b 100644
--- a/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java
+++ b/src/main/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.java
@@ -20,7 +20,6 @@ package org.apache.commons.math3.stat.inference;
import java.math.BigDecimal;
import java.util.Arrays;
import java.util.HashSet;
-import java.util.Iterator;
import org.apache.commons.math3.distribution.EnumeratedRealDistribution;
import org.apache.commons.math3.distribution.RealDistribution;
@@ -69,14 +68,9 @@ import org.apache.commons.math3.util.MathUtils;
* 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>For small samples (where the product of the sample sizes is less than
+ * {@value #LARGE_SAMPLE_PRODUCT}), the method presented in [4] is used to compute the
+ * exact p-value for the 2-sample test.</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>
@@ -97,8 +91,6 @@ import org.apache.commons.math3.util.MathUtils;
* 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
@@ -115,6 +107,8 @@ import org.apache.commons.math3.util.MathUtils;
* <li>[3] Jasjeet S. Sekhon. 2011. <a href="http://www.jstatsoft.org/article/view/v042i07">
* Multivariate and Propensity Score Matching Software with Automated Balance Optimization:
* The Matching package for R</a> Journal of Statistical Software, 42(7): 1-52.</li>
+ * <li>[4] Wilcox, Rand. 2012. Introduction to Robust Estimation and Hypothesis Testing,
+ * Chapter 5, 3rd Ed. Academic Press.</li>
* </ul>
* <br/>
* Note that [1] contains an error in computing h, refer to <a
@@ -136,16 +130,19 @@ public class KolmogorovSmirnovTest {
/** Convergence criterion for the sums in #pelzGood(double, double, int)} */
protected static final double PG_SUM_RELATIVE_ERROR = 1.0e-10;
- /** When product of sample sizes is less than this value, 2-sample K-S test is exact */
+ /** No longer used. */
+ @Deprecated
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.
+ * distribution to compute the p-value.
*/
protected static final int LARGE_SAMPLE_PRODUCT = 10000;
- /** Default number of iterations used by {@link #monteCarloP(double, int, int, boolean, int)} */
+ /** Default number of iterations used by {@link #monteCarloP(double, int, int, boolean, int)}.
+ * Deprecated as of version 3.6, as this method is no longer needed. */
+ @Deprecated
protected static final int MONTE_CARLO_ITERATIONS = 1000000;
/** Random data generator used by {@link #monteCarloP(double, int, int, boolean, int)} */
@@ -160,9 +157,12 @@ public class KolmogorovSmirnovTest {
/**
* Construct a KolmogorovSmirnovTest with the provided random data generator.
+ * The #monteCarloP(double, int, int, boolean, int) that uses the generator supplied to this
+ * constructor is deprecated as of version 3.6.
*
* @param rng random data generator used by {@link #monteCarloP(double, int, int, boolean, int)}
*/
+ @Deprecated
public KolmogorovSmirnovTest(RandomGenerator rng) {
this.rng = rng;
}
@@ -226,18 +226,9 @@ public class KolmogorovSmirnovTest {
* {@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. In the very
- * small sample case, if there are ties in the data, the actual sample values (including ties)
- * are used in generating the partitions (which are basically multi-set partitions in this
- * case).</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>For small samples (where the product of the sample sizes is less than
+ * {@value #LARGE_SAMPLE_PRODUCT}), the exact p-value is computed using the method presented
+ * in [4], implemented in {@link #exactP(double, int, int, boolean)}. </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>
@@ -274,11 +265,8 @@ public class KolmogorovSmirnovTest {
xa = x;
ya = y;
}
- if (lengthProduct < SMALL_SAMPLE_PRODUCT) {
- return integralExactP(integralKolmogorovSmirnovStatistic(xa, ya) + (strict?1l:0l), x.length, y.length);
- }
if (lengthProduct < LARGE_SAMPLE_PRODUCT) {
- return integralMonteCarloP(integralKolmogorovSmirnovStatistic(xa, ya) + (strict?1l:0l), x.length, y.length, MONTE_CARLO_ITERATIONS);
+ return exactP(kolmogorovSmirnovStatistic(xa, ya), x.length, y.length, strict);
}
return approximateP(kolmogorovSmirnovStatistic(x, y), x.length, y.length);
}
@@ -1000,16 +988,8 @@ public class KolmogorovSmirnovTest {
* 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}.
+ * The returned probability is exact, implemented by unwinding the recursive function
+ * definitions presented in [4] (class javadoc).
* </p>
*
* @param d D-statistic value
@@ -1020,50 +1000,10 @@ public class KolmogorovSmirnovTest {
* greater than (resp. greater than or equal to) {@code d}
*/
public double exactP(double d, int n, int m, boolean strict) {
- return integralExactP(calculateIntegralD(d, n, m, strict), n, m);
- }
-
- /**
- * Computes \(P(D_{n,m} >= d/(n*m))\), where \(D_{n,m}\) is the
- * 2-sample Kolmogorov-Smirnov statistic.
- * <p>
- * Here d is the D-statistic represented as long value.
- * The real D-statistic is obtained by dividing d by n*m.
- * See also {@link #exactP(double, int, int, boolean)}.
- *
- * @param d integral D-statistic
- * @param n first sample size
- * @param m second sample size
- * @return probability that a randomly selected m-n partition of m + n generates \(D_{n,m}\)
- * greater than or equal to {@code d/(n*m)}
- */
- private double integralExactP(long d, int n, int m) {
- 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 long curD = integralKolmogorovSmirnovStatistic(nSet, mSet);
- if (curD >= d) {
- tail++;
- }
- }
- return (double) tail / (double) CombinatoricsUtils.binomialCoefficient(n + m, n);
+ return 1 - n(m, n, m, n, calculateIntegralD(d, m, n, strict), strict) /
+ CombinatoricsUtils.binomialCoefficientDouble(n + m, m);
}
-
/**
* Uses the Kolmogorov-Smirnov distribution to approximate \(P(D_{n,m} > d)\) where \(D_{n,m}\)
* is the 2-sample Kolmogorov-Smirnov statistic. See
@@ -1270,4 +1210,61 @@ public class KolmogorovSmirnovTest {
data[i] += dist.sample();
}
}
+
+ /**
+ * The function C(i, j) defined in [4] (class javadoc), formula (5.5).
+ * defined to return 1 if |i/n - j/m| <= c; 0 otherwise. Here c is scaled up
+ * and recoded as a long to avoid rounding errors in comparison tests, so what
+ * is actually tested is |im - jn| <= cmn.
+ *
+ * @param i first path parameter
+ * @param j second path paramter
+ * @param m first sample size
+ * @param n second sample size
+ * @param cmn integral D-statistic (see {@link #calculateIntegralD(double, int, int, boolean)})
+ * @param strict whether or not the null hypothesis uses strict inequality
+ * @return C(i,j) for given m, n, c
+ */
+ private static int c(int i, int j, int m, int n, long cmn, boolean strict) {
+ if (strict) {
+ return FastMath.abs(i*(long)n - j*(long)m) <= cmn ? 1 : 0;
+ }
+ return FastMath.abs(i*(long)n - j*(long)m) < cmn ? 1 : 0;
+ }
+
+ /**
+ * The function N(i, j) defined in [4] (class javadoc).
+ * Returns the number of paths over the lattice {(i,j) : 0 <= i <= n, 0 <= j <= m}
+ * from (0,0) to (i,j) satisfying C(h,k, m, n, c) = 1 for each (h,k) on the path.
+ * The return value is integral, but subject to overflow, so it is maintained and
+ * returned as a double.
+ *
+ * @param i first path parameter
+ * @param j second path parameter
+ * @param m first sample size
+ * @param n second sample size
+ * @param cnm integral D-statistic (see {@link #calculateIntegralD(double, int, int, boolean)})
+ * @param strict whether or not the null hypothesis uses strict inequality
+ * @return number or paths to (i, j) from (0,0) representing D-values as large as c for given m, n
+ */
+ private static double n(int i, int j, int m, int n, long cnm, boolean strict) {
+ /*
+ * Unwind the recursive definition given in [4].
+ * Compute n(1,1), n(1,2)...n(2,1), n(2,2)... up to n(i,j), one row at a time.
+ * When n(i,*) are being computed, lag[] holds the values of n(i - 1, *).
+ */
+ final double[] lag = new double[n];
+ double last = 0;
+ for (int k = 0; k < n; k++) {
+ lag[k] = c(0, k + 1, m, n, cnm, strict);
+ }
+ for (int k = 1; k <= i; k++) {
+ last = c(k, 0, m, n, cnm, strict);
+ for (int l = 1; l <= j; l++) {
+ lag[l - 1] = c(k, l, m, n, cnm, strict) * (last + lag[l - 1]);
+ last = lag[l - 1];
+ }
+ }
+ return last;
+ }
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/2983ff81/src/test/R/KolmogorovSmirnovTestCases.R
----------------------------------------------------------------------
diff --git a/src/test/R/KolmogorovSmirnovTestCases.R b/src/test/R/KolmogorovSmirnovTestCases.R
index 3146325..b456291 100644
--- a/src/test/R/KolmogorovSmirnovTestCases.R
+++ b/src/test/R/KolmogorovSmirnovTestCases.R
@@ -156,6 +156,12 @@ uniform <- c(0.7930305, 0.6424382, 0.8747699, 0.7156518, 0.1845909, 0.2022326,
smallSample1 <- c(6, 7, 9, 13, 19, 21, 22, 23, 24)
smallSample2 <- c(10, 11, 12, 16, 20, 27, 28, 32, 44, 54)
+smallSample3 <- c(6, 7, 9, 13, 19, 21, 22, 23, 24, 29, 30, 34, 36, 41, 45, 47, 51, 63, 33, 91)
+smallSample4 <- c(10, 11, 12, 16, 20, 27, 28, 32, 44, 54, 56, 57, 64, 69, 71, 80, 81, 88, 90)
+smallSample5 <- c(-10, -5, 17, 21, 22, 23, 24, 30, 44, 50, 56, 57, 59, 67, 73, 75, 77, 78, 79, 80, 81, 83, 84, 85, 88, 90,
+ 92, 93, 94, 95, 98, 100, 101, 103, 105, 110)
+smallSample6 <- c(-2, -1, 0, 10, 14, 15, 16, 20, 25, 26, 27, 31, 32, 33, 34, 45, 47, 48, 51, 52, 53, 54, 60, 61, 62, 63,
+ 74, 82, 106, 107, 109, 11, 112, 113, 114)
bootSample1 <- c(0, 2, 4, 6, 8, 8, 10, 15, 22, 30, 33, 36, 38)
bootSample2 <- c(9, 17, 20, 33, 40, 51, 60, 60, 72, 90, 101)
roundingSample1 <- c(2,4,6,8,9,10,11,12,13)
@@ -185,7 +191,13 @@ verifyTwoSampleLargeSamples(gaussian, gaussian2, 0.0319983962391632, 0.202352941
"Two sample N(0, 1) vs N(0, 1.6)")
verifyTwoSampleSmallSamplesExact(smallSample1, smallSample2, 0.105577085453247, .5, tol,
-"Two sample small samples exact")
+"Two sample small samples exact 1")
+
+verifyTwoSampleSmallSamplesExact(smallSample3, smallSample4, 0.046298660942952, 0.426315789473684, tol,
+"Two sample small samples exact 2")
+
+verifyTwoSampleSmallSamplesExact(smallSample5, smallSample6, 0.00300743602233366, 0.41031746031746, tol,
+"Two sample small samples exact 3")
verifyTwoSampleBootstrap(bootSample1, bootSample2, 0.0059, 1E-3, "Two sample bootstrap - isolated failures possible")
verifyTwoSampleBootstrap(gaussian, gaussian2, 0.0237, 1E-2, "Two sample bootstrap - isolated failures possible")
http://git-wip-us.apache.org/repos/asf/commons-math/blob/2983ff81/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java b/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java
index 0a8a419..5c8df3f 100644
--- a/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java
+++ b/src/test/java/org/apache/commons/math3/stat/inference/KolmogorovSmirnovTestTest.java
@@ -25,7 +25,6 @@ import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.distribution.UniformRealDistribution;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.Well19937c;
-import org.apache.commons.math3.stat.inference.KolmogorovSmirnovTest;
import org.apache.commons.math3.util.CombinatoricsUtils;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
@@ -181,12 +180,48 @@ public class KolmogorovSmirnovTestTest {
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
+ // Reference values from R, version 3.2.0 - 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);
}
+ /** Small samples - exact p-value, checked against R */
+ @Test
+ public void testTwoSampleSmallSampleExact2() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final double[] smallSample1 = {
+ 6, 7, 9, 13, 19, 21, 22, 23, 24, 29, 30, 34, 36, 41, 45, 47, 51, 63, 33, 91
+ };
+ final double[] smallSample2 = {
+ 10, 11, 12, 16, 20, 27, 28, 32, 44, 54, 56, 57, 64, 69, 71, 80, 81, 88, 90
+ };
+ // Reference values from R, version 3.2.0 - R uses non-strict inequality in null hypothesis
+ Assert
+ .assertEquals(0.0462986609, test.kolmogorovSmirnovTest(smallSample1, smallSample2, false), TOLERANCE);
+ Assert.assertEquals(0.4263157895, test.kolmogorovSmirnovStatistic(smallSample1, smallSample2), TOLERANCE);
+ }
+
+ /** Small samples - exact p-value, checked against R */
+ @Test
+ public void testTwoSampleSmallSampleExact3() {
+ final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest();
+ final double[] smallSample1 = {
+ -10, -5, 17, 21, 22, 23, 24, 30, 44, 50, 56, 57, 59, 67, 73, 75, 77, 78, 79, 80, 81, 83, 84, 85, 88, 90,
+ 92, 93, 94, 95, 98, 100, 101, 103, 105, 110
+ };
+ final double[] smallSample2 = {
+ -2, -1, 0, 10, 14, 15, 16, 20, 25, 26, 27, 31, 32, 33, 34, 45, 47, 48, 51, 52, 53, 54, 60, 61, 62, 63,
+ 74, 82, 106, 107, 109, 11, 112, 113, 114
+ };
+ // Reference values from R, version 3.2.0 - R uses non-strict inequality in null hypothesis
+ Assert
+ .assertEquals(0.00300743602, test.kolmogorovSmirnovTest(smallSample1, smallSample2, false), TOLERANCE);
+ Assert.assertEquals(0.4103174603, test.kolmogorovSmirnovStatistic(smallSample1, smallSample2), TOLERANCE);
+ Assert
+ .assertEquals(0.00300743602, test.kolmogorovSmirnovTest(smallSample2, smallSample1, false), 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.
@@ -279,10 +314,9 @@ public class KolmogorovSmirnovTestTest {
}
/**
- * 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.
+ * Verifies that Monte Carlo simulation gives results close to exact p values.
*/
- // @Test
+ @Test
public void testTwoSampleMonteCarlo() {
final KolmogorovSmirnovTest test = new KolmogorovSmirnovTest(new Well19937c(1000));
final int sampleSize = 14;