Replace stdlib normal CDF dependency, make it more benchstat-like

ts-perf/packages/api/package.json

@@ -28,7 +28,6 @@
         "@esfx/collections-hashmap": "^1.0.2",
         "@esfx/collections-hashset": "^1.0.2",
         "@esfx/equatable": "^1.0.2",
-        "@stdlib/stats-base-dists-normal-cdf": "^0.3.1",
         "@ts-perf/core": "workspace:^",
         "iterable-query": "1.0.0-pre.16",
         "power-options": "^0.3.4",

ts-perf/packages/api/src/normalDistribution.ts

@@ -0,0 +1,105 @@
+// Rational approximations for erf/erfc, split by input range to preserve precision in the tails.
+const ERF_NUMERATOR = [
+    9.604973739870516e0,
+    9.002601972038427e1,
+    2.232005345946843e3,
+    7.003325141128051e3,
+    5.55923013010395e4,
+];
+
+const ERF_DENOMINATOR = [
+    3.356171416475031e1,
+    5.213579497801527e2,
+    4.594323829709801e3,
+    2.26290000613891e4,
+    4.92673942608636e4,
+];
+
+const ERFC_NUMERATOR = [
+    2.461969814735305e-10,
+    5.641895648310689e-1,
+    7.463210564422699e0,
+    4.863719709856814e1,
+    1.965208329560771e2,
+    5.264451949954774e2,
+    9.345285271719576e2,
+    1.027551886895157e3,
+    5.575353353693993e2,
+];
+
+const ERFC_DENOMINATOR = [
+    1.32281951154745e1,
+    86.707214088599,
+    354.93777888782,
+    9.757085017432055e2,
+    1.823909166879098e3,
+    2.24633760818711e3,
+    1.656663091941613e3,
+    5.575353408177277e2,
+];
+
+const ERFC_LARGE_NUMERATOR = [
+    5.641895835477551e-1,
+    1.275366707599781e0,
+    5.019050422511805e0,
+    6.160210979930536e0,
+    7.409742699504489e0,
+    2.978866653721002e0,
+];
+
+const ERFC_LARGE_DENOMINATOR = [
+    2.260528632201173e0,
+    9.39603524938001,
+    12.0489539808097,
+    1.708144507475659e1,
+    9.608968090632859e0,
+    3.369076451000815e0,
+];
+
+export function standardNormalCdf(x: number): number {
+    return complementaryErrorFunction(-x / Math.SQRT2) / 2;
+}
+
+function errorFunction(x: number): number {
+    if (Math.abs(x) > 1) {
+        return 1 - complementaryErrorFunction(x);
+    }
+
+    const xSquared = x * x;
+    return x * evaluatePolynomial(xSquared, ERF_NUMERATOR)
+        / evaluatePolynomialWithImplicitLeadingOne(xSquared, ERF_DENOMINATOR);
+}
+
+function complementaryErrorFunction(x: number): number {
+    const absoluteX = Math.abs(x);
+    if (absoluteX < 1) {
+        return 1 - errorFunction(x);
+    }
+
+    const exponent = -absoluteX * absoluteX;
+    if (exponent < -745) {
+        return x < 0 ? 2 : 0;
+    }
+
+    const numerator = absoluteX < 8 ? ERFC_NUMERATOR : ERFC_LARGE_NUMERATOR;
+    const denominator = absoluteX < 8 ? ERFC_DENOMINATOR : ERFC_LARGE_DENOMINATOR;
+    const result = Math.exp(exponent) * evaluatePolynomial(absoluteX, numerator)
+        / evaluatePolynomialWithImplicitLeadingOne(absoluteX, denominator);
+    return x < 0 ? 2 - result : result;
+}
+
+function evaluatePolynomial(x: number, coefficients: readonly number[]): number {
+    let result = coefficients[0];
+    for (let i = 1; i < coefficients.length; i++) {
+        result = result * x + coefficients[i];
+    }
+    return result;
+}
+
+function evaluatePolynomialWithImplicitLeadingOne(x: number, coefficients: readonly number[]): number {
+    let result = x + coefficients[0];
+    for (let i = 1; i < coefficients.length; i++) {
+        result = result * x + coefficients[i];
+    }
+    return result;
+}

