Mailman 3 value+error type (long) - SciPy-Dev

March 16, 2004


      I recently became fed up with performing error calculations on experimental data and decided to make a Python data type which I could use with Numeric. Now I can enter data with error along with the errors as in the following example and all the error calculations are automatically carried through for me. I thought there might be some interest amongst others for inclusion of this in the scipy distribution, although you might prefer it to be in C. I provide it as food for thought, but let me know if there is interest in including it as is and I can look at making it robust by adding testcases etc. and some nicer helper functions for constructors with relative errors etc. Note, the helper functions for extracting the limits would be better implemented as properties of an array type subclassed from the Numeric array type. I gather Numarray would allow this, but I couldn't think of a way to do it in Numeric and I wanted to maintain compatibility. Similarly, I couldn't think of a 
 nicer way to apply Ufuncs across the prime value and error extrema. It would be nice to not have to call ApplyUfuncToErr() to do this.
Here's a sample of usage of the class, followed by the class itself. Note I've used Gnuplot.py in this but you'll get the idea by looking at the code:

--
R_AB.py:

from Numeric import *
import Gnuplot, Gnuplot.funcutils
from ErrorVal import *
from scipy.interpolate import splrep,splev

# manometer readings [mmHg]
upperPressure = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0)
lowerPressure = ArrayOfErr([144., 246., 378., 469., 493., 505.], 1.0)

pressureDiff = upperPressure - lowerPressure

# Temperatures & Pressures from lookup table
lowerTablePressures = array([760. , 540. , 290. , 110. , 60.  , 40.  ])     # [mmHg]
lowerTableTemps     = array([4.210, 3.867, 3.332, 2.685, 2.373, 2.194])     # [K]
upperTablePressures = array([780. , 560. , 300. , 120. , 70.  , 50.  ])     # [mmHg]
upperTableTemps     = array([4.238, 3.902, 3.358, 2.735, 2.447, 2.290])     # [K]

# Linearly interpolate table values to produce temperatures
temps =  ( lowerTableTemps +
           (pressureDiff - lowerTablePressures) *
           (upperTableTemps - lowerTableTemps) /
           (upperTablePressures - lowerTablePressures) )

# Voltages across R_AB [V]
V_RStandard = ArrayOfErr([2.016, 2.016, 2.020, 2.017, 2.021, 2.019], 0.001)
R = 100.2                           # standard resistor value [Ohm]
I = V_RStandard / R                 # current [A]

# Voltages across Germanium resistor [V]
V_RAB = array([Err(6.167, 0.002),
               Err(6.168, 0.002),
               Err(6.170, 0.002),
               Err(6.160, 0.002),
               Err(6.153, 0.02),
               Err(5.894, 0.01)], PyObject)

R_AB = V_RAB / I
logRAB = ApplyUfuncToErr(R_AB, log10)       # This means log10(R_AB)
sqrtLogROnT = (logRAB/temps)**0.5

g = Gnuplot.Gnuplot(debug=1)

g.title('Allen-Bradley model log(R)=A+B[log(R)/T]^1/2')
g.ylabel('log(R)')
g.xlabel('[log(R)/T]^1/2')
g('set bar 1')                # set bar end length
g('set grid')

# spline fit
x1 = PrimeVals(sqrtLogROnT)[0]
x2 = PrimeVals(sqrtLogROnT)[5]
splineX = arange(x1,x2+.005,.01)
exactSpline = splrep(PrimeVals(sqrtLogROnT), PrimeVals(logRAB), s=0)
exactSplineY = splev(splineX, exactSpline)

g.plot(\
    Gnuplot.Data(PrimeVals(sqrtLogROnT), PrimeVals(logRAB), inline=0, with='linesp lt 1 pt 0'),
    Gnuplot.Data(splineX, exactSplineY, inline=0, with='linesp lt 2 pt 0'),     # exact spline
    Gnuplot.Data(PrimeVals(sqrtLogROnT),                        # x-points
                 PrimeVals(logRAB),                             # y-points
                 MinVals(sqrtLogROnT), MaxVals(sqrtLogROnT),    # x-axis error bars
                 MinVals(logRAB), MaxVals(logRAB),              # y-axis error bars
                 inline=0, with='xyerrorbars -1 0'),
    )

raw_input('Gary says copy to clipboard now\n')

--
ErrorVal.py:

from Numeric import *

