Package InversionTest ::
Module StatisticalTests
|
|
1 """
2 General statistical tests that operate on population samples,
3 such as the binomial test, sign test, and Wilcoxon signed rank test.
4 Also, contains relevant CDF, PMF, and PDF functions for related
5 distributions such as the binomial and normal distributions.
6
7 Note: If available, SciPy's optimized versions of binomial testing
8 and normal CDF calculations are utilized by default. These should
9 be marginally faster, as they have hooks into C code.
10
11 Author: Benjamin D. Nye
12 License: Apache License V2.0
13 """
14 import itertools
15 import random
16 from math import pi, erf, exp, floor, lgamma, log, sqrt, factorial
17 from InversionTest.Constants import (GREATER_THAN_HYPOTHESIS, LESS_THAN_HYPOTHESIS,
18 TWO_SIDED_HYPOTHESIS, TEST_HYPOTHESES)
19 from InversionTest.MonteCarloSampler import (MCBMeanDifferenceTest, MCBWilcoxonRankTest,
20 MCMeanDifferenceTest, MCWilcoxonRankTest,
21 meanDifference, convertSamplesToRanks)
22
23
24 try:
25 import scipy.stats
26 import scipy.stats.distributions
27 IS_SCIPY_AVAILABLE = True
28 except ImportError:
29 IS_SCIPY_AVAILABLE = False
30
31 SQRT_OF_TWO = sqrt(2.0)
32
33
34
38
39
40
42 """
43 Exact binomial test, where two-sided test uses a minlike formulation.
44 This two-sided approach was chosen to match frameworks like R and SciPy.
45 For reference, the minlike calculation is the sum of all p(k,n) where p(k,n)<=p(x,n)
46 @param x: Number of successes
47 @type x: int
48 @param n: Number of observations
49 @type n: int
50 @param p: Assumed probability of a success, in [0,1]
51 @type p: float
52 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
53 @type alternative: str
54 @return: Probability of the null hypothesis, in [0,1]
55 @rtype: float
56 """
57 raise NotImplementedError
58
60 """
61 A sign test, which works based on the counts that are greater or less
62 than the compared pairs or null hypothesis mean. Uses the binomial test
63 with p=0.5 to calculate the probability.
64 @param series: The series of values to test
65 @type series: list of float
66 @param series2: A series of comparison pairs, optionally. If None, mu is used instead.
67 @type series2: list of float or None
68 @param mu: A comparison value to compare all values in series against (used only if series2 is None)
69 @type mu: float
70 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
71 @type alternative: str
72 @return: Probability of the null hypothesis, in [0,1]
73 @rtype: float
74 """
75 raise NotImplementedError
76
78 """
79 A Wilcoxon Signed Rank test statistic, for distributions symmetric around the median
80 Returns the statistic value and the pValue
81 @param series: A series of values
82 @type series: list of float
83 @param series2: A second series of values, optionally (if None, mu is used instead)
84 @type series2: list of float or None
85 @param mu: The presumed median for values (used only if y is None)
86 @type mu: float
87 @return: Probability of the null hypothesis, given the alternative
88 @rtype: float
89 """
90 raise NotImplementedError
91
92
93
95 """
96 A generic permutation hypothesis test between two sample populations.
97 Runs all permutations of funct(x',y') where x' and y' are generated
98 from the data points in x and y, then finds where funct(x,y) falls
99 into the generated distribution.
100 NOTE: This is crushingly slow as len(x) + len(y) > 10. A Monte Carlo
101 or partial-coverage permutation test is recommended for larger N.
102 @param x: First sample set of data points
103 @type x: list of float
104 @param y: Second sample set of data points
105 @type y: list of float
106 @param funct: Statistic function, in form funct(x, y)
107 @type funct: callable
108 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
109 @type alternative: str
110 @return: Probability of the null hypothesis, given the alternative
111 @rtype: float
112 """
113 nX = len(x)
114 vObs = funct(x, y)
115 if alternative == GREATER_THAN_HYPOTHESIS:
116 counter = 0
117 for z in itertools.permutations(list(x)+list(y)):
118 if funct(z[:nX], z[nX:]) <= vObs:
119 counter += 1
120 pValue = counter/float(factorial(len(x)+len(y)))
121 elif alternative == LESS_THAN_HYPOTHESIS:
122 counter = 0
123 for z in itertools.permutations(list(x)+list(y)):
124 if funct(z[:nX], z[nX:]) >= vObs:
125 counter += 1
126 pValue = counter/float(factorial(len(x)+len(y)))
127 else:
128 counter = 0
129 tieCounter = 0
130 for z in itertools.permutations(list(x)+list(y)):
131 val = funct(z[:nX], z[nX:])
132 if val < vObs:
133 counter += 1
134 elif val == vObs:
135 tieCounter += 1
136 nSamples = factorial(len(x) + len(y))
137 cValue = counter/float(nSamples)
138 tieValue = tieCounter/float(nSamples)
139 pValue = 2.0*min(cValue+tieValue, 1-cValue)
140 return min(1.0, pValue)
141
144 """
145 Permutation test that tests for differences in the means of two samples
146 (e.g., a two-sample t-like statistic of mean(s1)-mean(s2)).
147 If nToApprox is defined, this sets a cutoff for exact estimation after
148 which a Monte Carlo approximation is used instead.
149 @param x: First sample set of data points
150 @type x: list of float
151 @param y: Second sample set of data points
152 @type y: list of float
153 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
154 @type alternative: str
155 @param pValue: The p-Value for the test to confirm, used for Monte Carlo early termination
156 @type pValue: float
157 @param iterations: The max number of iterations to run for Monte Carlo
158 @type iterations: int
159 @param useStoppingRule: If True, use version of MonteCarlo with an unbiased early stopping rule
160 @type useStoppingRule: bool
161 @param maxExactN: The largest N=len(x)+len(y) to calculate an exact test value.
162 For values higher than this, use a Monte Carlo approximation.
163 @type maxExactN: int
164 @return: Probability of the null hypothesis, given the alternative
165 @rtype: float
166 """
167 if maxExactN is None or (len(x) + len(y) <= maxExactN):
168 return permutationTest(x, y, meanDifference, alternative)
169 elif useStoppingRule:
170 return MCBMeanDifferenceTest(x, y, alternative, iterations, pValue).getEstimate()
171 else:
172 return MCMeanDifferenceTest(x, y, alternative, iterations).getEstimate()
173
176 """
177 Permutation test that tests for differences in the ranks of two samples.
178 @param x: First sample set of data points
179 @type x: list of float
180 @param y: Second sample set of data points
181 @type y: list of float
182 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
183 @type alternative: str
184 @param pValue: The p-Value for the test to confirm, used for Monte Carlo early termination
185 @type pValue: float
186 @param iterations: The max number of iterations to run for Monte Carlo
187 @type iterations: int
188 @param useStoppingRule: If True, use version of MonteCarlo with an unbiased early stopping rule
189 @type useStoppingRule: bool
190 @param maxExactN: The largest N=len(x)+len(y) to calculate an exact test value.
191 For values higher than this, use a Monte Carlo approximation.
192 @type maxExactN: int
193 @return: Probability of the null hypothesis, given the alternative
194 @rtype: float
195 """
196 if maxExactN is None or (len(x) + len(y) <= maxExactN):
197 x, y = convertSamplesToRanks(x, y)
198 def simpleWStat(x, y):
199 return sum(y)
200 return permutationTest(x, y, simpleWStat, alternative)
201 elif useStoppingRule:
202 return MCBWilcoxonRankTest(x, y, alternative, iterations, pValue).getEstimate()
203 else:
204 return MCWilcoxonRankTest(x, y, alternative, iterations).getEstimate()
205
206
208 """
209 A Wilcoxon Signed Rank test, for distributions symmetric around the median
210 @param x: A series of values
211 @type x: list of float
212 @param y: A second series of values, optionally (if None, mu is used instead)
213 @type y: list of float or None
214 @param mu: The presumed median for values (used only if y is None)
215 @type mu: float
216 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
217 @type alternative: str
218 @return: Probability of the null hypothesis, given the alternative
219 @rtype: float
220 """
221 wStat, pValue = wilcoxonSignedRankStatistic(x, y, mu)
222 if alternative != TWO_SIDED_HYPOTHESIS:
223 mean = wilcoxonMeanScore(x, y, mu)
224 transform = transformSymmetricPValueHypothesis
225 pValue = transform(mean-wStat, pValue, TWO_SIDED_HYPOTHESIS, alternative)
226 return pValue
227
228
229
231 """
232 Exact binomial test, where two-sided test uses a minlike formulation.
233 This two-sided approach was chosen to match frameworks like R and SciPy.
234 For reference, the minlike calculation is the sum of all p(k,n) where p(k,n)<=p(x,n)
235 @param x: Number of successes
236 @type x: int
237 @param n: Number of observations
238 @type n: int
239 @param p: Assumed probability of a success, in [0,1]
240 @type p: float
241 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
242 @type alternative: str
243 @param useMinlike: If True, calculate a minlike two-tail. Else, return 2*min(p_low, p_high).
244 @type useMinlike: bool
245 @return: Probability of the null hypothesis, in [0,1]
246 @rtype: float
247 """
248 x = int(floor(x))
249 n = int(floor(n))
250
251 if alternative == GREATER_THAN_HYPOTHESIS:
252 if x == 0:
253 pValue = 1.0
254 else:
255 pValue = 1.0-binomialCDF(max(0,x-1), n, p)
256 elif alternative == LESS_THAN_HYPOTHESIS:
257 pValue = binomialCDF(x, n, p)
258
259 elif x == n*p:
260 pValue = 1.0
261 elif useMinlike:
262
263 if x > n*p:
264 x = n-x
265 p = 1-p
266 tail1 = binomialCDF(x, n, p)
267 tail2 = 0
268 pStop = binomialPMF(x, n, p)
269 for k in xrange(0, int(floor(n*p))):
270 pmfVal = binomialPMF(n-k, n, p)
271 if pmfVal > pStop:
272 break
273 else:
274 tail2 += pmfVal
275 pValue = tail1 + tail2
276 else:
277
278 pValue = 2.0*min(pythonBinomialTest(x,n,p, LESS_THAN_HYPOTHESIS),
279 pythonBinomialTest(x,n,p, GREATER_THAN_HYPOTHESIS))
280 return max(0.0, min(1.0, pValue))
281
283 """
284 A sign test, which works based on the counts that are greater or less
285 than the compared pairs or null hypothesis mean. Uses the binomial test
286 with p=0.5 to calculate the probability.
287 @param series: The series of values to test
288 @type series: list of float
289 @param series2: A series of comparison pairs, optionally. If None, mu is used instead.
290 @type series2: list of float or None
291 @param mu: A comparison value to compare all values in series against (used only if series2 is None)
292 @type mu: float
293 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
294 @type alternative: str
295 @return: Probability of the null hypothesis, in [0,1]
296 @rtype: float
297 """
298 if series2:
299 if len(series) != len(series2):
300 raise SignTestInvalidPairLengthError("Sign test requires series of equal length")
301 signCount = [series2[i] < series[i] for i in xrange(len(series)) if series[i] != series2[i]]
302 else:
303 signCount = [mu < series[i] for i in xrange(len(series)) if series[i] != mu]
304 n = len(signCount)
305 x = sum(signCount)
306 if n == 0:
307 raise SignTestCancelationError("All terms in sign test canceled out, test invalid.")
308 pValue = binomialTest(x, n, 0.5, alternative)
309 return pValue
310
312 """
313 A Wilcoxon two-sided test. This uses a normal approximation and
314 adjusts for ties using a variance penalty of (t^3-t)/48.0 for each
315 tie of length t.
316 @param series: A series of values
317 @type series: list of float
318 @param series2: A second series of values, optionally (if None, mu is used instead)
319 @type series2: list of float or None
320 @param mu: The presumed median for values (used only if y is None)
321 @type mu: float
322 @return: Wilcoxon statistic value (W+), Probability of the null hypothesis
323 @rtype: float, float
324 """
325
326 if series2:
327 if len(series) != len(series2):
328 raise IndexError("Wilcoxon Signed Rank test requires series of equal length")
329 diffs = [series2[i]-series[i] for i in xrange(len(series))
330 if series2[i]-series[i] != 0]
331 else:
332 diffs = [mu-series[i] for i in xrange(len(series))
333 if mu-series[i] != 0]
334 N = len(diffs)
335 ties = []
336 wRankSum = 0
337 rank = 1
338
339 diffs.sort(key=absAndValKey)
340 for val, group in itertools.groupby(diffs, abs):
341 group = [x for x in group]
342 numItems = len(group)
343 if numItems > 1:
344 ties.append(numItems)
345 value = (rank + (numItems-1)/2.0)
346 wRankSum += value*len([x for x in group if x < 0])
347 rank += numItems
348
349 mean = N*(N+1)/4.0
350 variance = N*(N+1)*(2*N+1)/24.0 - sum([(t**3-t)/48.0 for t in ties])
351 minWRank = min(wRankSum, N*(N+1)/2.0 - wRankSum)
352 zScore = (minWRank-mean)/sqrt(variance)
353 cdfValue = normalCDFFunction(-abs(zScore), 0.0, 1.0)
354 pValue = 2.0*cdfValue
355 return wRankSum, pValue
356
358 """
359 Get the mean Wilcoxon score, given equality
360 @param series: A series of values
361 @type series: list of float
362 @param series2: A second series of values, optionally (if None, mu is used instead)
363 @type series2: list of float or None
364 @param mu: The presumed median for values (used only if y is None)
365 @type mu: float
366 @return: Mean Wilcoxon statistic value
367 @rtype: float, float
368 """
369 if series2:
370 if len(series) != len(series2):
371 raise IndexError("Wilcoxon Signed Rank test requires series of equal length")
372 diffs = [series2[i]-series[i] for i in xrange(len(series))
373 if series2[i]-series[i] != 0]
374 else:
375 diffs = [mu-series[i] for i in xrange(len(series))
376 if mu-series[i] != 0]
377 N = len(diffs)
378 return N*(N+1)/4.0
379
380
381
383 """
384 Wrapper for the SciPy binomial two-sided test and binomial CDF calculations
385 for one-tailed tests. This allows testing two-sided and both one-sided hypotheses.
386 @param x: Number of successes
387 @type x: int
388 @param n: Number of observations
389 @type n: int
390 @param p: Assumed probability of a success, in [0,1]
391 @type p: float
392 @param alternative: The alternate hypothesis for this test, from TEST_HYPOTHESES set
393 @type alternative: str
394 @return: Probability of the null hypothesis, in [0,1]
395 @rtype: float
396 """
397 if alternative == LESS_THAN_HYPOTHESIS:
398 pValue = 1.0-scipy.stats.binom.sf(x, n, p)
399 elif alternative == GREATER_THAN_HYPOTHESIS:
400 pValue = 1.0-scipy.stats.binom.cdf(max(0,x-1), n, p)
401 else:
402 pValue = scipy.stats.binom_test(x, n, p)
403 return pValue
404
406 """
407 The SciPy Wilcoxon two-sided test. This test always uses
408 a normal approximation and does not adjust for ties properly,
409 but should be (theoretically) faster than the pure python
410 versions here, as it uses C code under the hood.
411 @param series: A series of values
412 @type series: list of float
413 @param series2: A second series of values, optionally (if None, mu is used instead)
414 @type series2: list of float or None
415 @param mu: The presumed median for values (used only if y is None)
416 @type mu: float
417 @return: Wilcoxon statistic value of min(W+,W-), Probability of the null hypothesis
418 @rtype: float, float
419 """
420 if series2 is None:
421 series2 = [mu]*len(series)
422 wStatistic, pValue = scipy.stats.wilcoxon(series, series2)
423 return wStatistic, pValue
424
425
426
428 """
429 Binomial PMF, using log-gamma implementation
430 @param x: # of successes
431 @type x: int
432 @param n: Number of observations
433 @type n: int
434 @param p: Probability of a success
435 @type p: float
436 @return: Point mass for x in the binomial distribution
437 @rtype: float
438 """
439 return exp(lgamma(n+1) - lgamma(x+1) - lgamma(n-x+1) + x*log(p) + (n-x)*log(1-p))
440
442 """
443 Binomial CDF, using log-gamma implementation
444 @param x: # of successes
445 @type x: int
446 @param n: Number of observations
447 @type n: int
448 @param p: Probability of a success
449 @type p: float
450 @return: Cummulative distribution function for x in the binomial distribution
451 @rtype: float
452 """
453 lnMaxFact = lgamma(n+1)
454 q = 1.0-p
455 pValue = 0.0
456 if x <= n*p:
457 for k in xrange(0, int(x+1)):
458 pValue += exp(lnMaxFact - lgamma(k+1) - lgamma(n-k+1)
459 + k*log(p) + (n-k)*log(q))
460 else:
461 for k in xrange(int(x+1), n+1):
462 pValue += exp(lnMaxFact - lgamma(k+1) - lgamma(n-k+1)
463 + k*log(p) + (n-k)*log(q))
464 pValue = 1.0 - pValue
465 return pValue
466
468 """
469 A heuristic for when the normal approximation is reasonable
470 @param n: Number of observations
471 @type n: int
472 @param p: Probability of a success
473 @type p: float
474 @param threshold: Threshold for the heuristic, in [0, 1] where 0 is never and 1 is always
475 @type threshold: float
476 @param minN: Minimum n to have before using the approximation under any circumstances
477 @type minN: int
478 @return: True (use normal distribution) if heuristic < threshold and n > minN, else False
479 @rtype: bool
480 """
481 if n > minN:
482 heuristic = abs((sqrt(p/(1-p)) - sqrt((1-p)/p))/sqrt(n))
483 return heuristic < threshold
484 else:
485 return False
486
488 """
489 Normal CDF (Phi) implementation based on the error function in Python 2.7
490 @param x: Value for CDF, x in P(X<=x) where X ~ N(mu=loc, var=scale**2)
491 @type x: float
492 @param loc: Location parameter (mean)
493 @type loc: float
494 @param scale: Scale parameter (variance)
495 @type scale: float
496 @return: CDF value from the normal distribution
497 @rtype: float
498 """
499 x = (x-loc)/float(scale)
500 return 0.5 + 0.5*erf(x/SQRT_OF_TWO)
501
502
503
505 """
506 Return the absolute value and the value of a number
507 @param x: Some number
508 @type x: int or float
509 @return: abs(x), x
510 @rtype: tuple of (int or float, int or float)
511 """
512 return abs(x), x
513
515 """
516 Generic symmetric paired test, to convert between different-sided results
517 @param statistic: A statistic function, in the form f(x, y, mu)
518 @type statistic: callable
519 @param statisticAlternative: The hypothesis for the test alternative, from TEST_HYPOTHESES set
520 @type statisticAlternative: str
521 @param args: Variable arguments, to pass to the statistic
522 @type args: list of object
523 @param kwds: Variable keyword arguments, to pass to the statistic
524 @type kwds: dict of {str : object}
525 """
526 statisticValue, pValue = statistic(*args, **kwds)
527 pValue = transformSymmetricPValueHypothesis(statisticValue, pValue, statisticAlternative, alternative)
528 return pValue
529
558
559
560
562 """
563 Random number generator that randomly re-positions a standard MT
564 generator based on the system random entropy generator, so that
565 any permutation is theoretically possible.
566 """
567 RESET_COUNT = 2**1000
568 MAX_SKIP = 2**100
570 """ Initialize the random number generators and reset counter"""
571 self.count = 0
572 self.rng = random.Random()
573 self.sysRandom = random.SystemRandom()
574 self.random = self.rng.random
575
577 """ Generate a random number """
578 if self.count >= self.RESET_COUNT:
579 self.count = 0
580 self.rng.jumpahead(self.sysRandom.randrange(0,self.MAX_SKIP))
581 self.count += 1
582 return self.random()
583
584
585
586 signTest = pythonSignTestStatistic
587 wilcoxonSignedRankStatistic = pythonWilcoxonSignedRankStatistic
588 if IS_SCIPY_AVAILABLE:
589
590 binomialTest = scipyBinomialTestStatistic
591 normalCDFFunction = scipy.stats.norm.cdf
592 erfInverse = scipy.special.erfinv
593 else:
594 binomialTest = pythonBinomialTest
595 normalCDFFunction = pythonNormalCDF
596 erfInverse = pythonErfInverse
597