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