Skip to content

Instantly share code, notes, and snippets.

@jradavenport
Created September 10, 2021 19:24
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jradavenport/d3b5163f2e701f0729c504f48a77220b to your computer and use it in GitHub Desktop.
Save jradavenport/d3b5163f2e701f0729c504f48a77220b to your computer and use it in GitHub Desktop.
Example of a simple Monte Carlo simulation to estimate errors
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "fa6400ff",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib\n",
"\n",
"matplotlib.rcParams.update({'font.size':18})\n",
"matplotlib.rcParams.update({'font.family':'serif'})\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "30766a1e",
"metadata": {},
"outputs": [],
"source": [
"# let's do a \"monte carlo\" (aka random number) simulation!\n",
"# NOTE: this is NOT a Markov Chain Monte Carlo (MCMC) example\n",
"\n",
"# Here's our pretend distances\n",
"dist_1 = 12.3\n",
"dist_2 = 543.21\n",
"\n",
"# pretend errors\n",
"err_1 = 0.1\n",
"err_2 = 0.5\n",
"\n",
"# Here's our make-believe measurment:\n",
"meas = np.sqrt(dist_1**2 + dist_2**2)\n",
"# this is obviously very simple, but lets pretend it's the complicated ellipsoid thing"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "2f092dce",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f9d53606130>]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEACAYAAABPiSrXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXKklEQVR4nO3dfZBc1Xnn8e8DBiNhWbLRYJBgrbAggmoBv8hZCGveElHxym8Ufos3sSsQtLDZrHfDEkQ5OM7LrgTLOrHjdVGyjW0cvwVKJDiYBLOgQIRJeaiKwCUisGVs0EsYSUjgIN6sZ/+4d7xN06M53dN3ejTz/VR1Hc25p0+f6cv0j3vPvacjM5EkqcRBgx6AJOnAYWhIkooZGpKkYoaGJKmYoSFJKvaKQQ+gH+bPn5+LFi0a9DAk6YBy//3378jMoW6eMy1CY9GiRQwPDw96GJJ0QImIH3X7HE9PSZKKGRqSpGKGhiSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqZmhIkopNizvCpQPNGavvZMvuvcXtF86bxfqV5zY4IqmMoSENwJbde3l09fLi9otW3trgaKRynp6SJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxboKjYj4tYjYHRFf3E+bIyLicxGxLSKeiIh7IuLs/bR/T0TcX7d9LCKujYjZ3YxLkjQ5ikIjIuZHxE3A/wDm7qfdHOBu4CTgZOAo4FvAHRGxrEP7C4G/AD6RmUcCZwLvAv46Ig7u8neRJDWs9EjjBuAHwHnjtLscWAJcnJk7MnNfZq4CNgDXRcTPvvQpIl4DfAK4KTO/ApCZPwQuA84BPtTVbyJJalxpaKzIzCuA58ZqEBEBXARsysyNbZvXAsdRhcGo91Edtaxta3sbsBf4zcKxSZImSVFoZObjBc2OBxYAD3TYtqEuz2qpO7MuX9I+M18ANgKnRcQrS8YnSZoc/fyO8MV1ua3Dtq11eUIX7d9MdXTyUF9GJ00Ft62s//HWgQ5D6lU/L7kdnSB/psO20bp5E2j/EhGxIiKGI2J4ZGSki2FKA7T9weohHaCauE8jG25fPSlzTWYuzcylQ0NDvXQhSepSP0NjT10e3mHb7LY2vbSXJA1YP0Pj4bo8usO2BXX5SBft9wGb+zM0SVI/9DM0vk81gX1Kh22jdeta6u5u2wZARBxCdXPgfZn5bB/HJ0maoL6FRmYmcD1wYkQsadt8AdVRw10tdTcCTwHnt7V9G9Xpqc/3a2ySpP7o90T4NVSXyK6plx45KCKuBE4FLs3MF0cbZuYu4HeA90TEfwCIiEXAtVTh8qU+j02SNEGla099MCK2A9+tq94fEdsjov3GvKepbtrbBDwIbAeWA8sy8/b2fjPz88AHgMsi4gng74FvAm/PzJ/2+DtJkhpSdHNfZn4V+Gph2x1Uy4kUycwbqU5VSZKmOL9PQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQV6+f3aUgz1hmr72TL7r3jtvv6oTsBWDhvVtNDkhphaEh9sGX3Xh5dvXz8hl/4DADrf+PchkckNcPTU5KkYoaGJKmYoSFJKuachtRB6cT2KCe2NVMYGlIHxRPb0gxjaEjTVC9HS+tXelWX9s/QkKapbo+WFq28tcHRaLowNKQDwMJ5s7r+UHeeRU0wNKQDgKeNNFV4ya0kqZihIUkqZmhIkooZGpKkYoaGJKlYI1dPRcShwKXAhcBRwAvAA8AfZOY/tLU9ArgaWA4cDGwCrsrMdU2MTTOTy4JI/dHUJbdfAt4L/CpwE3A4sAZYHxG/kpl3AETEHOBuYDdwMrALuAK4IyLelpnfbmh8mmFcFkTqj76fnoqIY4EPADdn5o1Z+QnVkUcAv9vS/HJgCXBxZu7IzH2ZuQrYAFwXEd5HIklTSBNzGgvrcnNrZWbuAXYAxwJERAAXAZsyc2NbH2uB44BzGhifJKlHTYTGI8DzwOLWyoh4LTAf+Ke66nhgAdVcR7sNdXlWA+OTJPWo76GRmTupTkG9PSI+FBGHRsQQ1ZzGLuBjddPRUNnWoZutdXlCv8cnSepdI5fcZuYngUuAPwWeBp4Ajgb+XWY+WDebW5fPdOhitG7eWK8RESsiYjgihkdGRvoxbEnSOJqYCD84Ir4OXAt8GJhDFRg/orp6qv2UU/byOpm5JjOXZubSoaGhCY1ZklSmiSONC4H3A3+Umd/MzOczc3td/zxwQ30fx566/eEd+phdl3s6bJMkDUgTobGsLu9prczMZ4Fh4F8BPw88XG86ukMfC+rykQbGJ0nqUROh8aq67HTaaV9Lm+9TTXif0qHdaN26vo5MkjQhTYTGcF2e3loZEYcAbwaeA76XmQlcD5wYEUva+riA6j6PuxoYnySpR02Exp8BjwNXRcRZUZkDfBo4BliVmU/Vba8BHgLWRMT8iDgoIq4ETgUuzcwXGxifJKlHTdynMQL8AtVd3TcAT1KFyMnAhzPzD1raPg2cSbVI4YPAdqqFC5dl5u39HpskaWIaWdspM7dR3adR0nYH1XIikqQpzu/TkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxV7RVMcRMQf4KHA+MJcqoB4CPpuZf97S7gjgamA5cDCwCbgqM9c1NTZNLWesvpMtu/cWt184bxbrV57b4IgkjaWR0IiI+cA9wL3A6Zm5KyJOAP4GeDfw53W7OcDdwG7gZGAXcAVwR0S8LTO/3cT4NLVs2b2XR1cvL26/aOWtDY5G0v40dXrq/wDPABdn5i6AzHwE+D3gBy3tLgeW1O12ZOa+zFwFbACui4jGjoQkSd3re2hExM8B7wO+mJn7Wrdl5tcy84q6XQAXAZsyc2NbN2uB44Bz+j0+SVLvmvg/+XfU5fA47Y4HFgA3dti2oS7PAjxFpZdYOG9W16eoFs6b1dBopJmlidA4dfQfEbEGOA+YTTUJ/ieZ+Zf15sV1ua1DH1vr8oQGxqcDnJPg0uA0Mafxurq8Bfg+1QT3ScDDwM0RcUm9fW5dPtOhj9G6eWO9SESsiIjhiBgeGRmZ8KAlSeNrIjQOq8t/zMxrMvPpzNwJ/CfgceDqiHhVS/vs5UUyc01mLs3MpUNDQxMcsiSpRBOhMXqUcGdrZWa+UNe9Gjgd2FNvOrxDH7Prck+HbZKkAWkiNH5clzs7bHuiLoeoTlcBHN2h3YK6fKSP45IkTVAToXFfXR7ZYdvoeaQRqvmOrcApHdqN1q3r68gkSRPSRGjcAjxJddXUz0TEwVSX0D4J3JuZCVwPnBgRS9r6uADYDNzVwPgkST3qe2hk5lPAfwXeGhH/PSJeGRGzgU8Arwc+kpn/Uje/hupS3DURMT8iDoqIK6ku2700M1/s9/gkSb1rZBmRzLwBeBfVEcN2YAvVpbfnZeaXW9o9DZxJtUjhg3Xb5cCyzLy9ibFJknrX2NpOmXkL1amq8drtoFpORJI0xfl9GpKkYoaGJKmYS49L6olfnjUzGRqSeuKXZ81Mnp6SJBUzNCRJxTw9JQno/sut/GKrmcnQkAT45VYq4+kpSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVKzx0IiIhRGxJyJyjO1HRMTnImJbRDwREfdExNlNj0uS1L3JONL4DPDqThsiYg5wN3AScDJwFPAt4I6IWDYJY5MkdaHR0IiI91KFwXfHaHI5sAS4ODN3ZOa+zFwFbACuiwi/JEqSppDGQiMi5gGfAi4FnumwPYCLgE2ZubFt81rgOOCcpsYnSepek0ca1wJ3ZObfjrH9eGAB8ECHbRvq8qwmBiZJ6k0jp3/qiex3Uc1VjGVxXW7rsG1rXZ7Qv1FJkiaq70caEXEYsAa4LDN37Kfp3Lp82amrlrp5+3mdFRExHBHDIyMjPY1VktSdJk5PfQz4UWbeUNi+46W44z4pc01mLs3MpUNDQ710IUnqUl9PT0XEKcBvAW8saL6nLg/vsG12WxsdQM5YfSdbdu8tbr9w3qwGRyOpn/o9p7G8Lu+tLo76mdcCRMT2+udrgb+q/310h34W1OUjfR6fJsGW3Xt5dPXy8RtKOuD0NTTqeyxWtddHxDrgrMw8qqUuqCa8T+nQ1Wjdun6OT5I0MQNbeyozE7geODEilrRtvgDYDNw16QOTJI1p0AsWXgM8BKyJiPkRcVBEXAmcClyamS8OdniSpFZNLyPynXoe4xfrn7fXj6MBMvNp4ExgE/AgsJ1qXmRZZt7e5NgkSd1rdG2nzDy9oM0OquVEJElT3KBPT0mSDiCGhiSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqZmhIkooZGpKkYo3eES5JoxbOm8Wilbc22v/6lec21r8qhoakSdH0B3qTgaT/z9DQuPwmPkmjDA2Ny2/ikzTKiXBJUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSsb6HRkS8OiJ+OyLui4idEbEnIr4XEb8bEYd0aH9ERHwuIrZFxBMRcU9EnN3vcUmSJq6JI42vA9cAVwNDwHzgT4BVwNrWhhExB7gbOAk4GTgK+BZwR0Qsa2BskqQJaCI0DgL+NDNvzsx9mflCZn4e+Abw9rYwuBxYAlycmTvq9quADcB1EeHaWJI0hTQRGl8Fvtyh/jt1+RaAiAjgImBTZm5sa7sWOA44p4HxSZJ61PfQyMwbOoQAwKF1+WRdHg8sAB7o0HZDXZ7V5+FJkiZgMq+eegvwIvDN+ufFdbmtQ9utdXlC04OSJJWblNCIiGOBdwKfyszH6+q5dflMh6eM1s3bT58rImI4IoZHRkb6NlZJ0tgaD4167uI6YCPw0Q5Nspd+M3NNZi7NzKVDQ0MTGaIkqdBkXJ30v6iukDo9M59tqd9Tl4d3eM7stjaSpCmg0dCIiJXArwJnZub2ts0P1+XRHZ66oC4faWpskqTuNXZ6KiJ+G/hvwC9n5g/quiMiYlHd5PtUE96ndHj6aN26psYnSepeI6ERERcCvw+cl5kPtWx6B/BxgMxM4HrgxIhY0tbFBcBm4K4mxidJ6k3fT09FxAeAzwK3AudHxPktm98A7G75+RqqgFgTEe8GdgFXAKcC/z4zX+z3+CRJvWtiTmMl1RHMO+pHuy+N/iMzn46IM6nWqXoQOJhqrmNZZnqUIUlTTN9DIzPf0GX7HVTLiUiSpjgXBJyBzlh9J1t27y1uv3DerAZHI/XHwnmzWLTy1q6fs37luQ2NaHoyNGagLbv38ujq5YMehtRXvXz4dxsyMjSmBY8cJE0WQ2Ma8MhB0mQxNCTNWN3OgzgHYmhImsG6DQDnQCb3+zQkSQc4Q0OSVMzQkCQVc05Dkgo5cW5oSFIxJ849PSVJ6oKhIUkq5umpKabbJUHAZUEkTR5DY4pxSRBJU5mnpyRJxQwNSVIxQ0OSVMzQkCQVMzQkScW8eqpLvXxL3nRbRkBSmV6+t7zb/if788XQ6FK3l8ROx2UEJJVp+gN9EJ8vUyI0IuI9wJXAscBzwDeAj2XmMwMdWB/0ssCZJE1VAw+NiLgQ+Bzw65n5lYj4OeB24E0RsSwzfzrYEU6Mp6YkTScDnQiPiNcAnwBuysyvAGTmD4HLgHOADw1weJKkNoO+eup9wFxgbVv9bcBe4DcnfUSSpDEN+vTUmXX5QGtlZr4QERuB0yLilZn5XFMD6OVqKEmaqQYdGovrcluHbVuBNwPHAQ81NQAXCJSkcoMOjbl12ekqqdG6eZ2eGBErgBX1jz+JiE29DiKu7lg9H9jRa59qxPTZJxfGoEfQD9NnfxzAWj6/etkfr+/29QYdGqOy6ydkrgHWNDAWACJiODOXNtW/uuc+mVrcH1PLZO2PQU+E76nLwztsm93WRpI0YIMOjYfr8ugO2xYA+4DNkzccSdL+DDo07q7LU1orI+IQ4CTgvsx8dtJHVWns1Jd65j6ZWtwfU8uk7I/I7Ho6oX8vHvFa4IfA32Tm+1vq3wn8FXBRZl4/qPFJkl5qoKEBEBEXUSXkh+plRBZRLSPyOHDALyMiSdPJwEMDICLeS7Vg4THA81QLFl41HRYslKTpZNBzGgBk5o2Z+abMPDIzj8nMy/oRGBHxaxGxOyK+OMb20yLiSxHxWETsjIiRiFgbEW/s4bUuiYiNEfFERPwgIj4aEQdP9HeYTiZjf0TE6+r3/h/rPvZExHcj4uKImBY3R/TLZP59tPT5xoh4MSIe7bWP6WqSP6+OiojPRMTmup9tEfG3EfG2cZ+cmdPuQXWTy03Aj6juAflihza/UG+7BTiqrns9cC/wLHBGF6/3h1Q3I/5y/fMbgBHghkG/F1PhMZn7A9gI/DNwFhDAYcDH6r4/Pej3Yio8Jvvvo6XPg4H7634fHfT7MFUeA/i8Oh7YDvwx8Kq67i3ATuDacZ8/6DesoZ3wLeBq4MT97ITT6jd7blv94vo56wtfazHwInBNW/1/qfs5Z9Dvx6Afk7w//gn4jx3q19f9LB70+zHox2Tuj7bnXg58B/ixoTGY/VH/j9R3gJs7bLsM+Mi4fQz6DWtoJxxTl4v2sxOOAf7zGM/fCTxf+Fqr6tc4ra1+YV3/lUG/H4N+TPL++AiwoEP9/65f+4ODfj8G/ZjM/dHynOOA3cC/AR41NAazP4Cz69d4R6/jnSrLiPRVZj5e2ObTY2w+BHiy8OXGWql3S0TspDpNMqNN5v7IzE+OsenQuizdr9PWJP99jLoO+LPM/J5TSy81yfvjnXU5XNj+ZabERPhUEhE/D8zh5d/xMZbFwFPZeeJ+K7AwImZ32KYCPeyPsbwF2AX83YQHNYP1sj8i4sNU59//uKlxzVQ97I9TqY405kTE1+pJ9X+OiNsi4q0lHRgaL/dbVOtd/c/C9nPpvEovLfVzx9iu8XW7P14mIpYC/xb4+BjhrnJd7Y+IGAKuBVZkg9+LM4N1+/fxOqrQ+Duq++FOAN4I/BS4KyLePl4HhkaLiPhF4BKq/8Af6+Kpg7/ZZRqawP5o7eMwqptHb2Xsw3sV6HF/fBL4y8z0CK/Petwfh1F97v91Zn4hM5/NzK1UX639HPCp8TqYlnMavYiI44CbqW4q/IsunrqHzqv0giv19mwC+6O1jwBuoFr48gNZzwSqe73sj/qa/3OAJU2ObSaawN/H6JH2na2VmbkrIoaBMyPihMx8ZKwOPNIAImIB8G3gC5m5usunPwy8eox5iwXAVk+JdGeC+6PVZ6g+sH4lM3/Sl8HNQBPYH++i+j/bhyJi++gDOBY4tqXu/fvvRq0m+Pfx47rc2WHbE3U5tL8OZnxo1Odc7wBuzcyVLfUnR8ShYz/zZ8ZaqXcBcASwrk9DnRH6sD9G218L/BLVDZc76roF9X5RoYnsj8y8JDNfk5lHtT6Ax4DHWuq+0exvMX304e/jvro8ssO20bAY2V8HMzo0IuI1VIl9L9X1/a2+SXWk8JL2EfGqtnZfoJpEOr+t/j11+fn+jHb669P+ICL+kOr9/6XM3N6yqfUrgjWOfu0P9Uef9sdXgReA89razgXeBDy8v1NTMIPnNOo38zaqG2puAX6/7frxeW3tF1EtUfFURPzrzPwXgMx8OCJWAb8TEbdn5v+NiDcAVwFfzsyXnDtUZ/3aHxFxOfV7D1zU1sfZeORXpF/7Q/3Rx8+rzRHxR8BVEXE78DWqS3Y/C7wSuHTcwQzqLsgmH8AHqdZWGaG6smlv/fMDLW3eXW/b32NRS/sjqe5kvR84tMNrXgo8RHVecDPwe8ArBv1eTIXHZO4PqruO99fHxwf9fgz6MYi/j7rN4/Xr/LR+bAceH/T7MejHgD6vLgQ2UN0UuJPqSGVpyXinxNLokqQDw4ye05AkdcfQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JU7P8BGjIRQddamkQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Our measument may have no analytic way to propogate uncertainties, \n",
"# so now lets add Gaussian errors to the distances!\n",
"\n",
"N_sim = 1000 # how many simulated measurments to make\n",
"\n",
"dist_1_sim = np.random.normal(loc=dist_1, scale=err_1, size=N_sim)\n",
"dist_2_sim = np.random.normal(loc=dist_2, scale=err_2, size=N_sim)\n",
"\n",
"# you get a gaussian distribution of distances, centered at the true value, w/ the error as the stddev\n",
"_ = plt.hist(dist_1_sim, histtype='step', bins=25)\n",
"plt.plot([dist_1, dist_1], [0,100], c='C1')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "97df9c33",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"543.3492376915606 543.3496578861784 0.5127766990013339\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEACAYAAABCl1qQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVoElEQVR4nO3dfbBkdX3n8feHx4wsMIQZEEZhYoEGqkCjMwalJAMENQE0RA2ERYugzq4J7CYh6qVWZCoPZFZZkjKJy06iRhSiYDBCUCLWhDUSsBw28aEgoAgGBmeYQQchDMrDd//oHupy6Tunu6fP7Xvnvl9VXWf6nN85/e1f9fTnnvM753SqCkmStmeXcRcgSZr9DAtJUiPDQpLUyLCQJDUyLCRJjXYbdwGjsGjRolq6dOm4y5CkOeW2227bXFWL+2m7U4TF0qVLWbdu3bjLkKQ5Jcn3+m3rYShJUiPDQpLUyLCQJDUyLCRJjQwLSVIjw0KS1GigsEhyVpItSf56muXHJPl4kvuSPJRkU5JrkvzcNO3vTbKhx+P+Id6LJKklfV1nkWQRcBmwHNh3mjavBG4BrgOWV9WGJIcCfwPckuTEqrp56npV9fxhi5ckzYx+9ywuB+4GXtuwrR8Db62qDQBV9T3gbGBP4APDlylJGqd+r+BeWVX3J1m6nTb3A79XVQ9PnllVdyX5AZ29EknAsavXsn7L1r7bL1m4gJsnTmixImn7+gqLqmocQ+i2+fNpFu8O/HCAuqSd2votW7l39cl9t186cX2L1UjNWj8bKsnPAnsD10yz/OIk30qyMckdSS7tjpFIkmaJmTh19reAh4GLeywr4HHgWOAF3bZvAdYl2e7Ad5KVSdYlWbdp06YRlyxJmqzVsEjyauC/0hnzuK9Hk+VV9ftV9XBVPVFVa4HfBA4F/nB7266qNVW1rKqWLV7c1x12JUlDai0skrwI+CxwYVVd1atNVW3uMfvzwJPAKW3VJkkaTCthkeRg4EbgY1W1epB1q+op4CHggDZqkyQNbuRhkWQx8CXg+qqamDT/qCR7THq+IslJPdbfFdgf6LXXIUkag5GGRZL96OxR/DPw36csvg44eNLzFcB5PTbzOjqn9N4wytokScMb2c+qJvlPwBeApcC1wEVJJjdZ2GO1U5OcC6wBngCOAT4MbATeN6raJEk7pt97Q50JXArs2p11epLXAw9W1dHdeb8I/Hz33xf2sdm/oHNK7RnABcBewI/oBM7vV9X6vt6BJKl1/V7BfSVwZUObvwOyvTZT2m8C/qT7kCTNYv6ehSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhr19Rvckrbv2NVrWb9la9/tlyxc0GI10ugZFtIIrN+ylXtXnzzuMqTWeBhKktTIsJAkNTIsJEmNDAtJUiMHuKWd1DBnaN08cUKLFWkuGygskpwF/Dnwd1V19jRt9gf+J3AysCtwJ3BhVd00Tfs3AxcALwR+DHwaeH9VPTZIbZKebdAztJZOXN9iNZrr+joMlWRRks8AfwTsu512ewNfBo4AjgKeD3we+FKSk3q0Pwe4Cri0qg4AjgPeCPx9kl0HfC+SpJb0O2ZxOXA38NqGdu8GjgTeWVWbq+rpqvpj4OvAZUme2ZNJsh9wKfCZqroCoKruAc4HjgfeNtA7kSS1pt+wWFlV76VzmKinJAHeDtxZVbdPWXwN8CI6IbDNr9HZS7lmStsvAFuBd/RZmySpZX2FRVXd30ezw4CDgW/0WPb17vQXJs07rjt9VvuqegK4HTgmyZ791CdJatcoT519cXf6/R7LHuhODx+g/S509kakee/9u10OX5gYdxmax0Z56uy2ge9eZzFtm7dwB9o/S5KVwEqAQw45pN8apTnpZbvfx623fI8z/m//Zyx5s0KNUhvXWVTL7TsrVa0B1gAsW7ZsqG1Ic8XLD9kPgHt/w5sVajxGeRjq4e50rx7LnjelzTDtJUljMsqwuKs7PajHsoO7028P0P5p4LujKU2StCNGGRbfoTMwfXSPZdvm3TRp3penLAMgye50Luq7taoeH2F9kqQhjSwsqqqAjwIvSXLklMVvorOX8I+T5l0N/Ag4bUrbX6JzGOojo6pNkrRjRn3X2Q8AdwBrurcI2SXJBcBLgXdV1ZPbGlbVD4DfBd6c5D8DJFkKXEInVD4+4tokSUPq995QZybZAHytO+v0JBuSTL2g7hE6F9vdCXwT2EDnhoInVdUXp263qj4CnAGcn+RB4CvAdcApVfXUkO9JkjRifZ06W1VXAlf22XYzndt+9KWqrqZzSEqSNEv540eSpEaGhSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlq1MYtyiXNQUsWLmDpxGC/l3HzxAktVqTZxLCQBDDwF/8gwaK5z8NQkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhqNPCySrEryaJINPR4/SlJJDuy2PXs7bX971LVJkobT1i/lXVJVq6bOTHIFcGhVbWxqK0maPdoIi+/0mplkH+A04NwWXlOS1KKRh0VVfXKaRacDTwFXjfo1JUntmskB7rOBq6vq0Rl8TUnSCLQ1ZvEsSQ4HXg28t8fiZUluAI4A9gS+CfxZVV07E7VJvRy7ei3rt2ztu/2ShQtarEYavxkJCzp7FXdV1Vd6LPsZ4Jyq+mr3LKkLgc8luaCqVs9QfdKzrN+ylXtXnzzuMqRZo/XDUEl2Ad4KfKzH4quA5VX1VYCq2lhV5wK3AX+QZOl2trsyybok6zZt2tRC5ZKkbWZizOJE4GDg8qkLquqxqnqsxzrX0dnref10G62qNVW1rKqWLV68eGTFSpKeaybC4mzghqp6YIB1tl2HccDoy5EkDarVMYtJ11acNc3yVcAfVdUTUxYd2J1ubq86zScOWEs7pu0B7tOBR+kcVurlIuBa4P9Nmf/LwNPAje2VpvnEAWtpx7R9GOps4JM99hwmuyzJS6CzJ5LkEuCVwAeq6tst1ydJ6kNrexZJDqNzbcV/2U6zE+mcKXVtkoXAAjrXWZxVVVe0VZskaTCthUVVfQdIQ5u1wNq2apAkjYa/ZyFJamRYSJIaGRaSpEaGhSSpkWEhSWo0U3edlbSTWbJwAUsnrh+o/c0TJ7RYkdpkWEgayqBf/IMEi2YfD0NJkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpUSthkeTeJBt6PO7v0Xb/JH+V5PtJHkzyT0lWtFGXJGk4u7W14ap6flObJHsDXwa2AEcBPwDeC3wpyS9V1Y1t1SdJ6t+4D0O9GzgSeGdVba6qp6vqj4GvA5claS3MJEn9G1tYJAnwduDOqrp9yuJrgBcBx894YZKk5xjnnsVhwMHAN3os+3p3+gszV44kaTqthUWSi5N8K8nGJHckuTTJoklNXtydfr/H6g90p4e3VZ8kqX9tjQkU8DhwLPAY8Brg48CvJjmmqjYA+3bbPtZj/W3zFk73AklWAisBDjnkkNFUrTnh2NVrWb9l60DrLFm4oKVqpPmhrbBYXlWbJz1fm+Q3gWuBPwTeMWlZDfMCVbUGWAOwbNmyobahuWn9lq3cu/rkcZchzSuthMWUoNjm88CTwCnd5w93p3v1aPu8KW0kzXFLFi5g6cT1A7W/eeKEFivSIGbs1NSqeirJQ8AB3Vl3dacH9Wh+cHf67dYLkzQjBv3iHyRY1L6RD3AnWZHkpB7zdwX2B7btdXyHzkD20T02s23eTaOuT5I0uDbOhloBnNdj/uvo7MncAFBVBXwUeEmSI6e0fRPwXeAfW6hPkjSgtk6dPTXJuUn2SMergA8DG4H3TWr3AeAOYE2SRUl2SXIB8FLgXVX1ZEv1SZIG0EZY/AXwu8AZwD3AD4FPA/8AvKKq/n1bw6p6BDgOuBP4JrABOBk4qaq+2EJtkqQhjHyAu6o2AX/SffTTfjOd235I0jM8e2p28UZ9kmYlz56aXcZ911lJ0hxgWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaGRaSpEaGhSSpkWEhSWo08rBIsk+S85LcmuShJA8n+VaS9yTZfUrbs5M8mmRDj8dvj7o2SdJwdmthm58CjgfOBD4H7Aq8DVgDvAY4dUr7S6pqVQt1SJpHlixcwNKJ6wde5+aJE1qqaOfSRljsAvxpVX22+/xp4CNJTgR+PclJVXVjC68raR4b5kt/0HCZz9oYs7gS+ESP+bd0p8tbeE1JUotGvmdRVZdPs2iP7vSHo35NzW3Hrl7L+i1b+26/ZOGCFquR1Esbh6Gmsxx4ErhuyvxlSW4AjgD2BL4J/FlVXTuDtWmM1m/Zyr2rTx53GZK2Y0ZOnU3yQuANwIeq6v4pi38GuKiqDgVeCtwJfC7JRMM2VyZZl2Tdpk2bWqlbktTRelgkCXAZcDvwP6YsvgpYXlVfBaiqjVV1LnAb8AdJlk633apaU1XLqmrZ4sWL2ylekgTMzJ7FB4EjgVOq6vHJC6rqsap6rMc619E5RPb6GahPktSg1TGL7qGkXweOq6oNA6y6sTs9YPRVSZIG1dqeRZLzgN8BfrGq7u7O23/yoaUkq6Ze1d11YHe6ua36JEn9ayUskpwDXAS8tqrumLToVGDVpOcXAUf12MQv07mYz4v3JGkWGPlhqCRnAH8JXA+cluS0SYtfBmyZssplSd5aVXcm2Qd4P/BKYHVVfXvU9UmSBtfGmMUEnT2WU3nufaAAPj7p3ycCbwWuTbIQWEDnOouzquqKFmqTJA2hjSu4XzZA27XA2lHXIEkaLX/PQpLUaCZv96F5wns9STsfw0Ij572epJ2Ph6EkSY0MC0lSI8NCktTIsJAkNTIsJEmNDAtJUiNPnVUjr5uQZFiokddNSPIwlCSpkWEhSWpkWEiSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIaeW+oecgbA0rDGeb/zs0TJ7RY0cwxLOYhbwwoDWfQ/ztLJ65vsZqZZVhImreWLFww0Bf6fN7LNiwkzVs7yyGimTArBriTvDnJbUkeTHJfkkuSPG/cdUmSOsYeFknOAa4CLq2qA4DjgDcCf59k17EWJ0kCxhwWSfYDLgU+U1VXAFTVPcD5wPHA28ZYniSpa9xjFr8G7AtcM2X+F4CtwDuAj810UXONp8JKs9MwA+izdRxl3GFxXHf6jckzq+qJJLcDxyTZs6p+3FYBs+286UHrgU5NngorzT6Dflccu3rtrA2XcYfFi7vT7/dY9gDwCuBFwB1tFTDbzpv2Gghp/hr0i38mr+NIVc3Yiz3nxZO7gMOBn5q695DkU8DpwKur6pYe664EVnafvgS4s+VydxaLgM3jLmInYV+Ohv04GsP046FVtbifhuPes9hm4MSqqjXAmhZq2aklWVdVy8Zdx87AvhwN+3E02u7HcZ86+3B3ulePZc+b0kaSNCbjDou7utODeiw7GHga+O7MlSNJ6mXcYfHl7vToyTOT7A4cAdxaVY/PeFU7Nw/djY59ORr242i02o/jHuD+aeAe4IaqOn3S/DcAnwPeXlUfHVd9kqSOsYYFQJK300nEt1XVFUmWAl8E7gdOqqqnxlmfJGkWhAVAkrcAFwAvAH4CfBq4sKoeG2thkiRg/GMWAFTV1VX18qo6oKpeUFXnGxTTS3Jvkg09Hvc3rHdpkkqyqseyfZKcl+TWJA8leTjJt5K8pzuGtFNqqS/3TrIyyXVJ7k6yMck9ST6R5PDW3swYtdGPPdou6X4ux/8Xbova6stht7vNbLnOQgOqqucP0j7JMuC/bafJp+jcvPFMOuNFu9K5keMa4DXAqcNVOvu10JevAP4P8L+BM6vqkSQvBj4DfC3Jy6tqpzvLr4V+nOrDwD4DFTVHtdWXg253slmxZ6F2JdkN+Evgb7fTbBfgT6vqs1X1dFU9UVUfoXNI8JQkJ81ErbNdn30JnVvYnFdVjwBU1V3Ae+jcOPOcVoucAwbox23t3wIcBXytzbrmokH7cliGxfzwe8BjwGXbaXMl8Ike87fdamX5qIuao/rpy38BXtvj5Iz7utN92yhsjumnHwFIshD4EPCu7jp6tr77ckd4GGonl+QwYAJ4NXDAdO2q6vJpFu3Rnf5wxKXNOQP05cP0vvPAy7vTfxp9dXNHv/04ySXAl6rqH5Jc0Gpxc8wQfTk09yzmqCQXdwegNya5ozu4tahH0zV0Di/dPuRLLQeeBK4buthZru2+TLJXkjcCHwQ+Clw9grJnnTb6MckKOr+c+TujrXZ2a+szOcB2n8OwmJsKeBw4ls7pxr8FvAVYl+SZAazuNSwHARcP8yJJXgi8AfhQVfV1xsQc1GpfJvkknb2Mv6XzQ17n1Ww4X330Rt6PSX6Kzpfh+VU1n+5K29Znsq/tTr92lY859gAW9Zh3avfD8Ffd5wcCDwGvmdRmRbfNqj5eI8D1wDo6t5Af+/uew325B/Aq4F/p3Er/sHG/77nQj3S+BG+cMu+mztfW+N/zXOrLfre7vYdjFnNQ9f4r6/N0Dhed0n3+IeDqqhr2+PgHgSOBV9VOfH+umejLqvoJcEuSNwH/RufMleOH2dZsNep+THI0nb98f25kRc4RbX0m+9zutGbFFdwajSQbgAOqapckj9DZ5Zx8Rs4ewH7AfwCPQu/zrpNMAOcBx1XV3a0XPguNqi97bPff6PxC5N5V9R8jL3yWGbYfuwPZE8DU3xj+aWB3YGP3+SVVdUmLb2HWaPEz+cx2t9tw3LtcPgbeRV1B555ZU+fvCjwBPNiw7nYPndAJiY3AEZPm7Q8sHfd7nyt9Cfwq8PPTrPcv3fUOHvf7n+39OE37m9iJD0O1+JkcervbHg5wzz0r6HyhT/U6OqdC3zDshpOcA1xE5xqByb97fiqwatjtzmIraKcv30Dn6vdnSXIg8LPAhu5jZ7GClj6T89AK2unLHd6uYTE3nZrk3CR7pONVdG6FsBF43zAbTHIGnWPp/wyclmTVtgfwKyOqezYaeV92vTPJbyTZA545H/4qYE/g3VX19A5XPru01Y/zUVt9uWPbHfdul4+Bd1MX0znn/CvAemAL8O907kW0ZJp13kTnL9kf0NlNfbT7/H9NavOv3WXTPf563O99DvXlId3/fF8FHui23QhcCxw/7vc9V/pxSvtbust/0m2/be/soHG//7nQl8Nsd+rDAW5JUiMPQ0mSGhkWkqRGhoUkqZFhIUlqZFhIkhoZFpKkRoaFJKmRYSFJamRYSJIa/X+1qx1sR+6yyQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# now lets do do the final simulation, which prodcues 1000 estimates of the measurment\n",
"# in this case you can just do the math:\n",
"meas_sim = np.sqrt(dist_1_sim**2 + dist_2_sim**2)\n",
"\n",
"# but sometimes you need to loop over every simulation and collect the results into an array\n",
"\n",
"# lets plot the simulated results\n",
"_ = plt.hist(meas_sim, histtype='step', bins=25)\n",
"plt.plot([meas, meas], [0,100])\n",
"\n",
"meas_m = np.mean(meas_sim)\n",
"meas_err = np.std(meas_sim)\n",
"print(meas, meas_m, meas_err)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "d62a4cec",
"metadata": {},
"outputs": [],
"source": [
"# so our simulated resulting error for the thing is around 543.4 +/- 0.5"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8bc9d49a",
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment