Changeset 48

Timestamp:

05/07/09 00:46:04 (3 years ago)

Author:

dsimcha

Message:

Overhaul tests. Change kurtosis, skewness to more widely used, possibly less biased definitions. This only matters for small samples.

Files:

• trunk/distrib.d (modified) (5 diffs)
• trunk/gamma.d (modified) (3 diffs)
• trunk/random.d (modified) (15 diffs)
• trunk/summary.d (modified) (4 diffs)
• trunk/tests.d (modified) (36 diffs)

trunk/distrib.d

@@ -1254 +1254 @@
-
-///
-
+pure nothrow 
-    real toSquare = (X - X0) / gamma;
-    return 1.0L / (
@@ -1290 +1290 @@
-
-///
-
+pure nothrow 
-    return X0 + gamma * tan(PI * (p - 0.5L));
-}
@@ -1309 +1309 @@
-
-///
-
+pure nothrow
-    real diff = (X - mu);
-    real sign = (diff > 0) ? 1 : -1;
@@ -1323 +1323 @@
-
-///
-
+pure nothrow 
-    real diff = (mu - X);
-    real sign = (diff > 0) ? 1 : -1;
@@ -1352 +1352 @@
-
-/**Kolmogorov distribution.  Used in Kolmogorov-Smirnov testing.*/
-
+pure nothrow 
-    if(X == 0) {
-        //Handle as a special case.  Otherwise, get NAN b/c of divide by zero.

trunk/gamma.d

@@ -78 +78 @@
- * $(NAN) if sign is indeterminate.
- */
-
+ pure nothrow
-{
-    /* Author: Don Clugston. */
@@ -529 +529 @@
-// Continued fraction expansion #1 for incomplete beta integral
-// Use when x < (a+1)/(a+b+2)
-
+ pure nothrow
-{
-    real xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
@@ -612 +612 @@
-// Continued fraction expansion #2 for incomplete beta integral
-// Use when x > (a+1)/(a+b+2)
-
+ pure nothrow
-{
-    real  xk, pk, pkm1, pkm2, qk, qkm1, qkm2;

trunk/random.d

@@ -210 +210 @@
-    foreach(ref elem; observ)
-    elem = rNorm();
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 0  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 0  Observed:  ", median(observ));
@@ -229 +229 @@
-    foreach(ref elem; observ)
-    elem = rCauchy(2, 5);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: N/A  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 2  Observed:  ", median(observ));
@@ -251 +251 @@
-    foreach(ref elem; observ)
-    elem = rStudentT(5);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 0  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 0  Observed:  ", median(observ));
@@ -271 +271 @@
-    foreach(ref elem; observ)
-    elem = rFisher(5, 7);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: ",  7.0 / 5, "  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: ??  Observed:  ", median(observ));
@@ -290 +290 @@
-    foreach(ref elem; observ)
-    elem = rChiSquare(df);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: ", df, "  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: ", df - (2.0L / 3.0L), "  Observed:  ", median(observ));
@@ -715 +715 @@
-        foreach(ref elem; observ)
-        elem = rHypergeometric(n1, n2, n);
-Pval
+Test
-        writeln("100k samples from hypergeom.(", n1, ", ", n2, ", ", n, "):");
-        writeln("\tMean Expected: ", n * cast(real) n1 / (n1 + n2),
@@ -821 +821 @@
-    foreach(ref elem; observ)
-    elem = rLaplace();
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 0  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 0  Observed:  ", median(observ));
@@ -841 +841 @@
-    foreach(ref elem; observ)
-    elem = rExponential(2.0L);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 0.5  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 0.3465  Observed:  ", median(observ));
@@ -900 +900 @@
-    foreach(ref elem; observ)
-    elem = rGamma(2.0L, 3.0L);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 1.5  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: ??  Observed:  ", median(observ));
@@ -963 +963 @@
-        foreach(ref elem; observ)
-        elem = rBeta(a, b);
-Pval
+Test
-        auto summ = summary(observ);
-
+.p
-        writeln("\tMean Expected: ", a / (a + b), " Observed:  ", summ.mean);
-        writeln("\tMedian Expected: ??  Observed:  ", median(observ));
@@ -988 +988 @@
-    foreach(ref elem; observ)
-    elem = rLogistic(2.0L, 3.0L);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: 2  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: 2  Observed:  ", median(observ));
@@ -1007 +1007 @@
-    foreach(ref elem; observ)
-    elem = rLogNorm(-2.0L, 1.0L);
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: ", exp(-1.5), "  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: ", exp(-2.0L), "  Observed:  ", median(observ));
@@ -1028 +1028 @@
-    foreach(ref elem; observ)
-    elem = rWeibull(2.0L, 3.0L);
-Pval
-
+Test
+.p
-    delete observ;
-}
@@ -1057 +1057 @@
-        elem = rWald(4, 7);
-    }
-Pval
-
+Test
+.p
-    writeln("\tMean Expected: ", 4, "  Observed:  ", mean(observ));
-    writeln("\tMedian Expected: ??  Observed:  ", median(observ));
@@ -1077 +1077 @@
-        elem = rRayleigh(3);
-    }
-Pval
-
-
-
-
+Test
+.p
+
+
+

trunk/summary.d

@@ -497 +497 @@
-    ///
-    real skewness() const {
-var()
+_m2 - _mean * _mean
-        real numerator = _m3 - 3 * _mean * _m2 + 2 * _mean * _mean * _mean;
-        return numerator / pow(var, 1.5L);
@@ -506 +506 @@
-        real mean4 = mean * mean;
-        mean4 *= mean4;
-var()
+_m2 - _mean * _mean
-        return (_m4 - 4 * _mean * _m3 + 6 * _mean * _mean * _m2 - 3 * mean4) /
-               (vari * vari) - 3;
@@ -549 +549 @@
-
-unittest {
-Octave
--.92
-2.0775
-0.79523
+Matlab
+0.25
+1.36
+2.1657
-    writefln("Passed kurtosis unittest.");
-}
@@ -573 +573 @@
-    // Values from Octave.
-    assert(approxEqual(skewness([1,2,3,4,5].dup), 0));
-45618
-76783
+5443
+866
-
-    // Test handling of ranges that are not arrays.
-    string[] stringy = ["3", "1", "4", "1", "5", "9", "2", "6", "5"];
-    auto intified = map!(to!(int, string))(stringy);
-45618
+5443
-    writeln("Passed skewness test.");
-}

trunk/tests.d

@@ -67 +67 @@
-
-import dstats.base, dstats.distrib, dstats.alloc, dstats.summary, dstats.sort,
-
+
+
-
-version(unittest) {
@@ -91 +92 @@
-}
-
+/**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.
+    real testStat;
+
+    /**The P-value against the provided alternative.  This struct can
+     * be implicitly converted to just the P-value via alias this.*/
+    real 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.
+    real lowerBound;
+
+    ///  Upper bound of the confidence interval at the level specified.
+    real upperBound;
+
+    alias testRes this;
+
+    ///
+    string toString() {
+        return text("Test Statistic = ", testRes.testStat, "\nP = ", testRes.p,
+                "\nLower Confidence Bound = ", lowerBound,
+                "\nUpper Confidence Bound = ", upperBound);
+    }
+}
+
-/**One-sample Student's T-test for difference between mean of data and
- * a fixed value.  Alternatives are Alt.LESS, meaning mean(data) < mean,
@@ -96 +139 @@
- * != mean.
- *
-
-
+
+
+
+
-if(realIterable!(T)) {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    writeln("Passed 1-sample studentsTTest test.");
-}
@@ -127 +174 @@
- * mean(sample1) < mean(sample2), Alt.GREATER, meaning mean(sample1) >
- * mean(sample2), and Alt.TWOSIDE, meaning mean(sample1) != mean(sample2).
-
-
+
+
+
+
+
+
-if(realIterable!(T) && realIterable!(U)) {
-    size_t n1, n2;
@@ -143 +194 @@
-    real sx1x2 = sqrt(((n1 - 1) * s1summ.var + (n2 - 1) * s2summ.var) /
-                 (n1 + n2 - 2));
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
-unittest {
-    // Values from R.
-
-
+
+
+
+
+
+
+
+
-           0.1173));
-
+0, 
-           0.8827));
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    writeln("Passed 2-sample studentsTTest test.");
-}
@@ -174 +257 @@
- * meaning mean(sample1) > mean(sample2), and Alt.TWOSIDE, meaning mean(sample1)
- * != mean(sample2).
-
-
+
+
+
+
+
+
-if(realIterable!(T) && realIterable!(U)) {
-    size_t n1, n2;
@@ -190 +277 @@
-    auto v1 = s1summ.var, v2 = s2summ.var;
-    real sx1x2 = sqrt(v1 / n1 + v2 / n2);
-
+
+
-    real numerator = v1 / n1 + v2 / n2;
-    numerator *= numerator;
@@ -198 +286 @@
-    denom2 = denom2 * denom2 / (n2 - 1);
-    real df = numerator / (denom1 + denom2);
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    assert(approxEqual(welchTTest([1,3,5,7,9,11].dup, [2,2,1,3,4].dup), 0.06616));
-Alt.LESS)
-  
-Alt.GREATER)
-  
+0
+Alt.LESS),
+0
+Alt.GREATER),
-    writeln("Passed welchTTest test.");
-}
@@ -228 +341 @@
- * equal to testMean.
- *
-
-
+
+
+
+
+
+
-if(realInput!(T) && realInput!(U)) {
-    OnlineMeanSD msd;
@@ -240 +357 @@
-        len++;
-    }
-
-
+
+
+
+
+
+
+
-    if(alt == Alt.LESS) {
-
+
+
+
+
-    } else if(alt == Alt.GREATER) {
-
-
-
+
+
+
+
-    } else {
-
-
+
+
+
+
+
+
+
-}
-
-unittest {
-    // Values from R.
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-    writeln("Passed pairedTTest unittest.");
-}
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- *
- * Input ranges for this function must define a length.
- *
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-if(isInputRange!(T) && dstats.base.hasLength!(T)) {
-    ulong n1 = sample1.length, n2 = sample2.length, N = n1 + n2;
@@ -310 +613 @@
-}
-
-unittest {
-    assert(wilcoxonRankSum([1, 2, 3, 4, 5].dup, [2, 4, 6, 8, 10].dup) == 5);
-    assert(wilcoxonRankSum([2, 4, 6, 8, 10].dup, [1, 2, 3, 4, 5].dup) == 20);
-    assert(wilcoxonRankSum([3, 7, 21, 5, 9].dup, [2, 4, 6, 8, 10].dup) == 15);
-    writeln("Passed wilcoxonRankSum test.");
-}
-
-/**Computes a P-value for a Wilcoxon rank sum test score against the given
- * alternative. Alt.LESS means that mean rank(sample1) < mean rank(sample2).
- * Alt.GREATER means mean rank(sample1) > mean rank(sample2).  Alt.TWOSIDE means
- * mean rank(sample1) != mean rank(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.
- *
- * Note:  Exact p-value computation is never used when tieSum > 0, i.e. when
- * there were ties in the data, because it is not computationally feasible.
- * In these cases, exactThresh will be ignored.
- *
- * Returns:  The p-value against the given alternative.*/
-real wilcoxonRankSumPval(T)(T w, ulong n1, ulong n2, Alt alt = Alt.TWOSIDE,
-                           real tieSum = 0,  uint exactThresh = 50) {
@@ -361 +639 @@
-    assert(approxEqual(wilcoxonRankSumPval(1500, 50, 50, Alt.LESS), .957903));
-    assert(approxEqual(wilcoxonRankSumPval(8500, 100, 200, Alt.LESS), .01704));
-
+W
-    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6), 0.9273));
-    assert(approxEqual(wilcoxonRankSumPval(w, 5, 6, Alt.GREATER), 0.4636));
@@ -367 +645 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
- * 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.*/
@@ -561 +747 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-if(realInput!(T) && realInput!(U) &&
-   dstats.base.hasLength!(T) && dstats.base.hasLength!(U)) {
@@ -629 +873 @@
-}
-
-unittest {
-    assert(wilcoxonSignedRank([1,2,3,4,5].dup, [2,1,4,5,3].dup) == 7.5);
-    assert(wilcoxonSignedRank([3,1,4,1,5].dup, [2,7,1,8,2].dup) == 6);
-    assert(wilcoxonSignedRank([8,6,7,5,3].dup, [0,9,8,6,7].dup) == 5);
-    writeln("Passed wilcoxonSignedRank unittest.");
-}
-
-/**Computes a P-value for a Wilcoxon signed rank test score against the given
- * alternative. Alt.LESS means that elements of before are typically less
- * than corresponding elements of after.  Alt.GREATER means elements of
- * before are typically greater than corresponding elements of after.
- * Alt.TWOSIDE 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 tieSum > 0, i.e. when
- * there were ties in the data, because it is not computationally feasible.
- * In these cases, exactThresh will be ignored.
- *
- * May give a different answer than calling wilcoxonSignedRank, and then
- * passing the W value into wilcoxonSignedRankPval in the presence of ties
- * or equal pairs.
- *
- * The input ranges for this function must define a length and must be
- * forward ranges.
- *
- * Returns:  The p-value against the given alternative.*/
-real wilcoxonSignedRankPval(T)(T before, T after,
-      Alt alt = Alt.TWOSIDE, uint exactThresh = 50)
-if(realInput!(T) && dstats.base.hasLength!(T) && isForwardRange!(T))
-in {
-    assert(before.length == after.length);
-} body {
-    real tieSum;
-    real W = wilcoxonSignedRank(before, after, &tieSum);
-    ulong N = before.length;
-    while(!before.empty && !after.empty) {
-        if(before.front == after.front)
-            N--;
-        before.popFront;
-        after.popFront;
-
-    }
-    return wilcoxonSignedRankPval(W, N, alt, tieSum, exactThresh);
-}
-
-unittest {
-    // Values from R.
-    alias approxEqual ae;
-    // With ties, normal approx.
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,1,4,5,3].dup), 1));
-    assert(ae(wilcoxonSignedRankPval([3,1,4,1,5].dup, [2,7,1,8,2].dup), 0.7865));
-    assert(ae(wilcoxonSignedRankPval([8,6,7,5,3].dup, [0,9,8,6,7].dup), 0.5879));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.LESS), 0.5562));
-    assert(ae(wilcoxonSignedRankPval([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.LESS), 0.3932));
-    assert(ae(wilcoxonSignedRankPval([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.LESS), 0.2940));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,1,4,5,3].dup, Alt.GREATER), 0.5562));
-    assert(ae(wilcoxonSignedRankPval([3,1,4,1,5].dup, [2,7,1,8,2].dup, Alt.GREATER), 0.706));
-    assert(ae(wilcoxonSignedRankPval([8,6,7,5,3].dup, [0,9,8,6,7].dup, Alt.GREATER), 0.7918));
-
-    // Exact.
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,16,32].dup), 0.625));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.LESS), 0.3125));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,16,32].dup, Alt.GREATER), 0.7812));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup), 0.8125));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.LESS), 0.6875));
-    assert(ae(wilcoxonSignedRankPval([1,2,3,4,5].dup, [2,-4,-8,-16,32].dup, Alt.GREATER), 0.4062));
-    writeln("Passed wilcoxonSignedRankPval unittest.");
-}
-
-/**Computes Wilcoxon signed rank test P-value for
- * a set of observations against another set, using the given alternative.
- * Alt.LESS means that elements of before are typically less than corresponding
- * elements of after.  Alt.GREATER means elements of before are typically
- * greater than corresponding elements of after.  Alt.TWOSIDE means that
- * there is a significant difference in either direction.
- *
- * tieSum is a parameter that is used  internally to adjust for
- * ties in the data.  It is computed by wilcoxonSignedRank().
- * exactThresh is the threshold value of N at which this function
- * switches from exact to approximate computation of the p-value.  Exact
- * computation is very slow for large datasets, and for anything but very small
- * datasets, the asymptotic approximation is good enough for all practical
- * purposes.  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.
- *
- * Note:  Exact p-value computation is never used when tieSum > 0, i.e. when
- * there were ties in the data, because it is not computationally feasible.
- * In these cases, exactThresh will be ignored.
- *
- * Returns:  The p-value against the given alternative.*/
-real wilcoxonSignedRankPval(T)(T W, ulong N, Alt alt = Alt.TWOSIDE,
-     real tieSum = 0, uint exactThresh = 50)
@@ -835 +981 @@
- * Alt.GREATER, meaning elements of before are typically greater than
- * elements of after, and Alt.TWOSIDE, meaning that there is a significant
-
-
-
+
+
+
+
+
+
+
-    uint greater, less;
-    while(!before.empty && !after.empty) {
@@ -848 +998 @@
-        after.popFront;
-    }
+    real propGreater = cast(real) greater / (greater + less);
-    if(alt == Alt.LESS) {
-binomialCDF(greater, less + greater, 0.5
+TestRes(propGreater, binomialCDF(greater, less + greater, 0.5)
-    } else if(alt == Alt.GREATER) {
-binomialCDF(less, less + greater, 0.5
+TestRes(propGreater, binomialCDF(less, less + greater, 0.5)
-    } else if(less > greater) {
-2 * binomialCDF(greater, less + greater, 0.5
+TestRes(propGreater, 2 * binomialCDF(greater, less + greater, 0.5)
-    } else if(greater > less) {
-2 * binomialCDF(less, less + greater, 0.5
-1
+ TestRes(propGreater, 2 * binomialCDF(less, less + greater, 0.5)
+TestRes(propGreater, 1)
-}
-
@@ -867 +1018 @@
-    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));
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-    writeln("Passed signTest unittest.");
+}
+
+///For chiSqrFit, is expected value range counts or proportions?
+enum Expected {
+    ///
+    COUNT,
+
+    ///
+    PROPORTION
-}
-
@@ -906 +1037 @@
- * for testing whether a set of observations fits a discrete distribution.
- *
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
- * and is therefore not considered valid when expected values are below 5.
- *
@@ -922 +1059 @@
- *
- * uint[] observed = [980, 1028, 1001, 964, 1102];
-
-
-
+
+
+
+
- * ---
- */
-
-
+
+
+
+
+
+
+
-    uint len = 0;
-    real chiSq = 0;
+    real multiplier = 1;
+
+    if(countProp == Expected.PROPORTION) {
+        multiplier = 0;
+        foreach(elem; observed) {
+            multiplier += elem;
+        }
+    }
+
-    while(!observed.empty && !expected.empty) {
-
+ * multiplier
-        real diff = cast(real) observed.front - e;
-        chiSq += (diff * diff) / e;
@@ -939 +1091 @@
-        len++;
-    }
-chiSqrCDFR(chiSq, len - 1
+TestRes(chiSq, chiSqrCDFR(chiSq, len - 1)
-}
-
@@ -947 +1099 @@
-    uint[] observed = [980, 1028, 1001, 964, 1102];
-    auto expected = repeat(cast(real) sum(observed) / observed.length);
-
-
+
+
+
+
+
+
+
+
-    writeln("Passed chiSqrFit test.");
-}
@@ -1080 +1238 @@
-}
-
-
-
-/**Performs a Kolmogorov-Smirnov (K-S) 2-sample test and returns
- * the D value.  The K-S test is a non-parametric test for a difference between
- * two empirical distributions or between an empirical distribution and a
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-if(isInputRange!(T) && isInputRange!(U)) {
-    auto TAState = TempAlloc.getState;
@@ -1096 +1337 @@
-        TempAlloc.free(TAState);
-    }
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-    qsort(F);
-    qsort(Fprime);
@@ -1128 +1365 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-if(realInput!(T) && (is(Func == function) || is(Func == delegate))) {
-    scope(exit) TempAlloc.free;
-
-
-
-
-
-
-
-
+
+
+
+
+
-    qsort(Femp);
-    real D = 0;
@@ -1173 +1385 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-in {
-    assert(D >= -1);
@@ -1195 +1393 @@
-}
-
-
-
-
+
-in {
-    assert(D >= -1);
@@ -1205 +1401 @@
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
-/**Convenience function.  Converts a dynamic array to a static one, then
- * calls the overload.*/
-real
+TestRes
-if(isIntegral!(T))
-in {
@@ -1250 +1531 @@
-    newTable[1][0] = contingencyTable[1][0];
-    return fisherExact(newTable, alt);
-}
-
-
-/**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.TWOSIDE, meaning
- * true odds ratio != 1.
- *
- * Accepts a 2x2 contingency table as an array of arrays of uints.
- * For now, only does 2x2 contingency tables.
- * Examples:
- * ---
- * real res = fisherExact([[2u, 7], [8, 2]], Alt.LESS);
- * assert(approxEqual(res, 0.01852));  // Odds ratio is very small in this case.
- * ---
- * */
-real fisherExact(T)(const T[2][2] contingencyTable, Alt alt = Alt.TWOSIDE)
-if(isIntegral!(T)) {
-
-    static real fisherLower(const uint[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 real fisherUpper(const uint[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]);
-    }
-
-
-    if(alt == Alt.LESS)
-        return fisherLower(contingencyTable);
-    else if(alt == Alt.GREATER)
-        return fisherUpper(contingencyTable);
-
-    alias contingencyTable 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(real) (n + 1) * (n1 + 1)) / (n1 + n2 + 2));
-    immutable real pExact = hypergeometricPMF(c[0][0], n1, n2, n);
-    immutable real pMode = hypergeometricPMF(mode, n1, n2, n);
-
-    if(approxEqual(pExact, pMode, 1e-7)) {
-        return 1;
-    } else if(c[0][0] < mode) {
-        immutable real pLower = hypergeometricCDF(c[0][0], n1, n2, n);
-
-        // Special case to prevent binary search from getting stuck.
-        if(hypergeometricPMF(n, n1, n2, n) > pExact) {
-            return pLower;
-        }
-
-        // Binary search for where to begin upper half.
-        uint min = mode, max = n, guess = uint.max;
-        while(min != max) {
-            guess = (max == min + 1 && guess == min) ? max :
-                    (cast(ulong) max + cast(ulong) min) / 2UL;
-
-            immutable real 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 && min == max)
-            guess = min;
-
-        return std.algorithm.min(pLower +
-               hypergeometricCDFR(guess, n1, n2, n), 1.0L);
-    } else {
-        immutable real pUpper = hypergeometricCDFR(c[0][0], n1, n2, n);
-
-        // Special case to prevent binary search from getting stuck.
-        if(hypergeometricPMF(0, n1, n2, n) > pExact) {
-            return pUpper;
-        }
-
-        // Binary search for where to begin lower half.
-        uint min = 0, max = mode, guess = uint.max;
-        while(min != max) {
-            guess = (max == min + 1 && guess == min) ? max :
-                    (cast(ulong) max + cast(ulong) min) / 2UL;
-            real 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 && min == max)
-            guess = min;
-
-        return std.algorithm.min(pUpper +
-               hypergeometricCDF(guess, n1, n2, n), 1.0L);
-    }
-}
-
@@ -1394 +1567 @@
-
-    auto res = fisherExact([[19000u, 80000], [20000, 90000]]);
+    assert(approxEqual(res.testStat, 1.068731));
-    assert(approxEqual(res, 3.319e-9));
-    res = fisherExact([[18000u, 80000], [20000, 90000]]);
@@ -1528 +1702 @@
-}
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    switch(alt) {
-        case Alt.TWOSIDE:
-
+
+
+
+
+
-        case Alt.LESS:
-
+
+
+
+
+
-        case Alt.GREATER:
-
+
+
+
+
+
-        default:
-            assert(0);
-    }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
- *
- * Bugs:  Exact computation not yet implemented.  Uses asymptotic approximation
@@ -1565 +1828 @@
- * large N, but is not perfectly accurate.  Not valid for data with very large
- * amounts of ties.  */
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-    real cc = 2.0L / (N * (N - 1));  // Continuity correction.
-    real sd = sqrt(cast(real) (4 * N + 10) / (9 * N * (N - 1)));
@@ -1595 +1838 @@
-    switch(alt) {
-        case Alt.TWOSIDE:
-
-
+TestRes(tau, 
+)
-        case Alt.LESS:
-normalCDF(tau + cc, 0, sd
+TestRes(tau, normalCDF(tau + cc, 0, sd)
-        case Alt.GREATER:
-normalCDFR(tau - cc, 0, sd
+TestRes(tau, normalCDFR(tau - cc, 0, sd)
-        default:
-            assert(0);
@@ -1607 +1850 @@
-
-unittest {
-
-
-
+
+
-    // therefore quite bad.
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-
@@ -1621 +1875 @@
- * This is the most basic, intuitive version of the false discovery rate
- * statistic, and assumes all hypotheses are independent.
+ *
- * Returns:   An array of Q-values with indices
- * corresponding to the indices of the p-values passed in.*/
@@ -1627 +1882 @@
-    // Not optimized at all because I can't imagine anyone writing code where
-    // FDR calculations are the main bottleneck.
-
-
+
+
+
-    foreach(i, ref elem; perm)
-        elem = i;