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 }