# Uncertainty propagation using MonteCarloMeasurements gives weird results
**URL:** https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745
**Category:** Specific Domains
**Tags:** monte-carlo, measurements
**Created:** [January 13, 2025, 10:21pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745 "2025-01-13T22:21:44Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)
#### Post date: [January 13, 2025, 10:21pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/1 "2025-01-13T22:21:44Z")
Hi!
I am not familiar with MonteCarloMeasurements.jl, so I am probably doing something dumb, but I thought that it would be the good tool for what I want to do and I am getting some weird results.
I am trying to calculate a diffusion coefficient using an Arrhenius law and I simply want to propagate my uncertaities. As it involves an exponential in the expression, I thought it would be the occasion to try MonteCarloMeasurements.jl, as I think Measurements.jl should not be able to be applied here due to the high non-linearity with an exponential. I also want to combine it with Unitful, which should be fine.
I am basically calculating the following expression: D = D\_0 \exp(-\frac{E\_a}{RT}), where D\_0 and E\_a can have some uncertainty.
Here is what I tried and what I obtained, from a simple MWE:
```julia
using MonteCarloMeasurements
using Unitful
D0 = (2.5e-12±0)u"m^2/s"
Ea = (227±62)u"kJ/mol"
D = D0*exp(-Ea/(8.3145u"J/mol/K"*1100u"K"))
@show D
# D = 4.57403079910306e-16 m² s⁻¹ ± 1.6590917723694656e-14 m² s⁻¹
D02 = 2.5e-12±0
Ea2 = (227*1e3)±62000
D2 = D02*exp(-Ea2/((8.3145±0)*(1100±0)))
@show D2
# D2 = 4.57e-16 ± 1.7e-14
D03 = 2.5e-12u"m^2/s"
Ea3 = 227u"kJ/mol"
D3 = D03*exp(-Ea3/(8.3145u"J/mol/K"*1100u"K"))
@show D3
# D3 = 4.1578487110411763e-23 m² s⁻¹
D04 = 2.5e-12
Ea4 = (227*1e3)
D4 = D04*exp(-Ea4/(8.3145*1100))
@show D4
# D4 = 4.1578487110411763e-23
```
The correct results should be on the order of e-23 as obtained without MonteCarloMeasurements.jl. I assume I misunderstood something in the usage of this package because not only the uncertainty is probably wrong, but the main calculation is too. Playing with Measurements.jl gives the correct result, but a probably wrong uncertainty. I tried to play with the number of particles, but that didn’t change the results. Looking forward to understand what I missed!
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [January 13, 2025, 11:14pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/2 "2025-01-13T23:14:22Z")
This is all fine … the distribution just gets rather skewed due to the exponential and thus the mean (shown for uncertain values) is nowhere near the median (corresponding to a certain value):
```julia-repl
julia> D02 = 2.5e-12±0
2.5e-12 Particles{Float64, 2000}
julia> Ea2 = (227*1e3)±62000
227000.0 ± 62000.0 Particles{Float64, 2000}
julia> D2 = D02*exp(-Ea2/((8.3145±0)*(1100±0)))
4.57403e-16 ± 1.66e-14 Particles{Float64, 2000}
julia> @show D2
D2 = 4.57e-16 ± 1.7e-14
4.57403e-16 ± 1.66e-14 Particles{Float64, 2000}
julia> D2.particles |> mean
4.5740307991030655e-16
julia> D2.particles |> median
4.157886227774504e-23
```
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 14, 2025, 4:26am UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/3 "2025-01-14T04:26:37Z")
Hi there! 👋
I think the result is correct. As pointed out by @bertschi , the resulting distribution is highly skewed, here’s a histogram (note the logarithmic y-axis)
You can endow the particles with a nominal value like this, here I used the mean
```julia
D0 = with_nominal((2.5e-12±0)u"m^2/s", (2.5e-12)u"m^2/s")
Ea = with_nominal((227±62)u"kJ/mol", 227u"kJ/mol")
D = D0*exp(-Ea/(8.3145u"J/mol/K"*1100u"K"))
```
Then you can extract the nominal value of the result like this
```julia
nominal(D)
```
```julia
julia> nominal(D)
4.1578487110411763e-23 m² s⁻¹
```
and you see that it matches what you compute without uncertainty. With such a nonlinear function resulting in such a skewed distribution, propagating the mean of the expectation does not give you the expectation of the mean (see, e.g., [Jensen’s inequality](https://en.wikipedia.org/wiki/Jensen%27s_inequality)).
> [@Iddingsite](#):
>
> Playing with [Measurements.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Measurements) gives the correct result
probably not, it gives you the propagated mean of the input (like `nominal` above), which _is not_ the correct mean of the output 😉
Here’s a histogram of the output with only 10% of the uncertainty of above, in this plot it’s easier to see the much less extreme distribution
---
### Author:  [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)
#### Post date: [January 14, 2025, 10:19am UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/4 "2025-01-14T10:19:25Z")
Hi to you both, thx for your answers!
I understand more what’s happening behind the hood now and why I got this result. And I can definitely see the appeal of using MonteCarloMeasurements.jl here, as I can actually see what’s happening with the distributions and the propagations. It is really cool compared to Measurements.jl where you can’t really get a grasp of things for neophytes likes me and it will just give a result.
Now, to come back to my initial problem, I see that I can retrieve the correct value with `nominal` (which I don’t really understand how it is calculated, but fair enough), but what would be the way to estimate my uncertainty in this case? I am not a pro in statistics, but as my distribution is not normal, I guess I could calculate a median absolute deviation? I don’t really see the gain, apart from a better understanding of my problem and the knowledge that I can’t use linear error propagation theory for this equation, using MonteCarloMeasurements.jl here.
On the other hand, I could rewrite D = D\_0 \exp(-\frac{E\_a}{RT}) as log(D) = log(D\_0) - \frac{E\_a}{RT} to propagate the uncertainties, which gives something I can report:
```julia
using MonteCarloMeasurements
using Unitful
using Plots
D0 = (2.5e-12±0)u"m^2/s"
Ea = (227±62)u"kJ/mol"
D_log = (log(ustrip.(D0)) - Ea/(8.3145u"J/mol/K"*1100u"K"))
plot(D_log)
```
Which output
which looks good to me because it is normal. And Measurements.jl gives the same result. So I can get an uncertainty for log(D) but not really for D I assume?
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 14, 2025, 3:01pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/5 "2025-01-14T15:01:30Z")
> [@Iddingsite](#):
>
> I see that I can retrieve the correct value with `nominal`
The point of uncertainty propagation is that the _correct_ value is _unknown_. You are talking about the _expectation of the mean_, but the probability that this is the _correct_ value is actually 0.
> [@Iddingsite](#):
>
> but what would be the way to estimate my uncertainty in this case?
What information do you want to know? The standard deviation is one measure of uncertainty, even if the distribution is skewed. You can also compute some quantile of the output distribution if you prefer that, or several quantiles. Since the \log is normally distributed in this case, you can of course think in terms of this quantity like you suggest.
---
### Author:  [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)
#### Post date: [January 14, 2025, 5:05pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/6 "2025-01-14T17:05:06Z")
Thx for your answer.
> [@baggepinnen](#):
>
> The point of uncertainty propagation is that the _correct_ value is _unknown_. You are talking about the _expectation of the mean_, but the probability that this is the _correct_ value is actually 0.
I understand your point with the “correct” value. That makes sense.
> [@baggepinnen](#):
>
> What information do you want to know? The standard deviation is one measure of uncertainty, even if the distribution is skewed. You can also compute some quantile of the output distribution if you prefer that, or several quantiles. Since the \log is normally distributed in this case, you can of course think in terms of this quantity like you suggest.
Ok, what I am getting confused is what the number after \pm means here, in the case of D. So for instance
```julia
D = 4.57403079910306e-16 m² s⁻¹ ± 1.6590917723694656e-14 m² s⁻¹
```
which has a nominal value of
```julia
nominal(D) = 4.1578487110411763e-23 m² s⁻¹
```
My input in my uncertainty (`Ea = (227±62)u"kJ/mol"`) is for instance 2 \sigma. So, I guess ± 1.6590917723694656e-14 m² s⁻¹ is also twice the standard deviation? Which is huge compared to the order of magnitude of the nominal value but I guess it can be expected with this big uncertainty. I just want to make sure that I understand what this value represents and then I think I will have everything in my hands to take a decision and close this topic!
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 14, 2025, 5:15pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/7 "2025-01-14T17:15:01Z")
[The convention for \pm in MCM.jl is to work with standard deviations](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/api/#MonteCarloMeasurements.:%C2%B1), so \mu \pm \sigma creates a normally distributed number with mean \mu and standard deviation \sigma, this is also how the the uncertain numbers are displayed. See also the functions `pmean, pstd, pvar` etc. in the docs
- [Home · MonteCarloMeasurements Documentation](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/#Basic-Examples-1)
- [Home · MonteCarloMeasurements Documentation](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/#Constructors-1)
---
### Author:  [@Iddingsite](https://discourse.julialang.org/u/Iddingsite)
#### Post date: [January 14, 2025, 6:35pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/8 "2025-01-14T18:35:22Z")
Thank you for your response. So I believe this is an important distinction with Measurements.jl, as I believe it doesn’t matter in this case if we are using 1 \sigma or 2 \sigma as long as it is consistent.
I understand better the philosophy of the package now. Thanks a lot for your help and your patience!
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 14, 2025, 6:39pm UTC](https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745/9 "2025-01-14T18:39:47Z")
Measurements will linearize the function, in which case it doesn’t matter which multiple of standard deviation you use. When the function is nonlinear, this is very important!