Uncertainty Propagation in C++

Uncertainty Propagation in C++

Evan Manning

Real-world measurements have uncertainties. Sometimes it's handy to keep track of those uncertainties through subsequent calculations.

Introduction

Uncertainty is the one great certainty in life, but how many simulations take this into account? Consider a spacecraft trajectory simulator. Such a program might start with an initial position and velocity for the spacecraft, and project it forward in time. For each timestep the program moves the spacecraft in accordance with its current velocity and adjusts the velocity to account for the gravitational attraction of significant bodies. If the program is implemented perfectly (neglecting roundoff errors), then for a given set of initial conditions it will be able to predict perfectly where the spacecraft will be at any future moment.

But in the real world there will always be some uncertainty in our knowledge of the spacecraft's initial position and velocity, as well as the positions and masses of planets. What we really need is a simulator that can track the effect of initial uncertainties and provide an assessment of the uncertainty in the final position. When a spacecraft is headed for a close encounter with an asteroid, mission planners will need to know more than whether the most likely path of the spacecraft impacts the asteroid. They will need to know the probability of a collision.

If the example program were coded in C, the options for adding uncertainty appraisal to this application would be: (1) ignore it, (2) try a number of plausible sets of inputs and look at the output, or (3) rewrite the application to keep an assessment of uncertainty alongside the expected value of each uncertain value. C++ adds a fourth option: use an UncertainDouble variable in place of each potentially uncertain double, and let the UncertainDouble class take care of all the bookkeeping.

This article describes a set of classes that enable numbers with associated uncertainties to behave like builtins. In many cases, the new types can be dropped into calculations in place of their old, "certain" counterparts. The result is a simulation that accounts for uncertainty, but maintains badly needed notational simplicity.

Implementation

The code presented here follows the Gaussian model of uncertainty (see sidebar) and is built around two private double data members: value (also called mean) and uncertainty (also called sigma). I call this class UDouble, for "uncertain double." Listing 1 is the class definition. For the most part, it looks like the definition of any mathematical class. mean and deviation member functions give read-only access to the data members. Other differences are discussed below.

Listing 2 shows UDouble's constructors and destructors. The constructors initialize the value and uncertainty data members; the destructor does nothing.

C++'s builtin double data type can be seen as the degenerate case of UDouble with an uncertainty equal to zero. In fact, when uncertainty is zero, members of this class do behave exactly like the builtin double type. For this reason, I've provided one constructor that accepts a single argument of type double. This constructor uses the argument as the value and defaults the uncertainty to 0.0. This is the default conversion from double to UDouble.

The uncertainty associated with a number tells us how many digits are significant, which allows us to print that number more intelligently than is usual with doubles. The function uncertain_print (Listing 3) takes a mean and a deviation (and an optional ostream) and prints out mean +/- deviation, carefully printing only as many digits of the mean as correspond to the first two digits of the deviation. So 1.2345 +/- 0.2387 prints as "1.23 +/- 0.24" and 0.012345 +/- 5.321 prints as "0.0 +/- 5.3". uncertain_print is used by UDouble's operator<< (not shown). I've implemented it outside of the UDouble class so that it can be used by other UDoubleXX classes. operator>> is much simpler, and is similarly implemented in terms of uncertain_read.

Propagation of Single Uncertainties

