Uncertainty I: random and systematic errors, propagation

Lesson 1 introduced the two kinds of error; this lesson makes them quantitative. We describe random scatter with the normal distribution and its standard deviation, distinguish the scatter of single measurements from the uncertainty of their mean, show how to estimate uncertainties that cannot be found by repetition, derive the formula for propagating uncertainties through a calculation, test it against a Monte Carlo simulation, and learn how to combine several measurements of the same quantity with a weighted mean. CHEM-211 treated much of the same statistics from the analytical chemist's side, with t-tests and calibration curves; here the emphasis is on the physicist's problem of turning measured inputs into a derived result with an honest uncertainty.

Describing random scatter

Repeat a measurement many times and the results scatter. Very often the histogram of results approaches a normal, or Gaussian, distribution,

with mean \(\mu\) and standard deviation \(\sigma\). The central limit theorem explains why: when a measurement error is the sum of many small independent contributions, their sum tends to a normal distribution whatever the individual distributions look like. About 68 of every 100 results fall within \(\pm\sigma\) of the mean, 95 within \(\pm 2\sigma\) and 99.7 within \(\pm 3\sigma\). From \(N\) measurements \(x_i\) we estimate the mean and standard deviation by

The factor \(N - 1\) rather than \(N\) corrects for the fact that the deviations are measured from \(\bar{x}\), which is itself fitted to the data and so is closer to them than the true mean is. The standard deviation \(s\) describes the scatter of a single measurement. The uncertainty in the mean is smaller, the standard error

This is the basic reason repetition helps: the average of 100 measurements is ten times more certain than any one of them. With few measurements, \(s\) itself is poorly known, and a confidence interval for the mean should use Student's t distribution, introduced by William Sealy Gosset in 1908 while he worked on quality control at the Guinness brewery; with ten measurements, the 95 in 100 interval is \(\pm 2.26\,\sigma_{\bar{x}}\) rather than \(\pm 1.96\,\sigma_{\bar{x}}\).

Type A and Type B uncertainties

The international Guide to the Expression of Uncertainty in Measurement, published in 1993 and known as the GUM, classifies uncertainty estimates by how they are obtained rather than by the kind of error. Type A estimates come from statistical analysis of repeated measurements, as above. Type B estimates come from everything else: the resolution of an instrument, the manufacturer's accuracy specification, a calibration certificate, a handbook value or physical judgment. A ruler marked in millimeters has a reading uncertainty of order a fraction of a millimeter; a digital display that changes in steps of 0.01 V has a resolution uncertainty. If a quantity is equally likely to lie anywhere in an interval of half-width \(a\), as with rounding to the last digit, its standard uncertainty is \(a/\sqrt{3}\). A specification such as plus or minus 0.5 parts in 100 of reading plus two digits must be converted the same way. Systematic errors you know about should be corrected; the uncertainty of the correction is a Type B contribution.

Propagating uncertainties

Usually the quantity you want is calculated from several measured ones. If \(q = f(x, y, \dots)\) and the inputs have small, independent uncertainties \(\sigma_x, \sigma_y, \dots\), expand \(f\) to first order about the best estimates. A small change in \(x\) changes \(q\) by \((\partial f/\partial x)\,\delta x\). Independent changes add like perpendicular vectors, because the cross terms average to zero, so the variance of \(q\) is the sum of the variances contributed by each input:

To see why the cross terms vanish, write \(\delta q = a\,\delta x + b\,\delta y\) and square: \(\delta q^2 = a^2\delta x^2 + b^2\delta y^2 + 2ab\,\delta x\,\delta y\). Averaged over many repeats, \(\langle\delta x^2\rangle = \sigma_x^2\), \(\langle\delta y^2\rangle = \sigma_y^2\), and \(\langle\delta x\,\delta y\rangle = 0\) if the errors are independent. When the errors are correlated, the cross term remains as \(2ab\,\sigma_{xy}\), where \(\sigma_{xy}\) is the covariance; this matters, for example, when two quantities are computed from the same calibration. Several special cases follow directly. For sums and differences, \(q = x \pm y\), the absolute uncertainties add in quadrature. For products and quotients, \(q = x^m y^n\), the fractional uncertainties add in quadrature, weighted by the powers:

For the pendulum, \(g = 4\pi^2LT^{-2}\), so the fractional uncertainty of \(g\) combines the fractional uncertainty of \(L\) with twice that of \(T\). The factor of two makes timing more important than it first appears. Two warnings apply. A difference of nearly equal quantities can have a huge fractional uncertainty, even when each input is precise: \(1.003 \pm 0.002\) minus \(1.000 \pm 0.002\) is \(0.003 \pm 0.003\). And the linear formula assumes that \(f\) is close to linear over the range of the uncertainties; when it is not, the formula can mislead badly.

Monte Carlo propagation

A Monte Carlo calculation sidesteps the derivatives. Draw each input at random from its distribution, compute \(q\), repeat many thousands of times, and take the mean and standard deviation of the results. The method handles correlations, non-normal inputs and strong nonlinearity, at the cost of computing time, which is now cheap. It is also a valuable check on an analytic calculation. The program propagates the pendulum uncertainties both ways, then treats a deliberately nonlinear case: the range of a projectile launched at 12 m/s at 44 degrees, close to the 45-degree maximum, with an angle uncertainty of 3 degrees. Finally it combines three measurements of \(g\) with a weighted mean.

