Propagation of uncertainties
============================

The module :mod:`nbi_stat` provides a generic mechanism to propagate
uncertainties via the function :meth:`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

.. math::

   x &= 1.3 \pm 0.1\\
   y &= 0.22 \pm 0.05

and we want to calculate the uncertainty of

.. math::

   f(x,y) = x^y

Then we can do

.. runblock:: pycon

   >>> import sys; sys.path.insert(0,"..") # ignore
   >>> import nbi_stat # ignore
   >>> import numpy # ignore
   >>> 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)
   

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

.. math::

   \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.

.. runblock:: pycon

   >>> import sys; sys.path.insert(0,"..") # ignore
   >>> import nbi_stat # ignore
   >>> import numpy # ignore
   >>> 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)


Propagation of uncertainties on function parameters
---------------------------------------------------

Suppose we did a measurement which we fitted the expression

.. math::

   f(x;a,b) = ax^b

to, and found the parameter estimates 

.. math::

   a &= 0.9 \pm 0.2\\
   b &= 0.32 \pm 0.05

Now, we'd like to draw this function together with an uncertainty
band.  For this, we can also use :meth:`nbi_stat.propagate_uncertainty`
albeit with a bit more work

.. plot::
   :include-source:
   
   >>> 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 :meth:`nbi_stat.fit_plot` uses this to show an
uncertainty band around a fitted function.