root/trunk/tests.d

/**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.algorithm, std.functional, std.range, std.conv, std.math, std.traits,
       std.exception, std.typetuple;

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

private static import dstats.infotheory;

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

/**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,

    // These are kept for compatibility with the old style, are intentionally
    // lacking DDoc and may eventually be deprecated/removed.
    TWOSIDE = twoSided,
    LESS = less,
    GREATER = greater,
    NONE = 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 {
    ///
    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(T) 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].dup, [1,3,4,5,7,9].dup);
    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].dup, [1,3,4,5,7,9].dup, 0, Alt.less),
           0.1173));
    assert(approxEqual(studentsTTest([1,2,3,4,5].dup, [1,3,4,5,7,9].dup, 0, Alt.greater),
           0.8827));
    auto t2 = studentsTTest([1,3,5,7,9,11].dup, [2,2,1,3,4].dup, 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].dup, [2,2,1,3,4].dup, 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].dup, [2,2,1,3,4].dup, 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].dup), [1,3,4,5,7,9].dup, 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].dup, summary([1,3,4,5,7,9].dup), -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].dup, [1,3,4,5,7,9].dup, 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].dup, [2,2,1,3,4].dup).p, 0.06616));
    assert(approxEqual(welchTTest([1,3,5,7,9,11].dup, [2,2,1,3,4].dup, 0,
        Alt.less).p, 0.967));
    assert(approxEqual(welchTTest([1,3,5,7,9,11].dup, [2,2,1,3,4].dup, 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.
 *
 * 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) {
        double diff = cast(double) before.front - cast(double) after.front;
        before.popFront;
        after.popFront;
        msd.put(diff);
    }

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

/**Compute the test directly from summary statistics of the differences between
 * corresponding samples.
 *
 * Examples:
 * ---
 * float[] data1 = getSomeDataSet();
 * float[] data2 = getSomeOtherDataSet();
 * assert(data1.length == data2.length);
 *
 * // 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, 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].dup, [2,3,5,5,6].dup, 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].dup, [2,3,5,5,6].dup, 0, Alt.less).p, 0.0889));
    assert(approxEqual(pairedTTest([3,2,3,4,5].dup, [2,3,5,5,6].dup, 0, Alt.greater).p, 0.9111));
    assert(approxEqual(pairedTTest([3,2,3,4,5].dup, [2,3,5,5,6].dup, 0, Alt.twoSided).p, 0.1778));
    assert(approxEqual(pairedTTest([3,2,3,4,5].dup, [2,3,5,5,6].dup, 1, Alt.less).p, 0.01066));
    assert(approxEqual(pairedTTest([3,2,3,4,5].dup, [2,3,5,5,6].dup, 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) {
        mixin(newFrame);
        auto data = tempdup(dataIn[0]);
        auto withins = newStack!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 = newStack!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))) {
        mixin(newFrame);
        auto data = tempdup(dataIn[0]);
        auto withins = newStack!MeanSD(data.length);
        withins[] = MeanSD.init;
    } 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)) {
    mixin(newFrame);
    size_t N = 0;

    static if(dataIn.length == 1 && isInputRange!(typeof(dataIn[0].front))) {
        auto data = tempdup(dataIn[0]);
        alias ElementType!(typeof(data[0])) C;
        static if(dstats.base.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!(dstats.base.hasLength, T)) {
            enum bool useLength = true;
        } else {
            enum bool useLength = false;
        }
    }

    size_t[] lengths = newStack!size_t(data.length);
    static if(useLength) {
        foreach(i, rng; data) {
            auto rngLen = rng.length;
            lengths[i] = rngLen;
            N += rngLen;
        }
        auto dataArray = newStack!(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 = newStack!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))) {
        mixin(newFrame);
        auto data = tempdup(dataIn[0]);
        auto ranks = newStack!double(data.length);
        auto dataPoints = newStack!double(data.length);
        auto colMeans = newStack!Mean(data.length);
        colMeans[] = Mean.init;
    } 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))) {

    mixin(newFrame);
    alias Unqual!(CommonType!(ElementType!(T), ElementType!(U))) C;

    static if(dstats.base.hasLength!T && dstats.base.hasLength!U) {
        auto n1 = sample1.length, n2 = sample2.length, N = n1 + n2;
        auto combined = newStack!(C)(N);
        rangeCopy(combined[0..n1], sample1);
        rangeCopy(combined[n1..$], sample2);
    } 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 = newStack!(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;

    float* cache = (newStack!(float)((n1 + 1) * (W + 1))).ptr;
    float* cachePrev = (newStack!(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);
    }
    TempAlloc.free;
    TempAlloc.free;
    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;
    }

    mixin(newFrame);

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

        double[] diffRanks = newStack!(double)(before.length);
        byte[] signs = newStack!(byte)(before.length);
        double[] diffs = newStack!(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 = newStack!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, replicate(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);
    }

    float* cache = (newStack!(float)((N + 1) * (W + 1))).ptr;
    float* cachePrev = (newStack!(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);
    }
    TempAlloc.free;
    TempAlloc.free;
    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,

    // Kept for compatibility w/ old style, intentionally not documented, may
    // eventually be removed.
    COUNT = count,
    PROPORTION = 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.
 *
 * This is, for all practical purposes, an inherently non-directional test.
 * Therefore, the one-sided verses two-sided option is not provided.
 *
 * 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))) {
    mixin(newFrame);

    static if(isRandomAccessRange!U && dstats.base.hasLength!U) {
        alias countsIn counts;
    } else {
        auto counts = tempdup(countsIn);
    }

    uint N = sum(counts);

    double[] logPs;
    static if(std.range.hasLength!F) {
        logPs = newStack!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 = newStack!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 in some hypothetical experiment.
 * uint[] drug1 = [1000, 2000, 1500];
 * uint[] drug2 = [1500, 3000, 2300];
 * uint[] placebo = [500, 1100, 750];
 * assert(approxEqual(chiSquareContingency(drug1, drug2, placebo), 0.2397));
 * ---
 *
 * 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.  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));
}

// For converting between base e and base 2 logarithms.
private enum loge2 = 0.69314718055994530941723212145817656807550013436025525412;

// 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) {
    mixin(newFrame);
    static if(isForwardRange!(T[0]) && T.length == 1 &&
        isForwardRange!(typeof(rangesIn[0].front()))) {
        auto ranges = tempdup(rangesIn[0]);
    } else static if(allSatisfy!(isForwardRange, typeof(rangesIn))) {
        alias rangesIn ranges;
    } else {
        static assert(0, "Can only perform contingency table test" ~
            " on a tuple of ranges or a range of ranges.");
    }

    double[] colSums = newStack!(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 = chiSquareCDFR(chiSq, (nRows - 1) * (nCols - 1));

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

    // 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, an
exception 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);
---
*/
TestRes chiSquareObs(T, U)(T x, U y)
if(isInputRange!T && isInputRange!U) {
    uint xFreedom, yFreedom, n;
    typeof(return) ret;

    static if(!dstats.base.hasLength!T && !dstats.base.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));
}

