Welcome, guest | Sign In | My Account | Store | Cart

This class permits full error propagation for numeric values. It wraps a value and an associated error (standard deviation, measurement uncertainty...). Numeric operations are overloaded and permit use with other Uncertain objects, precise values without errors, generally speakin gmost other numeric (even array) objects. The "traitification" can easily be reverted by dumping all references to traits.

Python, 90 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#!/usr/bin/python
# -*- coding: utf8 -*-
# Uncertain quantities.
# (c) Robert Jordens <jordens@debian.org>
# Made available freely under the Python license

import numpy as np
from enthought.traits.api import HasStrictTraits
from enthought.traits.api import Str, Float, Bool, List, Dict, Enum

class Uncertain(HasStrictTraits):
    """
    Represents a numeric value with a known small uncertainty 
    (error, standard deviation...).
    Numeric operators are overloaded to work with other Uncertain or
    numeric objects.
    The uncertainty (error) must be small. Otherwise the linearization
    employed here becomes wrong.
    The usage of traits can easily be dumped.
    """
    value = Float
    error = Float(0.)
    def __init__(self, value=0., error=0., *a, **t):
        self.value = value
        self.error = abs(error)
        super(Uncertain, self).__init__(*a, **t)
    def __str__(self):
        return "%g+-%g" % (self.value, self.error)
    def __repr__(self):
        return "Uncertain(%s, %s)" % (self.value, self.error)
    def __float__(self):
        return self.value
    def assign(self, other):
        if isinstance(other, Uncertain):
            self.value = other.value
            self.error = other.error
        else:
            self.value = other
            self.error = 0.
    def __abs__(self):
        return Uncertain(abs(self.value), self.error)
    def __add__(self, other):
        if isinstance(other, Uncertain):
            v = self.value + other.value
            e = (self.error**2+other.error**2)**.5
            return Uncertain(v, e)
        else:
            return Uncertain(self.value+other, self.error)
    def __radd__(self, other):
        return self + other # __add__
    def __sub__(self, other):
        return self + (-other) # other.__neg__ and __add__
    def __rsub__(self, other):
        return -self + other # __neg__ and __add__
    def __mul__(self, other):
        if isinstance(other, Uncertain):
            v = self.value * other.value
            e = ((self.error*other.value)**2+(other.error*self.value)**2)**.5
            return Uncertain(v, e)
        else:
            return Uncertain(self.value*other,
                    self.error*other)
    def __rmul__(self, other):
        return self * other # __mul__
    def __neg__(self):
        return self*-1 # __mul__
    def __pos__(self):
        return self
    def __div__(self, other):
        return self*(1./other) # other.__div__ and __mul__
    def __rdiv__(self, other):
        return (self/other)**-1. # __pow__ and __div__
    def __pow__(self, other):
        if isinstance(other, Uncertain):
            v = self.value**other.value
            e = ((self.error*other.value*self.value**(other.value-1.))**2+
                (other.error*np.log(self.value)*self.value**other.value)**2)**.5
            return Uncertain(v, e)
        else:
            return Uncertain(self.value**other,
                    self.error*other*self.value**(other-1))
    def __rpow__(self, other):
        assert not isinstance(other, Uncertain)
            # otherwise other.__pow__ would have been called
        return Uncertain(other**self.value,
                self.error*np.log(other)*other**self.value)
    def exp(self):
        return np.e**self
    def log(self):
        return Uncertain(np.log(self.value), self.error/self.value)
  • error propagation formulas are probably correct
  • error vallues are in units of the value and not in percent
  • may operations are delegated to other simpler or known operations
  • IMO this class is a significant improvement over the classes presented in:

7 comments

Eric-Olivier LE BIGOT 14 years, 10 months ago  # | flag

