001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.stat.inference;
018
019import org.apache.commons.math3.distribution.NormalDistribution;
020import org.apache.commons.math3.exception.ConvergenceException;
021import org.apache.commons.math3.exception.MaxCountExceededException;
022import org.apache.commons.math3.exception.NoDataException;
023import org.apache.commons.math3.exception.NullArgumentException;
024import org.apache.commons.math3.stat.ranking.NaNStrategy;
025import org.apache.commons.math3.stat.ranking.NaturalRanking;
026import org.apache.commons.math3.stat.ranking.TiesStrategy;
027import org.apache.commons.math3.util.FastMath;
028
029/**
030 * An implementation of the Mann-Whitney U test (also called Wilcoxon rank-sum test).
031 *
032 */
033public class MannWhitneyUTest {
034
035    /** Ranking algorithm. */
036    private NaturalRanking naturalRanking;
037
038    /**
039     * Create a test instance using where NaN's are left in place and ties get
040     * the average of applicable ranks. Use this unless you are very sure of
041     * what you are doing.
042     */
043    public MannWhitneyUTest() {
044        naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
045                TiesStrategy.AVERAGE);
046    }
047
048    /**
049     * Create a test instance using the given strategies for NaN's and ties.
050     * Only use this if you are sure of what you are doing.
051     *
052     * @param nanStrategy
053     *            specifies the strategy that should be used for Double.NaN's
054     * @param tiesStrategy
055     *            specifies the strategy that should be used for ties
056     */
057    public MannWhitneyUTest(final NaNStrategy nanStrategy,
058                            final TiesStrategy tiesStrategy) {
059        naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
060    }
061
062    /**
063     * Ensures that the provided arrays fulfills the assumptions.
064     *
065     * @param x first sample
066     * @param y second sample
067     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
068     * @throws NoDataException if {@code x} or {@code y} are zero-length.
069     */
070    private void ensureDataConformance(final double[] x, final double[] y)
071        throws NullArgumentException, NoDataException {
072
073        if (x == null ||
074            y == null) {
075            throw new NullArgumentException();
076        }
077        if (x.length == 0 ||
078            y.length == 0) {
079            throw new NoDataException();
080        }
081    }
082
083    /** Concatenate the samples into one array.
084     * @param x first sample
085     * @param y second sample
086     * @return concatenated array
087     */
088    private double[] concatenateSamples(final double[] x, final double[] y) {
089        final double[] z = new double[x.length + y.length];
090
091        System.arraycopy(x, 0, z, 0, x.length);
092        System.arraycopy(y, 0, z, x.length, y.length);
093
094        return z;
095    }
096
097    /**
098     * Computes the <a
099     * href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney
100     * U statistic</a> comparing mean for two independent samples possibly of
101     * different length.
102     * <p>
103     * This statistic can be used to perform a Mann-Whitney U test evaluating
104     * the null hypothesis that the two independent samples has equal mean.
105     * </p>
106     * <p>
107     * Let X<sub>i</sub> denote the i'th individual of the first sample and
108     * Y<sub>j</sub> the j'th individual in the second sample. Note that the
109     * samples would often have different length.
110     * </p>
111     * <p>
112     * <strong>Preconditions</strong>:
113     * <ul>
114     * <li>All observations in the two samples are independent.</li>
115     * <li>The observations are at least ordinal (continuous are also ordinal).</li>
116     * </ul>
117     * </p>
118     *
119     * @param x the first sample
120     * @param y the second sample
121     * @return Mann-Whitney U statistic (maximum of U<sup>x</sup> and U<sup>y</sup>)
122     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
123     * @throws NoDataException if {@code x} or {@code y} are zero-length.
124     */
125    public double mannWhitneyU(final double[] x, final double[] y)
126        throws NullArgumentException, NoDataException {
127
128        ensureDataConformance(x, y);
129
130        final double[] z = concatenateSamples(x, y);
131        final double[] ranks = naturalRanking.rank(z);
132
133        double sumRankX = 0;
134
135        /*
136         * The ranks for x is in the first x.length entries in ranks because x
137         * is in the first x.length entries in z
138         */
139        for (int i = 0; i < x.length; ++i) {
140            sumRankX += ranks[i];
141        }
142
143        /*
144         * U1 = R1 - (n1 * (n1 + 1)) / 2 where R1 is sum of ranks for sample 1,
145         * e.g. x, n1 is the number of observations in sample 1.
146         */
147        final double U1 = sumRankX - ((long) x.length * (x.length + 1)) / 2;
148
149        /*
150         * It can be shown that U1 + U2 = n1 * n2
151         */
152        final double U2 = (long) x.length * y.length - U1;
153
154        return FastMath.max(U1, U2);
155    }
156
157    /**
158     * @param Umin smallest Mann-Whitney U value
159     * @param n1 number of subjects in first sample
160     * @param n2 number of subjects in second sample
161     * @return two-sided asymptotic p-value
162     * @throws ConvergenceException if the p-value can not be computed
163     * due to a convergence error
164     * @throws MaxCountExceededException if the maximum number of
165     * iterations is exceeded
166     */
167    private double calculateAsymptoticPValue(final double Umin,
168                                             final int n1,
169                                             final int n2)
170        throws ConvergenceException, MaxCountExceededException {
171
172        /* long multiplication to avoid overflow (double not used due to efficiency
173         * and to avoid precision loss)
174         */
175        final long n1n2prod = (long) n1 * n2;
176
177        // http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U#Normal_approximation
178        final double EU = n1n2prod / 2.0;
179        final double VarU = n1n2prod * (n1 + n2 + 1) / 12.0;
180
181        final double z = (Umin - EU) / FastMath.sqrt(VarU);
182
183        // No try-catch or advertised exception because args are valid
184        // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
185        final NormalDistribution standardNormal = new NormalDistribution(null, 0, 1);
186
187        return 2 * standardNormal.cumulativeProbability(z);
188    }
189
190    /**
191     * Returns the asymptotic <i>observed significance level</i>, or <a href=
192     * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
193     * p-value</a>, associated with a <a
194     * href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney
195     * U statistic</a> comparing mean for two independent samples.
196     * <p>
197     * Let X<sub>i</sub> denote the i'th individual of the first sample and
198     * Y<sub>j</sub> the j'th individual in the second sample. Note that the
199     * samples would often have different length.
200     * </p>
201     * <p>
202     * <strong>Preconditions</strong>:
203     * <ul>
204     * <li>All observations in the two samples are independent.</li>
205     * <li>The observations are at least ordinal (continuous are also ordinal).</li>
206     * </ul>
207     * </p><p>
208     * Ties give rise to biased variance at the moment. See e.g. <a
209     * href="http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf"
210     * >http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf</a>.</p>
211     *
212     * @param x the first sample
213     * @param y the second sample
214     * @return asymptotic p-value
215     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
216     * @throws NoDataException if {@code x} or {@code y} are zero-length.
217     * @throws ConvergenceException if the p-value can not be computed due to a
218     * convergence error
219     * @throws MaxCountExceededException if the maximum number of iterations
220     * is exceeded
221     */
222    public double mannWhitneyUTest(final double[] x, final double[] y)
223        throws NullArgumentException, NoDataException,
224        ConvergenceException, MaxCountExceededException {
225
226        ensureDataConformance(x, y);
227
228        final double Umax = mannWhitneyU(x, y);
229
230        /*
231         * It can be shown that U1 + U2 = n1 * n2
232         */
233        final double Umin = (long) x.length * y.length - Umax;
234
235        return calculateAsymptoticPValue(Umin, x.length, y.length);
236    }
237
238}