NIST StRD certified datasets

Part of the numerical-validation suite. The certified-value leg: every dataset NIST publishes, scored by digits of agreement.

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, all 58 of NIST's datasets (the per-dataset roster is below, guarded):

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 — and so is Nelson's response: its certified model is log[y] = b1 − b2·x1·exp(−b3·x2), so the loader's logResponse flag fits the NLS tool to a synthesized ln(y) column (the decision block below). Nonlinear fits start from NIST's Start 2; Start-1 robustness is an optimizer property, not a numerical-correctness one.

Coverage roster (all 58)

Every StRD dataset NIST publishes, suite by suite — covered ✅ with the engine path, or ❌ with the reason — so "all 58" can be checked at a glance rather than by counting files (and a future exclusion would have to say why, here). Guarded: an engine-free test in nist_strd_parser_test.dart parses this table and asserts that its ✅ set equals the committed files under test/validation/data/nist/ and the per-suite lists in nist_strd_test.dart, that each ✅ row's suite and difficulty match the file's header, and that every ❌ row carries a reason — so the roster cannot drift from the suite. (Added 2026-08-21, 0.21.0 item 11, with the ten nonlinear datasets that had been left unselected when the suite was built.)

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

• 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.

(Two datasets were excluded here until 2026-08-21. NoInt1 / NoInt2 certify regression through the origin, which the engine did not offer; the intercept option shipped as 0.21.0 item 12 and both now run under intercept: false — the conventions are in Regression through the origin. Nelson certifies a fit to log(y), a response transform the NLS tool does not apply; the decision, 0.21.0 item 13: harness-only — the StRD loader's logResponse flag synthesizes ln(y) exactly as it synthesizes the polynomial powers, and the dataset joined the suite under the average-difficulty floor (LRE 7.3). A product-side response_transform argument on nonlinear_regression is out of scope, not deferred — it would duplicate Compute Variable (LN(y)) for one case; the gap was in the harness, not the estimator.)

What the suite found (and fixed)

Building the suite surfaced two numerical-hygiene defects and one genuine limitation (since fixed), and widening it to the full nonlinear suite (2026-08-21) a third defect and a second 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 2026-07-02 → fixed 2026-08-21). The linear solver could not fit Filip's 10th-degree polynomial design: the condition number of the design is squared by the normal-equations path, and the returned coefficients were meaningless (LRE 0) — the classic result this dataset was designed to expose; historically several major packages failed it too. The fix (0.21.0 item 14): regression_linear now solves the weighted least-squares problem by Householder QR of the √W-scaled, column-equilibrated design (wls_qr in rust/src/analyze/regression.rs) — β by back-substitution on R, the covariance bread (XᵀWX)⁻¹ = R⁻¹R⁻ᵀ feeding the model-based SEs, the robust sandwich and the leverage exactly as before, so every downstream statistic moved to the new path at once (the plain fit, its through-origin variant, the stepwise oracle and the survey-design fit all go through it). Filip scores 7.2 (R's lm reaches the same class on this file), Longley moved 7.4 → 11.3, and the exact-fit Wampler family 4.8–8.6 → 6.1–12.7; the higher-difficulty linear floor rose from 4.5 to 5.0. Nothing else moved beyond last-digit noise (Norris 13.3 → 12.2 and Pontius 10.8 → 12.2 are rounding-level at 12+ digits; the cross-package reg_linear/* pins are unchanged). Pinned on the Rust side by a test that reads the committed Filip.dat verbatim (data and certified values) and asserts LRE ≥ 6 for every coefficient and SD, plus the solver's unit tests (closed-form weighted fit and bread, exact dependence refused, a degree-9 design the normal equations get wrong by a factor ~100). The conditioning note the record asked for came along — see Linear least squares by QR.
• NLS iteration cap (fixed, 2026-08-21). Adding the ten unselected nonlinear datasets (0.21.0 item 11) found MGH10 and Bennett5 (both higher difficulty) stalling unconverged at the Levenberg–Marquardt cap of 200 iterations — MGH10 at SSE ≈ 8600 against the certified 87.9 (LRE 0), Bennett5 at LRE 1.6. From NIST's Start 2 they need ~440 and ~450 iterations (Bennett5 converges from its distant Start 1 too, in ~490; MGH10 from Start 1 does not in 5000 and is not a target — it is the textbook case for it). The cap in rust/src/analyze/nonlinear.rs is now 1000, which lifts MGH10 to 8.5 and Bennett5 to 7.0; pinned by the mgh10_converges_from_start_2_within_the_iteration_cap Rust test (the 16 points embedded, certified values to 1e-7 relative). In the same pass the convergence flag was made honest: a fit that stalls because no damped step can decrease the SSE any further is at its floating-point floor (Lanczos2 reaches its certified minimum exactly this way, after 198 iterations) and now reports Converged: Yes; No means one thing — the cap was reached — and the Model Summary note says so.
• Lanczos1 (known limitation, bounded by the problem). An exact-data fit: the certified residual sum of squares is 1.4 × 10⁻²⁵ on data of magnitude 1, i.e. residuals of ~8 × 10⁻¹⁴ against a model that double precision evaluates to ~5 × 10⁻¹⁶ — each residual is good to two or three digits, so the SSE (LRE 2.9) and the parameter standard deviations (∝ √SSE, LRE 3.2) cannot reach the average-difficulty floor in double precision, while the estimates themselves sit at the 11-digit cap (10.6–11.0). Recorded in the test's knownLimits like Filip (still scored, asserted to stay below the floor); the same row is the reason McCullough's assessments treat Lanczos1's SDs separately.

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. One floor is 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). The higher-difficulty linear floor was 4.5 while the normal-equations solver held the exact-fit Wampler datasets (certified SDs of exactly 0, scored by absolute LRE) near it; the QR solver of 2026-08-21 (0.21.0 item 14) lifted the family to 6.1–12.7 and Filip to 7.2, and the floor now sits at 5.0, a digit under the new worst. The linear rows below are from that run. The NoInt1 / NoInt2 rows (regression through the origin, intercept: false) and the eleven nonlinear rows after MGH09 (Gauss2 … Nelson) were added from the watched runs of 2026-08-21 (Linux x64, the 0.21.0 tree — items 12, 11 and 13); the NoInt certified values are exact-arithmetic fits of generated integer data, and the engine reproduces them to the 15-digit cap less a bit.

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.