How do you propagate asymmetric errors? (on the practical way to...)

Question

I've used uncertainties's ufloatpackage to propagate uncertainties as I thought it's the standard way to do it. Now I want to propagate errors when I have different lower and upper uncertainties.

From this source I understand that

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

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.”

...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.

which is in the same line as this question in Physics.SE; the same source states that

if your method can’t handle large uncertainties, it’s pretty useless for astronomy.

and I do agree with the author. From this I guess known expressions (from first order linear expansions) such as

$$\delta f(x_1,..., x_n) \approx \sqrt{\left(\frac{\partial f}{\partial x_1}\right)^2 \delta x_1 ^2 + ... \left(\frac{\partial f}{\partial x_n}\right)^2 \delta x_n ^2 } = \sqrt{\sum\limits_{i=1}^n \left(\frac{\partial f}{\partial x_i} \delta x_i \right)^2}$$

won't work well.

Also there's this question in Astronomy.SE in which using simple error propagation won't work... so

Question: how do you manage this kind of problems when treating with asymmetric, or/and large uncertainties? What would be the "gold standard" for this kind of thing? Should I follow Floris's answer?:

In a Monte Carlo simulation, you sample the distributions of your input variables, and transform them according to the formula; you can then plot the resulting output distribution, and compute its shape etc.

do you have any example to understand this in deep?

Answer 1

I'm not sure this question really belongs here, but you mention the word "astronomy", and I'm an astronomer and I have an opinion on how to add numbers with asymmetric uncertainties:

Inherent ambiguity in a solution

If you have the full probability distribution functions (PDFs) of two asymmetric distributions, they should be added by a convolution. This will give you an exact result if you can do it analytically, or — in case this isn't possible — a numerical, almost-exact result. The Monte Carlo solution that you mention belongs to the latter.

However, if you're only handed out values such as $5^{+3}_{-2}$ and $2^{+2}_{-1}$, there are infinitely many distributions which can be described by three such numbers (i.e. a central value $x_0$, an upper error $\sigma_{+}$, and a lower error $\sigma_{-}$), and hence no "correct" way to add them, neither analytically nor numerically.

Additionally, you should be aware what the "central value" actually means; does it represent the mean, the mode, or the median? All three examples are common in the literature.

As an example, consider the distribution below.

Talking about a distribution, I think many consider the mean to be the "characteristic" value. But looking at a distribution like this, I think many will think of the peak as the characteristic value describing this distribution, i.e. the mode. On the other hand, if a confidence interval (CI) is defined as the area between the 16th and the 84th percentile (as for a normal distribution), only the median is guaranteed to lie inside the CI. Depending on your preferences, you could then describe the PDF in this example as $3.1_{-1.4}^{+1.4}$ (mean), $2.1_{-0.5}^{+2.4}$ (mode), or $2.7_{-1.1}^{+1.8}$ (median).

Asymmetric errors are not Gaussian

Nevertheless, whatever you do, don't add upper and lower errors in quadrature. Despite being common in the literature, this approach has no statistical justifications, and is only correct for normal (Gaussian) distributions, which are not asymmetric. That is, $$ 5^{+3}_{-2} + 2^{+2}_{-1} \ne 7^{+\sqrt{3^2+2^2}}_{-\sqrt{2^2+1^2}}. $$

That this is wrong, can be seen from the central limit theorem: In the limit of many distributions of the same asymmetry, the combined PDF should approach a Gaussian distribution. In contrast, errors added in quadrature never decrease in asymmetry.

An approximate solution

However, it turns out that, even if you don't know the full PDFs, there are ways too add numbers that do make… if not sense, then at least more sense than adding in quadrature.

One method, described by Barlow (2003), circumvents the problem by first transforming the asymmetric numbers to symmetric numbers in a clever way, then adding these in quadrature, then transforming back. Exactly how to do this transformation is, again, ambiguous, but it turns out that, for a wide range of different functional forms of asymmetric PDFs, this method is typically much better than just adding in quadrature in the "usual" way.

Based on this approach, I wrote an easy-to-use Python function add_asym, available on GitHub that does this, with the option of two different transformations. It is described and tested in detail in <blatant self-promotion> Laursen et al. (2019) (Appendix B) </blatant self-promotion>.

The following figure (from my paper) shows how the ratio between the upper and lower error declines, when correctly adding (i.e. convoluting) a large number (up to 25) of diverse PDFs, which are all described by a median ${x_0}_{-\sigma_{-}}^{+\sigma_{+}} = 0_{-2}^{+3}$.

