commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From t.@apache.org
Subject svn commit: r1242164 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/stat/inference/ test/java/org/apache/commons/math/stat/inference/
Date Wed, 08 Feb 2012 22:58:52 GMT
Author: tn
Date: Wed Feb  8 22:58:51 2012
New Revision: 1242164

URL: http://svn.apache.org/viewvc?rev=1242164&view=rev
Log:
Merged interface and implementation of OneWayAnova, MannWhitneyUTest and WilcoxonSignedRankTest.
JIRA: MATH-739

Removed:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTestImpl.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestImpl.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTest.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnova.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/TestUtils.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/MannWhitneyUTestTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/OneWayAnovaTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/TestUtilsTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTest.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/MannWhitneyUTest.java
Wed Feb  8 22:58:51 2012
@@ -16,17 +16,84 @@
  */
 package org.apache.commons.math.stat.inference;
 
+import org.apache.commons.math.distribution.NormalDistribution;
 import org.apache.commons.math.exception.ConvergenceException;
 import org.apache.commons.math.exception.MaxCountExceededException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.NullArgumentException;
+import org.apache.commons.math.stat.ranking.NaNStrategy;
+import org.apache.commons.math.stat.ranking.NaturalRanking;
+import org.apache.commons.math.stat.ranking.TiesStrategy;
+import org.apache.commons.math.util.FastMath;
 
 /**
- * An interface for Mann-Whitney U test (also called Wilcoxon rank-sum test).
+ * An implementation of the Mann-Whitney U test (also called Wilcoxon rank-sum test).
  *
  * @version $Id$
  */
