/**Hypothesis testing beyond simple CDFs.  All functions work with input
 * ranges with elements implicitly convertible to double unless otherwise noted.
 *
 * Author:  David Simcha*/
 /*
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */
module dstats.tests;

import std.functional, std.range, std.conv, std.math, std.traits,
       std.exception, std.typetuple;

import std.algorithm : reverse, copy, max, min, swap, map, reduce;

import dstats.base, dstats.distrib, dstats.alloc, dstats.summary, dstats.sort;

static import dstats.cor;

private static import dstats.infotheory;

version(unittest) {
    import std.stdio, dstats.random;
}

/**Alternative hypotheses.  Exact meaning varies with test used.*/
enum Alt {
    /// f(input1) != X
    twoSided,

    /// f(input1) < X
    less,

    /// f(input1) > X
    greater,

    /**
    Skip P-value computation (and confidence intervals if applicable)
    and just return the test statistic.
    */
    none
}

/**
A plain old data struct for returning the results of hypothesis tests.
*/
struct TestRes {
    /// The test statistic.  What exactly this is is specific to the test.
    double testStat;

    /**The P-value against the provided alternative.  This struct can
     * be implicitly converted to just the P-value via alias this.*/
    double p;

    /// Allow implicit conversion to the P-value.
    alias p this;

    ///
    string toString() {
        return text("Test Statistic = ", testStat, "\nP = ", p);
    }
}

/**
A plain old data struct for returning the results of hypothesis tests
that also produce confidence intervals.  Contains, can implicitly convert
to, a TestRes.
*/
struct ConfInt {
    ///  This is alias this'd.
    TestRes testRes;

    ///  Lower bound of the confidence interval at the level specified.
    double lowerBound;

    ///  Upper bound of the confidence interval at the level specified.
    double upperBound;

    alias testRes this;

    ///
    string toString() {
        return text("Test Statistic = ", testRes.testStat, "\nP = ", testRes.p,
                "\nLower Confidence Bound = ", lowerBound,
                "\nUpper Confidence Bound = ", upperBound);
    }
}

/**
Tests whether a struct/class has the necessary information for calculating
a T-test.  It must have a property .mean (mean), .stdev (stdandard deviation),
.var (variance), and .N (sample size).
*/
template isSummary(T) {
    enum bool isSummary = is(typeof(T.init.mean)) && is(typeof(T.init.stdev)) &&
        is(typeof(T.init.var)) && is(typeof(T.init.N));
}

/**
One-sample Student's T-test for difference between mean of data and
 a fixed value.  Alternatives are Alt.less, meaning mean(data) < testMean,
 Alt.greater, meaning mean(data) > testMean, and Alt.twoSided, meaning
 mean(data)!= testMean.

 data may be either an iterable with elements implicitly convertible to
 double or a summary struct (see isSummary).
 *
 Examples:
 ---
 uint[] data = [1,2,3,4,5];

 // Test the null hypothesis that the mean of data is >= 1 against the
 // alternative that the mean of data is < 1.  Calculate confidence
 // intervals at 90%.
 auto result1 = studentsTTest(data, 1, Alt.less, 0.9);

 // Do the same thing, only this time we've already calculated the summary
 // statistics explicitly before passing them to studensTTest.
 auto summary = meanStdev(data);
 writeln(summary.stdev);
 result2 = studentsTTest(summary, 1, Alt.less, 0.9);  // Same as result1.
 assert(result1 == result2);
 ---

 Returns:  A ConfInt containing T, the P-value and the boundaries of
 the confidence interval for mean(data) at the level specified.

 References:  http://en.wikipedia.org/wiki/Student%27s_t-test
 */
ConfInt studentsTTest(T)(
    T data,
    double testMean = 0,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) if( (isSummary!T || doubleIterable!T)) {
    enforceConfidence(confLevel);
    dstatsEnforce(isFinite(testMean), "testMean must not be infinite or NaN.");

    static if(isSummary!T) {
        return pairedTTest(data, testMean, alt, confLevel);
    } else static if(doubleIterable!T) {
        return pairedTTest(meanStdev(data), testMean, alt, confLevel);
    }
}

unittest {
    auto t1 = studentsTTest([1, 2, 3, 4, 5].dup, 2);
    assert(approxEqual(t1.testStat, 1.4142));
    assert(approxEqual(t1.p, 0.2302));
    assert(approxEqual(t1.lowerBound, 1.036757));
    assert(approxEqual(t1.upperBound, 4.963243));
    assert(t1 == studentsTTest( meanStdev([1,2,3,4,5].dup), 2));

    auto t2 = studentsTTest([1, 2, 3, 4, 5].dup, 2, Alt.less);
    assert(approxEqual(t2.p, .8849));
    assert(approxEqual(t2.testStat, 1.4142));
    assert(t2.lowerBound == -double.infinity);
    assert(approxEqual(t2.upperBound, 4.507443));

    auto t3 = studentsTTest( summary([1, 2, 3, 4, 5].dup), 2, Alt.greater);
    assert(approxEqual(t3.p, .1151));
    assert(approxEqual(t3.testStat, 1.4142));
    assert(approxEqual(t3.lowerBound, 1.492557));
    assert(t3.upperBound == double.infinity);
}

/**
Two-sample T test for a difference in means,
assumes variances of samples are equal.  Alteratives are Alt.less, meaning
mean(sample1) - mean(sample2) < testMean, Alt.greater, meaning
mean(sample1) - mean(sample2) > testMean, and Alt.twoSided, meaning
mean(sample1) - mean(sample2) != testMean.

sample1 and sample2 may be either iterables with elements implicitly
convertible to double or summary structs (see isSummary).

Returns:  A ConfInt containing the T statistic, the P-value, and the
boundaries of the confidence interval for the difference between means
of sample1 and sample2 at the specified level.

References:  http://en.wikipedia.org/wiki/Student%27s_t-test
 */
ConfInt studentsTTest(T, U)(
    T sample1,
    U sample2,
    double testMean = 0,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) if( (doubleIterable!T || isSummary!T) && (doubleIterable!U || isSummary!U)) {
    enforceConfidence(confLevel);
    dstatsEnforce(isFinite(testMean), "testMean must not be infinite or nan.");

    static if(isSummary!T) {
        alias sample1 s1summ;
    } else {
        immutable s1summ = meanStdev(sample1);
    }

    static if(isSummary!U) {
        alias sample2 s2summ;
    } else {
        immutable s2summ = meanStdev(sample2);
    }

    immutable n1 = s1summ.N, n2 = s2summ.N;

    immutable sx1x2 = sqrt((n1 * s1summ.mse + n2 * s2summ.mse) /
                 (n1 + n2 - 2));
    immutable normSd = (sx1x2 * sqrt((1.0 / n1) + (1.0 / n2)));
    immutable meanDiff = s1summ.mean - s2summ.mean;
    ConfInt ret;
    ret.testStat = (meanDiff - testMean) / normSd;
    if(alt == Alt.none) {
        return ret;
    } else if(alt == Alt.less) {
        ret.p = studentsTCDF(ret.testStat, n1 + n2 - 2);

        ret.lowerBound = -double.infinity;
        if(confLevel > 0) {
            immutable delta = invStudentsTCDF(1 - confLevel, n1 + n2 - 2)
                * normSd;
            ret.upperBound = meanDiff - delta;
        } else {
            ret.upperBound = meanDiff;
        }

    } else if(alt == Alt.greater) {
        ret.p = studentsTCDF(-ret.testStat, n1 + n2 - 2);

        ret.upperBound = double.infinity;
        if(confLevel > 0) {
            immutable delta = invStudentsTCDF(1 - confLevel, n1 + n2 - 2)
                * normSd;
            ret.lowerBound = meanDiff + delta;
        } else {
            ret.lowerBound = meanDiff;
        }

    } else {
        immutable t = ret.testStat;
        ret.p = 2 * ((t < 0) ?
                    studentsTCDF(t, n1 + n2 - 2) :
                    studentsTCDFR(t, n1 + n2 - 2));

        if(confLevel > 0) {
            immutable delta = invStudentsTCDF(
                0.5 * (1 - confLevel), n1 + n2 - 2) * normSd;
            ret.lowerBound = meanDiff + delta;
            ret.upperBound = meanDiff - delta;
        } else {
            ret.lowerBound = meanDiff;
            ret.upperBound = meanDiff;
        }
    }
    return ret;
}

unittest {
    // Values from R.
    auto t1 = studentsTTest([1,2,3,4,5], [1,3,4,5,7,9]);
    assert(approxEqual(t1.p, 0.2346));
    assert(approxEqual(t1.testStat, -1.274));
    assert(approxEqual(t1.lowerBound, -5.088787));
    assert(approxEqual(t1.upperBound, 1.422120));


    assert(approxEqual(studentsTTest([1,2,3,4,5], [1,3,4,5,7,9], 0, Alt.less),
           0.1173));
    assert(approxEqual(studentsTTest([1,2,3,4,5], [1,3,4,5,7,9], 0, Alt.greater),
           0.8827));
    auto t2 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 5);
    assert(approxEqual(t2.p, 0.44444));
    assert(approxEqual(t2.testStat, -0.7998));
    assert(approxEqual(t2.lowerBound, -0.3595529));
    assert(approxEqual(t2.upperBound, 7.5595529));


    auto t5 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 0, Alt.less);
    assert(approxEqual(t5.p, 0.965));
    assert(approxEqual(t5.testStat, 2.0567));
    assert(approxEqual(t5.upperBound, 6.80857));
    assert(t5.lowerBound == -double.infinity);

    auto t6 = studentsTTest([1,3,5,7,9,11], [2,2,1,3,4], 0, Alt.greater);
    assert(approxEqual(t6.p, 0.03492));
    assert(approxEqual(t6.testStat, 2.0567));
    assert(approxEqual(t6.lowerBound, 0.391422));
    assert(t6.upperBound == double.infinity);

    auto t7 = studentsTTest([1, 2, 4], [3]);
    assert(approxEqual(t7.p, 0.7418));
    assert(approxEqual(t7.testStat, 0.-.378));
    assert(approxEqual(t7.lowerBound, -8.255833));
    assert(approxEqual(t7.upperBound, 6.922499));

}

/**
Two-sample T-test for difference in means.  Does not assume variances are equal.
Alteratives are Alt.less, meaning mean(sample1) - mean(sample2) < testMean,
Alt.greater, meaning mean(sample1) - mean(sample2) > testMean, and
Alt.twoSided, meaning mean(sample1) - mean(sample2) != testMean.

sample1 and sample2 may be either iterables with elements implicitly
convertible to double or summary structs (see isSummary).

Returns:  A ConfInt containing the T statistic, the P-value, and the
boundaries of the confidence interval for the difference between means
of sample1 and sample2 at the specified level.

References:  http://en.wikipedia.org/wiki/Student%27s_t-test
 */
ConfInt welchTTest(T, U)(T sample1, U sample2, double testMean = 0,
    Alt alt = Alt.twoSided, double confLevel = 0.95)
if( (isSummary!T || doubleIterable!T) && (isSummary!U || doubleIterable!U)) {
    enforceConfidence(confLevel);
    dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or NaN.");

    static if(isSummary!T) {
        alias sample1 s1summ;
    } else {
        auto s1summ = meanStdev(sample1);
    }

    static if(isSummary!U) {
        alias sample2 s2summ;
    } else {
        auto s2summ = meanStdev(sample2);
    }

    immutable double n1 = s1summ.N,
                   n2 = s2summ.N;

    immutable v1 = s1summ.var, v2 = s2summ.var;
    immutable double sx1x2 = sqrt(v1 / n1 + v2 / n2);
    immutable double meanDiff = s1summ.mean - s2summ.mean - testMean;
    immutable double t = meanDiff / sx1x2;
    double numerator = v1 / n1 + v2 / n2;
    numerator *= numerator;
    double denom1 = v1 / n1;
    denom1 = denom1 * denom1 / (n1 - 1);
    double denom2 = v2 / n2;
    denom2 = denom2 * denom2 / (n2 - 1);
    immutable double df = numerator / (denom1 + denom2);

    ConfInt ret;
    ret.testStat = t;
    if(alt == Alt.none) {
        return ret;
    } else if(alt == Alt.less) {
        ret.p = studentsTCDF(t, df);
        ret.lowerBound = -double.infinity;

        if(confLevel > 0) {
            ret.upperBound = meanDiff +
                testMean - invStudentsTCDF(1 - confLevel, df) * sx1x2;
        } else {
            ret.upperBound = meanDiff + testMean;
        }

    } else if(alt == Alt.greater) {
        ret.p = studentsTCDF(-t, df);
        ret.upperBound = double.infinity;

        if(confLevel > 0) {
            ret.lowerBound = meanDiff +
                testMean + invStudentsTCDF(1 - confLevel, df) * sx1x2;
        } else {
            ret.lowerBound = meanDiff + testMean;
        }

    } else {
        ret.p = 2 * ((t < 0) ?
                     studentsTCDF(t, df) :
                     studentsTCDF(-t, df));

        if(confLevel > 0) {
            double delta = invStudentsTCDF(0.5 * (1 - confLevel), df) * sx1x2;
            ret.upperBound = meanDiff + testMean - delta;
            ret.lowerBound = meanDiff + testMean + delta;
        } else {
            ret.upperBound = meanDiff + testMean;
            ret.lowerBound = meanDiff + testMean;
        }
    }
    return ret;
}

unittest {
    // Values from R.
    auto t1 = welchTTest( meanStdev([1,2,3,4,5]), [1,3,4,5,7,9], 2);
    assert(approxEqual(t1.p, 0.02285));
    assert(approxEqual(t1.testStat, -2.8099));
    assert(approxEqual(t1.lowerBound, -4.979316));
    assert(approxEqual(t1.upperBound, 1.312649));

    auto t2 = welchTTest([1,2,3,4,5], summary([1,3,4,5,7,9]), -1, Alt.less);
    assert(approxEqual(t2.p, 0.2791));
    assert(approxEqual(t2.testStat, -0.6108));
    assert(t2.lowerBound == -double.infinity);
    assert(approxEqual(t2.upperBound, 0.7035534));

    auto t3 = welchTTest([1,2,3,4,5], [1,3,4,5,7,9], 0.5, Alt.greater);
    assert(approxEqual(t3.p, 0.9372));
    assert(approxEqual(t3.testStat, -1.7104));
    assert(approxEqual(t3.lowerBound, -4.37022));
    assert(t3.upperBound == double.infinity);

    assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4]).p, 0.06616));
    assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4], 0,
        Alt.less).p, 0.967));
    assert(approxEqual(welchTTest([1,3,5,7,9,11], [2,2,1,3,4], 0,
        Alt.greater).p, 0.03308));
}

/**
Paired T test.  Tests the hypothesis that the mean difference between
corresponding elements of before and after is testMean.  Alternatives are
Alt.less, meaning the that the true mean difference (before[i] - after[i])
is less than testMean, Alt.greater, meaning the true mean difference is
greater than testMean, and Alt.twoSided, meaning the true mean difference is not
equal to testMean.

before and after must be input ranges with elements implicitly convertible
to double and must have the same length.

Returns:  A ConfInt containing the T statistic, the P-value, and the
boundaries of the confidence interval for the mean difference between
corresponding elements of sample1 and sample2 at the specified level.

References:  http://en.wikipedia.org/wiki/Student%27s_t-test
*/
ConfInt pairedTTest(T, U)(
    T before,
    U after,
    double testMean = 0,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) if(doubleInput!(T) && doubleInput!(U) && isInputRange!T && isInputRange!U) {
    enforceConfidence(confLevel);
    dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");

    MeanSD msd;
    while(!before.empty && !after.empty) {
        immutable diff = cast(double) before.front - cast(double) after.front;
        before.popFront();
        after.popFront();
        msd.put(diff);
    }

    dstatsEnforce(before.empty && after.empty,
        "before and after have different lengths in pairedTTest.");

    return pairedTTest(msd, testMean, alt, confLevel);
}

/**
Compute a paired T test directly from summary statistics of the differences
between corresponding samples.

Examples:
---
float[] data1 = [8, 6, 7, 5, 3, 0, 9];
float[] data2 = [3, 6, 2, 4, 3, 6, 8];

// Calculate summary statistics on difference explicitly.
MeanSD summary;
foreach(i; 0..data1.length) {
    summary.put(data1[i] - data2[i]);
}

// Test the null hypothesis that the mean difference between corresponding
// elements (data1[i] - data2[i]) is greater than 5 against the null that it
// is <= 5.  Calculate confidence intervals at 99%.
auto result = pairedTTest(summary, 5, Alt.twoSided, 0.99);

// This is equivalent to:
auto result2 = pairedTTest(data1, data2, 5, Alt.twoSided, 0.99);
---

References:  http://en.wikipedia.org/wiki/Student%27s_t-test
 */
ConfInt pairedTTest(T)(
    T diffSummary,
    double testMean = 0,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) if(isSummary!T) {
    enforceConfidence(confLevel);
    dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");

    if(diffSummary.N < 2) {
        return ConfInt.init;
    }

    // Save typing.
    alias diffSummary msd;

    ConfInt ret;
    ret.testStat = (msd.mean - testMean) / msd.stdev * sqrt(msd.N);
    auto sampleMean = msd.mean;
    auto sampleSd = msd.stdev;
    double normSd = sampleSd / sqrt(msd.N);
    ret.testStat = (sampleMean - testMean) / normSd;

    if(alt == Alt.none) {
        return ret;
    } else if(alt == Alt.less) {
        ret.p = studentsTCDF(ret.testStat, msd.N - 1);
        ret.lowerBound = -double.infinity;

        if(confLevel > 0) {
            double delta = invStudentsTCDF(1 - confLevel, msd.N - 1) * normSd;
            ret.upperBound = sampleMean - delta;
        } else {
            ret.upperBound = sampleMean;
        }

    } else if(alt == Alt.greater) {
        ret.p = studentsTCDF(-ret.testStat, msd.N - 1);
        ret.upperBound = double.infinity;

        if(confLevel > 0) {
            double delta = invStudentsTCDF(1 - confLevel, msd.N - 1) * normSd;
            ret.lowerBound = sampleMean + delta;
        } else {
            ret.lowerBound = sampleMean;
        }

    } else {
        immutable double t = ret.testStat;
        ret.p = 2 * ((t < 0) ?
                      studentsTCDF(t, msd.N - 1) :
                      studentsTCDF(-t, msd.N - 1));

        if(confLevel > 0) {
            double delta = invStudentsTCDF(0.5 * (1 - confLevel), msd.N - 1) * normSd;
            ret.lowerBound = sampleMean + delta;
            ret.upperBound = sampleMean - delta;
        } else {
            ret.lowerBound = ret.upperBound = sampleMean;
        }

    }
    return ret;
}

