Numerical-validation suite

Every statistic ChakataStat prints is checked against an independently computed reference value before a release. The per-core unit tests (rust/src/stats.rs, rust/src/analyze/tests.rs) check a handful of known answers by hand and serve as spot checks. This suite is the complement: it runs the procedure catalog end to end through the shipped native engine and checks, for each procedure, the statistics a reader acts on — the test statistic, the p-value, the parameter estimates, the effect sizes — against values computed independently by an established reference package. Not every printed cell: the shape of every result is guarded separately, by analysis_test.dart. This is production-readiness gate 4 in production_readness.md.

The pages of this suite

This page is the entry point: how the suite works, and what it covers. The detail sits on its own pages so each can be read on its own.

  • NIST StRD certified datasets — the certified-value leg: every dataset NIST publishes, run end to end and scored by digits of agreement.
  • Cross-OS agreement — what the three operating systems' engines compute relative to each other, per pin: digits of agreement, worst pins first, regenerated from each release candidate's tier-2 run. The statement the tolerances cannot make.
  • Statistics-package file formats — how the SPSS, Stata and SAS readers are proven against the packages that wrote the files.
  • Tolerances and recorded conventions — what each tolerance means, and why a number is pinned as tightly as it is.

The recorded conventions, where a statistic has more than one defensible definition and ChakataStat's choice is recorded:

  • The convention inventory — every statistic with more than one defensible definition, and where ChakataStat's choice is written down: recorded, recorded now, single-definition, or pending a named item (0.23.0 item 9).
  • Descriptives, tests and the classical toolbox — the definitions the inventory found implicit, written once: quantiles and moments, continuity corrections, tie handling, effect-size denominators, CI families, post hoc tests, ties in Cox, lags in ADF, factor rotation, clustering defaults, chart geometry.
  • Regression — GLM quantile residuals per family, mediation follow-ups, regression through the origin, least squares by QR.
  • Latent variable — latent class analysis, item response theory, CFA follow-ups.
  • Clinical and survey — complex survey designs, equivalence tests (TOST), Bland–Altman method comparison.
  • Resampling and time series — the LOESS smoother, permutation tests, VAR / Granger causality.

And, for whoever maintains the suite: regenerating the reference values — the generator, the pinned reference environment, and what to do when a reference package changes its answer.

How it works

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

Artifact Produced by Holds
test/validation/reference_values.json tool/validation/generate_references.py for every checked statistic: the reference value, the tol the engine must meet, and the source (the exact reference call)
test/numerical_validation_test.dart hand-written for every id: the engine request and the table cell to read

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.

