Numerical and reproducibility repair

A successful run establishes that one procedure reached an output. It does not yet establish that the implementation is correct, the approximation is converged, correlated data have been analyzed properly, another person can repeat the result, or the underlying scientific model describes the intended regime.

This focused review builds the evidence needed between “it ran” and “the claim is supported.” It uses one deterministic derivative calculation and one correlated-mean calculation, then shows how tests, provenance, and independent checks limit the conclusion you can responsibly draw.

Helpful background. You should be able to distinguish an estimator’s bias from its variance and recognize that correlated observations contain less independent information than equally many independent observations. If that is the main obstacle, begin with Statistical ensembles and probability. No particular programming language is assumed; a calculator, spreadsheet, notebook, or short program is enough.

Separate the sources of error

The word “error” hides mechanisms that respond to different checks. Keep at least these five categories distinct:

Source	Typical cause	Evidence that addresses it
Truncation or discretization	Replacing a continuum or infinite process by a finite step, grid, basis, order, or volume	Refinement, observed order, extrapolation, and comparison with a bounded reference case
Roundoff	Finite precision, cancellation, overflow, underflow, or poor conditioning	Precision variation, reformulation, conditioning analysis, and scale-aware arithmetic checks
Statistical uncertainty	Finite random data, including autocorrelation and shared inputs	Covariance-aware standard errors, effective sample size, blocking or resampling, and repeated streams
Model or systematic error	Incorrect or incomplete physical model, regulator, calibration, fit family, boundary condition, or asymptotic assumption	Alternative models, external calibration, sensitivity studies, held-out observables, and theory arguments
Implementation defect	Wrong formula, sign, index, units, data handling, or program logic	Tests that can fail, analytic fixtures, invariants, and independent implementations

An implementation defect is not an uncertainty component to hide inside a wider error bar. Statistical uncertainty usually shrinks with more genuinely independent information; model error generally does not. Truncation can shrink under refinement until roundoff, instability, or computational cost takes over. When several uncertainties are combined, the covariance and the meaning of each component must be stated. Adding unlike quantities in quadrature is not justified merely because each has a percentage attached; the deeper comparison is developed in Complete lattice error budgets.

A convergent calculation with precision turnover

Approximate $f'(1)$ for $f(x)=e^x$ by the central difference

$$D(h) = \frac{e^{1+h}-e^{1-h}}{2h}.$$

Taylor expansion gives

$$D(h) = e+\frac e6h^2+\frac e{120}h^4+\mathcal O(h^6).$$

In exact arithmetic, the leading truncation error is therefore second order. If $E(h)=\lvert D(h)-e\rvert$ and the step is reduced by a factor $r>1$, the observed order is

$$p(h) = \frac{\log\!\left(E(h)/E(h/r)\right)}{\log r}.$$

It should approach $2$ only in the range where the $h^2$ term dominates. Here is one typical IEEE binary64 evaluation; the last digits can vary with the math library and evaluation path.

$h$	$D(h)$	$\lvert D(h)-e\rvert$
$10^{-1}$	$2.72281456394742$	$4.53\times10^{-3}$
$10^{-2}$	$2.71832713338271$	$4.53\times10^{-5}$
$10^{-3}$	$2.71828228150572$	$4.53\times10^{-7}$
$10^{-4}$	$2.71828183298961$	$4.53\times10^{-9}$
$10^{-5}$	$2.71828182851763$	$5.86\times10^{-11}$
$10^{-6}$	$2.71828182829559$	$1.63\times10^{-10}$
$10^{-8}$	$2.71828182185629$	$6.60\times10^{-9}$
$10^{-10}$	$2.71828115572248$	$6.73\times10^{-7}$

The first four refinements reduce the error by roughly $10^2$ per decade, as predicted. Near $h=10^{-5}$ the trend turns over. The numerator subtracts two nearby floating-point numbers, and division by $h$ amplifies their rounding error. A useful error model is

$$E(h) \lesssim C_{\mathrm t}h^2 + C_{\mathrm r}\frac{u}{h},$$

where $u$ is the unit roundoff and the constants depend on the function and implementation. Balancing these terms predicts an optimal scale of order $u^{1/3}$, not the smallest representable positive $h$. The exact constants and conditioning still have to be checked Higham 2002, chs. 1–2.

