{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Propagation of uncertainty (error)\n", "\n", "> Marcos Duarte \n", "> [Laboratory of Biomechanics and Motor Control](https://bmclab.pesquisa.ufabc.edu.br/) \n", "> Federal University of ABC, Brazil" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" }, "toc": true }, "source": [ "

Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Evaluation of measurement data\n", "\n", "> The result of a measurement is only an approximation or estimate of the value of the measurand and thus is complete only when accompanied by a statement of the uncertainty of that estimate. \n", "> \n", "> Uncertainty (of measurement) is a parameter, associated with the result of a measurement, that characterizes the dispersion of the values that could reasonably be attributed to the measurand. \n", "> \n", "> [Evaluation of measurement data - Guide to the expression of uncertainty in measurement (2008)](https://www.bipm.org/utils/common/documents/jcgm/JCGM_100_2008_E.pdf)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "For instance, to calculate the density of a material someone made the following measurements of mass and volume (m=d/V): \n", "\n", "$$ m = 4.0 \\pm 0.5 \\; g $$ \n", "$$ V = 2.0 \\pm 0.2 \\; cm^3 $$\n", "\n", "Where 0.5 g and 0.2 cm$^3$ each represent a value of one standard deviation and express the errors in the measured weight and volume. \n", "\n", "Based on these numbers, the material density is 2 g/cm$^3$; but how much is the error (uncertainty) of the density? " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "One way to estimate the error in the density is by propagation of uncertainty:\n", "\n", "> Propagation of uncertainty (or propagation of error) is the effect of variables' uncertainties (or errors) on the uncertainty of a function based on them. When the variables are the values of experimental measurements they have uncertainties due to measurement limitations (e.g., instrument precision and noise) which propagate to the combination of variables in the function ([Wikipedia](http://en.wikipedia.org/wiki/Propagation_of_uncertainty))." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Propagation of uncertainty by linear approximation\n", "\n", "The classic way for propagation of uncertainty is through linear approximation of error propagation theory. We start with the mathematical function that relates the physical quantities and estimate the contribution of their variations by means of partial derivatives (the rate of variation) of each of these primary quantities. The linearity of the mathematical function in the vicinity of the obtained value and the statistical correlations between these quantities are also taken into consideration in the propagation of uncertainty. \n", "\n", "If the measurand $f$ is defined in terms of the variables $x, y, z, ...$ by a general expression $f(x, y, z, ...)$, a first order approximation of the propagation of uncertainty ignoring correlations between these variables is given by the well known formula:\n", "\n", "$$ \\sigma _{f}\\;=\\;\\sqrt{\\left ( \\frac{\\partial f}{\\partial x} \\right)^2\\sigma _{x}^{2} + \\left ( \\frac{\\partial f}{\\partial y} \\right)^2 \\sigma _{y}^{2} + \\left ( \\frac{\\partial f}{\\partial z} \\right)^2 \\sigma _{z}^{2}\\: + \\ldots} $$\n", "\n", "Where $ \\partial f/\\partial i $ is the first partial derivative of $ f $ with respect to variable $ i $ and $ \\sigma_i $ is the uncertainty of the measurement of the variable $ i $ (one standard deviation). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Using the former formula for linear approximation of error propagation, the estimation of the uncertainty of the density is given by:\n", "\n", "$$ d=\\frac{m}{V} \\\\\n", "\\sigma _{d}\\;=\\;\\sqrt{ \\left ( \\frac{\\partial d}{\\partial m} \\right)^2 \\sigma _{m}^{2}+\\left ( \\frac{\\partial d}{\\partial V} \\right)^2 \\sigma _{V}^{2} } \\\\\n", "\\sigma _{d}\\;=\\;\\sqrt{ \\left ( \\frac{1}{V} \\right)^2 \\sigma _{m}^{2}+\\left ( -\\frac{m}{V^2} \\right)^2 \\sigma _{V}^{2} } $$\n", "\n", "And considering the values of mass and volume given earlier:\n", "\n", "$$ \\sigma_{d}\\;=\\;\\sqrt{\\left(\\frac{1}{2}\\right)^2\\times 0.5\\:^2+\\left(-\\frac{4}{2^2}\\right)^2\\times 0.2\\:^2}\\;\\;\\approx\\; 0.32 $$\n", "\n", "Then, the estimated value of the density can now be expressed in a complete form: d = 2.00±0.32 g/cm$^3$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Using the `uncertainties` package to calculate uncertainty\n", "\n", "It is essential that any result of a measurement is accompanied by the estimate of its uncertainty as above, this way we can better understand the value of the measure and how reliable it is. However, it is unlikely that anyone will repeat these calculations in the day-to-day work (although there are tables with partial derivatives for the most common expressions). Therefore, it is useful a software that can perform these calculations (at least for checking the results manually obtained). Fortunately, there are software for the propagation of uncertainty, [see this list on Wikipedia](https://en.wikipedia.org/wiki/List_of_uncertainty_propagation_software). \n", "\n", "One such software, free, open source, and cross-platform, is the library for [Python](https://python.org) called [uncertainties](https://pypi.org/project/uncertainties/). With [uncertainties](https://pypi.org/project/uncertainties/) installed, let's calculate again the uncertainty of the density:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:03:29.767571Z", "start_time": "2017-12-30T07:03:29.681534Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Import the necessary libraries\n", "from IPython.display import display, Math # IPython formatting\n", "from uncertainties import ufloat # to define variables with uncertainty" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Let's input the values:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:03:31.970156Z", "start_time": "2017-12-30T07:03:31.958144Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.00+/-0.32\n" ] } ], "source": [ "m = ufloat(nominal_value=4.0, std_dev=.5) # value with uncertainty\n", "V = ufloat(2.0, 0.2) # we can pass values directly\n", "d = m/V # calculate density\n", "\n", "print(d)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The result we obtained before, but without having to deduct the partial derivatives because [uncertainties](http://pythonhosted.org/uncertainties/) does this with an [algorithm for automatic differentiation](http://en.wikipedia.org/wiki/Automatic_differentiation)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Displaying the output in different formats" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:03:43.904812Z", "start_time": "2017-12-30T07:03:43.881866Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle d = 2.00 \\pm 0.32 g/cm^3$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle d = \\left(2.00 \\pm 0.32\\right) \\times 10^{0} g/cm^3$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle d = 2.00(32) g/cm^3$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(Math(f'd = {d:.3L} g/cm^3'))\n", "display(Math(f'd = {d:.2eL} g/cm^3'))\n", "display(Math(f'd = {d:.2uS} g/cm^3'))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The operation above could have been performed in a single line (after you imported the necessary libraries at the first time using them):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:04:37.275122Z", "start_time": "2017-12-30T07:04:37.269857Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.00+/-0.32\n" ] } ], "source": [ "print(ufloat(4.0, 0.5) / ufloat(2.0, 0.2))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "And how much each variable mass and volume contributed to the final uncertainty (i.e., the partial derivative of density in relation to the variable mass or volume times its standard deviation (each term on the right side of the equation above without squaring)) is given by:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:06:25.918229Z", "start_time": "2017-12-30T07:06:25.913124Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "dict_items([(2.0+/-0.2, 0.2), (4.0+/-0.5, 0.25)])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.error_components().items()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "By definition, according to the formula above, the square root of the quadratic sum of these values is the total uncertainty for the density. \n", "\n", "The value of each partial derivative of the density function with respect to the mass and volume variables is given by:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:08:09.720180Z", "start_time": "2017-12-30T07:08:09.712447Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "defaultdict(None, {2.0+/-0.2: -1.0, 4.0+/-0.5: 0.5})" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.derivatives" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Automatic deduction of the symbolic formula for error propagation \n", "\n", "Stepping back, we can even deduce the formula for the error propagation using the partial derivatives in case we want the actual formula using [Sympy](http://sympy.org/en/index.html):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:08:29.884899Z", "start_time": "2017-12-30T07:08:29.352608Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "from sympy import symbols, Symbol, sqrt, Add, Mul, Pow, init_printing, latex\n", "init_printing()\n", "\n", "def stdform(formula, formula_symb, *errors):\n", " \"\"\"\n", " Calculate the symbolic formula for error propagation using partial derivatives.\n", " \"\"\"\n", " from functools import reduce\n", "\n", " formula2 = sqrt(reduce(Add, (Mul(Pow(formula.diff(var), 2, evaluate=False),\n", " Symbol(rf'\\sigma_{var.name}', commutative=False)**2, evaluate=False)\n", " for var in formula.atoms(Symbol) if var.name in errors), 0.))\n", "\n", " formula3 = sqrt(reduce(Add, (Mul(Pow(formula.diff(var)/formula, 2, evaluate=False),\n", " Symbol(rf'\\sigma_{var.name}', commutative=False)**2, evaluate=False)\n", " for var in formula.atoms(Symbol) if var.name in errors), 0.))\n", "\n", " print('Uncertainty:')\n", " display(Math(rf'\\sigma_{formula_symb} = {latex(formula2)}'))\n", " print('Relative Uncertainty:')\n", " display(Math(rf'\\frac{{\\sigma_{formula_symb}}}{formula_symb} = {latex(formula3)}'))\n", "\n", " return formula2, formula3" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:08:32.163956Z", "start_time": "2017-12-30T07:08:32.138222Z" }, "slideshow": { "slide_type": "subslide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Example for uncertainty propagation for the density, d=m/v:\n", "Uncertainty:\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\sigma_d = \\sqrt{\\left(\\frac{1}{V}\\right)^{2} \\sigma_{m}^{2} + \\left(- \\frac{m}{V^{2}}\\right)^{2} \\sigma_{V}^{2}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Relative Uncertainty:\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\frac{\\sigma_d}d = \\sqrt{\\left(- \\frac{1}{V}\\right)^{2} \\sigma_{V}^{2} + \\left(\\frac{1}{m}\\right)^{2} \\sigma_{m}^{2}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print('Example for uncertainty propagation for the density, d=m/v:')\n", "m, V = symbols('m, V')\n", "d = m/V\n", "formula = stdform(d, 'd', 'm', 'V')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Same as before, but Sympy deduced for us." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Another example\n", "\n", "Let's estimate the uncertainty in the measurement of the volume of a cylinder given its diameter ($d=2.0\\pm0.1$) cm and height $h=10.0\\pm0.5$ cm. But now, let's go straight to the solution:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:09:02.652765Z", "start_time": "2017-12-30T07:09:02.611694Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Uncertainty:\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\sigma_V = \\sqrt{\\left(\\frac{\\pi d^{2}}{4}\\right)^{2} \\sigma_{h}^{2} + \\left(\\frac{\\pi d h}{2}\\right)^{2} \\sigma_{d}^{2}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Relative Uncertainty:\n" ] }, { "data": { "text/latex": [ "$\\displaystyle \\frac{\\sigma_V}V = \\sqrt{\\left(\\frac{2}{d}\\right)^{2} \\sigma_{d}^{2} + \\left(\\frac{1}{h}\\right)^{2} \\sigma_{h}^{2}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from sympy import pi\n", "d, h = symbols('d, h')\n", "V = h*pi*(d/2)**2\n", "\n", "formula = stdform(V, 'V', 'd', 'h')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2017-12-30T07:09:04.141161Z", "start_time": "2017-12-30T07:09:04.129675Z" }, "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle V = 31.4 \\pm 3.5 cm^3$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from math import pi\n", "h = ufloat(nominal_value=10.0, std_dev=0.5)\n", "d = ufloat(nominal_value=2.0, std_dev=0.1)\n", "V = h*pi*(d/2)**2\n", "display(Math(f'V = {V:.3L} cm^3'))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Other Python packages for uncertainty analysis\n", "\n", "The uncertainty package employs a linear approximation of error propagation, which is only valid for small errors, and it assumes that the errors follow normal distributions. There are other packages that can be used when these assumptions are violated: \n", " - [NIST Uncertainty Machine](https://uncertainty.nist.gov/): *A Web-based software application to evaluate the measurement uncertainty associated with an output quantity defined by a measurement model of the form $y=f(x_0,\t\\ldots,x_n)$.* The NIST Uncertainty Machine allows to select a probability distribution for each of the input quantities and can evaluate the measurement uncertainty using two methods: the linear approximation and the Monte Carlo method,\n", " - [mcerp](https://github.com/tisimst/mcerp): *A stochastic calculator for Monte Carlo methods that uses latin-hypercube sampling to perform non-order specific error propagation (or uncertainty analysis).* \n", " - [soerp](https://github.com/tisimst/soerp): *Python implementation of the original Fortran code SOERP by N. D. Cox to apply a second-order analysis to error propagation (or uncertainty analysis). The soerp package allows you to easily and transparently track the effects of uncertainty through mathematical calculations.* " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "And there are other packages using linear error propagation: \n", "\n", " - [scinum](https://github.com/riga/scinum): *Provides a simple Number class that wraps plain floats or NumPy arrays and adds support for multiple uncertainties, automatic (gaussian) error propagation, and scientific rounding.* \n", " - [GUM Tree Calculator](https://github.com/MSLNZ/GTC): *A Python package for processing data with measurement uncertainty.* " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## In conclusion\n", "\n", "There is no excuse to not perform propagation of uncertainty of the result of a measurement. Remember, your work is complete only when you report the uncertainty." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## References\n", "\n", "- [Evaluation of measurement data - Guide to the expression of uncertainty in measurement](https://www.bipm.org/utils/common/documents/jcgm/JCGM_100_2008_E.pdf). JCGM 100:2008. \n", "- [Avaliação de dados de medição - Guia para a expressão de incerteza de medição](http://www.inmetro.gov.br/inovacao/publicacoes/gum_final.pdf). GUM 2008 in Portuguese.\n", "- [uncertainties](https://pythonhosted.org/uncertainties/): a Python package for calculations with uncertainties, Eric O. Lebigot. " ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }