Prev
Unit 12 · Randomness and Estimation
Propagating Uncertainty
65m100 XP
Next
🎲Randomness and Estimation

Propagating Uncertainty

Covariance and the multivariate normal, error ellipses, linearised propagation J P Jᵀ, the unscented transform, collision probability

Over Siberia, 10 February 2009
Iridium 33, an active communications satellite, and Kosmos 2251, a Russian military satellite that had been out of service for years, are both in near-polar orbits about 790 km up over Siberia, crossing each other's paths almost at right angles. Every day an online service called SOCRATES, run by T. S. Kelso at CelesTrak, screens the public orbit catalogue for close approaches between satellites. Its report that morning predicts that these two will pass 584 m apart at 16:56 UTC. It is not an alarming number. Orbits computed from the public catalogue are uncertain by kilometres, and predicted passes closer than that happen every day; in the week's reports this pair never makes the top ten, and on the day it ranks 152nd. At 16:56 UTC Iridium 33 stops transmitting. The two satellites have collided at about 11.7 km/s, destroying both and filling their orbits with debris — more than 2,000 catalogued fragments by mid-2011. The prediction of the closest approach had been reasonable. What it did not carry was a statement of how uncertain it was, and therefore of how likely a collision was.
Every quantity an engineer computes from measured or uncertain inputs is itself uncertain. For a vector of uncertain quantities x\mathbf x the essential description has two parts: the mean μ=E[x]\boldsymbol\mu = E[\mathbf x] and the covariance matrix
P=E[(x−μ)(x−μ)T].P = E\big[(\mathbf x - \boldsymbol\mu)(\mathbf x - \boldsymbol\mu)^{\mathsf T}\big].
Its diagonal holds the variances σi2\sigma_i^2; the off-diagonal entries Pij=ρijσiσjP_{ij} = \rho_{ij}\sigma_i\sigma_j say how errors in different components move together. PP is symmetric and positive semidefinite. The multivariate normal distribution is the one fixed by these two alone, with density proportional to exp⁡ ⁣(−12(x−μ)TP−1(x−μ))\exp\!\big(-\tfrac12(\mathbf x - \boldsymbol\mu)^{\mathsf T}P^{-1}(\mathbf x - \boldsymbol\mu)\big). The quadratic form in the exponent is the squared Mahalanobis distance: distance measured in standard deviations, along the directions in which the errors actually lie.
💡Error ellipses are eigenvectors
The surfaces of constant density are ellipsoids. Their axes point along the eigenvectors of PP and their semi-axes are kλik\sqrt{\lambda_i} (Unit 6). For a 95% region in two dimensions, take k2=5.991k^2 = 5.991, the 95th percentile of the χ2\chi^2 distribution with two degrees of freedom — not k=2k = 2, which encloses only 86%. A GNSS receiver's horizontal accuracy, a surveyor's point error and an orbit's along-track/cross-track uncertainty are all reported as such ellipses. Least squares (Unit 5) produces one automatically: the estimate p^\hat{\mathbf p} of Ap≈bA\mathbf p \approx \mathbf b with independent errors σ\sigma has covariance σ2(ATA)−1\sigma^2(A^{\mathsf T}A)^{-1} — the matrix whose diagonal gave GNSS its dilution of precision in lesson 5-4.
Derivation
Ellipse axes from eigenvectors, and why 95% needs k2=5.991k^2 = 5.991
The axes. Write P=QΛQTP = Q\Lambda Q^{\mathsf T} with orthonormal eigenvectors in QQ (a symmetric matrix always allows this, lesson 5-3). In the rotated coordinates z=QT(x−μ)\mathbf z = Q^{\mathsf T}(\mathbf x - \boldsymbol\mu) the quadratic form loses its cross terms:
(x−μ)TP−1(x−μ)=zTΛ−1z=z12λ1+z22λ2.(\mathbf x - \boldsymbol\mu)^{\mathsf T}P^{-1}(\mathbf x - \boldsymbol\mu) = \mathbf z^{\mathsf T}\Lambda^{-1}\mathbf z = \frac{z_1^2}{\lambda_1} + \frac{z_2^2}{\lambda_2}.
Setting it equal to k2k^2 gives an ellipse with semi-axes kλ1k\sqrt{\lambda_1} and kλ2k\sqrt{\lambda_2} along the eigenvectors.
The probability inside. Each zi/λiz_i/\sqrt{\lambda_i} is an independent standard normal, so d2=w12+w22d^2 = w_1^2 + w_2^2 for two standard normals — exactly the radius of the Box–Muller derivation in lesson 12-1, where P(d2>s)=e−s/2P(d^2 > s) = e^{-s/2}. So
P(d2≤k2)=1−e−k2/2.P(d^2 \le k^2) = 1 - e^{-k^2/2}.
For 95%, e−k2/2=0.05e^{-k^2/2} = 0.05 and k2=−2ln⁡0.05=5.991k^2 = -2\ln0.05 = 5.991. With k=2k = 2 the formula gives 1−e−2=0.8651 - e^{-2} = 0.865: a '2-sigma' ellipse holds only 86% in two dimensions.
y=f(x)⟹Py≈J Px JT,J=∂f∂x∣μ\mathbf y = \mathbf f(\mathbf x) \quad\Longrightarrow\quad P_y \approx J\,P_x\,J^{\mathsf T}, \qquad J = \frac{\partial\mathbf f}{\partial\mathbf x}\Big|_{\boldsymbol\mu}

