<p><strong>Yes</strong>, 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.</p>
<h2>Example 1 with Parallax (<a href="https://en.wikipedia.org/wiki/Normal_distribution" rel="nofollow noreferrer">Normal</a> distribution to <a href="https://en.wikipedia.org/wiki/Inverse_distribution#Reciprocal_normal_distribution" rel="nofollow noreferrer">Reciprocal Normal</a> distribution)</h2>
<p>Suppose we measure the parallax angle of a star as <span class="math-container">$0.02\pm0.001$</span> 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, <a href="https://en.wikipedia.org/wiki/Parallax#Stellar_parallax" rel="nofollow noreferrer">wikipedia</a>:</p>
<p><a href="https://i.sstatic.net/VSA3b.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/VSA3b.png" alt="enter image description here" /></a></p>
<p>Then our estimate of distance to the star is <span class="math-container">$d=1/p = 50$</span> 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 <span class="math-container">$0.021$</span> and <span class="math-container">$0.019$</span> arcseconds, respectively. This gives us <span class="math-container">$47.619$</span> and <span class="math-container">$52.632$</span> parsecs.  So our distance estimate is <span class="math-container">$50^{+2.632}_{-2.381}$</span> parsecs.  Notice the upper bound is greater than the lower bound.</p>
<p>We can verify this result with a Monte Carlo run where we generate random samples from a normal distribution with mean <span class="math-container">$0.02$</span> and standard deviation of <span class="math-container">$0.001$</span>.  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 <span class="math-container">$50$</span>, with upper and lower bounds that each contain <span class="math-container">$34$</span>% of the samples between the bounds and the median also converging to our previously calculated bounds of <span class="math-container">$47.619$</span> and <span class="math-container">$52.632$</span> 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.</p>
<p><strong>Notes:</strong></p>
<ol>
<li>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.</li>
<li>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.</li>
<li>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 (<a href="https://en.wikipedia.org/wiki/Lutz%E2%80%93Kelker_bias" rel="nofollow noreferrer">Lutz-Kelker bias</a>).  For outstanding further reading, see: <a href="http://spiff.rit.edu/classes/phys440/lectures/lutz_kelker/lutz_kelker.html" rel="nofollow noreferrer">http://spiff.rit.edu/classes/phys440/lectures/lutz_kelker/lutz_kelker.html</a>.</li>
</ol>
<h2>Example 2 with Asteroid Orbital Inclination (<a href="https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution" rel="nofollow noreferrer">von-Mises Fisher</a> distribution to <a href="https://en.wikipedia.org/wiki/Rayleigh_distribution" rel="nofollow noreferrer">Rayleigh</a> distribution)</h2>
<p>Let's move on to a more complex example.  Swike asked the question: <a href="https://astronomy.stackexchange.com/questions/35471/why-are-asteroids-with-zero-orbital-inclination-rare/40518#40518">Why are asteroids with zero orbital inclination rare?</a>.</p>
<p>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 <span class="math-container">$p=3$</span> in <span class="math-container">$\mathbb{R}^3$</span>. I generated a 1000 3D points with a Monte Carlo run.</p>
<p><a href="https://i.sstatic.net/Gdo46.gif" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Gdo46.gif" alt="enter image description here" /></a></p>
<p>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.</p>
<p><a href="https://i.sstatic.net/P0b6B.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/P0b6B.png" alt="enter image description here" /></a></p>
<p>You could calculate asymmetric error bounds on either the input (3D) or output (1D) points.</p>

CC BY-SA 4.0
2021-11-22T23:39:48.943Z
2021-11-22T23:39:48.943Z
Connor Garcia
<p>I'm not sure this question really belongs here, but you mention the word &quot;astronomy&quot;, and I'm an astronomer and I have an opinion on how to <em>add</em> numbers with asymmetric uncertainties:</p>
<h3>Inherent ambiguity in a solution</h3>
<p>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.</p>
<p>However, if you're only handed out values such as <span class="math-container">$5^{+3}_{-2}$</span> and <span class="math-container">$2^{+2}_{-1}$</span>, there are infinitely many distributions which can be described by three such numbers (i.e. a central value <span class="math-container">$x_0$</span>, an upper error <span class="math-container">$\sigma_{+}$</span>, and a lower error <span class="math-container">$\sigma_{-}$</span>), and hence no &quot;correct&quot; way to add them, neither analytically nor numerically.</p>
<p>Additionally, you should be aware what the &quot;central value&quot; actually means; does it represent the <em>mean</em>, the <em>mode</em>, or the <em>median</em>? All three examples are common in the literature.</p>
<p>As an example, consider the distribution below.</p>
<p><a href="https://i.sstatic.net/mq2WY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/mq2WY.png" alt="mean-mode-median" /></a></p>
<p><em>Talking</em> about a distribution, I think many consider the mean to be the &quot;characteristic&quot; value. But <em>looking</em> 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
<span class="math-container">$3.1_{-1.4}^{+1.4}$</span> (mean),
<span class="math-container">$2.1_{-0.5}^{+2.4}$</span> (mode), or
<span class="math-container">$2.7_{-1.1}^{+1.8}$</span> (median).</p>
<h3>Asymmetric errors are not Gaussian</h3>
<p>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,
<span class="math-container">$$
5^{+3}_{-2} + 2^{+2}_{-1} \ne 7^{+\sqrt{3^2+2^2}}_{-\sqrt{2^2+1^2}}.
$$</span></p>
<p>That this is wrong, can be seen from the <a href="https://en.wikipedia.org/wiki/Central_limit_theorem" rel="nofollow noreferrer">central limit theorem</a>: 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.</p>
<h3>An approximate solution</h3>
<p>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 <em>sense</em>, then at least <em>more</em> sense than adding in quadrature.</p>
<p>One method, described by <a href="https://arxiv.org/abs/physics/0306138" rel="nofollow noreferrer">Barlow (2003)</a>, 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 <em>typically</em> much better than just adding in quadrature in the &quot;usual&quot; way.</p>
<p>Based on this approach, I wrote an easy-to-use Python function <code>add_asym</code>, <a href="https://github.com/anisotropela/add_asym" rel="nofollow noreferrer">available on GitHub</a> that does this, with the option of two different transformations. It is described and tested in detail in <code>&lt;blatant self-promotion&gt;</code> <a href="https://ui.adsabs.harvard.edu/abs/2019A%26A...627A..84L" rel="nofollow noreferrer">Laursen et al. (2019)</a> (Appendix B) <code>&lt;/blatant self-promotion&gt;</code>.</p>
<p>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 <span class="math-container">${x_0}_{-\sigma_{-}}^{+\sigma_{+}} = 0_{-2}^{+3}$</span>.</p>
<p>Specifically, from the bottom and up, the <span class="math-container">$\color{blue}{\mathsf{blue}}$</span> curves show the addition of <a href="https://en.wikipedia.org/wiki/Skew_normal_distribution" rel="nofollow noreferrer">skewed Gaussians</a>, <a href="https://en.wikipedia.org/wiki/Weibull_distribution" rel="nofollow noreferrer">Weibull</a>, <a href="https://en.wikipedia.org/wiki/Log-normal_distribution" rel="nofollow noreferrer">lognormal</a>, <a href="https://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution" rel="nofollow noreferrer">Fréchet</a>, and <a href="https://en.wikipedia.org/wiki/Log-logistic_distribution" rel="nofollow noreferrer">loglogistic</a> distributions.</p>
<p><a href="https://i.sstatic.net/REiKK.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/REiKK.png" alt="testCLT" /></a></p>
<p>The two different transformations (<span class="math-container">$\color{green}{\mathsf{green}}$</span> and <span class="math-container">$\color{olive}{\mathsf{olive\text{-}ish}}$</span>) are seen to be 1) quite consistent with each other, and 2) an acceptable match to (most of) the &quot;true&quot; results.</p>
<p>For comparison, the &quot;usual&quot;, but wrong, method of adding errors separately in quadrature is shown in <span class="math-container">$\color{red}{\mathsf{red}}$</span>, never departing from <span class="math-container">$\sigma_{+}/\sigma_{-}=1.5$</span>.</p>
<h3>Relation to astronomy… -ish</h3>
<p>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.</p>
<p>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 <em>logarithm</em> 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.</p>

CC BY-SA 4.0
2023-07-11T16:50:09.557Z
What needs to be done if we have to use division or multiplication instead?  Can we just use the logarithm and proceed as outlined for addition.
Christian Herenz
CC BY-SA 4.0
2023-07-12T08:56:09.163Z
@ChristianHerenz I'm a really bad statistician, but Peter "The Frequentist" Jakobsen will be visiting later today, so I'll ask him about his thoughts. But regarding using log, my gut feeling is that it's fine in principle, but the error (i.e. error compared to the unobtainable correct solution) will be correspondingly larger.
pela
CC BY-SA 4.0
2023-07-12T13:34:10.373Z
So PJ shuddered when I mentioned asymmetric uncertainties and had no real solution, but also suggested going to log space.
pela
CC BY-SA 4.0
2021-11-18T17:58:28.613Z
2022-11-17T11:00:21.273Z
2022-11-17T11:00:21.273Z
pela
<p>I've used <code>uncertainties</code>'s <code>ufloat</code>package to propagate uncertainties as I thought it's the standard way to do it. Now <strong>I want to propagate errors when I have different lower and upper uncertainties</strong>.</p>
<p>From <a href="https://newton.cx/%7Epeter/2013/04/propagating-uncertainties-the-lazy-and-absurd-way/" rel="noreferrer">this source</a> I understand that</p>
<blockquote>
<p>there’s the uncertainties package which is quite sophisticated and seems to be the <strong>gold standard</strong> for this kind of thing.</p>
</blockquote>
<p>but</p>
<blockquote>
<p>analytic propagation still fundamentally operates by ignoring nonlinear terms, which means, in the words of the uncertainties documentation, that <strong>“it is therefore important that uncertainties be small.”</strong></p>
</blockquote>
<blockquote>
<p>...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.</p>
</blockquote>
<p>which is in the same line as <a href="https://physics.stackexchange.com/questions/343151/error-propagation-with-asymmetric-uncertainties">this question</a> in Physics.SE; the same source states that</p>
<blockquote>
<p>if your method can’t handle large uncertainties, it’s pretty useless for astronomy.</p>
</blockquote>
<p>and I do agree with the author. From this I guess <a href="https://en.wikipedia.org/wiki/Propagation_of_uncertainty" rel="noreferrer">known</a> expressions (from first order linear expansions) such as</p>
<p><span class="math-container">$$\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}$$</span></p>
<p>won't work well.</p>
<p>Also there's <a href="https://astronomy.stackexchange.com/questions/38371/how-can-i-calculate-the-uncertainties-in-magnitude-like-the-cds-does/38372#38372">this question</a> in Astronomy.SE in which using simple error propagation won't work... so</p>
<p><strong>Question</strong>: how do you manage this kind of problems when treating with asymmetric, or/and large uncertainties? What would be the &quot;gold standard&quot; for this kind of thing? Should I follow <code>Floris</code>'s answer?:</p>
<blockquote>
<p>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.</p>
</blockquote>
<p>do you have any example to understand this in deep?</p>