ts-perf/packages/api/src/uDistribution.ts

@@ -0,0 +1,207 @@
+// Adapted from Go's golang.org/x/perf/internal/stats U distribution implementation.
+
+interface UEntry {
+    n1: number;
+    twoU: number;
+    value: number;
+}
+
+function choose(n: number, k: number): number {
+    if (k < 0 || k > n) {
+        return 0;
+    }
+
+    k = Math.min(k, n - k);
+    let result = 1;
+    for (let i = 1; i <= k; i++) {
+        result = (result * (n - k + i)) / i;
+    }
+    return result;
+}
+
+function probabilityMasses(n1: number, n2: number, maxU: number): Float64Array {
+    let n = n1;
+    let m = n2;
+    if (n > m) {
+        [n, m] = [m, n];
+    }
+
+    const memo = Array.from({ length: n + 1 }, () => new Float64Array(maxU + 1));
+    for (let currentM = 0; currentM <= m; currentM++) {
+        memo[0][0] = 1;
+
+        const nLimit = Math.min(n, currentM);
+        for (let currentN = 1; currentN <= nLimit; currentN++) {
+            const left = memo[currentN - 1];
+            const right = currentN <= currentM - 1 ? memo[currentN] : memo[currentM - 1];
+            const output = memo[currentN];
+            const uLimit = Math.min(maxU, currentN * currentM);
+
+            for (let u = uLimit; u >= 0; u--) {
+                const leftProbability = u >= currentM ? currentN * left[u - currentM] : 0;
+                const rightProbability = currentM * right[u];
+                output[u] = (leftProbability + rightProbability) / (currentN + currentM);
+            }
+        }
+    }
+
+    return memo[n];
+}
+
+function untiedCdf(n1: number, n2: number, u: number): number {
+    if (u < 0) {
+        return 0;
+    }
+
+    const maxU = n1 * n2;
+    if (u >= maxU) {
+        return 1;
+    }
+
+    let integerU = Math.floor(u);
+    const flip = integerU >= Math.floor((maxU + 1) / 2);
+    if (flip) {
+        integerU = maxU - integerU - 1;
+    }
+
+    const probabilities = probabilityMasses(n1, n2, integerU);
+    let probability = 0;
+    for (let i = 0; i <= integerU; i++) {
+        probability += probabilities[i];
+    }
+    return flip ? 1 - probability : probability;
+}
+
+function entryKey(n1: number, twoU: number): string {
+    return `${n1},${twoU}`;
+}
+
+function setEntry(entries: Map<string, UEntry>, n1: number, twoU: number, value: number): void {
+    entries.set(entryKey(n1, twoU), { n1, twoU, value });
+}
+
+function sum(values: readonly number[]): number {
+    return values.reduce((total, value) => total + value, 0);
+}
+
+function twoUMin(
+    n1: number,
+    ties: readonly number[],
+    tieCount: number,
+    coefficients: readonly number[],
+): number {
+    let twoU = -n1 * n1;
+    let remaining = n1;
+    for (let k = 1; k <= tieCount; k++) {
+        const count = Math.min(remaining, ties[k - 1]);
+        twoU += count * coefficients[k];
+        remaining -= count;
+    }
+    return twoU;
+}
+
+function twoUMax(
+    n1: number,
+    ties: readonly number[],
+    tieCount: number,
+    coefficients: readonly number[],
+): number {
+    let twoU = -n1 * n1;
+    let remaining = n1;
+    for (let k = tieCount; k > 0; k--) {
+        const count = Math.min(remaining, ties[k - 1]);
+        twoU += count * coefficients[k];
+        remaining -= count;
+    }
+    return twoU;
+}
+
+function tiedCumulativeCount(twoU: number, n1: number, ties: readonly number[]): number {
+    const rankCount = ties.length;
+    const coefficients = new Array<number>(rankCount + 1).fill(0);
+    coefficients[1] = ties[0];
+    for (let k = 2; k <= rankCount; k++) {
+        coefficients[k] = coefficients[k - 1] + ties[k - 2] + ties[k - 1];
+    }
+
+    const memo = Array.from({ length: rankCount + 1 }, () => new Map<string, UEntry>());
+    setEntry(memo[rankCount], n1, twoU, 0);
+
+    let tieSum = sum(ties);
+    for (let k = rankCount - 1; k >= 2; k--) {
+        tieSum -= ties[k];
+        for (const next of memo[k + 1].values()) {
+            const rankLow = Math.max(0, next.n1 - tieSum);
+            const rankHigh = Math.min(next.n1, ties[k]);
+            for (let rank = rankLow; rank <= rankHigh; rank++) {
+                const previousTwoU = next.twoU - rank * (coefficients[k + 1] - 2 * next.n1 + rank);
+                const previousN1 = next.n1 - rank;
+                if (
+                    twoUMin(previousN1, ties, k, coefficients) <= previousTwoU
+                    && previousTwoU <= twoUMax(previousN1, ties, k, coefficients)
+                ) {
+                    setEntry(memo[k], previousN1, previousTwoU, 0);
+                }
+            }
+        }
+    }
+
+    const firstTwoRanks = ties[0] + ties[1];
+    for (const entry of memo[2].values()) {
+        let cumulativeCount = 0;
+        const rankLow = Math.max(0, entry.n1 - ties[0]);
+        const rankHigh = Math.trunc((entry.twoU - entry.n1 * (ties[0] - entry.n1)) / firstTwoRanks);
+        for (let rank = rankLow; rank <= rankHigh; rank++) {
+            cumulativeCount += choose(ties[0], entry.n1 - rank) * choose(ties[1], rank);
+        }
+        entry.value = cumulativeCount;
+    }
+
+    tieSum = ties[0];
+    for (let k = 3; k <= rankCount; k++) {
+        tieSum += ties[k - 2];
+        for (const entry of memo[k].values()) {
+            let cumulativeCount = 0;
+            const rankLow = Math.max(0, entry.n1 - tieSum);
+            const rankHigh = Math.min(entry.n1, ties[k - 1]);
+            for (let rank = rankLow; rank <= rankHigh; rank++) {
+                const previousTwoU = entry.twoU - rank * (coefficients[k] - 2 * entry.n1 + rank);
+                const previousN1 = entry.n1 - rank;
+                const previous = memo[k - 1].get(entryKey(previousN1, previousTwoU));
+                let count = previous?.value ?? 0;
+                if (!previous && twoUMax(previousN1, ties, k - 1, coefficients) < previousTwoU) {
+                    count = choose(tieSum, previousN1);
+                }
+                cumulativeCount += count * choose(ties[k - 1], rank);
+            }
+            entry.value = cumulativeCount;
+        }
+    }
+
+    const result = memo[rankCount].get(entryKey(n1, twoU));
+    if (!result) {
+        throw new Error("Failed to compute the exact Mann-Whitney U distribution.");
+    }
+    return result.value;
+}
+
+function tiedCdf(n1: number, n2: number, ties: readonly number[], u: number): number {
+    if (u < 0) {
+        return 0;
+    }
+    if (u >= n1 * n2) {
+        return 1;
+    }
+
+    const twoU = Math.trunc(2 * u);
+    return tiedCumulativeCount(twoU, n1, ties) / choose(n1 + n2, n1);
+}
+
+export function uDistributionCdf(
+    n1: number,
+    n2: number,
+    ties: readonly number[],
+    u: number,
+): number {
+    return ties.some(tie => tie > 1) ? tiedCdf(n1, n2, ties, u) : untiedCdf(n1, n2, u);
+}