A precision repeat distinguishes the two mechanisms. In exact arithmetic the leading error at $h=10^{-10}$ is about $4.53\times10^{-21}$. If an arbitrary-precision evaluation with enough digits approaches that value while binary64 gives an error near $10^{-7}$, the turnover is numerical. Increasing precision should move the turnover to smaller $h$. If it does not, investigate the formula, inputs, conditioning, or reference rather than blaming ordinary roundoff.

Tests must be capable of failing

The derivative routine admits several complementary tests:

• Exact structural cases. A constant should give zero and a linear function should give its slope, within a tolerance justified by precision.
• Analytic fixtures. Compare several smooth functions and evaluation points with known derivatives, not only $e^x$ at one symmetric point.
• Convergence behavior. Verify the predicted order over a declared range; one close answer at one $h$ is not a convergence study.
• Metamorphic or invariant checks. The central formula obeys $D(-h)=D(h)$ when evaluated consistently. A broken symmetry can reveal an argument or sign defect even when no reference value is available.
• Deliberate defect. Introduce a plausible wrong denominator, shifted index, unit conversion, or sign and confirm that the intended test fails.
• Independent method. Compare with automatic differentiation, a complex step where applicable, symbolic differentiation, or a separately written finite-difference implementation. Record which assumptions the methods still share.

A tolerance belongs to the test design. Derive it from the analytic error model, precision, conditioning, sampling uncertainty, or a declared external reference before inspecting the result you hope to accept. A test suite can show sensitivity to specified defects; no finite suite proves that every input, regime, or physical model is correct.

Correlated data change the uncertainty, not the estimand

Let $X_1,\ldots,X_N$ be a stationary sequence with

$$\mathbb E[X_i]=\mu, \qquad \operatorname{Cov}(X_i,X_j) = \sigma^2\rho^{\lvert i-j\rvert}, \qquad 0\leq\rho<1,$$

and estimate the stationary mean by

$$\overline X=\frac1N\sum_{i=1}^N X_i.$$

Linearity gives $\mathbb E[\overline X]=\mu$, so the common-mean estimator is unbiased under the stated assumptions. Correlation changes its variance:

$$\operatorname{Var}(\overline X) = \frac{\sigma^2}{N^2} \left[ N+2\sum_{k=1}^{N-1}(N-k)\rho^k \right].$$

For large $N$ at fixed $\rho$,

$$\operatorname{Var}(\overline X) \sim \frac{\sigma^2}{N} \frac{1+\rho}{1-\rho}, \qquad N_{\mathrm{eff}} \sim N\frac{1-\rho}{1+\rho}.$$

Take a concrete marginal variance $\sigma^2=1$, $N=1000$, and $\rho=0.8$. The independent-sample formula gives

$$\operatorname{SE}_{\mathrm{naive}} = \frac1{\sqrt{1000}} \approx0.0316.$$

The exact finite-$N$ covariance sum instead gives

$$\operatorname{Var}(\overline X)=0.00896, \qquad \operatorname{SE}(\overline X)\approx0.0947, \qquad N_{\mathrm{eff}}\approx111.6.$$

The naive standard error is too small by nearly a factor of three. If this standard error is turned into a confidence or credible interval, the sampling distribution or posterior model and the intended coverage or probability interpretation must also be stated. For a sequence generated as a stationary AR(1) process $X_t=0.8X_{t-1}+\varepsilon_t$, unit marginal variance requires $\operatorname{Var}(\varepsilon_t)=1-0.8^2=0.36$; setting the innovation variance to one defines a different example.

In real data, $\rho$ and the covariance tail are unknown. Estimate them with a stated window, batch means, blocking, or another justified time-series method, and test stability as that choice changes. A stable sample mean alone does not establish stationarity or equilibration. Nor does the standard error include model bias, initialization bias, or uncertainty in a physical calibration. The correlation-aware variance framework is developed in Sokal 1997, pp. 131–192 and the practical window tradeoff in Wolff 2004, pp. 143–153. For deeper study, continue to Autocorrelation times and effective sample size.

A compact reproducibility packet

Reproducibility is not a property of the final number alone. It is a relation between a specified target, procedure, inputs, environment, outputs, and acceptance rule. Terminology differs across fields, so define the independence you actually tested:

Check	What is held fixed	What it can reveal
Rerun	Same method, implementation, inputs, and intended environment	Missing instructions, unstable dependencies, nondeterminism, or failure to execute again
Reanalysis	Same raw data or stored outputs, but a new analysis implementation or defensible analysis choices	Analysis-code defects and sensitivity to downstream choices; it does not test data generation
Independent reproduction	A separately written implementation and, when relevant, independently regenerated data from the mathematical and scientific specification	Shared-code defects and some generation failures; it can still share a wrong model or reference formula

These definitions are operational; report overlaps in code, libraries, data, calibration, formulas, and analyst choices. Broader terminology and its limits are discussed by National Academies 2019, ch. 3.

For a small result, the following plain-text packet is enough to expose most missing information:

question, estimand, and scientific regime:
formula, model, or source specification:
exact inputs, units, preprocessing, and checksums:
source revision and dependency versions:
platform, numeric type, and precision:
run command or ordered interaction steps:
tolerances, stopping rules, and failure behavior:
random generator, stream construction, and seed policy:
raw outputs, logs, and checksums:
analysis choices and changes made after seeing results:
expected tests, convergence model, and acceptance rule:
claim supported by this packet—and claims not tested:

A seed without the generator and stream construction is incomplete. A dependency lock without the inputs is incomplete. A notebook that depends on hidden execution order is incomplete. A checksum identifies bytes but does not explain what they mean. Test the packet by using it in a clean environment or handing it to another analyst, then add every missing step that had to be guessed. For a fuller treatment of independence, see Blinding, analysis choices, and independent reproduction.

Keep the scientific claim within the evidence

Different statements require different evidence:

Statement	Minimum supporting evidence	What remains outside it
The procedure executed	Recorded environment, command, inputs, and output	Correctness and scientific meaning
The implementation passed specified tests	Expected results, justified tolerances, and tests shown to detect intended defects	Untested inputs and shared defects
The numerical target is verified in a tested regime	Reference fixture, convergence and precision studies, and independent checks	Extrapolation beyond the tested regime
The estimator has the quoted uncertainty	Estimand, dependence model, covariance-aware analysis, and interval assumptions	Model bias and unmeasured systematics
The scientific conclusion is supported	All preceding evidence plus a justified physical model, observable, regulator or continuum argument, calibration, and systematic-error analysis	Claims outside the declared theory and regime

Agreement between two implementations can verify a bounded numerical target while both implement the same wrong physical model. More samples can narrow a Monte Carlo standard error while leaving finite-volume, cutoff, fit-model, or calibration errors unchanged. A sensitivity envelope is not automatically a probability distribution, and an independent reproduction is a test, not a new uncertainty term to add mechanically.

Exercises

1. Locate convergence and choose a robust test point

Using the derivative table, compute the observed order for the refinements $10^{-1}\to10^{-2}$ and $10^{-2}\to10^{-3}$. What happens for $10^{-5}\to10^{-6}$? Choose a robust binary64 step and justify a test tolerance.

Solution

For decade refinement, $r=10$. The first two estimates are

$$\begin{aligned} p(10^{-1}) &=\log_{10}\!\left( \frac{4.53\times10^{-3}}{4.53\times10^{-5}} \right) \approx2,\\ p(10^{-2}) &=\log_{10}\!\left( \frac{4.53\times10^{-5}}{4.53\times10^{-7}} \right) \approx2. \end{aligned}$$

For $10^{-5}\to10^{-6}$, the error increases from about $5.86\times10^{-11}$ to $1.63\times10^{-10}$, so the observed order is negative. The calculation has left the truncation-dominated range.

Although $h=10^{-5}$ happens to give the smallest error in this table, it is near the turnover and can be sensitive to evaluation details. A more robust fixture is $h=10^{-4}$, where the expected truncation error is about $4.53\times10^{-9}$ and second-order behavior is still clear. A tolerance of order $10^{-8}$ is defensible for this stated platform and formula after the preceding refinement and precision checks. It is not a universal tolerance for every implementation or function.

2. Quantify the cost of correlation

For the covariance model with $N=1000$, $\rho=0.8$, and $\sigma^2=1$, compute the large-$N$ effective sample size and standard error. Compare them with the independent-sample values and explain what additional checks real data require.

Solution

The large-$N$ effective sample size is

