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 }