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