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 }