Dataset Shape Exercises
clinic.csv 60×12: group sex age bp chol bmi s1 s2 s3 outcome visits grp descriptives, crosstabs, the t-test / ANOVA / Levene family, correlation, every regression family, the GLM family, factor/PCA, clustering, discriminant, MANOVA, most nonparametrics
paired.csv 30×3: pre post follow paired t, Wilcoxon, sign, Friedman
surv.csv 80×4: time event arm age Kaplan-Meier (+ log-rank), Cox regression
psm.csv 80×3: treat x1 x2 Propensity score (logistic score model, matching, IPTW, balance)
ts.csv 120×1: value (AR(1)) ACF / PACF / Box-Ljung, ARIMA
grunfeld.csv 200×5: firm year invest value capital (the classic 10-firm Grunfeld panel; statsmodels' copy minus American Steel) Panel regression (FE/RE/pooled, Hausman, clustered SEs)
grunfeld_ub.csv 185×5: the same, three firms' series deterministically shortened Panel regression on an unbalanced panel (harmonic T̄, per-entity θ)

The second oracle — R against Python, pin by pin

Until 0.23.0 every pin in reference_values.json had exactly one oracle: the Python package that produced it. An engine that shared that package's convention — or its mistake — passed without anyone learning anything. The second oracle closes that for the families a regulated user leans on most: tool/validation/generate_second_oracle_references.R recomputes the same statistics in R, from the same committed CSV bytes, with implementations that share no code with the Python ones, and writes test/validation/second_oracle_reference_values.json — the Python file's ids verbatim, the R value and the exact R call. The suite asserts the engine against both at the pin's own tolerance (the second oracle has no tolerance of its own), so the claim for a covered pin is two independent implementations agree, and the engine agrees with both.

267 of the 267 pins in reference_values.json have a second oracle, across 52 of its 52 families (asserted by a guard in numerical_validation_test.dart, like the coverage sentence above). The inventory, by the plan's lanes plus the fifth lane the maintainer added on 2026-08-28 so that nothing stays single-oracle:

Lane Families Pins R implementation
1 — linear models and GLMs reg_linear (16) · nls (5) · iv (6) · reg_logistic (4) · reg_probit (2) · reg_poisson (3) · reg_gamma (6) · reg_negbin (2) · glm_resid (17) · multinomial (3) · ordinal (2) · mediation (7) · mediation_mod (14) 87 base lm / glm / nls (partially linear), ivreg, MASS glm.nb / polr, nnet multinom; the quantile-residual brackets from R's ppois / pnbinom / pgamma on R's fitted means
2 — survival survival_km (2) · survival_cox (2) · ph_test (3) 7 survival survdiff, coxph(ties = "efron"), cox.zph
3 — meta-analysis meta_analysis (15) 15 metafor rma (FE, DL, REML), regtest(model = "lm")
4 — classical tests and descriptives explore (10) · crosstabs (5) · ttest_one (4) · ttest_ind (7) · ttest_paired (4) · anova (5) · levene (2) · means (2) · correlate (6) · partial (1) · glm (4) · manova (3) · rm_anova (3) · mann_whitney (4) · kruskal (2) · normality (3) · wilcoxon (2) · sign (1) · friedman (2) · binomial (1) · runs (2) 73 base t.test / cor.test / chisq.test / wilcox.test / kruskal.test / shapiro.test / ks.test / friedman.test / binom.test / manova / quantile(type = 5), car Anova(type = 3) / leveneTest, psych skew / kurtosi / partial.r
5 — the remainder cfa (27) · panel (27) · propensity (8) · acf pacf adf kpss arima ets (10) · factor discriminant hier reliability (9) · violin (4) 85 lavaan cfa; plm within / random(swar) / pooling / within_intercept / phtest / vcovHC; MatchIt matchit(nearest, linear.logit); base acf / pacf / Box.test / HoltWinters, urca ur.df / punitroot / ur.kpss; base eigen / dist / manova, psych alpha; the Gaussian KDE as a formula at the engine's grid

The rule for a disagreement. Where the two oracles differ by more than the pin's tolerance the R entry carries a divergence naming the convention behind it, the engine's choice stays the primary pin's, and the value is never averaged. The suite refuses an unexplained disagreement and a divergence recorded for a pair that agrees, so the record cannot go stale in either direction. Nine pins diverge today, three conventions:

  • reg_gamma/m2ll — the Gamma log-likelihood is evaluated at a different dispersion: R's logLik.glm uses deviance/n (the family's aic), statsmodels the Pearson χ²/df scale. Coefficients, deviance and the Pearson dispersion agree to 10+ digits. Recorded in conventions/regression.md.
  • ph_test/* (three pins) — two published forms of the Grambsch–Therneau test: survival ≥ 3.0's cox.zph uses the exact score test, lifelines and the engine the 1994 paper's simplified approximation. Recorded in conventions/clinical-and-survey.md.
  • panel/ub_* (five pins) — the unbalanced-panel Swamy–Arora σ²ᵤ: linearmodels and the engine use the harmonic-mean T̄ (Stata's xtreg, re), plm's swar follows Baltagi & Chang (1994)'s trace-based correction; θ, the RE coefficients and SEs and the Hausman statistic all move with it, while every balanced-panel pin agrees to 15 digits. Recorded in conventions/regression.md.

The agreement, where they agree. Closed-form statistics agree to 13–16 digits — including the balanced panel (15.4 median), the propensity block (13.8), factor/discriminant (14.3) and the violin density (14.4). The iteratively-estimated families are where the digits go — cfa 6.2 (two optimizers on F_ML), reg_negbin 6.1 (two estimators of θ), multinomial 6.4, ordinal 5.7, glm_resid 6.4 (brackets from two IRLS fits), survival_cox 7.6, meta_analysis 7.8 (REML), ets 7.1 (a one-parameter optimum), pacf 4.8, nls 3.3 on a flat objective where the pinned SSE and R² are optimizer-invariant but the parameters are not — and adf/p 2.3, two published approximations (MacKinnon 1994 and 1996) of the same Dickey–Fuller distribution, both inside the pin's 1e-3. The generator iterates every R fit to machine precision on purpose (R's default IRLS, REML and Holt–Winters thresholds leave 5–8 digits), and the suite's floor test is where a regeneration at the defaults would show.

The installation check — the suite's shipped subset

Everything above runs on the developer's machine and in CI. Since 0.23.0 (item 1 of that release) a subset of it ships inside the application and runs on the user's installed binary: Help → Verify this Installation in the app, ChakataStat --verify headless (exit 0 only on a full pass), the report rendered as ordinary output tables. The User Guide page Verify this Installation says what it proves and what it does not; this section records how the subset is chosen and how it is kept honest.

The manifest is a strict subset, never a fork. The shipped files — assets/verify/manifest.json and the CSVs under assets/verify/data/ — are generated from this suite by tool/validation/generate_verify_manifest.dart (run under flutter test), whose selection lives in test/validation/verify_manifest_source.dart. Every expected value, tolerance and provenance string is copied from its source row — reference_values.json, survey_reference_values.json, or the certified block of a NIST .dat — at build time, and test/verify_manifest_test.dart fails whenever the committed assets differ from a fresh build, so the manifest cannot be edited by hand and cannot drift when a reference moves. The same guard proves every catalog family is either represented by a check or excused by name with a reason recorded in the manifest — absence is written down, never implied.

The selection (195 checks over 19 datasets):

  • 63 pinned statistics from the package-referenced suites — 60 from the Python file, 3 from the R survey file — chosen so that every family with a pinned reference is represented and the families a regulated user runs weigh most: descriptive statistics (8), comparing means (12), correlation (5), regression and modeling (12: linear, logistic, probit, Poisson, negative binomial, gamma, IV, mediation), nonparametric tests (6), scale reduction and classification (6: reliability, factor, discriminant, propensity, CFA), survival (3), time series (4), meta-analysis (4) and complex samples (3, through a declared design). Each is judged by the same tolerance this suite pins it to.
  • The whole NIST linear least-squares suite — all eleven datasets, every certified parameter estimate and standard deviation plus R² and the residual standard deviation (132 checks), judged by the same LRE floor for the dataset's difficulty grade that nist_strd_test.dart asserts. The data ship as named columns generated from the .dat file, with the polynomial powers synthesized exactly as the suite synthesizes them.
  • Six families excused, each with its reason in the manifest: discovering the dataset and transforms print no statistic; charts return geometry rather than result tables; Bayesian tests, power and sample size and tables have no pinned reference row in this suite, and a subset cannot add one.

The verifier reads the engine's raw numbers off ToolSuccess.result (the seam added for it — printed cells carry three decimals, and digits of agreement cannot be counted from those) and scores each with the same logRelativeError this suite uses for NIST. On the developer's Linux build the shipped manifest passes in full; the worst digits are the negative-binomial coefficient (6.1, cross-implementation agreement with statsmodels, not a certification) and the Wampler5 estimates (6.1–6.5, the exact-fit polynomial family), both well inside their floors.

The adversarial battery — what every procedure does at the edge

Every defect 0.22.0 found was an edge case nobody had asked about. Since 0.23.0 item 6 the answer is asked and pinned for 84 rows of degenerate input, each driven through ToolDispatcher — the errors-as-data seam, so a crash is a failure — by test/adversarial_battery_test.dart over three small committed fixtures in test/fixtures/battery/, whose README.md is the list: each row names the fixture, the call, the outcome class and the reference or reason. A guard requires the README's row ids and the suite's to be the same set.

The rows cover four groups — degenerate data (a constant, n = 1 or 2, one case per group, an all-missing column), degenerate designs (exact collinearity, a one-level factor, complete and quasi-complete separation, weights that are zero, negative or huge, weighted N below the parameter count), numeric extremes (the far tail, 10¹⁵ and 10⁻¹⁵ scales, the limits of a double) and semantics edges (a filter selecting no case, an empty split group) — across every tool family; families with no applicable edge say so in the README rather than being skipped.

Each row lands in one of three classes: a number pinned against R, a closed form or an invariance; a named error; or an undefined cell whose reason is recorded and whose neighbours are still pinned. Ten behaviours were fixed when the battery was first run, because a row found a crash-free but wrong answer — a weighted regression printing a negative F on −2.2 df, a separated logistic fit printing B = 12,911 with blank standard errors, an all-zero Poisson outcome "fitted", a negative survival time dropped silently, a constant reliability item giving a standardized α of 0.95, a count of −0. Every decision a row forced is stated once in the conventions inventory's degenerate inputs section, which closed the inventory's last two waiting rows.

Coverage

267 statistics across 52 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. The website's copy of these numbers, and the NIST total, are guarded the same way by website_claims_test.dart, added 2026-08-23 after an outside reviewer found the homepage still saying 191/47/45 while this page correctly said 267/52/58.) 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) · moderated mediation (statsmodels OLS on the PROCESS model-7 and model-14 shapes: the moderated path coefficients and SEs, the conditional indirect effects at the mean ∓/± 1 SD probes, the index of moderated mediation for both single-path forms, and the conditional-slope SE from the coefficient covariance — bootstrap CIs again 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).