# Implementing error propagation

Is there a software package that supports, or could support, through operator overloading or extensions, code such as the following;

x = rand_arr(10) ; array 10 elements long
y = rand_arr(10)
z = x + y ; elemental addition (z[0]=x[0]+y[0];z[1]=x[1]+y[1],...)
print z

If the above, or something similar, is valid in this hypothetical language, does it also support, or could it support via operator overloading / custom classes / extensions / whatever, the following additional syntax:

x.err[*] = 0.1 % a fixed error
y.err = rand_arr(10) % a dynamic noisy error range
z = x + y
print z.err

Or whatever the implementation dictates, but the point being that the error is tracked for me through the z = x+y without work on my part.

I could imagine default support as per a standard (or customizable) function for +-/*e^ and log() operators, and support for extensions or hooks to add support for additional functions (system, and user) such as z = smooth(x,3) and z = my_func(x^2)

This would be a nice feature for a language.

## Error Explanation

There is a discussion of ERROR PROPAGATION here.

For the error tracking, I'm assuming that x.err is the error associated with each measurement, in this case, a constant value of 0.1. y.err is similar, but each measurement is part of some noise. Given a standard definition of error propagation, such as the following:

When adding x + y with errors dx and dy, the solution has error dx + dy

Or in code

z = x + y => dz = dx + dy

I can add lines of code to track dz, but is there a language that supports, or could support, dz without me doing extra work?

• On the SO version of this question (which I can no longer find), you had a link to an "error propagation" site. Could you provide that here?
– whuber
Dec 21 '10 at 18:46
• Found the SO version! stackoverflow.com/questions/4483266/…
– whuber
Dec 21 '10 at 19:03
• Copied SO version; I'll ask mmyers to remove it from SO.
– user88
Dec 21 '10 at 20:03
• K​a​b​l​a​m​m​o. Dec 21 '10 at 20:39

See Wikipedia's List of uncertainty propagation software, in particular Python uncertainties .
There's even a conference: http://probabilistic-programming.org/wiki/NIPS*2008_Workshop

(Beware -- random variables that propagate through a network and converge and correlate get messy.)

It's difficult to determine what you're asking for because you haven't specified the semantics of "x.err" etc., but it sounds like you might be interested in interval arithmetic. Implementations are available in Lisp, Numerica, Maple, Matlab, and Mathematica. Libraries are available for ADA, Fortran, C++, etc.

### Edit

In light of comments (scattered between here and SO), it appears the OP is asking about implementing automatic error propagation. In principle that's no harder than implementing, say, a complex number class or interval arithmetic, which are straightforward exercises in any OO system. The data structure would be a tuple $(x, \epsilon)$ where $x$ represents a real number and $\epsilon \ge 0$ quantifies its "error". The usual numbers would be embedded into this structure via $x \to (x, 0)$. The semantics would include

\eqalign{ (x, \epsilon) + (y, \delta) = &(x+y, (\epsilon^2 + \delta^2)^{1/2}) \cr (x, \epsilon) - (y, \delta) = &(x-y, (\epsilon^2 + \delta^2)^{1/2}) \cr (x, \epsilon) \times (y, \delta) = &(x y,((y \epsilon)^2 + (x \delta)^2)^{1/2}) \cr \text{etc.} \cr f((x, \epsilon), (y, \delta)) = &(f(x, y), || (\frac{\partial f}{\partial x}\epsilon, \frac{\partial f}{\partial y}\delta)||). \cr }

(The last generalizes everything preceding it and points the way to implementing powers, exponentials, logs, trig functions, etc.)

Such an implementation could truly be handy for those who understand deeply what's going on. For others--who might not be aware of or appreciate the importance of assumptions of independence, of relatively small errors, of differentiability of functions, etc.--it would be truly dangerous. As a simple example, when $(x, \epsilon)$ and $(y, \delta)$ are perfectly negatively correlated (implying $\epsilon = \delta$), the correct operation is $(x, \epsilon) + (y, \delta) = (x + y, 0)$.

• @whuber Well, it was just "Is there a programming language with built-in overloadable operators so I could implement small error propagation?", so I have moved it to StackOverflow and that's why it was deleted here -- also the OP edited it on SO, so I think it should stay deleted here.
– user88
Dec 20 '10 at 0:28
• @mbq Fair enough. Although the question was not well formulated, between it and its title it appeared there was interest in error propagation, which would seem to fit with the interests of a stats site.
– whuber
Dec 20 '10 at 15:17
• @mbq A clarification posted on the SO version makes it abundantly clear that the OP is about error propagation. There is no longer any question that this is relevant to stats. Apparently there's no longer any option for us to vote to migrate it back, but it definitely belongs here and should be re-opened.
– whuber
Dec 21 '10 at 16:12
• @whuber I still think it is more about operator overloading than about error propagation (not to mention what I think about error propagation), but you have it opened again.
– user88
Dec 21 '10 at 18:41
• @whuber and others... Thank you. Sorry for the poorly worded question/title. I debated here v. SO, as the question is both about error propagation AND software implementation of it. Anyway, this is an excellent answer. Thanks. Dec 21 '10 at 20:44

I recommend you look into R, the open source software for statistical computing. R supports many vectorized functions to deal with element by element computation of objects. Chapter 3 of The R Inferno provides a good outline on vectorization in R. From your post above, I'm assuming you have some background in computer programming, and R's command line interface will probably feel right at home to you.

Here's your pseudo code above implemented in R:

set.seed(42)    #Makes random numbers reproducible.
x <- sample(1:20, size = 10, replace = FALSE)
y <- sample(1:20, size = 10, replace = FALSE)
z <- x + y
z
> z
[1] 29 32 23 20 19 23 35  4 23 19

I'm not entirely sure what you are trying to do with the error coding, but R has the capability to support just about any mathematical or statistical computation you can dream up. If someone else hasn't already implemented a solution for your problem, you can leverage the fact R is a programming language and write your own code to solve it. It is relatively easy to interface with other languages such as C, C++, python from R so there are plenty of options there.

• MASS, the book (and package on CRAN) has a function doing error propagation. Jan 22 '16 at 11:32

The error propagation mentioned by @whuber in his edit to this answer is implemented in Maple (disclosure: I work for them) as well, in addition to the interval arithmetic he also mentions. The package that implements it is called ScientificErrorAnalysis (here is the help page). It might take some getting used to. An example:

with(ScientificErrorAnalysis):
a := Quantity(2, 0.1);
#      Quantity(2, 0.1)
b := Quantity(3, 0.2);
#      Quantity(3, 0.2)
c := Quantity(2, 0.1);
#      Quantity(2, 0.1)
a + b;
#      Quantity(2, 0.1) + Quantity(3, 0.2)
combine(a + b, errors);
#      Quantity(5., 0.2236067977)
combine(a * b, errors);
#      Quantity(6., 0.5000000000)
combine(log(a), errors);
#      Quantity(0.6931471806, 0.05000000000)
combine(a^2+a, errors);
#      Quantity(6., 0.5000000000)
combine(a^2+c, errors);
#      Quantity(6., 0.4123105626)

Basically, the idea is you do your computation with these Quantity objects and at the very end you run combine(..., errors) to get a single Quantity back. I won't claim that it's perfect, but it's not extremely naive: note that, while a and c are defined in the same way, a^2+a and a^2+c have different errors. This is because a^2 and a are correlated whereas a^2 and c are not.