import math, random
random.seed(397)
# 1. g = 4 pi^2 L / T^2 with L = 0.800 +/- 0.002 m and T = 1.7975 +/- 0.0015 s.
L, dL, T, dT = 0.800, 0.002, 1.7975, 0.0015
g = 4*math.pi**2*L/T**2
dg = g*math.sqrt((dL/L)**2 + (2*dT/T)**2)
print(f"analytic: g = {g:.4f} +/- {dg:.4f} m/s^2  (length term {g*dL/L:.4f}, period term {g*2*dT/T:.4f})")
N = 200000
samp = [4*math.pi**2*random.gauss(L, dL)/random.gauss(T, dT)**2 for _ in range(N)]
m = sum(samp)/N; s = math.sqrt(sum((x - m)**2 for x in samp)/(N - 1))
print(f"Monte Carlo ({N} draws): g = {m:.4f} +/- {s:.4f} m/s^2")
# 2. A strongly nonlinear case: the range of a projectile at a steep launch angle, where linear propagation fails.
v, th, dth = 12.0, math.radians(44.0), math.radians(3.0)
R = lambda a: v*v*math.sin(2*a)/9.81
lin = abs(2*v*v*math.cos(2*th)/9.81)*dth
mc = [R(random.gauss(th, dth)) for _ in range(N)]
mm = sum(mc)/N; ss = math.sqrt(sum((x - mm)**2 for x in mc)/(N - 1))
print(f"range at 44 +/- 3 degrees: linear propagation gives +/- {lin:.3f} m, Monte Carlo gives mean {mm:.3f} m (nominal {R(th):.3f}) and spread +/- {ss:.3f} m")
# 3. Weighted mean of three measurements of the same quantity.
meas = [(9.796, 0.012), (9.812, 0.006), (9.781, 0.020)]
w = [1/s**2 for _, s in meas]
xw = sum(wi*x for wi, (x, _) in zip(w, meas))/sum(w)
sw = 1/math.sqrt(sum(w))
chi2 = sum(((x - xw)/s)**2 for x, s in meas)
print(f"weighted mean = {xw:.4f} +/- {sw:.4f}, chi-squared = {chi2:.2f} for 2 degrees of freedom")

analytic: g = 9.7749 +/- 0.0294 m/s^2  (length term 0.0244, period term 0.0163)
Monte Carlo (200000 draws): g = 9.7749 +/- 0.0294 m/s^2
range at 44 +/- 3 degrees: linear propagation gives +/- 0.054 m, Monte Carlo gives mean 14.590 m (nominal 14.670) and spread +/- 0.125 m
weighted mean = 9.8069 +/- 0.0052, chi-squared = 3.22 for 2 degrees of freedom

For the pendulum the two methods agree perfectly, \(g = 9.7749 \pm 0.0294\ \text{m/s}^2\), because the function is close to linear over the small input uncertainties. The breakdown shows that the length contributes more than the period here, so a better length measurement is the next improvement. The projectile is different. Near 45 degrees the range is at a maximum, so its first derivative with respect to angle is almost zero and the linear formula gives only \(\pm 0.054\) m. But any error in angle, in either direction, can only decrease the range. The Monte Carlo result shows both effects: the spread is more than twice the linear estimate, and the mean range is shifted 8 cm below the nominal value. Near a maximum, the distribution of the result is lopsided and the linear rule fails.

Combining measurements: the weighted mean

When several independent measurements \(x_i \pm \sigma_i\) of the same quantity are available, the best combined estimate gives each a weight \(w_i = 1/\sigma_i^2\):

The formula follows from maximizing the joint likelihood of the measurements, or equivalently from minimizing \(\chi^2 = \sum_i (x_i - \bar{x})^2/\sigma_i^2\) with respect to \(\bar{x}\): setting the derivative to zero gives the weighted mean, and propagating the uncertainties gives its error. A measurement twice as precise counts four times as much. The value of \(\chi^2\) at the minimum is a consistency check. If the measurements agree within their errors, \(\chi^2\) should be about equal to the number of degrees of freedom, here the number of measurements minus one. In the example, \(\chi^2 = 3.22\) for two degrees of freedom, which happens about one time in five by chance: acceptable. A value of 20 would say that at least one measurement has an underestimated uncertainty or a systematic error, and the weighted mean should not be trusted until the disagreement is understood. The Particle Data Group, which compiles the world's measurements of particle properties, inflates the combined error by \(\sqrt{\chi^2/\nu}\) when measurements disagree, an honest admission of unknown systematics.

Estimating systematic uncertainties

Systematic errors are found by thinking about the apparatus and by varying the method. Calibrate instruments against standards. Reverse things that should not matter, such as the direction of a current or the orientation of a sample, and check that the result is unchanged. Measure the same quantity a second, independent way. Vary parameters that the theory says are irrelevant, such as the mass of a pendulum bob, and look for changes. Once a systematic effect is identified, correct for it and assign an uncertainty to the correction. Then combine all the contributions, random and systematic, in quadrature in the final error budget, and report them separately when they are comparable, as in \(g = 9.807 \pm 0.006\ (\text{stat}) \pm 0.010\ (\text{syst})\ \text{m/s}^2\), so that readers can see which kind dominates.