{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Error propagation: Three approaches\n", "- H&H 'Functional' Approach\n", "- H&H 'Calculus' Approach\n", "- Monte Carlo simulation" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import scipy as sp\n", "from scipy import stats\n", "import sympy as sym # for symbolic differentiation\n", "from sympy.interactive import init_printing # provides LaTex formatted cell output\n", "sym.init_printing()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In PHYS 211 we did an experiment to determine \"little $g$\" from measurement of the length and period of a pendulum. The following cell defines a function giving $g$ for known values of $l$ and $T$; the next cell it gives typical data and associated uncertainties." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def g(l,t):\n", " return 4*sp.pi**2*l/t**2" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "mean_l = 1.0 # Measured length of pendulum in meters\n", "sigma_l = 0.004 # Uncertainty in length - assumed to be \n", " #i.d. of normal distribution \n", "mean_t = 2.0; # Measured period\n", "sigma_t = 0.005 # Uncertainty in period - assumed to be \n", " # s.d. of normal distribution" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "### \"Functional approach\"\n", "\n", "\\begin{eqnarray*}\n", "(\\sigma_g)_L &=& g(L+\\Delta L, T) - g(L,T) \\\\\n", "(\\sigma_g)_T &=& g(L, T+\\Delta T) - g(L,T)\\\\\n", "\\sigma_g &=& \\sqrt{(\\sigma_g)_L^2 + (\\sigma_g)_T^2}\n", "\\end{eqnarray*}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAAUBAMAAAC9jfsOAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQ\nqzLsm4+cAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJUUlEQVRoBe1ZXWybZxV+/JPEjj875qeTpl3U\nnZBGNU3LqECABjObxBB01GhEYs0KQR0T0qiaSqsKGlCDxAViIymdUFcF5gsY4qpWS6gmr4svtlXV\nmGIqwihQ4ptyM6H+ha3d0pnnnPPa7/s5qbhJql3slfK+5/05z/Occ/z92AESJbzf3rMZ2CTKPgS8\n3Pizajz5xHaOjSfabkWtoy/KXqHiupdtzpla0b7ZGqLxE0UkGseqwP3Nc1UYDo+Ereuo4y2NT/Zo\nz6Gwa7zRgEI5JjraoqnpafAktr3hwc8CibEt7ZApZsdogXM9WjmVO/7i6rpVS/r7ha85LNMrMvyO\n1xIj5CQav70tazb2SGyqGTLdkrEgap9+pfv3QqNRzxeJ8zlEv8DNYqGCfBUv1QbKtiIWvoBXZG9w\n1HU2/2kTtjMA/A2nSokD2IDMJeDZzjsOh+djzQGpW6KCqVKXdieSnU7nGhTKmDIzcIt9GjyJ+kQt\nzNXwGnJvx6jCSUjL9Z3daCWC6Iv46krdsqNasp1OnS6cq94gYN3xWngo1tLV6JuyoKMnsWXJkOlW\ngCDqXvqN7k5mpRWVmQluTCNfIWKmhEQZH0dq1FbESk4jx1P4npSInc6jhcWmWfgx8GGcARawv4o3\ngS+9XnQ4PB+2LpCOQ1cwMuFohw5ioAi0DEpJXtuxzCzpYp8GT6Lb+Say0/gn8ExIFdoxWoBkSqsR\nIFXBln7dtqNhpU/sYRklVtVrUH7HawkJxd4M/KE3ehJblgyZbgVQPI3ap1/p8BUgDxznh7aIkTJy\nFwiZG0U0MXyRlq6olWohw3wlPsISSefm+5vOOga8giVgDr+vRRxbXRyRGDbnaG7DhzFfdrQfvA85\ncrehUMY0REpd7NcQkIhPdgL5d/FLYLEdcgV2jBYgmaNlBDjCT8FK3bKjWvKGw7nqNSi/47XYOd9/\nBpgiCHT0JLYsGTLdCuCj9ulXOtQB3lhPATcB86MYepeOw5ebuVJqmpauqDXSwvAVZiM9ap2bU7ZZ\nR+7G0egysL/Nq5A3OhGgOBxjzTk6N8ZQcrT1++Tc75gwQjkmKZEurtDQI9Ht9LKUaFvz+iXqoyWZ\no5VCfEpIVuiWHdXiS8RDUyUfMNxON2ABCds1AlS5oKMnsWXJkOlWZh+1Tz9PTJXYDU/yaA17WdOy\nu5kvXjqKkS8f324ras1XMPwWsElKJJ2bMwyz8p3bS3iaV1GNl1MdOH3b1/mZJk5/c45dgOgxR5ss\naokIr1DG5Eo0ampCDZ4EoA/vAlrN3c1+PjeP0wqZi1YKsfTG3hrP9emWHdWSv3Vvm9syB/UGAbud\nniMPBC36L0s0SR8deyRuWTPkdBMgiLqXfqUjYIp/uQo+xv4C0peEItfZivl7ka3qilob6yjw8TbJ\n6mhnc5HtrLmlJu9y2FbCLd9qAndjY01xBC9s7rhzS/zxYUd7ElKiXImdQCkJrERc7NcQkJgPjtTp\nOfwOu1VbnFbIXLSMIFoq4RACSInflUS0pJoJeQ2REoneIGDb8QHHqRNUc56ydHyhR+KWNUOsIQ8o\ngI+6l36lI+ZR/iXL+DaH+/Ez3syAp567WptfxsCMrqi1sSIlyoEl0s7mWiLdyfx992GkJzNzVeb5\nsICkJhRHzLA5R3TH52tGW9cSzfOoQimJKxEX+zUEJBAfYJd0gy3pV2tdOhuVzKKVEnWaeKgZQEr8\nVhLVQrxP8E9KBDxfM4hwpxewHPAtwatHSyTjCz2S7jJSEzyrugkQRO3TL3S8lC7wXGYCktbkvjfk\nWZRsYf/MyASSb+uKWnZ1/1xKpJ27cVC2WXdg6GoTt/1wrk2AZ5rsBi4qDq1Yc47dGx2yM0qSaWuJ\nHudZhVISVyIu9mvgMUcCiA8GKtKPSRdr0ecfZXukHaM1MotWEs830DtrdIvrlh0LC/hGsVui7EwY\nsO10HQkRtviNrkfyJ1ZM7n/MUFc3AXzUQfr58Jnh62eZhwsXtES89y5zli2icGWwjOQlTvLLao3w\nje5KVGeJtIPOVbZZB/hBqfL4YvMDwD3VfJmPtazgcC3WnKMBJIpIyYn88qvQEu3kRKCeFCaKEzn8\nDoM+DfAkug15lPK41kmsFS1G68g0WinEIyxRyUM63bIjWqp/sUcs56rXoPyOd+xn5XvBVJWLOnoS\nnWqGTLcCCJ5FrfQasEvPIDMhJZIbHd8bWuzm+Xc6Py1XkayoleJXoeXM2bOLv/qudJM61xKpJddi\nvsTjD6DTZImyZQwsKw7XYs058qsWAUeWrUTp1lNnz159HRHVK9RHlcRKJIt9GiY9ifogU5Gq/URy\nuHqL0RqZRSuFOCNXkYd0urljYT3Na4WnOFe9QcC64x37qbcQl4780sXxTI9Ep5oh0y0AT0oCLWqf\nfpeeeSnK8IS8LiQP4LyEmGV3V4HPomldMWsa/ImBe6PW8fubzCnbLH4IUsXNxcJFHAR2F/ldebCs\nODwea87R3Ig2eLFLS8cC7wHyyU1RgTDpVSSL/Ro8ifmcA/agUMEQHVdtMVqeONiNlhHgPJ9FK3XL\njmqp6HUmsZpeDd3vBFr6qDeDv1qw6ehJdKoZkt+h9hizj9qnX+kojwKQLGOWCdka3cvXGwxvRXIS\n/8Cpqq6oxavsVJsnR5g47WwuYaj1Gz6I8HjzpRJuReYt+aVkc9Fw5sv0CJocly+3MiZLmJo0EuAy\nn4lSIoUyEi2RLoqaQIMnUZ/EY40dE7ipMftXfhUOuAIzpOXyZYtWP2TI16NDK3Tr64Fq2SDltI+j\n6A0C1h2vRcIK22A1+rUI0tGT6FQzZLoVwEft06/p4c1SSjRQx3McZsdq+orxn/HtfL8Z+5dbUevV\nfZwid8+1knU63/HsAyWoVViYrdGFfsNjDxJlw4mHAcUZ/A79gqbHd8Hcfjv2aUeCY53TSBziQYVS\nkoEzS6f5ziqLfRoCEtlO82esCWzrdN5EP12XOUarZBqtRoDZHe0VunVHtSTGF5qAzlWvD1h3Ai2P\n8lzQouM/aMtbpo6exKaSIdOtAD5qS5sGrHS4uUrMNJGKAfaamj9aU7T/C3aD6WJ6MuuWRPn1hz+j\nrlOrrxPudWBvMF1MRS42W9MJf0bFXWuK6MGiSW/fAOsG08UjOhmfruFM/hkh/9Jbl5ZcF9Trgt5g\nuriOUny6hjN5x0WixO799h7NwCbgfwlODnJ00oAfAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left ( 9.869604401089358, \\quad 0.03947841760435722, \\quad -0.049163581851308535, \\quad 0.0630523848637\\right )$$" ], "text/plain": [ "(9.869604401089358, 0.03947841760435722, -0.049163581851308535, 0.063052384863\n", "7)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sigma_gl = g(mean_l+sigma_l,mean_t)-g(mean_l,mean_t)\n", "sigma_gt = g(mean_l,mean_t+sigma_t)- g(mean_l,mean_t)\n", "sigma_g = sp.sqrt(sigma_gl**2 + sigma_gt**2)\n", "g(mean_l,mean_t),sigma_gl,sigma_gt,sigma_g\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Result: $g = 9.87\\pm 0.06$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### \"Calculus approach\"\n", "\\begin{eqnarray*}\n", "(\\sigma_g)_L &=& \\left. \\frac{\\partial g}{\\partial L}\\right|_{\\bar{L},\\bar{T}}\\, \\sigma_L \\\\\n", "(\\sigma_g)_T &=& \\left. \\frac{\\partial g}{\\partial T}\\right|_{\\bar{L},\\bar{T}}\\, \\sigma_T \\\\\n", "\\sigma_g &=& \\sqrt{(\\sigma_g)_L^2 +(\\sigma_g)_T^2}\n", "\\end{eqnarray*}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Redefine $g(L,T)$ in terms of sympy floats so that we can do symbolic work. Will also use sym.sqrt() below.\n", "\n", "NOTE: sympy floats != regular python floats See http://docs.sympy.org/dev/gotchas.html\n", "Either\n", "sym.sqrt(sym.pi), or\n", "sp.sqrt(float(sym.py)), or\n", "sym.sqrt(sym.Float(sp.pi)), or\n", "sym.sqrt(sym.syimpify(sp.pi))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def g(l,t):\n", " return 4*l*sym.pi**2/t**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In sympy you must declare symbolic variable explicitly.
\n", "The method below allows you to control assumptions, as in
\n", "\n", "`n = sym.symbols('n',postive=True)`\n", "\n", "You can also import directly from a set of common symbols, e.g.,
\n", "\n", "`from sympy.abc import w`\n", "or\n", "`sym.var('z')`\n", "\n", "See http://docs.sympy.org/latest/gotchas.html#symbols" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "l, t = sym.symbols('l, t')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJUAAAA1BAMAAAC6kLSSAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiUSZq1TvELvdZiIy\nds1Wk1T5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD/0lEQVRIDaVXS4gTWRQ9lXQ+lY+pwY2oaAmC\nCIJpkRlxMUYQFKW1FyqKSrc/dKEYwc0sehKYAVEXRl2oNELAjQgzXUtxFp1B3doRQRHFro3owkVQ\n1MYP5X1VqfdJXvUnfRd555x77k29l7z3EmCeYWz4bZ4dRPkS7BRknmgUE47cwrBlNjd8FeM2zoma\nhQL2gIbKiLlhnVEKEY0jT9YBL75LihY+EOotgv+GdG85RDQW0TeITEtStPAGV40iwbQV8NQRrgOJ\nOrIu4q4k6WD+JFdvM5Q6HvA0axyGUUWqhmQz5BHj/XdhImbFGoT3BpzWTkRsqmzUMd4QihZZvNfj\nPftYvWn5vmuKe9fXxcCYInWTjEO9Jj3P+4gznsfyiW/sNeO/MuSH4b0FtoUsYhzBOxhLD91/6XCD\n/0RmlXMGdq/+bOO0InUTi3o9wpJ8Q6TG2EwnbCHQU5YweSInfb0W/sVMaiQa1Auo9kmpSpOUXxVf\n0kHuU36Ya4bT53ISgrNgvXIlMxRoXODSyzFJACrEtqdL/A1NJ6euJ3Pv3rLl8wDi9QK3AXmy5b4q\nvei5cKlQt0LRtLMfQiyPJ+ltywVbSFnqw/pJEXuLTLOwqimkxBeBBZqiTYskf0tK0PzSw8LA0MP+\ndcj0S5op9eXyIm87LqJP7nUBSFa5QQ9meyAtp1kX9S1CNTNDPvRhzMaEbgrcAKzAGolNA6nR0OA0\nefqY3+w5P62BJ2mCYw3OdCDpeS2d3q0VaMc43bKv5OkY8IoRSY1cqGJ5VK9fHr385z9NTZSUdLGp\nrE/mBg3nlT6lV5M17PIzf95k8R5gE2PRom8xtuqruCmwtlfUHG730hU9RE0nR2nUK2qOwEZod2Jk\nr5pu7Y2nbHKpL9CeEFIv5ZyktZ9sSMkAbsQGAvkaDrdTFbcN1EE9J2lj0zbqjH7/HkrXsbmdMa93\nWnyunpOFIoZsjY89lxwHZMKxek5WrK69nRqlpTrF/QGwOjijfzh0K0qfDu1tejQlWPp1U5GQ6uB+\nll338jlJh4TpqnVxxjvmmFEtAfudBvmcpM0Yb6lG0yI+PqiKOnaFLtOilKANlJGmTJmVl/9+fgrj\ndcmkhbFtUwPqOUk3U+6T6h0CRvGsrIoalmip52T2B5nUuxa0DP8/2aEp7pBoXZVz0r8bN6kmWoZZ\nhb+ukjNdI7LWlhSkZtqDobmilNF91qSMKX8YiLUw81qxfmuRYgOPSVamfpCGSz+IZhPLcFCxBct+\nVNbipRGZRuP16mYIfmOCbQYe2Yt1jqcF9+4q6QXBSrUHJTVn0v6/kZV+Uc65R7uA/9/Y32sHUZcO\npgjES0LsEd3hdQMc9QjEfz7E7R57hGXn/K/CT5nT+ekRLQimAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left ( - \\frac{8 l}{t^{3}} \\pi^{2}, \\quad \\frac{4 \\pi^{2}}{t^{2}}\\right )$$" ], "text/plain": [ "⎛ 2 2⎞\n", "⎜-8⋅π ⋅l 4⋅π ⎟\n", "⎜────────, ────⎟\n", "⎜ 3 2 ⎟\n", "⎝ t t ⎠" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sym.diff(g(l,t),t,1),sym.diff(g(l,t),l)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Evaluate the symbolic expressions at the values $\\bar{l}$ and $\\bar{t}$,\n", "and calculate the uncertainty." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAykAAAAUBAMAAAB7SG+3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQ\nqzLsm4+cAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJZElEQVRYCe1ZXYhcZxl+ZvZvZufM7PiTQOlF\nJkUooYhbg6JSzVgvKpqYkbqgSWNXqkWoJVtoiFI1o+CFaM3GijQl2rnQiiBkSFyDbJM9F21jaMuu\nwVijxp2belMkf2tM2o3j8z7vd+bMWbu92oZe9IN8P+95n+d5f+acM7MBcjW8Pd5KFdhowbwLeHb2\njwrr5CM7uM4+0gkW7Y6csGuVRpie9TNP2kV7Z1qIdh6vIjd7tEmrYVaMaOdtsmmVmwAirTy4c3YW\nLuwiBP/zzOxsOwTici/D+Qe/Wfk8HAPQuGJklMTvSjfPfsg8CZCSBzKxuYOcTZ6J0hF/X7ZKJ0SX\n0coqlY6dSPI3L69nz+i+IQZevit+ucmlV1LxC7Pu7o8B5SovfhzRj3GT7dBAuYlnWkN1t9gOn8Rz\ndm14PEx+/n4MvzIE/AWnarkDWIfCJQhNz8wYbEZfMoNWuQkg0ny3273uwi5SOAS8j8YFpwry98H5\ni91uG44BaFwxMkrik1Kugf01uhKgFOX2IkqvQpNnIiXx92WrEgjzRkrRp/C5kD9YGa9eapRciMEK\n92T3NWPrldT4hYkWMNdCVAfytE2j3KBboYZcHR/AwLhbbJefRole+IZ1hZPO0ZnF2Hf4LvBunAXO\nYF8TVyA0PTNjE/BbM2iVmwDiG2LLF1xY5xd3LQOf5QfGqYL8yGNw/sHjDwHCAGZcMTJK4pPSyFWM\nTQrgSnL7O3AQmiTiSuIX0CcrQRJdRiujNNDAZo9PlfHqpUb5KgZdxqdfYtKpiPiFKccoTgPH2LEq\nxuooXaBbaRzR5OhF7mTRbmABBdYp955xn8J5XwzfHQWewxIwh1+3oiUIY5L946PA/pgGrXITQKQl\nynYkHERGqNYGTjhVkHvnJ+D8bBYBhgHMuGJklMQvpdEnMF8XwJXk9lNgsaNJIq4k/jRbpROiy0pl\nlA5XeVGJ8cMTh3qmRvmGGHgZC6LqlVT8qnhxEuX/AqeA9cD8OEZ4wOjluFQbmOZOFu3GFjB6lVUY\nHPcpnE1bVw7fgSPRZQbT4b12CcJItG+6Dhxu8hzWwqUAECnwKxcOItYVhjLlVEGurQaQ36tmGPbu\n/7uSVSJ/EpqeYAQoRQ9ke8yuaJKIKzl/L1ul4xiLqm9klD7sF+wJbl3xeqbGxNeeorwcutITEb8w\ng8vqymALe4Bi3R6xHIuXjmDsM8d2uEW7+QZG/wNstK7YFM6mrSvl7m01/IT3Sos3TdvRRtU3on+z\nK1NAstLNASIFyGzCQcS7MsD71QJxkXxVDSCwfMuejlETE4x2SkaiEFbjD6FFDwSAKSVu2B0TuTuW\nSEhH/L1sFUOILhHRmlD4uvTKnhbNjE9l93r2jD1fxqCunL71i9z1RMSf9GCYH8pSA+/nfAGD1maU\nulsxfyeKTVm029BG5TVEU0xQk5+NPOzmlmI+vrC9hpu/wlYZ2qj6Ro7vtvMMOKzmJoD4KFpzYeeH\nd+UI71ejcpGTsK4YcCDO2efHMG7kIR1ZJfF7aLnf3RMAlmLiNmovXU4ScSXxp9l6OoZZMRIKrU8v\n1fC4x6eyq3pRzxh8FYO6cgc28J3eK6nKnvTgMAuVr+Or1LsLP+BTCnj0qWut+WUMHZJFuw0N60oJ\n7IomP6srulL46+4nMDhVmGuyVnx8C21c6cjxXlFXwko3AcTH6tPThJ3fuxJdoN2oXK6trhg/xwf5\nzzDBaKZkZJXEH0LD71sOMKXEbXiBQE4SCYkZf5qtp6PoEg1fEwqtT3djfCFW/n4zWD2jnjHxtRjU\nFWBgsk/Esw89eBB8E0zCMs3vfcXeK/kF7DvELyv5V2XRzm/sH1lXNPnZyH33Xoxci3Hrt+c6JDgY\nO5rbvtG7g9mVw1O8cNAB4gMeDsLO710ZqfMJZoFIpNDxrhgQuLcqTM9IUzKySs4fQiseEkApJm4T\nhuMkkZCY8afZKgZhEomwJhS+XuGX+RavWHysjNezZ3whybt4KOnK0EXPVkrO7z0YapCmckFd4RfR\nZZ6KVVSuDteRt8dZeVk7vgQLV6M2u6IJOovcdwd4HzTpvhi/A9jSTNC0pINvu/1NHrXKjYfFF4yU\nm/uCsPN7V4Z5IZV/HtYVAf/kbzBi3Eh8/+hX+mHgZ2i5KgauClC0FD0QjFgFbBqzHDXB+H+TZqsY\nHNMvw32/UvPL7EotJGZdUT1To3wVgwpXrvM1HqUiCb/1YI+BKxf0BAMGF3ia57/T5Wm7V8yi3QB/\noiwXzp1b/NnXbZrSWeTa2aOmXKP7NvCO3dJM0LSkYzODjnnUKjcetol0ChFjlnA4670yz3BS+UfP\nnbv2kvPz5X1vLIwbUxHt+pX+oHhNCWPL7IoAUvJA8D2wUjYpE0/M+AX0bBWDY95AKT5r90pIbF9s\nnoML9iMuGBWVYlDhinUM9Zc04Sem0LAP3+ikve3zB3C+Sqoip9srfJxPy+K7afAHPq+N+8QfW3am\ntu8O8CFZ3VStXMRj/DpTFYaembEJ/OnLoVVuAjhphfe3hOlgIurK+QY/MQokyBMlIO38DBqGg5YV\nI6MkPimRd/iiAK4kt0oDI1VNysTTcf5etoohiS6jlVE6z/eKx+eVUT1To3xDDCwc/xzCe9CzNSXx\new/4J6GH2I86ZliIrdGdsJ+AW5Gfwt9wqimLdryXTnVIMcaCafKzfSK0+wXfEXg4fqaGW1DgV2hD\nkyozhpvRz+1nnFa5CeB8BVbYhf2sruxndTLyl+H86+wjBMNwXLagMyOjJD4p5WvYPyWAK8lt/ezM\nn6HJM1E6zp9ma+k4hvH3j4xSuR09HvK3z6tXLzXKN8TAy7kGNlXJlZRU/MLkHpjdNcm/XbTxFK/P\nTLQAvv3/tXMHv0lM/CNYtHt+L48obble80nnXU9uq0G7ypmZFiHEjU7cTRZh7o8JSUd07Fsd+66j\nVW4COF+O3yldWOehs0ungZuaNPbJH+2edv7czjMxv2MbBma8n6e+kVESnyv9cuIj9DKAUpTb9m73\nCjR5JkrH+dNsFYMww1/r0+GPr/6cMLOrE/JXZVTP1ChfKAZdXnf8HnKlIuI3zCD//jfJxx+rZW1b\n41F4EzhfP8Qbp4TvvH4Eb4J1PR9ijbXnLa095SqMN04J7VVCWHsz/zqJ29ee9uTaU67CeOOUoqlV\nQlhzs/0l3/7Xa61Hba0JV+W7cUr5VWNY6wsD9vzP1Ti9Pd46FdgI/A/Ydu2j8eQ/7QAAAABJRU5E\nrkJggg==\n", "text/latex": [ "$$\\left ( 9.86960440108936, \\quad 0.0394784176043574, \\quad 0.0493480220054468, \\quad 0.0631963031544892\\right )$$" ], "text/plain": [ "(9.86960440108936, 0.0394784176043574, 0.0493480220054468, 0.0631963031544892)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sigma_gt = sym.diff(g(l,t),l).evalf(subs={l:mean_l,t:mean_t})*sigma_t\n", "sigma_gl = sym.diff(g(l,t),l).evalf(subs={l:mean_l,t:mean_t})*sigma_l\n", "sigma_g = sym.sqrt(sigma_gl**2 + sigma_gt**2).evalf()\n", "g(mean_l,mean_t).evalf(),sigma_gl,sigma_gt,sigma_g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Result: $g = 9.87 \\pm 0.06$\n", "\n", "This agrees with previous result." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Monte Carlo technique" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we are going to get the uncertainty in a different way, using simulated data. Our propagation of uncertainty formulas assume that the uncertainties are standard deviations of normal distribution. We can \"redo\" the experment, with new simulated \"measurements\" sampled from the assumed distributions.\n", "\n", "First we return to our definition of $g(L,T)$ that returns a regular python float." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def g(l,t):\n", " return 4*sp.pi**2*l/t**2" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 9.81457334, 9.87911409, 9.85254789, ..., 9.90789589,\n", " 9.9116262 , 9.77532495])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nEx = 10000 # Number of simulated experiments\n", "l = sp.random.normal(mean_l,sigma_l,nEx) # Pick values of l from normal distribution\n", "t = sp.random.normal(mean_t,sigma_t,nEx) # Pick values of t from normal distribution\n", "gg = g(l,t) # Calculate nEx values of g using random l's and t,s\n", "gg" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAAUBAMAAAAKIynxAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQ\nqzLsm4+cAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE6klEQVRIDc1WXWgcVRg9O/ubndkfLC0UHzJK\noPSpq0VREVzrQ6WmzQouaNPYgD4UtLhCQ5GKWR8URGq3KlJLpPtifXSJ1CJLmnlILaEtuwZjCRqz\nL/WlSGIbNY2J6/m+O7vdgHmQIPSDuffO+b5z5uy9c+8sEHBx94aVprdNwMXqd2ry0rH97KvHGi3k\nett7Mgenf0zq0QJDbyafB+zzF9pFrcHFUYOZXvTurT4C/DJdrVYMQVm7vetFQ9F0IL+z4deIkeTh\n/moVYuxJOB9iqz47h0QRE+VwtoW8xHxsRFQiGUy6gZMyFPA9D+hqNitw9uA5AdfE0/hW9bQXvUAO\nJ1zsaDabdSUY1pnmX6RRyqSvwl42NRAjFqtXMQFYvCshkWNtzEUgi4cQzPhI9CPg6sCKPP2NDGaA\naY4IOtNzHhAaex0I5rBT8p1hlWBnCZhe9KJLSA3iWSBhCIb1zLW0kTLpn4BPTY0aCfNH1mEXEEkj\nlYW9QEU7A2cw/htHPnLPUxxHxWGgJ4NFYJxDBYc9eRrjS50sGd2JYB0xIWmvevHTqGVRAS4YgmHV\nlUIpk/4EmGtojRqxaagBp4QtQC2D6N+sjt/ybDdY4shHKm2Hdijj3AKGG4CCbYeP6lPWNqk64kuE\ntFc93nCV+YQClGBYbYcm3efRodaoEYJf8NqDIb5QnMNl3mDu5ihS+87v9xEr3XZ4XyiDjzmHZRhQ\nHd4/1MDijaGyUDujlkP8TwLaqx7gvCIVQRiCYU1te5GYSPlpvMYxa9QIuwyvc3iQk7mA0E3ewG72\norYLXUWDXELLoVOgQy5xnwsDimzQCyw7iy5OCbUzuitIyh7QXvUQ+PqAVIxCCT7rcXSXjUM/HRfW\nKC8xwtZl04NDbHfjfVkVHD97u1xbQXjEIJW2Qxt0GCrExoswoP5w4GGn6eEFT7gd0Z3zHUpv9IBv\nypypBSjhSosVHPTnUNOIcN1ZY4xwDTjCuzjN1jp6Q95Dq47hEe45a1mRWKPt8ANxiG1vjTd80Hd4\nMP0HsKNMbmesXWXV46s0wk2XBZTQYoW5LY2UpJHnJTVqBDjCW8yrQ+7LFd50pZFcimRhyZInVi6j\n5dCpqEO+H54Piuz38l6+TIcuyzuDOyTm75SY0QukESTCQxVK0CaR1defUn46mqOK1KgRc/DSoawy\nz7Y6G5nVqURJ5lCQ47Ozt6/paRObnZ37rEBwL3xQHHLnHPRm/mUOgxlY8ou1V73Uijqs8SlK0IYb\nNMwySvnpd8Sq1KgROKvUwFnZKdZJzKd508XmgSTfw1IL4YltzkN0ZbA9nZSzEgJSlgc/J2Se76GA\nncGTOpwloL3qkR0hd54UJWgT5oSxjFImzQ9rNK01agTJ30W0h7sZ0V5nFzcZ4r2wCvgRk0VFmOYR\n6DtMZXDEm3B5r6A43Cw/LVFxTslpvCYOYbJBPa7PZEP1LBcnCjwT6VAJ2vBbtz2tDk16S/XcD6bG\nGImpw37OIs+cfBk4DPzav597PP+zj+Cr5hTCM4tTBJ9Yde08s1Bw4MxefiH7pz1WDjQQeVUSd+Ly\nUUpQT3vV+zz/GNNbi2yEYJrNYwcAldJ0X7PJHaQ1aiSgx9g++eptPN7euMQ6CoFB+eew8ahsXGId\nBbvAvbFO7j/ADmX+p5ig7qaNa1sbl1hPQeYv4LK5W8NK4x80cfIa4awungAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\left ( 9.86941494156, \\quad 0.0637665132787\\right )$$" ], "text/plain": [ "(9.86941494156, 0.0637665132787)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sp.mean(gg),sp.std(gg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Result: $g = 9.87 \\pm 0.06$ " ] }, { "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": { "collapsed": true }, "outputs": [], "source": [ "%load_ext version_information\n", "\n", "#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.5.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]" }, { "module": "IPython", "version": "4.2.0" }, { "module": "OS", "version": "Linux 3.10.0 327.el7.x86_64 x86_64 with redhat 7.3 Maipo" }, { "module": "numpy", "version": "1.11.0" }, { "module": "sympy", "version": "1.0" } ] }, "text/html": [ "
SoftwareVersion
Python3.5.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython4.2.0
OSLinux 3.10.0 327.el7.x86_64 x86_64 with redhat 7.3 Maipo
numpy1.11.0
sympy1.0
Thu Jun 22 10:52:15 2017 EDT
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.5.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] \\\\ \\hline\n", "IPython & 4.2.0 \\\\ \\hline\n", "OS & Linux 3.10.0 327.el7.x86\\_64 x86\\_64 with redhat 7.3 Maipo \\\\ \\hline\n", "numpy & 1.11.0 \\\\ \\hline\n", "sympy & 1.0 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Thu Jun 22 10:52:15 2017 EDT} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.5.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]\n", "IPython 4.2.0\n", "OS Linux 3.10.0 327.el7.x86_64 x86_64 with redhat 7.3 Maipo\n", "numpy 1.11.0\n", "sympy 1.0\n", "Thu Jun 22 10:52:15 2017 EDT" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%version_information numpy, sympy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [Root]", "language": "python", "name": "Python [Root]" }, "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }