Recipe18.14.Doing Arithmetic with Error Propagation


Recipe 18.14. Doing Arithmetic with Error Propagation

Credit: Mario Hilgemeier

Problem

You have numbers coming from measurements affected by known percentual uncertainties, and you want to perform arithmetic on these numbers while tracking the uncertainty in the results.

Solution

The simplest approach is to code a class that implements arithmetic operators applying the classical physicists' error-propagation formulas:

import math class Measurement(object):     ''' models a measurement with % uncertainty, provides arithmetic '''     def _ _init_ _(self, val, perc):         self.val = val                            # central value         self.perc = perc                          # % uncertainty         self.abs = self.val * self.perc / 100.0   # absolute error     def _ _repr_ _(self):         return "Measurement(%r, %r)" % (self.val, self.perc)     def _ _str_ _(self):         return "%g+-%g%%" % (self.val, self.perc)     def _addition_result(self, result, other_abs):         new_perc = 100.0 * (math.hypot(self.abs, other_abs) / result)         return Measurement(result, new_perc)     def _ _add_ _(self, other):         result = self.val + other.val         return self._addition_result(result, other.abs)     def _ _sub_ _(self, other):         result = self.val - other.val         return self._addition_result(result, other.abs)     def _multiplication_result(self, result, other_perc):         new_perc = math.hypot(self.perc, other_perc)         return Measurement(result, new_perc)     def _ _mul_ _(self, other):         result = self.val * other.val         return self._multiplication_result(result, other.perc)     def _ _div_ _(self, other):         result = self.val / other.val         return self._multiplication_result(result, other.perc)

Discussion

Here is a typical example of use for this Measurement class:

m1 = Measurement(100.0, 5.5)   # measured value of 100.0 with 5.5% error m2 = Measurement(50, 2)        # measured value of 50.0 with 2% error print "m1 = ", m1 print "m2 = ", m2 print "m1 + m2 = ", m1 + m2 print "m1 - m2 = ", m1 - m2 print "m1 * m2 = ", m1 * m2 print "m1 / m2 = ", m1 / m2 print "(m1+m2) * (m1-m2) = ", (m1+m2) * (m1-m2) print "(m1-m2) / (m1+m2) = ", (m1-m2) / (m1+m2) # emits: # m1 =  100+-5.5% # m2 =  50+-2% # m1 + m2 =  150+-3.72678% # m1 - m2 =  50+-11.1803% # m1 * m2 =  5000+-5.85235% # m1 / m2 =  2+-5.85235% # (m1+m2) * (m1-m2) =  7500+-11.7851% # (m1-m2) / (m1+m2) =  0.333333+-11.7851%

What is commonly known as a percentage error is of course best described as a percentage of uncertainty. For example, when we state that some measured quantity is 100 with an error of 5.5% (or, equivalently, 5.5%), we mean that we know, with a reasonable level of confidence, that the quantity lies somewhere between 94.5 and 105.5. The error-propagation formulas are somewhat heuristic, rather than rigorous, but they're quite traditional and have proven over the centuries that they perform acceptably in most large computations in physics or engineering.

Class Measurement, as shown in this recipe, does not support arithmetic with floatsonly arithmetic between instances of Measurement. For those rare occasions when I need, in such computations, numbers that are known "exactly", it is easiest to input them as "measurements with an error of 0%". For example, if I have measured some sphere's radius as 1 meter +- 3%, I can compute the sphere's volume (with the well-known formula, 4/3 pi times the cube of the radius) as follows:

r = Measurement(1, 3) v = Measurement(4/3.0*math.pi, 0) * r * r * r print v # emits: 4.18879+-5.19615%

Avoiding accidental operations with floats that are presumed to be exact, but in fact are not, is quite helpful: this way, when I need to state that a certain number has 0 error, I'm reminded to consider whether things are truly that way. If your applications are different, so that you do need operations between measurements and exact floats all over the place, you can insert, as the first line of every one of the arithmetic special methods, the following statement:

        if isinstance(other, float):             other = Measurement(other, 0)

Alternatively, you could perform this coercion in a special method named _ _coerce_ _, but that approach is considered obsolete and is discouraged in modern Python. If you do perform the coercion in the various arithmetic special methods (_ _add_ _, _ _sub_ _, etc.), don't forget to also add the _ _radd_ _, etc, equivalentsafter all, if you want to be able to code:

    some_measurement * 2.0

you will no doubt also want to be able to code:

    2.0 * some_measurement

and get exactly the same effects. For this purpose, in Python, your class needs to define the various _ _r... versions of the operator special methods. However, I'm not pursuing this line of reasoning further, because in most cases, you will be best served by not having the implicit ability to do arithmetic in an automatic way between measurements and floatsmuch like, with Python 2.4's decimal module, you can't implicitly do arithmetic in an automatic way between decimal numbers and floats.

See Also

Library Reference and Python in a Nutshell docs for module math.



Python Cookbook
Python Cookbook
ISBN: 0596007973
EAN: 2147483647
Year: 2004
Pages: 420

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net