I'd like to suggest my uncertainties module, which has many advantages over the above recipe:

  1. It gives a correct result for x-x, when x has an uncertainty: the correct result is exactly zero, contrary to what the above recipe gives. More generally, all the correlations between the various variables of a calculation are taken into account.

  2. uncertainties.py allows many more mathematical functions to be used (almost all the function from the standard math module), including matrix inversion with Numpy. Logical functions can still be used as well (x == y), etc.

Eric-Olivier LE BIGOT 14 years, 10 months ago  # | flag

… and I'd like to add that the uncertainties module does not depend on any other module, which makes it easier to deploy!

jordens+activestate Robert Jordens 14 years, 9 months ago  # | flag

Eric, you are absolutely right. Correlations (like in x-x) make the above module fail miserably. On the other hand people tend to forget that correlations also frequently appear hidden in e.g. the measuring process. Then all approaches that ignore the provenance of the quantities fail ;-)

I'll propose one other approach: symbolic calculation with sympy. Short and concise:

from sympy import *
from uncertain import Uncertain  

def stdev(term, *errors):
    return sqrt(reduce(Add, (
        (term.diff(var)*Symbol("sigma_"+var.name))**2
        for var in term.atoms(Symbol) if var.name in errors), 0.)) 

def stdev_for_uncertain(term, **vars):
    st =  stdev(term, *vars.keys())
    subs = {}
    for k,v in vars.iteritems():
        subs[Symbol(k)] = v.value
        subs[Symbol("sigma_"+k)] = v.error
    return st.n(subs=subs, maxprec=10) 

if __name__ == "__main__":
    x,y = symbols("xy")
    t = (x**2+y**2)/(x*y)
    print stdev(t, "x", "y")
    print stdev_for_uncertain(t, x=Uncertain(3,.1), y=Uncertain(5,.5))

That'll be as good as your module with regards to correlations. It can even do crazy stuff like integrals...

jordens+activestate Robert Jordens 14 years, 9 months ago  # | flag

Then I'd also like to propose a __str__ method for all uncertain implementations. It correctly truncates and rounds the representations of value error and also applies to large and small quantities that need exponential notation.

import numpy as np
def __str__(self):
    # number of significant figures to print after the error-digit
    pr = 0
    er, va = self.error, self.value
    if er == 0. or not np.isfinite(va):
        # exact
        return "%g" % va
    if not np.isfinite(er):
        return "%g(%g)" % (va, er)
    # 10^nsig < er < 10^(nsig+1)
    nsig = int(np.floor(np.log10(er).min()))
    # 10^vsig < |va| < 10^(vsig+1)
    if np.fabs(va) > 0:
        vsig = int(np.floor(np.log10(np.fabs(va)).max()))
    else:
        vsig = 0
    # is er or va reuire exp noation: do so and append es
    tsig = max(nsig, vsig)
    if abs(tsig) > 4:
        er /= 10**tsig
        nsig -= tsig
        va /= 10**tsig
        vsig -= tsig
        es = "e%d" % tsig
    else:
        es = ""
    # now the larger of the two is between 1e-5 and 1e5
    # valid: numer of digits that should appear befor the error
    valid = 1+pr+max(0, vsig-nsig)
    # digits after the comma
    nshown = pr-min(0, nsig)
    # keep as many digits as error has
    er = round(er, pr-nsig)
    va = round(va, pr-nsig)
    s = "%0.*f(%.0f)%s" % (nshown, va, er*10**nshown, es)
    return s
jordens+activestate Robert Jordens 14 years, 9 months ago  # | flag

A snippet from the testcases:

def test_str1(self):
    self.assertEqual(str(Uncertain(0.002534, 0.0)), "0.002534")
def test_str2(self):
    self.assertEqual(str(Uncertain(5, 2)), "5(2)")
def test_str4(self):
    self.assertEqual(str(Uncertain(3.12345, .3)), "3.1(3)")
def test_str5(self):
    self.assertEqual(str(Uncertain(3.589, .3)), "3.6(3)")
