Uncertainty Propagation: Method

GUM Law of Propagation and the Monte Carlo Method

The measurement model

A result is rarely read directly; it is computed from other quantities, Y = f(X₁, X₂, …, X_n). The density of a cylinder is its mass over its volume, a resistance is a voltage over a current, a strength allowable is a mean less a factor times a standard deviation. Each input quantity X_i is known imperfectly: the best estimate x_i is accompanied by a standard uncertainty u(x_i), the standard deviation of the probability distribution that expresses what is known about X_i. The question the tool answers is what uncertainty u(y) the output inherits, and what interval about y contains the value of Y with a stated probability. The model is typed in plain formula syntax with the input names; the tool does not track units, so every parameter of an input is entered in that input's unit and the result carries the unit the formula gives (mass in g and lengths in cm give g/cm³).

Standard uncertainty of the inputs (JCGM 100 clause 4, JCGM 101 clause 6)

A Type A evaluation uses repeated observations: the best estimate is the mean of n readings, u = s/√n from their sample standard deviation s, with ν = n − 1 degrees of freedom. JCGM 101 (§6.4.9) assigns such a quantity a scaled and shifted Student t distribution with ν = n − 1 and scale s/√n; the tool offers it as "Student t (mean, scale, ν)". A Type B evaluation uses other knowledge, a certificate, a specification, a resolution or experience, and assigns a distribution to it. The tool's distributions and the standard uncertainty of each are:

distribution	parameters entered	standard uncertainty u	typical use
Normal	mean, u	u	a Type A mean with many readings, or any quantity whose uncertainty is already a standard deviation
Normal from expanded	mean, U, k	U/k	a certificate that states an expanded uncertainty U with its coverage factor k (k = 2 when it only says "approximately 95 %")
Rectangular	mean, half-width a	a/√3	a quantity equally likely anywhere within ± a: a tolerance, a resolution (a = half the last digit), rounded data
Triangular	mean, half-width a	a/√6	limits ± a with values near the middle more likely than near the limits
U-shaped (arcsine)	mean, half-width a	a/√2	a quantity that cycles between its limits, such as a thermostatically controlled temperature (GUM H.1.3.4)
Lognormal	mean, sd	sd	a positive quantity with a long upper tail (a cost, a duration, a concentration); the tool converts the mean and sd to the log-scale parameters μln = ln(mean) − σln²/2, σln² = ln(1 + (sd/mean)²)
Student t	mean, scale s, ν	s (the scale; the distribution's own standard deviation is s √(ν/(ν − 2)))	a Type A mean of few readings, scale s/√n and ν = n − 1 (JCGM 101 §6.4.9); ν must exceed 2
PERT	min a, most likely m, max b	√((μ − a)(b − μ)/7), μ = (a + 4m + b)/6	a three-point expert estimate (project and risk analysis); the beta distribution with shape 1 + 4(m − a)/(b − a) and 1 + 4(b − m)/(b − a)
Constant	value	0	a defined constant or a quantity whose uncertainty is negligible

The degrees of freedom ν_i measure how reliable u(x_i) itself is: n − 1 for a Type A evaluation, the value stated on a certificate, and infinite for a Type B estimate whose limits are trusted. The GUM (G.4.2) allows a finite ν for a Type B estimate from the judged relative uncertainty of u itself, ν ≈ ½ (Δu/u)⁻²: limits believed reliable to 25 % give ν = 8, to 50 % give ν = 2 (its Annex H.1 uses both).

The law of propagation of uncertainty (JCGM 100 clause 5)

Expanding the model to first order about the best estimates, Y − y ≈ Σ c_i (X_i − x_i) with the sensitivity coefficients c_i = ∂f/∂x_i, and taking the variance of both sides,

u_c²(y) = Σ c_i² u²(x_i) + 2 Σ_i<j c_i c_j u(x_i) u(x_j) r(x_i, x_j),

the second sum vanishing for independent inputs. The tool differentiates the formula numerically: central differences (f(x + h) − f(x − h))/2h at two step sizes h and h/2, combined by Richardson extrapolation (4D(h/2) − D(h))/3, which is exact for polynomials up to the fourth degree and has an error of order h⁴ otherwise, with h = 10⁻³ of the input's magnitude (of its u when the best estimate is zero). The uncertainty budget tabulates x_i, u(x_i), c_i, the contribution |c_i| u(x_i) and its share (c_i u(x_i))²/u_c² of the variance; Figure 1 ranks the inputs by that share, which is the first thing to read when the uncertainty has to be reduced, since only the dominant inputs are worth improving. Correlation is entered as pairs "name name r" or as the whole matrix, which must be symmetric with a unit diagonal and positive definite; a positive correlation between two inputs that enter with sensitivities of the same sign increases u_c, of opposite signs decreases it (a difference of two readings from the same instrument).

Effective degrees of freedom and the coverage factor (JCGM 100 Annex G)

The expanded uncertainty U = k u_c defines the interval y ± U that is expected to contain Y with the coverage probability p. When every input has infinite degrees of freedom the output is taken as normal and k is the normal quantile, 1.960 for 95 %, 2.576 for 99 %. When some u(x_i) are themselves uncertain the output follows a t distribution with the effective degrees of freedom of the Welch–Satterthwaite formula,

ν_eff = u_c⁴(y) / Σ (c_i u(x_i))⁴/ν_i,

the terms with infinite ν_i contributing nothing. The GUM (G.4.1) truncates ν_eff to the next lower integer and reads k = t_p(ν_eff) from the t distribution, 2.26 for ν = 9 and 2.04 for ν = 30 at 95 %. The formula is derived for independent inputs; with correlation entered it is reported as indicative. Its Annex H.1 end-gauge calibration (u_c = 32 nm, ν_eff = 16.7 → 16, k₉₉ = 2.92, U = 93 nm) is one of the checks the tool's test script reproduces.

Why a Monte Carlo evaluation (JCGM 101)

The law of propagation rests on three things that may fail: the first-order expansion (a strongly nonlinear model, or an input whose best estimate is zero so the linear term vanishes, as in the mass-calibration example of JCGM 101 §9.3), the assumption that the output is normal or t distributed (inputs far from normal that dominate the budget, as the rectangular diameter does in the default case here), and symmetry of the output distribution (a lognormal input, a model like 1/x or e^x with a large relative uncertainty). JCGM 101 propagates the distributions themselves instead of their variances: each trial draws every input from its distribution, evaluates the model, and the M model values are a sample of the output distribution. Its estimate is their mean, its standard uncertainty their standard deviation, and the coverage intervals are read from the sorted values with no assumption about the shape. The tool draws rectangular, triangular and U-shaped values by their inverse distribution functions, normal and lognormal values by the polar Box–Muller method, Student t values as a normal over the root of a chi-square from a gamma variate (Marsaglia and Tsang), and PERT values as a beta variate from two gamma variates; the generator is the Mersenne Twister with the entered seed, so a run is repeatable.

Number of trials and the adaptive procedure (JCGM 101 §7.2, §7.9)

With a fixed count, M = 10⁵ trials give the 95 % interval endpoints to about one hundredth of u(y) and 10⁶ to a few thousandths; JCGM 101 recommends M at least 10⁴ times 1/(1 − p). The adaptive procedure (§7.9.4) instead runs blocks of max(10⁴, 100/(1 − p)) trials and after each block computes the standard deviation of the block estimates of y, u(y) and the two endpoints about their averages, s = √(Σ (y^(r) − ȳ)²/(h(h − 1))); when 2s is below the numerical tolerance δ for all four, the result has stabilized and all h × M values are used. The tolerance is defined from the number of significant digits wanted in u(y): writing u(y) = c × 10^l with c a two-digit integer, δ = ½ × 10^l (u = 2.0 gives δ = 0.05, u = 0.0074 gives δ = 0.00005). The tool caps the run at the entered trial count and at about eight seconds.

The two coverage intervals (JCGM 101 §7.7)

Any interval that contains 100p % of the model values is a coverage interval. The probabilistically symmetric one leaves (1 − p)/2 of the values below and above it; the shortest one is the narrowest of all such intervals, found by sliding a window of pM values along the sorted sample. For a symmetric output distribution the two coincide; for a skewed one the shortest interval moves toward the mode and is the one JCGM 101 recommends. The histogram of Figure 2 shows the shape against the normal distribution the GUM result implies, scaled to the same area; the skewness g₁ and excess kurtosis g₂ quantify the departure (a rectangular distribution has g₂ = −1.2; the sum of four of them −0.3). Figure 3 shows how the standard deviation settles as the trials accumulate.

Correlated inputs in the Monte Carlo (JCGM 101 §6.4.8, Annex C.5)

Correlated inputs are drawn from a multivariate normal distribution by the Cholesky factor L of the correlation matrix: a vector of independent standard normal scores w becomes z = L w with the required correlation, and each z_i is mapped through its input's quantile function. For normal and lognormal inputs this is exact. For the other distributions it is a Gaussian copula: the normal scores carry the entered correlation and the marginals keep their shapes, but the linear correlation of the drawn values is somewhat smaller than r (two rectangular inputs at r = 0.8 have a Pearson correlation of (6/π) arcsin(0.4) = 0.786, Kruskal 1958). Student t and PERT quantile functions are tabulated against the normal score and interpolated. The GUM form uses the entered r directly.

Validation of the GUM result (JCGM 101 §8.1)

Both methods are run on the same inputs, and the GUM interval y ± U is compared with the Monte Carlo probabilistically symmetric interval at the same coverage probability: d_low = |y − U − y_low| and d_high = |y + U − y_high|. When both are within the numerical tolerance δ, the GUM framework is validated for this model and these inputs and either result may be reported; when either exceeds δ, the first-order treatment is inadequate and the Monte Carlo result should be reported. JCGM 101 §9.2 shows the test on the additive model Y = X₁ + X₂ + X₃ + X₄: with four standard normal inputs both methods give [−3.92, 3.92]; with four rectangular inputs of the same standard deviation the GUM still gives ±3.92 but the true 95 % interval is ±2√3[2 − (3/5)^¼] = ±3.88 (its Annex E), a difference the two-digit tolerance δ = 0.05 just resolves. The test is strict by design: in the default case here the rectangular diameter carries 55 % of the variance, the output is flat-topped (g₂ ≈ −0.34), and the GUM interval is 4 % wider than the Monte Carlo one, which at δ = 0.00005 is a failure of validation even though the GUM result is the conservative one. Validation asks whether the two agree to the digits reported, not which is safer.

Reporting (JCGM 100 clause 7)

State the estimate and either the combined standard uncertainty u_c or the expanded uncertainty U with its coverage factor and probability, give the uncertainty to two significant digits and round the estimate to the same decimal place, and list the budget so that a reader can see what dominates. When the Monte Carlo result is reported, give the number of trials and the interval type (symmetric or shortest).

Limits

One output quantity (a vector output and its covariance are JCGM 102). The GUM evaluation is first order only; the second-order terms of JCGM 100 eq. 5.1.2 note are not added, which is why the Monte Carlo comparison matters for nonlinear models. The input distributions are those listed, assigned by the analyst; the tool does not form them from data (Bayesian updating of a prior with observations, JCGM 101 §6.1) and does not fit distributions to samples. Correlation among non-normal inputs is a Gaussian copula, not a joint distribution of another form. Welch–Satterthwaite is applied as published, for independent inputs. The Monte Carlo run is capped at 400,000 trials and about eight seconds, so a 99.9 % interval endpoint is resolved to a few hundredths of u only.