# 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()


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