unittest {
    // Values from R.
    auto t1 = pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1);
    assert(approxEqual(t1.p, 0.02131));
    assert(approxEqual(t1.testStat, -3.6742));
    assert(approxEqual(t1.lowerBound, -2.1601748));
    assert(approxEqual(t1.upperBound, 0.561748));

    assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.less).p, 0.0889));
    assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.greater).p, 0.9111));
    assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 0, Alt.twoSided).p, 0.1778));
    assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1, Alt.less).p, 0.01066));
    assert(approxEqual(pairedTTest([3,2,3,4,5], [2,3,5,5,6], 1, Alt.greater).p, 0.9893));
}

/**
Tests the null hypothesis that the variances of all groups are equal against
the alternative that heteroscedasticity exists.  data must be either a
tuple of ranges or a range of ranges.  central is an alias for the measure
of central tendency to be used.  This can be any function that maps a
forward range of numeric types to a numeric type.  The commonly used ones
are median (default) and mean (less robust).  Trimmed mean is sometimes
useful, but is currently not implemented in dstats.summary.

References:
Levene, Howard (1960). "Robust tests for equality of variances". in Ingram
Olkin, Harold Hotelling et al. Contributions to Probability and Statistics:
Essays in Honor of Harold Hotelling. Stanford University Press. pp. 278-292.

Examples:
---
int[] sample1 = [1,2,3,4,5];
int[] sample2 = [100,200,300,400,500];
auto result = levenesTest(sample1, sample2);

// Clearly the variances are different between these two samples.
assert( approxEqual(result.testStat, 10.08));
assert( approxEqual(result.p, 0.01310));
---
*/
TestRes levenesTest(alias central = median, T...)(T data) {
    return anovaLevene!(true, false, central, T)(data);
}

unittest {
    // Values from R's car package, which uses the median definition
    // exclusively.
    auto res1 = levenesTest([1,2,3,4,5][], [2,4,8,16,32][]);
    assert(approxEqual(res1.testStat, 3.0316));
    assert(approxEqual(res1.p, 0.1198), res1.toString());

    auto res2 = levenesTest([[1,2,3,4,5][], [100,200,300,400,500,600][]][]);
    assert(approxEqual(res2.testStat, 13.586));
    assert(approxEqual(res2.p, 0.005029));

    auto res3 = levenesTest([8,6,7,5,3,0,9][], [3,6,2,4,3,6][]);
    assert(approxEqual(res3.testStat, 1.1406));
    assert(approxEqual(res3.p, 0.3084));
}

/**
The F-test is a one-way ANOVA extension of the T-test to >2 groups.
It's useful when you have 3 or more groups with equal variance and want
to test whether their means are equal.  Data can be input as either a
tuple or a range.  This may contain any combination of ranges of numeric
types, MeanSD structs and Summary structs.

Note:  This test makes the assumption that all groups have equal variances,
also known as homoskedasticity.  For a similar test that does not make these
assumptions, see welchAnova.

Examples:
---
uint[] thing1 = [3,1,4,1],
       thing2 = [5,9,2,6,5,3],
       thing3 = [5,8,9,7,9,3];
auto result = fTest(thing1, meanStdev(thing2), summary(thing3));
assert(approxEqual(result.testStat, 4.9968));
assert(approxEqual(result.p, 0.02456));
---

References:  http://en.wikipedia.org/wiki/F-test

Returns:

A TestRes containing the F statistic and the P-value for the alternative
that the means of the groups are different against the null that they
are identical.
*/
TestRes fTest(T...)(T data) {
    return anovaLevene!(false, false, "dummy", T)(data);
}

/**
Same as fTest, except that this test does not require the assumption of
equal variances.  In exchange it's slightly less powerful.

References:

B.L. Welch. On the Comparison of Several Mean Values: An Alternative Approach
Biometrika, Vol. 38, No. 3/4 (Dec., 1951), pp. 330-336.
 */
TestRes welchAnova(T...)(T data) {
    return anovaLevene!(false, true, "dummy", T)(data);
}

unittest {
    // Values from R.
    uint[] thing1 = [3,1,4,1],
           thing2 = [5,9,2,6,5,3],
           thing3 = [5,8,9,7,9,3];
    auto result = fTest(thing1, meanStdev(thing2), summary(thing3));
    assert(approxEqual(result.testStat, 4.9968));
    assert(approxEqual(result.p, 0.02456));

    auto welchRes1 = welchAnova(thing1, thing2, thing3);
    assert( approxEqual(welchRes1.testStat, 6.7813));
    assert( approxEqual(welchRes1.p, 0.01706));

    // Test array case.
    auto res2 = fTest([thing1, thing2, thing3].dup);
    assert(approxEqual(result.testStat, res2.testStat));
    assert(approxEqual(result.p, res2.p));

    thing1 = [2,7,1,8,2];
    thing2 = [8,1,8];
    thing3 = [2,8,4,5,9];
    auto res3 = fTest(thing1, thing2, thing3);
    assert(approxEqual(res3.testStat, 0.377));
    assert(approxEqual(res3.p, 0.6953));

    auto res4 = fTest([summary(thing1), summary(thing2), summary(thing3)][]);
    assert(approxEqual(res4.testStat, res3.testStat));
    assert(approxEqual(res4.testStat, res3.testStat));

    auto welchRes2 = welchAnova(summary(thing1), thing2, thing3);
    assert( approxEqual(welchRes2.testStat, 0.342));
    assert( approxEqual(welchRes2.p, 0.7257));

    auto res5 = fTest([1, 2, 4], [3]);
    assert(approxEqual(res5.testStat, 0.1429));
    assert(approxEqual(res5.p, 0.7418));
}

// Levene's Test, Welch ANOVA and F test have massive overlap at the
// implementation level but less at the conceptual level, so I've combined
// the implementations into one horribly complicated but well-encapsulated
// templated function but left the interfaces as three unrelated functions.
private TestRes anovaLevene(bool levene, bool welch, alias central,  T...)
(T dataIn) {
    static if(dataIn.length == 1) {
        auto alloc = newRegionAllocator();
        auto data = alloc.array(dataIn[0]);
        auto withins = alloc.uninitializedArray!(MeanSD[])(data.length);
        withins[] = MeanSD.init;
    } else {
        enum len = dataIn.length;
        alias dataIn data;
        MeanSD[len] withins;
    }

    static if(levene) {
        static if(dataIn.length == 1) {
            auto centers = alloc.uninitializedArray!(double[])(data.length);
        } else {
            double[len] centers;
        }

        foreach(i, category; data) {
            static assert( isForwardRange!(typeof(category)) &&
                is(Unqual!(ElementType!(typeof(category))) : double),
                "Can only perform Levene's test on input ranges of elements " ~
                "implicitly convertible to doubles.");

            // The cast is to force conversions to double on alias this'd stuff
            // like the Mean struct.
            centers[i] = cast(double) central(category.save);
        }

        double preprocess(double dataPoint, size_t category) {
            return abs(dataPoint - centers[category]);
        }
    } else {
        static double preprocess(double dataPoint, size_t category) {
            return dataPoint;
        }
    }


    auto DFGroups = data.length - 1;
    ulong N = 0;

    foreach(category, range; data) {
        static if(isInputRange!(typeof(range)) &&
            is(Unqual!(ElementType!(typeof(range))) : double)) {
            foreach(elem; range) {
                double preprocessed = preprocess(elem, category);
                withins[category].put(preprocessed);
                N++;
            }
        } else static if(isSummary!(typeof(range))) {
            withins[category] = range.toMeanSD();
            N += roundTo!long(range.N);
        } else {
            static assert(0, "Can only perform ANOVA on input ranges of " ~
                "numeric types, MeanSD structs and Summary structs, not a " ~
                typeof(range).stringof ~ ".");
        }
    }

    static if(!welch) {
        immutable ulong DFDataPoints = N - data.length;
        double mu = 0;
        foreach(summary; withins) {
            mu += summary.mean * (summary.N / N);
        }

        double totalWithin = 0;
        double totalBetween = 0;
        foreach(group; withins) {
            totalWithin += group.mse * (group.N / DFDataPoints);
            immutable diffSq = (group.mean - mu) * (group.mean - mu);
            totalBetween += diffSq * (group.N / DFGroups);
        }

        immutable  F = totalBetween / totalWithin;
        if(isNaN(F)) {
            return TestRes.init;
        }

        return TestRes(F, fisherCDFR(F, DFGroups, DFDataPoints));
    } else {
        immutable double k = data.length;
        double sumW = 0;
        foreach(summary; withins) {
            sumW += summary.N / summary.var;
        }

        double sumFt = 0;
        foreach(summary; withins) {
            sumFt += ((1 - summary.N / summary.var / sumW) ^^ 2) / (summary.N - 1);
        }

        immutable kSqM1 = (k * k - 1.0);
        immutable df2 = 1.0 / (3.0 / kSqM1 * sumFt);
        immutable denom = 1 + 2 * (k - 2.0) / kSqM1 * sumFt;

        double yHat = 0;
        foreach(i, summary; withins) {
            yHat += summary.mean * (summary.N / summary.var);
        }
        yHat /= sumW;

        double numerator = 0;
        foreach(i, summary; withins) {
            numerator += summary.N / summary.var * ((summary.mean - yHat) ^^ 2);
        }
        numerator /= (k - 1);

        immutable F = numerator / denom;
        if(isNaN(F)) {
            return TestRes.init;
        }

        return TestRes(F, fisherCDFR(F, DFGroups, df2));
    }
}

/**Performs a correlated sample (within-subjects) ANOVA.  This is a
 * generalization of the paired T-test to 3 or more treatments.  This
 * function accepts data as either a tuple of ranges (1 for each treatment,
 * such that a given index represents the same subject in each range) or
 * similarly as a range of ranges.
 *
 * Returns:  A TestRes with the F-statistic and P-value for the null that
 * the the variable being measured did not vary across treatments against the
 * alternative that it did.
 *
 * Examples:
 * ---
 * // Test the hypothesis that alcohol, loud music, caffeine and sleep
 * // deprivation all have equivalent effects on programming ability.
 *
 * uint[] alcohol = [8,6,7,5,3,0,9];
 * uint[] caffeine = [3,6,2,4,3,6,8];
 * uint[] noSleep = [3,1,4,1,5,9,2];
 * uint[] loudMusic = [2,7,1,8,2,8,1];
 * // Subject 0 had ability of 8 under alcohol, 3 under caffeine, 3 under
 * // no sleep, 2 under loud music.  Subject 1 had ability of 6 under alcohol,
 * // 6 under caffeine, 1 under no sleep, and 7 under loud music, etc.
 * auto result = correlatedAnova(alcohol, caffeine, noSleep, loudMusic);
 * ---
 *
 * References:  "Concepts and Applications of Inferrential Statistics".
 *              Richard Lowry.  Vassar College.   version.
 *              http://faculty.vassar.edu/lowry/webtext.html
 */
TestRes correlatedAnova(T...)(T dataIn)
if(allSatisfy!(isInputRange, T)) {
    static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
        auto alloc = newRegionAllocator();
        auto data = alloc.array(dataIn[0]);
        auto withins = alloc.newArray!(MeanSD[])(data.length);
    } else {
        enum len = dataIn.length;
        alias dataIn data;
        MeanSD[len] withins;
    }
    MeanSD overallSumm;
    double nGroupNeg1 = 1.0 / data.length;

    bool someEmpty() {
        foreach(elem; data) {
            if(elem.empty) {
                return true;
            }
        }
        return false;
    }

    uint nSubjects = 0;
    double subjSum = 0;
    while(!someEmpty) {
        double subjSumInner = 0;
        foreach(i, elem; data) {
            auto dataPoint = elem.front;
            subjSumInner += dataPoint;
            overallSumm.put(dataPoint);
            withins[i].put(dataPoint);
            data[i].popFront;
        }
        nSubjects++;
        subjSum += subjSumInner * subjSumInner * nGroupNeg1;
    }
    double groupSum = 0;
    foreach(elem; withins) {
        groupSum += elem.mean * elem.N;
    }

    groupSum /= sqrt(cast(double) nSubjects * data.length);
    groupSum *= groupSum;
    immutable subjErr = subjSum - groupSum;

    double betweenDev = 0;
    immutable mu = overallSumm.mean;
    foreach(group; withins) {
        double diff = (group.mean - mu);
        diff *= diff;
        betweenDev += diff * (group.N / (data.length - 1));
    }

    size_t errDf = data.length * nSubjects - data.length - nSubjects + 1;
    double randError = -subjErr / errDf;
    foreach(group; withins) {
        randError += group.mse * (group.N / errDf);
    }

    immutable F = betweenDev / randError;
    if(!(F >= 0)) {
        return TestRes(double.nan, double.nan);
    }

    return TestRes(F, fisherCDFR(F, data.length - 1, errDf));
}

unittest {
    // Values from VassarStats utility at
    // http://faculty.vassar.edu/lowry/VassarStats.html, but they like to
    // round a lot, so the approxEqual tolerances are fairly wide.  I
    // think it's adequate to demonstrate the correctness of this function,
    // though.
    uint[] alcohol = [8,6,7,5,3,0,9];
    uint[] caffeine = [3,6,2,4,3,6,8];
    uint[] noSleep = [3,1,4,1,5,9,2];
    uint[] loudMusic = [2,7,1,8,2,8,1];
    auto result = correlatedAnova(alcohol, caffeine, noSleep, loudMusic);
    assert(approxEqual(result.testStat, 0.43, 0.0, 0.01));
    assert(approxEqual(result.p, 0.734, 0.0, 0.01));

    uint[] stuff1 = [3,4,2,6];
    uint[] stuff2 = [4,1,9,8];
    auto result2 = correlatedAnova([stuff1, stuff2].dup);
    assert(approxEqual(result2.testStat, 0.72, 0.0, 0.01));
    assert(approxEqual(result2.p, 0.4584, 0.0, 0.01));
}

/**The Kruskal-Wallis rank sum test.  Tests the null hypothesis that data in
 * each group is not stochastically ordered with respect to data in each other
 * groups.  This is a one-way non-parametric ANOVA and can be thought of
 * as either a generalization of the Wilcoxon rank sum test to >2 groups or
 * a non-parametric equivalent to the F-test.  Data can be input as either a
 * tuple of ranges (one range for each group) or a range of ranges
 * (one element for each group).
 *
 * Bugs:  Asymptotic approximation of P-value only, not exact.  In this case,
 * I'm not sure a practical way to compute the exact P-value even exists.
 *
 * Returns:  A TestRes with the K statistic and the P-value for the null that
 * no group is stochastically larger than any other against the alternative that
 * groups are stochastically ordered.
 *
 * References:  "Concepts and Applications of Inferrential Statistics".
 *              Richard Lowry.  Vassar College.   version.
 *              http://faculty.vassar.edu/lowry/webtext.html
 */
TestRes kruskalWallis(T...)(T dataIn)
if(doubleInput!(typeof(dataIn[0].front)) || allSatisfy!(doubleInput, T)) {
    auto alloc = newRegionAllocator();
    size_t N = 0;

    static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
        auto data = alloc.array(dataIn[0]);
        alias ElementType!(typeof(data[0])) C;
        static if(hasLength!(typeof(data[0]))) {
            enum bool useLength = true;
        } else {
            enum bool useLength = false;
        }
    } else {
        enum len = dataIn.length;
        alias dataIn data;
        alias staticMap!(ElementType, T) Es;
        alias CommonType!(Es) C;
        static if(allSatisfy!(hasLength, T)) {
            enum bool useLength = true;
        } else {
            enum bool useLength = false;
        }
    }

    size_t[] lengths = alloc.uninitializedArray!(size_t[])(data.length);
    static if(useLength) {
        foreach(i, rng; data) {
            auto rngLen = rng.length;
            lengths[i] = rngLen;
            N += rngLen;
        }
        auto dataArray = alloc.uninitializedArray!(Unqual!(C)[])(N);
        size_t pos = 0;
        foreach(rng; data) {
            foreach(elem; rng) {
                dataArray[pos++] = elem;
            }
        }
    } else {
        auto app = appender!(Unqual!(C)[])();
        foreach(i, rng; data) {
            size_t oldLen = dataArray.length;
            app.put(rng);
            lengths[i] = dataArray.length - oldLen;
            N += lengths[i];
        }
        auto dataArray = app.data;
    }

    double[] ranks = alloc.uninitializedArray!(double[])(dataArray.length);
    try {
        rankSort(dataArray, ranks);
    } catch(SortException) {
        return TestRes.init;
    }

    size_t index = 0;
    double denom = 0, numer = 0;
    double rBar = 0.5 * (N + 1);
    foreach(meanI, l; lengths) {
        Mean groupStats;
        foreach(i; index..index + l) {
            groupStats.put( ranks[i]);
            double diff = ranks[i] - rBar;
            diff *= diff;
            denom += diff;
        }
        index += l;
        double nDiff = groupStats.mean - rBar;
        nDiff *= nDiff;
        numer += l * nDiff;
    }
    double K = (N - 1) * (numer / denom);

    // Tie correction.
    double tieSum = 0;
    uint nTies = 1;
    foreach(i; 1..dataArray.length) {
        if(dataArray[i] == dataArray[i - 1]) {
            nTies++;
        } else if(nTies > 1) {
            double partialSum = nTies;
            partialSum = (partialSum * partialSum * partialSum) - partialSum;
            tieSum += partialSum;
            nTies = 1;
        }
    }
    if(nTies > 1) {
        double partialSum = nTies;
        partialSum = (partialSum * partialSum * partialSum) - partialSum;
        tieSum += partialSum;
    }
    double tieDenom = N;
    tieDenom = (tieDenom * tieDenom * tieDenom) - tieDenom;
    tieSum = 1 - (tieSum / tieDenom);
    K *= tieSum;

    if(isNaN(K)) {
        return TestRes(double.nan, double.nan);
    }

    return TestRes(K, chiSquareCDFR(K, data.length - 1));
}

unittest {
    // These values are from the VassarStat web tool at
    // http://faculty.vassar.edu/lowry/VassarStats.html .
    // R is actually wrong here because it apparently doesn't use a correction
    // for ties.
    auto res1 = kruskalWallis([3,1,4,1].idup, [5,9,2,6].dup, [5,3,5].dup);
    assert(approxEqual(res1.testStat, 4.15));
    assert(approxEqual(res1.p, 0.1256));

    // Test for other input types.
    auto res2 = kruskalWallis([[3,1,4,1].idup, [5,9,2,6].idup, [5,3,5].idup].dup);
    assert(res2 == res1);
    auto res3 = kruskalWallis(map!"a"([3,1,4,1].dup), [5,9,2,6].dup, [5,3,5].dup);
    assert(res3 == res1);
    auto res4 = kruskalWallis([map!"a"([3,1,4,1].dup),
                               map!"a"([5,9,2,6].dup),
                               map!"a"([5,3,5].dup)].dup);
    assert(res4 == res1);

    // Test w/ one more case, just with one input type.
    auto res5 = kruskalWallis([2,7,1,8,2].dup, [8,1,8,2].dup, [8,4,5,9,2].dup,
                              [7,1,8,2,8,1,8].dup);
    assert(approxEqual(res5.testStat, 1.06));
    assert(approxEqual(res5.p, 0.7867));
}

/**The Friedman test is a non-parametric within-subject ANOVA.  It's useful
 * when parametric assumptions cannot be made.  Usage is identical to
 * correlatedAnova().
 *
 * References:  "Concepts and Applications of Inferrential Statistics".
 *              Richard Lowry.  Vassar College.   version.
 *              http://faculty.vassar.edu/lowry/webtext.html
 *
 * Bugs:  No exact P-value calculation.  Asymptotic approx. only.
 */
TestRes friedmanTest(T...)(T dataIn)
if(doubleInput!(typeof(dataIn[0].front)) || allSatisfy!(doubleInput, T)) {
    static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
        auto alloc = newRegionAllocator();
        auto data = alloc.array(dataIn[0]);
        auto ranks = alloc.uninitializedArray!(double[])(data.length);
        auto dataPoints = alloc.uninitializedArray!(double[])(data.length);
        auto colMeans = alloc.newArray!(Mean[])(data.length);
    } else {
        enum len = dataIn.length;
        alias dataIn data;
        double[len] ranks;
        double[len] dataPoints;
        Mean[len] colMeans;
    }
    double rBar = cast(double) data.length * (data.length + 1.0) / 2.0;
    MeanSD overallSumm;

    bool someEmpty() {
        foreach(elem; data) {
            if(elem.empty) {
                return true;
            }
        }
        return false;
    }

    uint N = 0;
    while(!someEmpty) {
        foreach(i, range; data) {
            dataPoints[i] = data[i].front;
            data[i].popFront;
        }

        try {
            rank(dataPoints[], ranks[]);
        } catch(SortException) {
            return TestRes.init;
        }

        foreach(i, rank; ranks) {
            colMeans[i].put(rank);
            overallSumm.put(rank);
        }
        N++;
    }

    double between = 0;
    double mu = overallSumm.mean;
    foreach(mean; colMeans) {
        double diff = mean.mean - overallSumm.mean;
        between += diff * diff;
    }
    between *= N;
    double within = overallSumm.mse * (overallSumm.N / (overallSumm.N - N));
    double chiSq = between / within;
    double df = data.length - 1;
    return TestRes(chiSq, chiSquareCDFR(chiSq, df));
}

unittest {
    // Values from R
    uint[] alcohol = [8,6,7,5,3,0,9];
    uint[] caffeine = [3,6,2,4,3,6,8];
    uint[] noSleep = [3,1,4,1,5,9,2];
    uint[] loudMusic = [2,7,1,8,2,8,1];
    auto result = friedmanTest(alcohol, caffeine, noSleep, loudMusic);
    assert(approxEqual(result.testStat, 1.7463));
    assert(approxEqual(result.p, 0.6267));

    uint[] stuff1 = [3,4,2,6];
    uint[] stuff2 = [4,1,9,8];
    auto result2 = friedmanTest([stuff1, stuff2].dup);
    assert(approxEqual(result2.testStat, 1));
    assert(approxEqual(result2.p, 0.3173));
}

/**Computes Wilcoxon rank sum test statistic and P-value for
 * a set of observations against another set, using the given alternative.
 * Alt.less means that sample1 is stochastically less than sample2.
 * Alt.greater means sample1 is stochastically greater than sample2.
 * Alt.twoSided means sample1 is stochastically less than or greater than
 * sample2.
 *
 * exactThresh is the threshold value of (n1 + n2) at which this function
 * switches from exact to approximate computation of the p-value.  Do not set
 * exactThresh to more than 200, as the exact
 * calculation is both very slow and not numerically stable past this point,
 * and the asymptotic calculation is very good for N this large.  To disable
 * exact calculation entirely, set exactThresh to 0.
 *
 * Notes:  Exact p-value computation is never used when ties are present in the
 * data, because it is not computationally feasible.
 *
 * Input ranges for this function must define a length.
 *
 * This test is also known as the Mann-Whitney U test.
 *
 * Returns:  A TestRes containing the W test statistic and the P-value against
 * the given alternative.
 *
 * References:  http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U
 *
 * StackOverflow Question 376003  http://stackoverflow.com/questions/376003
 *
 * Loughborough University MLSC Statistics 2.3 The Mann-Whitney U Test
 * http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf
 */
TestRes wilcoxonRankSum(T, U)
(T sample1, U sample2, Alt alt = Alt.twoSided, uint exactThresh = 50)
if(isInputRange!T && isInputRange!U &&
is(typeof(sample1.front < sample2.front) == bool) &&
is(CommonType!(ElementType!T, ElementType!U))) {

    auto alloc = newRegionAllocator();
    alias Unqual!(CommonType!(ElementType!(T), ElementType!(U))) C;

    static if(hasLength!T && hasLength!U) {
        auto n1 = sample1.length, n2 = sample2.length, N = n1 + n2;
        auto combined = alloc.uninitializedArray!(C[])(N);
        copy(sample1, combined[0..n1]);
        copy(sample2, combined[n1..$]);
    } else {
        auto app = appender!(C[])();

        foreach(elem; sample1) {
            app.put(elem);
        }

        uint n1 = app.data.length;
        foreach(elem; sample2) {
            app.put(elem);
        }

        auto combined = app.data;
        uint N = combined.length;
        uint n2 = N - n1;
    }

    double[] ranks = alloc.uninitializedArray!(double[])(N);
    try {
        rankSort(combined, ranks);
    } catch(SortException) {
        return TestRes.init;
    }
    double w = reduce!("a + b")
        (0.0, ranks[0..n1]) - cast(ulong) n1 * (n1 + 1) / 2UL;

    if(alt == Alt.none) {
        return TestRes(w);
    }

    double tieSum = 0;
    // combined is sorted by rankSort.  Can use it to figure out how many
    // ties we have w/o another allocation or sorting.
    enum oneOverTwelve = 1.0 / 12.0;
    tieSum = 0;
    ulong nties = 1;
    foreach(i; 1..N) {
        if(combined[i] == combined[i - 1]) {
            nties++;
        } else {
            if(nties == 1)
                continue;
            tieSum += ((nties * nties * nties) - nties) * oneOverTwelve;
            nties = 1;
        }
    }
    // Handle last run.
    if(nties > 1) {
        tieSum += ((nties * nties * nties) - nties) * oneOverTwelve;
    }

    immutable p = wilcoxonRankSumPval(w, n1, n2, alt, tieSum, exactThresh);
    return TestRes(w, p);
}

 unittest {
     // Values from R.

    assert(wilcoxonRankSum([1, 2, 3, 4, 5].dup, [2, 4, 6, 8, 10].dup).testStat == 5);
    assert(wilcoxonRankSum([2, 4, 6, 8, 10].dup, [1, 2, 3, 4, 5].dup).testStat == 20);
    assert(wilcoxonRankSum([3, 7, 21, 5, 9].dup, [2, 4, 6, 8, 10].dup).testStat == 15);

     // Simple stuff (no ties) first.  Testing approximate
     // calculation first.
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
           Alt.twoSided, 0), 0.9273));
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
           Alt.less, 0), 0.6079));
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
           Alt.greater, 0).p, 0.4636));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
            Alt.twoSided, 0).p, 0.4113));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
            Alt.less, 0).p, 0.2057));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup,
        map!"a"([3,5,7,8,13,15].dup), Alt.greater, 0).p, 0.8423));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.twoSided, 0), .6745));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.less, 0), .3372));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.greater, 0), .7346));

    // Now, lots of ties.
    assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
           Alt.twoSided, 0), 0.3976));
    assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
           Alt.less, 0), 0.1988));
    assert(approxEqual(wilcoxonRankSum([1,2,3,4,5].dup, [2,3,4,5,6].dup,
           Alt.greater, 0), 0.8548));
    assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
           Alt.twoSided, 0), 0.9049));
    assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
           Alt.less, 0), 0.4524));
    assert(approxEqual(wilcoxonRankSum([1,2,1,1,2].dup, [1,2,3,1,1].dup,
           Alt.greater, 0), 0.64));

    // Now, testing the exact calculation on the same data.
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
       Alt.twoSided), 0.9307));
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
           Alt.less), 0.6039));
     assert(approxEqual(wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup,
           Alt.greater), 0.4654));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
            Alt.twoSided), 0.4286));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
            Alt.less), 0.2143));
     assert(approxEqual(wilcoxonRankSum([1,2,6,10,12].dup, [3,5,7,8,13,15].dup,
            Alt.greater), 0.8355));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.twoSided), .6905));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.less), .3452));
     assert(approxEqual(wilcoxonRankSum([1,3,5,7,9].dup, [2,4,6,8,10].dup,
            Alt.greater), .7262));
}

private
double wilcoxonRankSumPval(double w, ulong n1, ulong n2, Alt alt = Alt.twoSided,
                           double tieSum = 0,  uint exactThresh = 50) {
    if(alt == Alt.none) {
        return double.nan;
    }

    immutable double N = n1 + n2;

    if(N < exactThresh && tieSum == 0) {
        return wilcoxRSPExact(roundTo!uint(w), cast(uint) n1, cast(uint) n2, alt);
    }

    immutable sd = sqrt(cast(double) (n1 * n2) / (N * (N - 1)) *
             ((N * N * N - N) / 12 - tieSum));

    // Can happen if all samples are tied.
    if(!(sd > 0)) {
        return double.nan;
    }

    immutable mean = (n1 * n2) / 2.0;

    if(alt == Alt.twoSided) {
        if(abs(w - mean) < 0.5) {
            return 1;
        } else if(w < mean) {
            return 2 * normalCDF(w + 0.5, mean, sd);
        } else {
            assert(w > mean);
            return 2 * normalCDFR(w - 0.5, mean, sd);
        }
    } else if(alt == Alt.less) {
        return normalCDF(w + 0.5, mean, sd);
    } else if(alt == Alt.greater) {
        return normalCDFR(w - 0.5, mean, sd);
    }

    assert(0);
}

unittest {
    /* Values from R.  I could only get good values for Alt.less directly.
     * Using W-values to test Alt.twoSided, Alt.greater indirectly.*/
    assert(approxEqual(wilcoxonRankSumPval(1200, 50, 50, Alt.less), .3670));
    assert(approxEqual(wilcoxonRankSumPval(1500, 50, 50, Alt.less), .957903));
    assert(approxEqual(wilcoxonRankSumPval(8500, 100, 200, Alt.less), .01704));
    auto w = wilcoxonRankSum([2,4,6,8,12].dup, [1,3,5,7,11,9].dup).testStat;
    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6), 0.9273));
    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.greater), 0.4636));
    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.less), 0.6079));

    // Monte carlo unit testing:  Make sure that the exact and asymptotic
    // versions agree within a small epsilon;
    double maxEpsilon = 0;
    foreach(i; 0..1_000) {
        uint n1 = uniform(5U, 25U);
        uint n2 = uniform(5U, 25U);
        uint testStat = uniform!"[]"(0, (n1 * n2));

        foreach(alt; [Alt.less, Alt.greater, Alt.twoSided]) {
            double approxP = wilcoxonRankSumPval(testStat, n1, n2, alt, 0, 0);
            double exactP = wilcoxonRankSumPval(testStat, n1, n2, alt, 0, 50);
            double epsilon = abs(approxP - exactP);
            assert(epsilon < 0.02);
            maxEpsilon = max(maxEpsilon, epsilon);
        }
    }
}

/* Used internally by wilcoxonRankSum.  This function uses dynamic
 * programming to count the number of combinations of numbers [1..N] that sum
 * of length n1 that sum to <= W in O(N * W * n1) time.
 * Algorithm obtained from StackOverflow Question 376003
 * (http://stackoverflow.com/questions/376003).*/
private double wilcoxRSPExact(uint W, uint n1, uint n2, Alt alt = Alt.twoSided) {
    uint N = n1 + n2;
    immutable maxPossible = n1 * n2;

    switch(alt) {
        case Alt.less:
            if(W >= maxPossible)  { // Value impossibly large
                return 1;
            } else if(W * 2 <= maxPossible) {
                break;
            } else {
                return 1 - wilcoxRSPExact(maxPossible - W - 1, n1, n2, Alt.less);
            }
            assert(0);
        case Alt.greater:
            if(W > maxPossible)  { // Value impossibly large
                return 0;
            } else if(W * 2 >= maxPossible) {
                return wilcoxRSPExact(maxPossible - W, n1, n2, Alt.less);
            } else if(W <= 0) {
                return 1;
            } else {
                return 1 - wilcoxRSPExact(W - 1, n1, n2, Alt.less);
            }
            assert(0);
        case Alt.twoSided:
            if(W * 2 <= maxPossible) {
                return min(1, wilcoxRSPExact(W, n1, n2, Alt.less) +
                       wilcoxRSPExact(maxPossible - W, n1, n2, Alt.greater));
            } else {
                return min(1, wilcoxRSPExact(W, n1, n2, Alt.greater) +
                       wilcoxRSPExact(maxPossible - W, n1, n2, Alt.less));
            }
            assert(0);
        default:
            assert(0);
    }

    W += n1 * (n1 + 1) / 2UL;

    auto alloc = newRegionAllocator();
    float* cache = (alloc.uninitializedArray!(float[])((n1 + 1) * (W + 1))).ptr;
    float* cachePrev = (alloc.uninitializedArray!(float[])((n1 + 1) * (W + 1))).ptr;
    cache[0..(n1 + 1) * (W + 1)] = 0;
    cachePrev[0..(n1 + 1) * (W + 1)] = 0;

    /* Using doubles for the intermediate steps is too slow, but I didn't want to
     * lose too much precision.  Since my sums must be between 0 and 1, I am
     * using the entire bit space of a float to hold numbers between zero and
     * one.  This is precise to at least 1e-7.  This is good enough for a few
     * reasons:
     *
     * 1.  This is a p-value, and therefore will likely not be used in
     *     further calculations where rounding error would accumulate.
     * 2.  If this is too slow, the alternative is to use the asymptotic
     *     approximation.  This is can have relative errors of several orders
     *     of magnitude in the tails of the distribution, and is therefore
     *     clearly worse.
     * 3.  For very large N, where this function could give completely wrong
     *     answers, it would be so slow that any reasonable person would use the
     *     asymptotic approximation anyhow.*/


    // Algorithm based on StackOverflow question 376003.
    double comb = exp(-logNcomb(N, n1));
    double floatMax = cast(double) float.max;
    cache[0] = cast(float) (comb * floatMax);
    cachePrev[0] = cast(float) (comb * floatMax);

    foreach(i; 1..N + 1) {
        swap(cache, cachePrev);
        foreach(k; 1..min(i + 1, n1 + 1)) {

            uint minW = k * (k + 1) / 2;
            float* curK = cache + k * (W + 1);
            float* prevK = cachePrev + k * (W + 1);
            float* prevKm1 = cachePrev + (k - 1) * (W + 1);

            foreach(w; minW..W + 1) {
                curK[w] = prevK[w] + ((i <= w) ? prevKm1[w - i] : 0);
            }
        }
    }

    double sum = 0;
    float* lastLine = cache + n1 * (W + 1);
    foreach(w; 1..W + 1) {
        sum += (cast(double) lastLine[w] / floatMax);
    }

    return sum;
}

unittest {
    // Values from R.
    assert(approxEqual(wilcoxRSPExact(14, 5, 6), 0.9307));
    assert(approxEqual(wilcoxRSPExact(14, 5, 6, Alt.less), 0.4654));
    assert(approxEqual(wilcoxRSPExact(14, 5, 6, Alt.greater), 0.6039));
    assert(approxEqual(wilcoxRSPExact(16, 6, 5), 0.9307));
    assert(approxEqual(wilcoxRSPExact(16, 6, 5, Alt.less), 0.6039));
    assert(approxEqual(wilcoxRSPExact(16, 6, 5, Alt.greater), 0.4654));
    assert(approxEqual(wilcoxRSPExact(66, 10, 35, Alt.less), 0.001053));
    assert(approxEqual(wilcoxRSPExact(78, 13, 6, Alt.less), 1));

    // Mostly to make sure that underflow doesn't happen until
    // the N's are truly unreasonable:
    //assert(approxEqual(wilcoxRSPExact(6_000, 120, 120, Alt.less), 0.01276508));
}

/**Computes a test statistic and P-value for a Wilcoxon signed rank test against
 * the given alternative. Alt.less means that elements of before are stochastically
 * less than corresponding elements of after.  Alt.greater means elements of
 * before are stochastically greater than corresponding elements of after.
 * Alt.twoSided means there is a significant difference in either direction.
 *
 * exactThresh is the threshold value of before.length at which this function
 * switches from exact to approximate computation of the p-value.   Do not set
 * exactThresh to more than 200, as the exact calculation is both very slow and
 * not numerically stable past this point, and the asymptotic calculation is
 * very good for N this large.  To disable exact calculation entirely, set
 * exactThresh to 0.
 *
 * Notes:  Exact p-value computation is never used when ties are present,
 * because it is not computationally feasible.
 *
 * The input ranges for this function must define a length and must be
 * forward ranges.
 *
 * Returns:  A TestRes of the W statistic and the p-value against the given
 * alternative.
 *
 * References:  http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test
 *
 * StackOverflow Question 376003  http://stackoverflow.com/questions/376003
 *
 * Handbook of Parametric and nonparametric statistical procedures. David Sheskin.
 * Third Edition. (2004)  CRC Press. Pg. 616.
 */
TestRes wilcoxonSignedRank(T, U)(T before, U after, Alt alt = Alt.twoSided, uint exactThresh = 50)
if(doubleInput!(T) && doubleInput!(U) &&
is(typeof(before.front - after.front) : double)) {
    uint nZero = 0;
    byte sign(double input) {
        if(input < 0)
            return -1;
        if(input > 0)
            return 1;
        nZero++;
        return 0;
    }

    auto alloc = newRegionAllocator();

    static if(hasLength!T && hasLength!U) {
        dstatsEnforce(before.length == after.length,
            "Ranges must have same lengths for wilcoxonSignedRank.");

        double[] diffRanks = alloc.uninitializedArray!(double[])(before.length);
        byte[] signs = alloc.uninitializedArray!(byte[])(before.length);
        double[] diffs = alloc.uninitializedArray!(double[])(before.length);

        size_t ii = 0;
        while(!before.empty && !after.empty) {
            double diff = cast(double) before.front - cast(double) after.front;
            signs[ii] = sign(diff);
            diffs[ii] = abs(diff);
            ii++;
            before.popFront;
            after.popFront;
        }
    } else {
        double[] diffRanks;
        auto diffApp = appender!(double[])();
        auto signApp = appender!(byte[])();

        while(!before.empty && !after.empty) {
            double diff = cast(double) before.front - cast(double) after.front;
            signApp.put(sign(diff));
            diffApp.put(abs(diff));
            before.popFront;
            after.popFront;
        }

        auto diffs = diffApp.data;
        auto signs = signApp.data;
        diffRanks = alloc.uninitializedArray!(double[])(diffs.length);
    }
    try {
        rankSort(diffs, diffRanks);
    } catch(SortException) {
        return TestRes.init;
    }

    ulong N = diffs.length - nZero;

    double W = 0;
    foreach(i, dr; diffRanks) {
        if(signs[i] == 1) {
            W += dr - nZero;
        }
    }

    // Just a sanity check.  Should be mathematically impossible for this
    // assert to fail.  The 1e-5 is for round-off error.
    assert(W > -1e-5 && W <= (N * (N + 1) / 2) + 1e-5);

    if(alt == Alt.none) {
        return TestRes(W);
    }

    // Handle ties.
    double tieSum = 0;

    // combined is sorted by rankSort.  Can use it to figure out how many
    // ties we have w/o another allocation or sorting.
    enum denom = 1.0 / 48.0;
    ulong nties = 1;
    foreach(i; 1..diffs.length) {
        if(diffs[i] == diffs[i - 1] && diffs[i] != 0) {
            nties++;
        } else {
            if(nties == 1)
                continue;
            tieSum += ((nties * nties * nties) - nties) * denom;
            nties = 1;
        }
    }
    // Handle last run.
    if(nties > 1) {
        tieSum += ((nties * nties * nties) - nties) * denom;
    }
    if(nZero > 0 && tieSum == 0) {
        tieSum = double.nan;  // Signal that there were zeros and exact p-val can't be computed.
    }

    return TestRes(W, wilcoxonSignedRankPval(W, N, alt, tieSum, exactThresh));
}

unittest {
    // Values from R.
    alias approxEqual ae;
    assert(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup).testStat == 7.5);
    assert(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup).testStat == 6);
    assert(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup).testStat == 5);

    // With ties, normal approx.
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup), 1));
    assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, map!"a"([2,7,1,8,2].dup)), 0.7865));
    assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup), 0.5879));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.less), 0.5562));
    assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.less), 0.3932));
    assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.less), 0.2940));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.greater), 0.5562));
    assert(ae(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.greater), 0.706));
    assert(ae(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.greater), 0.7918));
    assert(ae(wilcoxonSignedRank(cast(int[]) [1,16,2,4,8], cast(int[]) [1,5,2,3,4]).testStat, 6));
    assert(ae(wilcoxonSignedRank(cast(int[]) [1,16,2,4,8], cast(int[]) [1,5,2,3,4]), 0.1814));

    // Exact.
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup), 0.625));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.less), 0.3125));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.greater), 0.7812));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup), 0.8125));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.less), 0.6875));
    assert(ae(wilcoxonSignedRank([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.greater), 0.4062));

    // Monte carlo unit testing.  Make sure exact, approx are really,
    // really close to each other.
    double maxEpsilon = 0;
    foreach(i; 0..1_000) {
        uint N = uniform(10U, 50U);
        uint testStat = uniform!"[]"(0, N * (N + 1) / 2);

        foreach(alt; [Alt.less, Alt.greater, Alt.twoSided]) {
            double approxP = wilcoxonSignedRankPval(testStat, N, alt, 0, 0);
            double exactP = wilcoxonSignedRankPval(testStat, N, alt, 0, 50);
            double epsilon = abs(approxP - exactP);
            assert(epsilon < 0.02);
            maxEpsilon = max(maxEpsilon, epsilon);
        }
    }
}

/**Same as the overload, but allows testing whether a range is stochastically
 * less than or greater than a fixed value mu rather than paired elements of
 * a second range.*/
TestRes wilcoxonSignedRank(T)(T data, double mu, Alt alt = Alt.twoSided, uint exactThresh = 50)
if(doubleInput!(T) && is(typeof(data.front - mu) : double)) {
    return wilcoxonSignedRank(data, repeat(mu, data.length), alt, exactThresh);
}

unittest {
    auto res = wilcoxonSignedRank([-8,-6,2,4,7].dup, 0);
    assert(approxEqual(res.testStat, 7));
    assert(approxEqual(res.p, 1));
}

private double wilcoxonSignedRankPval(double W, ulong N, Alt alt = Alt.twoSided,
     double tieSum = 0, uint exactThresh = 50)
in {
    assert(N > 0);
    assert(tieSum >= 0 || isNaN(tieSum));
} body {
    if(alt == Alt.none) {
        return double.nan;
    }

    if(tieSum == 0 && !isNaN(tieSum) && N <= exactThresh) {
        return wilcoxSRPExact(roundTo!uint(W), to!uint(N), alt);
    }

    if(isNaN(tieSum)) {
        tieSum = 0;
    }

    immutable expected = N * (N + 1) * 0.25;
    immutable sd = sqrt(N * (N + 1) * (2 * N + 1) / 24.0 - tieSum);

    if(alt == Alt.less) {
        return normalCDF(W + 0.5, expected, sd);
    } else if(alt == Alt.greater) {
        return normalCDFR(W - 0.5, expected, sd);
    } else {
        assert(alt == Alt.twoSided);
        if(abs(W - expected) <= 0.5) {
            return 1;
        } else if(W < expected) {
            return 2 * normalCDF(W + 0.5, expected, sd);
        } else {
            assert(W > expected);
            return 2 * normalCDFR(W - 0.5, expected, sd);
        }
    }
}
// Tested indirectly through other overload.

/* Yes, a little cut and paste coding was involved here from wilcoxRSPExact,
 * but this function and wilcoxRSPExact are just different enough that
 * it would be more trouble than it's worth to write one generalized
 * function.
 *
 * Algorithm adapted from StackOverflow question 376003
 * (http://stackoverflow.com/questions/376003).
 */
private double wilcoxSRPExact(uint W, uint N, Alt alt = Alt.twoSided) {
    immutable maxPossible = N * (N + 1) / 2;

    switch(alt) {
        case Alt.less:
            if(W >= maxPossible)  { // Value impossibly large
                return 1;
            } else if(W * 2 <= maxPossible) {
                break;
            } else {
                return 1 - wilcoxSRPExact(maxPossible - W - 1, N, Alt.less);
            }
        case Alt.greater:
            if(W > maxPossible)  { // Value impossibly large
                return 0;
            } else if(W == 0) {
                return 1;
            } else if(W * 2 >= maxPossible) {
                return wilcoxSRPExact(maxPossible - W, N, Alt.less);
            } else {
                return 1 - wilcoxSRPExact(W - 1, N, Alt.less);
            }
        case Alt.twoSided:
            if(W * 2 <= maxPossible) {
                return min(1, wilcoxSRPExact(W, N, Alt.less) +
                       wilcoxSRPExact(maxPossible - W, N, Alt.greater));
            } else {
                return min(1, wilcoxSRPExact(W, N, Alt.greater) +
                       wilcoxSRPExact(maxPossible - W, N, Alt.less));
            }
        default:
            assert(0);
    }

    auto alloc = newRegionAllocator();
    float* cache = (alloc.uninitializedArray!(float[])((N + 1) * (W + 1))).ptr;
    float* cachePrev = (alloc.uninitializedArray!(float[])((N + 1) * (W + 1))).ptr;
    cache[0..(N + 1) * (W + 1)] = 0;
    cachePrev[0..(N + 1) * (W + 1)] = 0;

    double comb = pow(2.0, -(cast(double) N));
    double floatMax = cast(double) float.max;
    cache[0] = cast(float) (comb * floatMax);
    cachePrev[0] = cast(float) (comb * floatMax);

    foreach(i; 1..N + 1) {
        swap(cache, cachePrev);
        foreach(k; 1..i + 1) {

            uint minW = k * (k + 1) / 2;
            float* curK = cache + k * (W + 1);
            float* prevK = cachePrev + k * (W + 1);
            float* prevKm1 = cachePrev + (k - 1) * (W + 1);

            foreach(w; minW..W + 1) {
                curK[w] = prevK[w] + ((i <= w) ? prevKm1[w - i] : 0);
            }
        }
    }

    double sum  = 0;
    foreach(elem; cache[0..(N + 1) * (W + 1)]) {
        sum += cast(double) elem / (cast(double) float.max);
    }

    return sum;
}

unittest {
    // Values from R.
    assert(approxEqual(wilcoxSRPExact(25, 10, Alt.less), 0.4229));
    assert(approxEqual(wilcoxSRPExact(25, 10, Alt.greater), 0.6152));
    assert(approxEqual(wilcoxSRPExact(25, 10, Alt.twoSided), 0.8457));
    assert(approxEqual(wilcoxSRPExact(31, 10, Alt.less), 0.6523));
    assert(approxEqual(wilcoxSRPExact(31, 10, Alt.greater), 0.3848));
    assert(approxEqual(wilcoxSRPExact(31, 10, Alt.twoSided), 0.7695));
}

/**Sign test for differences between paired values.  This is a very robust
 * but very low power test.  Alternatives are Alt.less, meaning elements
 * of before are typically less than corresponding elements of after,
 * Alt.greater, meaning elements of before are typically greater than
 * elements of after, and Alt.twoSided, meaning that there is a significant
 * difference in either direction.
 *
 * Returns:  A TestRes with the proportion of elements of before that were
 * greater than the corresponding element of after, and the P-value against
 * the given alternative.
 */
TestRes signTest(T, U)(T before, U after, Alt alt = Alt.twoSided)
if(doubleInput!(T) && doubleInput!(U) &&
is(typeof(before.front < after.front) == bool)) {
    ulong greater, less;
    while(!before.empty && !after.empty) {
        if(before.front < after.front) {
            less++;
        } else if(after.front < before.front) {
            greater++;
        }

        // Ignore equals.
        before.popFront;
        after.popFront;
    }

    double propGreater = to!double(greater) / (greater + less);

    final switch(alt) {
        case Alt.none:
            return TestRes(propGreater);
        case Alt.less:
            return TestRes(propGreater,
                binomialCDF(greater, less + greater, 0.5));
        case Alt.greater:
            return TestRes(propGreater,
                binomialCDF(less, less + greater, 0.5));
        case Alt.twoSided:
            if(less > greater) {
                return TestRes(propGreater,
                    2 * binomialCDF(greater, less + greater, 0.5));
            } else if(greater > less) {
                return  TestRes(propGreater,
                    2 * binomialCDF(less, less + greater, 0.5));
            } else {
                return TestRes(propGreater, 1);
            }
    }
}

unittest {
    alias approxEqual ae;
    assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup), 1));
    assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup, Alt.less), 0.5));
    assert(ae(signTest([1,3,4,2,5].dup, [1,2,4,8,16].dup, Alt.greater), 0.875));
    assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup, Alt.greater), 0.03125));
    assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup, Alt.less), 1));
    assert(ae(signTest([5,3,4,6,8].dup, [1,2,3,4,5].dup), 0.0625));

    assert(approxEqual(signTest([1,2,6,7,9].dup, 2), 0.625));
    assert(ae(signTest([1,2,6,7,9].dup, 2).testStat, 0.75));
}

/**Similar to the overload, but allows testing for a difference between a
 * range and a fixed value mu.*/
TestRes signTest(T)(T data, double mu, Alt alt = Alt.twoSided)
if(doubleInput!(T) && is(typeof(data.front < mu) == bool)) {
    return signTest(data, repeat(mu), alt);
}

/**
Two-sided binomial test for whether P(success) == p.  The one-sided
alternatives are covered by dstats.distrib.binomialCDF and binomialCDFR.
k is the number of successes observed, n is the number of trials, p
is the probability of success under the null.

Returns:  The P-value for the alternative that P(success) != p against
the null that P(success) == p.

Notes:  This test can also be performed using multinomialTest, but this
implementation is much faster and easier to use.
 */
double binomialTest(ulong k, ulong n, double p) {
    dstatsEnforce(k <= n, "k must be <= n for binomial test.");
    dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial test.");

    enum epsilon = 1 - 1e-6;  // Small but arbitrary constant to deal w/ rounding error.

    immutable mode = cast(long) ((n + 1) * p);
    if(k == mode ||
       approxEqual(binomialPMF(k, n, p), binomialPMF(mode, n, p), 1 - epsilon)) {
        return 1;
    } else if(k > mode) {
        immutable double upperPart = binomialCDFR(k, n, p);
        immutable pExact = binomialPMF(k, n, p);
        ulong ulim = mode, llim = 0, guess;
        while(ulim - llim > 1) {
            guess = (ulim + llim) / 2;
            immutable double pGuess = binomialPMF(guess, n, p);

            if(pGuess == pExact) {
                ulim = guess + 1;
                llim = guess;
                break;
            } else if(pGuess < pExact) {
                llim = guess;
            } else {
                ulim = guess;
            }
        }

        guess = ulim;
        while(binomialPMF(guess, n, p) < pExact * epsilon) {
            guess++;
        }
        while(guess > 0 && binomialPMF(guess, n, p) > pExact / epsilon) {
            guess--;
        }
        if(guess == 0 && binomialPMF(0, n, p) > pExact / epsilon) {
            return upperPart;
        }
        return upperPart + binomialCDF(guess, n, p);
    } else {
        static double myPMF(ulong k, ulong n, double p) {
            return k > n ? 0 : binomialPMF(k, n, p);
        }

        immutable lowerPart = binomialCDF(k, n, p);
        immutable pExact = binomialPMF(k, n, p);
        ulong ulim = n + 1, llim = mode, guess;
        while(ulim - llim > 1) {
            guess = (ulim + llim) / 2;
            immutable double pGuess = myPMF(guess, n, p);
            if(pGuess == pExact) {
                ulim = guess;
                llim = guess;
                break;
            } else if(pGuess < pExact) {
                ulim = guess;
            } else {
                llim = guess;
            }
        }

        // All this stuff is necessary to deal with round-off error properly.
        guess = llim;
        while(myPMF(guess, n, p) < pExact * epsilon && guess > 0) {
            guess--;
        }
        while(myPMF(guess, n, p) > pExact / epsilon) {
            guess++;
        }

        return lowerPart + ((guess > n) ? 0 : binomialCDFR(guess, n, p));
    }
}

unittest {
    // Values from R.
    assert(approxEqual(binomialTest(46, 96, 0.5), 0.759649));
    assert(approxEqual(binomialTest(44, 56, 0.5), 2.088e-5));
    assert(approxEqual(binomialTest(12, 56, 0.5), 2.088e-5));
    assert(approxEqual(binomialTest(0, 40, 0.25), 2.236e-5));
    assert(approxEqual(binomialTest(5, 16, 0.5), 0.2101));
    assert(approxEqual(binomialTest(0, 20, 0.4), 4.16e-5));
    assert(approxEqual(binomialTest(20, 20, 0.6), 4.16e-5));
    assert(approxEqual(binomialTest(6, 88, 0.1), 0.3784));
    assert(approxEqual(binomialTest(3, 4, 0.5), 0.625));
    assert(approxEqual(binomialTest(4, 7, 0.8), 0.1480));
    assert(approxEqual(binomialTest(3, 9, 0.8), 0.003066));
    assert(approxEqual(binomialTest(9, 9, 0.7), 0.06565));
    assert(approxEqual(binomialTest(2, 11, 0.1), 0.3026));
    assert(approxEqual(binomialTest(1, 11, 0.1), 1));
    assert(approxEqual(binomialTest(5, 11, 0.1), 0.002751));
    assert(approxEqual(binomialTest(5, 12, 0.5), 0.7744));
    assert(approxEqual(binomialTest(12, 12, 0.5), 0.0004883));
    assert(approxEqual(binomialTest(12, 13, 0.6), 0.02042));
    assert(approxEqual(binomialTest(0, 9, 0.1), 1));
}

///For chiSquareFit and gTestFit, is expected value range counts or proportions?
enum Expected {
    ///
    count,

    ///
    proportion
}

/**Performs a one-way Pearson's chi-square goodness of fit test between a range
of observed and a range of expected values.  This is a useful statistical
test for testing whether a set of observations fits a discrete distribution.

Returns:  A TestRes of the chi-square statistic and the P-value for the
alternative hypothesis that observed is not a sample from expected against
the null that observed is a sample from expected.

Notes:  By default, expected is assumed to be a range of expected proportions.
These proportions are automatically normalized, and can sum to any number.
By passing Expected.count in as the last parameter, calculating expected
counts will be skipped, and expected will assume to already be properly
normalized.  This is slightly faster, but more importantly
allows input ranges to be used.

The chi-square test relies on asymptotic statistical properties
and is therefore not considered valid, as a rule of thumb,  when expected
counts are below 5.  However, this rule is likely to be unnecessarily
stringent in most cases.

Examples:
---
// Test to see whether a set of categorical observations differs
// statistically from a discrete uniform distribution.

uint[] observed = [980, 1028, 1001, 964, 1102];
auto expected = repeat(1.0);
auto res2 = chiSquareFit(observed, expected);
assert(approxEqual(res2, 0.0207));
assert(approxEqual(res2.testStat, 11.59));
---
 *
References:  http://en.wikipedia.org/wiki/Pearson%27s_chi-square_test
 */
TestRes chiSquareFit(T, U)(
    T observed,
    U expected,
    Expected countProp = Expected.proportion
) if(doubleInput!(T) && doubleInput!(U)) {
    return goodnessFit!(pearsonChiSqElem, T, U)(observed, expected, countProp);
}

unittest {
    // Test to see whether a set of categorical observations differs
    // statistically from a discrete uniform distribution.
    uint[] observed = [980, 1028, 1001, 964, 1102];
    auto expected = repeat(cast(double) sum(observed) / observed.length);
    auto res = chiSquareFit(observed, expected, Expected.count);
    assert(approxEqual(res, 0.0207));
    assert(approxEqual(res.testStat, 11.59));

    auto expected2 = [5.0, 5, 5, 5, 5, 0];
    observed ~= 0;
    auto res2 = chiSquareFit(observed, expected2);
    assert(approxEqual(res2, 0.0207));
    assert(approxEqual(res2.testStat, 11.59));
}

// Alias for old name, for backwards compatibility.  Don't document it
// because it will be deprecated eventually.
alias chiSquareFit chiSqrFit;

/**The G or likelihood ratio chi-square test for goodness of fit.  Roughly
 * the same as Pearson's chi-square test (chiSquareFit), but may be more
 * accurate in certain situations and less accurate in others.  However, it is
 * still based on asymptotic distributions, and is not exact. Usage is is
 * identical to chiSquareFit.
 *
 * References:  http://en.wikipedia.org/wiki/G_test
 *
 */
TestRes gTestFit(T, U)(T observed, U expected, Expected countProp = Expected.proportion)
if(doubleInput!(T) && doubleInput!(U)) {
    return goodnessFit!(gTestElem, T, U)(observed, expected, countProp);
}
// No unittest because I can't find anything to test this against.  However,
// it's hard to imagine how it could be wrong, given that goodnessFit() and
// gTestElem() both work, and, as expected, this function produces roughly
// the same results as chiSquareFit.

private TestRes goodnessFit(alias elemFun, T, U)(T observed, U expected, Expected countProp)
if(doubleInput!(T) && doubleInput!(U)) {
    if(countProp == Expected.proportion) {
        dstatsEnforce(isForwardRange!(U),
            "Can't use expected proportions instead of counts with input ranges.");
    }

    uint len = 0;
    double chiSq = 0;
    double multiplier = 1;

    // Normalize proportions to add up to the sum of the data.
    if(countProp == Expected.proportion) {
        double expectSum = 0;
        multiplier = 0;
        auto obsCopy = observed.save;
        auto expCopy = expected.save;
        while(!obsCopy.empty && !expCopy.empty) {
            multiplier += obsCopy.front;
            expectSum += expCopy.front;
            obsCopy.popFront;
            expCopy.popFront;
        }
        multiplier /= expectSum;
    }

    while(!observed.empty && !expected.empty) {
        scope(exit) {
            observed.popFront();
            expected.popFront();
        }
        double e = expected.front * multiplier;

        // If e is zero, then we should just treat the cell as if it didn't
        // exist.
        if(e == 0) {
            dstatsEnforce(observed.front == 0,
                "Can't have non-zero observed value w/ zero expected value.");
            continue;
        }

        chiSq += elemFun(observed.front, e);
        len++;
    }

    if(isNaN(chiSq)) {
        return TestRes(double.nan, double.nan);
    }

    return TestRes(chiSq, chiSquareCDFR(chiSq, len - 1));
}

/**
The exact multinomial goodness of fit test for whether a set of counts
fits a hypothetical distribution.  counts is an input range of counts.
proportions is an input range of expected proportions.  These are normalized
automatically, so they can sum to any value.

Returns:  The P-value for the null that counts is a sample from proportions
against the alternative that it isn't.

Notes:  This test is EXTREMELY slow for anything but very small samples and
degrees of freedom.  The Pearson's chi-square (chiSquareFit()) or likelihood
ratio chi-square (gTestFit()) are good enough approximations unless sample
size is very small.
 */
double multinomialTest(U, F)(U countsIn, F proportions)
if(isInputRange!U && isInputRange!F &&
   isIntegral!(ElementType!U) && isFloatingPoint!(ElementType!(F))) {
    auto alloc = newRegionAllocator();

    static if(isRandomAccessRange!U && hasLength!U) {
        alias countsIn counts;
    } else {
        auto counts = alloc.array(countsIn);
    }

    uint N = sum(counts);

    double[] logPs;
    static if(std.range.hasLength!F) {
        logPs = alloc.uninitializedArray!(double[])(proportions.length);
        size_t pIndex;
        foreach(p; proportions) {
            logPs[pIndex++] = p;
        }
    } else {
        auto app = appender(logPs);
        foreach(p; proportions) {
            app.put(p);
        }
        logPs = app.data;
    }

    logPs[] /= reduce!"a + b"(0.0, logPs);
    foreach(ref elem; logPs) {
        elem = log(elem);
    }


    double[] logs = alloc.uninitializedArray!(double[])(N + 1);
    logs[0] = 0;
    foreach(i; 1..logs.length) {
        logs[i] = log(i);
    }

    double nFact = logFactorial(N);
    double pVal = 0;
    uint nLeft = N;
    double pSoFar = nFact;

    double pActual = nFact;
    foreach(i, count; counts) {
        pActual += logPs[i] * count - logFactorial(count);
    }
    pActual -= pActual * 1e-6;  // Epsilon to handle numerical inaccuracy.

    void doIt(uint pos) {
        if(pos == counts.length - 1) {
            immutable pOld = pSoFar;
            pSoFar += logPs[$ - 1] * nLeft - logFactorial(nLeft);

            if(pSoFar <= pActual) {
                pVal += exp(pSoFar);
            }
            pSoFar = pOld;
            return;
        }

        uint nLeftOld = nLeft;
        immutable pOld = pSoFar;
        double pAdd = 0;

        foreach(i; 0..nLeft + 1) {
            if(i > 0) {
                pAdd += logPs[pos] - logs[i];
            }
            pSoFar = pOld + pAdd;
            doIt(pos + 1);
            nLeft--;
        }
        nLeft = nLeftOld;
        pSoFar = pOld;
    }
    doIt(0);
    return pVal;
}

unittest {
    // Nothing to test this against for more than 1 df, but it matches
    // chi-square roughly and should take the same paths for 2 vs. N degrees
    // of freedom.
    for(uint n = 4; n <= 100; n += 4) {
        foreach(k; 0..n + 1) {
            for(double p = 0.05; p <= 0.95; p += 0.05) {
                double bino = binomialTest(k, n, p);
                double[] ps = [p, 1 - p];
                uint[] counts = [k, n - k];
                double multino = multinomialTest(counts, ps);
                //writeln(k, "\t", n, "\t", p, "\t", bino, "\t", multino);
                assert(approxEqual(bino, multino),
                    text(bino, '\t', multino, '\t', k, '\t', n, '\t', p));
            }
        }
    }
}

/**
Performs a Pearson's chi-square test on a contingency table of arbitrary
dimensions.  When the chi-square test is mentioned, this is usually the one
being referred to.  Takes a set of finite forward ranges, one for each column
in the contingency table.  These can be expressed either as a tuple of ranges
or a range of ranges.  Returns a P-value for the alternative hypothesis that
frequencies in each row of the contingency table depend on the column against
the null that they don't.

Notes:  The chi-square test relies on asymptotic statistical properties
and is therefore not exact.  The typical rule of thumb is that each cell
should have an expected value of at least 5.  However, this is likely to
be unnecessarily stringent.

Yates's continuity correction is never used in this implementation.  If
you want something that's guaranteed to be conservative, use fisherExact().

This is, for all practical purposes, an inherently non-directional test.
Therefore, the one-sided verses two-sided option is not provided.

For 2x2 contingency tables, fisherExact is a more conservative test, in that
the type I error rate is guaranteed to never be above the nominal P-value.
However, even for small sample sizes this test may produce results closer
to the true P-value, at the risk of possibly being non-conservative.

Examples:
---
// Test to see whether the relative frequency of outcome 0, 1, and 2
// depends on the treatment (drug1, drug2 or placebo) in some hypothetical
// experiment.  For example, 1500 people had outcome 2 if they were treated
// with drug1 and 1100 had outcome 1 if treated with placebo.
uint[] drug1 = [1000, 2000, 1500];
uint[] drug2 = [1500, 3000, 2300];
uint[] placebo = [500, 1100, 750];
auto result1 = chiSquareContingency(drug1, drug2, placebo);

// The following uses a range of ranges instead of an array of ranges,
// and will produce equivalent results.
auto rangeOfRanges = [drug1, drug2, placebo];
auto result2 = chiSquareContingency(rangeOfRanges);
---

References: http://en.wikipedia.org/wiki/Pearson%27s_chi-square_test
 *
 */
TestRes chiSquareContingency(T...)(T inputData) {
    return testContingency!(pearsonChiSqElem, T)(inputData);
}

unittest {
    // Test array version.  Using VassarStat's chi-square calculator.
    uint[][] table1 = [[60, 80, 70],
                       [20, 50, 40],
                       [10, 15, 11]];
    uint[][] table2 = [[60, 20, 10],
                       [80, 50, 15],
                       [70, 40, 11]];
    assert(approxEqual(chiSquareContingency(table1), 0.3449));
    assert(approxEqual(chiSquareContingency(table2), 0.3449));
    assert(approxEqual(chiSquareContingency(table1).testStat, 4.48));

    // Test tuple version.
    auto p1 = chiSquareContingency(cast(uint[]) [31, 41, 59],
                                cast(uint[]) [26, 53, 58],
                                cast(uint[]) [97, 93, 93]);
    assert(approxEqual(p1, 0.0059));

    auto p2 = chiSquareContingency(cast(uint[]) [31, 26, 97],
                                cast(uint[]) [41, 53, 93],
                                cast(uint[]) [59, 58, 93]);
    assert(approxEqual(p2, 0.0059));

    uint[] drug1 = [1000, 2000, 1500];
    uint[] drug2 = [1500, 3000, 2300];
    uint[] placebo = [500, 1100, 750];
    assert(approxEqual(chiSquareContingency(drug1, drug2, placebo), 0.2397));
}

// Alias for old name, for backwards compatibility.  Don't document it
// because it is deprecated and has been scheduled for deprecation for
// ages.
deprecated alias chiSquareContingency chiSqrContingency;

/**
This struct is a subtype of TestRes and is used to return the results of
gTestContingency and gTestObs.  Due to the information theoretic interpretation
of the G test, it contains an extra field to return the mutual information
in bits.
*/
struct GTestRes {
    ///
    TestRes testRes;

    ///
    alias testRes this;

    /**
    The mutual info of the two random variables in the joint distribution
    represented by the contingency table, in bits (base 2).
    */
    double mutualInfo;
}

/**
The G or likelihood ratio chi-square test for contingency tables.  Roughly
the same as Pearson's chi-square test (chiSquareContingency), but may be more
accurate in certain situations and less accurate in others.

Like Pearson's Chi-square, the G-test is based on asymptotic distributions,
and is not exact. Usage is is identical to chiSquareContingency.

Note:  This test can be thought of as a test for nonzero mutual information
between the random variables represented by the rows and the columns,
since the test statistic and P-value are strictly increasing
and strictly decreasing, respectively, in mutual information.  Therefore, this
function returns a GTestRes, which is a subtype of TestRes and also gives
the mutual information for use in information theoretic settings.

References:  http://en.wikipedia.org/wiki/G_test, last retrieved 1/22/2011
*/
GTestRes gTestContingency(T...)(T inputData) {
    return testContingency!(gTestElem, T)(inputData);
}

unittest {
    // Values from example at http://udel.edu/~mcdonald/statgtestind.html
    // Handbook of Biological Statistics.
    uint[] withoutCHD = [268, 199, 42];
    uint[] withCHD = [807, 759, 184];
    auto res = gTestContingency(withoutCHD, withCHD);
    assert(approxEqual(res.testStat, 7.3));
    assert(approxEqual(res.p, 0.026));
    assert(approxEqual(res.mutualInfo, 0.0023313));


    uint[] moringa = [127, 99, 264];
    uint[] vicinus = [116, 67, 161];
    auto res2 = gTestContingency(moringa, vicinus);
    assert(approxEqual(res2.testStat, 6.23));
    assert(approxEqual(res2.p, 0.044));
    assert(approxEqual(res2.mutualInfo, 0.00538613));
}

// Pearson and likelihood ratio code are pretty much the same.  Factor out
// the one difference into a function that's a template parameter.  However,
// for API simplicity, this is hidden and they look like two separate functions.
private GTestRes testContingency(alias elemFun, T...)(T rangesIn) {
    auto alloc = newRegionAllocator();
    static if(isInputRange!(T[0]) && T.length == 1 &&
        isForwardRange!(typeof(rangesIn[0].front()))) {
        auto ranges = alloc.array(rangesIn[0]);

        foreach(ref range; ranges) {
            range = range.save;
        }

    } else static if(allSatisfy!(isForwardRange, typeof(rangesIn))) {
        auto saved = saveAll(rangesIn);
        auto ranges = saved.expand;
    } else {
        static assert(0, "Can only perform contingency table test" ~
            " on a tuple of ranges or a range of ranges.");
    }

    double[] colSums = alloc.uninitializedArray!(double[])(ranges.length);
    colSums[] = 0;
    size_t nCols = 0;
    immutable size_t nRows = ranges.length;
    foreach(ri, range; ranges) {
        size_t curLen = 0;
        foreach(elem; range.save) {
            colSums[ri] += cast(double) elem;
            curLen++;
        }
        if(ri == 0) {
            nCols = curLen;
        } else {
            assert(curLen == nCols);
        }
    }

    bool noneEmpty() {
        foreach(range; ranges) {
            if(range.empty) {
                return false;
            }
        }
        return true;
    }

    void popAll() {
        foreach(i, range; ranges) {
            ranges[i].popFront;
        }
    }

    double sumRow() {
        double rowSum = 0;
        foreach(range; ranges) {
            rowSum += cast(double) range.front;
        }
        return rowSum;
    }

    double chiSq = 0;
    immutable double NNeg1 = 1.0 / sum(colSums);
    while(noneEmpty) {
        auto rowSum = sumRow();
        foreach(ri, range; ranges) {
            double expected = NNeg1 * rowSum * colSums[ri];
            chiSq += elemFun(range.front, expected);
        }
        popAll();
    }

    if(isNaN(chiSq)) {
        return GTestRes(TestRes(double.nan, double.nan), double.nan);
    }

    // This can happen in some cases due to numerical fuzz.
    if(chiSq > 1e-5 && chiSq <= 0) {
        return GTestRes(TestRes(0, 1), 0);
    }

    immutable pVal = (chiSq >= 0) ?
        chiSquareCDFR(chiSq, (nRows - 1) * (nCols - 1)) : double.nan;

    // 1 / (2 * LN2), for converting chiSq to mutualInfo.
    enum chiToMi = 1 / (2 * LN2);

    // This is the mutual information between the two random variables
    // represented by the contingency table, only if we're doing a G test.
    // If we're doing a Pearson's test, it's a completely meaningless quantity,
    // but never gets returned by any public function.
    immutable mutualInfo = chiSq * NNeg1 * chiToMi;

    return GTestRes(TestRes(chiSq, pVal), mutualInfo);
}

private double pearsonChiSqElem(double observed, double expected) pure nothrow {
    immutable diff = observed - expected;
    return diff * diff / expected;
}

private double gTestElem(double observed, double expected) pure nothrow {
    return (observed == 0 && expected > 0) ? 0 :
        (observed * log(observed / expected) * 2);
}

/**
Given two vectors of observations of jointly distributed variables x, y, tests
the null hypothesis that values in x are independent of the corresponding
values in y.  This is done using Pearson's Chi-Square Test.  For a similar test
that assumes the data has already been tabulated into a contingency table, see
chiSquareContingency.

x and y must both be input ranges.  If they are not the same length, a
DstatsArgumentException is thrown.

Examples:
---
// Test whether the appearance of "foo" vs. "bar" is independent of the
// appearance of "baz" vs. "xxx".
auto x = ["foo", "bar", "bar", "foo", "foo"];
auto y = ["xxx", "baz", "baz", "xxx", "baz"];
auto result = chiSquareObs(x, y);

// This is equivalent to:
auto contingency = new uint[][](2, 2);
foreach(i; 0..x.length) {
    immutable index1 = (x[i] == "foo");
    immutable index2 = (y[i] == "xxx");
    contingency[index1][index2]++;
}

auto result2 = chiSquareContingency(contingency);
---
*/
TestRes chiSquareObs(T, U)(T x, U y)
if(isInputRange!T && isInputRange!U) {
    uint xFreedom, yFreedom, n;
    typeof(return) ret;

    static if(!hasLength!T && !hasLength!U) {
        ret.testStat = toContingencyScore!(T, U, uint)
            (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
    } else {
        immutable minLen = min(x.length, y.length);
        if(minLen <= ubyte.max) {
            ret.testStat = toContingencyScore!(T, U, ubyte)
                (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
        } else if(minLen <= ushort.max) {
            ret.testStat = toContingencyScore!(T, U, ushort)
                (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
        } else {
            ret.testStat = toContingencyScore!(T, U, uint)
                (x, y, &pearsonChiSqElem, xFreedom, yFreedom, n);
        }
    }

    ret.p = chiSquareCDFR(ret.testStat, xFreedom * yFreedom);
    return ret;
}

unittest {
    // We know the chi-square contingency works, so test that the automatic
    // binning works, too.
    ubyte[] obs1 = [1, 2, 3, 1, 2, 3, 1, 2, 3];
    ubyte[] obs2 = [1, 3, 2, 1, 3, 2, 1, 3, 2];

    uint[][] cTable = [[3, 0, 0],
                       [0, 0, 3],
                       [0, 3, 0]];
    auto gRes = chiSquareContingency(cTable);
    auto miRes = chiSquareObs(obs1, obs2);
    foreach(ti, elem; miRes.tupleof) {
        assert(approxEqual(elem, gRes.tupleof[ti]));
    }

    auto x = ["foo", "bar", "bar", "foo", "foo"];
    auto y = ["xxx", "baz", "baz", "xxx", "baz"];
    auto result = chiSquareObs(x, y);
    assert(approxEqual(result.testStat, 2.22222222));
    assert(approxEqual(result.p, 0.136037));

    auto contingency = new uint[][](2, 2);
    foreach(i; 0..x.length) {
        immutable index1 = (x[i] == "foo");
        immutable index2 = (y[i] == "xxx");
        contingency[index1][index2]++;
    }

    auto result2 = chiSquareContingency(contingency);
    assert(approxEqual(result.testStat, result2.testStat),
        text(result.testStat, ' ', result2.testStat));
    assert(approxEqual(result.p, result2.p));
}

/**
Given two ranges of observations of jointly distributed variables x, y, tests
the null hypothesis that values in x are independent of the corresponding
values in y.  This is done using the Likelihood Ratio G test.  Usage is similar
to chiSquareObs.  For an otherwise identical test that assumes the data has
already been tabulated into a contingency table, see gTestContingency.

Note:  This test can be thought of as a test for nonzero mutual information
between x and y, since the test statistic and P-value are strictly increasing
and strictly decreasing, respectively, in mutual information.  Therefore, this
function returns a GTestRes, which is a subtype of TestRes and also gives
the mutual information for use in information theoretic settings.
*/
GTestRes gTestObs(T, U)(T x, U y)
if(isInputRange!T && isInputRange!U) {
    uint xFreedom, yFreedom, n;
    typeof(return) ret;

    static if(!hasLength!T && !hasLength!U) {
        ret.testStat = toContingencyScore!(T, U, uint)
            (x, y, &gTestElem, xFreedom, yFreedom, n);
    } else {
        immutable minLen = min(x.length, y.length);
        if(minLen <= ubyte.max) {
            ret.testStat = toContingencyScore!(T, U, ubyte)
                (x, y, &gTestElem, xFreedom, yFreedom, n);
        } else if(minLen <= ushort.max) {
            ret.testStat = toContingencyScore!(T, U, ushort)
                (x, y, &gTestElem, xFreedom, yFreedom, n);
        } else {
            ret.testStat = toContingencyScore!(T, U, uint)
                (x, y, &gTestElem, xFreedom, yFreedom, n);
        }
    }

    ret.p = chiSquareCDFR(ret.testStat, xFreedom * yFreedom);
    ret.mutualInfo = ret.testStat / (2 * LN2 * n);
    return ret;
}

unittest {
    // We know the g test contingency works, so test that the automatic binning
    // works, too.
    ubyte[] obs1 = [1, 2, 3, 1, 2, 3, 1, 2, 3];
    ubyte[] obs2 = [1, 3, 2, 1, 3, 2, 1, 3, 2];

    uint[][] cTable = [[3, 0, 0],
                       [0, 0, 3],
                       [0, 3, 0]];
    auto gRes = gTestContingency(cTable);
    auto miRes = gTestObs(obs1, obs2);
    foreach(ti, elem; miRes.tupleof) {
        assert(approxEqual(elem, gRes.tupleof[ti]));
    }

    auto x = ["foo", "bar", "bar", "foo", "foo"];
    auto y = ["xxx", "baz", "baz", "xxx", "baz"];
    auto result = gTestObs(x, y);
    assert(approxEqual(result.testStat, 2.91103));
    assert(approxEqual(result.p, 0.0879755));
    assert(approxEqual(result.mutualInfo, 0.419973));
}

package double toContingencyScore(T, U, Uint)
(T x, U y, double function(double, double) elemFun,
 out uint xFreedom, out uint yFreedom, out uint nPtr) {

    import dstats.infotheory : NeedsHeap;
    enum needsHeap = NeedsHeap!T ||
        NeedsHeap!U;
    alias dstats.infotheory.ObsEnt!(ElementType!T, ElementType!U) ObsType;

    static if(needsHeap) {
        Uint[ObsType] jointCounts;
        Uint[ElementType!T] xCounts;
        Uint[ElementType!U] yCounts;
    } else {
        auto alloc = newRegionAllocator();
        dstatsEnforce(x.length == y.length,
            "Can't calculate mutual info with different length vectors.");
        immutable len = x.length;
        auto jointCounts = StackHash!(ObsType, Uint)(max(20, len / 20), alloc);
        auto xCounts = StackHash!(ElementType!T, Uint)(max(10, len / 40), alloc);
        auto yCounts = StackHash!(ElementType!U, Uint)(max(10, len / 40), alloc);
    }

    uint n = 0;
    while(!x.empty && !y.empty) {
        n++;
        auto a = x.front;
        auto b = y.front;
        jointCounts[ObsType(a, b)]++;
        xCounts[a]++;
        yCounts[b]++;

        x.popFront();
        y.popFront();
    }

    dstatsEnforce(x.empty && y.empty,
        "Can't calculate mutual info with different length vectors.");

    xFreedom = cast(uint) xCounts.length - 1;
    yFreedom = cast(uint) yCounts.length - 1;
    nPtr = n;

    double ret = 0;
    immutable double nNeg1 = 1.0 / n;
    foreach(key1, marg1; xCounts) foreach(key2, marg2; yCounts) {
        immutable observed = jointCounts.get(
            ObsType(key1, key2), 0
        );
        immutable expected = marg1 * nNeg1 * marg2;
        ret += elemFun(observed, expected);
    }

    return ret;
}

/**
Fisher's Exact test for difference in odds between rows/columns
in a 2x2 contingency table.  Specifically, this function tests the odds
ratio, which is defined, for a contingency table c, as (c[0][0] * c[1][1])
 / (c[1][0] * c[0][1]).  Alternatives are Alt.less, meaning true odds ratio
< 1, Alt.greater, meaning true odds ratio > 1, and Alt.twoSided, meaning
true odds ratio != 1.

Accepts a 2x2 contingency table as an array of arrays of uints.
For now, only does 2x2 contingency tables.

Notes:  Although this test is "exact" in that it does not rely on asymptotic
approximations, it is very statistically conservative when the marginals
are not truly fixed in the experimental design in question.  If a
closer but possibly non-conservative approximation of the true P-value is
desired, Pearson's chi-square test (chiSquareContingency) may perform better,
even for small samples.

Returns:  A TestRes of the odds ratio and the P-value against the given
alternative.

Examples:
---
double res = fisherExact([[2u, 7], [8, 2]], Alt.less);
assert(approxEqual(res.p, 0.01852));  // Odds ratio is very small in this case.
assert(approxEqual(res.testStat, 4.0 / 56.0));
---

References:  http://en.wikipedia.org/wiki/Fisher%27s_Exact_Test
*/
TestRes fisherExact(T)(const T[2][2] contingencyTable, Alt alt = Alt.twoSided)
if(isIntegral!(T)) {
    foreach(range; contingencyTable) {
        foreach(elem; range) {
            dstatsEnforce(elem >= 0,
                "Cannot have negative elements in a contingency table.");
        }
    }

    static double fisherLower(const T[2][2] contingencyTable) {
        alias contingencyTable c;
        return hypergeometricCDF(c[0][0], c[0][0] + c[0][1], c[1][0] + c[1][1],
                                 c[0][0] + c[1][0]);
    }

    static double fisherUpper(const T[2][2] contingencyTable) {
        alias contingencyTable c;
        return hypergeometricCDFR(c[0][0], c[0][0] + c[0][1], c[1][0] + c[1][1],
                                 c[0][0] + c[1][0]);
    }


    alias contingencyTable c;  // Save typing.
    immutable oddsRatio = cast(double) c[0][0] * c[1][1] / c[0][1] / c[1][0];
    if(alt == Alt.none) {
        return TestRes(oddsRatio);
    } else if(alt == Alt.less) {
        return TestRes(oddsRatio, fisherLower(contingencyTable));
    } else if(alt == Alt.greater) {
        return TestRes(oddsRatio, fisherUpper(contingencyTable));
    }


    immutable uint n1 = c[0][0] + c[0][1],
                   n2 = c[1][0] + c[1][1],
                   n  = c[0][0] + c[1][0];

    immutable uint mode =
        cast(uint) ((cast(double) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
    immutable double pExact = hypergeometricPMF(c[0][0], n1, n2, n);
    immutable double pMode = hypergeometricPMF(mode, n1, n2, n);

    enum epsilon = 1 - 1e-5;
    if(approxEqual(pExact, pMode, 1 - epsilon)) {
        return TestRes(oddsRatio, 1);
    } else if(c[0][0] < mode) {
        immutable double pLower = hypergeometricCDF(c[0][0], n1, n2, n);

        if(hypergeometricPMF(n, n1, n2, n) > pExact / epsilon) {
            return TestRes(oddsRatio, pLower);
        }

        // Binary search for where to begin upper half.
        uint min = mode, max = n, guess = uint.max;
        while(max - min > 1) {
            guess = cast(uint) (
                    (max == min + 1 && guess == min) ? max :
                    (cast(ulong) max + cast(ulong) min) / 2UL);

            immutable double pGuess = hypergeometricPMF(guess, n1, n2, n);
            if(pGuess <= pExact &&
                hypergeometricPMF(guess - 1, n1, n2, n) > pExact) {
                break;
            } else if(pGuess < pExact) {
                max = guess;
            } else min = guess;
        }

        if(guess == uint.max) {
            guess = min;
        }

        while(guess > 0 && hypergeometricPMF(guess, n1, n2, n) < pExact * epsilon) {
            guess--;
        }

        while(hypergeometricPMF(guess, n1, n2, n) > pExact / epsilon) {
            guess++;
        }

        static import std.algorithm;
        double p = std.algorithm.min(pLower +
               hypergeometricCDFR(guess, n1, n2, n), 1.0);
        return TestRes(oddsRatio, p);
    } else {
        immutable double pUpper = hypergeometricCDFR(c[0][0], n1, n2, n);

        if(hypergeometricPMF(0, n1, n2, n) > pExact / epsilon) {
            return TestRes(oddsRatio, pUpper);
        }

        // Binary search for where to begin lower half.
        uint min = 0, max = mode, guess = uint.max;
        while(max - min > 1) {
            guess = cast(uint) (
                    (max == min + 1 && guess == min) ? max :
                    (cast(ulong) max + cast(ulong) min) / 2UL);
            immutable double pGuess = hypergeometricPMF(guess, n1, n2, n);

            if(pGuess <= pExact &&
                hypergeometricPMF(guess + 1, n1, n2, n) > pExact) {
                break;
            } else if(pGuess <= pExact) {
                min = guess;
            } else max = guess;
        }

        if(guess == uint.max) {
            guess = min;
        }

        while(hypergeometricPMF(guess, n1, n2, n) < pExact * epsilon) {
            guess++;
        }

        while(guess > 0 &&
            hypergeometricPMF(guess, n1, n2, n) > pExact / epsilon) {
            guess--;
        }

        static import std.algorithm;
        double p = std.algorithm.min(pUpper +
               hypergeometricCDF(guess, n1, n2, n), 1.0);
        return TestRes(oddsRatio, p);
    }
}

/**
Convenience function.  Converts a dynamic array to a static one, then
calls the overload.
*/
TestRes fisherExact(T)(const T[][] contingencyTable, Alt alt = Alt.twoSided)
if(isIntegral!(T)) {
    dstatsEnforce(contingencyTable.length == 2 &&
            contingencyTable[0].length == 2 &&
            contingencyTable[1].length == 2,
            "Fisher exact only supports 2x2 tables.");

    T[2][2] newTable;
    newTable[0][0] = contingencyTable[0][0];
    newTable[0][1] = contingencyTable[0][1];
    newTable[1][1] = contingencyTable[1][1];
    newTable[1][0] = contingencyTable[1][0];
    return fisherExact(newTable, alt);
}

unittest {
    // Simple, naive impl. of two-sided to test against.
    static double naive(const uint[][] c) {
        immutable uint n1 = c[0][0] + c[0][1],
                   n2 = c[1][0] + c[1][1],
                   n  = c[0][0] + c[1][0];
        immutable uint mode =
            cast(uint) ((cast(double) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
        immutable double pExact = hypergeometricPMF(c[0][0], n1, n2, n);
        immutable double pMode = hypergeometricPMF(mode, n1, n2, n);
        if(approxEqual(pExact, pMode, 1e-7))
            return 1;
        double sum = 0;
        foreach(i; 0..n + 1) {
            double pCur = hypergeometricPMF(i, n1, n2, n);
            if(pCur <= pExact / (1 - 1e-5))
                sum += pCur;
        }
        return sum;
    }

    uint[][] c = new uint[][](2, 2);

    foreach(i; 0..100_000) {
        c[0][0] = uniform(0U, 51U);
        c[0][1] = uniform(0U, 51U);
        c[1][0] = uniform(0U, 51U);
        c[1][1] = uniform(0U, 51U);
        double naiveAns = naive(c);
        double fastAns = fisherExact(c);
        assert(approxEqual(naiveAns, fastAns), text(c, naiveAns, fastAns));
    }

    auto res = fisherExact([[19000, 80000], [20000, 90000]]);
    assert(approxEqual(res.testStat, 1.068731));
    assert(approxEqual(res, 3.319e-9));
    res = fisherExact([[18000, 80000], [20000, 90000]]);
    assert(approxEqual(res, 0.2751));
    res = fisherExact([[14500, 20000], [30000, 40000]]);
    assert(approxEqual(res, 0.01106));
    res = fisherExact([[100, 2], [1000, 5]]);
    assert(approxEqual(res, 0.1301));
    res = fisherExact([[2, 7], [8, 2]]);
    assert(approxEqual(res, 0.0230141));
    res = fisherExact([[5, 1], [10, 10]]);
    assert(approxEqual(res, 0.1973244));
    res = fisherExact([[5, 15], [20, 20]]);
    assert(approxEqual(res, 0.0958044));
    res = fisherExact([[5, 16], [20, 25]]);
    assert(approxEqual(res, 0.1725862));
    res = fisherExact([[10, 5], [10, 1]]);
    assert(approxEqual(res, 0.1973244));
    res = fisherExact([[5, 0], [1, 4]]);
    assert(approxEqual(res.p, 0.04761904));
    res = fisherExact([[0, 1], [3, 2]]);
    assert(approxEqual(res.p, 1.0));
    res = fisherExact([[0, 2], [6, 4]]);
    assert(approxEqual(res.p, 0.4545454545));
    res = fisherExact([[2, 7], [8, 2]], Alt.less);
    assert(approxEqual(res, 0.01852));
    res = fisherExact([[5, 1], [10, 10]], Alt.less);
    assert(approxEqual(res, 0.9783));
    res = fisherExact([[5, 15], [20, 20]], Alt.less);
    assert(approxEqual(res, 0.05626));
    res = fisherExact([[5, 16], [20, 25]], Alt.less);
    assert(approxEqual(res, 0.08914));
    res = fisherExact([[2, 7], [8, 2]], Alt.greater);
    assert(approxEqual(res, 0.999));
    res = fisherExact([[5, 1], [10, 10]], Alt.greater);
    assert(approxEqual(res, 0.1652));
    res = fisherExact([[5, 15], [20, 20]], Alt.greater);
    assert(approxEqual(res, 0.985));
    res = fisherExact([[5, 16], [20, 25]], Alt.greater);
    assert(approxEqual(res, 0.9723));
}

/**
Performs a Kolmogorov-Smirnov (K-S) 2-sample test.  The K-S test is a
non-parametric test for a difference between two empirical distributions or
between an empirical distribution and a reference distribution.

Returns:  A TestRes with the K-S D value and a P value for the null that
FPrime is distributed identically to F against the alternative that it isn't.
This implementation uses a signed D value to indicate the direction of the
difference between distributions.  To get the D value used in standard
notation, simply take the absolute value of this D value.

Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.

References:  http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
 */
TestRes ksTest(T, U)(T F, U Fprime)
if(doubleInput!(T) && doubleInput!(U)) {
    double D = ksTestD(F, Fprime);
    return TestRes(D, ksPval(F.length, Fprime.length, D));
}

unittest {
    assert(approxEqual(ksTest([1,2,3,4,5], [1,2,3,4,5]).testStat, 0));
    assert(approxEqual(ksTestDestructive([1,2,3,4,5], [1,2,2,3,5]).testStat, -.2));
    assert(approxEqual(ksTest([-1,0,2,8, 6], [1,2,2,3,5]).testStat, .4));
    assert(approxEqual(ksTest([1,2,3,4,5], [1,2,2,3,5,7,8]).testStat, .2857));
    assert(approxEqual(ksTestDestructive([1, 2, 3, 4, 4, 4, 5],
           [1, 2, 3, 4, 5, 5, 5]).testStat, .2857));

    assert(approxEqual(ksTest([1, 2, 3, 4, 4, 4, 5], [1, 2, 3, 4, 5, 5, 5]),
           .9375));
    assert(approxEqual(ksTestDestructive([1, 2, 3, 4, 4, 4, 5],
        [1, 2, 3, 4, 5, 5, 5]), .9375));
}

template isArrayLike(T) {
    enum bool isArrayLike = hasSwappableElements!(T) && hasAssignableElements!(T)
        && hasLength!(T) && isRandomAccessRange!(T);
}

/**
One-sample Kolmogorov-Smirnov test against a reference distribution.
Takes a callable object for the CDF of refernce distribution.

Returns:  A TestRes with the Kolmogorov-Smirnov D value and a P value for the
null that Femp is a sample from F against the alternative that it isn't. This
implementation uses a signed D value to indicate the direction of the
difference between distributions.  To get the D value used in standard
notation, simply take the absolute value of this D value.

Bugs:  Exact calculation not implemented.  Uses asymptotic approximation.

Examples:
---
auto stdNormal = parametrize!(normalCDF)(0.0, 1.0);
auto empirical = [1, 2, 3, 4, 5];
auto res = ksTest(empirical, stdNormal);
---

References:  http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
 */
TestRes ksTest(T, Func)(T Femp, Func F)
if(doubleInput!(T) && is(ReturnType!(Func) : double)) {
    double D = ksTestD(Femp, F);
    return TestRes(D, ksPval(Femp.length, D));
}

unittest {
    auto stdNormal = paramFunctor!(normalCDF)(0.0, 1.0);
    assert(approxEqual(ksTest([1,2,3,4,5].dup, stdNormal).testStat, -.8413));
    assert(approxEqual(ksTestDestructive([-1,0,2,8, 6].dup, stdNormal).testStat, -.5772));
    auto lotsOfTies = [5,1,2,2,2,2,2,2,3,4].dup;
    assert(approxEqual(ksTest(lotsOfTies, stdNormal).testStat, -0.8772));

    assert(approxEqual(ksTest([0,1,2,3,4].dup, stdNormal), .03271));

    auto uniform01 = parametrize!(uniformCDF)(0, 1);
    assert(approxEqual(ksTestDestructive([0.1, 0.3, 0.5, 0.9, 1].dup, uniform01), 0.7591));

}

/**Same as ksTest, except sorts in place, avoiding memory allocations.*/
TestRes ksTestDestructive(T, U)(T F, U Fprime)
if(isArrayLike!(T) && isArrayLike!(U)) {
    double D = ksTestDDestructive(F, Fprime);
    return TestRes(D, ksPval(F.length, Fprime.length, D));
}

///Ditto.
TestRes ksTestDestructive(T, Func)(T Femp, Func F)
if(isArrayLike!(T) && is(ReturnType!Func : double)) {
    double D =  ksTestDDestructive(Femp, F);
    return TestRes(D, ksPval(Femp.length, D));
}

private double ksTestD(T, U)(T F, U Fprime)
if(isInputRange!(T) && isInputRange!(U)) {
    auto alloc = newRegionAllocator();
    return ksTestDDestructive(alloc.array(F), alloc.array(Fprime));
}

private double ksTestDDestructive(T, U)(T F, U Fprime)
if(isArrayLike!(T) && isArrayLike!(U)) {
    qsort(F);
    qsort(Fprime);
    double D = 0;
    size_t FprimePos = 0;
    foreach(i; 0..2) {  //Test both w/ Fprime x vals, F x vals.
        double diffMult = (i == 0) ? 1 : -1;
        foreach(FPos, Xi; F) {
            if(FPos < F.length - 1 && F[FPos + 1] == Xi)
                continue;  //Handle ties.
            while(FprimePos < Fprime.length && Fprime[FprimePos] <= Xi) {
                FprimePos++;
            }
            double diff = diffMult * (cast(double) (FPos + 1) / F.length -
                       cast(double) FprimePos / Fprime.length);
            if(abs(diff) > abs(D))
                D = diff;
        }
        swap(F, Fprime);
        FprimePos = 0;
    }
    return D;
}

private double ksTestD(T, Func)(T Femp, Func F)
if(doubleInput!(T) && is(ReturnType!Func : double)) {
    auto alloc = newRegionAllocator();
    return ksTestDDestructive(alloc.array(Femp), F);
}

private double ksTestDDestructive(T, Func)(T Femp, Func F)
if(isArrayLike!(T) && is(ReturnType!Func : double)) {
    qsort(Femp);
    double D = 0;

    foreach(FPos, Xi; Femp) {
        double diff = cast(double) FPos / Femp.length - F(Xi);
        if(abs(diff) > abs(D))
            D = diff;
    }

    return D;
}

private double ksPval(ulong N, ulong Nprime, double D)
in {
    assert(D >= -1);
    assert(D <= 1);
} body {
    return 1 - kolmogorovDistrib(sqrt(cast(double) (N * Nprime) / (N + Nprime)) * abs(D));
}

private double ksPval(ulong N, double D)
in {
    assert(D >= -1);
    assert(D <= 1);
} body {
    return 1 - kolmogorovDistrib(abs(D) * sqrt(cast(double) N));
}

/**
Wald-wolfowitz or runs test for randomness of the distribution of
elements for which positive() evaluates to true.  For example, given
a sequence of coin flips [H,H,H,H,H,T,T,T,T,T] and a positive() function of
"a == 'H'", this test would determine that the heads are non-randomly
distributed, since they are all at the beginning of obs.  This is done
by counting the number of runs of consecutive elements for which
positive() evaluates to true, and the number of consecutive runs for which
it evaluates to false.  In the example above, we have 2 runs.  These are the
block of 5 consecutive heads at the beginning and the 5 consecutive tails
at the end.

Alternatives are Alt.less, meaning that less runs than expected have been
observed and data for which positive() is true tends to cluster,
Alt.greater, which means that more runs than expected have been observed
and data for which positive() is true tends to not cluster even moreso than
expected by chance, and Alt.twoSided, meaning that elements for which
positive() is true cluster as much as expected by chance.

Bugs:  No exact calculation of the P-value.  Asymptotic approximation only.

References:  http://en.wikipedia.org/wiki/Runs_test
*/
double runsTest(alias positive = "a > 0", T)(T obs, Alt alt = Alt.twoSided)
if(isIterable!(T)) {
    RunsTest!(positive, ForeachType!(T)) r;
    foreach(elem; obs) {
        r.put(elem);
    }
    return r.p(alt);
}

unittest {
    // Values from R lawstat package, for which "a < median(data)" is
    // hard-coded as the equivalent to positive().  The median of this data
    // is 0.5, so everything works.
    immutable int[] data = [1,0,0,0,1,1,0,0,1,0,1,0,1,0,1,1,1,0,0,1].idup;
    assert(approxEqual(runsTest(data), 0.3581));
    assert(approxEqual(runsTest(data, Alt.less), 0.821));
    assert(approxEqual(runsTest(data, Alt.greater), 0.1791));
}

/**
Runs test as in runsTest(), except calculates online instead of from stored
array elements.
*/
struct RunsTest(alias positive = "a > 0", T) {
private:
    uint nPos;
    uint nNeg;
    uint nRun;
    bool lastPos;

    alias unaryFun!(positive) pos;

public:

    ///
    void put(T elem) {
        bool curPos = pos(elem);
        if(nRun == 0) {
            nRun = 1;
            if(curPos) {
                nPos++;
            } else {
                nNeg++;
            }
        } else if(pos(elem)) {
            nPos++;
            if(!lastPos) {
                nRun++;
            }
        } else {
            nNeg++;
            if(lastPos) {
                nRun++;
            }
        }
        lastPos = curPos;
    }

    ///
    uint nRuns() {
        return nRun;
    }

    ///
    double p(Alt alt = Alt.twoSided) {
        uint N = nPos + nNeg;
        double expected = 2.0 * nPos * nNeg / N + 1;
        double sd = sqrt((expected - 1) * (expected - 2) / (N - 1));
        if(alt == Alt.less) {
            return normalCDF(nRun, expected, sd);
        } else if(alt == Alt.greater) {
            return normalCDFR(nRun, expected, sd);
        } else {
            return 2 * ((nRun < expected) ?
                        normalCDF(nRun, expected, sd) :
                        normalCDFR(nRun, expected, sd));
        }
    }
}

deprecated {
    // Aliases for old names for correlation tests.
    alias pearsonCorTest pcorTest;
    alias spearmanCorTest scorTest;
    alias kendallCorTest kcorTest;
}

/**
Tests the hypothesis that the Pearson correlation between two ranges is
different from some 0.  Alternatives are Alt.less
(pearsonCor(range1, range2) < 0), Alt.greater (pearsonCor(range1, range2)
 0) and Alt.twoSided (pearsonCor(range1, range2) != 0).

Returns:  A ConfInt of the estimated Pearson correlation of the two ranges,
the P-value against the given alternative, and the confidence interval of
the correlation at the level specified by confLevel.

References:  http://en.wikipedia.org/wiki/Pearson_correlation
*/
ConfInt pearsonCorTest(T, U)(
    T range1,
    U range2,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) if(doubleInput!(T) && doubleInput!(U)) {
    enforceConfidence(confLevel);

    auto pearsonRes = dstats.cor.pearsonCor(range1, range2);
    if(isNaN(pearsonRes.cor)) {
        return ConfInt.init;
    }

    return pearsonCorTest(pearsonRes.cor, pearsonRes.N, alt, confLevel);
}

/**
Same as overload, but uses pre-computed correlation coefficient and sample
size instead of computing them.

Note:  This is a template only because of DMD Bug 2972.
 */
ConfInt pearsonCorTest()(
    double cor,
    double N,
    Alt alt = Alt.twoSided,
    double confLevel = 0.95
) {
    dstatsEnforce(N >= 0, "N must be >= 0 for pearsonCorTest.");
    dstatsEnforce(cor > -1.0 || approxEqual(cor, -1.0),
        "Correlation must be between 0, 1.");
    dstatsEnforce(cor < 1.0 || approxEqual(cor, 1.0),
         "Correlation must be between 0, 1.");
    enforceConfidence(confLevel);

    immutable double denom = sqrt((1 - cor * cor) / (N - 2));
    immutable double t = cor / denom;
    ConfInt ret;
    ret.testStat = cor;

    double sqN, z;
    if(confLevel > 0) {
        sqN = sqrt(N - 3.0);
        z = sqN * atanh(cor);
    }

    final switch(alt) {
        case Alt.none :
            return ret;
        case Alt.twoSided:
            ret.p = (abs(cor) >= 1) ? 0 :
                2 * ((t < 0) ? studentsTCDF(t, N - 2) : studentsTCDFR(t, N - 2));

            if(confLevel > 0) {
                double deltaZ = invNormalCDF(0.5 * (1 - confLevel));
                ret.lowerBound = tanh((z + deltaZ) / sqN);
                ret.upperBound = tanh((z - deltaZ) / sqN);
            } else {
                ret.lowerBound = cor;
                ret.upperBound = cor;
            }

            break;
        case Alt.less:
            if(cor >= 1) {
                ret.p = 1;
            } else if(cor <= -1) {
                ret.p = 0;
            } else {
                ret.p = studentsTCDF(t, N - 2);
            }

            if(confLevel > 0) {
                double deltaZ = invNormalCDF(1 - confLevel);
                ret.lowerBound = -1;
                ret.upperBound = tanh((z - deltaZ) / sqN);
            } else {
                ret.lowerBound = -1;
                ret.upperBound = cor;
            }

            break;
        case Alt.greater:
            if(cor >= 1) {
                ret.p = 0;
            } else if(cor <= -1) {
                ret.p = 1;
            } else {
                ret.p = studentsTCDFR(t, N - 2);
            }

            if(confLevel > 0) {
                double deltaZ = invNormalCDF(1 - confLevel);
                ret.lowerBound = tanh((z + deltaZ) / sqN);
                ret.upperBound = 1;
            } else {
                ret.lowerBound = cor;
                ret.upperBound = 1;
            }

            break;
    }
    return ret;
}

unittest {
    // Values from R.
    auto t1 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.twoSided);
    auto t2 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.less);
    auto t3 = pearsonCorTest([1,2,3,4,5].dup, [2,1,4,3,5].dup, Alt.greater);

    assert(approxEqual(t1.testStat, 0.8));
    assert(approxEqual(t2.testStat, 0.8));
    assert(approxEqual(t3.testStat, 0.8));

    assert(approxEqual(t1.p, 0.1041));
    assert(approxEqual(t2.p, 0.948));
    assert(approxEqual(t3.p, 0.05204));

    assert(approxEqual(t1.lowerBound, -0.2796400));
    assert(approxEqual(t3.lowerBound, -0.06438567));
    assert(approxEqual(t2.lowerBound, -1));

    assert(approxEqual(t1.upperBound, 0.9861962));
    assert(approxEqual(t2.upperBound, 0.9785289));
    assert(approxEqual(t3.upperBound, 1));

    // Test special case hack for cor = +- 1.
    uint[] myArr = [1,2,3,4,5];
    uint[] myArrReverse = myArr.dup;
    reverse(myArrReverse);

    auto t4 = pearsonCorTest(myArr, myArr, Alt.twoSided);
    auto t5 = pearsonCorTest(myArr, myArr, Alt.less);
    auto t6 = pearsonCorTest(myArr, myArr, Alt.greater);
    assert(approxEqual(t4.testStat, 1));
    assert(approxEqual(t4.p, 0));
    assert(approxEqual(t5.p, 1));
    assert(approxEqual(t6.p, 0));

    auto t7 = pearsonCorTest(myArr, myArrReverse, Alt.twoSided);
    auto t8 = pearsonCorTest(myArr, myArrReverse, Alt.less);
    auto t9 = pearsonCorTest(myArr, myArrReverse, Alt.greater);
    assert(approxEqual(t7.testStat, -1));
    assert(approxEqual(t7.p, 0));
    assert(approxEqual(t8.p, 0));
    assert(approxEqual(t9.p, 1));
}

/**
Tests the hypothesis that the Spearman correlation between two ranges is
different from some 0.  Alternatives are
Alt.less (spearmanCor(range1, range2) < 0), Alt.greater (spearmanCor(range1, range2)
> 0) and Alt.twoSided (spearmanCor(range1, range2) != 0).

Returns:  A TestRes containing the Spearman correlation coefficient and
the P-value for the given alternative.

Bugs:  Exact P-value computation not yet implemented.  Uses asymptotic
approximation only.  This is good enough for most practical purposes given
reasonably large N, but is not perfectly accurate.  Not valid for data with
very large amounts of ties.
*/
TestRes spearmanCorTest(T, U)(T range1, U range2, Alt alt = Alt.twoSided)
if(isInputRange!(T) && isInputRange!(U) &&
is(typeof(range1.front < range1.front) == bool) &&
is(typeof(range2.front < range2.front) == bool)) {

    static if(!hasLength!T) {
        auto alloc = newRegionAllocator();
        auto r1 = alloc.array(range1);
    } else {
        alias range1 r1;
    }
    immutable double N = r1.length;

    return pearsonCorTest(dstats.cor.spearmanCor(range1, range2), N, alt, 0);
}

unittest {
    // Values from R.
    int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
    int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];
    auto t1 = spearmanCorTest(arr1, arr2, Alt.twoSided);
    auto t2 = spearmanCorTest(arr1, arr2, Alt.less);
    auto t3 = spearmanCorTest(arr1, arr2, Alt.greater);

    assert(approxEqual(t1.testStat, -0.1769406));
    assert(approxEqual(t2.testStat, -0.1769406));
    assert(approxEqual(t3.testStat, -0.1769406));

    assert(approxEqual(t1.p, 0.4429));
    assert(approxEqual(t3.p, 0.7785));
    assert(approxEqual(t2.p, 0.2215));
}

/**
Tests the hypothesis that the Kendall Tau-b between two ranges is
different from 0.  Alternatives are Alt.less (kendallCor(range1, range2) < 0),
Alt.greater (kendallCor(range1, range2) > 0) and Alt.twoSided
(kendallCor(range1, range2) != 0).


exactThresh controls the maximum length of the range for which exact P-value
computation is used.  The default is 50.  Exact calculation is never used
when ties are present because it is not computationally feasible.
Do not set this higher than 100, as it will be very slow
and the asymptotic approximation is pretty good at even a fraction of this
size.

Note:

As an optimization, when a range is a SortedRange with predicate "a < b",
it is assumed already sorted and not sorted a second time by this function.
This is useful when applying this function multiple times with one of the
arguments the same every time.

Returns:  A TestRes containing the Kendall correlation coefficient and
the P-value for the given alternative.

References:  StackOverflow Question 948341 (http://stackoverflow.com/questions/948341)

The Variance of Tau When Both Rankings Contain Ties.  M.G. Kendall.
Biometrika, Vol 34, No. 3/4 (Dec., 1947), pp. 297-298
 */
TestRes kendallCorTest(T, U)(
    T range1,
    U range2,
    Alt alt = Alt.twoSided,
    uint exactThresh = 50
) if(isInputRange!(T) && isInputRange!(U)) {
    auto alloc = newRegionAllocator();

    static if(dstats.cor.isDefaultSorted!T) {
        alias range1 i1d;
    } else {
        auto i1d = alloc.array(range1);
    }

    static if(dstats.cor.isDefaultSorted!U) {
        alias range2 i2d;
    } else {
        auto i2d = alloc.array(range2);
    }

    immutable res = dstats.cor.kendallCorDestructiveLowLevel(i1d, i2d, true);
    immutable double n = i1d.length;

    immutable double var =
          (2.0 / 9) * n * (n - 1) * (2 * n + 5)
        - (2.0 / 9) * res.tieCorrectT1
        - (2.0 / 9) * res.tieCorrectU1
        + (4 / (9 * n * (n - 1) * (n - 2))) * res.tieCorrectT2 * res.tieCorrectU2
        + 2 / (n * (n - 1)) * res.tieCorrectT3 * res.tieCorrectU3;

    // Need the / 2 to change C, as used in Kendall's paper to S, as used here.
    immutable double sd = sqrt(var) / 2;

    enum double cc = 1;
    auto tau = res.tau;
    auto s = res.s;

    immutable bool noTies = res.tieCorrectT1 == 0 && res.tieCorrectU1 == 0;

    if(noTies && n <= exactThresh) {
        // More than uint.max data points for exact calculation is
        // not plausible.
        assert(i1d.length < uint.max);
        immutable N = cast(uint) i1d.length;
        immutable nSwaps = (N * (N - 1) / 2 - cast(uint) s) / 2;
        return TestRes(tau, kendallCorExactP(N, nSwaps, alt));
    }

    final switch(alt) {
        case Alt.none :
            return TestRes(tau);
        case Alt.twoSided:
            if(abs(s) <= cc) {
                return TestRes(tau, 1);
            } else if(s < 0) {
                return TestRes(tau, 2 * normalCDF(s + cc, 0, sd));
            } else {
                assert(s > 0);
                return TestRes(tau, 2 * normalCDFR(s - cc, 0, sd));
            }
            assert(0);

        case Alt.less:
            return TestRes(tau, normalCDF(s + cc, 0, sd));
        case Alt.greater:
            return TestRes(tau, normalCDFR(s - cc, 0, sd));
    }
}

// Dynamic programming algorithm for computing exact Kendall tau P-values.
// Thanks to ShreevatsaR from StackOverflow.
private double kendallCorExactP(uint N, uint swaps, Alt alt) {
    uint maxSwaps = N * (N - 1) / 2;
    assert(swaps <= maxSwaps);
    immutable expectedSwaps = cast(ulong) N * (N - 1) * 0.25;
    if(alt == Alt.greater) {
        if(swaps > expectedSwaps) {
            if(swaps == maxSwaps) {
                return 1;
            }
            return 1.0 - kendallCorExactP(N, maxSwaps - swaps - 1, Alt.greater);
        }
    } else if(alt == Alt.less) {
        if(swaps == 0) {
            return 1;
        }
        return kendallCorExactP(N, maxSwaps - swaps + 0, Alt.greater);
    } else if(alt == Alt.twoSided) {
        if(swaps < expectedSwaps) {
            return min(1, 2 * kendallCorExactP(N, swaps, Alt.greater));
        } else if(swaps > expectedSwaps) {
            return min(1, 2 * kendallCorExactP(N, swaps, Alt.less));
        } else {
            return 1;
        }
    } else {  // Alt.none
        return double.nan;
    }

    /* This algorithm was obtained from Question 948341 on stackoverflow.com
     * and works as follows:
     *
     * swaps is the number of swaps that would be necessary in a bubble sort
     * to sort one list in the same order as the other.  N is the sample size.
     * We want to find the number of ways that we could get a bubble sort
     * distance of at least swaps, and divide it by the total number of
     * permutations, pElem.
     *
     * The number of swaps necessary to sort a list is equivalent to the
     * number of inversions in the list, i.e. where i > j, but
     * list[i] < list[j].  This is a bottom-up dynamic programming algorithm
     * based on this principle.
     *
     * Assume c(N, k) is the number of permutations that require <= swaps
     * inversions.
     * We use the recurrence relation:
     * When k ≤ N - 1, c(N,k) = c(N,k-1) + c(N-1,k)
     * When k ≥ N,   c(N,k) = c(N,k-1) + c(N-1,k) - c(N-1,k-N)
     *
     * We also divide every value by the constant N! to turn this count into a
     * probability.
     */

    immutable double pElem = exp(-logFactorial(N));
    auto alloc = newRegionAllocator();
    double[] cur = alloc.uninitializedArray!(double[])(swaps + 1);
    double[] prev = alloc.uninitializedArray!(double[])(swaps + 1);

    prev[] = pElem;
    cur[0] = pElem;
    foreach(k; 1..N + 1) {
        immutable uint nSwapsPossible = k * (k - 1) / 2;
        immutable uint upTo = min(swaps, nSwapsPossible) + 1;
        foreach(j; 1..upTo) {
            if(j < k) {
                cur[j] = prev[j] + cur[j - 1];
            } else {
                cur[j] = prev[j] - prev[j - k] + cur[j - 1];
            }
        }
        cur[upTo..$] = cur[upTo - 1];
        swap(cur, prev);
    }

    return prev[$ - 1];
}

unittest {
    /* Values from R, with continuity correction enabled.  Note that large
     * one-sided inexact P-values are commented out because R seems to have a
     * slightly different interpretation of the proper continuity correction
     * than this library.  This library corrects the z-score in the direction
     * that would make the test more conservative.  R corrects towards zero.
     * I can't find a reference to support either one, but empirically it seems
     * like correcting towards more conservative results more accurately mirrors
     * the results of the exact test.  This isn't too big a deal anyhow since:
     *
     * 1.  The difference is small.
     * 2.  It only occurs on results that are very far from significance
     *     (P > 0.5).
     */
    int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
    int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];
    auto t1 = kendallCorTest(arr1, arr2, Alt.twoSided);
    auto t2 = kendallCorTest(arr1, arr2, Alt.less);
    auto t3 = kendallCorTest(arr1, arr2, Alt.greater);

    assert(approxEqual(t1.testStat, -.1448010));
    assert(approxEqual(t2.testStat, -.1448010));
    assert(approxEqual(t3.testStat, -.1448010));

    assert(approxEqual(t1.p, 0.3923745));
    //assert(approxEqual(t3.p, 0.8038127));
    assert(approxEqual(t2.p, 0.1961873));

    // Now, test the case of ties in both arrays.
    arr1 = [1,1,1,2,2,3,4,5,5,6];
    arr2 = [1,1,2,3,4,5,5,5,5,6];
    assert(approxEqual(kendallCorTest(arr1, arr2, Alt.twoSided).p, 0.001216776));
    //assert(approxEqual(kendallCorTest(arr1, arr2, Alt.less).p, 0.9993916));
    assert(approxEqual(kendallCorTest(arr1, arr2, Alt.greater).p, 0.0006083881));

    arr1 = [1,1,1,2,2,2,3,3,3,4,4,4,5,5,5];
    arr2 = [1,1,1,3,3,3,2,2,2,5,5,5,4,4,4];
    assert(approxEqual(kendallCorTest(arr1, arr2).p, 0.006123));
    assert(approxEqual(kendallCorTest(assumeSorted(arr1), arr2).p, 0.006123));

    // Test the exact stuff.  Still using values from R.
    uint[] foo = [1,2,3,4,5];
    uint[] bar = [1,2,3,5,4];
    uint[] baz = [5,3,1,2,4];

    assert(approxEqual(kendallCorTest(foo, foo).p, 0.01666666));
    assert(approxEqual(kendallCorTest(foo, foo, Alt.greater).p, 0.008333333));
    assert(approxEqual(kendallCorTest(foo, foo, Alt.less).p, 1));

    assert(approxEqual(kendallCorTest(foo, bar).p, 0.083333333));
    assert(approxEqual(kendallCorTest(foo, bar, Alt.greater).p, 0.041666667));
    assert(approxEqual(kendallCorTest(foo, bar, Alt.less).p, 0.9917));

    assert(approxEqual(kendallCorTest(foo, baz).p, 0.8167));
    assert(approxEqual(kendallCorTest(foo, baz, Alt.greater).p, 0.7583));
    assert(approxEqual(kendallCorTest(foo, baz, Alt.less).p, .4083));

    assert(approxEqual(kendallCorTest(bar, baz).p, 0.4833));
    assert(approxEqual(kendallCorTest(bar, baz, Alt.greater).p, 0.8833));
    assert(approxEqual(kendallCorTest(bar, baz, Alt.less).p, 0.2417));

    // A little monte carlo unittesting.  For large ranges, the deviation
    // between the exact and approximate version should be extremely small.
    foreach(i; 0..100) {
        uint nToTake = uniform(15, 65);
        auto lhs = array(take(randRange!rNormal(0, 1), nToTake));
        auto rhs = array(take(randRange!rNormal(0, 1), nToTake));
        if(i & 1) {
            lhs[] += rhs[] * 0.2;  // Make sure there's some correlation.
        } else {
            lhs[] -= rhs[] * 0.2;
        }
        double exact = kendallCorTest(lhs, rhs).p;
        double approx = kendallCorTest(lhs, rhs, Alt.twoSided, 0).p;
        assert(abs(exact - approx) < 0.01);

        exact = kendallCorTest(lhs, rhs, Alt.greater).p;
        approx = kendallCorTest(lhs, rhs, Alt.greater, 0).p;
        assert(abs(exact - approx) < 0.01);

        exact = kendallCorTest(lhs, rhs, Alt.less).p;
        approx = kendallCorTest(lhs, rhs, Alt.less, 0).p;
        assert(abs(exact - approx) < 0.01);
    }
}

/**A test for normality of the distribution of a range of values.  Based on
 * the assumption that normally distributed values will have a sample skewness
 * and sample kurtosis very close to zero.
 *
 * Returns:  A TestRes with the K statistic, which is Chi-Square distributed
 * with 2 degrees of freedom under the null, and the P-value for the alternative
 * that the data has skewness and kurtosis not equal to zero against the null
 * that skewness and kurtosis are near zero.  A normal distribution always has
 * skewness and kurtosis that converge to zero as sample size goes to infinity.
 *
 * Notes:  Contrary to popular belief, tests for normality should usually
 * not be used to deterimine whether T-tests are valid.  If the sample size is
 * large, T-tests are valid regardless of the distribution due to the central
 * limit theorem.  If the sample size is small, a test for normality will
 * likely not be very powerful, and a priori knowledge or simple inspection
 * of the data is often a better idea.
 *
 * References:
 * D'Agostino, Ralph B., Albert Belanger, and Ralph B. D'Agostino, Jr.
 * "A Suggestion for Using Powerful and Informative Tests of Normality",
 * The American Statistician, Vol. 44, No. 4. (Nov., 1990), pp. 316-321.
 */
TestRes dAgostinoK(T)(T range)
if(doubleIterable!(T)) {
    // You are not meant to understand this.  I sure don't.  I just implemented
    // these formulas off of Wikipedia, which got them from:

    // D'Agostino, Ralph B., Albert Belanger, and Ralph B. D'Agostino, Jr.
    // "A Suggestion for Using Powerful and Informative Tests of Normality",
    // The American Statistician, Vol. 44, No. 4. (Nov., 1990), pp. 316-321.

    // Amazing.  I didn't even realize things this complicated were possible
    // in 1990, before widespread computer algebra systems.

    // Notation from Wikipedia.  Keeping same notation for simplicity.
    double sqrtb1 = void, b2 = void, n = void;
    {
        auto summ = summary(range);
        sqrtb1 = summ.skewness;
        b2 = summ.kurtosis + 3;
        n = summ.N;
    }

    // Calculate transformed skewness.
    double Y = sqrtb1 * sqrt((n + 1) * (n + 3) / (6 * (n - 2)));
    double beta2b1Numer = 3 * (n * n + 27 * n - 70) * (n + 1) * (n + 3);
    double beta2b1Denom = (n - 2) * (n + 5) * (n + 7) * (n + 9);
    double beta2b1 = beta2b1Numer / beta2b1Denom;
    double Wsq = -1 + sqrt(2 * (beta2b1 - 1));
    double delta = 1.0 / sqrt(log(sqrt(Wsq)));
    double alpha = sqrt( 2.0 / (Wsq - 1));
    double Zb1 = delta * log(Y / alpha + sqrt(pow(Y / alpha, 2) + 1));

    // Calculate transformed kurtosis.
    double Eb2 = 3 * (n - 1) / (n + 1);
    double sigma2b2 = (24 * n * (n - 2) * (n - 3)) / (
        (n + 1) * (n + 1) * (n + 3) * (n + 5));
    double x = (b2 - Eb2) / sqrt(sigma2b2);

    double sqBeta1b2 = 6 * (n * n - 5 * n + 2) / ((n + 7) * (n + 9)) *
         sqrt((6 * (n + 3) * (n + 5)) / (n * (n - 2) * (n - 3)));
    double A = 6 + 8 / sqBeta1b2 * (2 / sqBeta1b2 + sqrt(1 + 4 / (sqBeta1b2 * sqBeta1b2)));
    double Zb2 = ((1 - 2 / (9 * A)) -
        cbrt((1 - 2 / A) / (1 + x * sqrt(2 / (A - 4)))) ) *
        sqrt(9 * A / 2);

    double K2 = Zb1 * Zb1 + Zb2 * Zb2;

    if(isNaN(K2)) {
        return TestRes(double.nan, double.nan);
    }

    return TestRes(K2, chiSquareCDFR(K2, 2));
}

unittest {
    // Values from R's fBasics package.
    int[] arr1 = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20];
    int[] arr2 = [8,6,7,5,3,0,9,8,6,7,5,3,0,9,3,6,2,4,3,6,8];

    auto r1 = dAgostinoK(arr1);
    auto r2 = dAgostinoK(arr2);

    assert(approxEqual(r1.testStat, 3.1368));
    assert(approxEqual(r1.p, 0.2084));

    assert(approxEqual(r2.testStat, 1.1816));
    assert(approxEqual(r2.p, 0.5539));
}

/**Fisher's method of meta-analyzing a set of P-values to determine whether
 * there are more significant results than would be expected by chance.
 * Based on a chi-square statistic for the sum of the logs of the P-values.
 *
 * Returns:  A TestRes containing the chi-square statistic and a P-value for
 * the alternative hypothesis that more small P-values than would be expected
 * by chance are present against the alternative that the distribution of
 * P-values is uniform or enriched for large P-values.
 *
 * References:  Fisher, R. A. (1948) "Combining independent tests of
 * significance", American Statistician, vol. 2, issue 5, page 30.
 * (In response to Question 14)
 */
TestRes fishersMethod(R)(R pVals)
if(doubleInput!R) {
    double chiSq = 0;
    uint df = 0;
    foreach(pVal; pVals) {
        dstatsEnforce(pVal >= 0 && pVal <= 1,
            "P-values must be between 0, 1 for Fisher's Method.");
        chiSq += log( cast(double) pVal);
        df += 2;
    }
    chiSq *= -2;
    return TestRes(chiSq, chiSquareCDFR(chiSq, df));
}

unittest {
    // First, basic sanity check.  Make sure w/ one P-value, we get back that
    // P-value.
    for(double p = 0.01; p < 1; p += 0.01) {
        assert(approxEqual(fishersMethod([p].dup).p, p));
    }
    float[] ps = [0.739, 0.0717, 0.01932, 0.03809];
    auto res = fishersMethod(ps);
    assert(approxEqual(res.testStat, 20.31));
    assert(res.p < 0.01);
}

/// For falseDiscoveryRate.
enum Dependency : bool {
    /**
    Assume that dependency among hypotheses may exist.  (More conservative,
    Benjamini-Yekutieli procedure.)
    */
    yes = true,

    /**
    Assume hypotheses are independent.  (Less conservative, Benjamine-
    Hochberg procedure.)
    */
    no = false
}

/**
Computes the false discovery rate statistic given a list of
p-values, according to Benjamini and Hochberg (1995) (independent) or
Benjamini and Yekutieli (2001) (dependent).  The Dependency parameter
controls whether hypotheses are assumed to be independent, or whether
the more conservative assumption that they are correlated must be made.

Returns:
An array of adjusted P-values with indices corresponding to the order of
the P-values in the input data.

References:
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate:
a practical and powerful approach to multiple testing. Journal of the Royal
Statistical Society Series B, 57, 289-200

Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery
rate in multiple testing under dependency. Annals of Statistics 29, 1165-1188.
 */
float[] falseDiscoveryRate(T)(T pVals, Dependency dep = Dependency.no)
if(doubleInput!(T)) {
    auto qVals = array(map!(to!float)(pVals));

    double C = 1;
    if(dep == Dependency.yes) {
        foreach(i; 2..qVals.length + 1) {
            C += 1.0 / i;
        }
    }

    auto alloc = newRegionAllocator();
    auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);
    foreach(i, ref elem; perm) {
        elem = i;
    }

    qsort(qVals, perm);

    foreach(i, ref q; qVals) {
        q = min(1.0f, q * C * cast(double) qVals.length / (cast(double) i + 1));
    }

    float smallestSeen = float.max;
    foreach(ref q; retro(qVals)) {
        if(q < smallestSeen) {
            smallestSeen = q;
        } else {
            q = smallestSeen;
        }
    }

    qsort(perm, qVals);  //Makes order of qVals correspond to input.
    return qVals;
}

unittest {
    // Comparing results to R.
    auto pVals = [.90, .01, .03, .03, .70, .60, .01].dup;
    auto qVals = falseDiscoveryRate(pVals);
    alias approxEqual ae;
    assert(ae(qVals[0], .9));
    assert(ae(qVals[1], .035));
    assert(ae(qVals[2], .052));
    assert(ae(qVals[3], .052));
    assert(ae(qVals[4], .816666666667));
    assert(ae(qVals[5], .816666666667));
    assert(ae(qVals[6], .035));

    auto p2 = [.1, .02, .6, .43, .001].dup;
    auto q2 = falseDiscoveryRate(p2);
    assert(ae(q2[0], .16666666));
    assert(ae(q2[1], .05));
    assert(ae(q2[2], .6));
    assert(ae(q2[3], .5375));
    assert(ae(q2[4], .005));

    // Dependent case.
    qVals = falseDiscoveryRate(pVals, Dependency.yes);
    assert(ae(qVals[0], 1));
    assert(ae(qVals[1], .09075));
    assert(ae(qVals[2], .136125));
    assert(ae(qVals[3], .136125));
    assert(ae(qVals[4], 1));
    assert(ae(qVals[5], 1));
    assert(ae(qVals[6], .09075));

    q2 = falseDiscoveryRate(p2, Dependency.yes);
    assert(ae(q2[0], .38055555));
    assert(ae(q2[1], .1141667));
    assert(ae(q2[2], 1));
    assert(ae(q2[3], 1));
    assert(ae(q2[4], .01141667));
}

/**Uses the Hochberg procedure to control the familywise error rate assuming
 * that hypothesis tests are independent.  This is more powerful than
 * Holm-Bonferroni correction, but requires the independence assumption.
 *
 * Returns:
 * An array of adjusted P-values with indices corresponding to the order of
 * the P-values in the input data.
 *
 * References:
 * Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of
 * significance. Biometrika, 75, 800-803.
 */
float[] hochberg(T)(T pVals)
if(doubleInput!(T)) {
    auto qVals = array(map!(to!float)(pVals));

    auto alloc = newRegionAllocator();
    auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);
    foreach(i, ref elem; perm)
        elem = i;

    qsort(qVals, perm);

    foreach(i, ref q; qVals) {
        dstatsEnforce(q >= 0 && q <= 1,
            "P-values must be between 0, 1 for hochberg.");
        q = min(1.0f, q * (cast(double) qVals.length - i));
    }

    float smallestSeen = float.max;
    foreach(ref q; retro(qVals)) {
        if(q < smallestSeen) {
            smallestSeen = q;
        } else {
            q = smallestSeen;
        }
    }

    qsort(perm, qVals);  //Makes order of qVals correspond to input.
    return qVals;
}

unittest {
    alias approxEqual ae;
    auto q = hochberg([0.01, 0.02, 0.025, 0.9].dup);
    assert(ae(q[0], 0.04));
    assert(ae(q[1], 0.05));
    assert(ae(q[2], 0.05));
    assert(ae(q[3], 0.9));

    auto p2 = [.1, .02, .6, .43, .001].dup;
    auto q2 = hochberg(p2);
    assert(ae(q2[0], .3));
    assert(ae(q2[1], .08));
    assert(ae(q2[2], .6));
    assert(ae(q2[3], .6));
    assert(ae(q2[4], .005));
}

/**Uses the Holm-Bonferroni method to adjust a set of P-values in a way that
 * controls the familywise error rate (The probability of making at least one
 * Type I error).  This is basically a less conservative version of
 * Bonferroni correction that is still valid for arbitrary assumptions and
 * controls the familywise error rate.  Therefore, there aren't too many good
 * reasons to use regular Bonferroni correction instead.
 *
 * Returns:
 * An array of adjusted P-values with indices corresponding to the order of
 * the P-values in the input data.
 *
 * References:
 * Holm, S. (1979). A simple sequentially rejective multiple test procedure.
 * Scandinavian Journal of Statistics, 6, 65-70.
 */
float[] holmBonferroni(T)(T pVals)
if(doubleInput!(T)) {
    auto alloc = newRegionAllocator();

    auto qVals = array(map!(to!float)(pVals));
    auto perm = alloc.uninitializedArray!(size_t[])(qVals.length);

    foreach(i, ref elem; perm) {
        elem = i;
    }

    qsort(qVals, perm);

    foreach(i, ref q; qVals) {
        dstatsEnforce(q >= 0 && q <= 1,
            "P-values must be between 0, 1 for holmBonferroni.");
        q = min(1.0, q * (cast(double) qVals.length - i));
    }

    foreach(i; 1..qVals.length) {
        if(qVals[i] < qVals[i - 1]) {
            qVals[i] = qVals[i - 1];
        }
    }

    qsort(perm, qVals);  //Makes order of qVals correspond to input.
    return qVals;
}

unittest {
    // Values from R.
    auto ps = holmBonferroni([0.001, 0.2, 0.3, 0.4, 0.7].dup);
    alias approxEqual ae;
    assert(ae(ps[0], 0.005));
    assert(ae(ps[1], 0.8));
    assert(ae(ps[2], 0.9));
    assert(ae(ps[3], 0.9));
    assert(ae(ps[4], 0.9));

    ps = holmBonferroni([0.3, 0.1, 0.4, 0.1, 0.5, 0.9].dup);
    assert(ps == [1f, 0.6f, 1f, 0.6f, 1f, 1f]);
}

// old unconstrained approxEqual to work around https://issues.dlang.org/show_bug.cgi?id=18287
static if (__VERSION__ == 2078) private
{
    /**
       Computes whether two values are approximately equal, admitting a maximum
       relative difference, and a maximum absolute difference.
       Params:
            lhs = First item to compare.
            rhs = Second item to compare.
            maxRelDiff = Maximum allowable difference relative to `rhs`.
            maxAbsDiff = Maximum absolute difference.
       Returns:
           `true` if the two items are approximately equal under either criterium.
           If one item is a range, and the other is a single value, then the result
           is the logical and-ing of calling `approxEqual` on each element of the
           ranged item against the single item. If both items are ranges, then
           `approxEqual` returns `true` if and only if the ranges have the same
           number of elements and if `approxEqual` evaluates to `true` for each
           pair of elements.
     */
    bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff, V maxAbsDiff = 1e-5)
    {
        import std.range.primitives : empty, front, isInputRange, popFront;
        static if (isInputRange!T)
        {
            static if (isInputRange!U)
            {
                // Two ranges
                for (;; lhs.popFront(), rhs.popFront())
                {
                    if (lhs.empty) return rhs.empty;
                    if (rhs.empty) return lhs.empty;
                    if (!approxEqual(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
                        return false;
                }
            }
            else static if (isIntegral!U)
            {
                // convert rhs to real
                return approxEqual(lhs, real(rhs), maxRelDiff, maxAbsDiff);
            }
            else
            {
                // lhs is range, rhs is number
                for (; !lhs.empty; lhs.popFront())
                {
                    if (!approxEqual(lhs.front, rhs, maxRelDiff, maxAbsDiff))
                        return false;
                }
                return true;
            }
        }
        else
        {
            static if (isInputRange!U)
            {
                // lhs is number, rhs is range
                for (; !rhs.empty; rhs.popFront())
                {
                    if (!approxEqual(lhs, rhs.front, maxRelDiff, maxAbsDiff))
                        return false;
                }
                return true;
            }
            else static if (isIntegral!T || isIntegral!U)
            {
                // convert both lhs and rhs to real
                return approxEqual(real(lhs), real(rhs), maxRelDiff, maxAbsDiff);
            }
            else
            {
                // two numbers
                //static assert(is(T : real) && is(U : real));
                if (rhs == 0)
                {
                    return fabs(lhs) <= maxAbsDiff;
                }
                static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
                {
                    if (lhs == lhs.infinity && rhs == rhs.infinity ||
                        lhs == -lhs.infinity && rhs == -rhs.infinity) return true;
                }
                return fabs((lhs - rhs) / rhs) <= maxRelDiff
                    || maxAbsDiff != 0 && fabs(lhs - rhs) <= maxAbsDiff;
            }
        }
    }

    /**
       Returns $(D approxEqual(lhs, rhs, 1e-2, 1e-5)).
     */
    bool approxEqual(T, U)(T lhs, U rhs)
    {
        return approxEqual(lhs, rhs, 1e-2, 1e-5);
    }
}