-public interface MannWhitneyUTest {
+public class MannWhitneyUTest {
+
+    /** Ranking algorithm. */
+    private NaturalRanking naturalRanking;
+
+    /**
+     * Create a test instance using where NaN's are left in place and ties get
+     * the average of applicable ranks. Use this unless you are very sure of
+     * what you are doing.
+     */
+    public MannWhitneyUTest() {
+        naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
+                TiesStrategy.AVERAGE);
+    }
+
+    /**
+     * Create a test instance using the given strategies for NaN's and ties.
+     * Only use this if you are sure of what you are doing.
+     *
+     * @param nanStrategy
+     *            specifies the strategy that should be used for Double.NaN's
+     * @param tiesStrategy
+     *            specifies the strategy that should be used for ties
+     */
+    public MannWhitneyUTest(final NaNStrategy nanStrategy,
+                            final TiesStrategy tiesStrategy) {
+        naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
+    }
+
+    /**
+     * Ensures that the provided arrays fulfills the assumptions.
+     *
+     * @param x first sample
+     * @param y second sample
+     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
+     * @throws NoDataException if {@code x} or {@code y} are zero-length.
+     */
+    private void ensureDataConformance(final double[] x, final double[] y)
+        throws NullArgumentException, NoDataException {
+
+        if (x == null ||
+            y == null) {
+            throw new NullArgumentException();
+        }
+        if (x.length == 0 ||
+            y.length == 0) {
+            throw new NoDataException();
+        }
+    }
+
+    /** Concatenate the samples into one array.
+     * @param x first sample
+     * @param y second sample
+     * @return concatenated array
+     */
+    private double[] concatenateSamples(final double[] x, final double[] y) {
+        final double[] z = new double[x.length + y.length];
+
+        System.arraycopy(x, 0, z, 0, x.length);
+        System.arraycopy(y, 0, z, x.length, y.length);
+
+        return z;
+    }
 
     /**
      * Computes the <a
@@ -56,8 +123,65 @@ public interface MannWhitneyUTest {
      * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
      * @throws NoDataException if {@code x} or {@code y} are zero-length.
      */
-    double mannWhitneyU(final double[] x, final double[] y)
-        throws NullArgumentException, NoDataException;
+    public double mannWhitneyU(final double[] x, final double[] y)
+        throws NullArgumentException, NoDataException {
+
+        ensureDataConformance(x, y);
+
+        final double[] z = concatenateSamples(x, y);
+        final double[] ranks = naturalRanking.rank(z);
+
+        double sumRankX = 0;
+
+        /*
+         * The ranks for x is in the first x.length entries in ranks because x
+         * is in the first x.length entries in z
+         */
+        for (int i = 0; i < x.length; ++i) {
+            sumRankX += ranks[i];
+        }
+
+        /*
+         * U1 = R1 - (n1 * (n1 + 1)) / 2 where R1 is sum of ranks for sample 1,
+         * e.g. x, n1 is the number of observations in sample 1.
+         */
+        final double U1 = sumRankX - (x.length * (x.length + 1)) / 2;
+
+        /*
+         * It can be shown that U1 + U2 = n1 * n2
+         */
+        final double U2 = x.length * y.length - U1;
+
+        return FastMath.max(U1, U2);
+    }
+
+    /**
+     * @param Umin smallest Mann-Whitney U value
+     * @param n1 number of subjects in first sample
+     * @param n2 number of subjects in second sample
+     * @return two-sided asymptotic p-value
+     * @throws ConvergenceException if the p-value can not be computed
+     * due to a convergence error
+     * @throws MaxCountExceededException if the maximum number of
+     * iterations is exceeded
+     */
+    private double calculateAsymptoticPValue(final double Umin,
+                                             final int n1,
+                                             final int n2)
+        throws ConvergenceException, MaxCountExceededException {
+
+        final int n1n2prod = n1 * n2;
+
+        // http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U#Normal_approximation
+        final double EU = (double) n1n2prod / 2.0;
+        final double VarU = (double) (n1n2prod * (n1 + n2 + 1)) / 12.0;
+
+        final double z = (Umin - EU) / FastMath.sqrt(VarU);
+
+        final NormalDistribution standardNormal = new NormalDistribution(0, 1);
+
+        return 2 * standardNormal.cumulativeProbability(z);
+    }
 
     /**
      * Returns the asymptotic <i>observed significance level</i>, or <a href=
@@ -76,7 +200,10 @@ public interface MannWhitneyUTest {
      * <li>All observations in the two samples are independent.</li>
      * <li>The observations are at least ordinal (continuous are also ordinal).</li>
      * </ul>
-     * </p>
+     * </p><p>
+     * Ties give rise to biased variance at the moment. See e.g. <a
+     * href="http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf"
+     * >http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf</a>.</p>
      *
      * @param x the first sample
      * @param y the second sample
@@ -88,7 +215,20 @@ public interface MannWhitneyUTest {
      * @throws MaxCountExceededException if the maximum number of iterations
      * is exceeded
      */
-    double mannWhitneyUTest(final double[] x, final double[] y)
+    public double mannWhitneyUTest(final double[] x, final double[] y)
         throws NullArgumentException, NoDataException,
-        ConvergenceException, MaxCountExceededException;
+        ConvergenceException, MaxCountExceededException {
+
+        ensureDataConformance(x, y);
+
+        final double Umax = mannWhitneyU(x, y);
+
+        /*
+         * It can be shown that U1 + U2 = n1 * n2
+         */
+        final double Umin = x.length * y.length - Umax;
+
+        return calculateAsymptoticPValue(Umin, x.length, y.length);
+    }
+
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnova.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnova.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnova.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/OneWayAnova.java
Wed Feb  8 22:58:51 2012
@@ -16,27 +16,46 @@
  */
 package org.apache.commons.math.stat.inference;
 
+import org.apache.commons.math.distribution.FDistribution;
 import org.apache.commons.math.exception.ConvergenceException;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.MaxCountExceededException;
 import org.apache.commons.math.exception.NullArgumentException;
 import org.apache.commons.math.exception.OutOfRangeException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.stat.descriptive.summary.Sum;
+import org.apache.commons.math.stat.descriptive.summary.SumOfSquares;
 
 import java.util.Collection;
 
 /**
- * An interface for one-way ANOVA (analysis of variance).
+ * Implements one-way ANOVA (analysis of variance) statistics.
  *
  * <p> Tests for differences between two or more categories of univariate data
  * (for example, the body mass index of accountants, lawyers, doctors and
  * computer programmers).  When two categories are given, this is equivalent to
  * the {@link org.apache.commons.math.stat.inference.TTest}.
- * </p>
+ * </p><p>
+ * Uses the {@link org.apache.commons.math.distribution.FDistribution
+ * commons-math F Distribution implementation} to estimate exact p-values.</p>
+ * <p>This implementation is based on a description at
+ * http://faculty.vassar.edu/lowry/ch13pt1.html</p>
+ * <pre>
+ * Abbreviations: bg = between groups,
+ *                wg = within groups,
+ *                ss = sum squared deviations
+ * </pre>
  *
  * @since 1.2
  * @version $Id$
  */
-public interface OneWayAnova {
+public class OneWayAnova {
+
+    /**
+     * Default constructor.
+     */
+    public OneWayAnova() {
+    }
 
     /**
      * Computes the ANOVA F-value for a collection of <code>double[]</code>
@@ -47,7 +66,15 @@ public interface OneWayAnova {
      * <code>double[]</code> arrays.</li>
      * <li> There must be at least two <code>double[]</code> arrays in
the
      * <code>categoryData</code> collection and each of these arrays must
-     * contain at least two values.</li></ul></p>
+     * contain at least two values.</li></ul></p><p>
+     * This implementation computes the F statistic using the definitional
+     * formula<pre>
+     *   F = msbg/mswg</pre>
+     * where<pre>
+     *  msbg = between group mean square
+     *  mswg = within group mean square</pre>
+     * are as defined <a href="http://faculty.vassar.edu/lowry/ch13pt1.html">
+     * here</a></p>
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
@@ -57,8 +84,13 @@ public interface OneWayAnova {
      * array is less than 2 or a contained <code>double[]</code> array does not
have
      * at least two values
      */
-    double anovaFValue(Collection<double[]> categoryData)
-        throws NullArgumentException, DimensionMismatchException;
+    public double anovaFValue(final Collection<double[]> categoryData)
+        throws NullArgumentException, DimensionMismatchException {
+
+        AnovaStats a = anovaStats(categoryData);
+        return a.F;
+
+    }
 
     /**
      * Computes the ANOVA P-value for a collection of <code>double[]</code>
@@ -69,7 +101,14 @@ public interface OneWayAnova {
      * <code>double[]</code> arrays.</li>
      * <li> There must be at least two <code>double[]</code> arrays in
the
      * <code>categoryData</code> collection and each of these arrays must
-     * contain at least two values.</li></ul></p>
+     * contain at least two values.</li></ul></p><p>
+     * This implementation uses the
+     * {@link org.apache.commons.math.distribution.FDistribution
+     * commons-math F Distribution implementation} to estimate the exact
+     * p-value, using the formula<pre>
+     *   p = 1 - cumulativeProbability(F)</pre>
+     * where <code>F</code> is the F value and <code>cumulativeProbability</code>
+     * is the commons-math implementation of the F distribution.</p>
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
@@ -81,9 +120,15 @@ public interface OneWayAnova {
      * @throws ConvergenceException if the p-value can not be computed due to a convergence
error
      * @throws MaxCountExceededException if the maximum number of iterations is exceeded
      */
-    double anovaPValue(Collection<double[]> categoryData)
+    public double anovaPValue(final Collection<double[]> categoryData)
         throws NullArgumentException, DimensionMismatchException,
-        ConvergenceException, MaxCountExceededException;
+        ConvergenceException, MaxCountExceededException {
+
+        AnovaStats a = anovaStats(categoryData);
+        FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
+        return 1.0 - fdist.cumulativeProbability(a.F);
+
+    }
 
     /**
      * Performs an ANOVA test, evaluating the null hypothesis that there
@@ -96,7 +141,15 @@ public interface OneWayAnova {
      * <code>categoryData</code> collection and each of these arrays must
      * contain at least two values.</li>
      * <li>alpha must be strictly greater than 0 and less than or equal to 0.5.
-     * </li></ul></p>
+     * </li></ul></p><p>
+     * This implementation uses the
+     * {@link org.apache.commons.math.distribution.FDistribution
+     * commons-math F Distribution implementation} to estimate the exact
+     * p-value, using the formula<pre>
+     *   p = 1 - cumulativeProbability(F)</pre>
+     * where <code>F</code> is the F value and <code>cumulativeProbability</code>
+     * is the commons-math implementation of the F distribution.</p>
+     * <p>True is returned iff the estimated p-value is less than alpha.</p>
      *
      * @param categoryData <code>Collection</code> of <code>double[]</code>
      * arrays each containing data for one category
@@ -111,8 +164,120 @@ public interface OneWayAnova {
      * @throws ConvergenceException if the p-value can not be computed due to a convergence
error
      * @throws MaxCountExceededException if the maximum number of iterations is exceeded
      */
-    boolean anovaTest(Collection<double[]> categoryData, double alpha)
-        throws NullArgumentException, DimensionMismatchException, OutOfRangeException,
-        ConvergenceException, MaxCountExceededException;
+    public boolean anovaTest(final Collection<double[]> categoryData,
+                             final double alpha)
+        throws NullArgumentException, DimensionMismatchException,
+        OutOfRangeException, ConvergenceException, MaxCountExceededException {
+
+        if ((alpha <= 0) || (alpha > 0.5)) {
+            throw new OutOfRangeException(
+                    LocalizedFormats.OUT_OF_BOUND_SIGNIFICANCE_LEVEL,
+                    alpha, 0, 0.5);
+        }
+        return anovaPValue(categoryData) < alpha;
+
+    }
+
+    /**
+     * This method actually does the calculations (except P-value).
+     *
+     * @param categoryData <code>Collection</code> of <code>double[]</code>
+     * arrays each containing data for one category
+     * @return computed AnovaStats
+     * @throws NullArgumentException if <code>categoryData</code> is <code>null</code>
+     * @throws DimensionMismatchException if the length of the <code>categoryData</code>
+     * array is less than 2 or a contained <code>double[]</code> array does not
contain
+     * at least two values
+     */
+    private AnovaStats anovaStats(final Collection<double[]> categoryData)
+        throws NullArgumentException, DimensionMismatchException {
+
+        if (categoryData == null) {
+            throw new NullArgumentException();
+        }
+
+        // check if we have enough categories
+        if (categoryData.size() < 2) {
+            throw new DimensionMismatchException(
+                    LocalizedFormats.TWO_OR_MORE_CATEGORIES_REQUIRED,
+                    categoryData.size(), 2);
+        }
+
+        // check if each category has enough data and all is double[]
+        for (double[] array : categoryData) {
+            if (array.length <= 1) {
+                throw new DimensionMismatchException(
+                        LocalizedFormats.TWO_OR_MORE_VALUES_IN_CATEGORY_REQUIRED,
+                        array.length, 2);
+            }
+        }
+
+        int dfwg = 0;
+        double sswg = 0;
+        Sum totsum = new Sum();
+        SumOfSquares totsumsq = new SumOfSquares();
+        int totnum = 0;
+
+        for (double[] data : categoryData) {
+
+            Sum sum = new Sum();
+            SumOfSquares sumsq = new SumOfSquares();
+            int num = 0;
+
+            for (int i = 0; i < data.length; i++) {
+                double val = data[i];
+
+                // within category
+                num++;
+                sum.increment(val);
+                sumsq.increment(val);
+
+                // for all categories
+                totnum++;
+                totsum.increment(val);
+                totsumsq.increment(val);
+            }
+            dfwg += num - 1;
+            double ss = sumsq.getResult() - sum.getResult() * sum.getResult() / num;
+            sswg += ss;
+        }
+        double sst = totsumsq.getResult() - totsum.getResult() *
+            totsum.getResult()/totnum;
+        double ssbg = sst - sswg;
+        int dfbg = categoryData.size() - 1;
+        double msbg = ssbg/dfbg;
+        double mswg = sswg/dfwg;
+        double F = msbg/mswg;
+
+        return new AnovaStats(dfbg, dfwg, F);
+    }
+
+    /**
+        Convenience class to pass dfbg,dfwg,F values around within AnovaImpl.
+        No get/set methods provided.
+    */
+    private static class AnovaStats {
+
+        /** Degrees of freedom in numerator (between groups). */
+        private final int dfbg;
+
+        /** Degrees of freedom in denominator (within groups). */
+        private final int dfwg;
+
+        /** Statistic. */
+        private final double F;
+
+        /**
+         * Constructor
+         * @param dfbg degrees of freedom in numerator (between groups)
+         * @param dfwg degrees of freedom in denominator (within groups)
+         * @param F statistic
+         */
+        private AnovaStats(int dfbg, int dfwg, double F) {
+            this.dfbg = dfbg;
+            this.dfwg = dfwg;
+            this.F = F;
+        }
+    }
 
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/TestUtils.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/TestUtils.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/TestUtils.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/TestUtils.java
Wed Feb  8 22:58:51 2012
@@ -18,6 +18,11 @@ package org.apache.commons.math.stat.inf
 
 import java.util.Collection;
 import org.apache.commons.math.MathException;
+import org.apache.commons.math.exception.ConvergenceException;
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.MaxCountExceededException;
+import org.apache.commons.math.exception.NullArgumentException;
+import org.apache.commons.math.exception.OutOfRangeException;
 import org.apache.commons.math.stat.descriptive.StatisticalSummary;
 
 /**
@@ -29,18 +34,25 @@ import org.apache.commons.math.stat.desc
  */
 public class TestUtils  {
 
-    /** Singleton TTest instance using default implementation. */
+    /** Singleton TTest instance. */
     private static final TTest T_TEST = new TTestImpl();
 
-    /** Singleton ChiSquareTest instance using default implementation. */
+    /** Singleton ChiSquareTest instance. */
     private static final ChiSquareTest CHI_SQUARE_TEST = new ChiSquareTestImpl();
 
-    /** Singleton ChiSquareTest instance using default implementation. */
+    /** Singleton ChiSquareTest instance. */
     private static final UnknownDistributionChiSquareTest UNKNOWN_DISTRIBUTION_CHI_SQUARE_TEST
=
         new ChiSquareTestImpl();
 
-    /** Singleton OneWayAnova instance using default implementation. */
-    private static final OneWayAnova ONE_WAY_ANANOVA = new OneWayAnovaImpl();
+    /** Singleton OneWayAnova instance. */
+    private static final OneWayAnova ONE_WAY_ANANOVA = new OneWayAnova();
+
+    /** Singleton MannWhitneyUTest instance using default ranking. */
+    private static final MannWhitneyUTest MANN_WHITNEY_U_TEST = new MannWhitneyUTest();
+
+    /** Singleton WilcoxonSignedRankTest instance. */
+    private static final WilcoxonSignedRankTest WILCOXON_SIGNED_RANK_TEST =
+        new WilcoxonSignedRankTest();
 
     /**
      * Prevent instantiation.
@@ -74,7 +86,7 @@ public class TestUtils  {
     public static boolean homoscedasticTTest(double[] sample1, double[] sample2,
             double alpha)
         throws IllegalArgumentException, MathException {
-        return T_TEST. homoscedasticTTest(sample1, sample2, alpha);
+        return T_TEST.homoscedasticTTest(sample1, sample2, alpha);
     }
 
     /**
@@ -174,7 +186,7 @@ public class TestUtils  {
     public static boolean tTest(double mu, StatisticalSummary sampleStats,
         double alpha)
         throws IllegalArgumentException, MathException {
-        return T_TEST. tTest(mu, sampleStats, alpha);
+        return T_TEST.tTest(mu, sampleStats, alpha);
     }
 
     /**
@@ -207,7 +219,7 @@ public class TestUtils  {
     public static boolean tTest(StatisticalSummary sampleStats1,
         StatisticalSummary sampleStats2, double alpha)
         throws IllegalArgumentException, MathException {
-        return T_TEST. tTest(sampleStats1, sampleStats2, alpha);
+        return T_TEST.tTest(sampleStats1, sampleStats2, alpha);
     }
 
     /**
@@ -257,7 +269,7 @@ public class TestUtils  {
      */
     public static boolean chiSquareTest(long[][] counts, double alpha)
         throws IllegalArgumentException, MathException {
-        return CHI_SQUARE_TEST. chiSquareTest(counts, alpha);
+        return CHI_SQUARE_TEST.chiSquareTest(counts, alpha);
     }
 
     /**
@@ -305,8 +317,8 @@ public class TestUtils  {
      *
      * @since 1.2
      */
-    public static double oneWayAnovaFValue(Collection<double[]> categoryData)
-    throws IllegalArgumentException, MathException {
+    public static double oneWayAnovaFValue(final Collection<double[]> categoryData)
+        throws NullArgumentException, DimensionMismatchException {
         return ONE_WAY_ANANOVA.anovaFValue(categoryData);
     }
 
@@ -315,8 +327,9 @@ public class TestUtils  {
      *
      * @since 1.2
      */
-    public static double oneWayAnovaPValue(Collection<double[]> categoryData)
-    throws IllegalArgumentException, MathException {
+    public static double oneWayAnovaPValue(final Collection<double[]> categoryData)
+        throws NullArgumentException, DimensionMismatchException,
+        ConvergenceException, MaxCountExceededException {
         return ONE_WAY_ANANOVA.anovaPValue(categoryData);
     }
 
@@ -325,8 +338,10 @@ public class TestUtils  {
      *
      * @since 1.2
      */
-    public static boolean oneWayAnovaTest(Collection<double[]> categoryData, double
alpha)
-    throws IllegalArgumentException, MathException {
+    public static boolean oneWayAnovaTest(final Collection<double[]> categoryData,
+                                          final double alpha)
+        throws NullArgumentException, DimensionMismatchException,
+        OutOfRangeException, ConvergenceException, MaxCountExceededException {
         return ONE_WAY_ANANOVA.anovaTest(categoryData, alpha);
     }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTest.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTest.java
Wed Feb  8 22:58:51 2012
@@ -16,19 +16,123 @@
  */
 package org.apache.commons.math.stat.inference;
 
+import org.apache.commons.math.distribution.NormalDistribution;
 import org.apache.commons.math.exception.ConvergenceException;
 import org.apache.commons.math.exception.DimensionMismatchException;
 import org.apache.commons.math.exception.MaxCountExceededException;
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.NullArgumentException;
 import org.apache.commons.math.exception.NumberIsTooLargeException;
+import org.apache.commons.math.stat.ranking.NaNStrategy;
+import org.apache.commons.math.stat.ranking.NaturalRanking;
+import org.apache.commons.math.stat.ranking.TiesStrategy;
+import org.apache.commons.math.util.FastMath;
 
 /**
- * An interface for Wilcoxon signed-rank test.
+ * An implementation of the Wilcoxon signed-rank test.
  *
  * @version $Id$
  */
-public interface WilcoxonSignedRankTest {
+public class WilcoxonSignedRankTest {
+
+    /** Ranking algorithm. */
+    private NaturalRanking naturalRanking;
+
+    /**
+     * Create a test instance where NaN's are left in place and ties get
+     * the average of applicable ranks. Use this unless you are very sure
+     * of what you are doing.
+     */
+    public WilcoxonSignedRankTest() {
+        naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
+                TiesStrategy.AVERAGE);
+    }
+
+    /**
+     * Create a test instance using the given strategies for NaN's and ties.
+     * Only use this if you are sure of what you are doing.
+     *
+     * @param nanStrategy
+     *            specifies the strategy that should be used for Double.NaN's
+     * @param tiesStrategy
+     *            specifies the strategy that should be used for ties
+     */
+    public WilcoxonSignedRankTest(final NaNStrategy nanStrategy,
+                                  final TiesStrategy tiesStrategy) {
+        naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
+    }
+
+    /**
+     * Ensures that the provided arrays fulfills the assumptions.
+     *
+     * @param x first sample
+     * @param y second sample
+     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
+     * @throws NoDataException if {@code x} or {@code y} are zero-length.
+     * @throws DimensionMismatchException if {@code x} and {@code y} do not
+     * have the same length.
+     */
+    private void ensureDataConformance(final double[] x, final double[] y)
+        throws NullArgumentException, NoDataException, DimensionMismatchException {
+
+        if (x == null ||
+            y == null) {
+                throw new NullArgumentException();
+        }
+        if (x.length == 0 ||
+            y.length == 0) {
+            throw new NoDataException();
+        }
+        if (y.length != x.length) {
+            throw new DimensionMismatchException(y.length, x.length);
+        }
+    }
+
+    /**
+     * Calculates y[i] - x[i] for all i
+     *
+     * @param x first sample
+     * @param y second sample
+     * @return z = y - x
+     */
+    private double[] calculateDifferences(final double[] x, final double[] y) {
+
+        final double[] z = new double[x.length];
+
+        for (int i = 0; i < x.length; ++i) {
+            z[i] = y[i] - x[i];
+        }
+
+        return z;
+    }
+
+    /**
+     * Calculates |z[i]| for all i
+     *
+     * @param z sample
+     * @return |z|
+     * @throws NullArgumentException if {@code z} is {@code null}
+     * @throws NoDataException if {@code z} is zero-length.
+     */
+    private double[] calculateAbsoluteDifferences(final double[] z)
+        throws NullArgumentException, NoDataException {
+
+        if (z == null) {
+            throw new NullArgumentException();
+        }
+
+        if (z.length == 0) {
+            throw new NoDataException();
+        }
+
+        final double[] zAbs = new double[z.length];
+
+        for (int i = 0; i < z.length; ++i) {
+            zAbs[i] = FastMath.abs(z[i]);
+        }
+
+        return zAbs;
+    }
 
     /**
      * Computes the <a
@@ -65,8 +169,94 @@ public interface WilcoxonSignedRankTest 
      * @throws DimensionMismatchException if {@code x} and {@code y} do not
      * have the same length.
      */
-    double wilcoxonSignedRank(final double[] x, final double[] y)
-        throws NullArgumentException, NoDataException, DimensionMismatchException;
+    public double wilcoxonSignedRank(final double[] x, final double[] y)
+        throws NullArgumentException, NoDataException, DimensionMismatchException {
+
+        ensureDataConformance(x, y);
+
+        // throws IllegalArgumentException if x and y are not correctly
+        // specified
+        final double[] z = calculateDifferences(x, y);
+        final double[] zAbs = calculateAbsoluteDifferences(z);
+
+        final double[] ranks = naturalRanking.rank(zAbs);
+
+        double Wplus = 0;
+
+        for (int i = 0; i < z.length; ++i) {
+            if (z[i] > 0) {
+                Wplus += ranks[i];
+            }
+        }
+
+        final int N = x.length;
+        final double Wminus = (((double) (N * (N + 1))) / 2.0) - Wplus;
+
+        return FastMath.max(Wplus, Wminus);
+    }
+
+    /**
+     * Algorithm inspired by
+     * http://www.fon.hum.uva.nl/Service/Statistics/Signed_Rank_Algorihms.html#C
+     * by Rob van Son, Institute of Phonetic Sciences & IFOTT,
+     * University of Amsterdam
+     *
+     * @param Wmax largest Wilcoxon signed rank value
+     * @param N number of subjects (corresponding to x.length)
+     * @return two-sided exact p-value
+     */
+    private double calculateExactPValue(final double Wmax, final int N) {
+
+        // Total number of outcomes (equal to 2^N but a lot faster)
+        final int m = 1 << N;
+
+        int largerRankSums = 0;
+
+        for (int i = 0; i < m; ++i) {
+            int rankSum = 0;
+
+            // Generate all possible rank sums
+            for (int j = 0; j < N; ++j) {
+
+                // (i >> j) & 1 extract i's j-th bit from the right
+                if (((i >> j) & 1) == 1) {
+                    rankSum += j + 1;
+                }
+            }
+
+            if (rankSum >= Wmax) {
+                ++largerRankSums;
+            }
+        }
+
+        /*
+         * largerRankSums / m gives the one-sided p-value, so it's multiplied
+         * with 2 to get the two-sided p-value
+         */
+        return 2 * ((double) largerRankSums) / ((double) m);
+    }
+
+    /**
+     * @param Wmin smallest Wilcoxon signed rank value
+     * @param N number of subjects (corresponding to x.length)
+     * @return two-sided asymptotic p-value
+     */
+    private double calculateAsymptoticPValue(final double Wmin, final int N) {
+
+        final double ES = (double) (N * (N + 1)) / 4.0;
+
+        /* Same as (but saves computations):
+         * final double VarW = ((double) (N * (N + 1) * (2*N + 1))) / 24;
+         */
+        final double VarS = ES * ((double) (2 * N + 1) / 6.0);
+
+        // - 0.5 is a continuity correction
+        final double z = (Wmin - ES - 0.5) / FastMath.sqrt(VarS);
+
+        final NormalDistribution standardNormal = new NormalDistribution(0, 1);
+
+        return 2*standardNormal.cumulativeProbability(z);
+    }
 
     /**
      * Returns the <i>observed significance level</i>, or <a href=
@@ -110,7 +300,25 @@ public interface WilcoxonSignedRankTest 
      * @throws MaxCountExceededException if the maximum number of iterations
      * is exceeded
      */
-    double wilcoxonSignedRankTest(final double[] x, final double[] y, boolean exactPValue)
+    public double wilcoxonSignedRankTest(final double[] x, final double[] y,
+                                         final boolean exactPValue)
         throws NullArgumentException, NoDataException, DimensionMismatchException,
-        NumberIsTooLargeException, ConvergenceException, MaxCountExceededException;
+        NumberIsTooLargeException, ConvergenceException, MaxCountExceededException {
+
+        ensureDataConformance(x, y);
+
+        final int N = x.length;
+        final double Wmax = wilcoxonSignedRank(x, y);
+
+        if (exactPValue && N > 30) {
+            throw new NumberIsTooLargeException(N, 30, true);
+        }
+
+        if (exactPValue) {
+            return calculateExactPValue(Wmax, N);
+        } else {
+            final double Wmin = ( (double)(N*(N+1)) / 2.0 ) - Wmax;
+            return calculateAsymptoticPValue(Wmin, N);
+        }
+    }
 }

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/MannWhitneyUTestTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/MannWhitneyUTestTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/MannWhitneyUTestTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/MannWhitneyUTestTest.java
Wed Feb  8 22:58:51 2012
@@ -30,7 +30,7 @@ import org.junit.Test;
 
 public class MannWhitneyUTestTest {
 
-    protected MannWhitneyUTest testStatistic = new MannWhitneyUTestImpl();
+    protected MannWhitneyUTest testStatistic = new MannWhitneyUTest();
 
     @Test
     public void testMannWhitneyUSimple() throws Exception {

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/OneWayAnovaTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/OneWayAnovaTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/OneWayAnovaTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/OneWayAnovaTest.java
Wed Feb  8 22:58:51 2012
@@ -32,7 +32,7 @@ import org.junit.Test;
 
 public class OneWayAnovaTest {
 
-    protected OneWayAnova testStatistic = new OneWayAnovaImpl();
+    protected OneWayAnova testStatistic = new OneWayAnova();
 
     private double[] emptyArray = {};
 

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/TestUtilsTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/TestUtilsTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/TestUtilsTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/TestUtilsTest.java
Wed Feb  8 22:58:51 2012
@@ -450,7 +450,7 @@ public class TestUtilsTest {
       {110.0, 115.0, 111.0, 117.0, 128.0};
 
     private List<double[]> classes = new ArrayList<double[]>();
-    private OneWayAnova oneWayAnova = new OneWayAnovaImpl();
+    private OneWayAnova oneWayAnova = new OneWayAnova();
 
     @Test
     public void testOneWayAnovaUtils() throws Exception {

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestTest.java?rev=1242164&r1=1242163&r2=1242164&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/stat/inference/WilcoxonSignedRankTestTest.java
Wed Feb  8 22:58:51 2012
@@ -25,14 +25,14 @@ import org.junit.Test;
 
 
 /**
- * Test cases for the ChiSquareTestImpl class.
+ * Test cases for the WilcoxonSignedRangTest class.
  *
  * @version $Id$
  */
 
 public class WilcoxonSignedRankTestTest {
 
-    protected WilcoxonSignedRankTest testStatistic = new WilcoxonSignedRankTestImpl();
+    protected WilcoxonSignedRankTest testStatistic = new WilcoxonSignedRankTest();
 
     @Test
     public void testWilcoxonSignedRankSimple() throws Exception {



Mime
View raw message