Recipe 18.14. Doing Arithmetic with Error PropagationCredit: Mario Hilgemeier ProblemYou 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. SolutionThe 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) DiscussionHere 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 AlsoLibrary Reference and Python in a Nutshell docs for module math. |