Numerical-validation suite

Production-readiness gate 4. A statistics package has to be able to defend every number it prints. The per-core unit tests (rust/src/stats.rs, rust/src/analyze/tests.rs) check a handful of known answers by hand; that is good spot-checking, but it only confirms the cases we thought to write down. This suite is the complement: it runs the procedure catalog end to end through the real native engine and checks each printed statistic against a value computed independently by an established reference package. See production_readness.md gate 4.

How it works

Two committed artifacts, joined only by a stable string id:

The Python generator computes the reference numbers with SciPy, statsmodels, scikit-learn, lifelines and pingouin. The Dart test loads the same data into the native engine, runs each procedure, pulls the matching cell out of the result table by label (not by raw index, so a column reorder fails loudly rather than silently reading the wrong number), and asserts it is within tol of the reference. Neither side encodes the other's specifics: the JSON never mentions a table layout, and the test never hard-codes a reference number.

A guard test — "every reference value is checked (no orphans either way)" — asserts the set of ids in the JSON and the set exercised by the probes are equal, so coverage cannot silently drift as either side changes.

Shared, byte-identical inputs

The reference tools and the engine must see exactly the same data, or a disagreement could be an input difference rather than a numeric one. So the four datasets under test/validation/data/ are committed CSVs, and both sides read those same bytes — the generator with csv, the test with a small typed-CSV loader (test/validation/validation_data.dart). The generator creates a dataset deterministically only if its file is absent; if the file exists it is read back verbatim, so a future NumPy RNG change can never drift the committed inputs out from under the committed references.

Coverage

191 statistics across 47 procedure families — essentially the full analytic catalog. (These two numbers are asserted, not transcribed: a guard in numerical_validation_test.dart parses this very sentence and fails if it disagrees with reference_values.json. It is worded this way because the count had already gone stale twice — quote it freely, the suite keeps it honest.) For each procedure the suite pins the statistics that matter (test statistics, p-values, parameter estimates, effect sizes), not every cell (the shape of every result is already guarded by test/analysis_test.dart).

Descriptives/Explore · Frequencies · Crosstabs (Pearson, likelihood-ratio, phi, Cramér's V) · one-sample / independent / paired t (+ Welch, Cohen's d) · one-way ANOVA (+ SS, η²) · Levene · Means · Pearson/Spearman/Kendall/partial correlation · linear / logistic / probit / Poisson / negative-binomial / gamma / multinomial / ordinal regression · nonlinear least squares (scipy curve_fit: the optimizer-invariant residual SS and R², plus well-identified parameters) · instrumental variables / 2SLS (statsmodels IV2SLS: the coefficients, the corrected standard error and the first-stage weak-instrument F) · mediation (statsmodels OLS, cross-checked against pingouin mediation_analysis: the a / b / c′ / c path coefficients and standard errors and the a·b indirect point estimate — the seed-dependent bootstrap CI is not pinned) · meta-analysis (statsmodels combine_effects(method_re="dl") on the dedicated heterogeneous meta.csv: fixed- and random-effects pooled estimates and SEs, Cochran's Q and its p, I², the DerSimonian–Laird τ², and Egger's regression test — intercept, SE, t and p via scipy linregress of the standardized effect on the precision; REML pooling on the committed BCG-trials bcg.csv, whose published metafor results — Viechtbauer 2010, JSS — are asserted at generation time, verifying the transcription and the estimator at once) · CFA (the canonical Holzinger–Swineford three-factor model on the committed hs1939.csv: free loadings and their expected-information SEs, the completely standardized solution, factor correlations, factor/residual variances with SEs, and the full fit block — χ², CFI, TLI, RMSEA with its 90% CI, SRMR — against an independent scipy implementation of the same ML problem, with semopy agreement and the published lavaan tutorial values asserted at generation time) · factorial GLM (Type III, + partial η²) · MANOVA · repeated-measures ANOVA (+ partial η²) · Mann-Whitney (+ rank-biserial r) · Kruskal-Wallis · Wilcoxon (+ rank-biserial r) · sign · Friedman · runs · binomial · Shapiro-Wilk / K-S normality · Cronbach's α · factor / PCA · hierarchical clustering · discriminant · Kaplan-Meier (log-rank) · Cox (+ Grambsch-Therneau proportional-hazards test) · ACF / PACF / Box-Ljung · ARIMA · ADF / KPSS stationarity · simple exponential smoothing.

Deliberately not value-checked (and why):

• k-means and two-step clustering — the partition depends on initialization and iteration order, so cluster labels and centers are implementation-specific; there is no stable cross-package reference. (Hierarchical clustering is checked, via its linkage-independent first-merge distance.)
• Trended / seasonal exponential smoothing (Holt, Holt-Winters) — the smoothing parameters sit on flat objective surfaces, so different optimizers land on different (α, β, γ) with near-identical fit; there is no stable cross-package value. (Simple exponential smoothing is checked — its single parameter has a stable optimum.) The trended/seasonal recursions are instead unit-tested at fixed parameters, where the filter is deterministic.
• Custom tables / multiple-response — these are layout generators over counts already covered by Frequencies/Crosstabs, not new statistics.
• Power & sample-size analysis — it takes no dataset, so it cannot ride the shared-CSV harness that joins the engine to reference_values.json. It is instead validated against the same independent reference (SciPy's noncentral nct / ncf / ncx2 distributions, evaluated at the central critical value): the numeric cores in rust/src/analyze/power.rs reproduce SciPy to ~1e-11 in the Rust unit tests, and test/power_analysis_test.dart re-checks those values end to end through the engine. The underlying noncentral CDFs added for it (noncentral_f_cdf, noncentral_chi_square_cdf in rust/src/stats.rs) carry their own identity cross-checks (χ²(1) against the shifted normal, F(1, ν) against the squared noncentral t).
• ROC curve / AUC — validated against the Mann-Whitney U identity rather than the shared-CSV harness: the area under the ROC curve equals U / (n₊·n₋) — the share of positive/negative score pairs the positive outranks, ties at half. That is an algorithm wholly independent of the engine's sort-and-sweep trapezoid, so it is an exact cross-reference, not a second copy of the same code. rust/src/analyze/tests.rs pins the trapezoid AUC to a by-hand pair count (and to a direct pair-enumeration helper) including the tied-score case, and test/roc_analysis_test.dart re-checks the AUC end to end through the engine. The Hanley-McNeil standard error (the CI) and the null-hypothesis standard error (the chance test) are pinned to their closed-form definitions in the same Rust test.
• Intraclass correlation (ICC) — validated against the published Shrout & Fleiss (1979) worked example (6 targets × 4 judges), whose six ICC values are the canonical reference reproduced by R's irr and pingouin. rust/src/analyze/tests.rs pins all six forms to those values and checks every confidence interval brackets its estimate (and that a frequency weight equals replication); test/icc_kendall_test.dart re-checks them through the engine. It rides this dedicated test rather than the shared-CSV harness because the reference is a fixed published dataset, not a value derived from the committed validation CSVs.
• Kendall's W — pinned to hand-computed values (perfect agreement W = 1, reversed rankings W = 0, a worked tie-corrected case) and to the identity that, without ties, its equivalent χ² equals the already-validated Friedman statistic exactly — an independent cross-check of the same machinery. (Since the 0.8.0 effect-size audit below, the Friedman test itself also reports this W inline, via the exact same rescaling χ²/(N(k−1)) of its own already-validated χ² — the Rust unit tests below therefore validate both call sites at once.)
• The 0.8.0 effect-size audit — a sweep of every test's output against the conventional effect size beside it (t tests → Cohen's d/Hedges' g; ANOVA → η²/partial η²/ω²; chi-square → Cramér's V/φ; nonparametric → rank-biserial/ε²) found and closed the gaps: rank-biserial r on Mann-Whitney and Wilcoxon (validated against pingouin.mwu/pingouin.wilcoxon's independent RBC, riding the existing shared-CSV harness above) and partial η² on factorial GLM and repeated-measures ANOVA (validated against pingouin.anova/pingouin.rm_anova's independent np2, likewise on the shared harness). Three statistics stay off the shared harness because no package in the pinned reference set computes them independently — each instead rides a self-checking identity in rust/src/analyze/nonparametric.rs / rust/src/analyze/glm.rs, re-checked end to end in test/analysis_test.dart and test/tool_dispatcher_test.dart: Kruskal-Wallis epsilon-squared (ε² = H/(N−1)) is a deterministic rescaling of the already-validated H, checked against a hand computation (H = 7.2, N = 9 ⇒ ε² = 0.9); Jonckheere-Terpstra's r (r = Z/√N) likewise rescales the already-validated Z (itself pinned to the exact null-variance enumeration recorded below), checked against a worked perfect-trend example; and MANOVA's partial η² (Pillai V/s; Wilks 1 − Λ^(1/s); Hotelling-Lawley T/(T+s); Roy's θ/(1+θ), the SPSS/GLM convention) is checked against the already-validated Pillai/Wilks/Hotelling values on a two-group worked example, where s = 1 collapses all four to the same number — an internal cross-consistency check the formulas must satisfy. Rank-biserial's sign convention (positive when the first-listed group/measurement ranks higher) is the Kerby (2014) favorable/unfavorable-pairs formula, cross-checked against pingouin's RBC sign on the same (group-order-preserving) data rather than assumed.
• Exact small-sample nonparametrics — the Mann-Whitney, Wilcoxon and Kruskal-Wallis exact significances are pinned to closed-form / hand-counted reference values that an enumerator must reproduce exactly: complete separation gives Mann-Whitney 2/C(m+n, m) and Kruskal-Wallis (g! relabellings)/(N!/Πnᵢ!) (6/90 for three pairs), the all-positive case gives Wilcoxon 2/2ⁿ, and these match R's wilcox.test exact mode. rust/src/analyze/tests.rs also checks the count vector sums to C(m+n, m) and is symmetric, and that ties disable the exact path; test/exact_nonparametric_test.dart re-checks the p's through the engine.
• Nonparametric family completion (Cochran's Q, Jonckheere-Terpstra, Dunn) — pinned to hand-computed / enumeration reference values rather than the shared-CSV harness, because the shared datasets carry no related dichotomous or ordered-trend variables and two of the three have no implementation in the pinned reference set (SciPy, statsmodels, scikit-learn, lifelines, pingouin ship neither Jonckheere-Terpstra nor Dunn). Each rides a self-checking identity: Cochran's Q is pinned to a worked 5×3 table (Q = 14/3) and to the identity that, for two treatments, it equals the uncorrected McNemar χ²; Jonckheere-Terpstra's tie-corrected null variance is checked against the exact null variance obtained by enumerating the equally-likely group assignments (untied {1}|{2}|{3} gives 11/12; a tied {1}|{1,2} gives 0.5), plus a perfect-trend worked value and the sign symmetry of a reversed trend; Dunn's pairwise z is pinned to a hand value (z = −4/√3.5 for the extreme pair of {1,2}|{3,4}|{5,6}) and its zero/symmetry properties. All are in rust/src/analyze/tests.rs and re-checked end to end through the engine in test/tool_dispatcher_test.dart.
• Hosmer-Lemeshow goodness-of-fit — the χ² on the logistic model's deciles of risk is pinned to a hand-computed value (four cases split into two groups give Ĥ = 12/7), with the deciles-of-risk grouping verified to scale a frequency weight exactly like replication (the χ² doubles when every weight doubles) and the degrees of freedom to be g − 2. It rides this dedicated test rather than the shared-CSV harness because the pinned reference set (SciPy, statsmodels, scikit-learn, lifelines, pingouin) has no Hosmer-Lemeshow implementation, and the statistic's group cutpoints are convention-dependent. rust/src/analyze/tests.rs pins the three properties; test/tool_dispatcher_test.dart re-checks it end to end through a fitted logistic model.
• Bootstrap confidence intervals (the bootstrap kind) — a resampling procedure cannot ride the shared-CSV harness: the interval depends on a pseudo-random resample, and no two packages agree on a bootstrap CI unless they share a PRNG and draw order, so there is no cross-package value to pin. It is instead validated against four properties that do hold independent of the randomness, in rust/src/analyze/resample.rs and re-checked end to end in test/bootstrap_test.dart: (1) the observed statistic equals the direct statistic on the full sample, and the interval brackets it; (2) byte-identical reproducibility — the same data, statistic, seed and replicate count give the same interval on every run and platform (the SplitMix64 contract documented in the module), while a different seed gives a different one; (3) the bootstrap standard error of the mean converges to the analytic s/√n at large B; and (4) the BCa interval reduces to the percentile interval for symmetric data (z₀ ≈ 0, acceleration ≈ 0). The pure statistic cores (mean, median, sd, variance, CV, Pearson, Spearman, mean difference) and the type-7 quantile are additionally pinned to hand-computed values, and Spearman is checked against the rank-of-a-monotone-transform identity (ρ = 1).
• Multiple imputation + Rubin pooling (the multiple_imputation kind) — like the bootstrap, a procedure built on pseudo-random draws cannot ride the shared-CSV harness: the chained-equations imputer redraws each missing entry from a Bayesian predictive distribution, and no two packages agree on the imputed values (and so on the pooled estimate) unless they share a PRNG and draw order. It is instead validated against properties that hold independent of the randomness, in rust/src/analyze/multiple_imputation.rs and re-checked end to end in test/multiple_imputation_test.dart: (1) Rubin's rules are exact closed forms — with identical per-imputation estimates the between-variance B = 0, so the pooled SE is exactly √Ū and the fraction of missing information is 0; with varying estimates the total variance is exactly Ū + (1 + 1/m)·B; (2) the no-missing identity — with no missing values the pooled linear regression equals plain OLS (the slope of y = 2x + 3 is recovered exactly and FMI ≈ 0); (3) byte-identical reproducibility — the same data, analysis, m, cycle count and seed give the same pooled result on every run and platform (the SplitMix64 contract documented in the module), while a different seed gives a different one; and (4) the imputer recovers a known linear dependence — missing y at a known x on the line y = 2x imputes to a mean near the true value. The Cholesky factor used for the posterior coefficient draw is pinned to L·Lᵀ = A (and to rejecting a non-positive-definite matrix).
• Dimension reduction — canonical correlation, correspondence analysis, classical MDS — validated against self-checking identities rather than the shared-CSV harness (the committed datasets carry no natural two-set or two-way categorical structure, and the pinned reference set does not ship a directly comparable canonical-correlation / correspondence routine). Each rides an exact identity in rust/src/analyze/reduction.rs, re-checked end to end in test/dimension_reduction_test.dart: the canonical correlation of one variable per set equals the absolute Pearson correlation between them (and Wilks' Λ = 1 − r²); the correspondence total inertia equals χ²/N on a worked 2×2 table (O = [[3,1],[1,3]], χ² = 2, inertia = 0.25); and classical MDS reproduces the exact pairwise Euclidean distances of a planar configuration (a rectangle embedded in two dimensions), with a truly planar set leaving no spurious third dimension.
• Parametric / AFT survival and competing risks — validated against closed-form identities the maximum-likelihood and cumulative-incidence cores must reproduce, in rust/src/analyze/parametric.rs and re-checked end to end in test/parametric_survival_test.dart. For the AFT fits the intercept-only MLE has a known answer: the exponential model's intercept is ln(mean time) and the log-normal model's intercept and scale are the ordinary normal MLE of log t (mean and population SD) — both reproduced to ~1e-3 by the Nelder-Mead fit — and replacing an event with right-censoring at the same time strictly increases the estimated survival time. The competing-risks CIF is pinned to a hand-worked example (times 1–4 with causes 1, 2, 1, censored give final cumulative incidence 0.5 and 0.25, summing to 1 − S(end) = 0.75), with the weight verified to act as case replication. They ride these dedicated identities rather than the shared-CSV harness because the committed datasets carry no censored multi-cause survival structure and the pinned reference set's parametric-survival routines (lifelines) would re-import a near-identical Nelder-Mead fit rather than an independent check.
• Quantile regression, GEE and loglinear (the C4 regression extensions) — validated against self-checking identities in rust/src/analyze/regression_ext.rs, re-checked end to end in test/regression_ext_test.dart. Quantile regression recovers an exact line y = 2x + 1 at every τ (zero residuals ⇒ all quantiles coincide) and an intercept-only fit equals the sample τ-quantile — its point estimates are exact even though the sparsity-based SEs are only asymptotic. GEE with an independence working correlation has point estimates identical to the corresponding GLM (Gaussian ⇒ OLS), pinned by comparing to the engine's own weighted least squares; the exchangeable/AR(1) structures only change the (robust) standard errors. Loglinear independence fitted counts equal the expected counts rowᵢ·colⱼ/n and its deviance equals the table's likelihood-ratio G², with the saturated model fitting perfectly (deviance 0). These ride identities rather than the shared-CSV harness because the committed datasets carry no clustered or multi-way categorical structure, and the identities are exact (not tolerance-limited) cross-checks.
• Bayesian t test (the bayesian_ttest kind, roadmap C4a) — the JZS Bayes factor is validated against a second, independent derivation rather than the shared-CSV harness, in rust/src/analyze/bayesian.rs. The reported BF₁₀ is the Rouder (2009) g-integral (the same integrand pingouin.bayesfactor_ttest evaluates); the test recomputes it a completely different way — as the ratio of marginal likelihoods, integrating the noncentral-t likelihood of the observed t against the Cauchy prior on δ — and requires the two to agree to 1e-3 across a range of (t, n, df). That dual agreement pins the math without a cross-package dependency (and one case is also pinned to the offline reference value BF₁₀(t=2.5, n=20) ≈ 2.7023, the quantity pingouin returns). The supporting pieces carry their own checks: the noncentral-t density integrates to 1 and equals a central difference of the shipped noncentral_t_cdf; the posterior of δ is sign-symmetric, brackets its median, and shrinks the MLE δ̂ = t/√n toward 0; and the evidence scale maps to the Jeffreys / Lee-&-Wagenmakers categories. It rides these identities rather than the shared harness because neither scipy nor pingouin is needed to defend the number — the two internal derivations are an exact cross-reference.
• Bayesian correlation (the bayesian_correlation kind, roadmap C4b) — the Bayes factor for Pearson's ρ is validated against two further independent derivations, in rust/src/analyze/bayesian.rs. The reported BF₁₀ is the Ly et al. (2016) closed form (the quantity pingouin's default bayesfactor_pearson returns), evaluated through an Euler-transformed Gaussian hypergeometric ₂F₁; the tests require it to agree (to 1e-3) with (1) a direct ρ-integral of the exact reduced-likelihood ratio L(ρ)/L(0) against the prior, and (2) the Savage-Dickey density ratio BF₁₀ = 1/[2·p(ρ=0 | data)] read off the normalized posterior — three routes to the same number through different math. The ₂F₁ series itself is pinned against its Euler integral representation (an independent algorithm), and the closed form against the offline value BF₁₀(r=0.5, n=20) ≈ 2.9159. The posterior of ρ is checked for shrinkage toward 0, sign-symmetry under r → −r, and bracketing of its median. It rides these identities rather than the shared harness for the same reason as the t test — the number is fully defensible from the three internal cross-references.
• Bayesian one-way ANOVA (the bayesian_anova kind, roadmap C4c) — the Zellner-Siow g-prior Bayes factor for the group-means model vs. the null is validated against independent derivations in rust/src/analyze/bayesian.rs. The reported BF₁₀ mixes the closed-form fixed-g ratio BF(g) = (1+g)^((N−k)/2)/(1+g(1−R²))^((N−1)/2) against a Cauchy prior on g; the tests pin that closed form against a brute-force design-matrix computation — |I + gX(XʼX)⁻¹Xʼ|^(−1/2)·[SST/zʼ(I+…)⁻¹z]^((N−1)/2) assembled from the actual sum-to-zero contrast matrix — agreeing to ~1e-16, and confirm the design-matrix R² equals the group-means SS_between/SS_total decomposition. The mixed (integrated) Bayes factor is checked across two integration substitutions (a log-g grid and t = g/(1+g)), pinned to the offline value BF₁₀(N=30, k=3, R²=0.30) ≈ 4.1553, and verified to favour the null at R²=0, increase with R² and N, and tend to 1 as the prior scale → 0. It rides these identities rather than the shared harness because R / BayesFactor::anovaBF is not available on CI; for a balanced design the g-prior BF coincides with Rouder et al.'s (2012) default-prior ANOVA up to the prior-scale convention.
• Bayesian linear regression (the bayesian_regression kind, roadmap C4d) — validated by independent property tests in rust/src/analyze/bayesian.rs. The coefficient posterior (the reference-prior NIG conjugate posterior) has means equal to the OLS estimate, pinned against an independent OLS fit; the residual variance posterior is Inverse-Gamma((N−k)/2, SSE/2), whose mean SSE/(N−k−2) and χ²-quantile credible interval are checked (the χ² quantile — by bisection of the shipped chi_square_p — is itself pinned to textbook values). The regression Bayes factor reuses the same g-prior integral as the Bayesian ANOVA (with k = p+1), already validated there, and is pinned to the offline reference BF₁₀(N=40, p=2, R²=0.25, r=√2/4) ≈ 4.4337 and checked to favour the null at R²=0. It rides these identities rather than the shared harness because R / BayesFactor::regressionBF is not available on CI, and the posterior means collapse to OLS by construction (the diffuse-prior limit).

NIST StRD certified datasets

The cross-package suite above proves the engine agrees with other trusted software. The StRD leg proves something stronger for the procedures it covers: agreement with values certified by a national metrology institute. The NIST Statistical Reference Datasets (StRD) were published by the U.S. National Institute of Standards and Technology specifically so statistical software can demonstrate its numerical accuracy: each dataset carries certified values computed at NIST in high-precision arithmetic — to 15 significant digits (11 for the nonlinear suite) — with graded numerical difficulty, including deliberately ill-behaved cases (Longley's extreme collinearity, Filip's 10th-degree polynomial, the NumAcc cancellation series). They are public domain and are the standard benchmark in the software-accuracy literature.

How it works

• The files are committed verbatim under test/validation/data/nist/ — byte for byte as published. Each .dat file carries both the data and the certified values in one ASCII document.
• Both are parsed out of the same file by test/validation/nist_strd.dart, so no certified number is hand-transcribed anywhere: if a check disagrees with NIST it can only be the engine, never a copying error. (The one transcription is the nonlinear model expressions — structure, not numbers — rewritten from each file's Model: block into the NLS tool's grammar.) test/nist_strd_parser_test.dart guards the files' integrity and the parser, engine-free.
• test/nist_strd_test.dart drives every dataset end to end through the real native engine (the same NativeEngine.open() path as the cross-package suite) and scores each statistic by the log relative error — LRE, the number of significant digits on which the engine and NIST agree (LRE = −log₁₀(|engine − certified| / |certified|)), the metric of McCullough, "Assessing the Reliability of Statistical Software", The American Statistician 52/53 (1998/1999). Every statistic must meet a minimum-LRE floor keyed on the suite and on NIST's published difficulty grade; the suite prints the per-dataset achieved LREs, recorded below.

Coverage

All four StRD suites, 45 datasets:

The polynomial linear datasets (Pontius, Filip, Wampler) publish only y and x; the higher powers their certified models include are synthesized as columns by the loader. Nonlinear fits start from NIST's Start 2; Start-1 robustness is an optimizer property, not a numerical-correctness one.

Deliberate exclusions (the same standard as the list above — documented, not silent):

• NoInt1 / NoInt2 (linear suite) — they certify regression through the origin, which the engine does not offer (every linear model fits an intercept). If a no-intercept option ever ships, these two files are the ready-made certification for it.
• Nelson (nonlinear suite) — its certified fit is to log(y), a response transform the NLS tool does not apply.
• Wampler F statistics — the exact-fit Wampler datasets certify a zero residual sum of squares, so no finite F exists to compare; their coefficients, SEs, R² and residual SD are all still checked.

What the suite found (and fixed)

Building the suite surfaced two numerical-hygiene defects and one genuine limitation — exactly the kind of discovery it exists for:

• One-way ANOVA cancellation (fixed, 2026-07-02). The between-groups sum of squares was accumulated on the raw values; on data with many constant leading digits (SmLs07-09 are graded higher difficulty for exactly this) the group spread rounded away before the means were ever differenced, leaving 0.5–2.7 digits of agreement. The fix in rust/src/analyze/compare_means.rs: the sums of squares are translation-invariant, so they are now accumulated about a data point (with Neumaier-compensated summation in moments), lifting SmLs07-09 to 3.9–4.0 — the double-precision ceiling for these files (see the floor note below) — and the average-difficulty tier from ~6.5–8.5 to ~9.9–10.2. Pinned by the one_way_anova_is_translation_invariant Rust test in rust/src/ffi.rs (a 2³⁰ offset of exactly representable data must not change SS/MS/F).
• NLS finite-difference step floor (fixed, 2026-07-02). The Jacobian's step had an absolute floor of 1e-6 — larger than some certified parameters themselves (Hahn1's b7 ≈ 1e-7) — so the optimizer steered on a meaningless derivative and stalled far from the minimum (Hahn1 0.7, Kirby2 4.5). The step in rust/src/analyze/nonlinear.rs is now relative to the parameter's own scale (floor only at exactly zero), lifting Hahn1 to 8.0 and Kirby2 to 7.6 and improving nearly every other fit.
• Filip (known limitation, open). The linear solver cannot fit Filip's 10th-degree polynomial design: its condition number (~1e15) is squared by the normal-equations path, and the returned coefficients are meaningless (LRE 0). This is the classic result this dataset was designed to expose — historically several major packages failed it too; solving it needs an orthogonal-decomposition (QR/SVD) least-squares path. Until then Filip is recorded in the test's knownLimits (still scored, and asserted to stay failing so an eventual solver improvement forces this record to be updated). Follow-up: a QR-based regression_linear core would close this and likely lift the Wampler family too; worth noting that on such designs today the engine also raises no conditioning warning — a gate-1 style surface worth adding alongside.

Achieved LREs

From the watched run of flutter test test/nist_strd_test.dart on 2026-07-03 (Windows x64, ChakataStat engine 0.5.0 — the in-house tabular engine). This run re-recorded the table after the Polars removal, and every achieved LRE came out identical to the 2026-07-02 Polars-era run — including ENSO's build-sensitive worst statistic — confirming the storage swap did not perturb a single scored digit. Min LRE is the worst statistic of that dataset; the cap is 15 (11 nonlinear). The floors sit about a digit under the achieved values because the iterative NLS fits amplify last-bit codegen/libm differences across builds and platforms (an engine rebuild alone was observed to move ENSO's worst statistic between 5.0 and 5.9); the closed-form suites are bit-stable. Two floors are bounded by the problem, not the algorithm: higher-difficulty ANOVA by data representation (nine constant leading digits leave an f64 only ~4 digits of the certified SS — NIST's decimal data simply cannot be represented more finely), and higher-difficulty linear by the exact-fit Wampler datasets, whose certified coefficient SDs are exactly 0 and score by absolute LRE.

For context, McCullough's published assessments graded major commercial packages of the era at similar or lower LREs on the same files (several failed Filip outright, and low single digits on SmLs07-09 were common). A current run always prints this table (the suite's tearDownAll), so the recorded numbers can be refreshed by pasting from any watched run.

Tolerances and recorded conventions

Most statistics agree to 1e-4–1e-6. A looser tolerance always names a reason in the reference's source string. Where a statistic has more than one defensible definition, the reference is computed in the engine's convention and the choice is recorded here — these are agreements about definitions, not discrepancies:

• Quantiles (Explore percentiles) — the engine uses the Hazen definition (plotting position (k − ½)/n); the reference uses numpy.percentile(method="hazen"). (NumPy's default "linear" is a different, equally valid convention.)
• Mann-Whitney — the engine reports the asymptotic test from a tie-corrected z with no continuity correction; the reference derives both z and p from that same z, rather than SciPy's continuity-corrected p.
• Type III SS (factorial GLM) — only well-defined under sum-to-zero contrasts; the reference fits C(group, Sum)*C(sex, Sum) so statsmodels' anova_lm(typ=3) reproduces the engine's main-effect SS in the presence of an interaction.
• ARIMA — the engine fits AR by conditional least squares; the AR(1) reference is the OLS lag-1 regression (the exact CLS estimator), not statsmodels' unconditional ML, which would differ slightly.
• Negative binomial — the dispersion α is ML-estimated, with a looser 5e-2 tolerance on the dispersion alone (the mean-model coefficients match tightly).
• CFA — five recorded conventions, each pinned against the references. (1) Scale setting: each factor's first indicator is the reference, its loading fixed at 1; factor variances are free. (2) χ² = N·F with divisor N (not N−1) — pinned by reproducing lavaan's published Holzinger–Swineford χ² = 85.306 = 301·F. (3) Standard errors use the expected (Fisher) information with divisor N — the lavaan/semopy default; the observed-Hessian alternative disagrees with both references exactly where the model misfits, so it is deliberately not used. (4) The baseline model behind CFI/TLI is the independence model (diagonal Σ at the ML variances, closed form F_b = Σ ln s_jj − ln|S|); a saturated (df = 0) model reports CFI = TLI = 1 and no χ² p-value by convention. (5) The RMSEA CI is 90% by convention (not the analysis CI level), from inverting the noncentral χ². The sample covariance uses divisor N (weighted: Σw), matching the estimates the references produce.
• Meta-analysis — three recorded conventions. (1) The between-study variance is the DerSimonian–Laird moment estimate, truncated at zero (as is I²) — the metafor/textbook convention; statsmodels' combine_effects reports the raw (possibly negative) τ² and pools with it, so meta.csv is built deliberately heterogeneous (the generator asserts raw τ² > 0) to keep the two sides comparable — a homogeneous dataset would make the engine's truncated RE collapse onto FE while statsmodels' would not. (2) Pooled and per-study intervals are large-sample Wald (z-based); the Knapp–Hartung adjustment is out of scope for v1. (3) Egger's test is the classic OLS of the standardized effect on the precision with the intercept tested at k − 2 df (equivalent to WLS of effect on SE under 1/SE² weights). The dataset's Weight Cases setting is meaningless by design for meta-analysis — rows are studies and the inverse-variance weights are the statistic's own — so, like time series, it carries the standard unweighted note and has no weighted reference.
• Case weighting — the engine treats Weight Cases as a frequency weight: a weight-w case behaves exactly like w coincident copies (the procedures weighted on 2026-06-18 — clustering, discriminant, RM-ANOVA/MANOVA, survival — are tested against that weight ≡ replication identity). One sub-case carries a convention: Cox with the Efron ties correction under non-integer weights. The tied-event discount steps over the weighted event total, which reproduces replication exactly for integer weights but is an approximation for fractional ones (the same approach lifelines takes); the default Breslow correction is exact for any weights. Time series (ACF/PACF, ARIMA) is deliberately unweighted — a frequency weight is undefined for a series taken in order — so it has no weighted reference. The REML τ² (0.11.0) is the standard fixed-point iteration on the restricted likelihood — metafor's default estimator — bounded at zero, iteration-capped with a hit cap reported in the output note; Q, its p and I² stay defined from the fixed-effect fit and are estimator-independent, so only τ², the random-effects pooling and the per-study random weights change with the tau_method choice.

Regenerating

The reference toolchain is needed only to author the values — CI just reads the committed JSON (so the suite runs everywhere the Dart suite runs, with no Python dependency on the CI machines). To regenerate after adding a procedure or changing a number:

python3 -m venv /tmp/refenv
/tmp/refenv/bin/pip install scipy statsmodels scikit-learn lifelines pingouin
/tmp/refenv/bin/python tool/validation/generate_references.py
dart run tool/cstool.dart test test/numerical_validation_test.dart

The package versions used for the committed values are recorded in the JSON header (_versions) — at the time of writing SciPy 1.18, statsmodels 0.14, scikit-learn 1.9, lifelines 0.30, pingouin 0.6, NumPy 2.5.

Adding a procedure

1. In the generator, add a add("<family>/<stat>", value, tol, "<source>") line computing the statistic with a reference package; re-run it.
2. In the Dart test, add a Probe (or an entry to an existing one) mapping each new id to the engine request and the result cell that holds it.
3. Run the test. The orphan-guard will tell you immediately if the two sides are not in step.

When the reference disagrees, decide whether it is a real engine bug or a convention difference (a different-but-valid definition). If the latter, compute the reference in the engine's convention and record the choice in the Tolerances and recorded conventions list above — never just loosen the tolerance to make a genuine discrepancy pass.