def test_str6(self):
    self.assertEqual(str(Uncertain(3.999, .3)), "4.0(3)")
def test_str7(self):
    self.assertEqual(str(Uncertain(.001, .3)), "0.0(3)")
def test_str8(self):
    self.assertEqual(str(Uncertain(1, 1e-6)), "1.000000(1)")
def test_str9(self):
    self.assertEqual(str(Uncertain(1.000001, 1e-6)), "1.000001(1)")
def test_str12(self):
    self.assertEqual(str(Uncertain(1e6, 1.49)), "1.000000(1)e6")
def test_str13(self):
    self.assertEqual(str(Uncertain(1e6, 1.51)), "1.000000(2)e6")
def test_str15(self):
    self.assertEqual(str(Uncertain(12345678, 209)), "1.23457(2)e7")
def test_str16(self):
    self.assertEqual(str(Uncertain(.0, .0)), "0")
def test_str18(self):
    self.assertEqual(str(Uncertain(0.002534, 0.00002)), "0.00253(2)")
def test_str19(self):
    self.assertEqual(str(Uncertain(1e4, 2)), "10000(2)")
def test_str20(self):
    self.assertEqual(str(Uncertain(1e5, 20)), "1.0000(2)e5")
def test_str21(self):
    self.assertEqual(str(Uncertain(1e-5, 2e-9)), "1.0000(2)e-5")
def test_str22(self):
    self.assertEqual(str(Uncertain(1e-4, 2e-8)), "0.00010000(2)")
def test_str23(self):
    self.assertEqual(str(Uncertain(0, 0.0000001)), "0.0000000(1)")
def test_str34(self):
    self.assertEqual(str(Uncertain(-1e-9, 0.04)), "-0.00(4)")
def test_str25(self):
    self.assertEqual(str(Uncertain(-1e-9, 4)), "-0(4)")
def test_str27(self):
    self.assertEqual(str(Uncertain(0, 20000)), "0(20000)")
def test_str29(self):
    self.assertEqual(str(Uncertain(0.001, 200000)), "0(2)e5")
def test_str30(self):
    self.assertEqual(str(Uncertain(1, 200000)), "0(2)e5")
def test_str31(self):
    self.assertEqual(str(Uncertain(0.001, 0.02)), "0.00(2)")
def test_str32(self):
    self.assertEqual(str(Uncertain(0.000001, 0.00002)), "0(2)e-5")
def test_str33(self):
    self.assertEqual(str(Uncertain(np.inf, 0)), "inf")
def test_str36(self):
    self.assertEqual(str(Uncertain(np.nan, 0)), "nan")
def test_str37(self):
    self.assertEqual(str(Uncertain(np.nan, np.nan)), "nan")
def test_str39(self):
    self.assertEqual(str(Uncertain(0, -np.inf)), "0(inf)")
def test_str41(self):
    self.assertEqual(str(Uncertain(.1, 234)), "0(200)")
def test_str42(self):
    self.assertEqual(str(Uncertain(432123, 72)), "4.3212(7)e5")
def test_str43(self):
    self.assertEqual(str(Uncertain(4321, 72)), "4320(70)")
Eric-Olivier LE BIGOT 13 years, 10 months ago  # | flag

Very interesting comments, Robert!

Your Sympy approach is quite simple. One problem is that it forces the user to keep track of variables. That's something that the uncertainties package does completely transparently (the equivalent to stdev(t, "x", "y") would directly be t.std_dev(), for instance).

Your printing routine produces really nice results! It's making me think that including something similar in the uncertainties package would be a good idea. :)

Eric-Olivier LE BIGOT 13 years, 10 months ago  # | flag

PS: The printing routine above yields surprising results for 0.1±1e-50:

0.10000000000000000555111512312578270211815834045410(1)

or 0.1±1e-200 (the printing fails).

These are admittedly difficult cases, though… In any way, your printing routine contains many interesting ideas. :)