Uncertainty calculation for given equation

Hi,

In MATLAB R2017a, is it possible to set uncertainties on variables (e.g. a=5±0.2), and then perform calculations with these variables and find the final uncertainty?

 a=5±0.2
 b=3±0.15
 c=19±0.01

 y=(a+b)/c^2

and find the uncertainty in y (e.g. y=x±0.2)?

I would like to know if MATLAB can do this automatically, and if it is possible to set the uncertainties. I have symbolic math toolbox and curve fitting toolbox available only.

Thank you for your help,

[Sterling Baird]

While not a MATLAB solution (hence a comment, not an answer) and mostly for other people who come across this, if you're looking for an analytic solution, Mathematica has this built-in. See Numbers with Uncertainty. You can even use a series expansion to order "n" (see AroundReplace). I think MATLAB FEX has some similar uncertainty implementations, at least for basic uncertainty propagation, and automatic differentiation (AD) may also be worth looking into.

Answer 1

[John D'Errico]

No. It is not possible to do in any automatic fashion. At least, not with any tools that I know of that are publically available. This is a problem known as statistical tolerancing. There are at least three methods that I can think of to solve the problem:

1. Monte Carlo simulation
2. Truncated Taylor series approximations
3. Taguchi style methods, including what we named modified Taguchi methods

1. Simulate the system using random number generation for the variables.

2. Approximate your function using a first order truncated Taylor series approximation. It is now linear, so the problem reduces to a linear combination of random variables with known distributions. If your random variables are assumed to be Gaussian (Normally distributed), then the solution is trivial.

This scheme can also be extended to a second order Taylor series approximation, although the computations get slightly more involved.

3. Taguchi methods reduce to a numerical integration. Modified Taguchi methods (as described by me and a colleague in a few papers) were shown to be equivalent to a Gauss-Hermite numerical integration. They allow you to compute approximations to the moments of the function. You can then use those estimated moments to infer the distribution of the response from its mean, variance, skewness, and Kurtosis.

You can look online under the phrase "statistical tolerancing modified Taguchi', and you should find some papers that describe the techniques.

https://www.jstor.org/stable/1269802?seq=1#page_scan_tab_contents

Taguchi methods even have a Wiki page. I suppose that were I feeling sufficiently industrious, I should add a page for the modified methods that we developed.

In fact, very long ago I did write some tools that do statistical tolerancing. In fact, it is a credit to MATLAB that my roughly 20 year old code worked fine. Sorry, but I cannot put the code into the public domain. Not really worth the effort anyway, since the solution is so easy without it.

Regardless, if I assume that your variables a,b,c are all uniformly distributed, it really is not that hard to do:

N = 1e8;
a = 5 + (rand(N,1)*2-1)*0.2;
b = 3 + (rand(N,1)*2-1)*0.15;
c = 19 + (rand(N,1)*2-1)*0.01;

y = (a+b)./(c.^2);

mean(y)
ans =
        0.0221605698565819

std(y)
ans =
       0.00040004350708596

skewness(y)
ans =
       0.00018202561148938

kurtosis(y)
ans =
           2.3545148004182

In fact, those estimates of the moments are pretty accurate. Using more accurate methods (fairly high order Gauss-Legendre numerical integration) they can be computed as:

      mean: 0.022160670958632
       std: 0.000400054003691638
  skewness: 0.00012265988869227
  kurtosis: 2.35442773957793

I'd trust the first 8 or 10 of those digits as being correct. So the Monte Carlo estimation was quite reasonable.

hist(y,1000)

See from the histogram that the result does not take on any simple distribution. Of course, the computations above were done using the presumption that your variables were uniformly distributed with the indicated limits. If you intended them to be Normally distributed instead, the results would be slightly different.

Answer 2

[Arpad Lenart]

Thank you for your reply. Your reply was very helpful. I think I will be using the Monte Carlo simulation for now, since I am just learning how to use MATLAB. So an accuracy of 8 digits will be more than enough.