DMelt:Units/3 Error Propagation

Table of contents

Error propagation

Propagation of uncertainty (or propagation of error) is the effect of variables' uncertainties (or errors, more specifically random errors) on the uncertainty of a function based on them. Error (or uncertainty) propagation for arbitrary (analytic and non-analytic) functions are supported using several methods:

• Using a Taylor expansion
• Using a Monte Carlo simulation

Error propagation using P1D

A typical measurement contains an uncertainty. A convenient way to keep data with errors is to use the  jhplot.P1D data container.

from jhplot import *
p=P1D("data with errors")
p.add(1, 2, 0.5)  # error on Y=2 is 0.5
p.add(2, 5, 0.9)  # error on Y=5 is 0.9

When you are using the method

"oper()" of this class, the errors are automatically propagated. This applies for subtraction, division, multiplication and addition.

Error propagation

DataMelt links the Python package uncertainties, therefore, you can program with Python as usual. Run this example in the DataMelt IDE and you can see this:

from uncertainties import ufloat
from uncertainties.umath import *  # sin(), etc.
x = ufloat(1, 0.1)  # x = 1+/-0.1

This approach calculates errors using a Taylor expansion. It also can differentiate a function.

Error propagation for complex functions

In this section we describe error propagation for an arbitrary transformation using the Java classes.

Error propagation for arbitrary transformation can be done by the class  jhplot.math.ValueErr. The approach is based on a Taylor expansion (linearization) and still have some limitations. Here is a simple example:

from jhplot.math import *
x=ValueErr(10,2)
ylog=x.log10(x)  # log10 and associated error
print ylog.toString()

Uncertainty propagation using a Monte Carlo method

Complex functions are difficult to deal with using the above approach. DataMelt includes Java classes for a Monte Carlo approach which:

• avoids any errors associated with the linearization of the model
• produces a distribution for the uncertain output as well as the mean and standard deviation.

Propagation of errors formula requires the computation of derivatives, which can be quite complicated for larger models. The Monte Carlo approach can handle correlated parameters as well as independent parameters.

Error propagation for arbitrary transformation using a Monte Carlo can be done by the class  jhpro.upropog.UPropogate. Here is a simple example:

from java.lang import Math 
from jhplot import *
from jhpro.upropog import UPropogate

class MyFunc(FNon):   
     def value(self, x):                
              y=x[0]*x[0]*Math.sqrt(x[0])*self.p[0]+self.p[1]
              if (x[0]<4): y=0
              return y 

pl = MyFunc("test",1,2)    
print "f(2)=  ",pl.value([2])
print pl.numberOfParameters() 
print pl.dimension()
print pl.parameterNames().tolist() 

# define parameters and their errors
p0=100
p1=20

p0_err=10
p1_err=2

# do calculations at x=10 (with zero error)
xval=10
xval_err=0

ff=UPropogate(pl,[p0,p1],[p0_err,p1_err],[xval],[xval_err],10000)
err=ff.getStd()
pl.setParameter("par0",p0) 
pl.setParameter("par1",p1)
val=pl.value([xval])
print "val=",val," +- ",err


hh=ff.getHisto()
c1=HPlot()
c1.visible()
c1.setAutoRange()
c1.draw(hh)

and here is a more custom approach without this class:

from java.lang import Math 
from jhplot import *
from java.util import Date;

class MyFunc(FNon):   
     def value(self, x):                
              y=x[0]*x[0]*Math.sqrt(x[0])*self.p[0]+self.p[1]
              if (x[0]<4): y=0
              return y 

pl = MyFunc("test",1,2)    
print "f(2)=  ",pl.value([2])
print pl.numberOfParameters() 
print pl.dimension()
print pl.parameterNames().tolist() 

p0=100
p1=20

p0_err=10
p1_err=2

from cern.jet.random import Normal
from cern.jet.random.engine import MersenneTwister
from  hep.aida.ref.histogram import Histogram1D

engine=MersenneTwister(Date())
rd0=Normal(p0, p0_err,engine)
rd1=Normal(p1, p1_err,engine)

