MATLAB: Uncertainty calculation for given equation

uncertainty

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,

Best Answer

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.

Related Solutions

MATLAB: Generating a population randomly starting from some parameters

This is sometimes known as statistical tolerancing, to identify the distribution of a function of some random variable, or a set of random variables.

Monte Carlo simulation is indeed one way you can solve the problem, at least approximately so. Thus, just evaluate the function Q at all of those sets of values, then plot a histogram, etc.

There are other methods you can use to compute or approximate the distribution of Q. One can use Taylor series approximations, or even Taguchi methods. (I wrote several papers on the topic of modified Taguchi methods for this exact purpose. Find one of them here in Technometrics .)

A problem is that it is often true that if your function Q is nonlinear, that no common distribution exists to describe the distribution of Q. In those cases, you might use a member of the Pearson or Johnson families of distributions to yield an approximate distribution.

Edit: I pointed out some ideas in my comment below. So lets look at what happens for a simple nonlinear function. How about Q(x) = x^2?

Q = @(x) x.^2;

Really, how simple can you get, and still be nonlinear?

Now, suppose that x is normally distributed, with mean 0 and variance 1. I think that anyone will see, even without any computation that Q(x) cannot look normally distributed, since Q(x) can NEVER take on negative values. It turns out that in this case, Q(X) will have a simple distribution, a Chi-squared random variable, with one degree of freedom.

ezplot(@(x) chi2pdf(x,1),[0,10])
grid on

Sometimes you get lucky. As you can see though, this distribution looks quite non-normmally distributed. The reason is Q(x) has the property that around zero, we cannot simply truncate the Taylor series of x^2 to only the linear term and below. That quadratic terms is hugely significant.

We can try to convince ourselves that I got the distribution right by a simple Monte Carlo simulation.

hist(Q(randn(10000,1)),100)

Yep. Looks close enough for government work to me.

But suppose we changed x slightly? Suppose that x was normally distributed, with mean 100, and variance 1? In fact, the true distribution here should be a non-central chi-square, but I'm feeling too lazy to generate that one. A Monte Carlo simulation will suffice.

hist(Q(100 + randn(100000,1)),100)

Locally, it turns out that the Taylor series of Q(x), expanded around the mean of x (i.e., x==100), is now pretty well approximated by a linear polynomial, so dropping the quadratic term does not hurt us too much. That is, we can write

Q(u + 100) = 10000 + 200*u + u^2

Thus, for moderately small values of u, we can truncate the polynomial to ignore the u^2 term. Q is indeed sufficiently locally linear out there when u is a Normal random variable with mean zero and variance 1. (Then u+100 will clearly have mean 100 and variance 1.)

There is much more we can do with the techniques of statistical tolerancing. I've just brushed the surface here. The multidimensional case, where Q is a function of several random variables is really no more complicated.

MATLAB: Sum of two random variables with different distributions

I'm sorry. There is no command in MATLaB that will give you the CDF of the sum of two general random variables. This is for good reason: there is NO simple way to write the CDF of the sum of two general, unrelated random variables, with arbitrary distributions.

There are many things we might wish to do that have no simple solutions. That is not to say this is impossible. In fact, there are several classical solution approaches. The general idea is called statistical tolerancing. I'll suggest a couple of ideas that you can use.

1. You can do a Monte Carlo simulation. Generate random samples from each component, then form the sum. You can then compute a sample CDF from the data points. Lots and lots of points here will yield a decent approximation to the CDF.

2. Compute the mean, variance, skewness, kurtosis, etc., of the sum. There are many ways this can be done, using Taylor series approximations, or using various Taguchi style methods. For example, as Walter points out, the mean of a sum is simple to compute. However, the higher order moments will take a little more work. Given those moments, then you can choose a member of a distribution family that matches those moments as well as possible. The usual families of distributions that are used are the Pearson and the Johnson families.

The above scheme (#2) has its flaws for two distributions that are so different from each other, especially when one of them is bounded. That suggests that no simple Pearson or Johnson member will be a good approximant. So I would suggest #1 as the best approach, to form a completely empirical CDF.

If you do choose option #2, then I would still strongly suggest using a Monte Carlo scheme to validate those results. This is especially true because a Monte Carlo simulation is so easy to do.

Best Answer

Related Solutions

MATLAB: Generating a population randomly starting from some parameters

MATLAB: Sum of two random variables with different distributions

Related Question