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
that “it is therefore important that uncertainties be small.” As far as I can
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
— I create an array that stores 1024 samples drawn from a normal distribution.
Yep. To do math on it, I just use
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
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
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,
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.