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.DimensionMismatchException;
022 import org.apache.commons.math3.exception.MaxCountExceededException;
023 import org.apache.commons.math3.exception.NoDataException;
024 import org.apache.commons.math3.exception.NullArgumentException;
025 import org.apache.commons.math3.exception.NumberIsTooLargeException;
026 import org.apache.commons.math3.stat.ranking.NaNStrategy;
027 import org.apache.commons.math3.stat.ranking.NaturalRanking;
028 import org.apache.commons.math3.stat.ranking.TiesStrategy;
029 import org.apache.commons.math3.util.FastMath;
030
031 /**
032 * An implementation of the Wilcoxon signed-rank test.
033 *
034 * @version $Id: WilcoxonSignedRankTest.java 1416643 2012-12-03 19:37:14Z tn $
035 */
036 public class WilcoxonSignedRankTest {
037
038 /** Ranking algorithm. */
039 private NaturalRanking naturalRanking;
040
041 /**
042 * Create a test instance where NaN's are left in place and ties get
043 * the average of applicable ranks. Use this unless you are very sure
044 * of what you are doing.
045 */
046 public WilcoxonSignedRankTest() {
047 naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
048 TiesStrategy.AVERAGE);
049 }
050
051 /**
052 * Create a test instance using the given strategies for NaN's and ties.
053 * Only use this if you are sure of what you are doing.
054 *
055 * @param nanStrategy
056 * specifies the strategy that should be used for Double.NaN's
057 * @param tiesStrategy
058 * specifies the strategy that should be used for ties
059 */
060 public WilcoxonSignedRankTest(final NaNStrategy nanStrategy,
061 final TiesStrategy tiesStrategy) {
062 naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
063 }
064
065 /**
066 * Ensures that the provided arrays fulfills the assumptions.
067 *
068 * @param x first sample
069 * @param y second sample
070 * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
071 * @throws NoDataException if {@code x} or {@code y} are zero-length.
072 * @throws DimensionMismatchException if {@code x} and {@code y} do not
073 * have the same length.
074 */
075 private void ensureDataConformance(final double[] x, final double[] y)
076 throws NullArgumentException, NoDataException, DimensionMismatchException {
077
078 if (x == null ||
079 y == null) {
080 throw new NullArgumentException();
081 }
082 if (x.length == 0 ||
083 y.length == 0) {
084 throw new NoDataException();
085 }
086 if (y.length != x.length) {
087 throw new DimensionMismatchException(y.length, x.length);
088 }
089 }
090
091 /**
092 * Calculates y[i] - x[i] for all i
093 *
094 * @param x first sample
095 * @param y second sample
096 * @return z = y - x
097 */
098 private double[] calculateDifferences(final double[] x, final double[] y) {
099
100 final double[] z = new double[x.length];
101
102 for (int i = 0; i < x.length; ++i) {
103 z[i] = y[i] - x[i];
104 }
105
106 return z;
107 }
108
109 /**
110 * Calculates |z[i]| for all i
111 *
112 * @param z sample
113 * @return |z|
114 * @throws NullArgumentException if {@code z} is {@code null}
115 * @throws NoDataException if {@code z} is zero-length.
116 */
117 private double[] calculateAbsoluteDifferences(final double[] z)
118 throws NullArgumentException, NoDataException {
119
120 if (z == null) {
121 throw new NullArgumentException();
122 }
123
124 if (z.length == 0) {
125 throw new NoDataException();
126 }
127
128 final double[] zAbs = new double[z.length];
129
130 for (int i = 0; i < z.length; ++i) {
131 zAbs[i] = FastMath.abs(z[i]);
132 }
133
134 return zAbs;
135 }
136
137 /**
138 * Computes the <a
139 * href="http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test">
140 * Wilcoxon signed ranked statistic</a> comparing mean for two related
141 * samples or repeated measurements on a single sample.
142 * <p>
143 * This statistic can be used to perform a Wilcoxon signed ranked test
144 * evaluating the null hypothesis that the two related samples or repeated
145 * measurements on a single sample has equal mean.
146 * </p>
147 * <p>
148 * Let X<sub>i</sub> denote the i'th individual of the first sample and
149 * Y<sub>i</sub> the related i'th individual in the second sample. Let
150 * Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
151 * </p>
152 * <p>
153 * <strong>Preconditions</strong>:
154 * <ul>
155 * <li>The differences Z<sub>i</sub> must be independent.</li>
156 * <li>Each Z<sub>i</sub> comes from a continuous population (they must be
157 * identical) and is symmetric about a common median.</li>
158 * <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
159 * ordered, so the comparisons greater than, less than, and equal to are
160 * meaningful.</li>
161 * </ul>
162 * </p>
163 *
164 * @param x the first sample
165 * @param y the second sample
166 * @return wilcoxonSignedRank statistic (the larger of W+ and W-)
167 * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
168 * @throws NoDataException if {@code x} or {@code y} are zero-length.
169 * @throws DimensionMismatchException if {@code x} and {@code y} do not
170 * have the same length.
171 */
172 public double wilcoxonSignedRank(final double[] x, final double[] y)
173 throws NullArgumentException, NoDataException, DimensionMismatchException {
174
175 ensureDataConformance(x, y);
176
177 // throws IllegalArgumentException if x and y are not correctly
178 // specified
179 final double[] z = calculateDifferences(x, y);
180 final double[] zAbs = calculateAbsoluteDifferences(z);
181
182 final double[] ranks = naturalRanking.rank(zAbs);
183
184 double Wplus = 0;
185
186 for (int i = 0; i < z.length; ++i) {
187 if (z[i] > 0) {
188 Wplus += ranks[i];
189 }
190 }
191
192 final int N = x.length;
193 final double Wminus = (((double) (N * (N + 1))) / 2.0) - Wplus;
194
195 return FastMath.max(Wplus, Wminus);
196 }
197
198 /**
199 * Algorithm inspired by
200 * http://www.fon.hum.uva.nl/Service/Statistics/Signed_Rank_Algorihms.html#C
201 * by Rob van Son, Institute of Phonetic Sciences & IFOTT,
202 * University of Amsterdam
203 *
204 * @param Wmax largest Wilcoxon signed rank value
205 * @param N number of subjects (corresponding to x.length)
206 * @return two-sided exact p-value
207 */
208 private double calculateExactPValue(final double Wmax, final int N) {
209
210 // Total number of outcomes (equal to 2^N but a lot faster)
211 final int m = 1 << N;
212
213 int largerRankSums = 0;
214
215 for (int i = 0; i < m; ++i) {
216 int rankSum = 0;
217
218 // Generate all possible rank sums
219 for (int j = 0; j < N; ++j) {
220
221 // (i >> j) & 1 extract i's j-th bit from the right
222 if (((i >> j) & 1) == 1) {
223 rankSum += j + 1;
224 }
225 }
226
227 if (rankSum >= Wmax) {
228 ++largerRankSums;
229 }
230 }
231
232 /*
233 * largerRankSums / m gives the one-sided p-value, so it's multiplied
234 * with 2 to get the two-sided p-value
235 */
236 return 2 * ((double) largerRankSums) / ((double) m);
237 }
238
239 /**
240 * @param Wmin smallest Wilcoxon signed rank value
241 * @param N number of subjects (corresponding to x.length)
242 * @return two-sided asymptotic p-value
243 */
244 private double calculateAsymptoticPValue(final double Wmin, final int N) {
245
246 final double ES = (double) (N * (N + 1)) / 4.0;
247
248 /* Same as (but saves computations):
249 * final double VarW = ((double) (N * (N + 1) * (2*N + 1))) / 24;
250 */
251 final double VarS = ES * ((double) (2 * N + 1) / 6.0);
252
253 // - 0.5 is a continuity correction
254 final double z = (Wmin - ES - 0.5) / FastMath.sqrt(VarS);
255
256 // No try-catch or advertised exception because args are valid
257 final NormalDistribution standardNormal = new NormalDistribution(0, 1);
258
259 return 2*standardNormal.cumulativeProbability(z);
260 }
261
262 /**
263 * Returns the <i>observed significance level</i>, or <a href=
264 * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
265 * p-value</a>, associated with a <a
266 * href="http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test">
267 * Wilcoxon signed ranked statistic</a> comparing mean for two related
268 * samples or repeated measurements on a single sample.
269 * <p>
270 * Let X<sub>i</sub> denote the i'th individual of the first sample and
271 * Y<sub>i</sub> the related i'th individual in the second sample. Let
272 * Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
273 * </p>
274 * <p>
275 * <strong>Preconditions</strong>:
276 * <ul>
277 * <li>The differences Z<sub>i</sub> must be independent.</li>
278 * <li>Each Z<sub>i</sub> comes from a continuous population (they must be
279 * identical) and is symmetric about a common median.</li>
280 * <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
281 * ordered, so the comparisons greater than, less than, and equal to are
282 * meaningful.</li>
283 * </ul>
284 * </p>
285 *
286 * @param x the first sample
287 * @param y the second sample
288 * @param exactPValue
289 * if the exact p-value is wanted (only works for x.length <= 30,
290 * if true and x.length > 30, this is ignored because
291 * calculations may take too long)
292 * @return p-value
293 * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
294 * @throws NoDataException if {@code x} or {@code y} are zero-length.
295 * @throws DimensionMismatchException if {@code x} and {@code y} do not
296 * have the same length.
297 * @throws NumberIsTooLargeException if {@code exactPValue} is {@code true}
298 * and {@code x.length} > 30
299 * @throws ConvergenceException if the p-value can not be computed due to
300 * a convergence error
301 * @throws MaxCountExceededException if the maximum number of iterations
302 * is exceeded
303 */
304 public double wilcoxonSignedRankTest(final double[] x, final double[] y,
305 final boolean exactPValue)
306 throws NullArgumentException, NoDataException, DimensionMismatchException,
307 NumberIsTooLargeException, ConvergenceException, MaxCountExceededException {
308
309 ensureDataConformance(x, y);
310
311 final int N = x.length;
312 final double Wmax = wilcoxonSignedRank(x, y);
313
314 if (exactPValue && N > 30) {
315 throw new NumberIsTooLargeException(N, 30, true);
316 }
317
318 if (exactPValue) {
319 return calculateExactPValue(Wmax, N);
320 } else {
321 final double Wmin = ( (double)(N*(N+1)) / 2.0 ) - Wmax;
322 return calculateAsymptoticPValue(Wmin, N);
323 }
324 }
325 }