Linearised propagation of covariance — the 'law of propagation of uncertainty' of the international measurement Guide (GUM). Exact when f is linear; the mean maps to f(μ).

Derivation
Where JPJTJPJ^{\mathsf T} comes from — and what it leaves out
The linear part. Expand f\mathbf f about the mean: y≈f(μ)+J(x−μ)\mathbf y \approx \mathbf f(\boldsymbol\mu) + J(\mathbf x - \boldsymbol\mu). The deviation of y\mathbf y from f(μ)\mathbf f(\boldsymbol\mu) is then JJ times the deviation of x\mathbf x, and lesson 5-1's rule for a linear map gives
Py=E[J(x−μ)(x−μ)TJT]=J Px JT.P_y = E\big[J(\mathbf x - \boldsymbol\mu)(\mathbf x - \boldsymbol\mu)^{\mathsf T}J^{\mathsf T}\big] = J\,P_x\,J^{\mathsf T}.
The part it leaves out. Keep the next Taylor term for one output: y≈f(μ)+J(x−μ)+12(x−μ)T∇2f (x−μ)y \approx f(\boldsymbol\mu) + J(\mathbf x - \boldsymbol\mu) + \tfrac12(\mathbf x - \boldsymbol\mu)^{\mathsf T}\nabla^2f\,(\mathbf x - \boldsymbol\mu). The linear term averages to zero, but the quadratic one does not:
E[y]≈f(μ)+12∑i,j∂2f∂xi∂xjPij.E[y] \approx f(\boldsymbol\mu) + \tfrac12\sum_{i,j}\frac{\partial^2f}{\partial x_i\partial x_j}P_{ij}.
Linearisation drops this shift of the mean entirely.
The radar banana, exactly. Take the target due east at r=1000r = 1000 m with a bearing error ϵ\epsilon of standard deviation σb=15∘=0.2618\sigma_b = 15^\circ = 0.2618 rad. East is rcos⁡ϵr\cos\epsilon, and for a normal ϵ\epsilon the exact average is E[cos⁡ϵ]=e−σb2/2E[\cos\epsilon] = e^{-\sigma_b^2/2}:
E[east]=1000 e−0.26182/2=966.3 m,E[\text{east}] = 1000\,e^{-0.2618^2/2} = 966.3 \text{ m},
the million-sample answer (966.35 m), while linearisation keeps 1000 m. The second-order formula above gives 1000(1−σb2/2)=965.71000(1 - \sigma_b^2/2) = 965.7 m — most of the way there.
📐From radar measurement to position
A radar measures a target at range 10 km (σ = 10 m) and bearing 30° (σ = 0.1°). Position east =rsin⁡b= r\sin b, north =rcos⁡b= r\cos b, with Jacobian J=(sin⁡brcos⁡bcos⁡b−rsin⁡b)J = \begin{pmatrix} \sin b & r\cos b \\ \cos b & -r\sin b \end{pmatrix}. Then JPJTJPJ^{\mathsf T} gives σeast=15.92\sigma_{\text{east}} = 15.92 m, σnorth=12.29\sigma_{\text{north}} = 12.29 m and a correlation of −0.453-0.453: the errors are not independent any more. The eigenvectors show why: the ellipse has semi-axes 17.45 m across the line of sight (rσbr\sigma_b) and 10.00 m along it (σr\sigma_r). A Monte Carlo run of 200,000 samples agrees to the second decimal, and 95.03% of the samples fall inside the 95% ellipse.
💡When linearisation fails
Linearisation replaces f\mathbf f by its tangent plane at the mean. That is accurate only if f\mathbf f is nearly linear over the spread of the inputs. Give the radar a poor bearing — 15° of noise at 1000 m — and the true positions lie along an arc, a curved 'banana'. Monte Carlo with a million samples puts the mean at 966.35 m east, because points spread around an arc are on average closer than its radius; linearisation still says 1000 m, 33.7 m too far, and predicts an east–west spread of 2.00 m where the true one is 46.81 m. The covariance matrix is not wrong because of rounding or bad code; the method's assumption has failed.
The unscented transform. In 1997 Simon Julier and Jeffrey Uhlmann proposed a middle way between linearisation and Monte Carlo: choose a small set of sigma points that have exactly the mean and covariance of the input, push each one through the exact nonlinear function, and compute the weighted mean and covariance of the results. For nn inputs:
χ0=μ,χ±i=μ±((n+κ)P)i,w0=κn+κ,w±i=12(n+κ),\boldsymbol\chi_0 = \boldsymbol\mu, \qquad \boldsymbol\chi_{\pm i} = \boldsymbol\mu \pm \big(\sqrt{(n + \kappa)P}\big)_i, \qquad w_0 = \frac{\kappa}{n + \kappa}, \quad w_{\pm i} = \frac{1}{2(n + \kappa)},
where (M)i(\sqrt{M})_i is the ii-th column of a Cholesky factor (lesson 3-3). Five points for the radar give a mean of 966.31 m and an east–west spread of 47.68 m — against the million-sample answer of 966.35 m and 46.81 m. No derivatives are needed. The idea is a quadrature rule: the sigma points and weights integrate the function against a Gaussian, like the Gauss points of lesson 8-4.
Derivation
The sigma points reproduce the mean and covariance exactly
Let si\mathbf s_i be the columns of S=(n+κ)PS = \sqrt{(n + \kappa)P}, so that SST=∑isisiT=(n+κ)PSS^{\mathsf T} = \sum_i\mathbf s_i\mathbf s_i^{\mathsf T} = (n + \kappa)P.
Mean. The pairs μ±si\boldsymbol\mu \pm \mathbf s_i are symmetric about μ\boldsymbol\mu, and the weights add to 1 (w0+2n⋅12(n+κ)=κ+nn+κw_0 + 2n \cdot \frac{1}{2(n + \kappa)} = \frac{\kappa + n}{n + \kappa}):
∑w χ=w0μ+12(n+κ)∑i[(μ+si)+(μ−si)]=μ.\sum w\,\boldsymbol\chi = w_0\boldsymbol\mu + \frac{1}{2(n + \kappa)}\sum_i\big[(\boldsymbol\mu + \mathbf s_i) + (\boldsymbol\mu - \mathbf s_i)\big] = \boldsymbol\mu.
Covariance. The centre point contributes nothing (χ0−μ=0\boldsymbol\chi_0 - \boldsymbol\mu = \mathbf 0), and each pair contributes sisiT\mathbf s_i\mathbf s_i^{\mathsf T} twice:
∑w (χ−μ)(χ−μ)T=22(n+κ)∑isisiT=P.\sum w\,(\boldsymbol\chi - \boldsymbol\mu)(\boldsymbol\chi - \boldsymbol\mu)^{\mathsf T} = \frac{2}{2(n + \kappa)}\sum_i\mathbf s_i\mathbf s_i^{\mathsf T} = P.
So the 2n+12n + 1 points carry exactly the input's mean and covariance (and, by symmetry, zero skewness). Whatever the nonlinear function does to them, it does to a set of points that is statistically correct to second order — which is why the transform captures the banana's mean shift that linearisation misses.
📐Collision probability, and a paradox
At closest approach, project the relative position onto the plane perpendicular to the relative velocity. A collision means the true miss lies inside a small disc — here a combined radius of 10 m — around zero. With a predicted miss of 584 m and an isotropic position uncertainty σ\sigma, the probability is the integral of the 2-D Gaussian over that disc:
  • σ=50\sigma = 50 m: 9×10−329 \times 10^{-32} — a confident miss;
  • σ=413\sigma = 413 m: 1.1×10−41.1 \times 10^{-4} — the maximum, at σ=miss/2\sigma = \text{miss}/\sqrt2;
  • σ=1\sigma = 1 km: 4.2×10−54.2 \times 10^{-5} (Monte Carlo with 2×1072 \times 10^7 samples: 4.2×10−5±0.3×10−54.2 \times 10^{-5} \pm 0.3 \times 10^{-5});
  • σ=5\sigma = 5 km: 2.0×10−62.0 \times 10^{-6}.
Beyond the maximum, worse tracking makes the encounter look safer — the probability is diluted over a huge region. Operators therefore treat a low probability computed from a large covariance with suspicion. (These numbers are illustrative, not a reconstruction of the 2009 encounter.)
Derivation
Why the probability peaks at σ=d/2\sigma = d/\sqrt2
The collision disc (radius R=10R = 10 m) is tiny compared with the miss d=584d = 584 m, so the Gaussian density is nearly constant over it. The probability is then the disc's area times the density at its centre:
Pc≈πR2⋅12πσ2 e−d2/(2σ2)=R22σ2 e−d2/(2σ2).P_c \approx \pi R^2 \cdot \frac{1}{2\pi\sigma^2}\,e^{-d^2/(2\sigma^2)} = \frac{R^2}{2\sigma^2}\,e^{-d^2/(2\sigma^2)}.
Two effects compete as σ\sigma grows: the factor 1/σ21/\sigma^2 spreads the probability more thinly, while the exponential stops suppressing the distant disc. Write s=σ2s = \sigma^2 and differentiate:
dds(e−d2/(2s)s)=e−d2/(2s)s2(d22s−1)=0⟹s=d22,σ=d2=413 m.\frac{d}{ds}\Big(\frac{e^{-d^2/(2s)}}{s}\Big) = \frac{e^{-d^2/(2s)}}{s^2}\Big(\frac{d^2}{2s} - 1\Big) = 0 \quad\Longrightarrow\quad s = \frac{d^2}{2}, \quad \sigma = \frac{d}{\sqrt2} = 413 \text{ m}.
At the peak Pc=R2d2e−1=1005842×0.368=1.08×10−4P_c = \frac{R^2}{d^2}e^{-1} = \frac{100}{584^2} \times 0.368 = 1.08 \times 10^{-4} — the lab's maximum. Beyond it, more uncertainty means a smaller computed probability: dilution.
Given for one instant an intelligence which could comprehend all the forces by which nature is animated and the respective situation of the beings who compose it — an intelligence sufficiently vast to submit these data to analysis — … for it, nothing would be uncertain and the future, as the past, would be present to its eyes.
– Pierre-Simon Laplace, A Philosophical Essay on Probabilities (1814)
💡Where you will use this
Orbital Mechanics: conjunction assessment and covariance propagation along an orbit. GNSS Navigation: dilution of precision and error ellipses from the least-squares fix. Rocket Propulsion and Guidance: dispersion analysis of landing points and staging conditions. Geotechnical and Structural Engineering: first-order reliability methods (FORM) are linearised propagation; Monte Carlo and point-estimate methods are its sampled cousins. Every laboratory: measurement uncertainty budgets follow the same J P Jᵀ.
📚Further reading
  • Joint Committee for Guides in Metrology, Evaluation of measurement data — Guide to the expression of uncertainty in measurement (JCGM 100:2008, the 'GUM') — the international standard for JPJTJPJ^{\mathsf T} in laboratories; free from the BIPM.
  • Simon Julier and Jeffrey Uhlmann, 'Unscented filtering and nonlinear estimation', Proceedings of the IEEE 92(3), 401–422 (2004) — the unscented transform by its inventors.
  • Yaakov Bar-Shalom, X. Rong Li and Thiagalingam Kirubarajan, Estimation with Applications to Tracking and Navigation (Wiley, 2001) — covariance, error ellipses and consistency tests for engineers.
↓But there is a problem
Propagation tells you how uncertain a prediction is: carry the covariance through the model, and watch it grow. It cannot make the prediction any better. A navigation system has to do better. It has a model of how the vehicle moves, whose predictions grow more uncertain with every step, and measurements of where it is, which are noisy — and it must produce, every fraction of a second, the best estimate it can from both, together with an honest statement of how good that estimate is. In 1960 a Hungarian-born engineer, Rudolf Kálmán, published a recursive way to do exactly that. Later that year he visited NASA's Ames Research Center, where Stanley Schmidt's group was working out how a spacecraft could navigate to the Moon.