If x is 0.5 +/- 0.1, what is f(x)? We could make a good guess by looking at f(0.5) for the mean and then at f(0.4) and f(0.6) to estimate the deviation. The Gaussian approximation requires the value of f() at 0.5 (f(0.5)) and the slope of f() in the neighborhood of 0.5 (f'(0.5)) (The slope or derivative is the ratio of small changes in f(x) to small change in x.) The mean of f(x) is f() applied to the mean of x (<x>) and the deviation of f(x) is the deviation of x scaled by the slope of f() at the mean of x: f(x) = f(<x>) +/- f'(<x>)dx. This formula may look daunting but it is really quite simple to use when the slope of f() is known. The slope is known for all functions needed to make UDouble act like a double: unary +, unary -, acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod, frexp, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, and tanh.

Unary + and - have slopes of 1.0 and -1.0 respectively and so are easily implemented (Listing 4) . Slopes for the rest of these functions are calculable in terms of each other and a single value of x. For example, the slope of sin(x) is cos(x), the slope of exp(x) is exp(x), the slope of ceil(x) is 0, except where x is an integer, where the slope is infinite. (Listing 5) .

When the a function's derivative form is not known, the slope can be approximated by taking the difference between f(mean + sigma) and f(mean - sigma). Listing 6 presents the single-argument version of PropagateUncertaintiesBySlope, which does exactly this.

Multiple Sources of Uncertainty

Many operations combine two potentially uncertain inputs into one uncertain output. The simplest of these is binary "+": if a is 1.0 +/- 0.1 and b is 1.0 +/- 0.1 what is c = a + b? The value (or mean) of c follows familiar rules: <c> = <a> + <b> = 1.0 + 1.0 = 2.0. But the uncertainty of c will depend on whether or not a and b are correlated, that is, whether or not their uncertainties share a common source. I consider first two cases when a and b are completely correlated, then a case when they are not.

Adding Correlated Variables

Consider the case when b is a. This case might arise when we have two blocks known to be identical and we measure only one of them. Adding a + b is equivalent to asking how long two of the blocks are laid end-to-end, and of course, this is exactly twice the length of a single block: c = a + b = a + a = 2.0 * a = 2.0 * (1.0 +/- 0.1) = 2.0 +/- 0.2. Here the uncertainties of a and b have the same source and the same sign, so they simply add. You might think of them as parallel vectors. To carry out the same operation with an ideal uncertain double class, we would write

UDouble a(1.0, 0.1), b;
b = a;
cout << a << " + " << b << " = " << (a + b) << endl;

which prints:

1.00 +/- 0.10  +  1.00 +/- 0.10  =  2.00 +/- 0.20

For negative correlation, consider the case where b is (2 - a). This is somewhat harder to imagine, but perhaps we have a box known to be two meters long and two blocks that together fill it perfectly. Even though our measurement of block a is imperfect we know that b is (2 - a). In this case c = a + b = a + (2 - a) = 2 + (a - a) = 2.0 +/- 0.0 Since the uncertainties of a and b have a common source but opposite sense they subtract. You might think of them as parallel vectors pointing in opposite directions (antiparallel vectors). This case might be coded as follows:

UDouble a(1.0, 0.1), b;
b = 2.0 - a;
cout << a << " + " << b << " = " << (a + b) << endl;

which prints:

1.00 +/- 0.10  +  1.00 +/- 0.10  =
2.0000000 +/- 0.0

UDouble is a class template, which allows it to expand to two almost identical classes, parameterized by the variable is_correlated. The int parameter is_correlated is conceptually a boolean, but I didn't use the new boolean type because it is not yet widely available. typedefs hide all this complexity from the users. Users simply work with UDoubleCorr and UDoubleUncorr. For class UDoubleCorr uncertainties add simply; for UDoubleUncorr uncertainties add by hypotenuse, as I show below. Anti-correlated uncertainties can cancel when added, so the private uncertainty data member can be negative in UDoubleUncorr. Listing 7 shows operator+= and binary operator+.

Adding Uncorrelated Variables

The uncorrelated case arises when the uncertainties in a and b come from independent sources. In this case the uncertainty vectors would be (on average) at right angles and their sum would be their hypotenuse, or the square root of the sum of their squares: sqrt(0.1² + 0.1²) = 0.1 * sqrt(2) ~= 0.14 so c = a + b = 2.00 +/- 0.14. This case is an application of the uncorrelated version of the UDouble<is_correlated> class, UDoubleUncorr. It can be coded as follows:

UDoubleUncorr a(1.0, 0.1), b(1.0, 0.1);
cout << a << " + " << b << " = " << (a + b) << endl;

which prints:

>1.00 +/- 0.10 + 1.00 +/- 0.10 = 2.00 +/- 0.14

Adding Partially Correlated Variables

The most complicated kind of uncertainty addition occurs when the two operands are partially correlated. If a and b are independent and c is their sum, as in the previous case, then a and c are neither perfectly correlated nor perfectly independent. a is correlated with the a portion of c but uncorrelated with the b portion of c. The code presented here cannot trace such partial correlations, but more advanced methods can do so. One such case might be coded:

UDouble a(1.0, 0.1), b(1.0, 0.1), c;
//make "c" be half correlated with "a" and half with "b"
//renormalized to be 1.00 +/- 0.10
c = (a + b) / sqrt(2.0) + 1.0 - sqrt(2.0);
cout << a <<& quot; + " << c << " = " << (a + c) << endl;

which prints:

1.00 +/- 0.10  +  1.00 +/- 0.10  = 2.00 +/- 0.18

All binary functions use slope to figure out the uncertainty from each source and then add the two uncertainties either simply (if correlated) or as hypotenuse (not correlated). The operators which work this way are +=, -=, *=, /=, binary +, -, *, and /, and fmod, atan2, and pow. Some examples are given in Listings 7 & 8.

Listing 9 shows how uncertainties are propagated by slope through an unknown function of two uncertain variables.

The only operations defined for builtin type double which UDouble lacks are casting to other numerical types and relational operators. I've left these operations out because UDoubles are conceptually multi-valued. What should

UDoubleUncorr = ud(100.0, 3.0);
int i = ud;

yield? ud is 100.0 +/- 3.0 and so is likely to be near 100, but could well be anywhere from 90 to 110, and in theory might be -10,230. The most sensible single value is the mean, and if this is wanted it is available through int i = ud.mean();. Similarly, there is no single answer to the question of whether 100.0 +/- 3.0 is greater than 101.0 +/- 6.0. It is possible to assign a probability to this value but if that probability is expressed as a simple floating-point number between 0.0 and 1.0, then expressions like if (uda > udb) will almost always evaluate as true. If a comparison of means is wanted then the appropriate expression is if (uda.mean() > udb.mean()).

Demo Program

These classes, as well as many others which I describe briefly later are available from the JPL ftp site: airs1.jpl.nasa.gov in the directory pub/evan/c++, or by writing to me at [email protected]. The ftp site also also contains a program that puts UDouble through its paces and explains the results. This demo program is implemented mostly in terms of another class, UDoubleTest, which is composed of one private UDoubleUncorr member and one UDoubleCorr member, and distributes most operations to those classes.

Using UDouble

An application that uses UDouble must include header udms.h (also available from JPL) and must change all double variables that can hold uncertain values to UDouble. The application can input data via operator>> if the data is formatted as "mean +/- sigma"; otherwise you'll need custom input routines. Output should be formatted correctly by the overloaded operator<< without modification.

Use the class UDoubleCorr in cases with only one source of uncertainty. You can use UDoubleUncorr when there are multiple uncertainties and each independent uncertainty gets mixed with other uncertainties exactly once. Many applications will fit neither of these sets of restrictions, and so will need the more advanced classes in this collection.

This package contains multiple and very different implementations of the same functionality. This collection of classes now contains a UDoubleTest class that has members of the other UDoubleXX types and distributes all operations to the members. The multiple distribution makes it is easy to compare the output of various methods and see at a glance where they differ.

In the context of the complete set of uncertainty classes, the classes presented here are known as UDoubleMS, for uncertain double mean-sigma.

UDoubleMS offers the simplest possible model of uncertainty. The other UDouble classes on the ftp site offer more accurate, but more computationally expensive, solutions to the problem of modeling uncertainties. The correlation tracking class, UDoubleCT, uses the same underlying Gaussian model as the UDoubleMS class, but keeps track of uncertainties from multiple sources correctly. The ensemble class, UDoubleEnsemble, does not depend on the Gaussian model but instead models each uncertain variable with an ensemble of possible values.

I hope someday to add another class which uses some knowledge of the second derivative of functions (curve) to improve accuracy and (optionally) to warn when curves begin to break down the applicability of the Gaussian model. It would also warn when discontinuities threaten the applicability of this model.

Conclusion

I like to think of these classes as adding intelligence to an application. A carefully designed application using UDouble is not only making the basic calculations that a more primitive application would make but also "thinking" about the accuracy of its results. With Gaussian breakdown checking, the application could even check the accuracy of the first-order check on accuracy.

Acknowledgements

This work was partly done at the Jet Propulsion Laboratory, California Institute of Technology, under a contract from the National Aeronautics and Space Administration.

Evan Manning has degrees in Applied Physics from Caltech and Stanford. He has been working as a self-taught C programmer in defense and space applications for the past 8 years. Currently he works for Telos Information Systems as a consultant at NASA's Jet Propulsion Laboratory. He can be reached at [email protected].