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:
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!
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):
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)
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).
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
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:
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 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?
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.
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.
I understand your point with the ācorrectā value. That makes sense.
Ok, what I am getting confused is what the number after \pm means here, in the case of D. So for instance
D = 4.57403079910306e-16 mĀ² sā»Ā¹ Ā± 1.6590917723694656e-14 mĀ² sā»Ā¹
which has a nominal value of
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!
The convention for \pm in MCM.jl is to work with standard deviations, 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
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!
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!
