Propagation of uncertainties

The module nbi_stat provides a generic mechanism to propagate uncertainties via the function propagate_uncertainty(). It takes a expression (defined as a Python call-able), variable values and uncertainties and calculates the uncertainty on the expression. It does so by way of numerical differentiation of the passed expression.

Simple propagation of uncertainties

Suppose we measured

\[\begin{split}x &= 1.3 \pm 0.1\\ y &= 0.22 \pm 0.05\end{split}\]

and we want to calculate the uncertainty of

\[f(x,y) = x^y\]

Then we can do

>>> def f(x,y):
...     return x**y
... 
>>> xy = numpy.array((1.3,0.22))
>>> exy = numpy.array((0.1,0.05))
>>> z   = f(*xy)
>>> ez  = numpy.sqrt(nbi_stat.propagate_uncertainty(lambda p:f(*p),xy,exy))
>>> nbi_stat.print_result(z, [ez], 1)
    1.06 +/-     0.02

By default, the uncertainties on the independent variables (\(x\), and \(y\) above) are used as the step size for the numerical differentiation

\[\left.\frac{\mathrm{d}f}{\mathrm{d}x}\right|_{x=x_0} \approx \frac{f(x_0+\Delta x) - f(x_0-\Delta x)}{2\Delta x}\]

However, if the uncertainties are large, or the function is not in some otherway, well approximated by this, one can pass step-sizes for each component.

>>> def f(x,y):
...     return x**y
... 
>>> xy = numpy.array((1.3,0.8))
>>> exy = numpy.array((0.1,0.09))
>>> z   = f(*xy)
>>> ez  = numpy.sqrt(nbi_stat.propagate_uncertainty(lambda p:f(*p),xy,exy,0.2*exy))
>>> nbi_stat.print_result(z, [ez], 1)
    1.23 +/-     0.08

Propagation of uncertainties on function parameters

Suppose we did a measurement which we fitted the expression

\[f(x;a,b) = ax^b\]

to, and found the parameter estimates

\[\begin{split}a &= 0.9 \pm 0.2\\ b &= 0.32 \pm 0.05\end{split}\]

Now, we’d like to draw this function together with an uncertainty band. For this, we can also use nbi_stat.propagate_uncertainty() albeit with a bit more work

>>> import nbi_stat
>>>
>>> def f(x,a,b):
>>>     return a*x**b
>>>
>>> p = np.array((0.9,0.32))
>>> ep = np.array((0.2,0.05))
>>>
>>> x = np.linspace(0,10)
>>> y = f(x,*p)
>>> ey = np.sqrt([nbi_stat.propagate_uncertainty(lambda p:f(xi,*p),p,ep)
>>>               for xi in x])
>>>
>>> plt.fill_between(x,y-ey,y+ey,
>>>                  label=r"$\delta f$",alpha=0.5,color='y')
>>> plt.plot(x,y,label=r"$f$")
>>> plt.legend()

(Source code, png, hires.png, pdf)

_images/propagation-1.png

Note, the function nbi_stat.fit_plot() uses this to show an uncertainty band around a fitted function.