github.com/jgbaldwinbrown/perf@v0.1.1/pkg/stats/utest.go (about) 1 // Copyright 2015 The Go Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 package stats 6 7 import ( 8 "math" 9 "sort" 10 ) 11 12 // A LocationHypothesis specifies the alternative hypothesis of a 13 // location test such as a t-test or a Mann-Whitney U-test. The 14 // default (zero) value is to test against the alternative hypothesis 15 // that they differ. 16 type LocationHypothesis int 17 18 //go:generate stringer -type LocationHypothesis 19 20 const ( 21 // LocationLess specifies the alternative hypothesis that the 22 // location of the first sample is less than the second. This 23 // is a one-tailed test. 24 LocationLess LocationHypothesis = -1 25 26 // LocationDiffers specifies the alternative hypothesis that 27 // the locations of the two samples are not equal. This is a 28 // two-tailed test. 29 LocationDiffers LocationHypothesis = 0 30 31 // LocationGreater specifies the alternative hypothesis that 32 // the location of the first sample is greater than the 33 // second. This is a one-tailed test. 34 LocationGreater LocationHypothesis = 1 35 ) 36 37 // A MannWhitneyUTestResult is the result of a Mann-Whitney U-test. 38 type MannWhitneyUTestResult struct { 39 // N1 and N2 are the sizes of the input samples. 40 N1, N2 int 41 42 // U is the value of the Mann-Whitney U statistic for this 43 // test, generalized by counting ties as 0.5. 44 // 45 // Given the Cartesian product of the two samples, this is the 46 // number of pairs in which the value from the first sample is 47 // greater than the value of the second, plus 0.5 times the 48 // number of pairs where the values from the two samples are 49 // equal. Hence, U is always an integer multiple of 0.5 (it is 50 // a whole integer if there are no ties) in the range [0, N1*N2]. 51 // 52 // U statistics always come in pairs, depending on which 53 // sample is "first". The mirror U for the other sample can be 54 // calculated as N1*N2 - U. 55 // 56 // There are many equivalent statistics with slightly 57 // different definitions. The Wilcoxon (1945) W statistic 58 // (generalized for ties) is U + (N1(N1+1))/2. It is also 59 // common to use 2U to eliminate the half steps and Smid 60 // (1956) uses N1*N2 - 2U to additionally center the 61 // distribution. 62 U float64 63 64 // AltHypothesis specifies the alternative hypothesis tested 65 // by this test against the null hypothesis that there is no 66 // difference in the locations of the samples. 67 AltHypothesis LocationHypothesis 68 69 // P is the p-value of the Mann-Whitney test for the given 70 // null hypothesis. 71 P float64 72 } 73 74 // MannWhitneyExactLimit gives the largest sample size for which the 75 // exact U distribution will be used for the Mann-Whitney U-test. 76 // 77 // Using the exact distribution is necessary for small sample sizes 78 // because the distribution is highly irregular. However, computing 79 // the distribution for large sample sizes is both computationally 80 // expensive and unnecessary because it quickly approaches a normal 81 // approximation. Computing the distribution for two 50 value samples 82 // takes a few milliseconds on a 2014 laptop. 83 var MannWhitneyExactLimit = 50 84 85 // MannWhitneyTiesExactLimit gives the largest sample size for which 86 // the exact U distribution will be used for the Mann-Whitney U-test 87 // in the presence of ties. 88 // 89 // Computing this distribution is more expensive than computing the 90 // distribution without ties, so this is set lower. Computing this 91 // distribution for two 25 value samples takes about ten milliseconds 92 // on a 2014 laptop. 93 var MannWhitneyTiesExactLimit = 25 94 95 // MannWhitneyUTest performs a Mann-Whitney U-test [1,2] of the null 96 // hypothesis that two samples come from the same population against 97 // the alternative hypothesis that one sample tends to have larger or 98 // smaller values than the other. 99 // 100 // This is similar to a t-test, but unlike the t-test, the 101 // Mann-Whitney U-test is non-parametric (it does not assume a normal 102 // distribution). It has very slightly lower efficiency than the 103 // t-test on normal distributions. 104 // 105 // Computing the exact U distribution is expensive for large sample 106 // sizes, so this uses a normal approximation for sample sizes larger 107 // than MannWhitneyExactLimit if there are no ties or 108 // MannWhitneyTiesExactLimit if there are ties. This normal 109 // approximation uses both the tie correction and the continuity 110 // correction. 111 // 112 // This can fail with ErrSampleSize if either sample is empty or 113 // ErrSamplesEqual if all sample values are equal. 114 // 115 // This is also known as a Mann-Whitney-Wilcoxon test and is 116 // equivalent to the Wilcoxon rank-sum test, though the Wilcoxon 117 // rank-sum test differs in nomenclature. 118 // 119 // [1] Mann, Henry B.; Whitney, Donald R. (1947). "On a Test of 120 // Whether one of Two Random Variables is Stochastically Larger than 121 // the Other". Annals of Mathematical Statistics 18 (1): 50–60. 122 // 123 // [2] Klotz, J. H. (1966). "The Wilcoxon, Ties, and the Computer". 124 // Journal of the American Statistical Association 61 (315): 772-787. 125 func MannWhitneyUTest(x1, x2 []float64, alt LocationHypothesis) (*MannWhitneyUTestResult, error) { 126 n1, n2 := len(x1), len(x2) 127 if n1 == 0 || n2 == 0 { 128 return nil, ErrSampleSize 129 } 130 131 // Compute the U statistic and tie vector T. 132 x1 = append([]float64(nil), x1...) 133 x2 = append([]float64(nil), x2...) 134 sort.Float64s(x1) 135 sort.Float64s(x2) 136 merged, labels := labeledMerge(x1, x2) 137 138 R1 := 0.0 139 T, hasTies := []int{}, false 140 for i := 0; i < len(merged); { 141 rank1, nx1, v1 := i+1, 0, merged[i] 142 // Consume samples that tie this sample (including itself). 143 for ; i < len(merged) && merged[i] == v1; i++ { 144 if labels[i] == 1 { 145 nx1++ 146 } 147 } 148 // Assign all tied samples the average rank of the 149 // samples, where merged[0] has rank 1. 150 if nx1 != 0 { 151 rank := float64(i+rank1) / 2 152 R1 += rank * float64(nx1) 153 } 154 T = append(T, i-rank1+1) 155 if i > rank1 { 156 hasTies = true 157 } 158 } 159 U1 := R1 - float64(n1*(n1+1))/2 160 161 // Compute the smaller of U1 and U2 162 U2 := float64(n1*n2) - U1 163 Usmall := math.Min(U1, U2) 164 165 var p float64 166 if !hasTies && n1 <= MannWhitneyExactLimit && n2 <= MannWhitneyExactLimit || 167 hasTies && n1 <= MannWhitneyTiesExactLimit && n2 <= MannWhitneyTiesExactLimit { 168 // Use exact U distribution. U1 will be an integer. 169 if len(T) == 1 { 170 // All values are equal. Test is meaningless. 171 return nil, ErrSamplesEqual 172 } 173 174 dist := UDist{N1: n1, N2: n2, T: T} 175 switch alt { 176 case LocationDiffers: 177 if U1 == U2 { 178 // The distribution is symmetric about 179 // Usmall. Since the distribution is 180 // discrete, the CDF is discontinuous 181 // and if simply double CDF(Usmall), 182 // we'll double count the 183 // (non-infinitesimal) probability 184 // mass at Usmall. What we want is 185 // just the integral of the whole CDF, 186 // which is 1. 187 p = 1 188 } else { 189 p = dist.CDF(Usmall) * 2 190 } 191 192 case LocationLess: 193 p = dist.CDF(U1) 194 195 case LocationGreater: 196 p = 1 - dist.CDF(U1-1) 197 } 198 } else { 199 // Use normal approximation (with tie and continuity 200 // correction). 201 t := tieCorrection(T) 202 N := float64(n1 + n2) 203 μ_U := float64(n1*n2) / 2 204 σ_U := math.Sqrt(float64(n1*n2) * ((N + 1) - t/(N*(N-1))) / 12) 205 if σ_U == 0 { 206 return nil, ErrSamplesEqual 207 } 208 numer := U1 - μ_U 209 // Perform continuity correction. 210 switch alt { 211 case LocationDiffers: 212 numer -= mathSign(numer) * 0.5 213 case LocationLess: 214 numer += 0.5 215 case LocationGreater: 216 numer -= 0.5 217 } 218 z := numer / σ_U 219 switch alt { 220 case LocationDiffers: 221 p = 2 * math.Min(StdNormal.CDF(z), 1-StdNormal.CDF(z)) 222 case LocationLess: 223 p = StdNormal.CDF(z) 224 case LocationGreater: 225 p = 1 - StdNormal.CDF(z) 226 } 227 } 228 229 return &MannWhitneyUTestResult{N1: n1, N2: n2, U: U1, 230 AltHypothesis: alt, P: p}, nil 231 } 232 233 // labeledMerge merges sorted lists x1 and x2 into sorted list merged. 234 // labels[i] is 1 or 2 depending on whether merged[i] is a value from 235 // x1 or x2, respectively. 236 func labeledMerge(x1, x2 []float64) (merged []float64, labels []byte) { 237 merged = make([]float64, len(x1)+len(x2)) 238 labels = make([]byte, len(x1)+len(x2)) 239 240 i, j, o := 0, 0, 0 241 for i < len(x1) && j < len(x2) { 242 if x1[i] < x2[j] { 243 merged[o] = x1[i] 244 labels[o] = 1 245 i++ 246 } else { 247 merged[o] = x2[j] 248 labels[o] = 2 249 j++ 250 } 251 o++ 252 } 253 for ; i < len(x1); i++ { 254 merged[o] = x1[i] 255 labels[o] = 1 256 o++ 257 } 258 for ; j < len(x2); j++ { 259 merged[o] = x2[j] 260 labels[o] = 2 261 o++ 262 } 263 return 264 } 265 266 // tieCorrection computes the tie correction factor Σ_j (t_j³ - t_j) 267 // where t_j is the number of ties in the j'th rank. 268 func tieCorrection(ties []int) float64 { 269 t := 0 270 for _, tie := range ties { 271 t += tie*tie*tie - tie 272 } 273 return float64(t) 274 }