Specifically, from the bottom and up, the $\color{blue}{\mathsf{blue}}$ curves show the addition of skewed Gaussians, Weibull, lognormal, Fréchet, and loglogistic distributions.

The two different transformations ($\color{green}{\mathsf{green}}$ and $\color{olive}{\mathsf{olive\text{-}ish}}$) are seen to be 1) quite consistent with each other, and 2) an acceptable match to (most of) the "true" results.

For comparison, the "usual", but wrong, method of adding errors separately in quadrature is shown in $\color{red}{\mathsf{red}}$, never departing from $\sigma_{+}/\sigma_{-}=1.5$.

Relation to astronomy… -ish

So far I've only been talking mathematics. The reason I wrote the code mentioned above is that I needed to add asymmetric uncertainties in my work with galaxies. Galaxy masses, like galaxy luminosities, galaxy radii, gas cloud sizes, city populations, forest areas, and many, many other entities in nature, are distributed highly asymmetrically.

Various scaling relations tend to make astrophysical quantities (and other quantities in nature) well described by power laws, and thus it is often better to consider the distribution of the logarithm of some quantity, although there's not always a physical justification to do this. For this reason I prefer referring to a PDF using the median, with the 16th and the 84th percentiles as the lower and upper limits for the CI, because this is same in linear and logarithmic space.

Answer 2

Yes, Floris's answer is the gold standard when an analytic statistical approach isn't apparent. Here are two astronomy examples of how to do Monte Carlo runs to characterize your resulting distributions.

Example 1 with Parallax (Normal distribution to Reciprocal Normal distribution)

Suppose we measure the parallax angle of a star as $0.02\pm0.001$ arcseconds. Further suppose the errors on our measurement of parallax are normally distributed and expressed as one-sigma standard deviation from the mean. Image credit, wikipedia:

Then our estimate of distance to the star is $d=1/p = 50$ parsecs. But what is the error in our distance estimate? If we take the upper and lower bounds for the angle measurement errors, we have $0.021$ and $0.019$ arcseconds, respectively. This gives us $47.619$ and $52.632$ parsecs. So our distance estimate is $50^{+2.632}_{-2.381}$ parsecs. Notice the upper bound is greater than the lower bound.

We can verify this result with a Monte Carlo run where we generate random samples from a normal distribution with mean $0.02$ and standard deviation of $0.001$. Each one of these samples in arcseconds can be translated to distance with the above formula. You will see the median of these distance samples converge to $50$, with upper and lower bounds that each contain $34$% of the samples between the bounds and the median also converging to our previously calculated bounds of $47.619$ and $52.632$ parsecs. This is the median, 16th, and 84th percentile distribution that is Pela's favorite in the other answer. Note that the mean of the distance sample will creep up to 50.14 parsecs. As Pela mentions, asymmetric errors occur only in non-Gaussian distributions in which we don't expect the mean and median to coincide.

Notes:

1. This example is for a single variable propagating symmetric error bounds through a transformation that results in asymmetric errors bounds. Since our distance as a function of parallax is monotonic, we can pass our upper and lower bounds through the transformation to get the new bounds.
2. For multi-variate inputs and outputs, the technique is the same. We generate inputs according to their known distributions, calculate the outputs according to our functions, and then perform statistical analysis of the outputs to characterize their distributions.
3. Here, parallax is an excellent example for Monte Carlo runs, because large uncertainties in parallax will lead to large asymmetric bounds in distance estimates. We could even end up with a bias contributing to the belief that stars are closer than they actually are (Lutz-Kelker bias). For outstanding further reading, see: http://spiff.rit.edu/classes/phys440/lectures/lutz_kelker/lutz_kelker.html.

Example 2 with Asteroid Orbital Inclination (von-Mises Fisher distribution to Rayleigh distribution)

Let's move on to a more complex example. Swike asked the question: Why are asteroids with zero orbital inclination rare?.

In my answer, I assume that the axes of revolution of the main belt asteroids are normally distributed with highest concentration perpendicular to the ecliptic. Then their intersection with the Unit Sphere could be expressed as a von-Mises Fisher distribution with $p=3$ in $\mathbb{R}^3$. I generated a 1000 3D points with a Monte Carlo run.

Then I calculated their angle (inclination) from the reference plane normal. The result looks like a Rayleigh distribution, similar to the graphic Swike posted in his original question.

You could calculate asymmetric error bounds on either the input (3D) or output (1D) points.