class Err:
    def __init__(self, value, posErr=0., negErr=None):
        self.value = float(value)
        self.posErr = float(posErr)
        if negErr:
            self.negErr = float(negErr)
        else:
            self.negErr = float(posErr)

    def __repr__(self):
        return "Err(%s,%s,%s)" % (self.value, self.posErr, self.negErr)

    def __str__(self):
        return "%s +%s/-%s" % (self.value, self.posErr, self.negErr)

    def __abs__(self):
        value = abs(self.value)
        if value == self.value:
            posErr = self.posErr
            negErr = self.negErr
        else:
            negErr = self.posErr
            posErr = self.negErr
        return Err(value, posErr, negErr)

    def __neg__(self):
        return Err(-self.value, self.negErr, self.posErr)

    def __pos__(self):
        return Err(self.value, self.posErr, self.negErr)

    def __add__(self, other):
        value = self.value + other.value
        posErr = self.posErr + other.posErr
        negErr = self.negErr + other.negErr
        return Err(value, posErr, negErr)

    def __radd__(self, other):
        value = self.value + other.value
        posErr = self.posErr + other.posErr
        negErr = self.negErr + other.negErr
        return Err(value, posErr, negErr)

    def __sub__(self, other):
        value = self.value - other.value
        posErr = self.posErr + other.negErr
        negErr = self.negErr + other.posErr
        return Err(value, posErr, negErr)

    def __rsub__(self, other):
        value = other.value - self.value
        posErr = other.posErr + self.negErr
        negErr = other.negErr + self.posErr
        return Err(value, posErr, negErr)

    def __mul__(self, other):
        value = self.value * other.value
        posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
        negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
        return Err(value, posErr, negErr)

    def __rmul__(self, other):
        value = self.value * other.value
        posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
        negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
        return Err(value, posErr, negErr)

    def __div__(self, other):
        value = self.value / other.value
        posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
        negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
        return Err(value, posErr, negErr)

    def __rdiv__(self, other):
        value = other.value / self.value
        posErr = value * ((self.posErr/self.value) + (other.posErr/other.value))
        negErr = value * ((self.negErr/self.value) + (other.negErr/other.value))
        return Err(value, posErr, negErr)

    def __pow__(self, y):
        value = self.value**y.value
        posErr = (self.value + self.posErr)**(y.value + y.posErr) - value
        negErr = value - (self.value - self.negErr)**(y.value - y.negErr)
        return Err(value, posErr, negErr)

    def __rpow__(self, y):
        value = y.value**self.value
        posErr = (y.value + y.posErr)**(self.value + self.posErr) - value
        negErr = value - (y.value - y.negErr)**(self.value - self.negErr)
        return Err(value, posErr, negErr)

    def __coerce__(self, other):
        if isinstance(other, Err):
            return self, other
        try:
            return self, Err(float(other))
        except ValueError:
            pass

    def ApplyFunc(self, func):
        value = func(self.value)
        posErr = func(self.value + self.posErr) - value
        negErr = value - func(self.value - self.negErr)
        return Err(value, posErr, negErr)

def ArrayOfErr(errArray, posErr=0., negErr=None):
    '''
    Builds a Numeric array of Err() objects given an array and
    a +ve or both a +ve and -ve errorval
    Example:
    egErr = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0)
    eg2Err = ArrayOfErr([909., 802., 677., 585., 560., 548.], 1.0, 1.5)
    '''
    valArray = asarray(errArray)
    arrayLength = valArray.shape[0]
    errs = repeat(array([Err(0.0, posErr, negErr)], PyObject), arrayLength)
    return valArray + errs

def PrimeVals(errArray):
    '''
    Extract a Numeric array with just the prime value of the Err() objects
    '''
    li = []
    for i in errArray:
        li.append(i.value)
    return asarray(li)

def PosErrs(errArray):
    '''
    Return a Numeric array with the posErr values of the Err() objects
    '''
    li = []
    for i in errArray:
        li.append(i.posErr)
    return asarray(li)

def MaxVals(errArray):
    '''
    Return a Numeric array with the prime+posErr values of the Err() objects
    '''
    li = []
    for i in errArray:
        li.append(i.value + i.posErr)
    return asarray(li)

def NegErrs(errArray):
    '''
    Return a Numeric array with the negErr values of the Err() objects
    '''
    li = []
    for i in errArray:
        li.append(i.negErr)
    return asarray(li)

def MinVals(errArray):
    '''
    Return a Numeric array with the prime-negArr values of the Err() objects
    '''
    li = []
    for i in errArray:
        li.append(i.value - i.negErr)
    return asarray(li)

def ErrToFloatArrays(errArray):
    '''
    Return a Numeric array with 3 rows:
    row 1: the prime value of the Err() objects
    row 2: the prime+posErr values of the Err() objects
    row 3: the prime-negErr values of the Err() objects
    '''
    r1 = PrimeVals(errArray)
    r2 = MaxVals(errArray)
    r3 = MinVals(errArray)
    return array([r1, r2, r3])

def FloatToErrArrays(floatArray):
    '''
    Convert a Numeric array of the form produced by ErrToFloatArrays, which contains 3 rows,
    to an errArray
    '''
    r1 = floatArray[0]
    r2 = floatArray[1]-r1
    r3 = r1-floatArray[2]
    return asarray(map(Err,r1,r2,r3), PyObject)

def ApplyUfuncToErr(errArray, func):
    '''
    Example:
    a = ArrayOfErr([99., 82., 67., 58., 50., 54.], 1.0)
    b = ApplyUfuncToErr(a, log)         # This means log(a)
    '''
    a = ErrToFloatArrays(errArray)
    b = func(a)
    c = FloatToErrArrays(b)
    return c

--
-- 
___________________________________________________________
Sign-up for Ads Free at Mail.com
http://promo.mail.com/adsfreejump.htm

value+error type (long)

Gary Ruben

tags

participants (1)