The OpenD Programming Language

1 /**Hypothesis testing beyond simple CDFs.  All functions work with input
2  * ranges with elements implicitly convertible to double unless otherwise noted.
3  *
4  * Author:  David Simcha*/
5  /*
6  * License:
7  * Boost Software License - Version 1.0 - August 17th, 2003
8  *
9  * Permission is hereby granted, free of charge, to any person or organization
10  * obtaining a copy of the software and accompanying documentation covered by
11  * this license (the "Software") to use, reproduce, display, distribute,
12  * execute, and transmit the Software, and to prepare derivative works of the
13  * Software, and to permit third-parties to whom the Software is furnished to
14  * do so, all subject to the following:
15  *
16  * The copyright notices in the Software and this entire statement, including
17  * the above license grant, this restriction and the following disclaimer,
18  * must be included in all copies of the Software, in whole or in part, and
19  * all derivative works of the Software, unless such copies or derivative
20  * works are solely in the form of machine-executable object code generated by
21  * a source language processor.
22  *
23  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
24  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25  * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
26  * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
27  * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
28  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29  * DEALINGS IN THE SOFTWARE.
30  */
31 module dstats.tests;
32 
33 import std.functional, std.range, std.conv, std.math, std.traits,
34        std.exception, std.typetuple;
35 
36 import std.algorithm : reverse, copy, max, min, swap, map, reduce;
37 
38 import dstats.base, dstats.distrib, dstats.alloc, dstats.summary, dstats.sort;
39 
40 static import dstats.cor;
41 
42 private static import dstats.infotheory;
43 
44 version(unittest) {
45     import std.stdio, dstats.random;
46 }
47 
48 /**Alternative hypotheses.  Exact meaning varies with test used.*/
49 enum Alt {
50     /// f(input1) != X
51     twoSided,
52 
53     /// f(input1) < X
54     less,
55 
56     /// f(input1) > X
57     greater,
58 
59     /**
60     Skip P-value computation (and confidence intervals if applicable)
61     and just return the test statistic.
62     */
63     none
64 }
65 
66 /**
67 A plain old data struct for returning the results of hypothesis tests.
68 */
69 struct TestRes {
70     /// The test statistic.  What exactly this is is specific to the test.
71     double testStat;
72 
73     /**The P-value against the provided alternative.  This struct can
74      * be implicitly converted to just the P-value via alias this.*/
75     double p;
76 
77     /// Allow implicit conversion to the P-value.
78     alias p this;
79 
80     ///
81     string toString() {
82         return text("Test Statistic = ", testStat, "\nP = ", p);
83     }
84 }
85 
86 /**
87 A plain old data struct for returning the results of hypothesis tests
88 that also produce confidence intervals.  Contains, can implicitly convert
89 to, a TestRes.
90 */
91 struct ConfInt {
92     ///  This is alias this'd.
93     TestRes testRes;
94 
95     ///  Lower bound of the confidence interval at the level specified.
96     double lowerBound;
97 
98     ///  Upper bound of the confidence interval at the level specified.
99     double upperBound;
100 
101     alias testRes this;
102 
103     ///
104     string toString() {
105         return text("Test Statistic = ", testRes.testStat, "\nP = ", testRes.p,
106                 "\nLower Confidence Bound = ", lowerBound,
107                 "\nUpper Confidence Bound = ", upperBound);
108     }
109 }
110 
111 /**
112 Tests whether a struct/class has the necessary information for calculating
113 a T-test.  It must have a property .mean (mean), .stdev (stdandard deviation),
114 .var (variance), and .N (sample size).
115 */
116 template isSummary(T) {
117     enum bool isSummary = is(typeof(T.init.mean)) && is(typeof(T.init.stdev)) &&
118         is(typeof(T.init.var)) && is(typeof(T.init.N));
119 }
120 
121 /**
122 One-sample Student's T-test for difference between mean of data and
123  a fixed value.  Alternatives are Alt.less, meaning mean(data) < testMean,
124  Alt.greater, meaning mean(data) > testMean, and Alt.twoSided, meaning
125  mean(data)!= testMean.
126 
127  data may be either an iterable with elements implicitly convertible to
128  double or a summary struct (see isSummary).
129  *
130  Examples:
131  ---
132  uint[] data = [1,2,3,4,5];
133 
134  // Test the null hypothesis that the mean of data is >= 1 against the
135  // alternative that the mean of data is < 1.  Calculate confidence
136  // intervals at 90%.
137  auto result1 = studentsTTest(data, 1, Alt.less, 0.9);
138 
139  // Do the same thing, only this time we've already calculated the summary
140  // statistics explicitly before passing them to studensTTest.
141  auto summary = meanStdev(data);
142  writeln(summary.stdev);
143  result2 = studentsTTest(summary, 1, Alt.less, 0.9);  // Same as result1.
144  assert(result1 == result2);
145  ---
146 
147  Returns:  A ConfInt containing T, the P-value and the boundaries of
148  the confidence interval for mean(data) at the level specified.
149 
150  References:  http://en.wikipedia.org/wiki/Student%27s_t-test
151  */
152 ConfInt studentsTTest(T)(
153     T data,
154     double testMean = 0,
155     Alt alt = Alt.twoSided,
156     double confLevel = 0.95
157 ) if( (isSummary!T || doubleIterable!T)) {
158     enforceConfidence(confLevel);
159     dstatsEnforce(isFinite(testMean), "testMean must not be infinite or NaN.");
160 
161     static if(isSummary!T) {
162         return pairedTTest(data, testMean, alt, confLevel);
163     } else static if(doubleIterable!T) {
164         return pairedTTest(meanStdev(data), testMean, alt, confLevel);
165     }
166 }
167 
168 unittest {
169     auto t1 = studentsTTest([1, 2, 3, 4, 5].dup, 2);
170     assert(approxEqual(t1.testStat, 1.4142));
171     assert(approxEqual(t1.p, 0.2302));
172     assert(approxEqual(t1.lowerBound, 1.036757));
173     assert(approxEqual(t1.upperBound, 4.963243));
174     assert(t1 == studentsTTest( meanStdev([1,2,3,4,5].dup), 2));
175 
176     auto t2 = studentsTTest([1, 2, 3, 4, 5].dup, 2, Alt.less);
177     assert(approxEqual(t2.p, .8849));
178     assert(approxEqual(t2.testStat, 1.4142));
179     assert(t2.lowerBound == -double.infinity);
180     assert(approxEqual(t2.upperBound, 4.507443));
181 
182     auto t3 = studentsTTest( summary([1, 2, 3, 4, 5].dup), 2, Alt.greater);
183     assert(approxEqual(t3.p, .1151));
184     assert(approxEqual(t3.testStat, 1.4142));
185     assert(approxEqual(t3.lowerBound, 1.492557));
186     assert(t3.upperBound == double.infinity);
187 }
188 
189 /**
190 Two-sample T test for a difference in means,
191 assumes variances of samples are equal.  Alteratives are Alt.less, meaning
192 mean(sample1) - mean(sample2) < testMean, Alt.greater, meaning
193 mean(sample1) - mean(sample2) > testMean, and Alt.twoSided, meaning
194 mean(sample1) - mean(sample2) != testMean.
195 
196 sample1 and sample2 may be either iterables with elements implicitly
197 convertible to double or summary structs (see isSummary).
198 
199 Returns:  A ConfInt containing the T statistic, the P-value, and the
200 boundaries of the confidence interval for the difference between means
201 of sample1 and sample2 at the specified level.
202 
203 References:  http://en.wikipedia.org/wiki/Student%27s_t-test
204  */
205 ConfInt studentsTTest(T, U)(
206     T sample1,
207     U sample2,
208     double testMean = 0,
209     Alt alt = Alt.twoSided,
210     double confLevel = 0.95
211 ) if( (doubleIterable!T || isSummary!T) && (doubleIterable!U || isSummary!U)) {
212     enforceConfidence(confLevel);
213     dstatsEnforce(isFinite(testMean), "testMean must not be infinite or nan.");
214 
215     static if(isSummary!T) {
216         alias sample1 s1summ;
217     } else {
218         immutable s1summ = meanStdev(sample1);
219     }
220 
221     static if(isSummary!U) {
222         alias sample2 s2summ;
223     } else {
224         immutable s2summ = meanStdev(sample2);
225     }
226 
227     immutable n1 = s1summ.N, n2 = s2summ.N;
228 
229     immutable sx1x2 = sqrt((n1 * s1summ.mse + n2 * s2summ.mse) /
230                  (n1 + n2 - 2));
231     immutable normSd = (sx1x2 * sqrt((1.0 / n1) + (1.0 / n2)));
232     immutable meanDiff = s1summ.mean - s2summ.mean;
233     ConfInt ret;
234     ret.testStat = (meanDiff - testMean) / normSd;
235     if(alt == Alt.none) {
236         return ret;
237     } else if(alt == Alt.less) {
238         ret.p = studentsTCDF(ret.testStat, n1 + n2 - 2);
239 
240         ret.lowerBound = -double.infinity;
241         if(confLevel > 0) {
242             immutable delta = invStudentsTCDF(1 - confLevel, n1 + n2 - 2)
243                 * normSd;
244             ret.upperBound = meanDiff - delta;
245         } else {
246             ret.upperBound = meanDiff;
247         }
248 
249     } else if(alt == Alt.greater) {
250         ret.p = studentsTCDF(-ret.testStat, n1 + n2 - 2);
251 
252         ret.upperBound = double.infinity;
253         if(confLevel > 0) {
254             immutable delta = invStudentsTCDF(1 - confLevel, n1 + n2 - 2)
255                 * normSd;
256             ret.lowerBound = meanDiff + delta;
257         } else {
258             ret.lowerBound = meanDiff;
259         }
260 
261     } else {
262         immutable t = ret.testStat;
263         ret.p = 2 * ((t < 0) ?
264                     studentsTCDF(t, n1 + n2 - 2) :
265                     studentsTCDFR(t, n1 + n2 - 2));
266 
267         if(confLevel > 0) {
268             immutable delta = invStudentsTCDF(
269                 0.5 * (1 - confLevel), n1 + n2 - 2) * normSd;
270             ret.lowerBound = meanDiff + delta;
271             ret.upperBound = meanDiff - delta;
272         } else {
273             ret.lowerBound = meanDiff;
274             ret.upperBound = meanDiff;
275         }
276     }
277     return ret;
278 }
279 
280 unittest {
281     // Values from R.
282     auto t1 = studentsTTest([1,2,3,4,5], [1,3,4,5,7,9]);
283     assert(approxEqual(t1.p, 0.2346));
284     assert(approxEqual(t1.testStat, -1.274));
285     assert(approxEqual(t1.lowerBound, -5.088787));
286     assert(approxEqual(t1.upperBound, 1.422120));
287 
288 
289     assert(approxEqual(studentsTTest([1,2,3,4,5], [1,3,4,5,7,9], 0, Alt.less),
290            0.1173));
291     assert(approxEqual(studentsTTest([1,2,3,4,5], [1,3,4,5,7,9], 0, Alt.greater),
292            0.8827));
293     auto t2 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 5);
294     assert(approxEqual(t2.p, 0.44444));
295     assert(approxEqual(t2.testStat, -0.7998));
296     assert(approxEqual(t2.lowerBound, -0.3595529));
297     assert(approxEqual(t2.upperBound, 7.5595529));
298 
299 
300     auto t5 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 0, Alt.less);
301     assert(approxEqual(t5.p, 0.965));
302     assert(approxEqual(t5.testStat, 2.0567));
303     assert(approxEqual(t5.upperBound, 6.80857));
304     assert(t5.lowerBound == -double.infinity);
305 
306     auto t6 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 0, Alt.greater);
307     assert(approxEqual(t6.p, 0.03492));
308     assert(approxEqual(t6.testStat, 2.0567));
309     assert(approxEqual(t6.lowerBound, 0.391422));
310     assert(t6.upperBound == double.infinity);
311 
312     auto t7 = studentsTTest([1, 2, 4], [3]);
313     assert(approxEqual(t7.p, 0.7418));
314     assert(approxEqual(t7.testStat, 0.-.378));
315     assert(approxEqual(t7.lowerBound, -8.255833));
316     assert(approxEqual(t7.upperBound, 6.922499));
317 
318 }
319 
320 /**
321 Two-sample T-test for difference in means.  Does not assume variances are equal.
322 Alteratives are Alt.less, meaning mean(sample1) - mean(sample2) < testMean,
323 Alt.greater, meaning mean(sample1) - mean(sample2) > testMean, and
324 Alt.twoSided, meaning mean(sample1) - mean(sample2) != testMean.
325 
326 sample1 and sample2 may be either iterables with elements implicitly
327 convertible to double or summary structs (see isSummary).
328 
329 Returns:  A ConfInt containing the T statistic, the P-value, and the
330 boundaries of the confidence interval for the difference between means
331 of sample1 and sample2 at the specified level.
332 
333 References:  http://en.wikipedia.org/wiki/Student%27s_t-test
334  */
335 ConfInt welchTTest(T, U)(T sample1, U sample2, double testMean = 0,
336     Alt alt = Alt.twoSided, double confLevel = 0.95)
337 if( (isSummary!T || doubleIterable!T) && (isSummary!U || doubleIterable!U)) {
338     enforceConfidence(confLevel);
339     dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or NaN.");
340 
341     static if(isSummary!T) {
342         alias sample1 s1summ;
343     } else {
344         auto s1summ = meanStdev(sample1);
345     }
346 
347     static if(isSummary!U) {
348         alias sample2 s2summ;
349     } else {
350         auto s2summ = meanStdev(sample2);
351     }
352 
353     immutable double n1 = s1summ.N,
354                    n2 = s2summ.N;
355 
356     immutable v1 = s1summ.var, v2 = s2summ.var;
357     immutable double sx1x2 = sqrt(v1 / n1 + v2 / n2);
358     immutable double meanDiff = s1summ.mean - s2summ.mean - testMean;
359     immutable double t = meanDiff / sx1x2;
360     double numerator = v1 / n1 + v2 / n2;
361     numerator *= numerator;
362     double denom1 = v1 / n1;
363     denom1 = denom1 * denom1 / (n1 - 1);
364     double denom2 = v2 / n2;
365     denom2 = denom2 * denom2 / (n2 - 1);
366     immutable double df = numerator / (denom1 + denom2);
367 
368     ConfInt ret;
369     ret.testStat = t;
370     if(alt == Alt.none) {
371         return ret;
372     } else if(alt == Alt.less) {
373         ret.p = studentsTCDF(t, df);
374         ret.lowerBound = -double.infinity;
375 
376         if(confLevel > 0) {
377             ret.upperBound = meanDiff +
378                 testMean - invStudentsTCDF(1 - confLevel, df) * sx1x2;
379         } else {
380             ret.upperBound = meanDiff + testMean;
381         }
382 
383     } else if(alt == Alt.greater) {
384         ret.p = studentsTCDF(-t, df);
385         ret.upperBound = double.infinity;
386 
387         if(confLevel > 0) {
388             ret.lowerBound = meanDiff +
389                 testMean + invStudentsTCDF(1 - confLevel, df) * sx1x2;
390         } else {
391             ret.lowerBound = meanDiff + testMean;
392         }
393 
394     } else {
395         ret.p = 2 * ((t < 0) ?
396                      studentsTCDF(t, df) :
397                      studentsTCDF(-t, df));
398 
399         if(confLevel > 0) {
400             double delta = invStudentsTCDF(0.5 * (1 - confLevel), df) * sx1x2;
401             ret.upperBound = meanDiff + testMean - delta;
402             ret.lowerBound = meanDiff + testMean + delta;
403         } else {
404             ret.upperBound = meanDiff + testMean;
405             ret.lowerBound = meanDiff + testMean;
406         }
407     }
408     return ret;
409 }
410 
411 unittest {
412     // Values from R.
413     auto t1 = welchTTest( meanStdev([1,2,3,4,5]), [1,3,4,5,7,9], 2);
414     assert(approxEqual(t1.p, 0.02285));
415     assert(approxEqual(t1.testStat, -2.8099));
416     assert(approxEqual(t1.lowerBound, -4.979316));
417     assert(approxEqual(t1.upperBound, 1.312649));
418 
419     auto t2 = welchTTest([1,2,3,4,5], summary([1,3,4,5,7,9]), -1, Alt.less);
420     assert(approxEqual(t2.p, 0.2791));
421     assert(approxEqual(t2.testStat, -0.6108));
422     assert(t2.lowerBound == -double.infinity);
423     assert(approxEqual(t2.upperBound, 0.7035534));
424 
425     auto t3 = welchTTest([1,2,3,4,5], [1,3,4,5,7,9], 0.5, Alt.greater);
426     assert(approxEqual(t3.p, 0.9372));
427     assert(approxEqual(t3.testStat, -1.7104));
428     assert(approxEqual(t3.lowerBound, -4.37022));
429     assert(t3.upperBound == double.infinity);
430 
431     assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4]).p, 0.06616));
432     assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4], 0,
433         Alt.less).p, 0.967));
434     assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4], 0,
435         Alt.greater).p, 0.03308));
436 }
437 
438 /**
439 Paired T test.  Tests the hypothesis that the mean difference between
440 corresponding elements of before and after is testMean.  Alternatives are
441 Alt.less, meaning the that the true mean difference (before[i] - after[i])
442 is less than testMean, Alt.greater, meaning the true mean difference is
443 greater than testMean, and Alt.twoSided, meaning the true mean difference is not
444 equal to testMean.
445 
446 before and after must be input ranges with elements implicitly convertible
447 to double and must have the same length.
448 
449 Returns:  A ConfInt containing the T statistic, the P-value, and the
450 boundaries of the confidence interval for the mean difference between
451 corresponding elements of sample1 and sample2 at the specified level.
452 
453 References:  http://en.wikipedia.org/wiki/Student%27s_t-test
454 */
455 ConfInt pairedTTest(T, U)(
456     T before,
457     U after,
458     double testMean = 0,
459     Alt alt = Alt.twoSided,
460     double confLevel = 0.95
461 ) if(doubleInput!(T) && doubleInput!(U) && isInputRange!T && isInputRange!U) {
462     enforceConfidence(confLevel);
463     dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");
464 
465     MeanSD msd;
466     while(!before.empty && !after.empty) {
467         immutable diff = cast(double) before.front - cast(double) after.front;
468         before.popFront();
469         after.popFront();
470         msd.put(diff);
471     }
472 
473     dstatsEnforce(before.empty && after.empty,
474         "before and after have different lengths in pairedTTest.");
475 
476     return pairedTTest(msd, testMean, alt, confLevel);
477 }
478 
479 /**
480 Compute a paired T test directly from summary statistics of the differences
481 between corresponding samples.
482 
483 Examples:
484 ---
485 float[] data1 = [8, 6, 7, 5, 3, 0, 9];
486 float[] data2 = [3, 6, 2, 4, 3, 6, 8];
487 
488 // Calculate summary statistics on difference explicitly.
489 MeanSD summary;
490 foreach(i; 0..data1.length) {
491     summary.put(data1[i] - data2[i]);
492 }
493 
494 // Test the null hypothesis that the mean difference between corresponding
495 // elements (data1[i] - data2[i]) is greater than 5 against the null that it
496 // is <= 5.  Calculate confidence intervals at 99%.
497 auto result = pairedTTest(summary, 5, Alt.twoSided, 0.99);
498 
499 // This is equivalent to:
500 auto result2 = pairedTTest(data1, data2, 5, Alt.twoSided, 0.99);
501 ---
502 
503 References:  http://en.wikipedia.org/wiki/Student%27s_t-test
504  */
505 ConfInt pairedTTest(T)(
506     T diffSummary,
507     double testMean = 0,
508     Alt alt = Alt.twoSided,
509     double confLevel = 0.95
510 ) if(isSummary!T) {
511     enforceConfidence(confLevel);
512     dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");
513 
514     if(diffSummary.N < 2) {
515         return ConfInt.init;
516     }
517 
518     // Save typing.
519     alias diffSummary msd;
520 
521     ConfInt ret;
522     ret.testStat = (msd.mean - testMean) / msd.stdev * sqrt(msd.N);
523     auto sampleMean = msd.mean;
524     auto sampleSd = msd.stdev;
525     double normSd = sampleSd / sqrt(msd.N);
526     ret.testStat = (sampleMean - testMean) / normSd;
527 
528     if(alt == Alt.none) {
529         return ret;
530     } else if(alt == Alt.less) {
531         ret.p = studentsTCDF(ret.testStat, msd.N - 1);
532         ret.lowerBound = -double.infinity;
533 
534         if(confLevel > 0) {
535             double delta = invStudentsTCDF(1 - confLevel, msd.N - 1) * normSd;
536             ret.upperBound = sampleMean - delta;
537         } else {
538             ret.upperBound = sampleMean;
539         }
540 
541     } else if(alt == Alt.greater) {
542         ret.p = studentsTCDF(-ret.testStat, msd.N - 1);
543         ret.upperBound = double.infinity;
544 
545         if(confLevel > 0) {
546             double delta = invStudentsTCDF(1 - confLevel, msd.N - 1) * normSd;
547             ret.lowerBound = sampleMean + delta;
548         } else {
549             ret.lowerBound = sampleMean;
550         }
551 
552     } else {
553         immutable double t = ret.testStat;
554         ret.p = 2 * ((t < 0) ?
555                       studentsTCDF(t, msd.N - 1) :
556                       studentsTCDF(-t, msd.N - 1));
557 
558         if(confLevel > 0) {
559             double delta = invStudentsTCDF(0.5 * (1 - confLevel), msd.N - 1) * normSd;
560             ret.lowerBound = sampleMean + delta;
561             ret.upperBound = sampleMean - delta;
562         } else {
563             ret.lowerBound = ret.upperBound = sampleMean;
564         }
565 
566     }
567     return ret;
568 }
569 
570 unittest {
571     // Values from R.
572     auto t1 = pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1);
573     assert(approxEqual(t1.p, 0.02131));
574     assert(approxEqual(t1.testStat, -3.6742));
575     assert(approxEqual(t1.lowerBound, -2.1601748));
576     assert(approxEqual(t1.upperBound, 0.561748));
577 
578     assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.less).p, 0.0889));
579     assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.greater).p, 0.9111));
580     assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.twoSided).p, 0.1778));
581     assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1, Alt.less).p, 0.01066));
582     assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1, Alt.greater).p, 0.9893));
583 }
584 
585 /**
586 Tests the null hypothesis that the variances of all groups are equal against
587 the alternative that heteroscedasticity exists.  data must be either a
588 tuple of ranges or a range of ranges.  central is an alias for the measure
589 of central tendency to be used.  This can be any function that maps a
590 forward range of numeric types to a numeric type.  The commonly used ones
591 are median (default) and mean (less robust).  Trimmed mean is sometimes
592 useful, but is currently not implemented in dstats.summary.
593 
594 References:
595 Levene, Howard (1960). "Robust tests for equality of variances". in Ingram
596 Olkin, Harold Hotelling et al. Contributions to Probability and Statistics:
597 Essays in Honor of Harold Hotelling. Stanford University Press. pp. 278-292.
598 
599 Examples:
600 ---
601 int[] sample1 = [1,2,3,4,5];
602 int[] sample2 = [100,200,300,400,500];
603 auto result = levenesTest(sample1, sample2);
604 
605 // Clearly the variances are different between these two samples.
606 assert( approxEqual(result.testStat, 10.08));
607 assert( approxEqual(result.p, 0.01310));
608 ---
609 */
610 TestRes levenesTest(alias central = median, T...)(T data) {
611     return anovaLevene!(true, false, central, T)(data);
612 }
613 
614 unittest {
615     // Values from R's car package, which uses the median definition
616     // exclusively.
617     auto res1 = levenesTest([1,2,3,4,5][], [2,4,8,16,32][]);
618     assert(approxEqual(res1.testStat, 3.0316));
619     assert(approxEqual(res1.p, 0.1198), res1.toString());
620 
621     auto res2 = levenesTest([[1,2,3,4,5][], [100,200,300,400,500,600][]][]);
622     assert(approxEqual(res2.testStat, 13.586));
623     assert(approxEqual(res2.p, 0.005029));
624 
625     auto res3 = levenesTest([8,6,7,5,3,0,9][], [3,6,2,4,3,6][]);
626     assert(approxEqual(res3.testStat, 1.1406));
627     assert(approxEqual(res3.p, 0.3084));
628 }
629 
630 /**
631 The F-test is a one-way ANOVA extension of the T-test to >2 groups.
632 It's useful when you have 3 or more groups with equal variance and want
633 to test whether their means are equal.  Data can be input as either a
634 tuple or a range.  This may contain any combination of ranges of numeric
635 types, MeanSD structs and Summary structs.
636 
637 Note:  This test makes the assumption that all groups have equal variances,
638 also known as homoskedasticity.  For a similar test that does not make these
639 assumptions, see welchAnova.
640 
641 Examples:
642 ---
643 uint[] thing1 = [3,1,4,1],
644        thing2 = [5,9,2,6,5,3],
645        thing3 = [5,8,9,7,9,3];
646 auto result = fTest(thing1, meanStdev(thing2), summary(thing3));
647 assert(approxEqual(result.testStat, 4.9968));
648 assert(approxEqual(result.p, 0.02456));
649 ---
650 
651 References:  http://en.wikipedia.org/wiki/F-test
652 
653 Returns:
654 
655 A TestRes containing the F statistic and the P-value for the alternative
656 that the means of the groups are different against the null that they
657 are identical.
658 */
659 TestRes fTest(T...)(T data) {
660     return anovaLevene!(false, false, "dummy", T)(data);
661 }
662 
663 /**
664 Same as fTest, except that this test does not require the assumption of
665 equal variances.  In exchange it's slightly less powerful.
666 
667 References:
668 
669 B.L. Welch. On the Comparison of Several Mean Values: An Alternative Approach
670 Biometrika, Vol. 38, No. 3/4 (Dec., 1951), pp. 330-336.
671  */
672 TestRes welchAnova(T...)(T data) {
673     return anovaLevene!(false, true, "dummy", T)(data);
674 }
675 
676 unittest {
677     // Values from R.
678     uint[] thing1 = [3,1,4,1],
679            thing2 = [5,9,2,6,5,3],
680            thing3 = [5,8,9,7,9,3];
681     auto result = fTest(thing1, meanStdev(thing2), summary(thing3));
682     assert(approxEqual(result.testStat, 4.9968));
683     assert(approxEqual(result.p, 0.02456));
684 
685     auto welchRes1 = welchAnova(thing1, thing2, thing3);
686     assert( approxEqual(welchRes1.testStat, 6.7813));
687     assert( approxEqual(welchRes1.p, 0.01706));
688 
689     // Test array case.
690     auto res2 = fTest([thing1, thing2, thing3].dup);
691     assert(approxEqual(result.testStat, res2.testStat));
692     assert(approxEqual(result.p, res2.p));
693 
694     thing1 = [2,7,1,8,2];
695     thing2 = [8,1,8];
696     thing3 = [2,8,4,5,9];
697     auto res3 = fTest(thing1, thing2, thing3);
698     assert(approxEqual(res3.testStat, 0.377));
699     assert(approxEqual(res3.p, 0.6953));
700 
701     auto res4 = fTest([summary(thing1), summary(thing2), summary(thing3)][]);
702     assert(approxEqual(res4.testStat, res3.testStat));
703     assert(approxEqual(res4.testStat, res3.testStat));
704 
705     auto welchRes2 = welchAnova(summary(thing1), thing2, thing3);
706     assert( approxEqual(welchRes2.testStat, 0.342));
707     assert( approxEqual(welchRes2.p, 0.7257));
708 
709     auto res5 = fTest([1, 2, 4], [3]);
710     assert(approxEqual(res5.testStat, 0.1429));
711     assert(approxEqual(res5.p, 0.7418));
712 }
713 
714 // Levene's Test, Welch ANOVA and F test have massive overlap at the
715 // implementation level but less at the conceptual level, so I've combined
716 // the implementations into one horribly complicated but well-encapsulated
717 // templated function but left the interfaces as three unrelated functions.
718 private TestRes anovaLevene(bool levene, bool welch, alias central,  T...)
719 (T dataIn) {
720     static if(dataIn.length == 1) {
721         auto alloc = newRegionAllocator();
722         auto data = alloc.array(dataIn[0]);
723         auto withins = alloc.uninitializedArray!(MeanSD[])(data.length);
724         withins[] = MeanSD.init;
725     } else {
726         enum len = dataIn.length;
727         alias dataIn data;
728         MeanSD[len] withins;
729     }
730 
731     static if(levene) {
732         static if(dataIn.length == 1) {
733             auto centers = alloc.uninitializedArray!(double[])(data.length);
734         } else {
735             double[len] centers;
736         }
737 
738         foreach(i, category; data) {
739             static assert( isForwardRange!(typeof(category)) &&
740                 is(Unqual!(ElementType!(typeof(category))) : double),
741                 "Can only perform Levene's test on input ranges of elements " ~
742                 "implicitly convertible to doubles.");
743 
744             // The cast is to force conversions to double on alias this'd stuff
745             // like the Mean struct.
746             centers[i] = cast(double) central(category.save);
747         }
748 
749         double preprocess(double dataPoint, size_t category) {
750             return abs(dataPoint - centers[category]);
751         }
752     } else {
753         static double preprocess(double dataPoint, size_t category) {
754             return dataPoint;
755         }
756     }
757 
758 
759     auto DFGroups = data.length - 1;
760     ulong N = 0;
761 
762     foreach(category, range; data) {
763         static if(isInputRange!(typeof(range)) &&
764             is(Unqual!(ElementType!(typeof(range))) : double)) {
765             foreach(elem; range) {
766                 double preprocessed = preprocess(elem, category);
767                 withins[category].put(preprocessed);
768                 N++;
769             }
770         } else static if(isSummary!(typeof(range))) {
771             withins[category] = range.toMeanSD();
772             N += roundTo!long(range.N);
773         } else {
774             static assert(0, "Can only perform ANOVA on input ranges of " ~
775                 "numeric types, MeanSD structs and Summary structs, not a " ~
776                 typeof(range).stringof ~ ".");
777         }
778     }
779 
780     static if(!welch) {
781         immutable ulong DFDataPoints = N - data.length;
782         double mu = 0;
783         foreach(summary; withins) {
784             mu += summary.mean * (summary.N / N);
785         }
786 
787         double totalWithin = 0;
788         double totalBetween = 0;
789         foreach(group; withins) {
790             totalWithin += group.mse * (group.N / DFDataPoints);
791             immutable diffSq = (group.mean - mu) * (group.mean - mu);
792             totalBetween += diffSq * (group.N / DFGroups);
793         }
794 
795         immutable  F = totalBetween / totalWithin;
796         if(isNaN(F)) {
797             return TestRes.init;
798         }
799 
800         return TestRes(F, fisherCDFR(F, DFGroups, DFDataPoints));
801     } else {
802         immutable double k = data.length;
803         double sumW = 0;
804         foreach(summary; withins) {
805             sumW += summary.N / summary.var;
806         }
807 
808         double sumFt = 0;
809         foreach(summary; withins) {
810             sumFt += ((1 - summary.N / summary.var / sumW) ^^ 2) / (summary.N - 1);
811         }
812 
813         immutable kSqM1 = (k * k - 1.0);
814         immutable df2 = 1.0 / (3.0 / kSqM1 * sumFt);
815         immutable denom = 1 + 2 * (k - 2.0) / kSqM1 * sumFt;
816 
817         double yHat = 0;
818         foreach(i, summary; withins) {
819             yHat += summary.mean * (summary.N / summary.var);
820         }
821         yHat /= sumW;
822 
823         double numerator = 0;
824         foreach(i, summary; withins) {
825             numerator += summary.N / summary.var * ((summary.mean - yHat) ^^ 2);
826         }
827         numerator /= (k - 1);
828 
829         immutable F = numerator / denom;
830         if(isNaN(F)) {
831             return TestRes.init;
832         }
833 
834         return TestRes(F, fisherCDFR(F, DFGroups, df2));
835     }
836 }
837 
838 /**Performs a correlated sample (within-subjects) ANOVA.  This is a
839  * generalization of the paired T-test to 3 or more treatments.  This
840  * function accepts data as either a tuple of ranges (1 for each treatment,
841  * such that a given index represents the same subject in each range) or
842  * similarly as a range of ranges.
843  *
844  * Returns:  A TestRes with the F-statistic and P-value for the null that
845  * the the variable being measured did not vary across treatments against the
846  * alternative that it did.
847  *
848  * Examples:
849  * ---
850  * // Test the hypothesis that alcohol, loud music, caffeine and sleep
851  * // deprivation all have equivalent effects on programming ability.
852  *
853  * uint[] alcohol = [8,6,7,5,3,0,9];
854  * uint[] caffeine = [3,6,2,4,3,6,8];
855  * uint[] noSleep = [3,1,4,1,5,9,2];
856  * uint[] loudMusic = [2,7,1,8,2,8,1];
857  * // Subject 0 had ability of 8 under alcohol, 3 under caffeine, 3 under
858  * // no sleep, 2 under loud music.  Subject 1 had ability of 6 under alcohol,
859  * // 6 under caffeine, 1 under no sleep, and 7 under loud music, etc.
860  * auto result = correlatedAnova(alcohol, caffeine, noSleep, loudMusic);
861  * ---
862  *
863  * References:  "Concepts and Applications of Inferrential Statistics".
864  *              Richard Lowry.  Vassar College.   version.
865  *              http://faculty.vassar.edu/lowry/webtext.html
866  */
867 TestRes correlatedAnova(T...)(T dataIn)
868 if(allSatisfy!(isInputRange, T)) {
869     static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
870         auto alloc = newRegionAllocator();
871         auto data = alloc.array(dataIn[0]);
872         auto withins = alloc.newArray!(MeanSD[])(data.length);
873     } else {
874         enum len = dataIn.length;
875         alias dataIn data;
876         MeanSD[len] withins;
877     }
878     MeanSD overallSumm;
879     double nGroupNeg1 = 1.0 / data.length;
880 
881     bool someEmpty() {
882         foreach(elem; data) {
883             if(elem.empty) {
884                 return true;
885             }
886         }
887         return false;
888     }
889 
890     uint nSubjects = 0;
891     double subjSum = 0;
892     while(!someEmpty) {
893         double subjSumInner = 0;
894         foreach(i, elem; data) {
895             auto dataPoint = elem.front;
896             subjSumInner += dataPoint;
897             overallSumm.put(dataPoint);
898             withins[i].put(dataPoint);
899             data[i].popFront;
900         }
901         nSubjects++;
902         subjSum += subjSumInner * subjSumInner * nGroupNeg1;
903     }
904     double groupSum = 0;
905     foreach(elem; withins) {
906         groupSum += elem.mean * elem.N;
907     }
908 
909     groupSum /= sqrt(cast(double) nSubjects * data.length);
910     groupSum *= groupSum;
911     immutable subjErr = subjSum - groupSum;
912 
913     double betweenDev = 0;
914     immutable mu = overallSumm.mean;
915     foreach(group; withins) {
916         double diff = (group.mean - mu);
917         diff *= diff;
918         betweenDev += diff * (group.N / (data.length - 1));
919     }
920 
921     size_t errDf = data.length * nSubjects - data.length - nSubjects + 1;
922     double randError = -subjErr / errDf;
923     foreach(group; withins) {
924         randError += group.mse * (group.N / errDf);
925     }
926 
927     immutable F = betweenDev / randError;
928     if(!(F >= 0)) {
929         return TestRes(double.nan, double.nan);
930     }
931 
932     return TestRes(F, fisherCDFR(F, data.length - 1, errDf));
933 }
934 
935 unittest {
936     // Values from VassarStats utility at
937     // http://faculty.vassar.edu/lowry/VassarStats.html, but they like to
938     // round a lot, so the approxEqual tolerances are fairly wide.  I
939     // think it's adequate to demonstrate the correctness of this function,
940     // though.
941     uint[] alcohol = [8,6,7,5,3,0,9];
942     uint[] caffeine = [3,6,2,4,3,6,8];
943     uint[] noSleep = [3,1,4,1,5,9,2];
944     uint[] loudMusic = [2,7,1,8,2,8,1];
945     auto result = correlatedAnova(alcohol, caffeine, noSleep, loudMusic);
946     assert(approxEqual(result.testStat, 0.43, 0.0, 0.01));
947     assert(approxEqual(result.p, 0.734, 0.0, 0.01));
948 
949     uint[] stuff1 = [3,4,2,6];
950     uint[] stuff2 = [4,1,9,8];
951     auto result2 = correlatedAnova([stuff1, stuff2].dup);
952     assert(approxEqual(result2.testStat, 0.72, 0.0, 0.01));
953     assert(approxEqual(result2.p, 0.4584, 0.0, 0.01));
954 }
955 
956 /**The Kruskal-Wallis rank sum test.  Tests the null hypothesis that data in
957  * each group is not stochastically ordered with respect to data in each other
958  * groups.  This is a one-way non-parametric ANOVA and can be thought of
959  * as either a generalization of the Wilcoxon rank sum test to >2 groups or
960  * a non-parametric equivalent to the F-test.  Data can be input as either a
961  * tuple of ranges (one range for each group) or a range of ranges
962  * (one element for each group).
963  *
964  * Bugs:  Asymptotic approximation of P-value only, not exact.  In this case,
965  * I'm not sure a practical way to compute the exact P-value even exists.
966  *
967  * Returns:  A TestRes with the K statistic and the P-value for the null that
968  * no group is stochastically larger than any other against the alternative that
969  * groups are stochastically ordered.
970  *
971  * References:  "Concepts and Applications of Inferrential Statistics".
972  *              Richard Lowry.  Vassar College.   version.
973  *              http://faculty.vassar.edu/lowry/webtext.html
974  */
975 TestRes kruskalWallis(T...)(T dataIn)
976 if(doubleInput!(typeof(dataIn[0].front)) || allSatisfy!(doubleInput, T)) {
977     auto alloc = newRegionAllocator();
978     size_t N = 0;
979 
980     static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
981         auto data = alloc.array(dataIn[0]);
982         alias ElementType!(typeof(data[0])) C;
983         static if(hasLength!(typeof(data[0]))) {
984             enum bool useLength = true;
985         } else {
986             enum bool useLength = false;
987         }
988     } else {
989         enum len = dataIn.length;
990         alias dataIn data;
991         alias staticMap!(ElementType, T) Es;
992         alias CommonType!(Es) C;
993         static if(allSatisfy!(hasLength, T)) {
994             enum bool useLength = true;
995         } else {
996             enum bool useLength = false;
997         }
998     }
999 
1000     size_t[] lengths = alloc.uninitializedArray!(size_t[])(data.length);
1001     static if(useLength) {
1002         foreach(i, rng; data) {
1003             auto rngLen = rng.length;
1004             lengths[i] = rngLen;
1005             N += rngLen;
1006         }
1007         auto dataArray = alloc.uninitializedArray!(Unqual!(C)[])(N);
1008         size_t pos = 0;
1009         foreach(rng; data) {
1010             foreach(elem; rng) {
1011                 dataArray[pos++] = elem;
1012             }
1013         }
1014     } else {
1015         auto app = appender!(Unqual!(C)[])();
1016         foreach(i, rng; data) {
1017             size_t oldLen = dataArray.length;
1018             app.put(rng);
1019             lengths[i] = dataArray.length - oldLen;
1020             N += lengths[i];
1021         }
1022         auto dataArray = app.data;
1023     }
1024 
1025     double[] ranks = alloc.uninitializedArray!(double[])(dataArray.length);
1026     try {
1027         rankSort(dataArray, ranks);
1028     } catch(SortException) {
1029         return TestRes.init;
1030     }
1031 
1032     size_t index = 0;
1033     double denom = 0, numer = 0;
1034     double rBar = 0.5 * (N + 1);
1035     foreach(meanI, l; lengths) {
1036         Mean groupStats;
1037         foreach(i; index..index + l) {
1038             groupStats.put( ranks[i]);
1039             double diff = ranks[i] - rBar;
1040             diff *= diff;
1041             denom += diff;
1042         }
1043         index += l;
1044         double nDiff = groupStats.mean - rBar;
1045         nDiff *= nDiff;
1046         numer += l * nDiff;
1047     }
1048     double K = (N - 1) * (numer / denom);
1049 
1050     // Tie correction.
1051     double tieSum = 0;
1052     uint nTies = 1;
1053     foreach(i; 1..dataArray.length) {
1054         if(dataArray[i] == dataArray[i - 1]) {
1055             nTies++;
1056         } else if(nTies > 1) {
1057             double partialSum = nTies;
1058             partialSum = (partialSum * partialSum * partialSum) - partialSum;
1059             tieSum += partialSum;
1060             nTies = 1;
1061         }
1062     }
1063     if(nTies > 1) {
1064         double partialSum = nTies;
1065         partialSum = (partialSum * partialSum * partialSum) - partialSum;
1066         tieSum += partialSum;
1067     }
1068     double tieDenom = N;
1069     tieDenom = (tieDenom * tieDenom * tieDenom) - tieDenom;
1070     tieSum = 1 - (tieSum / tieDenom);
1071     K *= tieSum;
1072 
1073     if(isNaN(K)) {
1074         return TestRes(double.nan, double.nan);
1075     }
1076 
1077     return TestRes(K, chiSquareCDFR(K, data.length - 1));
1078 }
1079 
1080 unittest {
1081     // These values are from the VassarStat web tool at
1082     // http://faculty.vassar.edu/lowry/VassarStats.html .
1083     // R is actually wrong here because it apparently doesn't use a correction
1084     // for ties.
1085     auto res1 = kruskalWallis([3,1,4,1].idup, [5,9,2,6].dup, [5,3,5].dup);
1086     assert(approxEqual(res1.testStat, 4.15));
1087     assert(approxEqual(res1.p, 0.1256));
1088 
1089     // Test for other input types.
1090     auto res2 = kruskalWallis([[3,1,4,1].idup, [5,9,2,6].idup, [5,3,5].idup].dup);
1091     assert(res2 == res1);
1092     auto res3 = kruskalWallis(map!"a"([3,1,4,1].dup), [5,9,2,6].dup, [5,3,5].dup);
1093     assert(res3 == res1);
1094     auto res4 = kruskalWallis([map!"a"([3,1,4,1].dup),
1095                                map!"a"([5,9,2,6].dup),
1096                                map!"a"([5,3,5].dup)].dup);
1097     assert(res4 == res1);
1098 
1099     // Test w/ one more case, just with one input type.
1100     auto res5 = kruskalWallis([2,7,1,8,2].dup, [8,1,8,2].dup, [8,4,5,9,2].dup,
1101                               [7,1,8,2,8,1,8].dup);
1102     assert(approxEqual(res5.testStat, 1.06));
1103     assert(approxEqual(res5.p, 0.7867));
1104 }
1105 
1106 /**The Friedman test is a non-parametric within-subject ANOVA.  It's useful
1107  * when parametric assumptions cannot be made.  Usage is identical to
1108  * correlatedAnova().
1109  *
1110  * References:  "Concepts and Applications of Inferrential Statistics".
1111  *              Richard Lowry.  Vassar College.   version.
1112  *              http://faculty.vassar.edu/lowry/webtext.html
1113  *
1114  * Bugs:  No exact P-value calculation.  Asymptotic approx. only.
1115  */
1116 TestRes friedmanTest(T...)(T dataIn)
1117 if(doubleInput!(typeof(dataIn[0].front)) || allSatisfy!(doubleInput, T)) {
1118     static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
1119         auto alloc = newRegionAllocator();
1120         auto data = alloc.array(dataIn[0]);
1121         auto ranks = alloc.uninitializedArray!(double[])(data.length);
1122         auto dataPoints = alloc.uninitializedArray!(double[])(data.length);
1123         auto colMeans = alloc.newArray!(Mean[])(data.length);
1124     } else {
1125         enum len = dataIn.length;
1126         alias dataIn data;
1127         double[len] ranks;
1128         double[len] dataPoints;
1129         Mean[len] colMeans;
1130     }
1131     double rBar = cast(double) data.length * (data.length + 1.0) / 2.0;
1132     MeanSD overallSumm;
1133 
1134     bool someEmpty() {
1135         foreach(elem; data) {
1136             if(elem.empty) {
1137                 return true;
1138             }
1139         }
1140         return false;
1141     }
1142 
1143     uint N = 0;
1144     while(!someEmpty) {
1145         foreach(i, range; data) {
1146             dataPoints[i] = data[i].front;
1147             data[i].popFront;
1148         }
1149 
1150         try {
1151             rank(dataPoints[], ranks[]);
1152         } catch(SortException) {
1153             return TestRes.init;
1154         }
1155 
1156         foreach(i, rank; ranks) {
1157             colMeans[i].put(rank);
1158             overallSumm.put(rank);
1159         }
1160         N++;
1161     }
1162 
1163     double between = 0;
1164     double mu = overallSumm.mean;
1165     foreach(mean; colMeans) {
1166         double diff = mean.mean - overallSumm.mean;
1167         between += diff * diff;
1168     }
1169     between *= N;
1170     double within = overallSumm.mse * (overallSumm.N / (overallSumm.N - N));
1171     double chiSq = between / within;
1172     double df = data.length - 1;
1173     return TestRes(chiSq, chiSquareCDFR(chiSq, df));
1174 }
1175 
1176 unittest {
1177     // Values from R
1178     uint[] alcohol = [8,6,7,5,3,0,9];
1179     uint[] caffeine = [3,6,2,4,3,6,8];
1180     uint[] noSleep = [3,1,4,1,5,9,2];
1181     uint[] loudMusic = [2,7,1,8,2,8,1];
1182     auto result = friedmanTest(alcohol, caffeine, noSleep, loudMusic);
1183     assert(approxEqual(result.testStat, 1.7463));
1184     assert(approxEqual(result.p, 0.6267));
1185 
1186     uint[] stuff1 = [3,4,2,6];
1187     uint[] stuff2 = [4,1,9,8];
1188     auto result2 = friedmanTest([stuff1, stuff2].dup);
1189     assert(approxEqual(result2.testStat, 1));
1190     assert(approxEqual(result2.p, 0.3173));
1191 }
1192 
1193 /**Computes Wilcoxon rank sum test statistic and P-value for
1194  * a set of observations against another set, using the given alternative.
1195  * Alt.less means that sample1 is stochastically less than sample2.
1196  * Alt.greater means sample1 is stochastically greater than sample2.
1197  * Alt.twoSided means sample1 is stochastically less than or greater than
1198  * sample2.
1199  *
1200  * exactThresh is the threshold value of (n1 + n2) at which this function
1201  * switches from exact to approximate computation of the p-value.  Do not set
1202  * exactThresh to more than 200, as the exact
1203  * calculation is both very slow and not numerically stable past this point,
1204  * and the asymptotic calculation is very good for N this large.  To disable
1205  * exact calculation entirely, set exactThresh to 0.
1206  *
1207  * Notes:  Exact p-value computation is never used when ties are present in the
1208  * data, because it is not computationally feasible.
1209  *
1210  * Input ranges for this function must define a length.
1211  *
1212  * This test is also known as the Mann-Whitney U test.
1213  *
1214  * Returns:  A TestRes containing the W test statistic and the P-value against
1215  * the given alternative.
1216  *
1217  * References:  http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U
1218  *
1219  * StackOverflow Question 376003  http://stackoverflow.com/questions/376003
1220  *
1221  * Loughborough University MLSC Statistics 2.3 The Mann-Whitney U Test
1222  * http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf
1223  */
1224 TestRes wilcoxonRankSum(T, U)
1225 (T sample1, U sample2, Alt alt = Alt.twoSided, uint exactThresh = 50)
1226 if(isInputRange!T && isInputRange!U &&
1227 is(typeof(sample1.front < sample2.front) == bool) &&
1228 is(CommonType!(ElementType!T, ElementType!U))) {
1229 
1230     auto alloc = newRegionAllocator();
1231     alias Unqual!(CommonType!(ElementType!(T), ElementType!(U))) C;
1232 
1233     static if(hasLength!T && hasLength!U) {
1234         auto n1 = sample1.length, n2 = sample2.length, N = n1 + n2;
1235         auto combined = alloc.uninitializedArray!(C[])(N);
1236         copy(sample1, combined[0..n1]);
1237         copy(sample2, combined[n1..$]);
1238     } else {
1239         auto app = appender!(C[])();
1240 
1241         foreach(elem; sample1) {
1242             app.put(elem);
1243         }
1244 
1245         uint n1 = app.data.length;
1246         foreach(elem; sample2) {
1247             app.put(elem);
1248         }
1249 
1250         auto combined = app.data;
1251         uint N = combined.length;
1252         uint n2 = N - n1;
1253     }
1254 
1255     double[] ranks = alloc.uninitializedArray!(double[])(N);
1256     try {
1257         rankSort(combined, ranks);
1258     } catch(SortException) {
1259         return TestRes.init;
1260     }
1261     double w = reduce!("a + b")
1262         (0.0, ranks[0..n1]) - cast(ulong) n1 * (n1 + 1) / 2UL;
1263 
1264     if(alt == Alt.none) {
1265         return TestRes(w);
1266     }
1267 
1268     double tieSum = 0;
1269     // combined is sorted by rankSort.  Can use it to figure out how many
1270     // ties we have w/o another allocation or sorting.
1271     enum oneOverTwelve = 1.0 / 12.0;
1272     tieSum = 0;
1273     ulong nties = 1;
1274     foreach(i; 1..N) {
1275         if(combined[i] == combined[i - 1]) {
1276             nties++;
1277         } else {
1278             if(nties == 1)
1279                 continue;
1280             tieSum += ((nties * nties * nties) - nties) * oneOverTwelve;
1281             nties = 1;
1282         }
1283     }
1284     // Handle last run.
1285     if(nties > 1) {
1286         tieSum += ((nties * nties * nties) - nties) * oneOverTwelve;
1287     }
1288 
1289     immutable p = wilcoxonRankSumPval(w, n1, n2, alt, tieSum, exactThresh);
1290     return TestRes(w, p);
1291 }
1292 
1293  unittest {
1294      // Values from R.
1295 
1296     assert(wilcoxonRankSum([1, 2, 3, 4, 5].dup, [2, 4, 6, 8, 10].dup).testStat == 5);
1297     assert(wilcoxonRankSum([2, 4, 6, 8, 10].dup, [1, 2, 3, 4, 5].dup).testStat == 20);
1298     assert(wilcoxonRankSum([3, 7, 21, 5, 9].dup, [2, 4, 6, 8, 10].dup).testStat == 15);
1299 
1300      // Simple stuff (no ties) first.  Testing approximate
1301      // calculation first.
1302      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1303            Alt.twoSided, 0), 0.9273));
1304      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1305            Alt.less, 0), 0.6079));
1306      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1307            Alt.greater, 0).p, 0.4636));
1308      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
1309             Alt.twoSided, 0).p, 0.4113));
1310      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
1311             Alt.less, 0).p, 0.2057));
1312      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup,
1313         map!"a"([3,5,7,8,13,15].dup), Alt.greater, 0).p, 0.8423));
1314      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1315             Alt.twoSided, 0), .6745));
1316      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1317             Alt.less, 0), .3372));
1318      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1319             Alt.greater, 0), .7346));
1320 
1321     // Now, lots of ties.
1322     assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
1323            Alt.twoSided, 0), 0.3976));
1324     assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
1325            Alt.less, 0), 0.1988));
1326     assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
1327            Alt.greater, 0), 0.8548));
1328     assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
1329            Alt.twoSided, 0), 0.9049));
1330     assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
1331            Alt.less, 0), 0.4524));
1332     assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
1333            Alt.greater, 0), 0.64));
1334 
1335     // Now, testing the exact calculation on the same data.
1336      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1337        Alt.twoSided), 0.9307));
1338      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1339            Alt.less), 0.6039));
1340      assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
1341            Alt.greater), 0.4654));
1342      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
1343             Alt.twoSided), 0.4286));
1344      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
1345             Alt.less), 0.2143));
1346      assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
1347             Alt.greater), 0.8355));
1348      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1349             Alt.twoSided), .6905));
1350      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1351             Alt.less), .3452));
1352      assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
1353             Alt.greater), .7262));
1354 }
1355 
1356 private
1357 double wilcoxonRankSumPval(double w, ulong n1, ulong n2, Alt alt = Alt.twoSided,
1358                            double tieSum = 0,  uint exactThresh = 50) {
1359     if(alt == Alt.none) {
1360         return double.nan;
1361     }
1362 
1363     immutable double N = n1 + n2;
1364 
1365     if(N < exactThresh && tieSum == 0) {
1366         return wilcoxRSPExact(roundTo!uint(w), cast(uint) n1, cast(uint) n2, alt);
1367     }
1368 
1369     immutable sd = sqrt(cast(double) (n1 * n2) / (N * (N - 1)) *
1370              ((N * N * N - N) / 12 - tieSum));
1371 
1372     // Can happen if all samples are tied.
1373     if(!(sd > 0)) {
1374         return double.nan;
1375     }
1376 
1377     immutable mean = (n1 * n2) / 2.0;
1378 
1379     if(alt == Alt.twoSided) {
1380         if(abs(w - mean) < 0.5) {
1381             return 1;
1382         } else if(w < mean) {
1383             return 2 * normalCDF(w + 0.5, mean, sd);
1384         } else {
1385             assert(w > mean);
1386             return 2 * normalCDFR(w - 0.5, mean, sd);
1387         }
1388     } else if(alt == Alt.less) {
1389         return normalCDF(w + 0.5, mean, sd);
1390     } else if(alt == Alt.greater) {
1391         return normalCDFR(w - 0.5, mean, sd);
1392     }
1393 
1394     assert(0);
1395 }
1396 
1397 unittest {
1398     /* Values from R.  I could only get good values for Alt.less directly.
1399      * Using W-values to test Alt.twoSided, Alt.greater indirectly.*/
1400     assert(approxEqual(wilcoxonRankSumPval(1200, 50, 50, Alt.less), .3670));
1401     assert(approxEqual(wilcoxonRankSumPval(1500, 50, 50, Alt.less), .957903));
1402     assert(approxEqual(wilcoxonRankSumPval(8500, 100, 200, Alt.less), .01704));
1403     auto w = wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup).testStat;
1404     assert(approxEqual(wilcoxonRankSumPval(w, 5, 6), 0.9273));
1405     assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.greater), 0.4636));
1406     assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.less), 0.6079));
1407 
1408     // Monte carlo unit testing:  Make sure that the exact and asymptotic
1409     // versions agree within a small epsilon;
1410     double maxEpsilon = 0;
1411     foreach(i; 0..1_000) {
1412         uint n1 = uniform(5U, 25U);
1413         uint n2 = uniform(5U, 25U);
1414         uint testStat = uniform!"[]"(0, (n1 * n2));
1415 
1416         foreach(alt; [Alt.less, Alt.greater, Alt.twoSided]) {
1417             double approxP = wilcoxonRankSumPval(testStat, n1, n2, alt, 0, 0);
1418             double exactP = wilcoxonRankSumPval(testStat, n1, n2, alt, 0, 50);
1419             double epsilon = abs(approxP - exactP);
1420             assert(epsilon < 0.02);
1421             maxEpsilon = max(maxEpsilon, epsilon);
1422         }
1423     }
1424 }
1425 
1426 /* Used internally by wilcoxonRankSum.  This function uses dynamic
1427  * programming to count the number of combinations of numbers [1..N] that sum
1428  * of length n1 that sum to <= W in O(N * W * n1) time.
1429  * Algorithm obtained from StackOverflow Question 376003
1430  * (http://stackoverflow.com/questions/376003).*/
1431 private double wilcoxRSPExact(uint W, uint n1, uint n2, Alt alt = Alt.twoSided) {
1432     uint N = n1 + n2;
1433     immutable maxPossible = n1 * n2;
1434 
1435     switch(alt) {
1436         case Alt.less:
1437             if(W >= maxPossible)  { // Value impossibly large
1438                 return 1;
1439             } else if(W * 2 <= maxPossible) {
1440                 break;
1441             } else {
1442                 return 1 - wilcoxRSPExact(maxPossible - W - 1, n1, n2, Alt.less);
1443             }
1444             assert(0);
1445         case Alt.greater:
1446             if(W > maxPossible)  { // Value impossibly large
1447                 return 0;
1448             } else if(W * 2 >= maxPossible) {
1449                 return wilcoxRSPExact(maxPossible - W, n1, n2, Alt.less);
1450             } else if(W <= 0) {
1451                 return 1;
1452             } else {
1453                 return 1 - wilcoxRSPExact(W - 1, n1, n2, Alt.less);
1454             }
1455             assert(0);
1456         case Alt.twoSided:
1457             if(W * 2 <= maxPossible) {
1458                 return min(1, wilcoxRSPExact(W, n1, n2, Alt.less) +
1459                        wilcoxRSPExact(maxPossible - W, n1, n2, Alt.greater));
1460             } else {
1461                 return min(1, wilcoxRSPExact(W, n1, n2, Alt.greater) +
1462                        wilcoxRSPExact(maxPossible - W, n1, n2, Alt.less));
1463             }
1464             assert(0);
1465         default:
1466             assert(0);
1467     }
1468 
1469     W += n1 * (n1 + 1) / 2UL;
1470 
1471     auto alloc = newRegionAllocator();
1472     float* cache = (alloc.uninitializedArray!(float[])((n1 + 1) * (W + 1))).ptr;
1473     float* cachePrev = (alloc.uninitializedArray!(float[])((n1 + 1) * (W + 1))).ptr;
1474     cache[0..(n1 + 1) * (W + 1)] = 0;
1475     cachePrev[0..(n1 + 1) * (W + 1)] = 0;
1476 
1477     /* Using doubles for the intermediate steps is too slow, but I didn't want to
1478      * lose too much precision.  Since my sums must be between 0 and 1, I am
1479      * using the entire bit space of a float to hold numbers between zero and
1480      * one.  This is precise to at least 1e-7.  This is good enough for a few
1481      * reasons:
1482      *
1483      * 1.  This is a p-value, and therefore will likely not be used in
1484      *     further calculations where rounding error would accumulate.
1485      * 2.  If this is too slow, the alternative is to use the asymptotic
1486      *     approximation.  This is can have relative errors of several orders
1487      *     of magnitude in the tails of the distribution, and is therefore
1488      *     clearly worse.
1489      * 3.  For very large N, where this function could give completely wrong
1490      *     answers, it would be so slow that any reasonable person would use the
1491      *     asymptotic approximation anyhow.*/
1492 
1493 
1494     // Algorithm based on StackOverflow question 376003.
1495     double comb = exp(-logNcomb(N, n1));
1496     double floatMax = cast(double) float.max;
1497     cache[0] = cast(float) (comb * floatMax);
1498     cachePrev[0] = cast(float) (comb * floatMax);
1499 
1500     foreach(i; 1..N + 1) {
1501         swap(cache, cachePrev);
1502         foreach(k; 1..min(i + 1, n1 + 1)) {
1503 
1504             uint minW = k * (k + 1) / 2;
1505             float* curK = cache + k * (W + 1);
1506             float* prevK = cachePrev + k * (W + 1);
1507             float* prevKm1 = cachePrev + (k - 1) * (W + 1);
1508 
1509             foreach(w; minW..W + 1) {
1510                 curK[w] = prevK[w] + ((i <= w) ? prevKm1[w - i] : 0);
1511             }
1512         }
1513     }
1514 
1515     double sum = 0;
1516     float* lastLine = cache + n1 * (W + 1);
1517     foreach(w; 1..W + 1) {
1518         sum += (cast(double) lastLine[w] / floatMax);
1519     }
1520 
1521     return sum;
1522 }
1523 
1524 unittest {
1525     // Values from R.
1526     assert(approxEqual(wilcoxRSPExact(14, 5, 6), 0.9307));
1527     assert(approxEqual(wilcoxRSPExact(14, 5, 6, Alt.less), 0.4654));
1528     assert(approxEqual(wilcoxRSPExact(14, 5, 6, Alt.greater), 0.6039));
1529     assert(approxEqual(wilcoxRSPExact(16, 6, 5), 0.9307));
1530     assert(approxEqual(wilcoxRSPExact(16, 6, 5, Alt.less), 0.6039));
1531     assert(approxEqual(wilcoxRSPExact(16, 6, 5, Alt.greater), 0.4654));
1532     assert(approxEqual(wilcoxRSPExact(66, 10, 35, Alt.less), 0.001053));
1533     assert(approxEqual(wilcoxRSPExact(78, 13, 6, Alt.less), 1));
1534 
1535     // Mostly to make sure that underflow doesn't happen until
1536     // the N's are truly unreasonable:
1537     //assert(approxEqual(wilcoxRSPExact(6_000, 120, 120, Alt.less), 0.01276508));
1538 }
1539 
1540 /**Computes a test statistic and P-value for a Wilcoxon signed rank test against
1541  * the given alternative. Alt.less means that elements of before are stochastically
1542  * less than corresponding elements of after.  Alt.greater means elements of
1543  * before are stochastically greater than corresponding elements of after.
1544  * Alt.twoSided means there is a significant difference in either direction.
1545  *
1546  * exactThresh is the threshold value of before.length at which this function
1547  * switches from exact to approximate computation of the p-value.   Do not set
1548  * exactThresh to more than 200, as the exact calculation is both very slow and
1549  * not numerically stable past this point, and the asymptotic calculation is
1550  * very good for N this large.  To disable exact calculation entirely, set
1551  * exactThresh to 0.
1552  *
1553  * Notes:  Exact p-value computation is never used when ties are present,
1554  * because it is not computationally feasible.
1555  *
1556  * The input ranges for this function must define a length and must be
1557  * forward ranges.
1558  *
1559  * Returns:  A TestRes of the W statistic and the p-value against the given
1560  * alternative.
1561  *
1562  * References:  http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test
1563  *
1564  * StackOverflow Question 376003  http://stackoverflow.com/questions/376003
1565  *
1566  * Handbook of Parametric and nonparametric statistical procedures. David Sheskin.
1567  * Third Edition. (2004)  CRC Press. Pg. 616.
1568  */
1569 TestRes wilcoxonSignedRank(T, U)(T before, U after, Alt alt = Alt.twoSided, uint exactThresh = 50)
1570 if(doubleInput!(T) && doubleInput!(U) &&
1571 is(typeof(before.front - after.front) : double)) {
1572     uint nZero = 0;
1573     byte sign(double input) {
1574         if(input < 0)
1575             return -1;
1576         if(input > 0)
1577             return 1;
1578         nZero++;
1579         return 0;
1580     }
1581 
1582     auto alloc = newRegionAllocator();
1583 
1584     static if(hasLength!T && hasLength!U) {
1585         dstatsEnforce(before.length == after.length,
1586             "Ranges must have same lengths for wilcoxonSignedRank.");
1587 
1588         double[] diffRanks = alloc.uninitializedArray!(double[])(before.length);
1589         byte[] signs = alloc.uninitializedArray!(byte[])(before.length);
1590         double[] diffs = alloc.uninitializedArray!(double[])(before.length);
1591 
1592         size_t ii = 0;
1593         while(!before.empty && !after.empty) {
1594             double diff = cast(double) before.front - cast(double) after.front;
1595             signs[ii] = sign(diff);
1596             diffs[ii] = abs(diff);
1597             ii++;
1598             before.popFront;
1599             after.popFront;
1600         }
1601     } else {
1602         double[] diffRanks;
1603         auto diffApp = appender!(double[])();
1604         auto signApp = appender!(byte[])();
1605 
1606         while(!before.empty && !after.empty) {
1607             double diff = cast(double) before.front - cast(double) after.front;
1608             signApp.put(sign(diff));
1609             diffApp.put(abs(diff));
1610             before.popFront;
1611             after.popFront;
1612         }
1613 
1614         auto diffs = diffApp.data;
1615         auto signs = signApp.data;
1616         diffRanks = alloc.uninitializedArray!(double[])(diffs.length);
1617     }
1618     try {
1619         rankSort(diffs, diffRanks);
1620     } catch(SortException) {
1621         return TestRes.init;
1622     }
1623 
1624     ulong N = diffs.length - nZero;
1625 
1626     double W = 0;
1627     foreach(i, dr; diffRanks) {
1628         if(signs[i] == 1) {
1629             W += dr - nZero;
1630         }
1631     }
1632 
1633     // Just a sanity check.  Should be mathematically impossible for this
1634     // assert to fail.  The 1e-5 is for round-off error.
1635     assert(W > -1e-5 && W <= (N * (N + 1) / 2) + 1e-5);
1636 
1637     if(alt == Alt.none) {
1638         return TestRes(W);
1639     }
1640 
1641     // Handle ties.
1642     double tieSum = 0;
1643 
1644     // combined is sorted by rankSort.  Can use it to figure out how many
1645     // ties we have w/o another allocation or sorting.
1646     enum denom = 1.0 / 48.0;
1647     ulong nties = 1;
1648     foreach(i; 1..diffs.length) {
1649         if(diffs[i] == diffs[i - 1] && diffs[i] != 0) {
1650             nties++;
1651         } else {
1652             if(nties == 1)
1653                 continue;
1654             tieSum += ((nties * nties * nties) - nties) * denom;
1655             nties = 1;
1656         }
1657     }
1658     // Handle last run.
1659     if(nties > 1) {
1660         tieSum += ((nties * nties * nties) - nties) * denom;
1661     }
1662     if(nZero > 0 && tieSum == 0) {
1663         tieSum = double.nan;  // Signal that there were zeros and exact p-val can't be computed.
1664     }
1665 
1666     return TestRes(W, wilcoxonSignedRankPval(W, N, alt, tieSum, exactThresh));
1667 }
1668 
1669 unittest {
1670     // Values from R.
1671     alias approxEqual ae;
1672     assert(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup).testStat == 7.5);
1673     assert(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup).testStat == 6);
1674     assert(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup).testStat == 5);
1675 
1676     // With ties, normal approx.
1677     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup), 1));
1678     assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, map!"a"([2,7,1,8,2].dup)), 0.7865));
1679     assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup), 0.5879));
1680     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.less), 0.5562));
1681     assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.less), 0.3932));
1682     assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.less), 0.2940));
1683     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.greater), 0.5562));
1684     assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.greater), 0.706));
1685     assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.greater), 0.7918));
1686     assert(ae(wilcoxonSignedRank(cast(int[]) [1,16,2,4,8], cast(int[]) [1,5,2,3,4]).testStat, 6));
1687     assert(ae(wilcoxonSignedRank(cast(int[]) [1,16,2,4,8], cast(int[]) [1,5,2,3,4]), 0.1814));
1688 
1689     // Exact.
1690     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup), 0.625));
1691     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.less), 0.3125));
1692     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.greater), 0.7812));
1693     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup), 0.8125));
1694     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.less), 0.6875));
1695     assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.greater), 0.4062));
1696 
1697     // Monte carlo unit testing.  Make sure exact, approx are really,
1698     // really close to each other.
1699     double maxEpsilon = 0;
1700     foreach(i; 0..1_000) {
1701         uint N = uniform(10U, 50U);
1702         uint testStat = uniform!"[]"(0, N * (N + 1) / 2);
1703 
1704         foreach(alt; [Alt.less, Alt.greater, Alt.twoSided]) {
1705             double approxP = wilcoxonSignedRankPval(testStat, N, alt, 0, 0);
1706             double exactP = wilcoxonSignedRankPval(testStat, N, alt, 0, 50);
1707             double epsilon = abs(approxP - exactP);
1708             assert(epsilon < 0.02);
1709             maxEpsilon = max(maxEpsilon, epsilon);
1710         }
1711     }
1712 }
1713 
1714 /**Same as the overload, but allows testing whether a range is stochastically
1715  * less than or greater than a fixed value mu rather than paired elements of
1716  * a second range.*/
1717 TestRes wilcoxonSignedRank(T)(T data, double mu, Alt alt = Alt.twoSided, uint exactThresh = 50)
1718 if(doubleInput!(T) && is(typeof(data.front - mu) : double)) {
1719     return wilcoxonSignedRank(data, repeat(mu, data.length), alt, exactThresh);
1720 }
1721 
1722 unittest {
1723     auto res = wilcoxonSignedRank([-8,-6,2,4,7].dup, 0);
1724     assert(approxEqual(res.testStat, 7));
1725     assert(approxEqual(res.p, 1));
1726 }
1727 
1728 private double wilcoxonSignedRankPval(double W, ulong N, Alt alt = Alt.twoSided,
1729      double tieSum = 0, uint exactThresh = 50)
1730 in {
1731     assert(N > 0);
1732     assert(tieSum >= 0 || isNaN(tieSum));
1733 } do {
1734     if(alt == Alt.none) {
1735         return double.nan;
1736     }
1737 
1738     if(tieSum == 0 && !isNaN(tieSum) && N <= exactThresh) {
1739         return wilcoxSRPExact(roundTo!uint(W), to!uint(N), alt);
1740     }
1741 
1742     if(isNaN(tieSum)) {
1743         tieSum = 0;
1744     }
1745 
1746     immutable expected = N * (N + 1) * 0.25;
1747     immutable sd = sqrt(N * (N + 1) * (2 * N + 1) / 24.0 - tieSum);
1748 
1749     if(alt == Alt.less) {
1750         return normalCDF(W + 0.5, expected, sd);
1751     } else if(alt == Alt.greater) {
1752         return normalCDFR(W - 0.5, expected, sd);
1753     } else {
1754         assert(alt == Alt.twoSided);
1755         if(abs(W - expected) <= 0.5) {
1756             return 1;
1757         } else if(W < expected) {
1758             return 2 * normalCDF(W + 0.5, expected, sd);
1759         } else {
1760             assert(W > expected);
1761             return 2 * normalCDFR(W - 0.5, expected, sd);
1762         }
1763     }
1764 }
1765 // Tested indirectly through other overload.
1766 
1767 /* Yes, a little cut and paste coding was involved here from wilcoxRSPExact,
1768  * but this function and wilcoxRSPExact are just different enough that
1769  * it would be more trouble than it's worth to write one generalized
1770  * function.
1771  *
1772  * Algorithm adapted from StackOverflow question 376003
1773  * (http://stackoverflow.com/questions/376003).
1774  */
1775 private double wilcoxSRPExact(uint W, uint N, Alt alt = Alt.twoSided) {
1776     immutable maxPossible = N * (N + 1) / 2;
1777 
1778     switch(alt) {
1779         case Alt.less:
1780             if(W >= maxPossible)  { // Value impossibly large
1781                 return 1;
1782             } else if(W * 2 <= maxPossible) {
1783                 break;
1784             } else {
1785                 return 1 - wilcoxSRPExact(maxPossible - W - 1, N, Alt.less);
1786             }
1787         case Alt.greater:
1788             if(W > maxPossible)  { // Value impossibly large
1789                 return 0;
1790             } else if(W == 0) {
1791                 return 1;
1792             } else if(W * 2 >= maxPossible) {
1793                 return wilcoxSRPExact(maxPossible - W, N, Alt.less);
1794             } else {
1795                 return 1 - wilcoxSRPExact(W - 1, N, Alt.less);
1796             }
1797         case Alt.twoSided:
1798             if(W * 2 <= maxPossible) {
1799                 return min(1, wilcoxSRPExact(W, N, Alt.less) +
1800                        wilcoxSRPExact(maxPossible - W, N, Alt.greater));
1801             } else {
1802                 return min(1, wilcoxSRPExact(W, N, Alt.greater) +
1803                        wilcoxSRPExact(maxPossible - W, N, Alt.less));
1804             }
1805         default:
1806             assert(0);
1807     }
1808 
1809     auto alloc = newRegionAllocator();
1810     float* cache = (alloc.uninitializedArray!(float[])((N + 1) * (W + 1))).ptr;
1811     float* cachePrev = (alloc.uninitializedArray!(float[])((N + 1) * (W + 1))).ptr;
1812     cache[0..(N + 1) * (W + 1)] = 0;
1813     cachePrev[0..(N + 1) * (W + 1)] = 0;
1814 
1815     double comb = pow(2.0, -(cast(double) N));
1816     double floatMax = cast(double) float.max;
1817     cache[0] = cast(float) (comb * floatMax);
1818     cachePrev[0] = cast(float) (comb * floatMax);
1819 
1820     foreach(i; 1..N + 1) {
1821         swap(cache, cachePrev);
1822         foreach(k; 1..i + 1) {
1823 
1824             uint minW = k * (k + 1) / 2;
1825             float* curK = cache + k * (W + 1);
1826             float* prevK = cachePrev + k * (W + 1);
1827             float* prevKm1 = cachePrev + (k - 1) * (W + 1);
1828 
1829             foreach(w; minW..W + 1) {
1830                 curK[w] = prevK[w] + ((i <= w) ? prevKm1[w - i] : 0);
1831             }
1832         }
1833     }
1834 
1835     double sum  = 0;
1836     foreach(elem; cache[0..(N + 1) * (W + 1)]) {
1837         sum += cast(double) elem / (cast(double) float.max);
1838     }
1839 
1840     return sum;
1841 }
1842 
1843 unittest {
1844     // Values from R.
1845     assert(approxEqual(wilcoxSRPExact(25, 10, Alt.less), 0.4229));
1846     assert(approxEqual(wilcoxSRPExact(25, 10, Alt.greater), 0.6152));
1847     assert(approxEqual(wilcoxSRPExact(25, 10, Alt.twoSided), 0.8457));
1848     assert(approxEqual(wilcoxSRPExact(31, 10, Alt.less), 0.6523));
1849     assert(approxEqual(wilcoxSRPExact(31, 10, Alt.greater), 0.3848));
1850     assert(approxEqual(wilcoxSRPExact(31, 10, Alt.twoSided), 0.7695));
1851 }
1852 
1853 /**Sign test for differences between paired values.  This is a very robust
1854  * but very low power test.  Alternatives are Alt.less, meaning elements
1855  * of before are typically less than corresponding elements of after,
1856  * Alt.greater, meaning elements of before are typically greater than
1857  * elements of after, and Alt.twoSided, meaning that there is a significant
1858  * difference in either direction.
1859  *
1860  * Returns:  A TestRes with the proportion of elements of before that were
1861  * greater than the corresponding element of after, and the P-value against
1862  * the given alternative.
1863  */
1864 TestRes signTest(T, U)(T before, U after, Alt alt = Alt.twoSided)
1865 if(doubleInput!(T) && doubleInput!(U) &&
1866 is(typeof(before.front < after.front) == bool)) {
1867     ulong greater, less;
1868     while(!before.empty && !after.empty) {
1869         if(before.front < after.front) {
1870             less++;
1871         } else if(after.front < before.front) {
1872             greater++;
1873         }
1874 
1875         // Ignore equals.
1876         before.popFront;
1877         after.popFront;
1878     }
1879 
1880     double propGreater = to!double(greater) / (greater + less);
1881 
1882     final switch(alt) {
1883         case Alt.none:
1884             return TestRes(propGreater);
1885         case Alt.less:
1886             return TestRes(propGreater,
1887                 binomialCDF(greater, less + greater, 0.5));
1888         case Alt.greater:
1889             return TestRes(propGreater,
1890                 binomialCDF(less, less + greater, 0.5));
1891         case Alt.twoSided:
1892             if(less > greater) {
1893                 return TestRes(propGreater,
1894                     2 * binomialCDF(greater, less + greater, 0.5));
1895             } else if(greater > less) {
1896                 return  TestRes(propGreater,
1897                     2 * binomialCDF(less, less + greater, 0.5));
1898             } else {
1899                 return TestRes(propGreater, 1);
1900             }
1901     }
1902 }
1903 
1904 unittest {
1905     alias approxEqual ae;
1906     assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup), 1));
1907     assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup, Alt.less), 0.5));
1908     assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup, Alt.greater), 0.875));
1909     assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup, Alt.greater), 0.03125));
1910     assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup, Alt.less), 1));
1911     assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup), 0.0625));
1912 
1913     assert(approxEqual(signTest([1,2,6,7,9].dup, 2), 0.625));
1914     assert(ae(signTest([1,2,6,7,9].dup, 2).testStat, 0.75));
1915 }
1916 
1917 /**Similar to the overload, but allows testing for a difference between a
1918  * range and a fixed value mu.*/
1919 TestRes signTest(T)(T data, double mu, Alt alt = Alt.twoSided)
1920 if(doubleInput!(T) && is(typeof(data.front < mu) == bool)) {
1921     return signTest(data, repeat(mu), alt);
1922 }
1923 
1924 /**
1925 Two-sided binomial test for whether P(success) == p.  The one-sided
1926 alternatives are covered by dstats.distrib.binomialCDF and binomialCDFR.
1927 k is the number of successes observed, n is the number of trials, p
1928 is the probability of success under the null.
1929 
1930 Returns:  The P-value for the alternative that P(success) != p against
1931 the null that P(success) == p.
1932 
1933 Notes:  This test can also be performed using multinomialTest, but this
1934 implementation is much faster and easier to use.
1935  */
1936 double binomialTest(ulong k, ulong n, double p) {
1937     dstatsEnforce(k <= n, "k must be <= n for binomial test.");
1938     dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial test.");
1939 
1940     enum epsilon = 1 - 1e-6;  // Small but arbitrary constant to deal w/ rounding error.
1941 
1942     immutable mode = cast(long) ((n + 1) * p);
1943     if(k == mode ||
1944        approxEqual(binomialPMF(k, n, p), binomialPMF(mode, n, p), 1 - epsilon)) {
1945         return 1;
1946     } else if(k > mode) {
1947         immutable double upperPart = binomialCDFR(k, n, p);
1948         immutable pExact = binomialPMF(k, n, p);
1949         ulong ulim = mode, llim = 0, guess;
1950         while(ulim - llim > 1) {
1951             guess = (ulim + llim) / 2;
1952             immutable double pGuess = binomialPMF(guess, n, p);
1953 
1954             if(pGuess == pExact) {
1955                 ulim = guess + 1;
1956                 llim = guess;
1957                 break;
1958             } else if(pGuess < pExact) {
1959                 llim = guess;
1960             } else {
1961                 ulim = guess;
1962             }
1963         }
1964 
1965         guess = ulim;
1966         while(binomialPMF(guess, n, p) < pExact * epsilon) {
1967             guess++;
1968         }
1969         while(guess > 0 && binomialPMF(guess, n, p) > pExact / epsilon) {
1970             guess--;
1971         }
1972         if(guess == 0 && binomialPMF(0, n, p) > pExact / epsilon) {
1973             return upperPart;
1974         }
1975         return upperPart + binomialCDF(guess, n, p);
1976     } else {
1977         static double myPMF(ulong k, ulong n, double p) {
1978             return k > n ? 0 : binomialPMF(k, n, p);
1979         }
1980 
1981         immutable lowerPart = binomialCDF(k, n, p);
1982         immutable pExact = binomialPMF(k, n, p);
1983         ulong ulim = n + 1, llim = mode, guess;
1984         while(ulim - llim > 1) {
1985             guess = (ulim + llim) / 2;
1986             immutable double pGuess = myPMF(guess, n, p);
1987             if(pGuess == pExact) {
1988                 ulim = guess;
1989                 llim = guess;
1990                 break;
1991             } else if(pGuess < pExact) {
1992                 ulim = guess;
1993             } else {
1994                 llim = guess;
1995             }
1996         }
1997 
1998         // All this stuff is necessary to deal with round-off error properly.
1999         guess = llim;
2000         while(myPMF(guess, n, p) < pExact * epsilon && guess > 0) {
2001             guess--;
2002         }
2003         while(myPMF(guess, n, p) > pExact / epsilon) {
2004             guess++;
2005         }
2006 
2007         return lowerPart + ((guess > n) ? 0 : binomialCDFR(guess, n, p));
2008     }
2009 }
2010 
2011 unittest {
2012     // Values from R.
2013     assert(approxEqual(binomialTest(46, 96, 0.5), 0.759649));
2014     assert(approxEqual(binomialTest(44, 56, 0.5), 2.088e-5));
2015     assert(approxEqual(binomialTest(12, 56, 0.5), 2.088e-5));
2016     assert(approxEqual(binomialTest(0, 40, 0.25), 2.236e-5));
2017     assert(approxEqual(binomialTest(5, 16, 0.5), 0.2101));
2018     assert(approxEqual(binomialTest(0, 20, 0.4), 4.16e-5));
2019     assert(approxEqual(binomialTest(20, 20, 0.6), 4.16e-5));
2020     assert(approxEqual(binomialTest(6, 88, 0.1), 0.3784));
2021     assert(approxEqual(binomialTest(3, 4, 0.5), 0.625));
2022     assert(approxEqual(binomialTest(4, 7, 0.8), 0.1480));
2023     assert(approxEqual(binomialTest(3, 9, 0.8), 0.003066));
2024     assert(approxEqual(binomialTest(9, 9, 0.7), 0.06565));
2025     assert(approxEqual(binomialTest(2, 11, 0.1), 0.3026));
2026     assert(approxEqual(binomialTest(1, 11, 0.1), 1));
2027     assert(approxEqual(binomialTest(5, 11, 0.1), 0.002751));
2028     assert(approxEqual(binomialTest(5, 12, 0.5), 0.7744));
2029     assert(approxEqual(binomialTest(12, 12, 0.5), 0.0004883));
2030     assert(approxEqual(binomialTest(12, 13, 0.6), 0.02042));
2031     assert(approxEqual(binomialTest(0, 9, 0.1), 1));
2032 }
2033 
2034 ///For chiSquareFit and gTestFit, is expected value range counts or proportions?
2035 enum Expected {
2036     ///
2037     count,
2038 
2039     ///
2040     proportion
2041 }
2042 
2043 /**Performs a one-way Pearson's chi-square goodness of fit test between a range
2044 of observed and a range of expected values.  This is a useful statistical
2045 test for testing whether a set of observations fits a discrete distribution.
2046 
2047 Returns:  A TestRes of the chi-square statistic and the P-value for the
2048 alternative hypothesis that observed is not a sample from expected against
2049 the null that observed is a sample from expected.
2050 
2051 Notes:  By default, expected is assumed to be a range of expected proportions.
2052 These proportions are automatically normalized, and can sum to any number.
2053 By passing Expected.count in as the last parameter, calculating expected
2054 counts will be skipped, and expected will assume to already be properly
2055 normalized.  This is slightly faster, but more importantly
2056 allows input ranges to be used.
2057 
2058 The chi-square test relies on asymptotic statistical properties
2059 and is therefore not considered valid, as a rule of thumb,  when expected
2060 counts are below 5.  However, this rule is likely to be unnecessarily
2061 stringent in most cases.
2062 
2063 Examples:
2064 ---
2065 // Test to see whether a set of categorical observations differs
2066 // statistically from a discrete uniform distribution.
2067 
2068 uint[] observed = [980, 1028, 1001, 964, 1102];
2069 auto expected = repeat(1.0);
2070 auto res2 = chiSquareFit(observed, expected);
2071 assert(approxEqual(res2, 0.0207));
2072 assert(approxEqual(res2.testStat, 11.59));
2073 ---
2074  *
2075 References:  http://en.wikipedia.org/wiki/Pearson%27s_chi-square_test
2076  */
2077 TestRes chiSquareFit(T, U)(
2078     T observed,
2079     U expected,
2080     Expected countProp = Expected.proportion
2081 ) if(doubleInput!(T) && doubleInput!(U)) {
2082     return goodnessFit!(pearsonChiSqElem, T, U)(observed, expected, countProp);
2083 }
2084 
2085 unittest {
2086     // Test to see whether a set of categorical observations differs
2087     // statistically from a discrete uniform distribution.
2088     uint[] observed = [980, 1028, 1001, 964, 1102];
2089     auto expected = repeat(cast(double) sum(observed) / observed.length);
2090     auto res = chiSquareFit(observed, expected, Expected.count);
2091     assert(approxEqual(res, 0.0207));
2092     assert(approxEqual(res.testStat, 11.59));
2093 
2094     auto expected2 = [5.0, 5, 5, 5, 5, 0];
2095     observed ~= 0;
2096     auto res2 = chiSquareFit(observed, expected2);
2097     assert(approxEqual(res2, 0.0207));
2098     assert(approxEqual(res2.testStat, 11.59));
2099 }
2100 
2101 // Alias for old name, for backwards compatibility.  Don't document it
2102 // because it will be deprecated eventually.
2103 alias chiSquareFit chiSqrFit;
2104 
2105 /**The G or likelihood ratio chi-square test for goodness of fit.  Roughly
2106  * the same as Pearson's chi-square test (chiSquareFit), but may be more
2107  * accurate in certain situations and less accurate in others.  However, it is
2108  * still based on asymptotic distributions, and is not exact. Usage is is
2109  * identical to chiSquareFit.
2110  *
2111  * References:  http://en.wikipedia.org/wiki/G_test
2112  *
2113  */
2114 TestRes gTestFit(T, U)(T observed, U expected, Expected countProp = Expected.proportion)
2115 if(doubleInput!(T) && doubleInput!(U)) {
2116     return goodnessFit!(gTestElem, T, U)(observed, expected, countProp);
2117 }
2118 // No unittest because I can't find anything to test this against.  However,
2119 // it's hard to imagine how it could be wrong, given that goodnessFit() and
2120 // gTestElem() both work, and, as expected, this function produces roughly
2121 // the same results as chiSquareFit.
2122 
2123 private TestRes goodnessFit(alias elemFun, T, U)(T observed, U expected, Expected countProp)
2124 if(doubleInput!(T) && doubleInput!(U)) {
2125     if(countProp == Expected.proportion) {
2126         dstatsEnforce(isForwardRange!(U),
2127             "Can't use expected proportions instead of counts with input ranges.");
2128     }
2129 
2130     uint len = 0;
2131     double chiSq = 0;
2132     double multiplier = 1;
2133 
2134     // Normalize proportions to add up to the sum of the data.
2135     if(countProp == Expected.proportion) {
2136         double expectSum = 0;
2137         multiplier = 0;
2138         auto obsCopy = observed.save;
2139         auto expCopy = expected.save;
2140         while(!obsCopy.empty && !expCopy.empty) {
2141             multiplier += obsCopy.front;
2142             expectSum += expCopy.front;
2143             obsCopy.popFront;
2144             expCopy.popFront;
2145         }
2146         multiplier /= expectSum;
2147     }
2148 
2149     while(!observed.empty && !expected.empty) {
2150         scope(exit) {
2151             observed.popFront();
2152             expected.popFront();
2153         }
2154         double e = expected.front * multiplier;
2155 
2156         // If e is zero, then we should just treat the cell as if it didn't
2157         // exist.
2158         if(e == 0) {
2159             dstatsEnforce(observed.front == 0,
2160                 "Can't have non-zero observed value w/ zero expected value.");
2161             continue;
2162         }
2163 
2164         chiSq += elemFun(observed.front, e);
2165         len++;
2166     }
2167 
2168     if(isNaN(chiSq)) {
2169         return TestRes(double.nan, double.nan);
2170     }
2171 
2172     return TestRes(chiSq, chiSquareCDFR(chiSq, len - 1));
2173 }
2174 
2175 /**
2176 The exact multinomial goodness of fit test for whether a set of counts
2177 fits a hypothetical distribution.  counts is an input range of counts.
2178 proportions is an input range of expected proportions.  These are normalized
2179 automatically, so they can sum to any value.
2180 
2181 Returns:  The P-value for the null that counts is a sample from proportions
2182 against the alternative that it isn't.
2183 
2184 Notes:  This test is EXTREMELY slow for anything but very small samples and
2185 degrees of freedom.  The Pearson's chi-square (chiSquareFit()) or likelihood
2186 ratio chi-square (gTestFit()) are good enough approximations unless sample
2187 size is very small.
2188  */
2189 double multinomialTest(U, F)(U countsIn, F proportions)
2190 if(isInputRange!U && isInputRange!F &&
2191    isIntegral!(ElementType!U) && isFloatingPoint!(ElementType!(F))) {
2192     auto alloc = newRegionAllocator();
2193 
2194     static if(isRandomAccessRange!U && hasLength!U) {
2195         alias countsIn counts;
2196     } else {
2197         auto counts = alloc.array(countsIn);
2198     }
2199 
2200     uint N = sum(counts);
2201 
2202     double[] logPs;
2203     static if(std.range.hasLength!F) {
2204         logPs = alloc.uninitializedArray!(double[])(proportions.length);
2205         size_t pIndex;
2206         foreach(p; proportions) {
2207             logPs[pIndex++] = p;
2208         }
2209     } else {
2210         auto app = appender(logPs);
2211         foreach(p; proportions) {
2212             app.put(p);
2213         }
2214         logPs = app.data;
2215     }
2216 
2217     logPs[] /= reduce!"a + b"(0.0, logPs);
2218     foreach(ref elem; logPs) {
2219         elem = log(elem);
2220     }
2221 
2222 
2223     double[] logs = alloc.uninitializedArray!(double[])(N + 1);
2224     logs[0] = 0;
2225     foreach(i; 1..logs.length) {
2226         logs[i] = log(cast(double)i);
2227     }
2228 
2229     double nFact = logFactorial(N);
2230     double pVal = 0;
2231     uint nLeft = N;
2232     double pSoFar = nFact;
2233 
2234     double pActual = nFact;
2235     foreach(i, count; counts) {
2236         pActual += logPs[i] * count - logFactorial(count);
2237     }
2238     pActual -= pActual * 1e-6;  // Epsilon to handle numerical inaccuracy.
2239 
2240     void doIt(uint pos) {
2241         if(pos == counts.length - 1) {
2242             immutable pOld = pSoFar;
2243             pSoFar += logPs[$ - 1] * nLeft - logFactorial(nLeft);
2244 
2245             if(pSoFar <= pActual) {
2246                 pVal += exp(pSoFar);
2247             }
2248             pSoFar = pOld;
2249             return;
2250         }
2251 
2252         uint nLeftOld = nLeft;
2253         immutable pOld = pSoFar;
2254         double pAdd = 0;
2255 
2256         foreach(i; 0..nLeft + 1) {
2257             if(i > 0) {
2258                 pAdd += logPs[pos] - logs[i];
2259             }
2260             pSoFar = pOld + pAdd;
2261             doIt(pos + 1);
2262             nLeft--;
2263         }
2264         nLeft = nLeftOld;
2265         pSoFar = pOld;
2266     }
2267     doIt(0);
2268     return pVal;
2269 }
2270 
2271 unittest {
2272     // Nothing to test this against for more than 1 df, but it matches
2273     // chi-square roughly and should take the same paths for 2 vs. N degrees
2274     // of freedom.
2275     for(uint n = 4; n <= 100; n += 4) {
2276         foreach(k; 0..n + 1) {
2277             for(double p = 0.05; p <= 0.95; p += 0.05) {
2278                 double bino = binomialTest(k, n, p);
2279                 double[] ps = [p, 1 - p];
2280                 uint[] counts = [k, n - k];
2281                 double multino = multinomialTest(counts, ps);
2282                 //writeln(k, "\t", n, "\t", p, "\t", bino, "\t", multino);
2283                 assert(approxEqual(bino, multino),
2284                     text(bino, '\t', multino, '\t', k, '\t', n, '\t', p));
2285             }
2286         }
2287     }
2288 }
2289 
2290 /**
2291 Performs a Pearson's chi-square test on a contingency table of arbitrary
2292 dimensions.  When the chi-square test is mentioned, this is usually the one
2293 being referred to.  Takes a set of finite forward ranges, one for each column
2294 in the contingency table.  These can be expressed either as a tuple of ranges
2295 or a range of ranges.  Returns a P-value for the alternative hypothesis that
2296 frequencies in each row of the contingency table depend on the column against
2297 the null that they don't.
2298 
2299 Notes:  The chi-square test relies on asymptotic statistical properties
2300 and is therefore not exact.  The typical rule of thumb is that each cell
2301 should have an expected value of at least 5.  However, this is likely to
2302 be unnecessarily stringent.
2303 
2304 Yates's continuity correction is never used in this implementation.  If
2305 you want something that's guaranteed to be conservative, use fisherExact().
2306 
2307 This is, for all practical purposes, an inherently non-directional test.
2308 Therefore, the one-sided verses two-sided option is not provided.
2309 
2310 For 2x2 contingency tables, fisherExact is a more conservative test, in that
2311 the type I error rate is guaranteed to never be above the nominal P-value.
2312 However, even for small sample sizes this test may produce results closer
2313 to the true P-value, at the risk of possibly being non-conservative.
2314 
2315 Examples:
2316 ---
2317 // Test to see whether the relative frequency of outcome 0, 1, and 2
2318 // depends on the treatment (drug1, drug2 or placebo) in some hypothetical
2319 // experiment.  For example, 1500 people had outcome 2 if they were treated
2320 // with drug1 and 1100 had outcome 1 if treated with placebo.
2321 uint[] drug1 = [1000, 2000, 1500];
2322 uint[] drug2 = [1500, 3000, 2300];
2323 uint[] placebo = [500, 1100, 750];
2324 auto result1 = chiSquareContingency(drug1, drug2, placebo);
2325 
2326 // The following uses a range of ranges instead of an array of ranges,
2327 // and will produce equivalent results.
2328 auto rangeOfRanges = [drug1, drug2, placebo];
2329 auto result2 = chiSquareContingency(rangeOfRanges);
2330 ---
2331 
2332 References: http://en.wikipedia.org/wiki/Pearson%27s_chi-square_test
2333  *
2334  */
2335 TestRes chiSquareContingency(T...)(T inputData) {
2336     return testContingency!(pearsonChiSqElem, T)(inputData);
2337 }
2338 
2339 unittest {
2340     // Test array version.  Using VassarStat's chi-square calculator.
2341     uint[][] table1 = [[60, 80, 70],
2342                        [20, 50, 40],
2343                        [10, 15, 11]];
2344     uint[][] table2 = [[60, 20, 10],
2345                        [80, 50, 15],
2346                        [70, 40, 11]];
2347     assert(approxEqual(chiSquareContingency(table1), 0.3449));
2348     assert(approxEqual(chiSquareContingency(table2), 0.3449));
2349     assert(approxEqual(chiSquareContingency(table1).testStat, 4.48));
2350 
2351     // Test tuple version.
2352     auto p1 = chiSquareContingency(cast(uint[]) [31, 41, 59],
2353                                 cast(uint[]) [26, 53, 58],
2354                                 cast(uint[]) [97, 93, 93]);
2355     assert(approxEqual(p1, 0.0059));
2356 
2357     auto p2 = chiSquareContingency(cast(uint[]) [31, 26, 97],
2358                                 cast(uint[]) [41, 53, 93],
2359                                 cast(uint[]) [59, 58, 93]);
2360     assert(approxEqual(p2, 0.0059));
2361 
2362     uint[] drug1 = [1000, 2000, 1500];
2363     uint[] drug2 = [1500, 3000, 2300];
2364     uint[] placebo = [500, 1100, 750];
2365     assert(approxEqual(chiSquareContingency(drug1, drug2, placebo), 0.2397));
2366 }
2367 
2368 // Alias for old name, for backwards compatibility.  Don't document it
2369 // because it is deprecated and has been scheduled for deprecation for
2370 // ages.
2371 deprecated alias chiSquareContingency chiSqrContingency;
2372 
2373 /**
2374 This struct is a subtype of TestRes and is used to return the results of
2375 gTestContingency and gTestObs.  Due to the information theoretic interpretation
2376 of the G test, it contains an extra field to return the mutual information
2377 in bits.
2378 */
2379 struct GTestRes {
2380     ///
2381     TestRes testRes;
2382 
2383     ///
2384     alias testRes this;
2385 
2386     /**
2387     The mutual info of the two random variables in the joint distribution
2388     represented by the contingency table, in bits (base 2).
2389     */
2390     double mutualInfo;
2391 }
2392 
2393 /**
2394 The G or likelihood ratio chi-square test for contingency tables.  Roughly
2395 the same as Pearson's chi-square test (chiSquareContingency), but may be more
2396 accurate in certain situations and less accurate in others.
2397 
2398 Like Pearson's Chi-square, the G-test is based on asymptotic distributions,
2399 and is not exact. Usage is is identical to chiSquareContingency.
2400 
2401 Note:  This test can be thought of as a test for nonzero mutual information
2402 between the random variables represented by the rows and the columns,
2403 since the test statistic and P-value are strictly increasing
2404 and strictly decreasing, respectively, in mutual information.  Therefore, this
2405 function returns a GTestRes, which is a subtype of TestRes and also gives
2406 the mutual information for use in information theoretic settings.
2407 
2408 References:  http://en.wikipedia.org/wiki/G_test, last retrieved 1/22/2011
2409 */
2410 GTestRes gTestContingency(T...)(T inputData) {
2411     return testContingency!(gTestElem, T)(inputData);
2412 }
2413 
2414 unittest {
2415     // Values from example at http://udel.edu/~mcdonald/statgtestind.html
2416     // Handbook of Biological Statistics.
2417     uint[] withoutCHD = [268, 199, 42];
2418     uint[] withCHD = [807, 759, 184];
2419     auto res = gTestContingency(withoutCHD, withCHD);
2420     assert(approxEqual(res.testStat, 7.3));
2421     assert(approxEqual(res.p, 0.026));
2422     assert(approxEqual(res.mutualInfo, 0.0023313));
2423 
2424 
2425     uint[] moringa = [127, 99, 264];
2426     uint[] vicinus = [116, 67, 161];
2427     auto res2 = gTestContingency(moringa, vicinus);
2428     assert(approxEqual(res2.testStat, 6.23));
2429     assert(approxEqual(res2.p, 0.044));
2430     assert(approxEqual(res2.mutualInfo, 0.00538613));
2431 }
2432 
2433 // Pearson and likelihood ratio code are pretty much the same.  Factor out
2434 // the one difference into a function that's a template parameter.  However,
2435 // for API simplicity, this is hidden and they look like two separate functions.
2436 private GTestRes testContingency(alias elemFun, T...)(T rangesIn) {
2437     auto alloc = newRegionAllocator();
2438     static if(isInputRange!(T[0]) && T.length == 1 &&
2439         isForwardRange!(typeof(rangesIn[0].front()))) {
2440         auto ranges = alloc.array(rangesIn[0]);
2441 
2442         foreach(ref range; ranges) {
2443             range = range.save;
2444         }
2445 
2446     } else static if(allSatisfy!(isForwardRange, typeof(rangesIn))) {
2447         auto saved = saveAll(rangesIn);
2448         auto ranges = saved.expand;
2449     } else {
2450         static assert(0, "Can only perform contingency table test" ~
2451             " on a tuple of ranges or a range of ranges.");
2452     }
2453 
2454     double[] colSums = alloc.uninitializedArray!(double[])(ranges.length);
2455     colSums[] = 0;
2456     size_t nCols = 0;
2457     immutable size_t nRows = ranges.length;
2458     foreach(ri, range; ranges) {
2459         size_t curLen = 0;
2460         foreach(elem; range.save) {
2461             colSums[ri] += cast(double) elem;
2462             curLen++;
2463         }
2464         if(ri == 0) {
2465             nCols = curLen;
2466         } else {
2467             assert(curLen == nCols);
2468         }
2469     }
2470 
2471     bool noneEmpty() {
2472         foreach(range; ranges) {
2473             if(range.empty) {
2474                 return false;
2475             }
2476         }
2477         return true;
2478     }
2479 
2480     void popAll() {
2481         foreach(i, range; ranges) {
2482             ranges[i].popFront;
2483         }
2484     }
2485 
2486     double sumRow() {
2487         double rowSum = 0;
2488         foreach(range; ranges) {
2489             rowSum += cast(double) range.front;
2490         }
2491         return rowSum;
2492     }
2493 
2494     double chiSq = 0;
2495     immutable double NNeg1 = 1.0 / sum(colSums);
2496     while(noneEmpty) {
2497         auto rowSum = sumRow();
2498         foreach(ri, range; ranges) {
2499             double expected = NNeg1 * rowSum * colSums[ri];
2500             chiSq += elemFun(range.front, expected);
2501         }
2502         popAll();
2503     }
2504 
2505     if(isNaN(chiSq)) {
2506         return GTestRes(TestRes(double.nan, double.nan), double.nan);
2507     }
2508 
2509     // This can happen in some cases due to numerical fuzz.
2510     if(chiSq > 1e-5 && chiSq <= 0) {
2511         return GTestRes(TestRes(0, 1), 0);
2512     }
2513 
2514     immutable pVal = (chiSq >= 0) ?
2515         chiSquareCDFR(chiSq, (nRows - 1) * (nCols - 1)) : double.nan;
2516 
2517     // 1 / (2 * LN2), for converting chiSq to mutualInfo.
2518     enum chiToMi = 1 / (2 * LN2);
2519 
2520     // This is the mutual information between the two random variables
2521     // represented by the contingency table, only if we're doing a G test.
2522     // If we're doing a Pearson's test, it's a completely meaningless quantity,
2523     // but never gets returned by any public function.
2524     immutable mutualInfo = chiSq * NNeg1 * chiToMi;
2525 
2526     return GTestRes(TestRes(chiSq, pVal), mutualInfo);
2527 }
2528 
2529 private double pearsonChiSqElem(double observed, double expected) pure nothrow {
2530     immutable diff = observed - expected;
2531     return diff * diff / expected;
2532 }
2533 
2534 private double gTestElem(double observed, double expected) pure nothrow {
2535     return (observed == 0 && expected > 0) ? 0 :
2536         (observed * log(observed / expected) * 2);
2537 }
2538 
2539 /**
2540 Given two vectors of observations of jointly distributed variables x, y, tests
2541 the null hypothesis that values in x are independent of the corresponding
2542 values in y.  This is done using Pearson's Chi-Square Test.  For a similar test
2543 that assumes the data has already been tabulated into a contingency table, see
2544 chiSquareContingency.
2545 
2546 x and y must both be input ranges.  If they are not the same length, a
2547 DstatsArgumentException is thrown.
2548 
2549 Examples:
2550 ---
2551 // Test whether the appearance of "foo" vs. "bar" is independent of the
2552 // appearance of "baz" vs. "xxx".
2553 auto x = ["foo", "bar", "bar", "foo", "foo"];
2554 auto y = ["xxx", "baz", "baz", "xxx", "baz"];
2555 auto result = chiSquareObs(x, y);
2556 
2557 // This is equivalent to:
2558 auto contingency = new uint[][](2, 2);
2559 foreach(i; 0..x.length) {
2560     immutable index1 = (x[i] == "foo");
2561     immutable index2 = (y[i] == "xxx");
2562     contingency[index1][index2]++;
2563 }
2564 
2565 auto result2 = chiSquareContingency(contingency);
2566 ---
2567 */
2568 TestRes chiSquareObs(T, U)(T x, U y)
2569 if(isInputRange!T && isInputRange!U) {
2570     uint xFreedom, yFreedom, n;
2571     typeof(return) ret;
2572 
2573     static if(!hasLength!T && !hasLength!U) {
2574         ret.testStat = toContingencyScore!(T, U, uint)
2575             (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
2576     } else {
2577         immutable minLen = min(x.length, y.length);
2578         if(minLen <= ubyte.max) {
2579             ret.testStat = toContingencyScore!(T, U, ubyte)
2580                 (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
2581         } else if(minLen <= ushort.max) {
2582             ret.testStat = toContingencyScore!(T, U, ushort)
2583                 (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
2584         } else {
2585             ret.testStat = toContingencyScore!(T, U, uint)
2586                 (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
2587         }
2588     }
2589 
2590     ret.p = chiSquareCDFR(ret.testStat, xFreedom * yFreedom);
2591     return ret;
2592 }
2593 
2594 unittest {
2595     // We know the chi-square contingency works, so test that the automatic
2596     // binning works, too.
2597     ubyte[] obs1 = [1, 2, 3, 1, 2, 3, 1, 2, 3];
2598     ubyte[] obs2 = [1, 3, 2, 1, 3, 2, 1, 3, 2];
2599 
2600     uint[][] cTable = [[3, 0, 0],
2601                        [0, 0, 3],
2602                        [0, 3, 0]];
2603     auto gRes = chiSquareContingency(cTable);
2604     auto miRes = chiSquareObs(obs1, obs2);
2605     foreach(ti, elem; miRes.tupleof) {
2606         assert(approxEqual(elem, gRes.tupleof[ti]));
2607     }
2608 
2609     auto x = ["foo", "bar", "bar", "foo", "foo"];
2610     auto y = ["xxx", "baz", "baz", "xxx", "baz"];
2611     auto result = chiSquareObs(x, y);
2612     assert(approxEqual(result.testStat, 2.22222222));
2613     assert(approxEqual(result.p, 0.136037));
2614 
2615     auto contingency = new uint[][](2, 2);
2616     foreach(i; 0..x.length) {
2617         immutable index1 = (x[i] == "foo");
2618         immutable index2 = (y[i] == "xxx");
2619         contingency[index1][index2]++;
2620     }
2621 
2622     auto result2 = chiSquareContingency(contingency);
2623     assert(approxEqual(result.testStat, result2.testStat),
2624         text(result.testStat, ' ', result2.testStat));
2625     assert(approxEqual(result.p, result2.p));
2626 }
2627 
2628 /**
2629 Given two ranges of observations of jointly distributed variables x, y, tests
2630 the null hypothesis that values in x are independent of the corresponding
2631 values in y.  This is done using the Likelihood Ratio G test.  Usage is similar
2632 to chiSquareObs.  For an otherwise identical test that assumes the data has
2633 already been tabulated into a contingency table, see gTestContingency.
2634 
2635 Note:  This test can be thought of as a test for nonzero mutual information
2636 between x and y, since the test statistic and P-value are strictly increasing
2637 and strictly decreasing, respectively, in mutual information.  Therefore, this
2638 function returns a GTestRes, which is a subtype of TestRes and also gives
2639 the mutual information for use in information theoretic settings.
2640 */
2641 GTestRes gTestObs(T, U)(T x, U y)
2642 if(isInputRange!T && isInputRange!U) {
2643     uint xFreedom, yFreedom, n;
2644     typeof(return) ret;
2645 
2646     static if(!hasLength!T && !hasLength!U) {
2647         ret.testStat = toContingencyScore!(T, U, uint)
2648             (x, y, &gTestElem, xFreedom, yFreedom, n);
2649     } else {
2650         immutable minLen = min(x.length, y.length);
2651         if(minLen <= ubyte.max) {
2652             ret.testStat = toContingencyScore!(T, U, ubyte)
2653                 (x, y, &gTestElem, xFreedom, yFreedom, n);
2654         } else if(minLen <= ushort.max) {
2655             ret.testStat = toContingencyScore!(T, U, ushort)
2656                 (x, y, &gTestElem, xFreedom, yFreedom, n);
2657         } else {
2658             ret.testStat = toContingencyScore!(T, U, uint)
2659                 (x, y, &gTestElem, xFreedom, yFreedom, n);
2660         }
2661     }
2662 
2663     ret.p = chiSquareCDFR(ret.testStat, xFreedom * yFreedom);
2664     ret.mutualInfo = ret.testStat / (2 * LN2 * n);
2665     return ret;
2666 }
2667 
2668 unittest {
2669     // We know the g test contingency works, so test that the automatic binning
2670     // works, too.
2671     ubyte[] obs1 = [1, 2, 3, 1, 2, 3, 1, 2, 3];
2672     ubyte[] obs2 = [1, 3, 2, 1, 3, 2, 1, 3, 2];
2673 
2674     uint[][] cTable = [[3, 0, 0],
2675                        [0, 0, 3],
2676                        [0, 3, 0]];
2677     auto gRes = gTestContingency(cTable);
2678     auto miRes = gTestObs(obs1, obs2);
2679     foreach(ti, elem; miRes.tupleof) {
2680         assert(approxEqual(elem, gRes.tupleof[ti]));
2681     }
2682 
2683     auto x = ["foo", "bar", "bar", "foo", "foo"];
2684     auto y = ["xxx", "baz", "baz", "xxx", "baz"];
2685     auto result = gTestObs(x, y);
2686     assert(approxEqual(result.testStat, 2.91103));
2687     assert(approxEqual(result.p, 0.0879755));
2688     assert(approxEqual(result.mutualInfo, 0.419973));
2689 }
2690 
2691 package double toContingencyScore(T, U, Uint)
2692 (T x, U y, double function(double, double) elemFun,
2693  out uint xFreedom, out uint yFreedom, out uint nPtr) {
2694 
2695     import dstats.infotheory : NeedsHeap;
2696     enum needsHeap = NeedsHeap!T ||
2697         NeedsHeap!U;
2698     alias dstats.infotheory.ObsEnt!(ElementType!T, ElementType!U) ObsType;
2699 
2700     static if(needsHeap) {
2701         Uint[ObsType] jointCounts;
2702         Uint[ElementType!T] xCounts;
2703         Uint[ElementType!U] yCounts;
2704     } else {
2705         auto alloc = newRegionAllocator();
2706         dstatsEnforce(x.length == y.length,
2707             "Can't calculate mutual info with different length vectors.");
2708         immutable len = x.length;
2709         auto jointCounts = StackHash!(ObsType, Uint)(max(20, len / 20), alloc);
2710         auto xCounts = StackHash!(ElementType!T, Uint)(max(10, len / 40), alloc);
2711         auto yCounts = StackHash!(ElementType!U, Uint)(max(10, len / 40), alloc);
2712     }
2713 
2714     uint n = 0;
2715     while(!x.empty && !y.empty) {
2716         n++;
2717         auto a = x.front;
2718         auto b = y.front;
2719         jointCounts[ObsType(a, b)]++;
2720         xCounts[a]++;
2721         yCounts[b]++;
2722 
2723         x.popFront();
2724         y.popFront();
2725     }
2726 
2727     dstatsEnforce(x.empty && y.empty,
2728         "Can't calculate mutual info with different length vectors.");
2729 
2730     xFreedom = cast(uint) xCounts.length - 1;
2731     yFreedom = cast(uint) yCounts.length - 1;
2732     nPtr = n;
2733 
2734     double ret = 0;
2735     immutable double nNeg1 = 1.0 / n;
2736     foreach(key1, marg1; xCounts) foreach(key2, marg2; yCounts) {
2737         immutable observed = jointCounts.get(
2738             ObsType(key1, key2), 0
2739         );
2740         immutable expected = marg1 * nNeg1 * marg2;
2741         ret += elemFun(observed, expected);
2742     }
2743 
2744     return ret;
2745 }
2746 
2747 /**
2748 Fisher's Exact test for difference in odds between rows/columns
2749 in a 2x2 contingency table.  Specifically, this function tests the odds
2750 ratio, which is defined, for a contingency table c, as (c[0][0] * c[1][1])
2751  / (c[1][0] * c[0][1]).  Alternatives are Alt.less, meaning true odds ratio
2752 < 1, Alt.greater, meaning true odds ratio > 1, and Alt.twoSided, meaning
2753 true odds ratio != 1.
2754 
2755 Accepts a 2x2 contingency table as an array of arrays of uints.
2756 For now, only does 2x2 contingency tables.
2757 
2758 Notes:  Although this test is "exact" in that it does not rely on asymptotic
2759 approximations, it is very statistically conservative when the marginals
2760 are not truly fixed in the experimental design in question.  If a
2761 closer but possibly non-conservative approximation of the true P-value is
2762 desired, Pearson's chi-square test (chiSquareContingency) may perform better,
2763 even for small samples.
2764 
2765 Returns:  A TestRes of the odds ratio and the P-value against the given
2766 alternative.
2767 
2768 Examples:
2769 ---
2770 double res = fisherExact([[2u, 7], [8, 2]], Alt.less);
2771 assert(approxEqual(res.p, 0.01852));  // Odds ratio is very small in this case.
2772 assert(approxEqual(res.testStat, 4.0 / 56.0));
2773 ---
2774 
2775 References:  http://en.wikipedia.org/wiki/Fisher%27s_Exact_Test
2776 */
2777 TestRes fisherExact(T)(const T[2][2] contingencyTable, Alt alt = Alt.twoSided)
2778 if(isIntegral!(T)) {
2779     foreach(range; contingencyTable) {
2780         foreach(elem; range) {
2781             dstatsEnforce(elem >= 0,
2782                 "Cannot have negative elements in a contingency table.");
2783         }
2784     }
2785 
2786     static double fisherLower(const T[2][2] contingencyTable) {
2787         alias contingencyTable c;
2788         return hypergeometricCDF(c[0][0], c[0][0] + c[0][1], c[1][0] + c[1][1],
2789                                  c[0][0] + c[1][0]);
2790     }
2791 
2792     static double fisherUpper(const T[2][2] contingencyTable) {
2793         alias contingencyTable c;
2794         return hypergeometricCDFR(c[0][0], c[0][0] + c[0][1], c[1][0] + c[1][1],
2795                                  c[0][0] + c[1][0]);
2796     }
2797 
2798 
2799     alias contingencyTable c;  // Save typing.
2800     immutable oddsRatio = cast(double) c[0][0] * c[1][1] / c[0][1] / c[1][0];
2801     if(alt == Alt.none) {
2802         return TestRes(oddsRatio);
2803     } else if(alt == Alt.less) {
2804         return TestRes(oddsRatio, fisherLower(contingencyTable));
2805     } else if(alt == Alt.greater) {
2806         return TestRes(oddsRatio, fisherUpper(contingencyTable));
2807     }
2808 
2809 
2810     immutable uint n1 = c[0][0] + c[0][1],
2811                    n2 = c[1][0] + c[1][1],
2812                    n  = c[0][0] + c[1][0];
2813 
2814     immutable uint mode =
2815         cast(uint) ((cast(double) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
2816     immutable double pExact = hypergeometricPMF(c[0][0], n1, n2, n);
2817     immutable double pMode = hypergeometricPMF(mode, n1, n2, n);
2818 
2819     enum epsilon = 1 - 1e-5;
2820     if(approxEqual(pExact, pMode, 1 - epsilon)) {
2821         return TestRes(oddsRatio, 1);
2822     } else if(c[0][0] < mode) {
2823         immutable double pLower = hypergeometricCDF(c[0][0], n1, n2, n);
2824 
2825         if(hypergeometricPMF(n, n1, n2, n) > pExact / epsilon) {
2826             return TestRes(oddsRatio, pLower);
2827         }
2828 
2829         // Binary search for where to begin upper half.
2830         uint min = mode, max = n, guess = uint.max;
2831         while(max - min > 1) {
2832             guess = cast(uint) (
2833                     (max == min + 1 && guess == min) ? max :
2834                     (cast(ulong) max + cast(ulong) min) / 2UL);
2835 
2836             immutable double pGuess = hypergeometricPMF(guess, n1, n2, n);
2837             if(pGuess <= pExact &&
2838                 hypergeometricPMF(guess - 1, n1, n2, n) > pExact) {
2839                 break;
2840             } else if(pGuess < pExact) {
2841                 max = guess;
2842             } else min = guess;
2843         }
2844 
2845         if(guess == uint.max) {
2846             guess = min;
2847         }
2848 
2849         while(guess > 0 && hypergeometricPMF(guess, n1, n2, n) < pExact * epsilon) {
2850             guess--;
2851         }
2852 
2853         while(hypergeometricPMF(guess, n1, n2, n) > pExact / epsilon) {
2854             guess++;
2855         }
2856 
2857         static import std.algorithm;
2858         double p = std.algorithm.min(pLower +
2859                hypergeometricCDFR(guess, n1, n2, n), 1.0);
2860         return TestRes(oddsRatio, p);
2861     } else {
2862         immutable double pUpper = hypergeometricCDFR(c[0][0], n1, n2, n);
2863 
2864         if(hypergeometricPMF(0, n1, n2, n) > pExact / epsilon) {
2865             return TestRes(oddsRatio, pUpper);
2866         }
2867 
2868         // Binary search for where to begin lower half.
2869         uint min = 0, max = mode, guess = uint.max;
2870         while(max - min > 1) {
2871             guess = cast(uint) (
2872                     (max == min + 1 && guess == min) ? max :
2873                     (cast(ulong) max + cast(ulong) min) / 2UL);
2874             immutable double pGuess = hypergeometricPMF(guess, n1, n2, n);
2875 
2876             if(pGuess <= pExact &&
2877                 hypergeometricPMF(guess + 1, n1, n2, n) > pExact) {
2878                 break;
2879             } else if(pGuess <= pExact) {
2880                 min = guess;
2881             } else max = guess;
2882         }
2883 
2884         if(guess == uint.max) {
2885             guess = min;
2886         }
2887 
2888         while(hypergeometricPMF(guess, n1, n2, n) < pExact * epsilon) {
2889             guess++;
2890         }
2891 
2892         while(guess > 0 &&
2893             hypergeometricPMF(guess, n1, n2, n) > pExact / epsilon) {
2894             guess--;
2895         }
2896 
2897         static import std.algorithm;
2898         double p = std.algorithm.min(pUpper +
2899                hypergeometricCDF(guess, n1, n2, n), 1.0);
2900         return TestRes(oddsRatio, p);
2901     }
2902 }
2903 
2904 /**
2905 Convenience function.  Converts a dynamic array to a static one, then
2906 calls the overload.
2907 */
2908 TestRes fisherExact(T)(const T[][] contingencyTable, Alt alt = Alt.twoSided)
2909 if(isIntegral!(T)) {
2910     dstatsEnforce(contingencyTable.length == 2 &&
2911             contingencyTable[0].length == 2 &&
2912             contingencyTable[1].length == 2,
2913             "Fisher exact only supports 2x2 tables.");
2914 
2915     T[2][2] newTable;
2916     newTable[0][0] = contingencyTable[0][0];
2917     newTable[0][1] = contingencyTable[0][1];
2918     newTable[1][1] = contingencyTable[1][1];
2919     newTable[1][0] = contingencyTable[1][0];
2920     return fisherExact(newTable, alt);
2921 }
2922 
2923 unittest {
2924     // Simple, naive impl. of two-sided to test against.
2925     static double naive(const uint[][] c) {
2926         immutable uint n1 = c[0][0] + c[0][1],
2927                    n2 = c[1][0] + c[1][1],
2928                    n  = c[0][0] + c[1][0];
2929         immutable uint mode =
2930             cast(uint) ((cast(double) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
2931         immutable double pExact = hypergeometricPMF(c[0][0], n1, n2, n);
2932         immutable double pMode = hypergeometricPMF(mode, n1, n2, n);
2933         if(approxEqual(pExact, pMode, 1e-7))
2934             return 1;
2935         double sum = 0;
2936         foreach(i; 0..n + 1) {
2937             double pCur = hypergeometricPMF(i, n1, n2, n);
2938             if(pCur <= pExact / (1 - 1e-5))
2939                 sum += pCur;
2940         }
2941         return sum;
2942     }
2943 
2944     uint[][] c = new uint[][](2, 2);
2945 
2946     foreach(i; 0..100_000) {
2947         c[0][0] = uniform(0U, 51U);
2948         c[0][1] = uniform(0U, 51U);
2949         c[1][0] = uniform(0U, 51U);
2950         c[1][1] = uniform(0U, 51U);
2951         double naiveAns = naive(c);
2952         double fastAns = fisherExact(c);
2953         assert(approxEqual(naiveAns, fastAns), text(c, " - ", naiveAns, " - ", fastAns));
2954     }
2955 
2956     auto res = fisherExact([[19000, 80000], [20000, 90000]]);
2957     assert(approxEqual(res.testStat, 1.068731));
2958     assert(approxEqual(res, 3.319e-9));
2959     res = fisherExact([[18000, 80000], [20000, 90000]]);
2960     assert(approxEqual(res, 0.2751));
2961     res = fisherExact([[14500, 20000], [30000, 40000]]);
2962     assert(approxEqual(res, 0.01106));
2963     res = fisherExact([[100, 2], [1000, 5]]);
2964     assert(approxEqual(res, 0.1301));
2965     res = fisherExact([[2, 7], [8, 2]]);
2966     assert(approxEqual(res, 0.0230141));
2967     res = fisherExact([[5, 1], [10, 10]]);
2968     assert(approxEqual(res, 0.1973244));
2969     res = fisherExact([[5, 15], [20, 20]]);
2970     assert(approxEqual(res, 0.0958044));
2971     res = fisherExact([[5, 16], [20, 25]]);
2972     assert(approxEqual(res, 0.1725862));
2973     res = fisherExact([[10, 5], [10, 1]]);
2974     assert(approxEqual(res, 0.1973244));
2975     res = fisherExact([[5, 0], [1, 4]]);
2976     assert(approxEqual(res.p, 0.04761904));
2977     res = fisherExact([[0, 1], [3, 2]]);
2978     assert(approxEqual(res.p, 1.0));
2979     res = fisherExact([[0, 2], [6, 4]]);
2980     assert(approxEqual(res.p, 0.4545454545));
2981     res = fisherExact([[2, 7], [8, 2]], Alt.less);
2982     assert(approxEqual(res, 0.01852));
2983     res = fisherExact([[5, 1], [10, 10]], Alt.less);
2984     assert(approxEqual(res, 0.9783));
2985     res = fisherExact([[5, 15], [20, 20]], Alt.less);
2986     assert(approxEqual(res, 0.05626));
2987     res = fisherExact([[5, 16], [20, 25]], Alt.less);
2988     assert(approxEqual(res, 0.08914));
2989     res = fisherExact([[2, 7], [8, 2]], Alt.greater);
2990     assert(approxEqual(res, 0.999));
2991     res = fisherExact([[5, 1], [10, 10]], Alt.greater);
2992     assert(approxEqual(res, 0.1652));
2993     res = fisherExact([[5, 15], [20, 20]], Alt.greater);
2994     assert(approxEqual(res, 0.985));
2995     res = fisherExact([[5, 16], [20, 25]], Alt.greater);
2996     assert(approxEqual(res, 0.9723));
2997 }
2998 
2999 /**
3000 Performs a Kolmogorov-Smirnov (K-S) 2-sample test.  The K-S test is a
3001 non-parametric test for a difference between two empirical distributions or
3002 between an empirical distribution and a reference distribution.
3003 
3004 Returns:  A TestRes with the K-S D value and a P value for the null that
3005 FPrime is distributed identically to F against the alternative that it isn't.
3006 This implementation uses a signed D value to indicate the direction of the
3007 difference between distributions.  To get the D value used in standard
3008 notation, simply take the absolute value of this D value.
3009 
3010 Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.
3011 
3012 References:  http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
3013  */
3014 TestRes ksTest(T, U)(T F, U Fprime)
3015 if(doubleInput!(T) && doubleInput!(U)) {
3016     double D = ksTestD(F, Fprime);
3017     return TestRes(D, ksPval(F.length, Fprime.length, D));
3018 }
3019 
3020 unittest {
3021     assert(approxEqual(ksTest([1,2,3,4,5], [1,2,3,4,5]).testStat, 0));
3022     assert(approxEqual(ksTestDestructive([1,2,3,4,5], [1,2,2,3,5]).testStat, -.2));
3023     assert(approxEqual(ksTest([-1,0,2,8, 6], [1,2,2,3,5]).testStat, .4));
3024     assert(approxEqual(ksTest([1,2,3,4,5], [1,2,2,3,5,7,8]).testStat, .2857));
3025     assert(approxEqual(ksTestDestructive([1, 2, 3, 4, 4, 4, 5],
3026            [1, 2, 3, 4, 5, 5, 5]).testStat, .2857));
3027 
3028     assert(approxEqual(ksTest([1, 2, 3, 4, 4, 4, 5], [1, 2, 3, 4, 5, 5, 5]),
3029            .9375));
3030     assert(approxEqual(ksTestDestructive([1, 2, 3, 4, 4, 4, 5],
3031         [1, 2, 3, 4, 5, 5, 5]), .9375));
3032 }
3033 
3034 template isArrayLike(T) {
3035     enum bool isArrayLike = hasSwappableElements!(T) && hasAssignableElements!(T)
3036         && hasLength!(T) && isRandomAccessRange!(T);
3037 }
3038 
3039 /**
3040 One-sample Kolmogorov-Smirnov test against a reference distribution.
3041 Takes a callable object for the CDF of refernce distribution.
3042 
3043 Returns:  A TestRes with the Kolmogorov-Smirnov D value and a P value for the
3044 null that Femp is a sample from F against the alternative that it isn't. This
3045 implementation uses a signed D value to indicate the direction of the
3046 difference between distributions.  To get the D value used in standard
3047 notation, simply take the absolute value of this D value.
3048 
3049 Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.
3050 
3051 Examples:
3052 ---
3053 auto stdNormal = parametrize!(normalCDF)(0.0, 1.0);
3054 auto empirical = [1, 2, 3, 4, 5];
3055 auto res = ksTest(empirical, stdNormal);
3056 ---
3057 
3058 References:  http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
3059  */
3060 TestRes ksTest(T, Func)(T Femp, Func F)
3061 if(doubleInput!(T) && is(ReturnType!(Func) : double)) {
3062     double D = ksTestD(Femp, F);
3063     return TestRes(D, ksPval(Femp.length, D));
3064 }
3065 
3066 unittest {
3067     auto stdNormal = paramFunctor!(normalCDF)(0.0, 1.0);
3068     assert(approxEqual(ksTest([1,2,3,4,5].dup, stdNormal).testStat, -.8413));
3069     assert(approxEqual(ksTestDestructive([-1,0,2,8, 6].dup, stdNormal).testStat, -.5772));
3070     auto lotsOfTies = [5,1,2,2,2,2,2,2,3,4].dup;
3071     assert(approxEqual(ksTest(lotsOfTies, stdNormal).testStat, -0.8772));
3072 
3073     assert(approxEqual(ksTest([0,1,2,3,4].dup, stdNormal), .03271));
3074 
3075     auto uniform01 = parametrize!(uniformCDF)(0, 1);
3076     assert(approxEqual(ksTestDestructive([0.1, 0.3, 0.5, 0.9, 1].dup, uniform01), 0.7591));
3077 
3078 }
3079 
3080 /**Same as ksTest, except sorts in place, avoiding memory allocations.*/
3081 TestRes ksTestDestructive(T, U)(T F, U Fprime)
3082 if(isArrayLike!(T) && isArrayLike!(U)) {
3083     double D = ksTestDDestructive(F, Fprime);
3084     return TestRes(D, ksPval(F.length, Fprime.length, D));
3085 }
3086 
3087 ///Ditto.
3088 TestRes ksTestDestructive(T, Func)(T Femp, Func F)
3089 if(isArrayLike!(T) && is(ReturnType!Func : double)) {
3090     double D =  ksTestDDestructive(Femp, F);
3091     return TestRes(D, ksPval(Femp.length, D));
3092 }
3093 
3094 private double ksTestD(T, U)(T F, U Fprime)
3095 if(isInputRange!(T) && isInputRange!(U)) {
3096     auto alloc = newRegionAllocator();
3097     return ksTestDDestructive(alloc.array(F), alloc.array(Fprime));
3098 }
3099 
3100 private double ksTestDDestructive(T, U)(T F, U Fprime)
3101 if(isArrayLike!(T) && isArrayLike!(U)) {
3102     qsort(F);
3103     qsort(Fprime);
3104     double D = 0;
3105     size_t FprimePos = 0;
3106     foreach(i; 0..2) {  //Test both w/ Fprime x vals, F x vals.
3107         double diffMult = (i == 0) ? 1 : -1;
3108         foreach(FPos, Xi; F) {
3109             if(FPos < F.length - 1 && F[FPos + 1] == Xi)
3110                 continue;  //Handle ties.
3111             while(FprimePos < Fprime.length && Fprime[FprimePos] <= Xi) {
3112                 FprimePos++;
3113             }
3114             double diff = diffMult * (cast(double) (FPos + 1) / F.length -
3115                        cast(double) FprimePos / Fprime.length);
3116             if(abs(diff) > abs(D))
3117                 D = diff;
3118         }
3119         swap(F, Fprime);
3120         FprimePos = 0;
3121     }
3122     return D;
3123 }
3124 
3125 private double ksTestD(T, Func)(T Femp, Func F)
3126 if(doubleInput!(T) && is(ReturnType!Func : double)) {
3127     auto alloc = newRegionAllocator();
3128     return ksTestDDestructive(alloc.array(Femp), F);
3129 }
3130 
3131 private double ksTestDDestructive(T, Func)(T Femp, Func F)
3132 if(isArrayLike!(T) && is(ReturnType!Func : double)) {
3133     qsort(Femp);
3134     double D = 0;
3135 
3136     foreach(FPos, Xi; Femp) {
3137         double diff = cast(double) FPos / Femp.length - F(Xi);
3138         if(abs(diff) > abs(D))
3139             D = diff;
3140     }
3141 
3142     return D;
3143 }
3144 
3145 private double ksPval(ulong N, ulong Nprime, double D)
3146 in {
3147     assert(D >= -1);
3148     assert(D <= 1);
3149 } do {
3150     return 1 - kolmogorovDistrib(sqrt(cast(double) (N * Nprime) / (N + Nprime)) * abs(D));
3151 }
3152 
3153 private double ksPval(ulong N, double D)
3154 in {
3155     assert(D >= -1);
3156     assert(D <= 1);
3157 } do {
3158     return 1 - kolmogorovDistrib(abs(D) * sqrt(cast(double) N));
3159 }
3160 
3161 /**
3162 Wald-wolfowitz or runs test for randomness of the distribution of
3163 elements for which positive() evaluates to true.  For example, given
3164 a sequence of coin flips [H,H,H,H,H,T,T,T,T,T] and a positive() function of
3165 "a == 'H'", this test would determine that the heads are non-randomly
3166 distributed, since they are all at the beginning of obs.  This is done
3167 by counting the number of runs of consecutive elements for which
3168 positive() evaluates to true, and the number of consecutive runs for which
3169 it evaluates to false.  In the example above, we have 2 runs.  These are the
3170 block of 5 consecutive heads at the beginning and the 5 consecutive tails
3171 at the end.
3172 
3173 Alternatives are Alt.less, meaning that less runs than expected have been
3174 observed and data for which positive() is true tends to cluster,
3175 Alt.greater, which means that more runs than expected have been observed
3176 and data for which positive() is true tends to not cluster even moreso than
3177 expected by chance, and Alt.twoSided, meaning that elements for which
3178 positive() is true cluster as much as expected by chance.
3179 
3180 Bugs:  No exact calculation of the P-value.  Asymptotic approximation only.
3181 
3182 References:  http://en.wikipedia.org/wiki/Runs_test
3183 */
3184 double runsTest(alias positive = "a > 0", T)(T obs, Alt alt = Alt.twoSided)
3185 if(isIterable!(T)) {
3186     RunsTest!(positive, ForeachType!(T)) r;
3187     foreach(elem; obs) {
3188         r.put(elem);
3189     }
3190     return r.p(alt);
3191 }
3192 
3193 unittest {
3194     // Values from R lawstat package, for which "a < median(data)" is
3195     // hard-coded as the equivalent to positive().  The median of this data
3196     // is 0.5, so everything works.
3197     immutable int[] data = [1,0,0,0,1,1,0,0,1,0,1,0,1,0,1,1,1,0,0,1].idup;
3198     assert(approxEqual(runsTest(data), 0.3581));
3199     assert(approxEqual(runsTest(data, Alt.less), 0.821));
3200     assert(approxEqual(runsTest(data, Alt.greater), 0.1791));
3201 }
3202 
3203 /**
3204 Runs test as in runsTest(), except calculates online instead of from stored
3205 array elements.
3206 */
3207 struct RunsTest(alias positive = "a > 0", T) {
3208 private:
3209     uint nPos;
3210     uint nNeg;
3211     uint nRun;
3212     bool lastPos;
3213 
3214     alias unaryFun!(positive) pos;
3215 
3216 public:
3217 
3218     ///
3219     void put(T elem) {
3220         bool curPos = pos(elem);
3221         if(nRun == 0) {
3222             nRun = 1;
3223             if(curPos) {
3224                 nPos++;
3225             } else {
3226                 nNeg++;
3227             }
3228         } else if(pos(elem)) {
3229             nPos++;
3230             if(!lastPos) {
3231                 nRun++;
3232             }
3233         } else {
3234             nNeg++;
3235             if(lastPos) {
3236                 nRun++;
3237             }
3238         }
3239         lastPos = curPos;
3240     }
3241 
3242     ///
3243     uint nRuns() {
3244         return nRun;
3245     }
3246 
3247     ///
3248     double p(Alt alt = Alt.twoSided) {
3249         uint N = nPos + nNeg;
3250         double expected = 2.0 * nPos * nNeg / N + 1;
3251         double sd = sqrt((expected - 1) * (expected - 2) / (N - 1));
3252         if(alt == Alt.less) {
3253             return normalCDF(nRun, expected, sd);
3254         } else if(alt == Alt.greater) {
3255             return normalCDFR(nRun, expected, sd);
3256         } else {
3257             return 2 * ((nRun < expected) ?
3258                         normalCDF(nRun, expected, sd) :
3259                         normalCDFR(nRun, expected, sd));
3260         }
3261     }
3262 }
3263 
3264 deprecated {
3265     // Aliases for old names for correlation tests.
3266     alias pearsonCorTest pcorTest;
3267     alias spearmanCorTest scorTest;
3268     alias kendallCorTest kcorTest;
3269 }
3270 
3271 /**
3272 Tests the hypothesis that the Pearson correlation between two ranges is
3273 different from some 0.  Alternatives are Alt.less
3274 (pearsonCor(range1, range2) < 0), Alt.greater (pearsonCor(range1, range2)
3275  0) and Alt.twoSided (pearsonCor(range1, range2) != 0).
3276 
3277 Returns:  A ConfInt of the estimated Pearson correlation of the two ranges,
3278 the P-value against the given alternative, and the confidence interval of
3279 the correlation at the level specified by confLevel.
3280 
3281 References:  http://en.wikipedia.org/wiki/Pearson_correlation
3282 */
3283 ConfInt pearsonCorTest(T, U)(
3284     T range1,
3285     U range2,
3286     Alt alt = Alt.twoSided,
3287     double confLevel = 0.95
3288 ) if(doubleInput!(T) && doubleInput!(U)) {
3289     enforceConfidence(confLevel);
3290 
3291     auto pearsonRes = dstats.cor.pearsonCor(range1, range2);
3292     if(isNaN(pearsonRes.cor)) {
3293         return ConfInt.init;
3294     }
3295 
3296     return pearsonCorTest(pearsonRes.cor, pearsonRes.N, alt, confLevel);
3297 }
3298 
3299 /**
3300 Same as overload, but uses pre-computed correlation coefficient and sample
3301 size instead of computing them.
3302 
3303 Note:  This is a template only because of DMD Bug 2972.
3304  */
3305 ConfInt pearsonCorTest()(
3306     double cor,
3307     double N,
3308     Alt alt = Alt.twoSided,
3309     double confLevel = 0.95
3310 ) {
3311     dstatsEnforce(N >= 0, "N must be >= 0 for pearsonCorTest.");
3312     dstatsEnforce(cor > -1.0 || approxEqual(cor, -1.0),
3313         "Correlation must be between 0, 1.");
3314     dstatsEnforce(cor < 1.0 || approxEqual(cor, 1.0),
3315          "Correlation must be between 0, 1.");
3316     enforceConfidence(confLevel);
3317 
3318     immutable double denom = sqrt((1 - cor * cor) / (N - 2));
3319     immutable double t = cor / denom;
3320     ConfInt ret;
3321     ret.testStat = cor;
3322 
3323     double sqN, z;
3324     if(confLevel > 0) {
3325         sqN = sqrt(N - 3.0);
3326         z = sqN * atanh(cor);
3327     }
3328 
3329     final switch(alt) {
3330         case Alt.none :
3331             return ret;
3332         case Alt.twoSided:
3333             ret.p = (abs(cor) >= 1) ? 0 :
3334                 2 * ((t < 0) ? studentsTCDF(t, N - 2) : studentsTCDFR(t, N - 2));
3335 
3336             if(confLevel > 0) {
3337                 double deltaZ = invNormalCDF(0.5 * (1 - confLevel));
3338                 ret.lowerBound = tanh((z + deltaZ) / sqN);
3339                 ret.upperBound = tanh((z - deltaZ) / sqN);
3340             } else {
3341                 ret.lowerBound = cor;
3342                 ret.upperBound = cor;
3343             }
3344 
3345             break;
3346         case Alt.less:
3347             if(cor >= 1) {
3348                 ret.p = 1;
3349             } else if(cor <= -1) {
3350                 ret.p = 0;
3351             } else {
3352                 ret.p = studentsTCDF(t, N - 2);
3353             }
3354 
3355             if(confLevel > 0) {
3356                 double deltaZ = invNormalCDF(1 - confLevel);
3357                 ret.lowerBound = -1;
3358                 ret.upperBound = tanh((z - deltaZ) / sqN);
3359             } else {
3360                 ret.lowerBound = -1;
3361                 ret.upperBound = cor;
3362             }
3363 
3364             break;
3365         case Alt.greater:
3366             if(cor >= 1) {
3367                 ret.p = 0;
3368             } else if(cor <= -1) {
3369                 ret.p = 1;
3370             } else {
3371                 ret.p = studentsTCDFR(t, N - 2);
3372             }
3373 
3374             if(confLevel > 0) {
3375                 double deltaZ = invNormalCDF(1 - confLevel);
3376                 ret.lowerBound = tanh((z + deltaZ) / sqN);
3377                 ret.upperBound = 1;
3378             } else {
3379                 ret.lowerBound = cor;
3380                 ret.upperBound = 1;
3381             }
3382 
3383             break;
3384     }
3385     return ret;
3386 }
3387 
3388 unittest {
3389     // Values from R.
3390     auto t1 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.twoSided);
3391     auto t2 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.less);
3392     auto t3 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.greater);
3393 
3394     assert(approxEqual(t1.testStat, 0.8));
3395     assert(approxEqual(t2.testStat, 0.8));
3396     assert(approxEqual(t3.testStat, 0.8));
3397 
3398     assert(approxEqual(t1.p, 0.1041));
3399     assert(approxEqual(t2.p, 0.948));
3400     assert(approxEqual(t3.p, 0.05204));
3401 
3402     assert(approxEqual(t1.lowerBound, -0.2796400));
3403     assert(approxEqual(t3.lowerBound, -0.06438567));
3404     assert(approxEqual(t2.lowerBound, -1));
3405 
3406     assert(approxEqual(t1.upperBound, 0.9861962));
3407     assert(approxEqual(t2.upperBound, 0.9785289));
3408     assert(approxEqual(t3.upperBound, 1));
3409 
3410     // Test special case hack for cor = +- 1.
3411     uint[] myArr = [1,2,3,4,5];
3412     uint[] myArrReverse = myArr.dup;
3413     reverse(myArrReverse);
3414 
3415     auto t4 = pearsonCorTest(myArr, myArr, Alt.twoSided);
3416     auto t5 = pearsonCorTest(myArr, myArr, Alt.less);
3417     auto t6 = pearsonCorTest(myArr, myArr, Alt.greater);
3418     assert(approxEqual(t4.testStat, 1));
3419     assert(approxEqual(t4.p, 0));
3420     assert(approxEqual(t5.p, 1));
3421     assert(approxEqual(t6.p, 0));
3422 
3423     auto t7 = pearsonCorTest(myArr, myArrReverse, Alt.twoSided);
3424     auto t8 = pearsonCorTest(myArr, myArrReverse, Alt.less);
3425     auto t9 = pearsonCorTest(myArr, myArrReverse, Alt.greater);
3426     assert(approxEqual(t7.testStat, -1));
3427     assert(approxEqual(t7.p, 0));
3428     assert(approxEqual(t8.p, 0));
3429     assert(approxEqual(t9.p, 1));
3430 }
3431 
3432 /**
3433 Tests the hypothesis that the Spearman correlation between two ranges is
3434 different from some 0.  Alternatives are
3435 Alt.less (spearmanCor(range1, range2) < 0), Alt.greater (spearmanCor(range1, range2)
3436 > 0) and Alt.twoSided (spearmanCor(range1, range2) != 0).
3437 
3438 Returns:  A TestRes containing the Spearman correlation coefficient and
3439 the P-value for the given alternative.
3440 
3441 Bugs:  Exact P-value computation not yet implemented.  Uses asymptotic
3442 approximation only.  This is good enough for most practical purposes given
3443 reasonably large N, but is not perfectly accurate.  Not valid for data with
3444 very large amounts of ties.
3445 */
3446 TestRes spearmanCorTest(T, U)(T range1, U range2, Alt alt = Alt.twoSided)
3447 if(isInputRange!(T) && isInputRange!(U) &&
3448 is(typeof(range1.front < range1.front) == bool) &&
3449 is(typeof(range2.front < range2.front) == bool)) {
3450 
3451     static if(!hasLength!T) {
3452         auto alloc = newRegionAllocator();
3453         auto r1 = alloc.array(range1);
3454     } else {
3455         alias range1 r1;
3456     }
3457     immutable double N = r1.length;
3458 
3459     return pearsonCorTest(dstats.cor.spearmanCor(range1, range2), N, alt, 0);
3460 }
3461 
3462 unittest {
3463     // Values from R.
3464     int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
3465     int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];
3466     auto t1 = spearmanCorTest(arr1, arr2, Alt.twoSided);
3467     auto t2 = spearmanCorTest(arr1, arr2, Alt.less);
3468     auto t3 = spearmanCorTest(arr1, arr2, Alt.greater);
3469 
3470     assert(approxEqual(t1.testStat, -0.1769406));
3471     assert(approxEqual(t2.testStat, -0.1769406));
3472     assert(approxEqual(t3.testStat, -0.1769406));
3473 
3474     assert(approxEqual(t1.p, 0.4429));
3475     assert(approxEqual(t3.p, 0.7785));
3476     assert(approxEqual(t2.p, 0.2215));
3477 }
3478 
3479 /**
3480 Tests the hypothesis that the Kendall Tau-b between two ranges is
3481 different from 0.  Alternatives are Alt.less (kendallCor(range1, range2) < 0),
3482 Alt.greater (kendallCor(range1, range2) > 0) and Alt.twoSided
3483 (kendallCor(range1, range2) != 0).
3484 
3485 
3486 exactThresh controls the maximum length of the range for which exact P-value
3487 computation is used.  The default is 50.  Exact calculation is never used
3488 when ties are present because it is not computationally feasible.
3489 Do not set this higher than 100, as it will be very slow
3490 and the asymptotic approximation is pretty good at even a fraction of this
3491 size.
3492 
3493 Note:
3494 
3495 As an optimization, when a range is a SortedRange with predicate "a < b",
3496 it is assumed already sorted and not sorted a second time by this function.
3497 This is useful when applying this function multiple times with one of the
3498 arguments the same every time.
3499 
3500 Returns:  A TestRes containing the Kendall correlation coefficient and
3501 the P-value for the given alternative.
3502 
3503 References:  StackOverflow Question 948341 (http://stackoverflow.com/questions/948341)
3504 
3505 The Variance of Tau When Both Rankings Contain Ties.  M.G. Kendall.
3506 Biometrika, Vol 34, No. 3/4 (Dec., 1947), pp. 297-298
3507  */
3508 TestRes kendallCorTest(T, U)(
3509     T range1,
3510     U range2,
3511     Alt alt = Alt.twoSided,
3512     uint exactThresh = 50
3513 ) if(isInputRange!(T) && isInputRange!(U)) {
3514     auto alloc = newRegionAllocator();
3515 
3516     static if(dstats.cor.isDefaultSorted!T) {
3517         alias range1 i1d;
3518     } else {
3519         auto i1d = alloc.array(range1);
3520     }
3521 
3522     static if(dstats.cor.isDefaultSorted!U) {
3523         alias range2 i2d;
3524     } else {
3525         auto i2d = alloc.array(range2);
3526     }
3527 
3528     immutable res = dstats.cor.kendallCorDestructiveLowLevel(i1d, i2d, true);
3529     immutable double n = i1d.length;
3530 
3531     immutable double var =
3532           (2.0 / 9) * n * (n - 1) * (2 * n + 5)
3533         - (2.0 / 9) * res.tieCorrectT1
3534         - (2.0 / 9) * res.tieCorrectU1
3535         + (4 / (9 * n * (n - 1) * (n - 2))) * res.tieCorrectT2 * res.tieCorrectU2
3536         + 2 / (n * (n - 1)) * res.tieCorrectT3 * res.tieCorrectU3;
3537 
3538     // Need the / 2 to change C, as used in Kendall's paper to S, as used here.
3539     immutable double sd = sqrt(var) / 2;
3540 
3541     enum double cc = 1;
3542     auto tau = res.tau;
3543     auto s = res.s;
3544 
3545     immutable bool noTies = res.tieCorrectT1 == 0 && res.tieCorrectU1 == 0;
3546 
3547     if(noTies && n <= exactThresh) {
3548         // More than uint.max data points for exact calculation is
3549         // not plausible.
3550         assert(i1d.length < uint.max);
3551         immutable N = cast(uint) i1d.length;
3552         immutable nSwaps = (N * (N - 1) / 2 - cast(uint) s) / 2;
3553         return TestRes(tau, kendallCorExactP(N, nSwaps, alt));
3554     }
3555 
3556     final switch(alt) {
3557         case Alt.none :
3558             return TestRes(tau);
3559         case Alt.twoSided:
3560             if(abs(s) <= cc) {
3561                 return TestRes(tau, 1);
3562             } else if(s < 0) {
3563                 return TestRes(tau, 2 * normalCDF(s + cc, 0, sd));
3564             } else {
3565                 assert(s > 0);
3566                 return TestRes(tau, 2 * normalCDFR(s - cc, 0, sd));
3567             }
3568             assert(0);
3569 
3570         case Alt.less:
3571             return TestRes(tau, normalCDF(s + cc, 0, sd));
3572         case Alt.greater:
3573             return TestRes(tau, normalCDFR(s - cc, 0, sd));
3574     }
3575 }
3576 
3577 // Dynamic programming algorithm for computing exact Kendall tau P-values.
3578 // Thanks to ShreevatsaR from StackOverflow.
3579 private double kendallCorExactP(uint N, uint swaps, Alt alt) {
3580     uint maxSwaps = N * (N - 1) / 2;
3581     assert(swaps <= maxSwaps);
3582     immutable expectedSwaps = cast(ulong) N * (N - 1) * 0.25;
3583     if(alt == Alt.greater) {
3584         if(swaps > expectedSwaps) {
3585             if(swaps == maxSwaps) {
3586                 return 1;
3587             }
3588             return 1.0 - kendallCorExactP(N, maxSwaps - swaps - 1, Alt.greater);
3589         }
3590     } else if(alt == Alt.less) {
3591         if(swaps == 0) {
3592             return 1;
3593         }
3594         return kendallCorExactP(N, maxSwaps - swaps + 0, Alt.greater);
3595     } else if(alt == Alt.twoSided) {
3596         if(swaps < expectedSwaps) {
3597             return min(1, 2 * kendallCorExactP(N, swaps, Alt.greater));
3598         } else if(swaps > expectedSwaps) {
3599             return min(1, 2 * kendallCorExactP(N, swaps, Alt.less));
3600         } else {
3601             return 1;
3602         }
3603     } else {  // Alt.none
3604         return double.nan;
3605     }
3606 
3607     /* This algorithm was obtained from Question 948341 on stackoverflow.com
3608      * and works as follows:
3609      *
3610      * swaps is the number of swaps that would be necessary in a bubble sort
3611      * to sort one list in the same order as the other.  N is the sample size.
3612      * We want to find the number of ways that we could get a bubble sort
3613      * distance of at least swaps, and divide it by the total number of
3614      * permutations, pElem.
3615      *
3616      * The number of swaps necessary to sort a list is equivalent to the
3617      * number of inversions in the list, i.e. where i > j, but
3618      * list[i] < list[j].  This is a bottom-up dynamic programming algorithm
3619      * based on this principle.
3620      *
3621      * Assume c(N, k) is the number of permutations that require <= swaps
3622      * inversions.
3623      * We use the recurrence relation:
3624      * When k ≤ N - 1, c(N,k) = c(N,k-1) + c(N-1,k)
3625      * When k ≥ N,   c(N,k) = c(N,k-1) + c(N-1,k) - c(N-1,k-N)
3626      *
3627      * We also divide every value by the constant N! to turn this count into a
3628      * probability.
3629      */
3630 
3631     immutable double pElem = exp(-logFactorial(N));
3632     auto alloc = newRegionAllocator();
3633     double[] cur = alloc.uninitializedArray!(double[])(swaps + 1);
3634     double[] prev = alloc.uninitializedArray!(double[])(swaps + 1);
3635 
3636     prev[] = pElem;
3637     cur[0] = pElem;
3638     foreach(k; 1..N + 1) {
3639         immutable uint nSwapsPossible = k * (k - 1) / 2;
3640         immutable uint upTo = min(swaps, nSwapsPossible) + 1;
3641         foreach(j; 1..upTo) {
3642             if(j < k) {
3643                 cur[j] = prev[j] + cur[j - 1];
3644             } else {
3645                 cur[j] = prev[j] - prev[j - k] + cur[j - 1];
3646             }
3647         }
3648         cur[upTo..$] = cur[upTo - 1];
3649         swap(cur, prev);
3650     }
3651 
3652     return prev[$ - 1];
3653 }
3654 
3655 unittest {
3656     /* Values from R, with continuity correction enabled.  Note that large
3657      * one-sided inexact P-values are commented out because R seems to have a
3658      * slightly different interpretation of the proper continuity correction
3659      * than this library.  This library corrects the z-score in the direction
3660      * that would make the test more conservative.  R corrects towards zero.
3661      * I can't find a reference to support either one, but empirically it seems
3662      * like correcting towards more conservative results more accurately mirrors
3663      * the results of the exact test.  This isn't too big a deal anyhow since:
3664      *
3665      * 1.  The difference is small.
3666      * 2.  It only occurs on results that are very far from significance
3667      *     (P > 0.5).
3668      */
3669     int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
3670     int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];
3671     auto t1 = kendallCorTest(arr1, arr2, Alt.twoSided);
3672     auto t2 = kendallCorTest(arr1, arr2, Alt.less);
3673     auto t3 = kendallCorTest(arr1, arr2, Alt.greater);
3674 
3675     assert(approxEqual(t1.testStat, -.1448010));
3676     assert(approxEqual(t2.testStat, -.1448010));
3677     assert(approxEqual(t3.testStat, -.1448010));
3678 
3679     assert(approxEqual(t1.p, 0.3923745));
3680     //assert(approxEqual(t3.p, 0.8038127));
3681     assert(approxEqual(t2.p, 0.1961873));
3682 
3683     // Now, test the case of ties in both arrays.
3684     arr1 = [1,1,1,2,2,3,4,5,5,6];
3685     arr2 = [1,1,2,3,4,5,5,5,5,6];
3686     assert(approxEqual(kendallCorTest(arr1, arr2, Alt.twoSided).p, 0.001216776));
3687     //assert(approxEqual(kendallCorTest(arr1, arr2, Alt.less).p, 0.9993916));
3688     assert(approxEqual(kendallCorTest(arr1, arr2, Alt.greater).p, 0.0006083881));
3689 
3690     arr1 = [1,1,1,2,2,2,3,3,3,4,4,4,5,5,5];
3691     arr2 = [1,1,1,3,3,3,2,2,2,5,5,5,4,4,4];
3692     assert(approxEqual(kendallCorTest(arr1, arr2).p, 0.006123));
3693     assert(approxEqual(kendallCorTest(assumeSorted(arr1), arr2).p, 0.006123));
3694 
3695     // Test the exact stuff.  Still using values from R.
3696     uint[] foo = [1,2,3,4,5];
3697     uint[] bar = [1,2,3,5,4];
3698     uint[] baz = [5,3,1,2,4];
3699 
3700     assert(approxEqual(kendallCorTest(foo, foo).p, 0.01666666));
3701     assert(approxEqual(kendallCorTest(foo, foo, Alt.greater).p, 0.008333333));
3702     assert(approxEqual(kendallCorTest(foo, foo, Alt.less).p, 1));
3703 
3704     assert(approxEqual(kendallCorTest(foo, bar).p, 0.083333333));
3705     assert(approxEqual(kendallCorTest(foo, bar, Alt.greater).p, 0.041666667));
3706     assert(approxEqual(kendallCorTest(foo, bar, Alt.less).p, 0.9917));
3707 
3708     assert(approxEqual(kendallCorTest(foo, baz).p, 0.8167));
3709     assert(approxEqual(kendallCorTest(foo, baz, Alt.greater).p, 0.7583));
3710     assert(approxEqual(kendallCorTest(foo, baz, Alt.less).p, .4083));
3711 
3712     assert(approxEqual(kendallCorTest(bar, baz).p, 0.4833));
3713     assert(approxEqual(kendallCorTest(bar, baz, Alt.greater).p, 0.8833));
3714     assert(approxEqual(kendallCorTest(bar, baz, Alt.less).p, 0.2417));
3715 
3716     // A little monte carlo unittesting.  For large ranges, the deviation
3717     // between the exact and approximate version should be extremely small.
3718     foreach(i; 0..100) {
3719         uint nToTake = uniform(15, 65);
3720         auto lhs = array(take(randRange!rNormal(0, 1), nToTake));
3721         auto rhs = array(take(randRange!rNormal(0, 1), nToTake));
3722         if(i & 1) {
3723             lhs[] += rhs[] * 0.2;  // Make sure there's some correlation.
3724         } else {
3725             lhs[] -= rhs[] * 0.2;
3726         }
3727         double exact = kendallCorTest(lhs, rhs).p;
3728         double approx = kendallCorTest(lhs, rhs, Alt.twoSided, 0).p;
3729         assert(abs(exact - approx) < 0.01);
3730 
3731         exact = kendallCorTest(lhs, rhs, Alt.greater).p;
3732         approx = kendallCorTest(lhs, rhs, Alt.greater, 0).p;
3733         assert(abs(exact - approx) < 0.01);
3734 
3735         exact = kendallCorTest(lhs, rhs, Alt.less).p;
3736         approx = kendallCorTest(lhs, rhs, Alt.less, 0).p;
3737         assert(abs(exact - approx) < 0.01);
3738     }
3739 }
3740 
3741 /**A test for normality of the distribution of a range of values.  Based on
3742  * the assumption that normally distributed values will have a sample skewness
3743  * and sample kurtosis very close to zero.
3744  *
3745  * Returns:  A TestRes with the K statistic, which is Chi-Square distributed
3746  * with 2 degrees of freedom under the null, and the P-value for the alternative
3747  * that the data has skewness and kurtosis not equal to zero against the null
3748  * that skewness and kurtosis are near zero.  A normal distribution always has
3749  * skewness and kurtosis that converge to zero as sample size goes to infinity.
3750  *
3751  * Notes:  Contrary to popular belief, tests for normality should usually
3752  * not be used to deterimine whether T-tests are valid.  If the sample size is
3753  * large, T-tests are valid regardless of the distribution due to the central
3754  * limit theorem.  If the sample size is small, a test for normality will
3755  * likely not be very powerful, and a priori knowledge or simple inspection
3756  * of the data is often a better idea.
3757  *
3758  * References:
3759  * D'Agostino, Ralph B., Albert Belanger, and Ralph B. D'Agostino, Jr.
3760  * "A Suggestion for Using Powerful and Informative Tests of Normality",
3761  * The American Statistician, Vol. 44, No. 4. (Nov., 1990), pp. 316-321.
3762  */
3763 TestRes dAgostinoK(T)(T range)
3764 if(doubleIterable!(T)) {
3765     // You are not meant to understand this.  I sure don't.  I just implemented
3766     // these formulas off of Wikipedia, which got them from:
3767 
3768     // D'Agostino, Ralph B., Albert Belanger, and Ralph B. D'Agostino, Jr.
3769     // "A Suggestion for Using Powerful and Informative Tests of Normality",
3770     // The American Statistician, Vol. 44, No. 4. (Nov., 1990), pp. 316-321.
3771 
3772     // Amazing.  I didn't even realize things this complicated were possible
3773     // in 1990, before widespread computer algebra systems.
3774 
3775     // Notation from Wikipedia.  Keeping same notation for simplicity.
3776     double sqrtb1 = void, b2 = void, n = void;
3777     {
3778         auto summ = summary(range);
3779         sqrtb1 = summ.skewness;
3780         b2 = summ.kurtosis + 3;
3781         n = summ.N;
3782     }
3783 
3784     // Calculate transformed skewness.
3785     double Y = sqrtb1 * sqrt((n + 1) * (n + 3) / (6 * (n - 2)));
3786     double beta2b1Numer = 3 * (n * n + 27 * n - 70) * (n + 1) * (n + 3);
3787     double beta2b1Denom = (n - 2) * (n + 5) * (n + 7) * (n + 9);
3788     double beta2b1 = beta2b1Numer / beta2b1Denom;
3789     double Wsq = -1 + sqrt(2 * (beta2b1 - 1));
3790     double delta = 1.0 / sqrt(log(sqrt(Wsq)));
3791     double alpha = sqrt( 2.0 / (Wsq - 1));
3792     double Zb1 = delta * log(Y / alpha + sqrt(pow(Y / alpha, 2) + 1));
3793 
3794     // Calculate transformed kurtosis.
3795     double Eb2 = 3 * (n - 1) / (n + 1);
3796     double sigma2b2 = (24 * n * (n - 2) * (n - 3)) / (
3797         (n + 1) * (n + 1) * (n + 3) * (n + 5));
3798     double x = (b2 - Eb2) / sqrt(sigma2b2);
3799 
3800     double sqBeta1b2 = 6 * (n * n - 5 * n + 2) / ((n + 7) * (n + 9)) *
3801          sqrt((6 * (n + 3) * (n + 5)) / (n * (n - 2) * (n - 3)));
3802     double A = 6 + 8 / sqBeta1b2 * (2 / sqBeta1b2 + sqrt(1 + 4 / (sqBeta1b2 * sqBeta1b2)));
3803     double Zb2 = ((1 - 2 / (9 * A)) -
3804         cbrt((1 - 2 / A) / (1 + x * sqrt(2 / (A - 4)))) ) *
3805         sqrt(9 * A / 2);
3806 
3807     double K2 = Zb1 * Zb1 + Zb2 * Zb2;
3808 
3809     if(isNaN(K2)) {
3810         return TestRes(double.nan, double.nan);
3811     }
3812 
3813     return TestRes(K2, chiSquareCDFR(K2, 2));
3814 }
3815 
3816 unittest {
3817     // Values from R's fBasics package.
3818     int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
3819     int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];
3820 
3821     auto r1 = dAgostinoK(arr1);
3822     auto r2 = dAgostinoK(arr2);
3823 
3824     assert(approxEqual(r1.testStat, 3.1368));
3825     assert(approxEqual(r1.p, 0.2084));
3826 
3827     assert(approxEqual(r2.testStat, 1.1816));
3828     assert(approxEqual(r2.p, 0.5539));
3829 }
3830 
3831 /**Fisher's method of meta-analyzing a set of P-values to determine whether
3832  * there are more significant results than would be expected by chance.
3833  * Based on a chi-square statistic for the sum of the logs of the P-values.
3834  *
3835  * Returns:  A TestRes containing the chi-square statistic and a P-value for
3836  * the alternative hypothesis that more small P-values than would be expected
3837  * by chance are present against the alternative that the distribution of
3838  * P-values is uniform or enriched for large P-values.
3839  *
3840  * References:  Fisher, R. A. (1948) "Combining independent tests of
3841  * significance", American Statistician, vol. 2, issue 5, page 30.
3842  * (In response to Question 14)
3843  */
3844 TestRes fishersMethod(R)(R pVals)
3845 if(doubleInput!R) {
3846     double chiSq = 0;
3847     uint df = 0;
3848     foreach(pVal; pVals) {
3849         dstatsEnforce(pVal >= 0 && pVal <= 1,
3850             "P-values must be between 0, 1 for Fisher's Method.");
3851         chiSq += log( cast(double) pVal);
3852         df += 2;
3853     }
3854     chiSq *= -2;
3855     return TestRes(chiSq, chiSquareCDFR(chiSq, df));
3856 }
3857 
3858 unittest {
3859     // First, basic sanity check.  Make sure w/ one P-value, we get back that
3860     // P-value.
3861     for(double p = 0.01; p < 1; p += 0.01) {
3862         assert(approxEqual(fishersMethod([p].dup).p, p));
3863     }
3864     float[] ps = [0.739, 0.0717, 0.01932, 0.03809];
3865     auto res = fishersMethod(ps);
3866     assert(approxEqual(res.testStat, 20.31));
3867     assert(res.p < 0.01);
3868 }
3869 
3870 /// For falseDiscoveryRate.
3871 enum Dependency : bool {
3872     /**
3873     Assume that dependency among hypotheses may exist.  (More conservative,
3874     Benjamini-Yekutieli procedure.)
3875     */
3876     yes = true,
3877 
3878     /**
3879     Assume hypotheses are independent.  (Less conservative, Benjamine-
3880     Hochberg procedure.)
3881     */
3882     no = false
3883 }
3884 
3885 /**
3886 Computes the false discovery rate statistic given a list of
3887 p-values, according to Benjamini and Hochberg (1995) (independent) or
3888 Benjamini and Yekutieli (2001) (dependent).  The Dependency parameter
3889 controls whether hypotheses are assumed to be independent, or whether
3890 the more conservative assumption that they are correlated must be made.
3891 
3892 Returns:
3893 An array of adjusted P-values with indices corresponding to the order of
3894 the P-values in the input data.
3895 
3896 References:
3897 Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate:
3898 a practical and powerful approach to multiple testing. Journal of the Royal
3899 Statistical Society Series B, 57, 289-200
3900 
3901 Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery
3902 rate in multiple testing under dependency. Annals of Statistics 29, 1165-1188.
3903  */
3904 float[] falseDiscoveryRate(T)(T pVals, Dependency dep = Dependency.no)
3905 if(doubleInput!(T)) {
3906     auto qVals = array(map!(to!float)(pVals));
3907 
3908     double C = 1;
3909     if(dep == Dependency.yes) {
3910         foreach(i; 2..qVals.length + 1) {
3911             C += 1.0 / i;
3912         }
3913     }
3914 
3915     auto alloc = newRegionAllocator();
3916     auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);
3917     foreach(i, ref elem; perm) {
3918         elem = i;
3919     }
3920 
3921     qsort(qVals, perm);
3922 
3923     foreach(i, ref q; qVals) {
3924         q = min(1.0f, q * C * cast(double) qVals.length / (cast(double) i + 1));
3925     }
3926 
3927     float smallestSeen = float.max;
3928     foreach(ref q; retro(qVals)) {
3929         if(q < smallestSeen) {
3930             smallestSeen = q;
3931         } else {
3932             q = smallestSeen;
3933         }
3934     }
3935 
3936     qsort(perm, qVals);  //Makes order of qVals correspond to input.
3937     return qVals;
3938 }
3939 
3940 unittest {
3941     // Comparing results to R.
3942     auto pVals = [.90, .01, .03, .03, .70, .60, .01].dup;
3943     auto qVals = falseDiscoveryRate(pVals);
3944     alias approxEqual ae;
3945     assert(ae(qVals[0], .9));
3946     assert(ae(qVals[1], .035));
3947     assert(ae(qVals[2], .052));
3948     assert(ae(qVals[3], .052));
3949     assert(ae(qVals[4], .816666666667));
3950     assert(ae(qVals[5], .816666666667));
3951     assert(ae(qVals[6], .035));
3952 
3953     auto p2 = [.1, .02, .6, .43, .001].dup;
3954     auto q2 = falseDiscoveryRate(p2);
3955     assert(ae(q2[0], .16666666));
3956     assert(ae(q2[1], .05));
3957     assert(ae(q2[2], .6));
3958     assert(ae(q2[3], .5375));
3959     assert(ae(q2[4], .005));
3960 
3961     // Dependent case.
3962     qVals = falseDiscoveryRate(pVals, Dependency.yes);
3963     assert(ae(qVals[0], 1));
3964     assert(ae(qVals[1], .09075));
3965     assert(ae(qVals[2], .136125));
3966     assert(ae(qVals[3], .136125));
3967     assert(ae(qVals[4], 1));
3968     assert(ae(qVals[5], 1));
3969     assert(ae(qVals[6], .09075));
3970 
3971     q2 = falseDiscoveryRate(p2, Dependency.yes);
3972     assert(ae(q2[0], .38055555));
3973     assert(ae(q2[1], .1141667));
3974     assert(ae(q2[2], 1));
3975     assert(ae(q2[3], 1));
3976     assert(ae(q2[4], .01141667));
3977 }
3978 
3979 /**Uses the Hochberg procedure to control the familywise error rate assuming
3980  * that hypothesis tests are independent.  This is more powerful than
3981  * Holm-Bonferroni correction, but requires the independence assumption.
3982  *
3983  * Returns:
3984  * An array of adjusted P-values with indices corresponding to the order of
3985  * the P-values in the input data.
3986  *
3987  * References:
3988  * Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of
3989  * significance. Biometrika, 75, 800-803.
3990  */
3991 float[] hochberg(T)(T pVals)
3992 if(doubleInput!(T)) {
3993     auto qVals = array(map!(to!float)(pVals));
3994 
3995     auto alloc = newRegionAllocator();
3996     auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);
3997     foreach(i, ref elem; perm)
3998         elem = i;
3999 
4000     qsort(qVals, perm);
4001 
4002     foreach(i, ref q; qVals) {
4003         dstatsEnforce(q >= 0 && q <= 1,
4004             "P-values must be between 0, 1 for hochberg.");
4005         q = min(1.0f, q * (cast(double) qVals.length - i));
4006     }
4007 
4008     float smallestSeen = float.max;
4009     foreach(ref q; retro(qVals)) {
4010         if(q < smallestSeen) {
4011             smallestSeen = q;
4012         } else {
4013             q = smallestSeen;
4014         }
4015     }
4016 
4017     qsort(perm, qVals);  //Makes order of qVals correspond to input.
4018     return qVals;
4019 }
4020 
4021 unittest {
4022     alias approxEqual ae;
4023     auto q = hochberg([0.01, 0.02, 0.025, 0.9].dup);
4024     assert(ae(q[0], 0.04));
4025     assert(ae(q[1], 0.05));
4026     assert(ae(q[2], 0.05));
4027     assert(ae(q[3], 0.9));
4028 
4029     auto p2 = [.1, .02, .6, .43, .001].dup;
4030     auto q2 = hochberg(p2);
4031     assert(ae(q2[0], .3));
4032     assert(ae(q2[1], .08));
4033     assert(ae(q2[2], .6));
4034     assert(ae(q2[3], .6));
4035     assert(ae(q2[4], .005));
4036 }
4037 
4038 /**Uses the Holm-Bonferroni method to adjust a set of P-values in a way that
4039  * controls the familywise error rate (The probability of making at least one
4040  * Type I error).  This is basically a less conservative version of
4041  * Bonferroni correction that is still valid for arbitrary assumptions and
4042  * controls the familywise error rate.  Therefore, there aren't too many good
4043  * reasons to use regular Bonferroni correction instead.
4044  *
4045  * Returns:
4046  * An array of adjusted P-values with indices corresponding to the order of
4047  * the P-values in the input data.
4048  *
4049  * References:
4050  * Holm, S. (1979). A simple sequentially rejective multiple test procedure.
4051  * Scandinavian Journal of Statistics, 6, 65-70.
4052  */
4053 float[] holmBonferroni(T)(T pVals)
4054 if(doubleInput!(T)) {
4055     auto alloc = newRegionAllocator();
4056 
4057     auto qVals = array(map!(to!float)(pVals));
4058     auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);
4059 
4060     foreach(i, ref elem; perm) {
4061         elem = i;
4062     }
4063 
4064     qsort(qVals, perm);
4065 
4066     foreach(i, ref q; qVals) {
4067         dstatsEnforce(q >= 0 && q <= 1,
4068             "P-values must be between 0, 1 for holmBonferroni.");
4069         q = min(1.0, q * (cast(double) qVals.length - i));
4070     }
4071 
4072     foreach(i; 1..qVals.length) {
4073         if(qVals[i] < qVals[i - 1]) {
4074             qVals[i] = qVals[i - 1];
4075         }
4076     }
4077 
4078     qsort(perm, qVals);  //Makes order of qVals correspond to input.
4079     return qVals;
4080 }
4081 
4082 unittest {
4083     // Values from R.
4084     auto ps = holmBonferroni([0.001, 0.2, 0.3, 0.4, 0.7].dup);
4085     alias approxEqual ae;
4086     assert(ae(ps[0], 0.005));
4087     assert(ae(ps[1], 0.8));
4088     assert(ae(ps[2], 0.9));
4089     assert(ae(ps[3], 0.9));
4090     assert(ae(ps[4], 0.9));
4091 
4092     ps = holmBonferroni([0.3, 0.1, 0.4, 0.1, 0.5, 0.9].dup);
4093     assert(ps == [1f, 0.6f, 1f, 0.6f, 1f, 1f]);
4094 }
4095 
4096 // old unconstrained approxEqual to work around https://issues.dlang.org/show_bug.cgi?id=18287
4097 static if (__VERSION__ == 2078) private
4098 {
4099     /**
4100        Computes whether two values are approximately equal, admitting a maximum
4101        relative difference, and a maximum absolute difference.
4102        Params:
4103             lhs = First item to compare.
4104             rhs = Second item to compare.
4105             maxRelDiff = Maximum allowable difference relative to `rhs`.
4106             maxAbsDiff = Maximum absolute difference.
4107        Returns:
4108            `true` if the two items are approximately equal under either criterium.
4109            If one item is a range, and the other is a single value, then the result
4110            is the logical and-ing of calling `approxEqual` on each element of the
4111            ranged item against the single item. If both items are ranges, then
4112            `approxEqual` returns `true` if and only if the ranges have the same
4113            number of elements and if `approxEqual` evaluates to `true` for each
4114            pair of elements.
4115      */
4116     bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff, V maxAbsDiff = 1e-5)
4117     {
4118         import std.range.primitives : empty, front, isInputRange, popFront;
4119         static if (isInputRange!T)
4120         {
4121             static if (isInputRange!U)
4122             {
4123                 // Two ranges
4124                 for (;; lhs.popFront(), rhs.popFront())
4125                 {
4126                     if (lhs.empty) return rhs.empty;
4127                     if (rhs.empty) return lhs.empty;
4128                     if (!approxEqual(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
4129                         return false;
4130                 }
4131             }
4132             else static if (isIntegral!U)
4133             {
4134                 // convert rhs to real
4135                 return approxEqual(lhs, real(rhs), maxRelDiff, maxAbsDiff);
4136             }
4137             else
4138             {
4139                 // lhs is range, rhs is number
4140                 for (; !lhs.empty; lhs.popFront())
4141                 {
4142                     if (!approxEqual(lhs.front, rhs, maxRelDiff, maxAbsDiff))
4143                         return false;
4144                 }
4145                 return true;
4146             }
4147         }
4148         else
4149         {
4150             static if (isInputRange!U)
4151             {
4152                 // lhs is number, rhs is range
4153                 for (; !rhs.empty; rhs.popFront())
4154                 {
4155                     if (!approxEqual(lhs, rhs.front, maxRelDiff, maxAbsDiff))
4156                         return false;
4157                 }
4158                 return true;
4159             }
4160             else static if (isIntegral!T || isIntegral!U)
4161             {
4162                 // convert both lhs and rhs to real
4163                 return approxEqual(real(lhs), real(rhs), maxRelDiff, maxAbsDiff);
4164             }
4165             else
4166             {
4167                 // two numbers
4168                 //static assert(is(T : real) && is(U : real));
4169                 if (rhs == 0)
4170                 {
4171                     return fabs(lhs) <= maxAbsDiff;
4172                 }
4173                 static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
4174                 {
4175                     if (lhs == lhs.infinity && rhs == rhs.infinity ||
4176                         lhs == -lhs.infinity && rhs == -rhs.infinity) return true;
4177                 }
4178                 return fabs((lhs - rhs) / rhs) <= maxRelDiff
4179                     || maxAbsDiff != 0 && fabs(lhs - rhs) <= maxAbsDiff;
4180             }
4181         }
4182     }
4183 
4184     /**
4185        Returns $(D approxEqual(lhs, rhs, 1e-2, 1e-5)).
4186      */
4187     bool approxEqual(T, U)(T lhs, U rhs)
4188     {
4189         return approxEqual(lhs, rhs, 1e-2, 1e-5);
4190     }
4191 }