Skip to content

Instantly share code, notes, and snippets.

@koepferl
Last active July 18, 2021 09:10
Show Gist options
  • Save koepferl/a4ef10e0e3eb21a8b489c1fc5d119d70 to your computer and use it in GitHub Desktop.
Save koepferl/a4ef10e0e3eb21a8b489c1fc5d119d70 to your computer and use it in GitHub Desktop.
uncertainties
name: example-environment
dependencies:
- python=3.7
- numpy
- matplotlib
- ipywidgets
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
{
"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": null,
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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: ???"
]
},
{
"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": null,
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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
}
Display the source blob
Display the rendered blob
Raw
{
"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": null,
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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: ???"
]
},
{
"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": null,
"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": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment