Correctness

A statistics program is a measuring instrument, and an instrument is only as good as its calibration. This chapter explains how Dioptra’s numbers are checked, then gives the full record: every comparison, every NIST data set, and the cases where a difference was investigated.

The short version: all 58 NIST StRD reference data sets used (univariate, linear regression, one-way ANOVA and nonlinear regression) pass, and 7,841 comparisons in 300 oracle cases across 28 families of methods agree with SciPy, statsmodels, lifelines and scikit-learn within the tolerance stored with each case, with no failures.

What “correct” means for a number

A computer stores a number in 64 bits (IEEE 754 double precision), which holds about 15 to 16 significant decimal digits. Every arithmetic step can lose some of them. Most of the time the loss is invisible. Sometimes it is catastrophic: subtract two nearly equal large numbers and the leading digits cancel, leaving mostly rounding noise.

A classic example is the variance. The textbook “one-pass” formula

s² = (Σx² − (Σx)²/n) / (n − 1)

is exact in algebra and dangerous in a computer. Take the values 100000001, 100000002, 100000003. The answer is 1. But Σx² is about 3 × 10¹⁶, beyond the range in which doubles hold every integer, and the subtraction of two such numbers leaves almost nothing of the right digits. A program using that formula can report a variance of 0, or a negative one. The two-pass formula (compute the mean, then sum squared deviations from it) gets 1 exactly. Both formulas are “correct”; only one is correct in a computer.

So the useful question is not “is the formula right?” but “how many digits of the answer are right?”

LRE: counting correct digits

The log relative error (LRE) answers that question with one number. If the certified value is c and the program reports x,

LRE = −log₁₀( |x − c| / |c| )

An LRE of 15 means the first 15 significant digits agree; 10 means 10 digits; 0 means none. If c = 1.2345678901 and x = 1.2345678312, the relative error is about 4.8 × 10⁻⁸ and the LRE is 7.3: seven digits right. (When c is 0, the absolute error is used instead.) McCullough (1998) introduced LRE for exactly this purpose: judging statistical software on problems with known answers. As a rough scale, 4 digits is enough for any figure or table in a paper, 7 is what a sound algorithm reaches on ordinary data, and 15 is all a double can hold.

LRE is capped at the number of digits the reference certifies: 15 for most NIST sets, 11 for nonlinear regression. A cap means “as right as can be checked”.

NIST StRD: certified answers to hard problems

The US National Institute of Standards and Technology publishes the Statistical Reference Datasets (StRD): data sets with answers certified to 15 digits by computations in multiple-precision arithmetic. They are graded by difficulty. The lower sets are ordinary data. The higher ones are built to break naive algorithms: values like 1000000000000.4 that differ only in the last digits (ANOVA SmLs07–09), a tenth-degree polynomial on nearly collinear predictors (Filip), nonlinear models whose parameters are almost redundant (MGH09, Bennett5).

Dioptra’s test suite reads each NIST file, runs the same routine the app uses, and computes the LRE of every certified quantity. The test fails if any falls below a set bar. The tables below show the result for every data set.

Some LREs are lower than 15 for a reason no program can beat. The certified value of a data set with entries like 1000000.1 is computed from the decimal numbers, but a computer cannot hold 1000000.1 in binary; it stores the nearest double, which is slightly different. The exact answer for the numbers the computer actually holds therefore differs from the certified answer in the ninth or tenth digit. To separate the algorithm’s error from this input rounding, the tables have starred columns: the LRE against the exact answer for the doubles, computed in rational arithmetic. There Dioptra reaches 15 digits, or close to it, nearly everywhere.

Oracle fixtures: agreeing with the reference implementations

NIST covers only four families of methods. For everything else (t tests, rank tests, the studentised range, Fisher’s exact test, Kaplan–Meier, ROC, PCA and so on) the reference is the established open-source implementation: SciPy, statsmodels, lifelines and scikit-learn. A Python script in the repository (scripts/oracle/) runs those libraries on fixed inputs and writes the inputs and outputs to JSON fixture files with a tolerance for each case. Dioptra’s tests read the fixtures and compare.

The tolerance is chosen per quantity: tight (10⁻¹² relative) for exact rational results such as exact p values, looser where the reference itself is approximate (SciPy’s randomised Dunnett integration, for example). Where Dioptra and a reference disagree beyond rounding, the disagreement was investigated to its source, often in high-precision arithmetic. Several of the notes below record cases where the reference library, not Dioptra, was the less accurate of the two.

A method with no reference check does not ship. The few methods without an oracle fixture (intraclass correlation, Cohen’s kappa, the bootstrap, permutation tests) are tested against published worked examples and mathematical identities, listed at the end of the record.

Reproducing the record

Everything below is regenerated by

cargo run -p dioptra-stats --example correctness_report

and the same comparisons run in cargo test -p dioptra-stats on every change.

Generated from cargo run -p dioptra-stats --example correctness_report (the same comparisons as cargo test -p dioptra-stats), then annotated. This page covers special functions, distributions, descriptive statistics, t tests, Mann–Whitney, Wilcoxon, normality tests, outliers, ANOVA (one-way, two-way, repeated measures, Kruskal–Wallis, Friedman), multiple comparisons, contingency tables, correlation, linear regression, nonlinear regression (curve fitting), survival analysis, ROC curves, Bland–Altman and PCA. Intraclass correlation, Cohen’s kappa, the bootstrap and permutation tests have no oracle fixture; they are tested against published worked examples and identities (last section).

Every comparison uses the tolerance stored with its fixture case (research/05-nist-and-oracle-notes.md §3); exact discrete p-values (Mann–Whitney, Wilcoxon) are tightened to 1e-12 relative as that policy asks.

Oracle fixtures

Worst error / tolerance: the largest |got - expected| / max(abs, rel |expected|) over the family, with each case’s own tolerance (< 1 passes). Max rel. error is over expected values with |value| > 1e-290, including tiny values where the abs tolerance governs; min digits = -log10(max rel. error).

NIST StRD univariate (LRE, McCullough 1998)

LRE capped at 15 (NIST certifies 15 digits). Columns without a star score against the NIST certified values; starred columns against the exact statistics of the inputs as parsed into doubles (rational arithmetic), which separates the algorithm from input rounding: the NumAcc3/NumAcc4 values (1e6 + 0.1 and 1e7 + 0.1 steps) are not representable in binary, which caps any double-precision SD at 9.5 and 8.3 digits against the certified 0.1. Minimum LRE: 8.25 against certified values, 14.74 against the exact statistics of the double inputs.

NIST StRD linear regression (LRE)

Parameters and SDs: minimum over the coefficients. ANOVA table: minimum over the regression and residual SS and MS and F (Wampler1’s F is infinite and not scored). Min*: minimum LRE against the exact least-squares solution for the inputs as parsed into doubles (rational arithmetic; values below 1e-12 scored by absolute error).

NIST StRD one-way ANOVA (LRE)

Min*: minimum LRE against the exact statistics of the inputs as parsed into doubles.

On the linear sets the refined QR (equilibrated, pivoted Householder plus double-double refinement of the augmented system) returns the correctly rounded least-squares solution of the double inputs everywhere (Min* = 15) except Wampler5’s R^2 = 0.0022, which is formed as a difference (13.8). Against the certified values the remaining gap is input rounding: Wampler2’s y values and Pontius’s x values are decimal fractions that binary cannot hold. Filip reaches 14.0 (McCullough’s “acceptable” bar is 7; normal equations score 0, and plain Householder QR about 7-8). The simple-regression routine reaches 13.9 on Norris and 14.7 on NoInt1/NoInt2 through the same solver (tests/nist_linear.rs::nist_simple_path).

On the ANOVA sets every statistic is correctly rounded for the double inputs (Min* = 15: sums of squares on data shifted by the first observation, two-pass means, compensated sums). Against the certified values, SmLs04-09 are limited by their inputs: 1000000.4 and 1000000000000.4 are not representable in binary. On the parsed doubles the exact SS_between of SmLs05 is 16.080000001825 (certified 16.08, so at most LRE 9.94) and of SmLs07 1.6801563 (certified 1.68, at most 4.03); no double-precision program can do better on these inputs. AtmWtAg (7 constant digits) reaches 10.2 for the same reason. The tests hold the lower sets to 13, the average sets to 10 (9.9 for the input-capped SmLs05/06) and every set to 13 against the exact statistics of its double inputs.

NIST StRD nonlinear regression (LRE, both starting points)

LRE capped at 11 (NIST certifies 11 digits). Parameters and SDs: minimum over the parameters. Status “Ambiguous” is Prism’s flag for a parameter with dependency > 0.9999 (the fit converged; the parameters are nearly redundant). Starts that did not converge: 0. Worst parameter LRE over all 54 fits: 9.33.

Every one of the 27 sets converges from both NIST starting points with plain Levenberg–Marquardt (MINPACK lmder port: scaled trust region, pivoted QR, Moré’s lambda update) followed by Newton polishing. No restarts and no geodesic acceleration are needed; geodesic acceleration (available as FitOptions::geodesic) was tried and is not the default, because from MGH09’s first start it converges to a different stationary point. The asked bars were lower >= 7, average >= 6, higher >= 4; the worst parameter LRE is 9.3 (Lanczos3, whose data carry five significant digits), and 52 of the 54 fits have every parameter at 10.3 or better. MGH09, MGH10, Bennett5, Eckerle4, Rat43 and BoxBOD reach 10.4 to 11 from both starts (MGH09 start 1 takes 395 iterations, Bennett5 start 1 745).

Why the polish matters: LM stops on a relative sum-of-squares test, and near the optimum SS(b + d) - SS(b) is of order d^2, so a test on SS can only resolve the parameters to about sqrt(eps). Without the polish, ENSO, MGH09 and Thurber stop at LRE 6.8-7.4. The Newton step solves (J’J + sum r_i H_i) d = -J’r with J’r accumulated in compensated arithmetic and the second derivatives H_i from central differences of the exact Jacobian; a step is accepted while the equilibrated gradient norm falls. That resolves the parameters to about eps x condition, which is what the table shows.

Lanczos1 is the exception, and the limit is its data, not the algorithm. Its certified residual SS is 1.4307867721e-25: the y values are the model rounded to 13 digits, so the residuals are about 7.7e-14, roughly 100 units in the last place of y. Parsing the y values into doubles moves each by up to 2.2e-16, which changes the minimum SS by 6.5e-4 relative. Evaluated in 60-digit decimal arithmetic at Dioptra’s parameters, the SS is 1.43078687e-25 with the decimal data (certified: 1.43078677e-25) but 1.42986e-25 with the data as parsed into doubles: an LRE of 3.2 is the best any double-precision program can reach on these inputs. Dioptra’s 2.9 adds the rounding of the model evaluation in double. The parameters themselves reach 10.6, and the parameter SDs and residual SD inherit the SS’s 3.2 because they scale with sqrt(SS). tests/nist_nonlinear.rs holds every other set to 9 on all four quantities from both starts, and Lanczos1’s SS, SDs and residual SD to 2.5.

Notes on the fixtures

• Studentized range, far upper tail. The worst relative errors in that row (8e-6) are in the SciPy-generated reference, not in Dioptra. For k = 2 the studentized range has a closed form, P(Q > q) = P(|T_df| > q/sqrt 2). At q = 8, df = 1000 that is 2.0119785627859e-8 (40-digit mpmath and our t distribution agree); the fixture holds 2.0119947042e-8 (SciPy studentized_range.sf, off by 8.0e-6 relative), and Dioptra gives 2.01197856277e-8 (1e-11). The case passes on the fixture’s absolute tolerance (1e-11). tests/distributions.rs::quadrature_identities pins k = 2 to the t distribution at 1e-10 relative, and the ignored test studentized_range_vs_nested_quadrature checks k = 2..20, df = 1..1000 and inf against nested adaptive quadrature (worst: 1e-13 absolute on the cdf, 1e-11 relative on the upper tail).

• Dunnett. The fixture’s own double integral is accurate to about 2e-9 relative (p_greater = 1.3658528033e-5 for the unequal-n case); nested adaptive quadrature gives 1.36585280593e-5, which Dioptra matches to 3e-14. Dunnett with one treatment equals the t test exactly and is tested that way.

• ROUT robust centre. The fixture’s centre comes from SciPy’s minimize_scalar; its merit gradient there is about 1e-8, ours about 1e-13, and the two centres differ by up to 3e-10 relative (tolerance 1e-9). Flagged indices agree everywhere.

• Noncentral t. All six tables pass at 1e-10, including x = -3, df = 5, nc = 8, where the reference is 1.899970783e-20 and SciPy returns 8.0e-20.

• Shapiro–Wilk. W and p agree with SciPy to 1.3e-13 even though the fixture tolerance is 1e-7 / 1e-6 (the AS R94 polynomials are single precision; we use the full-precision normal scores, as R does, and so does SciPy).

• NIST NumAcc3/NumAcc4. The SD LREs of 9.5 and 8.3 against the certified 0.1 are the limit set by the inputs: 1000000.1 and 10000000.1 are not representable in binary, and the exact SD of the parsed doubles is 0.10000000003492 and 0.10000000055879. Dioptra reproduces those exact values to 15 digits (starred columns).

• Spearman, perfect correlation (correlation row). For n5 the ranks agree exactly, so rho = 1 and the t-approximation p is 0. SciPy reports 1.40e-24 because its rho rounds to 0.9999999999999999; Dioptra’s rounds to a slightly different neighbour of 1 and gives a p of the same order. The case passes on the absolute tolerance (1e-10); it is the only reason the correlation row shows a relative error above 1. Every other correlation comparison, including the exact Spearman and Kendall p-values, agrees to 2e-13 relative or better.

• Tukey–Kramer, Dunnett, p-value adjustment. The 2.0e-9 relative error is the fixture’s own Dunnett double integral (see the Dunnett note above); Tukey p-values agree to 8e-12 (the worst is a tail p of 5.3e-5), Tukey critical values and CIs to 1e-13, and all seven p-adjustment methods to 1e-13. The loose library cross-checks stored in the same fixture (statsmodels’ qsturng interpolation, SciPy’s randomized-QMC dunnett p and CI, tolerance 2e-3) are reported on their own “sanity” row so that they do not mask the exact comparisons; SciPy’s QMC Dunnett p for the unequal-n case is 3.8e-5 against the exact 2.73e-5.

• Two-way ANOVA trap. The fixture’s type3_treatment_coding_WRONG_for_reference (type III with dummy coding) is checked to differ from Dioptra’s sum-to-zero type III on the unbalanced case; matching it would be the classic error.

• Fisher r x c. The network algorithm is checked against plain enumeration of every table on 60 pseudo-random tables (p to 1e-10, table counts exactly, contingency::tests::network_matches_enumeration) and matches all five fixture tables, including the transposed 3 x 2 and the 2 x 2 that equals fisher_exact. The conditional MLE odds ratio and its exact CI agree with SciPy to 1e-13.

• Exact discrete p-values. Runs-test p (Swed–Eisenhart), Kendall and Spearman exact p and Fisher p are compared at 1e-12 relative or tighter where the fixture tolerance is looser, as the tolerance policy asks for rational quantities.

• Nonlinear regression. The fixture optima come from SciPy curve_fit and least_squares, polished with central-difference Jacobians; Dioptra uses exact derivatives (hand-written for the built-in library, forward automatic differentiation for parsed equations). Parameters, SEs, CIs, SS, AICc, R^2 and the covariance agree to 4e-8 relative or better, except two covariance entries of the Gaussian case (Amplitude-Mean and Mean-SD, 2e-10 and 1e-11: zero but for rounding, the X grid being symmetric about the mean), which differ by 3e-4 relative and pass on the absolute tolerance; that is SciPy’s finite-difference Jacobian. Every fit/ case is run twice: through the parsed fixture equation from the fixture’s start, and through the built-in Prism equation (4PL, one phase decay, Michaelis–Menten, specific binding with Hill slope, Gaussian) from the Prism-style initial-value rules; both reach the same optimum. The robust cases (fixed-scale Cauchy and soft-L1) agree to 6e-8.

• Cox regression. lifelines stops its Newton iterations early: its score at its own estimate is -1.6e-6, -2.1e-7 and -4.6e-6 for the three cases, and for small_censor_at_end its coefficient is 2.8e-6 relative from the maximum (the fixture asks 1e-6). Dioptra’s coefficient equals the maximum found by Newton iterations in 50-digit decimal arithmetic (-0.84574684893160798) to 1e-15, and its log-likelihood is higher than lifelines’ (-17.930671038121 against -17.930671038126). The test therefore checks the coefficient, SE and log-likelihood against the 50-digit maxima at 1e-11, checks that Dioptra’s likelihood is at least lifelines’, and compares the lifelines values at 1e-5; that last comparison is the 4.7e-6 in the survival row. Every KM, Greenwood, confidence-interval, log-rank, Gehan–Breslow–Wilcoxon, Tarone–Ware and hazard-ratio value agrees to 1e-11.

• ROC. The relative error of 1.0 in the ROC row is a Wilson lower limit at x = 0: the fixture holds 1.39e-17 (rounding in NumPy) where the exact value, and Dioptra’s, is 0; it passes on the absolute tolerance. Everything else agrees to 1e-14.

• Hand-tested (no fixture). ICC: the six Shrout & Fleiss (1979) coefficients of their worked example with R psych::ICC’s F tests and McGraw–Wong CIs, to the printed digits. Cohen’s kappa and its SEs: a hand-computed 2 x 2 table. Bootstrap: the BCa interval of a normal mean matches the t interval, BCa shifts a skewed variance interval the right way, and results are reproducible from the seed. Permutation test: exact enumeration against a hand count, Monte Carlo against the t test. xoshiro256**: the reference outputs from state (1, 2, 3, 4). DeLong’s paired comparison: identities (a monotone transform gives zero difference and zero SE; swapping the tests flips z). Profile-likelihood CIs equal the t CIs on a straight line to 1e-7 and are right-skewed for a rate constant; a global 4PL fit with shared Top and Bottom recovers noiseless parameters to 1e-9; ROUT on a curve finds the two planted outliers and none in clean data.

Performance (release build, one core)

From cargo test -p dioptra-stats --release --test benches -- --ignored --nocapture:

The r x c limits are the FISHER_RXC_WORK_LIMIT (5e7 path-child expansions) and FISHER_RXC_PATH_LIMIT (4e6 stored paths per stage). Tables near independence finish well inside them up to 5 x 6 with N = 200; strongly associated tables of that size do not (R’s FEXACT also needs an enlarged workspace or simulation there), and the chi-square test is the fallback. Prism offers Fisher’s test for 2 x 2 tables only. The Tukey time is almost all the studentized-range quantile for the simultaneous CI.