xval=10;


pl.setParameter("par0",p0+4*p0_err) 
pl.setParameter("par1",p1-4*p1_err)
x1=pl.value([xval])

pl.setParameter("par0",p0+4*p0_err) 
pl.setParameter("par1",p1+4*p1_err)
x2=pl.value([xval])

pl.setParameter("par0",p0-4*p0_err) 
pl.setParameter("par1",p1-4*p1_err)
x3=pl.value([xval])

pl.setParameter("par0",p0-4*p0_err) 
pl.setParameter("par1",p1+4*p1_err)
x4=pl.value([xval])

z=P0D("data")
z.add(x1)
z.add(x2)
z.add(x3)
z.add(x4)
hh=H1D("error",100,z.getMin(),z.getMax())

for i in range(10000):
    p0n=rd0.nextDouble()
    p1n=rd1.nextDouble()
    pl.setParameter("par0",p0n) 
    pl.setParameter("par1",p1n) 
    #print pl.parameters() 
    ff=pl.value([xval])
    hh.fill(ff)


all=hh.getStat()
err=all['standardDeviation']
pl.setParameter("par0",p0) 
pl.setParameter("par1",p1)
val=pl.value([xval])
print "val=",val," +- ",err


c1=HPlot()
c1.visible()
c1.setAutoRange()
c1.draw(hh)

Error propagation using VisAd

The Java package  VisAd contains powerful classes  visad.Real and  visad.Data to perform calculations with units and errors.

In VisAD the values of a Scalar may be either Real which is used for numeric quantities. The form of constructor for Real:

Real(double value, double error)

allows us to provide an error estimate for the value. This estimate can then be propagated through mathematical operations to provide an error estimate. Here is Java example:

import visad.Real;
Real a,b,c;
a = new Real(10.,1.);
b = new Real(255.,10.);
c = (Real) a.add(b);
System.out.println("sum = "+c);

When you run this example, you still get: sum = 265.0

The VisAD default for most math operations, however, is to not propagate errors, so to make use of this, we must explicitly indicate how we want error estimate to propagate. This is done by using an alternate signature of the add method (and all other math operators):

import visad.Real;
Real a,b,c;
a = new Real(10.,1.);
b = new Real(255.,10.);
c = (Real) a.add(b, Data.NEAREST_NEIGHBOR, Data.INDEPENDENT);
System.out.println("sum = "+c);
System.out.println("error of sum is="+c.getError().getErrorValue());

When you run this example, you get: "sum = 265.0" and "error of sum is=10.04987562112089"

The constants supplied to the add method are the type of interpolation and the type of error propagation. In this simple case, the type of interpolation is not really relevant, but as you will see later, VisAD Data types may contain finite approximations to continuous functions and when these are combined mathematically, may need to be resampled in order to match domains.

This example shows how to define values with errors and how to perform mathematical transformation while propagating errors (independently) using Jython/Python:

from visad import SI,Data,Real,RealType

def printReal(x):
   print x.toString(), "+-", (x.getError()).getErrorValue()

a=Real(10,2) # 10 +- 2 (error)
b=Real(30,5) # 30 +- 5 (error)

c=a.add(b,Data.NEAREST_NEIGHBOR, Data.INDEPENDENT)
print "a+b=", printReal(c)

c=a.divide(b,Data.NEAREST_NEIGHBOR, Data.INDEPENDENT)
print "a+b=", printReal(c)

a=Real(0.9,0.1) # 10 +- 2 (error)
c=a.cos(Data.NEAREST_NEIGHBOR, Data.INDEPENDENT)
print "acos(a)=", printReal(c)

a=Real(0.9,0.1) # 10 +- 2 (error)
c=a.pow(c,Data.NEAREST_NEIGHBOR, Data.INDEPENDENT)
print "a^4=", printReal(c)

print "use units:"
k=RealType("kelvin",SI.kelvin)
t2=Real(k,273.);
t1= Real(k,255.);
xsum =t1.add(t2);
print "sum = ",xsum," ",xsum.getUnit()