{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Error propagation: Three approaches\n", "- H&H 'Functional' Approach\n", "- Monte Carlo simulation\n", "- H&H 'Calculus' Approach\n", "\n", "We will use a simple example of a formula for a physical quantitity $g$ which \n", "depends on $x$, $y$, and $z$:\n", "\n", "$$\n", "f(x,y,z) = \\frac{xy}{z^2}.\n", "$$\n", "\n", "Let's assume that we have measured the following values:\n", " + $x = 1.00 \\pm 0.05$\n", " + $y = 2.00 \\pm 0.07$\n", " + $z = 3.00 \\pm 0.12$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def f(x,y,z):\n", " '''Model function for error propagation examples'''\n", " return x*y/z**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assign values to variables $x$, $y$, and $z$, with associated uncertainties " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "x = 1.\n", "alpha_x = 0.05\n", "\n", "y = 2.\n", "alpha_y = 0.07\n", "\n", "z = 3.\n", "alpha_z = 0.12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### \"Best\" value of $f$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.2222222222222222\n" ] } ], "source": [ "best = f(x,y,z)\n", "print(best)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## \"Functional approach\"\n", "\n", "\\begin{eqnarray*}\n", "(\\alpha_f)_x &=& f(x+\\alpha_x,y,z) - f(x,y,z) \\\\\n", "(\\alpha_f)_y &=& f(x, y+\\alpha_y,z) - f(x,y,z)\\\\\n", "(\\alpha_f)_z &=& f(x,y, z+\\alpha_z) - f(x,y,z)\\\\\n", "\\mbox{}\\\\\n", "\\alpha_f &=& \\sqrt{(\\alpha_f)_x^2 + (\\alpha_f)_y^2 +(\\alpha_f)_z^2 }\n", "\\end{eqnarray*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an array`u` to hold the values $(\\alpha_f)_x$, $(\\alpha_f)_x$, and $(\\alpha_f)_z$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "u = np.zeros(3)\n", "u[0] = f(x+alpha_x, y, z)\n", "u[1] = f(x, y+alpha_y, z)\n", "u[2] = f(x, y, z + alpha_z)\n", "u = u-best" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate total uncertainty\n", "\n", "$$\n", "\\alpha_f = \\sqrt{\\left(\\alpha_f\\right)^2_x + \\left(\\alpha_f\\right)^2_y + \\left(\\alpha_f\\right)^2_z}\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.021564448330840195\n" ] } ], "source": [ "unc = np.sqrt(u@u)\n", "print(unc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Final Result: $0.22\\pm 0.02$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Fractional uncertainties" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fractional uncertainties [ 0.05 0.035 -0.07544379]\n" ] } ], "source": [ "print(\"Fractional uncertainties\", u/best)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "$$ \\frac{\\alpha_x}{f} = 5\\%, \\quad\\quad \\frac{\\alpha_y}{f}= 3.5\\%, \\quad\\mbox{and}\\quad \n", "\\frac{\\alpha_z}{f} = 7.5\\% $$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo technique\n", "\n", "When we write something like $x = 1.00 \\pm 0.05$ it means that we believe\n", "that $x$ is a random variable selected from a normal distribution with \n", "a mean of $1.00$ and a standard deviation of $0.05$. In this technique we \n", "simulate very many experiments with values drawn from the appropriate \n", "distributions." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "nEx = 10000 #Number of simulated experiments\n", "x_sim = stats.norm.rvs(x,alpha_x,nEx) # array for simulated values of x\n", "y_sim = stats.norm.rvs(y,alpha_y,nEx) # array for simulated values of y\n", "z_sim = stats.norm.rvs(z,alpha_z,nEx) # array for simulated values of z " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate values of $f$ using simulated data" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.22355457533480017 0.022514402883885898\n" ] } ], "source": [ "f_sim = f(x_sim, y_sim, z_sim)\n", "print(np.mean(f_sim), np.std(f_sim))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Result: $0.22 \\pm 0.02$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## \"Calculus approach\"\n", "\\begin{eqnarray*}\n", "(\\alpha_f)_x &=& \\left. \\frac{\\partial f}{\\partial x}\\right|_{x = 1.0\\, y = 2.0\\, z=3.0}\\, \\alpha_x \\\\\n", "(\\alpha_f)_y &=& \\left. \\frac{\\partial f}{\\partial y}\\right|_{x = 1.0\\, y = 2.0\\, z=3.0}\\, \\alpha_y \\\\\n", "(\\alpha_f)_z &=& \\left. \\frac{\\partial f}{\\partial z}\\right|_{x = 1.0\\, y = 2.0\\, z=3.0}\\, \\alpha_z\\\\\n", "\\mbox{}\\\\\n", "\\alpha_f &=& \\sqrt{(\\alpha_f)_x^2 +(\\alpha_f)_y^2 + (\\alpha_f)_z^2}\n", "\\end{eqnarray*}\n", "\n", "In this example the partial derivatives are easy:\n", "\n", "$$\n", "\\frac{\\partial f}{\\partial x} = \\frac{y}{z^2}\\quad\\quad \n", "\\frac{\\partial f}{\\partial y} = \\frac{x}{z^2}\\quad\\quad\n", "\\frac{\\partial f}{\\partial y} = -\\frac{2xy}{z^3}.\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.023726840560069584\n" ] } ], "source": [ "alpha_f_x = alpha_x*y/z**2\n", "alpha_f_y = alpha_y*x/z**2\n", "alpha_f_z = -2*alpha_z*x*y/z**3\n", "\n", "alpha_f = np.sqrt(alpha_f_x**2 + alpha_f_x**2 + alpha_f_z**2 )\n", "\n", "print(alpha_f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Result: $0.22 \\pm 0.2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Version Information\n", "\n", "`version_information` is from J.R. Johansson (jrjohansson at gmail.com)
\n", "See Introduction to scientific computing with Python:
\n", "http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb
\n", "for more information and instructions for package installation.
\n", "\n", "If `version_information` has been installed system wide (as it has been on Bucknell linux computers with shared file systems), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "%load_ext version_information" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.7.7 64bit [GCC 7.3.0]" }, { "module": "IPython", "version": "7.16.1" }, { "module": "OS", "version": "Linux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core" }, { "module": "numpy", "version": "1.18.5" }, { "module": "scipy", "version": "1.5.2" } ] }, "text/html": [ "
SoftwareVersion
Python3.7.7 64bit [GCC 7.3.0]
IPython7.16.1
OSLinux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core
numpy1.18.5
scipy1.5.2
Wed Aug 05 20:18:56 2020 EDT
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.7.7 64bit [GCC 7.3.0] \\\\ \\hline\n", "IPython & 7.16.1 \\\\ \\hline\n", "OS & Linux 3.10.0 1062.9.1.el7.x86\\_64 x86\\_64 with centos 7.7.1908 Core \\\\ \\hline\n", "numpy & 1.18.5 \\\\ \\hline\n", "scipy & 1.5.2 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Wed Aug 05 20:18:56 2020 EDT} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.7.7 64bit [GCC 7.3.0]\n", "IPython 7.16.1\n", "OS Linux 3.10.0 1062.9.1.el7.x86_64 x86_64 with centos 7.7.1908 Core\n", "numpy 1.18.5\n", "scipy 1.5.2\n", "Wed Aug 05 20:18:56 2020 EDT" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%version_information numpy, scipy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "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.7.8" } }, "nbformat": 4, "nbformat_minor": 1 }