Gaussian Error Propagation with SymPy

Physikalisches Fortgeschrittenen-Praktikum Heidelberg, Klaus Reygers

This notebook shows how to use symbolic derivatives to propagate independent measurement uncertainties. The example is the volume of a cylinder,

\[V = \pi r^2 h,\]

where the measured radius \(r\) and height \(h\) have uncertainties \(\sigma_r\) and \(\sigma_h\).

Gaussian error propagation

For a quantity \(f(x_1, x_2, \ldots)\) that depends on independent measured variables, the variance is approximated by

\[ \sigma_f^2 = \sum_i \left(\frac{\partial f}{\partial x_i}\sigma_{x_i}\right)^2. \]

The following helper implements this formula symbolically with SymPy.

import sympy as sp
from IPython.display import Math, display

def gaussian_uncertainty(expr, variables):
    """Return the Gaussian uncertainty of expr for independent variables.

    Parameters
    ----------
    expr : sympy expression
        Quantity calculated from the measured variables.
    variables : sequence of (sympy.Symbol, sympy.Symbol)
        Pairs of measured variables and their standard uncertainties,
        e.g. [(x, sigma_x), (y, sigma_y)].
    """
    if not variables:
        raise ValueError("At least one variable-uncertainty pair is required.")

    variance = sp.S.Zero
    for variable, uncertainty in variables:
        derivative = sp.diff(expr, variable)
        variance += derivative**2 * uncertainty**2

    return sp.sqrt(sp.simplify(variance))

Example: volume of a cylinder

First define the symbols and the formula. The symbols are declared positive so that SymPy can simplify square roots in a physically meaningful way.

r, h, sigma_r, sigma_h = sp.symbols("r h sigma_r sigma_h", positive=True)
volume = sp.pi * r**2 * h

volume

\(\displaystyle \pi h r^{2}\)

Now apply the propagation formula. The uncertainty is calculated from the derivatives with respect to \(r\) and \(h\).

sigma_volume = gaussian_uncertainty(volume, [(r, sigma_r), (h, sigma_h)])

display(Math(rf"\sigma_V = {sp.latex(sigma_volume)}"))

\(\displaystyle \sigma_V = \pi r \sqrt{4 h^{2} \sigma_{r}^{2} + r^{2} \sigma_{h}^{2}}\)

It is often useful to look at the relative uncertainty. For the cylinder volume, the radius uncertainty enters with a factor of two because \(V\) depends on \(r^2\).

relative_uncertainty = sp.simplify(sigma_volume / volume)

display(Math(rf"\frac{{\sigma_V}}{{V}} = {sp.latex(relative_uncertainty)}"))

\(\displaystyle \frac{\sigma_V}{V} = \frac{\sqrt{4 h^{2} \sigma_{r}^{2} + r^{2} \sigma_{h}^{2}}}{h r}\)

Insert measured values

Use a dictionary for substitutions. This keeps the numerical calculation readable and reduces the chance of mixing up a variable and its uncertainty.

measurements = {
    r: 3.0,          # cm
    sigma_r: 0.1,    # cm
    h: 5.0,          # cm
    sigma_h: 0.1,    # cm
}

volume_value = float(volume.subs(measurements))
sigma_value = float(sigma_volume.subs(measurements))

display(Math(rf"V = ({volume_value:.1f} \pm {sigma_value:.1f})\,\mathrm{{cm}}^3"))

\(\displaystyle V = (141.4 \pm 9.8)\,\mathrm{cm}^3\)

Which measurement matters most?

The variance is a sum of squared contributions. Inspecting the contributions is a good habit: it tells you where a more precise measurement would improve the final result most.

contributions = {
    "radius": (sp.diff(volume, r) * sigma_r)**2,
    "height": (sp.diff(volume, h) * sigma_h)**2,
}

numeric_contributions = {
    name: float(term.subs(measurements))
    for name, term in contributions.items()
}
total_variance = sum(numeric_contributions.values())

for name, variance_part in numeric_contributions.items():
    fraction = variance_part / total_variance
    print(f"{name:>6}: {100 * fraction:5.1f}% of the variance")

radius:  91.7% of the variance
height:   8.3% of the variance