Error propagation in Python with SciPy: linear, Monte Carlo, uncertainties (tested)

To propagate errors with SciPy, take numerical partial derivatives with scipy.optimize.approx_fprime and combine them with the input sigmas (or the covariance matrix), or sample the inputs with scipy.stats and run a Monte Carlo. For everyday arithmetic, the uncertainties package does the linear propagation for you. All three agree on the example below.

1. Three ways, one example: g from a pendulum

b1.py propagates the uncertainty of a length (2.00 ± 0.05 m) and a period (2.84 ± 0.03 s) into g = 4π²L/T²:

import numpy as np
from scipy import stats
from scipy.optimize import approx_fprime
import uncertainties
from uncertainties import ufloat

# measured: length L = 2.00 +/- 0.05 m, period T = 2.84 +/- 0.03 s -> g = 4 pi^2 L / T^2
L, sL, T, sT = 2.00, 0.05, 2.84, 0.03
def g(p):
    return 4 * np.pi**2 * p[0] / p[1]**2

# 1. first-order (linear) propagation with a numerical gradient
p = np.array([L, T])
J = approx_fprime(p, g, 1e-8)
sg_lin = np.sqrt((J[0] * sL)**2 + (J[1] * sT)**2)
print("g            ", round(g(p), 4))
print("linear sigma ", round(sg_lin, 4))

# 2. Monte Carlo with scipy.stats
rng = np.random.default_rng(0)
Ls = stats.norm(L, sL).rvs(200_000, random_state=rng)
Ts = stats.norm(T, sT).rvs(200_000, random_state=rng)
gs = g([Ls, Ts])
print("MC mean, std ", round(gs.mean(), 4), round(gs.std(ddof=1), 4))

# 3. the uncertainties package
gu = 4 * np.pi**2 * ufloat(L, sL) / ufloat(T, sT)**2
print("uncertainties", gu, "| version", uncertainties.__version__)

# 4. correlated inputs: pass the covariance, not two sigmas
cov = np.array([[sL**2, 0.8 * sL * sT], [0.8 * sL * sT, sT**2]])
print("correlated   ", round(float(np.sqrt(J @ cov @ J)), 4))

Output (SciPy 1.18.1, NumPy 2.5.4, uncertainties 3.2.3, Python 3.13.5):

g             9.7893
linear sigma  0.3204
MC mean, std  9.7926 0.3215
uncertainties 9.79+/-0.32 | version 3.2.3
correlated    0.1473

• Linear (first-order) propagation: σ_g² = Σ (∂g/∂x_i · σ_i)². scipy.optimize.approx_fprime gives the partial derivatives numerically, so you don't have to derive them by hand.
• Monte Carlo: sample the inputs with scipy.stats.norm(...).rvs, push the samples through the formula, take the standard deviation. 0.3215 vs 0.3204 here: they agree because the relative errors are small. When they disagree, the function is too non-linear for the linear formula and the Monte Carlo number is the one to trust.
• uncertainties package: ufloat(value, sigma) objects carry the linear propagation through ordinary arithmetic and print as 9.79+/-0.32.
• Correlated inputs: with a correlation of 0.8 the result drops from 0.32 to 0.15, because L and T enter with opposite signs. Use sqrt(J @ cov @ J) with the full covariance matrix; adding the sigmas in quadrature assumes independent inputs.

2. Fitted parameters

If your inputs come from curve_fit, use the returned pcov as the covariance in point 4; the parameters of a fit are usually correlated. Standard errors themselves: SciPy standard error (stats.sem, linregress, curve_fit).

Plot the result with error bars: matplotlib errorbar, plotly error_y.

Written with AI assistance. Measured on 2026-10-23 with SciPy 1.18.1, NumPy 2.5.4, uncertainties 3.2.3 and Python 3.13.5 on Linux. The script ran verbatim and the output was pasted unedited. PlainTables is our extension; SciPy is free.