$$N_{\mathrm{eff}} \sim 1000\frac{1-0.8}{1+0.8} = 111.1.$$

Thus

$$\operatorname{SE}(\overline X) \approx \frac1{\sqrt{111.1}} \approx0.0949,$$

whereas the independent formula gives $N_{\mathrm{eff}}=1000$ and $\operatorname{SE}=0.0316$. Correlation has reduced the information by about a factor of nine and increased the standard error by about a factor of three. The exact finite-$N$ values are $N_{\mathrm{eff}}\approx111.6$ and $\operatorname{SE}\approx0.0947$.

Real data require evidence for equilibration or stationarity, an estimated autocorrelation tail or block length, stability under reasonable window or block changes, and enough independent streams or history to expose slow modes. The algebra with known $\rho$ does not supply those checks.

3. Classify three repeat calculations

Classify each case and state what it does not test:

1. A colleague uses the recorded command, same source, same inputs, and locked dependencies on a clean machine.
2. A second analyst receives the same raw output but writes a new fitting script and varies the justified fit window.
3. A separate group implements the published algorithm from the mathematical specification and regenerates the synthetic data with an independent random stream.

Solution

The first is a rerun. It tests whether the preserved procedure and environment description are sufficient and exposes some dependency or nondeterminism problems. It shares the implementation, inputs, and model, so it does not test code defects or data generation.

The second is a reanalysis. It tests downstream analysis code and fit-window sensitivity. It still shares the raw data and every upstream generation or measurement defect.

The third is an independent reproduction with both implementation and synthetic generation independence. It adds strong protection against shared code and random-stream defects. The groups can nevertheless agree because they share the same incorrect mathematical specification or scientific model, so agreement alone does not establish a physical conclusion.

4. Decide what more sampling can repair

A result quotes a $1.2\%$ correlation-aware statistical error, a $0.8\%$ discretization estimate, a shared $0.5\%$ calibration uncertainty, and a $3\%$ spread across defensible physical models. May these four percentages be added in quadrature automatically? Which component is reduced by generating four times as much effectively independent data?

Solution

No. Quadrature represents the variance of a sum only under an appropriate probabilistic model with the relevant covariances. The statistical term may have such an interpretation. The discretization number might be a bound or an extrapolation sensitivity; the calibration can be shared across all samples; and the model spread is not automatically a random variable with known coverage. Each needs its own meaning and dependence structure before any combination.

Four times as much effectively independent data reduces the statistical standard error by a factor of two, from $1.2\%$ to about $0.6\%$. It does not automatically reduce the discretization, calibration, or model components. Those require refinement, better calibration, alternative observables or models, and scientific validation.

Re-check and return

Repeat the computational and evidence diagnostic with a numerical integral whose exact value is known and a newly generated correlated sequence. Your result is ready to carry forward when it contains:

• a precise target, units, regime, and bounded reference case;
• tests that demonstrably fail under at least one plausible defect;
• a refinement study with an observed order and a precision repeat;
• correlation-aware uncertainty and its stationarity or equilibration assumptions;
• a complete reproducibility packet tested in a clean rerun or independent analysis;
• explicit shared dependencies for every repeat; and
• separate execution, numerical-verification, uncertainty, and scientific claims, with unsupported claims marked as such.

If all seven elements are explicit, continue to Computational field theory onboarding or Reproduce and validate a result. If only the covariance analysis remains weak, return to Statistical ensembles and probability and repeat the correlated-mean exercise with a new value of $\rho$. If the refinement or precision behavior is unexplained, do not enlarge the tolerance; isolate the deterministic calculation before attaching it to a scientific result.

References

• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, second edition, Society for Industrial and Applied Mathematics, 2002, doi:10.1137/1.9780898718027.
• National Academies of Sciences, Engineering, and Medicine, Reproducibility and Replicability in Science, National Academies Press, 2019, doi:10.17226/25303.
• Alan D. Sokal, “Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms,” in Cécile DeWitt-Morette, Pierre Cartier, and Antoine Folacci, eds., Functional Integration: Basics and Applications, Springer, 1997, pp. 131–192, doi:10.1007/978-1-4899-0319-8_6.
• Ulli Wolff, “Monte Carlo Errors with Less Errors,” Computer Physics Communications 156 (2004): 143–153, doi:10.1016/S0010-4655(03)00467-3.