Propagating Uncertainties: The Lazy and Absurd Way

2013 April 3

I needed to write some code that does calculations and propagates
uncertainties under a fairly generic set of conditions. A well-explored
problem, surely? And indeed, in Python there’s the
`uncertainties`

package which is
quite sophisticated and seems to be the gold standard for this kind of thing.

Being eminently reasonable, `uncertainties`

represents uncertain variables
with a mean and standard deviation, propagating errors analytically. It does
this quite robustly,
computing the needed derivatives magically,
but analytic propagation still fundamentally operates by ignoring nonlinear
terms, which means, in the words of the
`uncertainties`

documentation,
that “it is therefore important that uncertainties be small.” As far as I can
tell, `uncertainties`

does analytic propagation as well as anything out there,
but honestly, if your method can’t handle large uncertainties, it’s pretty
useless for astronomy.

Well, if analytic error propagation doesn’t work, I guess we have to do it empirically.

So I wrote a little Python module. To
represent 5 ± 3 I don’t create a variable that stores `mean=5`

and `stddev=3`

— I create an array that stores 1024 samples drawn from a normal distribution.
Yep. To do math on it, I just use `numpy`

‘s
vectorized operations. When I report a result, I look at the 16th, 50th, and
84th percentile points of the resulting distribution.

Ridiculous? Yes. Inefficient? Oh yes. Effective? *Also* yes, in many cases.

For instance: the `uncertainties`

package doesn’t support asymmetric error
bars or upper limits. My understanding is that these could be implemented, but
they badly break the assumptions of analytic error propagation — an asymmetric
error bar *by definition* cannot be represented by a simple mean and standard
deviation, and an upper limit measurement *by definition* has a large
uncertainty compared to its best value. But I can do math with these values
simply by drawing my 1024 samples from the right distribution —
skew normal or
uniform between zero and the limit. I can mix perfectly-known values,
“standard” (i.e. normally-distributed) uncertain values, upper limits, and
anything else, and everything Just Works. (It might be hard to define the
“uncertainty” on a complex function of a mixture of all of these, but that’s
because it’s genuinely poorly-defined — analytic propagation is just
misleading you!)

Another example: `uncertainties`

spends a lot of effort tracking correlations,
so that if `x = 5 ± 3`

, then `x - x = 0`

precisely, not 0 ± 4.2. My approach
gets this for free.

I’ve found that approaching uncertainties this way helps clarify your thinking
too. You worry: is 1024 samples big enough? Well, did you actually measure 5 ±
3 by taking 1024 samples? Probably not. As
Boss Hogg points out,
the uncertainties on your uncertainties are large.
I’m pretty sure that only in *extreme* circumstances would the number of
samples actually limit your ability to understand your uncertainties.

Likewise: what if you’re trying to compute `log(x)`

for `x = 9 ± 3`

? With 1024
samples, you’ll quite likely end up trying to take the logarithm of a negative
number. Well, *that’s telling you something*. In many such cases, `x`

is
something like a luminosity, and while you might not be confident that it’s
much larger than zero, I can *guarantee* you it’s not actually less than zero.
The assumption that x is drawn from a normal distribution is failing. Now,
living in the real world, you want to try to handle these corner cases, but if
they happen persistently, you’re being told that the breakdown of the
assumption is a significant enough effect that you need to figure out what to
do about it.

Now obviously this approach has some severe drawbacks. But it was *super* easy
to implement and has Just Worked remarkably well. Those are big deals.