/**
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 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(!dstats.base.hasLength!T && !dstats.base.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 * loge2 * 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) {

    enum needsHeap = dstats.infotheory.NeedsHeap!T ||
        dstats.infotheory.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 {
        mixin(newFrame);
        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));
        auto xCounts = StackHash!(ElementType!T, Uint)(max(10, len / 40));
        auto yCounts = StackHash!(ElementType!U, Uint)(max(10, len / 40));
    }

    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++;
        }

        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--;
        }

        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].dup, [1,2,3,4,5].dup).testStat, 0));
    assert(approxEqual(ksTestDestructive([1,2,3,4,5].dup, [1,2,2,3,5].dup).testStat, -.2));
    assert(approxEqual(ksTest([-1,0,2,8, 6].dup, [1,2,2,3,5].dup).testStat, .4));
    assert(approxEqual(ksTest([1,2,3,4,5].dup, [1,2,2,3,5,7,8].dup).testStat, .2857));
    assert(approxEqual(ksTestDestructive([1, 2, 3, 4, 4, 4, 5].dup,
           [1, 2, 3, 4, 5, 5, 5].dup).testStat, .2857));

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

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

/**One-sample KS test against a reference distribution, doesn't modify input
 * data.  Takes a function pointer or delegate for the CDF of refernce
 * distribution.
 *
 * Returns:  A TestRes with the K-S 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];
 * double 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 TAState = TempAlloc.getState;
    scope(exit) {
        TempAlloc.free(TAState);
        TempAlloc.free(TAState);
    }
    return ksTestDDestructive(tempdup(F), tempdup(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)) {
    scope(exit) TempAlloc.free;
    return ksTestDDestructive(tempdup(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 - kolmDist(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 - kolmDist(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, IterType!(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));
        }
    }
}

// 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);

    PearsonCor pearsonRes = 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:  T must be a numeric type.  The only reason this is a template and
 * not a plain old function is DMD bug 2972.
 */
ConfInt pearsonCorTest(T)(double cor, T N, Alt alt = Alt.twoSided, double confLevel = 0.95)
if(isNumeric!(T)) {
    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(t4.p == 0);
    assert(t5.p == 1);
    assert(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(t7.p == 0);
    assert(t8.p == 0);
    assert(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(!dstats.base.hasLength!T) {
        auto r1 = tempdup(range1);
        scope(exit) TempAlloc.free();
    } else {
        alias range1 r1;
    }
    double N = r1.length;

    return pearsonCorTest(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 some 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.
 *
 * 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)) {
    mixin(newFrame);
    auto i1d = tempdup(range1);
    auto i2d = tempdup(range2);
    immutable res = 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));
    double[] cur = newStack!double(swaps + 1);
    double[] prev = newStack!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);
    }
    TempAlloc.free;
    TempAlloc.free;
    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));

    // 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 = toArray(take(randRange!rNorm(0, 1), nToTake));
        auto rhs = toArray(take(randRange!rNorm(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 {
    /// Assume that dependency among hypotheses may exist.  (More conservative.)
    yes,

    /// Assume hypotheses are independent.  (Less conservative.)
    no,

    // Kept for compatibility with old style, intentionally not documented,
    // may eventually be removed.
    TRUE = yes,
    FALSE = no
}

/**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;
        }
    }

    mixin(newFrame);
    auto perm = newStack!(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.TRUE);
    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.TRUE);
    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));

    mixin(newFrame);
    auto perm = newStack!(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)) {
    mixin(newFrame);

    auto qVals = array(map!(to!float)(pVals));
    auto perm = newStack!(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]);
}


// Verify that there are no TempAlloc memory leaks anywhere in the code covered
// by the unittest.  This should always be the last unittest of the module.
unittest {
    auto TAState = TempAlloc.getState;
    assert(TAState.used == 0);
    assert(TAState.nblocks < 2);
}