Python Cookbook 2Nd Edition Jun 1002005 [Electronic resources] نسخه متنی

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.