uncertainties

environment.yml

name: example-environment

dependencies:
  - python=3.7
  - numpy
  - matplotlib
  - ipywidgets

Uncertainties_and_Monte_Carlo-compiled.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Übung zu Unsicherheiten\n",
    "\n",
    "#### 1) Messwerte mit Unsicherheiten\n",
    "Gegeben sind die Messwerte mit Messunsicherheiten eines Zylinders aus unbekanntem Material. Ermittelt werden soll die Dichte $\\rho$ des Zylinders und die zugehörigen Unsicherheiten $\\Delta\\rho$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# load essential libraries\n",
    "import numpy as np \n",
    "\n",
    "\n",
    "m = 2670 #g\n",
    "D = 52.1 #mm\n",
    "H = 160. #mm\n",
    "\n",
    "dm = 2   #g\n",
    "dD = 0.1 #mm\n",
    "dH = 0.5 #mm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "#### 2) Erwartungswert $\\overline{\\rho}$ der Dichte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rho [g/mm^3] 0.00782755180785807\n",
      "rho [g/cm^3] 7.827551807858071\n"
     ]
    }
   ],
   "source": [
    "# expectation value\n",
    "rho = 4 * m / (np.pi * D**2 * H)\n",
    "print('rho [g/mm^3]', rho)\n",
    "print('rho [g/cm^3]', rho * 1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 3) Messunsicherheit $\\Delta \\rho$ der Dichte "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "drho [g/mm^3] 3.91869545965466e-05\n",
      "drho [g/cm^3] 0.0391869545965466\n"
     ]
    }
   ],
   "source": [
    "# uncertainty\n",
    "drho = rho * ((dm / m)**2 + (2 * dD / D)**2 + (dH / H)**2)**0.5\n",
    "print('drho [g/mm^3]', drho)\n",
    "print('drho [g/cm^3]', drho * 1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "#### 4) Relative Unsicherheit der Dichte $\\frac{\\Delta \\rho}{\\overline{\\rho}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relative [.] 0.0050062849225995356\n",
      "relative [%] 0.5006284922599535\n"
     ]
    }
   ],
   "source": [
    "# relative uncertainty\n",
    "rel = drho / rho\n",
    "print('relative [.]', rel)\n",
    "print('relative [%]', rel * 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 5) Monte Carlo Betrachtung der Unsicherheiten $\\left(\\Delta m = 2g,~\\Delta D = 0.1 mm,~\\Delta H = 0.5mm\\right)$\n",
    "\n",
    "Es wird nun \"numerisch\" 1 Mio. mal gemessen. Welchen Wert erhalten wir dann für die Unsicherheit $\\Delta\\rho_{MonteCarlo}$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monte Carlo: rho [g/mm^3]: 0.007827747464370816 +- 3.9194483550048165e-05\n",
      "\n",
      "Monte Carlo: rho [g/cm^3]: 7.827747464370816 +- 0.03919448355004816\n",
      "Gauss Error: rho [g/cm^3]: 7.827551807858071 +- 0.0391869545965466\n"
     ]
    }
   ],
   "source": [
    "# Monte Carlo\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(seed=2)\n",
    "\n",
    "# number of measurements\n",
    "N = 1000000\n",
    "\n",
    "# values for m, D, H after N measurements\n",
    "m_MC = np.random.normal(m, dm, N)\n",
    "D_MC = np.random.normal(D, dD, N)\n",
    "H_MC = np.random.normal(H, dH, N)\n",
    "\n",
    "# calculated possible N values for rho \n",
    "rho_MC = 4 * m_MC / (np.pi * D_MC**2 * H_MC)\n",
    "\n",
    "# mean and std of rho\n",
    "rho_MC_mean = np.mean(rho_MC)\n",
    "rho_MC_std = np.std(rho_MC)\n",
    "\n",
    "print('Monte Carlo: rho [g/mm^3]:', rho_MC_mean, '+-', rho_MC_std)\n",
    "print('')\n",
    "print('Monte Carlo: rho [g/cm^3]:', rho_MC_mean * 1e3, '+-', rho_MC_std * 1e3)\n",
    "print('Gauss Error: rho [g/cm^3]:', rho * 1e3, '+-', drho * 1e3)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Fazit: Hier ändern sich Erwartungswert und Unsicherheit der Dichte nur geringfügig. Der Unterschied verschwindet hinter den gültigen Stellen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 5) Monte Carlo Betrachtung der Unsicherheiten $\\left(\\Delta m = [2;20]g,~\\Delta D = [0.1;10]mm,~\\Delta H = [0.5;10]mm\\right)$\n",
    "\n",
    "Es wird nun \"numerisch\" 1 Mio. mal gemessen. Welchen Wert erhalten wir dann für die Unsicherheit $\\Delta\\rho_{MonteCarlo}$ wenn es größere Ausgangsunsicherheiten $\\Delta m$, $\\Delta D$ oder $\\Delta H$ gibt?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Monte Carlo widget\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(seed=2)\n",
    "\n",
    "# for general plotting\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# interactive plots\n",
    "from ipywidgets import interactive\n",
    "import ipywidgets as widgets\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.dpi']= 200\n",
    "\n",
    "def update(dm = widgets.IntSlider(min=2., max=20., step=6, value=2.,     description='dm [g]:'), \n",
    "           dD = widgets.FloatSlider(min=0.1, max=10., step=2, value=0.1, description='dD [mm]:'), \n",
    "           dH = widgets.FloatSlider(min=0.5, max=10., step=2, value=0.5, description='dH [mm]:')):\n",
    "    \n",
    "    # Gaussian Error Propagation\n",
    "    rho_Gauss = 4 * m / (np.pi * D**2 * H)\n",
    "    drho_Gauss = rho_Gauss * ((dm / m)**2 + (2 * dD / D)**2 + (dH / H)**2)**0.5\n",
    "\n",
    "    \n",
    "    # Monte Carlo Error Propagation\n",
    "    N = 1000000\n",
    "\n",
    "    m_MC = np.random.normal(m, dm, N)\n",
    "    D_MC = np.random.normal(D, dD, N)\n",
    "    H_MC = np.random.normal(H, dH, N)\n",
    "    \n",
    "    rho_MC = 4 * m_MC / (np.pi * D_MC**2 * H_MC)\n",
    "    \n",
    "    # mean and standard deviation\n",
    "    rho_MC_mean = np.mean(rho_MC)\n",
    "    rho_MC_std = np.std(rho_MC)\n",
    "    \n",
    "    # results\n",
    "    #print('Monte Carlo: rho [g/mm^3]:', '%1.5f'%rho_MC_mean, '+-', '%1.5f'%rho_MC_std)\n",
    "    #print('')\n",
    "    print('Monte Carlo: rho [g/cm^3]:', '%1.2f'%(rho_MC_mean * 1e3), '+-', '%1.2f'%(rho_MC_std * 1e3))\n",
    "    print('Gauss Error: rho [g/cm^3]:', '%1.2f'%(rho_Gauss * 1e3), '+-', '%1.2f'%(drho_Gauss * 1e3))\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,4))\n",
    "    \n",
    "    # Diagram: histogram, bins = 500\n",
    "    h = plt.hist(rho_MC * 1e3, bins=500, density=True, label='MC: dm=' + str(dm) + 'g, dD=' + str(dD) + 'mm, dH=' + str(dH) + 'mm')   \n",
    "    \n",
    "    # for plot of gauss function using dF and F_mean from above\n",
    "    xmin = np.percentile(rho_MC * 1e3, 0., axis=0)\n",
    "    xmax = np.percentile(rho_MC * 1e3, 99.9, axis=0)\n",
    "    xfine = np.linspace(xmin, xmax, 100)\n",
    "    gauss = 1. / (2 * np.pi * (drho_Gauss*1e3)**2)**0.5 * np.exp(-(xfine - rho_Gauss*1e3)**2/(2 * (drho_Gauss*1e3) **2))\n",
    "    plt.plot(xfine, gauss, 'r-', alpha=0.5, label=r'Gauss $\\frac{1}{\\sqrt{2 \\pi \\sigma^2}}\\ \\exp{-\\frac{(x-\\overline{\\rho})^2}{2\\sigma^2}}$')\n",
    "    \n",
    "    \n",
    "    plt.xlabel(r'Dichte $\\rho (g/cm^3)$')\n",
    "    plt.ylabel('genormte Anzahl')\n",
    "    plt.legend(loc='best',prop={'size': 12})\n",
    "    plt.xlim(xmin, xmax)\n",
    "    \n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### 6) Interaktives Diagramm\n",
    "\n",
    "Frage 1: Wie verändert sich das Histogramm für große Werte in $\\Delta m$, $\\Delta D$, $\\Delta H$?\n",
    "\n",
    "Frage 2: Wie verhält es sich mit Kombinationen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb10181a41234810b4d6904610c612a0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=2, description='dm [g]:', max=20, min=2, step=6), FloatSlider(value=0.1,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Frage 1: Wie verändert sich das Histogramm für große Werte in dm, dD, dH?\n",
    "#Frage 2: Wie verhält es sich mit Kombinationen?\n",
    "\n",
    "interactive_plot = interactive(update)\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

Uncertainties_and_Monte_Carlo-empty.ipynb

Uncertainties_and_Monte_Carlo.ipynb