<data-analysis>
How do you propagate asymmetric errors? (on the practical way to...)
nuwe
# How do you propagate asymmetric errors? (on the practical way to...)

(Asked by: username_0 on 2021-11-18)

I've used `uncertainties`'s `ufloat`package 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](https://newton.cx/%7Epeter/2013/04/propagating-uncertainties-the-lazy-and-absurd-way/) 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](https://physics.stackexchange.com/questions/343151/error-propagation-with-asymmetric-uncertainties) 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](https://en.wikipedia.org/wiki/Propagation_of_uncertainty) 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](https://astronomy.stackexchange.com/questions/38371/how-can-i-calculate-the-uncertainties-in-magnitude-like-the-cds-does/38372#38372) 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?

**Tags:** data-analysis

---

## 2 answers

### Answer by username_2 on 2021-11-18. Score: 17 | Accepted answer

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.

[![mean-mode-median](https://i.sstatic.net/mq2WY.png)](https://i.sstatic.net/mq2WY.png)

*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](https://en.wikipedia.org/wiki/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)](https://arxiv.org/abs/physics/0306138), 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](https://github.com/anisotropela/add_asym) that does this, with the option of two different transformations. It is described and tested in detail in  [Laursen et al. (2019)](https://ui.adsabs.harvard.edu/abs/2019A%26A...627A..84L) (Appendix B) .

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](https://en.wikipedia.org/wiki/Skew_normal_distribution), [Weibull](https://en.wikipedia.org/wiki/Weibull_distribution), [lognormal](https://en.wikipedia.org/wiki/Log-normal_distribution), [Fréchet](https://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution), and [loglogistic](https://en.wikipedia.org/wiki/Log-logistic_distribution) distributions.

[![testCLT](https://i.sstatic.net/REiKK.png)](https://i.sstatic.net/REiKK.png)

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.

> **Comments:**
> **username_3**: What needs to be done if we have to use division or multiplication instead?  Can we just use the logarithm and proceed as outlined for addition.  
> **username_2**: @ChristianHerenz I'm a really bad statistician, but Peter "The Frequentist" Jakobsen will be visiting later today, so I'll ask him about his thoughts. But regarding using log, my gut feeling is that it's fine in principle, but the error (i.e. error compared to the unobtainable correct solution) will be correspondingly larger.  
> **username_2**: So PJ shuddered when I mentioned asymmetric uncertainties and had no real solution, but also suggested going to log space.  

---

### Answer by username_1 on 2021-11-22. Score: 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](https://en.wikipedia.org/wiki/Normal_distribution) distribution to [Reciprocal Normal](https://en.wikipedia.org/wiki/Inverse_distribution#Reciprocal_normal_distribution) 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](https://en.wikipedia.org/wiki/Parallax#Stellar_parallax):

[![enter image description here](https://i.sstatic.net/VSA3b.png)](https://i.sstatic.net/VSA3b.png)

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](https://en.wikipedia.org/wiki/Lutz%E2%80%93Kelker_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](https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution) distribution to [Rayleigh](https://en.wikipedia.org/wiki/Rayleigh_distribution) distribution)

Let's move on to a more complex example. Swike asked the question: [Why are asteroids with zero orbital inclination rare?](https://astronomy.stackexchange.com/questions/35471/why-are-asteroids-with-zero-orbital-inclination-rare/40518#40518).

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.

[![enter image description here](https://i.sstatic.net/Gdo46.gif)](https://i.sstatic.net/Gdo46.gif)

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.

[![enter image description here](https://i.sstatic.net/P0b6B.png)](https://i.sstatic.net/P0b6B.png)

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


---

