{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Error Propagation using Random Gaussian Numbers\n", "\n", "Example calculation of propagating uncertainties, both when adding and multiplying number, and also in the general case. The propagation can be done both analytically (using the error propagation formula) and also using simulation.\n", "\n", "The example is based on FIRST doing the error propagation analytically, and then verifying it by running a so-called Monte-Carlo (MC) program, which uses random numbers for propagating errors.\n", "\n", "## References:\n", "- Barlow: page 48-61\n", "- Bevington: page 36-48\n", "\n", "## Author(s), contact(s), and dates:\n", "- Author: Troels C. Petersen (NBI)\n", "- Email: petersen@nbi.dk\n", "- Date: 12th of November 2020" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "DO THE FOLLOWING ANALYTICAL EXERCISE FIRST!!!\n", "\n", "1. A class of students estimate by eye, that the length of the table in Auditorium A is $L = (3.5\\pm 0.4)$m, and that the width is $W = (0.8\\pm 0.2)$m.\n", "\n", " Assuming that there is no correlation between these two measurements, calculate ANALYTICALLY what the Perimeter (P), area (A), and diagonal (D) length is including (propagated) uncertainties. Repeat the calculation, given that the correlation between length and width is $\\rho(L,W) = 0.5$ - not an unreasonable number, given that they are estimated by the same (uncertain) scale.\n", " \n", "NOTE: This is a complete standard problem, that you will be asked to solve again and again in the course. For this reason, make sure that you understand how to do it, and become good at doing it reasonably fast." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Defining the parameters:\n", "mu1 = 3.5\n", "sig1 = 0.4\n", "mu2 = 0.8\n", "sig2 = 0.2\n", "rho12 = 0.5 # Correlation parameter!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "if not (-1.0 <= rho12 <= 1.0): \n", " raise ValueError(f\"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Note on analytic solutions:\n", "\n", "Python includes symbolic algebra in the package *SymPy*, which seems to be both simple and powerful (and in Python). In addition, printing with Latex can also be included (see below), which (in combination) is very nice.\n", "\n", "Below is a SymPy and Latex example with the hope that it will wet your appetite." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from IPython.core.display import Latex\n", "\n", "def lprint(*args,**kwargs):\n", " \"\"\"Pretty print arguments as LaTeX using IPython display system \n", " \n", " Parameters\n", " ----------\n", " args : tuple \n", " What to print (in LaTeX math mode)\n", " kwargs : dict \n", " optional keywords to pass to `display` \n", " \"\"\"\n", " display(Latex('$$'+' '.join(args)+'$$'),**kwargs)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$P = 2 L + 2 W$$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$$\\sigma_{P} = \\sqrt{4 \\sigma_{L}^{2} + 4 \\sigma_{W}^{2}}$$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$$P = (8.6 \\pm 0.9)\\,\\mathrm{m}$$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Import SymPy: \n", "from sympy import * \n", " \n", "# Define variables:\n", "L,W,P,A,D = symbols(\"L, W, P, A, D\")\n", "dL,dW,dP,dA,dD = symbols(\"sigma_L, sigma_W, sigma_P, sigma_A, sigma_D\")\n", "\n", "# Perimeter:\n", "# Define relation, and print:\n", "P = 2*L + 2*W\n", "lprint(latex(Eq(symbols('P'),P)))\n", "\n", "# Calculate uncertainty and print:\n", "dP = sqrt((P.diff(L) * dL)**2 + (P.diff(W) * dW)**2)\n", "lprint(latex(Eq(symbols('sigma_P'), dP)))\n", "\n", "# Turn expression into numerical functions \n", "fP = lambdify((L,W),P)\n", "fdP = lambdify((L,dL,W,dW),dP)\n", "\n", "# Define values and their errors\n", "vL, vdL = mu1,sig1\n", "vW, vdW = mu2,sig2\n", "\n", "# Numerically evaluate expressions and print \n", "vP = fP(vL,vW)\n", "vdP = fdP(vL,vdL,vW,vdW)\n", "lprint(fr'P = ({vP:.1f} \\pm {vdP:.1f})\\,\\mathrm{{m}}')\n", "\n", "\n", "\n", "\n", "# NOTE: Do the above analytical calculation before you continue below! Possibly use SymPy for the differentiations." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import numpy as np # Matlab like syntax for linear algebra and functions\n", "import matplotlib.pyplot as plt # Plots and figures like you know them from Matlab" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Error propagation - Simulation\n", "\n", "Now we want to try to see, if we can solve the above error propagation problem using simulation. The method is relatively straight forward: You simply take \"realistic\" values of the input parameters x (here Length (x1) and Width (x2)), calculate the resulting value y (here Perimeter, Area, and Diagonal), and do this many times. The resulting distribution of y should be centered around the value y(x1,x2), and the standard deviation should reflect the uncertainty in y from the uncertainties in the input variables.\n", "\n", "This is a much more clumsy of calculating the uncertainty, but comes with the advantage, that if the resulting uncertainty is not Gaussian, then one can actually see this (i.e. it avoids the assumptions used in the usual error propagation formula)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# First we set the parameters of the program:\n", "N_exp = 10000 # Number of \"experiments\" (i.e. drawing from random distributions)\n", "save_plots = False\n", "r = np.random\n", "r.seed(42)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define parameters for two random numbers (Gaussianly distributed):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "mu1 = 3.5\n", "sig1 = 0.4\n", "mu2 = 0.8\n", "sig2 = 0.2\n", "rho12 = 0.5 # Correlation parameter!" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "if not (-1.0 <= rho12 <= 1.0): \n", " raise ValueError(f\"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we calculate numbers that allows the transform from uncorrelated variables `u` and `v` to correlated random numbers `x1` and `x2` below (see Barlow page 42-44 for method):" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# NOTE: The following part of the code is for generating two Gaussian random numbers that are correlated.\n", "\n", "# Define the parameters needed for the transformation:\n", "theta = 0.5 * np.arctan( 2.0 * rho12 * sig1 * sig2 / ( np.square(sig1) - np.square(sig2) ) )\n", "sigu = np.sqrt( np.abs( (((sig1*np.cos(theta))**2) - (sig2*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n", "sigv = np.sqrt( np.abs( (((sig2*np.cos(theta))**2) - (sig1*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n", "\n", "# Produce random numbers with the (possible) correlation:\n", "u = r.normal(0.0, sigu, N_exp)\n", "v = r.normal(0.0, sigv, N_exp)\n", "x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v\n", "x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v\n", "x12_all = np.array([x1_all, x2_all])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Now YOU have to define what you want to consider (Perimeter, Area, Diagonal):\n", "y_all = x1_all - 2*x2_all # Silly formula - you have to put this in yourself!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the above is nothing more than a matrix multiplication written out! Also note that the absolute value is taken before the square root to avoid `np.sqrt(x)` with `x<0`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "Plot both input distribution and resulting 2D-histogram on screen:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4kAAAIzCAYAAACtCFzyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZylV10n/s+3tySdzt5NEkjCImERhQA9DOBrNBp2QXAEBASDg0ZRVBz4jcIAOsogzLjg/FAgCBKQ1bgADqtRRGXtYFgDEpCQJvvS2bqT7uo688e9Xal0P09VV93qutW33++87qvqnmc55z73Vrq+9X3O91RrLQAAAJAkq8Y9AAAAAFYOQSIAAAAzBIkAAADMECQCAAAwQ5AIAADADEEiAAAAMwSJwCGrqr5SVWeOexzjVFU/XlWXVdUtVfXgju2tqu49hnGdWVVbl+hcm6rq61V1+H7s+wdV9QtL0e/+qKqPV9XPjniON1TVy+fYPud7WFXfrqpHjTKGWed6V1U9ZT/2e2BVfXIp+gTgDlV136q6aNbjpqp6YVUdX1Ufq6pvDL8eN9d5BInAROr6xbeqnltV/7zneWvtAa21j89znnsMf8lec4CGOm6/l+QFrbUNrbV/HdcgDnAw+htJ/qy1dtt+7Pu/k/z3qlp3gMay5Fprv9Ba+5392beq3lpVrzwQ46iqByZ5UJL3zbdva+2LSbZV1ZMOxFgADlWtta+31s5orZ2R5KFJtif56wz+LbygtXZ6kguGz3sJEgHGaAUEn3dP8pUxj+GAqarDkpyd5M/3Z//W2hVJvpbkxw7kuCbUzyd5R2ut7ef+7xgeA8CBcVaSb7bWLk3y5CTnDdvPSzLnXR+CROCQNTvbWFUPq6otw9syrqqqPxju9onh123DWzIfUVWrquplVXVpVV1dVW+rqmNmnfenh9uuq6qX79XPb1XV+VX151V1U5LnDvv+VFVtq6orqup1szNZwyzbLw5vEbm5qn6nqr5neMxNVfXevsxX31ir6rCquiXJ6iRfqKpv7sf1Oqyqfq+qvjO8Rm+oqiOG286sqq1V9aJhP1dU1c/MOvaEqvrAcLyfq6pX7snqVtWea/yF4TX+yVnH9Z3vCVX11eH1+G5Vvbhn2P8xybbW2tbhcccPx/mk4fMNVXVJVf30rGM+nuRH57gOf1FVV1bVjVX1iap6wKxtb62qP66q/zsc22eq6ntmbX90VX1teOzrklRPH4dX1Y6q2jh8/rKqmqqqo4fPX1lVr53V5ytnHfv/Da/X5VX1X2a1n5Pkp5L8t+F1/sCsLs+oqi8Ox/WeGt6aW1Ubq+pvh5/N66vqn6qq73eHxyf5x+Fxhw33//5Z/d9l+Jo2zbrOZ9UgkAdg6T0jybuG3584/EPonj+I3mWuA8f9F2yAleKPkvxRa+3tVbUhyfcN238wyb8nOba1NpUkw1+8n5vkh5NcneRtSV6X5DlV9b1J/iTJ45J8Nsmrktxtr76enORpSX46yWFJvjfJryXZkuSUJB9K8otJXjvrmMdlcNvIqUk+n+SRGfzCf12STyV5Zu74C+Fsz+0aa2vtOUk2VFVL8qDW2iX7cY1ek+ReSc5IsivJO5O8IslLhttPSnLM8PU+Osn5VfU3rbUbkvxxkluH+9wjyUeSXJokrbUf3HscNZgrOtf53pzk6a21f6rBvIp79oz5+5N8fc+T1tr1w/fvbTW4PfJ/Jrmotfa2WcdcnOQn5rgOH0ryX5LsHF6TdwyvyR7PzOD9+nwG78n/TPKMYcD3l8Nj35fkBUl+Icnb9+6gtXZbVX0uyQ8Nj/nB4fX6gWH/P5jkD/c+rqoel+TFGfz1+N+TvGnWOc+tqkcm2dpae9lehz59OObbkvxLBp+ZNyR5UZKtSfYEdg9Psk+msKqOzOA9+Pqwr9ur6t1Jnp3k12ddl79rrV0z3Oe7VbUryX2TfHHvcwIcLB77w0e2667fvSx9XfjF27+Swf+r9zi3tXbu3vsN/3j8Y7nj3+gFESQCk+xvqmpq1vN1Gfzi3mVXkntX1cbW2rVJPj3HeX8qyR+01r6VJFX1kiRfHma6nprkA621PVmyVyT5lb2O/1Rr7W+G3+9IcuGsbd+uqjdmEBzMDhJf01q7KclXqurLST46q/8PJXlwuoPE3rHuCXr3R1VVkp9L8sDW2vXDtldlECju+QdoV5LfHp73g8NM5X2Hwc5PJPm+1tr2JF+tqvOSnDlPt53ny+C92ZXke6vqC8Og8Yaecxyb5ObZDa21j1bVX2QwJ+OEDALJ2W4eHteptfaWPd9X1W8luaGqjmmt3Ths/qvW2meH29+RZE9W+glJvtpaO3+47bUZBGF9/jHJD1XV+5I8MMnvDp//Q5L/kOSfOo55egbzL788a3zPnKOPPf5Pa+3y4TEfyB1B764kJye5+zCA7+ozueN6zb7W52UQ2L+ktTad5DlJ/tdex815rQEOBtddvzuf/chpy9LX6pO/cVtrbfN+7Pr4JJ9vrV01fH5VVZ3cWruiqk7O4A/HvdxuCkyyp7TWjt3zyCA71+d5Se6T5GvD2yGfOMe+d80wCzZ0aQZ/dDtxuO2yPRuGQdF1ex1/2ewnVXWf4S19V9bgFtRXJdm41zFXzfp+R8fzDYsY60JsSrI+yYXDWw+3Jflw7sgwJcl1ewWe24fj2jTsc/brvtM16NF3vmQQdD4hyaVV9Y9V9Yiec9yQ5KiO9nMzyBb/WWtt7/fnqCTbuk5WVaur6tVV9c3he/Xt4abZ79eVPWPe+7PRMvd1+McMAumHJPlSko9l8MeDhye5ZPjHjL3dqY/c+b2fS9+Y/3eSS5J8tKq+VVV9hQ72XK+Za91a+0wG2eMfqqr7Jbl3kvfvdVzvtQZgJM/MHbeaJoP//549/P7szFNkTJAIkKS19o3W2jMzuEf/NRlkQI5Mx611SS7PoODLHqclmcogcLsig1tGkyQ1mLN3wt7d7fX89RkUSzm9tXZ0kpemZ67aIsw11oW4NoNg9AGzAu9jWmt9wels1wz7PGVW26kL7P9OWmufa609OYP362+SvLdn1y9mEPzPqKrVSd6Ywa23z699q6reP8kXes73rAxuF35UBrfC3mPPafdj2Fdk1useZmfnug6fzCBz+uNJ/rG19tUM3r8fzXDu33x9DPefbX+Lygx2bu3m1tqLWmv3SvKkJP+1qs7q2O/WJN/MXtc6g2ziszPIIp4/u8JsVd01g+z+1wNwEGtJppfpv/1RVeszmKbxV7OaX53k0VX1jeG2V891DkEiQJKqenZVbRreFrcns7E7gwBnOoO5eHu8K8mvVdU9h/MXX5XkPcOs1/lJnlRVjxzOB/gfmT+AOCrJTUluGWZcnr9kL2zuse634XV5U5I/rKq7JElV3a2qHrsfx+7O4B+q36qq9cPX+NN77XZV7nyNe1XVuqr6qeEtnrsyuHZ9k0E+m+TYqpo9L/Slw6//JYMlQN42DBz3+KEM5v11OSrJ7Rlkh9dncD331/9N8oCq+s81qGr7KxnMu+w0zEJfmOSXckdQ+MkMKoL2BYnvzaAY0vcOf0n4zb227/d1TpKqemJV3XsY0O65zn3X+oMZXLvZ3p5BkPvsDILy2c5M8vettdv3dzwAzK+1tr21dsKsaRBprV3XWjurtXb68Ov1c51DkAgw8LgM5vvdkkERm2e01m4b/qL+P5P8y/A2y4cneUsGv/x+IoPiILcl+eUkaa19Zfj9uzPI6tycwX3/c/0i/OIMMlQ3ZxCIvWcJX1fvWBfh1zO49fDTw1st/y6DTNf+eEEGmbcrh+N5V+58TX4ryXnDa/z0/TjfczKYv3lTBsVfnt21U2ttZ5K37tleVQ9N8l+T/PQweH1NBn8E/o3h9pMzKCT0N13nyyDQuTTJd5N8NXPPXd17LNdmULDo1RkEmadnUCRmLv+YZG0Gwe6e50fljqq7e/fxoQzmsv59Bu/V3++1y5szmMu5rar6XuNsp2fwPt+SQYGkP5ljbdFzk/zUMKDcM56tGcwDbtl3PuNPZVAcB+Ag17K7TS/LY7nU/i9nBMBCDbN32zK4lfTfxz2elaKqXpPkpNba2fPuPHpfmzIIUB7cWtsxz76/n8GaUn9yoMc1iarqnUneO6swU6rqLUkun11RtQZLY5zbWuubSwpw0Hjogw5rn/zw3oXMD4zD7/rvF+5n4ZqRqG4KsMRqsAbfBRncZvp7GRQd+fY4xzRuw1tM12VwLf5DBoWCfnY5+h4uuXC//dx3rmqjzKO19qzZz6vqHkn+cwbVd2fv96UkAkRgIgzmJE5W4m1st5vWYKHgz1bVF6rqK1X1Pzr2OawGi/peUoMFie+x/CMFWLAnZ1Aw5vIMbtd7RnPbxlEZzEu8NYN5c7+feSqrcXCrqt9J8uUk/1sWHeDgMrbbTYdzFo5srd1SVWuT/HOSX22tfXrWPr+YwZpcv1BVz0jy4621nxzLgAEAAPbykAcd1v7pw711yJbUhrt+Z1luNx1bJrEN3DJ8unb42DtifXLuWBz6/CRnzZ4QDwAAwNIa65zEYcnxCzNYYPePhwvvzna3DBcFbq1NVdWNGaw31rWAMAAAwLJqadk9YbNKxhokDsuPn1FVxyb566r6vtbal2ft0pU13OcdqKpzkpyTJEceeeRD73e//apPAAAArAAXXnjhta21TeMeBwMrorppa21bVX08g3XKZgeJW5OcmmTrcOHhY5Lss/Bja+3cDNZnyubNm9uWLVsO+JgBAIClUVWXjnsMo1DddIlU1aZhBjFVdUSSRyX52l67vT/JnjW0nprk71UIBAAAOHDGmUk8Ocl5w3mJqzJYfPdvq+q3k2xprb0/yZuTvL2qLskgg/iM8Q0XAADgzlqS3ROWSRxbkNha+2L2Wlx32P6KWd/fluRpyzkuAACAQ9mKmJMIAABwsDInEQAAgIklkwgAALBILZm4dRJlEgEAAJghkwgAADCC6XEPYInJJAIAADBDJhEAAGCRWtrErZMokwgAAMAMmUQAAIDFasnuyUokyiQCAABwB5lEAACARWpR3RQAAIAJJpMIAACwaJXdqXEPYknJJAIAADBDJhEAAGCRWpJp1U0BAACYVDKJAAAAIzAnEQAAgIklkwgAALBILTKJAAAATDCZRAAAgBFMN5lEAAAAJpQgEQAAgBluNwUAAFgkhWsAAACYaDKJAAAAi9RS2T1hubfJejUAAACMRCYRAABgBJbAAAAAYGLJJAIAACyS6qYAAABMNJlEAACARavsbpOVe5usVwMAAMBIZBIBAAAWqSWZnrDc22S9GgAAAEYikwgAADAC1U0BAACYWDKJAAAAi9Sa6qYAAABMMJlEAACAEUybkwgAAMCkkkkEAABYpJZk94Tl3ibr1QAAADASmUQAAIBFU90UAACACSaTCAAAsEgtyfSE5d4m69UAAAAwEplEAACAEexu1kkEAABgQskkAgAALFJLWScRAACAySWTCAAAMIJp6yQCAAAwqWQSAQAAFqkl5iQulao6tar+oaourqqvVNWvduxzZlXdWFUXDR+vGMdYAQAADhXjzCROJXlRa+3zVXVUkgur6mOtta/utd8/tdaeOIbxAQAAzKmlrJO4VFprV7TWPj/8/uYkFye527jGAwAAcDCrqmOr6vyq+trwjs1HVNXxVfWxqvrG8Otx851nRdw8W1X3SPLgJJ/p2PyIqvpCVX2oqh6wrAMDAACYx3RWLctjP/xRkg+31u6X5EEZJOJ+I8kFrbXTk1wwfD6nsQeJVbUhyV8meWFr7aa9Nn8+yd1baw9K8v8n+Zuec5xTVVuqass111xzYAcMAACwwlTV0Ul+MMmbk6S1trO1ti3Jk5OcN9ztvCRPme9cYw0Sq2ptBgHiO1prf7X39tbaTa21W4bffzDJ2qra2LHfua21za21zZs2bTrg4wYAAEiS1pLdbdWyPJJs3JMcGz7OmTWUeyW5JsmfVdW/VtWfVtWRSU5srV0xGGu7Isld5ntNYytcU1WVQZR7cWvtD3r2OSnJVa21VlUPyyCovW4ZhwkAALBSXNta29yzbU2ShyT55dbaZ6rqj7Ift5b2nWhcfiDJc5J8qaouGra9NMlpSdJae0OSpyZ5flVNJdmR5BmttTaOwQIAAKxgW5Nsba3tqfNyfgZB4lVVdXJr7YqqOjnJ1fOdaGxBYmvtn5PMWSu2tfa6JK9bnhEBAAAsVGV67rBmWbTWrqyqy6rqvq21ryc5K8lXh4+zk7x6+PV9851rnJlEAAAAls4vJ3lHVa1L8q0kP5PBlL33VtXzknwnydPmO4kgEQAAYJFasqeozNi11i5K0jVn8ayFnGdlvBoAAABWBJlEAACAEeyesNzbZL0aAAAARiKTCAAAsEgtlek2/uqmS0kmEQAAgBkyiQAAACMwJxEAAICJJZMIAACwSC3J9ApZJ3GpTNarAQAAYCQyiQAAAItW2R3VTQEAAJhQMokAAACLZE4iAAAAE00mEQAAYATmJAIAADCxZBIBAAAWqbUyJxEAAIDJJZMIAAAwgt0yiQAAAEwqmUQAAIBFakmmVTcFAABgUskkAgAALFqZkwgAAMDkkkkEAABYpJZkupmTCAAAwISSSQQAABjB7gnLvU3WqwEAAGAkMokAAACL1FLmJAIAADC5ZBIBAABGMD1hubfJejUAAACMRCYRAABgkVpLdpuTCAAAwKSSSQQAABiB6qYAAABMLEEiAAAAM9xuCgAAsEgtlek2Wbm3yXo1AAAAjEQmEQAAYAS7o3ANAAAAE0omEQAAYJFaLIEBAADABJNJBAAAWDTVTQEAAJhgMokAAAAjmFbdFAAAgEklkwgAALBIrSW7VTcFAABgUskkAgAAjEB1UwAAACaWTCIAAMAitVSmzUkEAABgUo0tSKyqU6vqH6rq4qr6SlX9asc+VVX/p6ouqaovVtVDxjFWAACAPtOpZXksl3HebjqV5EWttc9X1VFJLqyqj7XWvjprn8cnOX34+I9JXj/8CgAAwAEwtiCxtXZFkiuG399cVRcnuVuS2UHik5O8rbXWkny6qo6tqpOHxwIAAIxVS8xJPBCq6h5JHpzkM3ttuluSy2Y93zps2/v4c6pqS1Vtueaaaw7UMAEAACbe2KubVtWGJH+Z5IWttZv23txxSNunobVzk5ybJJs3b95nOwAAwIFincQlVFVrMwgQ39Fa+6uOXbYmOXXW81OSXL4cYwMAADgUjS2TWFWV5M1JLm6t/UHPbu9P8oKqencGBWtuNB8RAABYMdrkrZM4zttNfyDJc5J8qaouGra9NMlpSdJae0OSDyZ5QpJLkmxP8jNjGCcAAMAhY5zVTf853XMOZ+/TkvzS8owIAABgYVqyrGsYLofJmmEJAADASMZe3RQAAOBgNmlzEmUSAQAAmCGTCAAAsEgtMokAAABMMJlEAACAEUxaJlGQCAAAMCGq6ttJbk6yO8lUa21zVR2f5D1J7pHk20me3lq7oe8cbjcFAABYpJbKdFuexwL8cGvtjNba5uHz30hyQWvt9CQXDJ/3EiQCAABMticnOW/4/XlJnjLXzm43BQAAGMF0lm1O4saq2jLr+bmttXP32qcl+WhVtSRvHG4/sbV2RZK01q6oqrvM1YkgEQAA4OBw7axbSPv8QGvt8mEg+LGq+tpCO3G7KQAAwIRorV0+/Hp1kr9O8rAkV1XVyUky/Hr1XOcQJAIAACxWy4opXFNVR1bVUXu+T/KYJF9O8v4kZw93OzvJ++Y6j9tNAQAAJsOJSf66qpJBrPfO1tqHq+pzSd5bVc9L8p0kT5vrJIJEAACARWrJQpenOGBaa99K8qCO9uuSnLW/53G7KQAAADNkEgEAAEawUjKJS0UmEQAAgBkyiQAAAIvUsn+VRw8mMokAAADMkEkEAAAYQZNJBAAAYFLJJAIAAIxgOjKJAAAATCiZRAAAgEVqzTqJAAAATDCZRAAAgBGobgoAAMDEkkkEAABYtDInEQAAgMklkwgAADACcxIBAACYWDKJAAAAi9RinUQAAAAmmEwiAADAYrWktXEPYmnJJAIAADBDJhEAAGAE0zEnEQAAgAklkwgAALBILdZJBAAAYILJJAIAACxaWScRAACAySWTCAAAMALrJAIAADCxZBIBAABGoLopAAAAE0smEQAAYJFak0kEAABggskkAgAAjMA6iQAAAEwsQSIAAAAz3G4KAAAwgtbGPYKlJZMIAADAjLEGiVX1lqq6uqq+3LP9zKq6saouGj5esdxjBAAAmEtrtSyP5TLu203fmuR1Sd42xz7/1Fp74vIMBwAA4NA21iCxtfaJqrrHOMcAAACwWC3Lm+VbDgfDnMRHVNUXqupDVfWAcQ8GAABgko37dtP5fD7J3Vtrt1TVE5L8TZLT996pqs5Jck6SnHbaacs7QgBgxXj0qqct6riPTf/FEo9k+c312ifh9cFKNmHFTVd2JrG1dlNr7Zbh9x9MsraqNnbsd25rbXNrbfOmTZuWfZwAAACTYkVnEqvqpCRXtdZaVT0sg6D2ujEPCwAAYKBl4uYkjjVIrKp3JTkzycaq2prkN5OsTZLW2huSPDXJ86tqKsmOJM9obdKWqgQAAFg5xl3d9JnzbH9dBktkAAAArEwTlsZa0XMSAQAAWF4rek4iALC0JqX652Jfxzj7X6prOO7XfqDN9/pW2mcREnMSAQBWtEM5iDiUXzuwdASJAAAAI5i00prmJAIAADBDJhEAAGCRWiZvTqJMIgAAADNkEgHgIDVXFchas/aA9//YDWd3tk9v397ZPt+YPrrznfu0LaaS52PWPauzvU3tWtC4+vZfSstRyVO1UDjAWpIJyyQKEgHgEFJr1nYGYweb5QiCV6pJD+om/fXBwUCQCAAAMALVTQEAAJhYMokAAACjkEkEAABgUskkAnBQOJQrNC6mwudSnm/V+vVL2n+fvmqpk2Ix7+NCj1nMz0FfHws911xjneSfT0hq4tZJFCQCwARqU7s6fzFfyoBz1ZELCx77lsZI+gPRvva5ztVnoUtdLGbJjj5LHegvxFwB2jjHBaxcgkQAAIBRmJMIAADApJJJBAAAWKyWiZuTKJMIAADADJlEAACAUUzYnERBIsAYKRl/YC115caFvCePWfesBZ+/r8rmYi1mDIxm3Ne8r/+FVGIFDm5VtTrJliTfba09sarumeTdSY5P8vkkz2mt7ZzrHIJEAA4Kh3LQvPrYYzvbd2/b1ntM3/INfe2rjzum91zTty5suYlat657w1xLYPQsp7H7hhu791/E2o19y2bMt9TFUqg1azsDtccd/3O9x8z1/nb9PCwmQF3Kn6tD+WcUkhU1J/FXk1yc5Ojh89ck+cPW2rur6g1Jnpfk9XOdwJxEAACACVBVpyT50SR/OnxeSX4kyfnDXc5L8pT5ziOTCAAAMIqVMyfxtUn+W5Kjhs9PSLKttTY1fL41yd3mO4lMIgAAwMFhY1VtmfU4Z8+GqnpikqtbaxfO2r/rPth5Q1qZRAAAgFEsXybx2tba5p5tP5Dkx6rqCUkOz2BO4muTHFtVa4bZxFOSXD5fJ4JEAMZioZVd5yvKsRzVG7vGvBxFT1iZ5io6sxT7z2WhRWr6ft76iiIlyYevf9OC+gDGq7X2kiQvSZKqOjPJi1trP1VVf5HkqRlUOD07yfvmO5cgEWCMVAM8sPqu72ICzoUup1HrFh489i2BMX3LrZ3tc/2C31cZc81dT+7uY1t3FdEkWXVaz/SV627oPWYhfSdJjji8s3n1+iM629uNN3e37+yv6t5b8bXnveqruJokbddUd3vPe7iUFWrnshSfXWDi/HqSd1fVK5P8a5I3z3eAIBEAAGCxWpK2opbASGvt40k+Pvz+W0ketpDjFa4BAABghkwiAADACNrKWQJjScgkAgAAMEMmEYCJ1VewY76CIAp9QP/PwTgLbi20KjIsmwnLJAoSARiLcf5C11eBcjFWbTiys716qnImSdu+o7N9zbHHLGj/Wtv/z/iau5/au63Lrkfcv3fb2htu62xfleO6z3XiUZ3t667qrkg6OKinWujRPdd3Tc9rn6Pial+V1qmvfaN7/573I0nveOeqONvZx/o5Kqju7P+cLmTJl7n+KDJ9y60qogL7ECQCAACMYoVVNx2VOYkAAADMkEkEAAAYQU3YnESZRAAAAGbIJAIAACxWi+qmALAUFlo9cb5lK+BQspCfH0u+AAslSATgoNCmdi142Yy+X47nWgKj75jVx82xHMJCnXpSd/t1N3Y273rQvTrbV+2a7u1ix10O62w/7Ibu1z51xOrec63tWVViZ89SF31uvXf3khlJcsTlt3a233LP7j7WX9G9LMecv9js6FnK4yEP6Gxvu3b3nqp2dV/H3Td0v4d9n5+5ljFp67qX2di9bVvvMZ3nmeNn52Bbd3AljgmSUt0UAACAySWTCAAAMIoJm5MokwgAAMAMmUQAAIBRTFgmUZAIwMjmq464VMUm+vpZs2njkpwfJtViKpj2HaN4DEw+QSIAK8qq9es726e3b1/wufqqmPb1kSS1bl33hvVHdDbvPvHYzvadxx3e28cRl1zb2T51jxM722/b2D2mW0/qnzWyu+dlHDfV/efu6+4/x68E9++uzHn0Zd3VVdfe0t2+44T+Cqq7NnRXMT36327ubJ8+vLsK7c3ff5fePo78zi2927rUTd0VV5MkJ3RXal29tntc09u6q57Ope3s/vzOtaTFR3e+c582S1zAMpiwTKI5iQAAAMyQSQQAAFisFuskAgAAMLlkEgEAAEZQEzYnUZAIwH5bbAEMhTNgefhZA5aCIBGAkc1VbbGvwmjSXUr/Meue1bnvmrufuuBx9R1z+737K2Cuu7q7ouWOu3ZX31x3w22d7Tef0n9NrnjEyZ3tG7/UUxV0Y/fskNZfLDTH/PvuzvZrv6/7n/7q3j1JMn1Yd/tVD+ueg3PkZd197O45T5IctbW7/dbTNnS291VKPXxb9zVMkts3dVe1Peyansq5R/RXqM2O7ve995ie6qZzVe2dqwrvR245b5+2xQaIXT+Hgk1YoAnLJJqTCAAAwIyxBolV9Zaqurqqvtyzvarq/1TVJVX1xap6yHKPEQAA4FAy7lUxlNsAACAASURBVEziW5M8bo7tj09y+vBxTpLXL8OYAAAADlnzBolV9YKqOu5AdN5a+0SS6+fY5clJ3tYGPp3k2KrqnsgBAAAwBtWW57Fc9qdwzUlJPldVn0/yliQfaa0t1xDvluSyWc+3DtuuWKb+AVacuQpKdBWgWMx55ipEA4zfwVRYZqn+nwUsn3mDxNbay6rq5Ukek+Rnkryuqt6b5M2ttW8e4PF1lU3bJ0CtqnMyuB01p5122gEeEsChq6/aYtvZX8F0wX1sOLK7j6O725Nk+ojuoLam+itd9tl+96M722+6e/c/maumuvu+bY57cNZ1F7rM1Q/tvsFn6sju1zG9vr8k6U3f0139s62eo4xpjyMv7T7XEVd0Vzed6nmrbtvU/360Nd2vfdXt3X1vuKL779VTR3SPKUnWX9H9OV112VXdY9o11XuuOqH7DZ765r93tq8+9tju88zxszNX5dOFWmgwJniDBWr9/+85GO3XnMRh5vDK4WMqyXFJzq+q/3UAx5YMMoez65efkuTyjvGd21rb3FrbvGnTpgM8JAAAgMm1P3MSf6WqLkzyv5L8S5Lvb609P8lDk/zEAR7f+5P89LDK6cOT3Nhac6spAACwMrRlfCyT/ZmTuDHJf26tXTq7sbU2XVVPHKXzqnpXkjOTbKyqrUl+M8na4fnfkOSDSZ6Q5JIk2zO43RUAAIADZH/mJL5ijm0Xj9J5a+2Z82xvSX5plD4AAAAOqGXM8i2Hca+TCAAAwAqyP7ebAnCQ6Cs1v1IrFT7+nv913EOAQ4JlKODAWs41DJeDIBHgINP3C91C102bay3EvmUokuTD179pn7bHn/Irvfu37Tu6+z/mqO4DTj2ps3n3UYf39rH65ts623fepft11FT/v+a3ntj9T+PNd1/YbwCrb+8vh3778d3nus/DLu1sv+Sq7srdU1ce0dvHLz3qo53tf3zhmZ3t67/Uf323n7KwJTiO/Gb3Z+u4r/bfwHTzj9za2X78B7qXXVmzo/sarr+i+7OQJNc8uPtcJ1/V/Vmsm27uPVefvqUupm/pfn21rv/nsE0tfGmZrv8/jHtNRUEoHHzcbgoAAMAMmUQAAIBRTNjtpjKJAAAAzJBJBAAAGMWEZRIFiQDLYCkrC467CMVjN5y9T9uqY48Zw0gAgANBkAjAPtrOnQs74OgN/dtO6Akgd3VXxuyrSLp6R/f+SXLzfborSt560sJnVdx+XHf7Ufe/obN9x4XHd7b/8I9+vrePG3f1VyXt8msPvKCz/bv37Rlskvde+tDO9p8541Od7Zffvz/Q/5fv3quz/eat3VVB+6q31u7+iq9H/X33+37ENd0VPndtWN3dflR/tdCTL7ims336mO6qp6vmqG7a1h/W3X5F989OX8Xg3du29faRLE1lUNVF4cCqNnlLYJiTCAAAwAyZRAAAgFG0/jslDkYyiQAAAMyQSQQAABjFhM1JFCQCHAK6KqLWmv4CH8DKtdBqyY9Z96ze/dtUd2GgvnMBhwZBIsAyWMpftpZqyYxa1x8kTm/f3tm+an13Fcjp73y3fwD3u2f3MUd0999XxfSaB3f3nSQ7u4tsZvsp3edad0N3Zcwkafe5tbP9fsde39n+rxu7q4L+6HFf6O3j/Gv/Q2f7Uzd+rrP907fcu7P9ul3dFTOT5P7HX9XZfsq67texcW1/Jc8dJ63rbP/HG+7T2b665/pO9Q83rectWXtr968qR39rR3ffN9/W30lf32sWPvumrry2e8O67ms1fUv358ofa2AyqG4KAADAxJJJBAAAGIVMIgAAAJNKJhEAAGCx2uTNSRQkAoeEuaoBJktXWGauKoIf3fnOJemj77WsPvbYJTk/sHLM9/8ugANBkAgw4dbc/dTO9ulrrus/5nu6K5Lmpp4KmEf3lBdNcvtxh3e2Tx3RXc7y9mO72w+7sf/PtLcfV53tR17afa6+SppJ0rcgwNFru6tm/vgjt3S291UwTZI/O+2fOttfdOVDOtv/+rMP7Wx/wX/6u94+3vjl/9TZvvXW7mqs37n2+N5znXjcTZ3tdVP3rxG7D+t+rzZ+sbeLHLatuxJtn1W3db9TbW3/m9tbUXdrd6XStmuq/1w91UpXbegu4bqYpSYWutRFn6X6AxUwhwnLJJqTCAAAwAyZRAAAgFHIJAIAADCpZBIBAABGMGnVTWUSAQAAmCGTCNBjMaXna0139cSke3mMWte//0duOW/B/QMHn3Euc7Ecy/YABx9BInBIWKp1EOfTpnZ19jXXL2ILtWr9+u6+d+7sPuDo7pL8q+Yo758d3cs95ITjune/e/8ajatv617a4NYTu/8JOnzbdHcfG/tvflnTvRpBdnav9pDbT+p/7af0LPdw067upTw+f+Upne27pvqXYnj54Td0tn/71hM624866ZbO9i/d3N13kuy8aV1n+3fSvdTFrhsO6z3X1hs2dbavu6n7PbnrP3cv97Duhp7PVZLdPctTrNvafa0y1fMeHtH9PiXp/Vzvvrp7CYy5/ojTt63t3Nn5Bx7rHQIHE0EiAADAKMxJBAAAYFLJJAIAACxWU90UAACACSaTCLBCPXbD2eMeAgCwPyYskyhIBOjRVxF1MVUKF1pKfjEB4qoju6ue7r74m53t9cD79J/rsqs623eeeFRn++GX39x7ruse0l0RtU9fFdPb5jjNjpO7K6Ku2nR7Z/uxR+3oPddVNxzd31GH1au7+z5mfX8fW3tezBcvu1tn++7buyulfmbq7r19rDmyu/rn1JVHdLavv6L/5qK127vb11/d/dr7KtrOZe313Z3sPmFDdx9buyuSthuv6e1j+paeMriLtJBlaparwvJCWeYCJkdVHZ7kE0kOyyDOO7+19ptVdc8k705yfJLPJ3lOa62nJPqA200BAABG0ZbpMbfbk/xIa+1BSc5I8riqeniS1yT5w9ba6UluSPK8+U4kSAQAADjItYE9C+uuHT5akh9Jcv6w/bwkT5nvXIJEAACARaoMqpsux2PesVStrqqLklyd5GNJvplkW2ttzxyErUm65zbMIkgEAAA4OGysqi2zHufM3tha291aOyPJKUkeluT+HeeYN9xUuAZgiS2msM3eVq3vLkIDAKxAy1fd9NrW2ub5dmqtbauqjyd5eJJjq2rNMJt4SpLL5ztekAiwAk1v7yknmWT1scd2tu++4cbO9lUbjuxsrxv7+2gnbexsX3dVdxXTtra7+maSHP3t7gqj209e19k+dVh1tu86urs9SabXd1fT3PC57kqetz68uypnkpx1r693tn/4wgf2HtNl3Wk39G775088oLN9w9bu1zjdc3l3H9Z9DZNkffdl77Xpi/0HTK/pvvFo1VT3dVz37e4Ko+3G/iq4tbb7V5LVa9d2j2lb9+d9Lm1qV++2ruqjlqEBDiZVtSnJrmGAeESSR2VQtOYfkjw1gwqnZyd533znEiQCAAAs1n7OF1wGJyc5r6pWZzCt8L2ttb+tqq8meXdVvTLJvyZ583wnEiQCAAAc5FprX0zy4I72b2UwP3G/KVwDAADADJlEAACAUayM202XjEwiAAAAM2QSgUNe35IVtaa7qiLAQs21NE5XZVXgIDNhmURBIsAymCvg/OjOd+7T9rjjf653/+lbbl1Q3337rz7mqN5j6qapzvZ2dPdyGnNZc3P30gpHrOle7mHbvQ/rbD/5k91jSpLdh3WvEXF91xLCSVZd3P86Pvrdfeb8J0nW3tY93t2Hd/9mcO3XupcRSZK7frJ76Yi2uruPW0/qvvFn3U29XWTdrd19TPdc975lLpLkiC9t7d6wq2dJiaO7P1u1vntJkkEnh3c27760p+9l0Hb2L5kBMMkEiQAAACNYIUtgLBlzEgEAAJghkwgAADAKmUQAAAAmlUwiAHDI66o+ulwVjlU+hYNcy8RlEgWJAD3a1K4l+wXtMeuetbC+d+7s3zbVXXFxzabuaprTt27vPlFPNckkyY7b+rd1qO3dFUyT5Jbvv0tn+/X3665IuvHL3VVMd23ov/llxwnd20748u7uMd2tu+8kOXFL9zE7Tug+Zs3t3b8Z9FURnav/w27sPtexlyy8yubhl9/c2V439VTHnevz0KPt6qk4e90Nnc1zVeZdfZfuz2/f573PXIFdrVnbWU24L0hbyv8H9J1nrgARYFwEiQAAACNQ3XQJVdXjqurrVXVJVf1Gx/bnVtU1VXXR8PGz4xgnAADAoWJsmcSqWp3kj5M8OsnWJJ+rqve31r66167vaa29YNkHCAAAsD9kEpfMw5Jc0lr7VmttZ5J3J3nyGMcDAABwyBvnnMS7Jbls1vOtSf5jx34/UVU/mOTfkvxaa+2yjn2Ag9h8hRsOpup+fa9lrmIajzv+5w7UcACAZTBpcxLHGSR2lX3b+/J+IMm7Wmu3V9UvJDkvyY/sc6Kqc5KckySnnXbaUo8TmHDLEYT2VUnsCxCnt/dUJE2yav36zvapa65d0P59FSiTJOuP6GyuXd2VP2+/x/G9pzrqs9/pbv+X7qqV7aTuKpd9fSfJ+hOP6mxfvaO7j3U3r+s9V5+NH/xWZ/vtD7x7Z/thX7y0/2Rre/5ocPSG7vae96rvWiVJLruyp+/uf/rb0Uf2nqptu7Gzve9zupilI6Z7+li1fn0+cst5Cz7fQozzD1EH0x/BgEPHOG833Zrk1FnPT0ly+ewdWmvXtdb21FV/U5KHdp2otXZua21za23zpk2bDshgAQAAOrVleiyTcQaJn0tyelXds6rWJXlGkvfP3qGqTp719MeSXLyM4wMAADjkjO1209baVFW9IMlHkqxO8pbW2leq6reTbGmtvT/Jr1TVjyWZSnJ9kueOa7wAAAD7WOYs33IY55zEtNY+mOSDe7W9Ytb3L0nykuUeFwAAwKFqrEEiwFKbq1LquAtEPGbds/ZpW7Whv1gIALDyVborch7MBInA2I07eFsOfdUep2+5NR/d+c592ucKdtvO7oqdazZ1V7qcvrW/UmqvXVPd7d1FT7Pu6lv7z9VXsbNHW7u6s72u7K7emiS9tUpvuqW7j3uc2Huu+uxXOtun13W/h4d996buE/VVME1/Jc/avqO7vafa7PTF3+zto+8PEH2fhzbHuZJ0fk77PHbD2d19z1G1t011f64BWH6CRAAAgFFM2JzEcVY3BQAAYIWRSQQAABhBySQCAAAwqWQSAQAARjFhmURBIsCYdS2NAStNV8Xdvqq91VMJFoCDgyARmCjjXk6jr/++JS16f8nuaU/6lwpY6FIXU9f0Lymx5q4nd/d95TWd7bu3bes/1/1O7z7mkm93tveuHXnCcb199I2rb+mIVRf9W++50hPgrDr2mM72vtcxlwUHUbu63/O5lo3YvW1b5+dxMX+UWOjyFG3nrgUv7ZKM/+cXYNEmLJNoTiIAAAAzZBIBAAAWq6luCgAAwASTSQQAABjFhGUSBYkAPeYrstFH8Q1Wisef+PxxDwGAg5AgEWAZ9FUrnatq5EKDzb6qlX19rFq/vvdcu6/ur3zaZc3dT+3d1lt5tKfCZ9u5s/tE193Q20cdc1Rn++7vXtl7zEJNXX5FZ/vqY49d1Pk+fP2b9nvfvvd2rvdw1ZHd21Yfd0w+dNXr97vvpP8PJgutegowqcxJBAAAYGIJEgEAAJjhdlMAAIBRuN0UAACASSWTCHCQ6Ssi0lcch5WprxjNR3e+c5lHMrDYar4ATF7hGkEiQI+lXMqir5JnrVubj9xy3pL00aZ2HfDlN/oCm6lLL+s9pi94XbXhyM723du2dZ9o+/b+gfUc01d5tLeP9FcM7avk2XeupQza+/qe6z1f6qDP0i4Ahw5BIgAAwGK1mJMIAADA5JJJBAAAGIVMIgAAAJNKJhGYKHMV61iphTf6xrxU413MNXnshrOXpG8WTpVRgINLRXVTABahr4LpYgKCcQa7fVU2F3PM7m3bOl/LYoLaxx3/c73HfPj6N80zwtH679t/vmu1VK99qSzH52ql/qEGgDsTJAIAAIxiwjKJ5iQCAAAwQyYRAABgBNUmK5UokwgAAMAMmUQAAIDFapm4OYmCRIAVqq/SZa1Z23vMR3e+80AN54Cw3AMArDzVJuz+2c2bN7ctW7aMexgAI5kreBpnkLiYJRoOtmU+AFh+VXVha23zuMexGEduPLV974/92rL0teXPXrQs18mcRAAAAGa43RQAAGAUk3VzpkwiAAAAd5BJBAAAGEFNWCZRkAisaIsplAIHgs8iAIcKQSLACrRSg46lHtdKfZ0AsCATlkk0JxEAAIAZMokAAACL1SZvTqJMIgAAADNkEgEAAEYhkwgAAMCkkkkEVjTVL1kpfBYB6FIxJxEAAIAVpqpOrap/qKqLq+orVfWrw/bjq+pjVfWN4dfj5juXIBEAAGAUrS3PY25TSV7UWrt/kocn+aWq+t4kv5Hkgtba6UkuGD6fkyARAADgINdau6K19vnh9zcnuTjJ3ZI8Ocl5w93OS/KU+c5lTiIAAMAIVtqcxKq6R5IHJ/lMkhNba1ckg0Cyqu4y3/GCRAAAgIPDxqraMuv5ua21c2fvUFUbkvxlkhe21m6qqgV3IkgEYB+PXvW03m2qfALALC3LuU7ita21zX0bq2ptBgHiO1prfzVsvqqqTh5mEU9OcvV8nQgSATigBJUAcODVIGX45iQXt9b+YNam9yc5O8mrh1/fN9+5BIkAAAAjqOlxjyBJ8gNJnpPkS1V10bDtpRkEh++tqucl+U6S/tuFhsYaJFbV45L8UZLVSf60tfbqvbYfluRtSR6a5LokP9la+/ZyjxMAAGAla639c5K+CYhnLeRcY1sCo6pWJ/njJI9P8r1Jnjlcx2O25yW5obV27yR/mOQ1yztKAACAQ8s4M4kPS3JJa+1bSVJV785gDY+vztrnyUl+a/j9+UleV1XVWv9Kkl//+tdz5pln3qnt6U9/en7xF38x27dvzxOe8IR9jnnuc5+b5z73ubn22mvz1Kc+dZ/tz3/+8/OTP/mTueyyy/Kc5zxnn+0vetGL8qQnPSlf//rX8/M///P7bH/Zy16WRz3qUbnooovywhe+cJ/tr3rVq/LIRz4yn/zkJ/PSl750n+2vfe1rc8YZZ+Tv/u7v8spXvnKf7W984xtz3/veNx/4wAfy+7//+/tsf/vb355TTz0173nPe/L6179+n+3nn39+Nm7cmLe+9a1561vfus/2D37wg1m/fn3+5E/+JO9973v32f7xj388SfJ7v/d7+du//ds7bTviiCPyoQ99KEnyO7/zO7ngggvutP2EE07IX/7lXyZJXvKSl+RTn/rUnbafcsop+fM///MkyQtf+MJcdNFFd9p+n/vcJ+eeOyjodM455+Tf/u3f7rT9jDPOyGtf+9okybOf/exs3br1Ttsf8YhH5Hd/93eTJD/xEz+R66677k7bzzrrrLz85S9Pkjz+8Y/Pjh077rT9iU98Yl784hcnyT6fu8Rnz2fv4P3sXdduy0l1am5r2/PlfPZO284880yfPZ+9JP6/57Pnszebz95on72D3gpbAmNUY8skZrCw42Wznm8dtnXu01qbSnJjkhP2PlFVnVNVW6pqy65duw7QcAEAACZfzZGUO7AdVz0tyWNbaz87fP6cJA9rrf3yrH2+Mtxn6/D5N4f7XNd1ziTZvHlz27JlS99mAABghamqC+da2mEl23D8qe1BZ+2bPT4QPnn+i5flOo0zk7g1yamznp+S5PK+fapqTZJjkly/LKMDAAA4BI0zSPxcktOr6p5VtS7JMzJYw2O2PWt6JMlTk/z9XPMRAQAAllVL0tryPJbJ2ArXtNamquoFST6SwRIYb2mtfaWqfjvJltba+zNYDPLtVXVJBhnEZ4xrvAAAAIeCsa6T2Fr7YJIP7tX2ilnf35b9WOwRAABgXGrC7nUc5+2mAAAArDBjzSQCAAAc9GQSAQAAmFQyiQAAAItUMScRAACACSaTCAAAsFjLvIbhcpBJBAAAYIZMIgAAwAjMSQQAAGBiySQCAACMQiYRAACASSWTCAAAMAJzEgEAAJhYMokAAACL1ZJMT1YqUSYRAACAGTKJAAAAo5isRKJMIgAAAHeQSQQAABiB6qYAAABMLJlEAACAUbTJSiXKJAIAADBDJhEAAGAE5iQCAAAwsWQSAQAAFqvFOokAAABMLplEAACARaokpbopAAAAk0omEQAAYBTT4x7A0pJJBAAAYIZMIgAAwAjMSQQAAGBiCRIBAACY4XZTAACAxWrDxwSRSQQAAGCGTCIAAMCitUThGgAAACaVTCIAAMAIarISiTKJAAAA3EEmEQAAYBTmJAIAADCpZBIBAAAWqyU1Pe5BLC2ZRAAAAGbIJAIAAIzCnEQAAAAmlUwiAADAKCYrkSiTCAAAwB1kEgEAAEZQ5iQCAAAwqWQSAQAARiGTCAAAwKSSSQQAAFislmR63INYWjKJAAAAzJBJBAAAWKRKU910KVTV8VX1sar6xvDrcT377a6qi4aP9y/3OAEAAA4147rd9DeSXNBaOz3JBcPnXXa01s4YPn5s+YYHAACwn1pbnscyGVeQ+OQk5w2/Py/JU8Y0DgAAAGYZV5B4YmvtiiQZfr1Lz36HV9WWqvp0VQkkAQCAlWfCMokHrHBNVf1dkpM6Nv33BZzmtNba5VV1ryR/X1Vfaq19s6Ovc5KckySnnXbaosYLAADAAQwSW2uP6ttWVVdV1cmttSuq6uQkV/ec4/Lh129V1ceTPDjJPkFia+3cJOcmyebNmyertBAAALByWSdxybw/ydnD789O8r69d6iq46rqsOH3G5P8QJKvLtsIAQAADiJV9Zaqurqqvjyrbb9WlphtXEHiq5M8uqq+keTRw+epqs1V9afDfe6fZEtVfSHJPyR5dWtNkAgAAKwo1dqyPPbDW5M8bq+2/V1ZYsYBu910Lq2165Kc1dG+JcnPDr//ZJLvX+ahAQAAHJRaa5+oqnvs1fzkJGcOvz8vyceT/Ppc5xlLkAgAADAxlrHy6CLcaWWJqupbWWKGIBEAAODgsLGqtsx6fu6wiOeSEiQCAAAs2rKuYXhta23zAo/Zr5UlZhtX4RoAAAAOvHlXltibTCIAAMBitayYOYlV9a4MitRsrKqtSX4zg5Uk3ltVz0vynSRPm+88gkQAAIAJ0Fp7Zs+mfVaWmIsgEQAAYBTT4x7A0jInEQAAgBmCRAD4f+3db6j2d10H8Pd795T5b0Qt0rZBC0QQH6TYyAZiZjUtkqLApIKe7EkTBSPUHkgRRE+iIAmG2h9Sh8xGUss1qRgS1ZxZejsHYwi7XbGipCn4Z55PD87F1U3dW9s593X97vt3vV5wwfld5+K63vDlcM7nfH6f7xcA2HK7KQAAwCn0Etm45mLRSQQAAGBLJxEAAOA0dBIBAABYK51EAACAk5okRzqJAAAArJROIgAAwImNmUQAAADWSycRAADgNHQSAQAAWCudRAAAgNPQSQQAAGCtdBIBAABOyjmJAAAArJlOIgAAwIlNMkdLh7iodBIBAADY0kkEAAA4DbubAgAAsFY6iQAAACdld1MAAADWTCcRAADgNMwkAgAAsFY6iQAAAKehkwgAAMBa6SQCAACc2OgkAgAAsF46iQAAACc1SY6Olk5xUekkAgAAsKWTCAAAcBpmEgEAAFgrnUQAAIDT0EkEAABgrXQSAQAATmySI51EAAAAVkonEQAA4KQmmXFOIgAAACulSAQAAGDL7aYAAACnYeMaAAAA1konEQAA4DRGJxEAAICV0kkEAAA4qZnkyBEYAAAArJROIgAAwGmYSQQAAGCtFikS2/5027Ntj9q+8iled3PbB9s+1PYd+8wIAADwdMzR0V4e+7JUJ/GzSX4yyb1P9oK2Z5K8J8nrk7w0yc+0fel+4gEAABymRWYSZ+aBJGn7VC+7MclDM/Pw5rW3J3ljks/tPCAAAMDTMmYS9+jaJI+cd31u8xwAAAA7srNOYtuPJ3nhBb71KzPzp0/nLS7w3AVL9La3JLllc/m1tp99eilZkWuS/PvSIdgra36YrPvhseaHybofnpcsHeDEJsnRujqJOysSZ+Z1p3yLc0muP+/6uiSPPsln3ZbktiRp+8mZedLNcFgn6354rPlhsu6Hx5ofJut+eNp+cukM/I9L+ZzE+5K8uO0NSb6Y5E1J3rxsJAAAgP9l9rfz6D4sdQTGT7Q9l+RVSf687d2b57+z7V1JMjNPJLk1yd1JHkjy4Zk5u0ReAACAQ7HU7qZ3JrnzAs8/muQN513fleSuZ/j2t50uHZcp6354rPlhsu6Hx5ofJut+eC7bNZ8ks7KZxM7KtmsFAADYl6uv+Lb5vit/ZC+fdc83PnT/PuZ1L+WZRAAAgEvbjJnES1nbm9s+2Pahtu9YOg+71/b9bR9z7MnhaHt9279u+0Dbs23funQmdq/tVW3/oe0/bdb9V5fOxH60PdP2H9v+2dJZ2I+2X2j7mbaftuPlYWj7LW3vaPv5ze/3Vy2d6dCtppPY9kyS9yT5oRwfn3Ff24/OzOeWTcaO/UGS303yRwvnYH+eSPL2mflU2xckub/tPX7WV+9rSV47M19u+6wkn2j7FzPzd0sHY+femuMN7K5eOgh79QMz45zEw/E7ST42Mz/V9tlJnrt0oGdqbTOJa+ok3pjkoZl5eGa+nuT2JG9cOBM7NjP3JvmPpXOwPzPzLzPzqc3Xj+f4j8drl03Frs2xL28un7V5rOs3Mv9H2+uS/GiS9y6dBdiNtlcneXWS9yXJzHx9Zr60bCrWVCRem+SR867PxR+OsGptvyvJy5P8/bJJ2IfNbYefTvJYkntmxrqv328n+eUk6xr24f8zSf6y7f1tb1k6DDv33Un+Lcnvb24tf2/b5y0d6hmbo/089mRNRWIv8Jz/MsNKtX1+ko8kedvM/NfSedi9mfnmzHxPkuuS3Nj2ZUtnYnfa/liSx2bm/qWzsHc3PFXEhQAAAnBJREFUzcwrkrw+yS+2ffXSgdipK5O8IsnvzczLk3wlib1FFraamcQcdw6vP+/6uiSPLpQF2KHNTNpHknxgZv5k6Tzs18x8qe3fJLk5iU2r1uumJD/e9g1Jrkpydds/npmfXTgXO7Y5Nzsz81jbO3M8UnTvsqnYoXNJzp13d8gducyKxMfzn3d/fO64Zk8ft5dZ3TUVifcleXHbG5J8Mcmbkrx52UjAxda2OZ5beGBmfmvpPOxH229P8o1NgficJK9L8psLx2KHZuadSd6ZJG1fk+SXFIjrt7nN8IqZeXzz9Q8n+bWFY7FDM/OvbR9p+5KZeTDJDya5rDajm5mbl85wsa2mSJyZJ9remuTuJGeSvH9mzi4cix1r+6Ekr0lyTdtzSd49M+9bNhU7dlOSn0vymc18WpK8a2buWjATu/eiJH+42cn6iiQfnhlHIsD6fEeSO4//H5grk3xwZj62bCT24C1JPrDZ2fThJL+wcJ6D1xljewAAABxb08Y1AAAAnJIiEQAAgC1FIgAAAFuKRAAAALYUiQAAAGwpEgEAANhSJAIAALClSATgktb2e9v+c9ur2j6v7dm2L1s6FwCsVWdm6QwA8JTa/nqSq5I8J8m5mfmNhSMBwGopEgG45LV9dpL7knw1yffPzDcXjgQAq+V2UwAuB9+a5PlJXpDjjiIAsCM6iQBc8tp+NMntSW5I8qKZuXXhSACwWlcuHQAAnkrbn0/yxMx8sO2ZJH/b9rUz81dLZwOANdJJBAAAYMtMIgAAAFuKRAAAALYUiQAAAGwpEgEAANhSJAIAALClSAQAAGBLkQgAAMCWIhEAAICt/wavgVIPcZz8xgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(14, 8))\n", "counts, xedges, yedges, im = ax.hist2d(x1_all, x2_all, bins=[120, 80], range= [[0.0, 6.0], [-1.0, 3.0]], cmin=1)\n", "ax.plot([0.0, 6.0], [0.0, 0.0], \"--k\") # NOTE: This draws a line from [x1, x2], [y1, y2] with dashed line (\"--\") and in black (\"k\")\n", "fig.colorbar(im) # ticks=[-1, 0, 1]\n", "\n", "ax.set(title='Histogram of lengths (x) and widths (y)',\n", " xlabel='x', \n", " ylabel='y',\n", " aspect='equal', # NOTE: This forces the x- and y-axis to have the SAME scale!!!\n", " )\n", "\n", "fig.tight_layout()\n", "fig\n", "\n", "if save_plots :\n", " fig.savefig(\"Dist_2Dgauss.pdf\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we bin `y_all` and fit it with a Gaussian distribution:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# In case you want to see if it is really Gaussian, you can plot (or fit) with this function:\n", "def gaussian(x, N, mu, sigma):\n", " return N * 1.0 / (sigma*np.sqrt(2*np.pi)) * np.exp(-0.5* (x-mu)**2/sigma**2)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAFlCAYAAAApo6aBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAacUlEQVR4nO3df7Cld10f8PfHJMBOod5CVpreZLt0TFvRqYFuYzo706HBVkwcglNwYlsIGGdtiy2Otrr4R9VaZtaZKv4sNhpqsCqkqCUlsS2yMI6ZEkwwICFSV9zKbjIkSDbAsNJJ/PSP+0Tu3j039+ze5+w599zXa+bMfc73+Z5zv988u0/e+3y/z/ep7g4AANv3ZfNuAADAshCsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkF8+7AUly6aWX9v79++fdDACALd13332f7u69k/YtRLDav39/7r333nk3AwBgS1X1fzfbZygQAGAkghUAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADASwQoAYCSCFQDASAQrAICRCFYAACMRrAAARnLxvBsALL+DR47m5KnTZ5StruzJ3YevnVOLAGZDsAJm7uSp0zl+5PozyvYfvnNOrQGYHcEKWFqulAEXmmAFLC1XyoALzeR1AICRCFYAACMRrAAARiJYAQCMRLACABiJuwKBhWbJBGAnEayAhWbJBGAnmXoosKouqqrfrap3D+9fUFX3VNUfVNU7quoZQ/kzh/fHhv37Z9N0AIDFci5zrN6Q5MF1738kyZu7+8okjyW5eSi/Oclj3f2VSd481AMAWHpTBauqujzJ9Ul+fnhfSa5N8s6hym1JXjFs3zC8z7D/pUN9AIClNu0cqx9P8r1JnjO8f16SU939xPD+RJLVYXs1ySeTpLufqKrHh/qfXv+FVXUoyaEk2bdv3/m2H9ihVlf2nDVXyqR0YKfbMlhV1TcleaS776uqlzxVPKFqT7HvSwXdtyS5JUkOHDhw1n5guU0KUCalAzvdNFesDiZ5eVVdl+RZSf5i1q5grVTVxcNVq8uTPDTUP5HkiiQnquriJF+e5DOjtxzYtVztAhbVlsGqu9+Y5I1JMlyx+tfd/U+q6r8meWWStye5Kcm7ho/cMbz/38P+o93tihQwGle7gEW1nZXXvy/Jd1fVsazNobp1KL81yfOG8u9Ocnh7TQQA2BnOaYHQ7n5/kvcP259IcvWEOn+a5FUjtA0AYEfxrEAAgJEIVgAAIxGsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkghUAwEgEKwCAkQhWAAAjOadH2gAsqtWVPWc9iHl1Zc+cWgPsVoIVMKqDR47m5KnTZ5RdiIBz9+FrZ/47ALYiWAGjOnnqdI4fuX7ezQCYC3OsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkghUAwEgstwDseputvWVtLOBcCVbArjdp7a2Nq7gDTMNQIADASAQrAICRCFYAACMRrAAARmLyOrAwVlf2nDVpfHVlz5xaA3DuBCtgYVjeANjpDAUCAIxky2BVVc+qqg9W1Yer6oGq+qGh/Beq6o+q6v7hddVQXlX1k1V1rKo+UlUvnnUnAAAWwTRDgV9Mcm13f76qLkny21X1G8O+f9Pd79xQ/xuTXDm8vi7JW4afAABLbcsrVr3m88PbS4ZXP81HbkjytuFzH0iyUlWXbb+pAACLbao5VlV1UVXdn+SRJO/p7nuGXW8ahvveXFXPHMpWk3xy3cdPDGUbv/NQVd1bVfc++uij2+gCAMBimCpYdfeT3X1VksuTXF1VX5PkjUn+ZpK/k+S5Sb5vqF6TvmLCd97S3Qe6+8DevXvPq/EAAIvknO4K7O5TSd6f5GXd/fAw3PfFJP85ydVDtRNJrlj3scuTPDRCWwEAFto0dwXuraqVYXtPkq9P8vtPzZuqqkryiiQfHT5yR5LXDHcHXpPk8e5+eCatBwBYINPcFXhZktuq6qKsBbHbu/vdVXW0qvZmbejv/iT/bKh/V5LrkhxL8oUkrxu/2QCztdkq8BYxBZ7OlsGquz+S5EUTyieeXbq7k7x++00DGN+0j82ZFKA2fg5gI4+0AaZy8MjRnDx1+oyynXgFZ6e1F9hZBCtgKidPnc7xI9efUeYKDsCZBCvgvE07rAawWwhWwHkzrAZwpnNaxwoAgM0JVgAAIxGsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkghUAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADASwQoAYCSCFQDASAQrAICRCFYAACO5eN4NANgpVlf2ZP/hOyeW33342jm0CFg0ghXAlDYLT5PCFrA7GQoEABiJYAUAMBLBCgBgJIIVAMBItgxWVfWsqvpgVX24qh6oqh8ayl9QVfdU1R9U1Tuq6hlD+TOH98eG/ftn2wUAgMUwzRWrLya5tru/NslVSV5WVdck+ZEkb+7uK5M8luTmof7NSR7r7q9M8uahHgDA0ttyuYXu7iSfH95eMrw6ybVJ/vFQfluSH0zyliQ3DNtJ8s4kP11VNXwPwNKZtL6Vta1gd5pqHauquijJfUm+MsnPJPnDJKe6+4mhyokkq8P2apJPJkl3P1FVjyd5XpJPb/jOQ0kOJcm+ffu21wuAOZoUoKxtBbvTVJPXu/vJ7r4qyeVJrk7yVZOqDT/rafat/85buvtAdx/Yu3fvtO0FAFhY53RXYHefSvL+JNckWamqp654XZ7koWH7RJIrkmTY/+VJPjNGYwEAFtk0dwXuraqVYXtPkq9P8mCS9yV55VDtpiTvGrbvGN5n2H/U/CoAYDeYZo7VZUluG+ZZfVmS27v73VX1sSRvr6p/n+R3k9w61L81yS9W1bGsXam6cQbtBgBYONPcFfiRJC+aUP6JrM232lj+p0leNUrrAAB2ECuvAwCMRLACABiJYAUAMJKpFggFdr6DR47m5KnTZ5RZHRxgXIIV7BInT53O8SPXn1FmdXCAcRkKBAAYiWAFADASwQoAYCSCFQDASAQrAICRCFYAACOx3AJwhknrXSVra14B8PQEK+AMk9a7AmA6hgIBAEYiWAEAjESwAgAYiWAFADASwQoAYCTuCoRdbHVlT/YfvvOsMgDOj2AFu9jdh6+ddxMAloqhQACAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADASwQoAYCSCFQDASLYMVlV1RVW9r6oerKoHquoNQ/kPVtXJqrp/eF237jNvrKpjVfXxqvqGWXYAAGBRTLPy+hNJvqe7P1RVz0lyX1W9Z9j35u7+D+srV9ULk9yY5KuT/JUkv1lVf727nxyz4QAAi2bLK1bd/XB3f2jY/lySB5OsPs1Hbkjy9u7+Ynf/UZJjSa4eo7EAAIvsnJ4VWFX7k7woyT1JDib5zqp6TZJ7s3ZV67Gsha4PrPvYiUwIYlV1KMmhJNm3b995NB1gZzl45GhOnjp9Rtnqyh7PbIQlMnWwqqpnJ/nVJN/V3Z+tqrck+eEkPfz80STflqQmfLzPKui+JcktSXLgwIGz9gMsm5OnTuf4kevPKNt/+M45tQaYhanuCqyqS7IWqn6pu38tSbr7U939ZHf/WZKfy5eG+04kuWLdxy9P8tB4TQYAWEzT3BVYSW5N8mB3/9i68svWVfvmJB8dtu9IcmNVPbOqXpDkyiQfHK/JAACLaZqhwINJXp3k96rq/qHs+5N8a1VdlbVhvuNJviNJuvuBqro9yceydkfh690RCADsBlsGq+7+7UyeN3XX03zmTUnetI12AQDsOFZeBwAYiWAFADASwQoAYCSCFQDASM5p5XVgZ9hshW8AZkuwgiU0aYVvLqzVlT1nraou3MLyE6wAZsDz/2B3MscKAGAkghUAwEgMBQLM0WZzsQwlws4kWAHM0aQAtTFoATuHoUAAgJEIVgAAIxGsAABGIlgBAIxEsAIAGIm7AmGH81xAgMUhWMEO57mAAIvDUCAAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADASwQoAYCSCFQDASAQrAICRbBmsquqKqnpfVT1YVQ9U1RuG8udW1Xuq6g+Gn39pKK+q+smqOlZVH6mqF8+6EwAAi2CaK1ZPJPme7v6qJNckeX1VvTDJ4STv7e4rk7x3eJ8k35jkyuF1KMlbRm81AMAC2jJYdffD3f2hYftzSR5MsprkhiS3DdVuS/KKYfuGJG/rNR9IslJVl43ecgCABXNOc6yqan+SFyW5J8nzu/vhZC18JfmKodpqkk+u+9iJoQwAYKlNHayq6tlJfjXJd3X3Z5+u6oSynvB9h6rq3qq699FHH522GQAAC2uqYFVVl2QtVP1Sd//aUPypp4b4hp+PDOUnklyx7uOXJ3lo43d29y3dfaC7D+zdu/d82w8AsDCmuSuwktya5MHu/rF1u+5IctOwfVOSd60rf81wd+A1SR5/asgQAGCZXTxFnYNJXp3k96rq/qHs+5McSXJ7Vd2c5I+TvGrYd1eS65IcS/KFJK8btcUAAAtqy2DV3b+dyfOmkuSlE+p3ktdvs10AADuOldcBAEYyzVAgsCAOHjmak6dOn1G2urJnTq0BYCPBCnaQk6dO5/iR6+fdDAA2YSgQAGAkrljBgjLsB7DzCFawoAz7Aew8hgIBAEYiWAEAjESwAgAYiTlWAAtmdWVP9h++86yyuw9fO6cWAdMSrAAWzKQAtTFoAYvJUCAAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYieUWYAF44DLAchCsYAF44DLAcjAUCAAwEsEKAGAkghUAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADCSLYNVVb21qh6pqo+uK/vBqjpZVfcPr+vW7XtjVR2rqo9X1TfMquEAAItmmitWv5DkZRPK39zdVw2vu5Kkql6Y5MYkXz185j9W1UVjNRYAYJFtGay6+7eSfGbK77shydu7+4vd/UdJjiW5ehvtAwDYMbbzEObvrKrXJLk3yfd092NJVpN8YF2dE0MZANuwurIn+w/feVbZ3YevnVOLgEnON1i9JckPJ+nh548m+bYkNaFuT/qCqjqU5FCS7Nu37zybAbA7TApQG4MWMH/ndVdgd3+qu5/s7j9L8nP50nDfiSRXrKt6eZKHNvmOW7r7QHcf2Lt37/k0AwBgoZxXsKqqy9a9/eYkT90xeEeSG6vqmVX1giRXJvng9poIALAzbDkUWFW/kuQlSS6tqhNJfiDJS6rqqqwN8x1P8h1J0t0PVNXtST6W5Ikkr+/uJ2fTdACAxbJlsOrub51QfOvT1H9Tkjdtp1EAADvRdu4KBGCO3CkIi0ewAtih3CkIi8ezAgEARiJYAQCMRLACABiJYAUAMBLBCgBgJIIVAMBIBCsAgJEIVgAAIxGsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkghUAwEgEKwCAkQhWAAAjEawAAEYiWAEAjESwAgAYiWAFADASwQoAYCSCFQDASC6edwNgmR08cjQnT53est7qyp4L0Bp2g9WVPdl/+M6zyu4+fO2cWgS7i2AFM3Ty1OkcP3L9vJvBLjIpQG0MWsDsGAoEABiJYAUAMJIthwKr6q1JvinJI939NUPZc5O8I8n+JMeTfEt3P1ZVleQnklyX5AtJXtvdH5pN02GxTJpPZe4UwO4yzRyrX0jy00netq7scJL3dveRqjo8vP++JN+Y5Mrh9XVJ3jL8hKVnPhWLatKE9s3qmeQO27NlsOru36qq/RuKb0jykmH7tiTvz1qwuiHJ27q7k3ygqlaq6rLufnisBgNwbqYNSya5w/ad7xyr5z8VloafXzGUryb55Lp6J4YyAIClN/bk9ZpQ1hMrVh2qqnur6t5HH3105GYAAFx45xusPlVVlyXJ8PORofxEkivW1bs8yUOTvqC7b+nuA919YO/evefZDACAxXG+weqOJDcN2zclede68tfUmmuSPG5+FQCwW0yz3MKvZG2i+qVVdSLJDyQ5kuT2qro5yR8nedVQ/a6sLbVwLGvLLbxuBm0GAFhI09wV+K2b7HrphLqd5PXbbRQAwE5k5XUAgJEIVgAAIxGsAABGIlgBAIxEsAIAGIlgBQAwEsEKAGAkghUAwEi2XCAUONvBI0dz8tTpM8pWV/bMqTUALArBCrawWYg6fuT6ObUIgEUlWMEWTp46LUQBMBVzrAAARiJYAQCMRLACABiJYAUAMBLBCgBgJIIVAMBIBCsAgJEIVgAAI7FAKABJ1p4osP/wnWeV3X342jm1CHYewQqAJJkYoDYGLeDpCVYAnJNJz89MXN2CRLCCM2z2wGXgSzZ7fqarWyBYwRk8cBmA7XBXIADASAQrAICRCFYAACMRrAAARrKtyetVdTzJ55I8meSJ7j5QVc9N8o4k+5McT/It3f3Y9poJALD4xrhi9fe7+6ruPjC8P5zkvd19ZZL3Du8BAJbeLIYCb0hy27B9W5JXzOB3AAAsnO2uY9VJ/ldVdZL/1N23JHl+dz+cJN39cFV9xXYbCcB8bPb8QGCy7Qarg9390BCe3lNVvz/tB6vqUJJDSbJv375tNgOAWfCIGjg32xoK7O6Hhp+PJPn1JFcn+VRVXZYkw89HNvnsLd19oLsP7N27dzvNAABYCOcdrKrqL1TVc57aTvIPk3w0yR1Jbhqq3ZTkXdttJADATrCdocDnJ/n1qnrqe365u/9HVf1Oktur6uYkf5zkVdtvJgDA4jvvYNXdn0jytRPK/yTJS7fTKACAnWi7k9cBIMnmdxCaAM9uIlixKx08cjQnT50+q9xt5HD+JgWojUELlp1gxdKbFKJWV/bk+JHr59Qi2D1cxWK3EaxYeidPnRaiYE5cxWK3mcUjbQAAdiXBCgBgJIIVAMBIBCsAgJEIVgAAIxGsAABGIlgBAIxEsAIAGIlgBQAwEiuvA7CQNnsclcfhsMgEKwAW0qTHUXkcDotOsGKpbPYvXAC4EAQrlooHLgMwTyavAwCMxBUrdoRJQ3yTGPaDxbe6suesuVImpbMsBCt2BEN8sDwmBaiDR45ODFuw0whWAMydq1UsC8GKubJODXAuDCOy6AQr5so6NcC5mBSgnDNYJO4KBAAYiStWLJzNLvUDTDLpnHEunzWMyJgEKxaOkxxwLrZzzjCMyNgEKy4Yj5sBYNkJVszEZiHKWlTAInGXIWObWbCqqpcl+YkkFyX5+e4+MqvfxeKxoCewE7jLkLHNJFhV1UVJfibJP0hyIsnvVNUd3f2xWfw+ZsMaUwBrnA+Z1qyuWF2d5Fh3fyJJqurtSW5IIlgtgHN57t7Gq06THjux2WcBlsWkq/CbPYZns0f2CGa7w6yC1WqST657fyLJ183ody2Nsf/ibRagtjPXyUkAYM20zzxMJp93px1yPJd/DDtHz1919/hfWvWqJN/Q3d8+vH91kqu7+1+uq3MoyaHh7d9I8vHRG3K2S5N8+gL8nkWk77vXbu7/bu57srv7r++714Xo/1/t7r2TdszqitWJJFese395kofWV+juW5LcMqPfP1FV3dvdBy7k71wU+r47+57s7v7v5r4nu7v/+r47+57Mv/+zeqTN7yS5sqpeUFXPSHJjkjtm9LsAABbCTK5YdfcTVfWdSf5n1pZbeGt3PzCL3wUAsChmto5Vd9+V5K5Zff95uqBDjwtG33ev3dz/3dz3ZHf3X993r7n2fyaT1wEAdqNZzbECANh1li5YVdXLqurjVXWsqg5P2P/MqnrHsP+eqtp/4Vs5O1P0/7VV9WhV3T+8vn0e7ZyFqnprVT1SVR/dZH9V1U8O/20+UlUvvtBtnJUp+v6Sqnp83XH/txe6jbNSVVdU1fuq6sGqeqCq3jChzlIe+yn7vszH/llV9cGq+vDQ/x+aUGcpz/lT9n1pz/fJ2lNequp3q+rdE/bN77h399K8sjZR/g+T/LUkz0jy4SQv3FDnXyT52WH7xiTvmHe7L3D/X5vkp+fd1hn1/+8leXGSj26y/7okv5GkklyT5J55t/kC9v0lSd4973bOqO+XJXnxsP2cJP9nwp/7pTz2U/Z9mY99JXn2sH1JknuSXLOhzlKe86fs+9Ke74f+fXeSX57053uex33Zrlj9+aN0uvv/JXnqUTrr3ZDktmH7nUleWlV1Ads4S9P0f2l1928l+czTVLkhydt6zQeSrFTVZRemdbM1Rd+XVnc/3N0fGrY/l+TBrD39Yb2lPPZT9n1pDcfz88PbS4bXxonDS3nOn7LvS6uqLk9yfZKf36TK3I77sgWrSY/S2XiS+fM63f1EkseTPO+CtG72pul/kvyjYTjknVV1xYT9y2ra/z7L6u8Owwa/UVVfPe/GzMJwuf9FWfvX+3pLf+yfpu/JEh/7YTjo/iSPJHlPd2967JftnD9F35PlPd//eJLvTfJnm+yf23FftmA1KY1uTPDT1Nmppunbf0+yv7v/VpLfzJcS/W6wzMd+Kx/K2iMYvjbJTyX5b3Nuz+iq6tlJfjXJd3X3ZzfunvCRpTn2W/R9qY99dz/Z3Vdl7QkfV1fV12yosrTHfoq+L+X5vqq+Kckj3X3f01WbUHZBjvuyBastH6Wzvk5VXZzky7M8QyjTPEroT7r7i8Pbn0vyty9Q2xbBNH8+llJ3f/apYYNeW2Pukqq6dM7NGk1VXZK1YPFL3f1rE6os7bHfqu/Lfuyf0t2nkrw/ycs27Frmc36Szfu+xOf7g0leXlXHszbl5dqq+i8b6sztuC9bsJrmUTp3JLlp2H5lkqM9zG5bAlv2f8O8kpdnbU7GbnFHktcMd4hdk+Tx7n543o26EKrqLz81v6Cqrs7a3/0/mW+rxjH069YkD3b3j21SbSmP/TR9X/Jjv7eqVobtPUm+Psnvb6i2lOf8afq+rOf77n5jd1/e3fuz9v+5o939TzdUm9txn9nK6/PQmzxKp6r+XZJ7u/uOrJ2EfrGqjmUtvd44vxaPa8r+/6uqenmSJ7LW/9fOrcEjq6pfydodUJdW1YkkP5C1CZ3p7p/N2pMArktyLMkXkrxuPi0d3xR9f2WSf15VTyQ5neTGZfify+Bgklcn+b1hvkmSfH+SfcnSH/tp+r7Mx/6yJLdV1UVZC4y3d/e7d8k5f5q+L+35fpJFOe5WXgcAGMmyDQUCAMyNYAUAMBLBCgBgJIIVAMBIBCsAgJEIVgAAIxGsAABGIlgBAIzk/wNpMQl8ORcciwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xmin, xmax = 0.0, 4.0\n", "\n", "fig2, ax2 = plt.subplots(figsize=(10, 6));\n", "counts, bin_edges, _ = ax2.hist(y_all, 100, range=(xmin, xmax), histtype='step');\n", "bin_centers = (bin_edges[1:] + bin_edges[:-1])/2\n", "s_counts = np.sqrt(counts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the distribution of \"whatever you put into it\" (initially x1-2*x2), which shows what output you get and what uncertainty to expect (given by the width - think about this!). We can thus get the result by simply recording the mean and width (SD):" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Mean = 1.895, SD = 0.402\n" ] } ], "source": [ "print(f\" Mean = {y_all.mean():5.3f}, SD = {y_all.std(ddof=1):5.3f}\")" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "if save_plots:\n", " fig2.savefig(\"Dist_ErrorProp.pdf\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---------------------------------------------------------------------------------- \n", "\n", "# Questions:\n", "\n", "0. First solve the problem of obtaining the Perimeter, Area & Diagonal with uncertainty ANALYTICALLY.\n", "\n", "1. Now look at the program, and assure yourself that you understand what is going on. Put in the correct expression for y in terms of x1=L and x2=W in order to calculate the circumference, area, and diagonal length, and run the program. Does the output correspond well with the results you expected from your analytical calculations to begin with?\n", "\n", "2. Imagine that you wanted to know the central value and uncertainty of y1 and y2, given the\n", " same above PDFs for `x1`=$L$ and `x2`=$W$:\n", " \n", " `y1 = log(square(x1*tan(x2))+sqrt((x1-x2)/(cos(x2)+1.0+x1)))`\n", " \n", " `y2 = 1.1+sin(20*x1)`\n", "\n", " Get the central value of y, and see if you can quickly differentiate this with\n", " respect to `x1` and `x2`, and thus predict what uncertainty to expect for y using\n", " the error propagation formula. It is (for once) OK to give up on the first expression :-)\n", " Next, try to estimate the central value and uncertainty using random numbers\n", " like above - do you trust this result more? And are the distributions Gaussian?\n", "\n", "\n", "### Advanced questions:\n", "\n", "3. Try to generate `x1` and `x2` with e.g. a quadratic correlation, and see that despite\n", " not having any linear correlation, the result on perimeter, area, and diagonal\n", " length is still affected.\n", "\n", "---------------------------------------------------------------------------------- " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "executable": "/usr/bin/env python", "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }