content/error-propagation/investigating-propagation-of-errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating propagation of errors\n",
    "\n",
    "This notebook was made to investigate the propagation of errors formula.\n",
    "We imagine that we have a function $q(x,y)$ and we want to propagate the\n",
    "uncertainty on $x$ and $y$ (denoted $\\sigma_x$ and $\\sigma_y$, respectively) through to the quantity $q$.\n",
    "\n",
    "The most straight forward way to do this is just randomly sample $x$ and $y$, evaluate $q$ and look at it's distribution. This is really the definition of what we mean by propagation of uncertianty. It's very easy to do with some simply python code.\n",
    "\n",
    "The calculus formula for the propagation of errors is really an approximation. This is the formula for a general $q(x,y)$\n",
    "\\begin{equation}\n",
    "\\sigma_q^2 = \\left( \\frac{\\partial q}{\\partial x} \\sigma_x \\right)^2 + \\left( \\frac{\\partial q}{\\partial y}\\sigma_y \\right)^2\n",
    "\\end{equation}\n",
    "\n",
    "In the special case of addition  $q(x,y) = x\\pm y$ we have $\\sigma_q^2 = \\sigma_x^2 + \\sigma_y^2$.\n",
    "\n",
    "In the special case of multiplication $q(x,y) = x y$ and division $q(x,y) = x / y$ we have $(\\sigma_q/q)^2 = (\\sigma_x/x)^2 + (\\sigma_y/y)^2$, which we can rewrite as $\\sigma_q = (x/y) \\sqrt{(\\sigma_x/x)^2 + (\\sigma_y/y)^2}$\n",
    "\n",
    "Let's try out these formulas and compare the direct approach of making the distribution to the prediction from these formulas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%pylab inline --no-import-all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import widgets  \n",
    "from ipywidgets import interact, interactive, fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup repeated observations of two variables x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_x = .8\n",
    "std_x = .15\n",
    "mean_y = 3.\n",
    "std_y = .9\n",
    "N = 1000\n",
    "x = np.random.normal(mean_x, std_x, N)\n",
    "y = np.random.normal(mean_y, std_y, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check propagation of errors under addition "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'p(x)')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VNXWx/HvSiiBFLq0BBIQEaKgELAg6lXkRUXACwgqCEq5iCgWioqiInYUpIkKKNJUQIqKDSuChdCUKj2EIhh6SIAk+/3jTELAQEKYkz0zWZ/nyZMzM2cmv0GTNfucs9cWYwxKKaUUQJDtAEoppXyHFgWllFJZtCgopZTKokVBKaVUFi0KSimlsmhRUEoplUWLglJKqSxaFJRSSmXRoqCUUipLEdsBzlX58uVNdHS07RhKKeVXli5d+o8xpkJu+/ldUYiOjiY+Pt52DKWU8isisi0v++nhI6WUUlm0KCillMqiRUEppVQWvzunoJRS+XHixAkSExNJTU21HcVVISEhREZGUrRo0Xw9X4uCUqpQSExMJDw8nOjoaETEdhxXGGNISkoiMTGRmJiYfL2GHj5SShUKqamplCtXLmALAoCIUK5cufMaDWlRUEoVGoFcEDKd73vUoqCUUiqLFgWllPJTYWFhXn9NPdGsVD4sWPN3nvdtVreii0lUoElPTyc4ONjaz9eRglJKFZCtW7dy8cUX06VLF+rVq0e7du04evQo0dHRDBkyhGuuuYYZM2awadMmWrRoQcOGDWnatCnr1q0DYMuWLVx11VU0atSIp59+2pWMWhSUUoWPiDtfebB+/Xp69uzJH3/8QUREBGPHjgWc+QU///wzHTt2pGfPnowaNYqlS5cybNgwevfuDUDfvn25//77WbJkCZUqVXLln0YPHymlVAGKioqiSZMmAHTq1ImRI0cC0KFDBwCOHDnC4sWLad++fdZzjh07BsCiRYuYNWsWAJ07d2bgwIFez6dFQSlV+Bhj7Ueffslo5u3Q0FAAMjIyKF26NCtWrMjT871NDx8ppVQBSkhI4JdffgFg+vTpXHPNNac8HhERQUxMDDNmzACcWcorV64EoEmTJnz44YcATJ061ZV8WhSUUqoA1alTh0mTJlGvXj327dvH/fff/699pk6dyoQJE6hfvz6xsbHMnTsXgDfffJMxY8bQqFEjDh486Eo+PXyklFIFKCgoiHHjxp1y39atW0+5HRMTw5dffvmv58bExGSNMgAef/xx7+fz+isqpZTyWzpSUMpleZ3oppPcAl90dDSrVq2yHeOsdKSglFIqixYFpZRSWbQoKKWUyqJFQSmlVBY90axUPgSlpuRtVmxQEBnFQ9wPpM7ZuXS6zYu8XCiwdetWWrZs+a+TzYMHD+baa6+lWbNmOT5vzpw5XHTRRdStW9crWc9Gi4JS52LPHujShRtyuIb8jE9pdgtrnnudtNJlXAym/NmQIUPO+vicOXNo2bJlgRQFPXykVF4tXAiXXQZffklGcDDpISVy/coIDuaCBfO5ov1NRPy53PY7UD4gPT2dHj16EBsbS/PmzUlJSaFr167MnDkTcCak1a1bl3r16tGvXz8WL17MvHnz6N+/P5dddhmbNm1yNZ+OFJTKjTHw2mvw5JOQng5Nm/LzkFEcvyD31sUhO7dz6SM9KLVqBXGdW/PXgGdJvPPePLdZVoFnw4YNTJ8+nXfffZc77rgjq+spwL59+5g9ezbr1q1DRDhw4AClS5emVatWtGzZknbt2rmeT4uCCnjntUra/v3QtSvMm+fcHjgQhg7l+F9JeXq91CpRxE+ey0WvPUfUtIlc/MKTlF72O2ufG0Z6qPeXUlS+LyYmhssuuwyAhg0bntLiIiIigpCQELp3786tt95Ky5YtCzyfHj5S6kyWLoUGDZyCULq08/3ll6HIuX2WMsWKs37Qi/w5bBxpJUOp9MUcGndoQeiGtS4FV76sePHiWdvBwcGkpaVl3S5SpAi///47bdu2Zc6cObRo0aLA82lRUOp0xsBbb8HVV8PWrRAXB8uWwW23ndfL/n1zG37/+CuO1LqY0C0badzxFirP+cg7mVVAOHLkCAcPHuSWW25hxIgRWWsqhIeHc/jw4QLJoIePlMruyBHo2ROmT3du9+4Nb7wB2T7dnY+jMRfy+/T5XDz0CarM+YjYQX0pvew31j/5gldeX+WdL/aaOnz4MK1btyY1NRVjDMOHDwegY8eO9OjRg5EjRzJz5kxq1qzpWgYxFlcgyo+4uDgTHx9vO4byI3k9pxC6cR1XDewF69ZBaCiMHw8dO57Xa55NlU+mUXvokwQfS+Vw7VjCP50NtWqd9+uqnK1du5Y6derYjlEgcnqvIrLUGBOX23P18JFSQPjaP2nc8WanIMTGQnz8GQuCt+z8710smfY5R6vFEL5+NTRq5Px8pSzSoqAKvaDUFGIH9CY4JQXatYPffoOLLy6Qn33k4lh+m/E1e69rBgcPQqdOcPx4gfxspXKiRUEVehcOf4GwzRtIrlELJk1yDh0VoPSwcFa9MhaqV3eueMpldqvKP387XJ4f5/setSioQq3soh+oNmU8GUWKsOqVMVCypJUc6eERMHmyM6ntpZdg0SIrOQJZSEgISUlJAV0YjDEkJSUREpL/flt69ZEqtIoe2EfsoL4AbO4zgMN169kN1LSpMznu5Zehc2dYsQIiIuxmCiCRkZEkJiayd+9e21FcFRISQmRkZL6f72pREJEWwJtAMDDeGPPyGfZrB8wAGhlj9NIi5T5juPi5ARTf+zf7G1zB1vseALzfOfOcPfccfPUVLF8ODz8MEyfazRNAihYtSkxMjO0YPs+1w0ciEgyMAW4G6gJ3isi/WvyJSDjwEPCbW1mUOl3leTOo+PVnpIWGsfqlkRAcbDuSo1gxmDIFQkLgvffgk09sJ1KFjJvnFBoDG40xm40xx4EPgdY57Pc88CqQ6mIWpbKE7Eig9gtPArD+iaGkRla3nOg0devCK6842z17wq5ddvOoQsXNolAV2J7tdqLnviwicjkQZYz57GwvJCI9RSReROID/Xigcll6OrFPPESR5CPsaXYru9p0sJ0oZ336wE03QVISdOuWtwV9lPICN4tCTr2Bs/7PFpEgYDjwWG4vZIx5xxgTZ4yJq1ChghcjqsKm+ntjKbP0V46Vv4C1z77quy2sg4Kcw0dlysAXXzi9mJQqAG4WhUQgKtvtSGBnttvhwCXADyKyFbgSmCciuU7DVio/wtf+Sc1RrwKwZugITpQpZzlRLqpWhbffdrb79YP16+3mUYWCm0VhCVBLRGJEpBjQEZiX+aAx5qAxprwxJtoYEw38CrTSq4+UGzJnLQelnWD7nfeS1PQG25Hypn175/LUlBRntvOJE7YTqQDnWlEwxqQBfYCvgLXAx8aY1SIyRERaufVzlcpJ9lnLGx572nacczNqlDPbOT5eZzsr12mXVBXwlo3/mAY9OpBRpAhLpn3O4dj6tiPl6KytnH/6Ca6/3jkHsnChs9aDUucgr11SdUaz8lt5mWhW5MB+rsyctdy7n88WhFxdey307w+vvnpytnN4uO1UKgBp7yMV0Gq//DQhe3Zz4LJGbOvWx3ac8zNkCNSvD5s3w+OP206jApQWBRWwwtatpvKnM8koWozVL43CnOPayj6neHFntnNQkHNV0saNthOpAKRFQQWsmqOcWcGJHbqQUi3abhhvueQS6NIF0tPh2Wdtp1EBSIuCCkgRfyyjwg9fk16iBFt7PGQ7jncNHgxFi8K0abB6te00KsD4+XhaqZzVHOk05E24uzvHy/vHLPi8d2gtQbMePWDsWKdAzJrlai5VuOhIQQWcMr8votwvP5EWFs62e3vbjuOOQYOcTqqffOKs1qaUl2hRUIHFGGq+6YwStnW9n7TSZSwHckmVKk7TPICnnrKbRQUULQoqoJRb+C2lVyzheJmyJNzT03Ycdw0c6MxV+PJL+Pln22lUgNCioAJHRgY1RzpXHG3t/iDpoWGWA7msfHl45BFne9Agba+tvEKLggoYFyyYT8TaPzlWoSKJHbvajlMwHn3Uaa/900+wYIHtNCoAaFFQgSE9nRqeeQlb/vcwGSElLAcqIKVKwYABzraOFpQX6CWpKiBU+vwTwjZvIKVqFDva3m07juuyX74a1OwOmgx7g+JLlrByzAfsvaFF1mNnbbKnVA50pKD8npw4QY2xwwDYfP9jmGLFLCcqWBklQ9na02n6V2PUK5CRYTmR8mdaFJTfqzJ7OiW3byM55kJ239bOdhwrEu/oTGqlqoT/tZaKX8y1HUf5MS0Kyq8FHUslZtxwADb16e//Te/yyRQrzub7HwWgxpjXkLQ0y4mUv9KioPxa5IeTCPl7F4drx7Kn+W2241i1q/UdHI2KJnTbZirP/dh2HOWntCgovxWcnEz0+JEAbHpwoNNSuhAzRYuyuU9/AGLeeh05fsxyIuWPCvdvkfJrUVPHU2xfEgfrNeCf62+yHccn7L65DUcurE2JXTuoOmOK7TjKD2lRUP7pwAGqvzcWgI19n3DWLlYQHOyMmoCYd96Eo0ctB1L+RouC8k/DhlH00EH2NW7C/iub2k7jU/beeDMHY+tT/J89MHq07TjKz2hRUP5n714YMQKATQ/pWsX/InLy3+WVV+DQIbt5lF/RoqD8z5gxkJzMP01v4ODljWyn8Un7mlzPgcsbw759MH687TjKj2hRUP7l6FGnKOB0QlVnIMLW7p71FkaMgBMn7OZRfkOLgvIvH3wA//wDcXEcaHil7TQ+7Z9rm0Ht2rB9O8yYYTuO8hNaFJT/yMiAN95wtvv10yuOchMUBI895my//rp2UFV5okVB+Y9PP4UNG6B6dWjb1nYa/9C5M1SoAMuWwQ8/2E6j/IAWBeU/hjmdUHn4YSikPY7OWUjIybWcM//9lDoLLQrKP/z2m7MOcalS0K2b7TT+pXdvpzjMnw9r1thOo3ycFgXlH15/3fn+v/85i9WrvCtfHrp2dbYzz8kodQZaFJTv27IFZs1yDhk99JDtNP7pkUecE/OTJ8Pu3bbTKB+mRUH5vhEjnCuP7roLqla1ncY/XXQRtG4Nx49nzfNQKidaFJRv278fJkxwth991G4Wf5d5eerYsZCcbDeL8llaFJRve/tt5w/YTTdB/fq20/i3Jk3giiuc1heTJtlOo3yUq0VBRFqIyHoR2Sgi/+pcJiK9RORPEVkhIj+LSF038yg/c+wYjHQW0aFfP7tZAoHIyX/HN96A9HS7eZRPcq0oiEgwMAa4GagL3JnDH/1pxphLjTGXAa8CemmEOmn6dNi1Cy691BkpqPN3++0QEwObNsHcubbTKB/k5kihMbDRGLPZGHMc+BBonX0HY0z2nr6hgM7DVw5jTl6G+thj2tLCW4KDnSuR4OS/r1LZuFkUqgLbs91O9Nx3ChF5QEQ24YwU9HpD5fj6a1i1CqpUgTvvtJ0msNx7L5QpA4sXO19KZeNmUcjpo92/RgLGmDHGmJrAQOCpHF9IpKeIxItI/N69e70cU/mkzJYMDz4IxYrZzRJowsKgVy9nW0cL6jRuFoVEICrb7Uhg51n2/xBok9MDxph3jDFxxpi4ChUqeDGi8kkrV8KCBRAa6sxgVt734INQtCjMnu2cX1DKw82isASoJSIxIlIM6AjMy76DiNTKdvNWYIOLeZS/yPz02r27c5hDeV/lytCpk3PuZvhw22mUD3GtKBhj0oA+wFfAWuBjY8xqERkiIq08u/URkdUisgJ4FOjiVh7lJxITnauOgoKgb1/baQJb5mTAiRMhKcluFuUzXO0/bIyZD8w/7b7B2bb1t16dYuszLxOdlsbuFq1YlVIS1vxtO1LguuQSaNECvvwSxo2DQYNsJ1I+QGc0K99x6BCRH38AQELX+y2HKSQyJ7ONGgWpqXazKJ+gRUH5jgkTKHLkMPsbXsmhSy+3naZwuOEGp33I33/D1Km20ygfoEVB+Ya0tKyWFglde1kOU4hkb30xfLiu46y0KCgfMW8ebN3K0aho9l7f3HaawuWOO5yrkVavhm+/tZ1GWaZFQfmGESMA2N65u3PlkSo4xYrBAw84257/Dqrw0t8+Zd/SpbBwIUREsLNNR9tpCqf//c9Zx/nzz2H9ettplEVaFJR9mZ9Oe/QgPTTMbpbCqnx56NzZ2c5sV64KJS0Kyq6dO+Gjj5xDRn362E5TuGVOFnz/fWfFO1UouTp5TalcvfUWnDgBbdtCdLROVrNogZTn8quvo9ziH9nwwnC23fdAjvs1q1uxgJOpgqQjBWVPSoozkxbg4YftZlEAJHTuAUDU1IlIWprlNMqGcyoKIhLqWVFNqfM3bRr88w80bOisH6ysS7rmBpKjaxKyewcVvp2f+xNUwDlrURCRIBG5S0Q+F5E9wDpgl6eJ3WundTlVKu+MOXmC+eGHdWU1XxEUxPZO3QGo9sG7lsMoG3IbKXwP1ASeACoZY6KMMRcATYFfgZdFpJPLGVUg+u47Z2W1SpWcyVPKZ+xqdQcnIkpResUSIv5YZjuOKmC5FYVmxpjnjTF/GGMyMu80xuwzxswyxrQFPnI3ogpImaOEBx7QldV8THpoKDva3g1A1JTxltOognbWomCMOQEgIs1Of0xEumTfR6k827ABPvsMihfXldV81Pa77iMjOJiKX82j+N+7bMdRBSivl6QOFpG2QD8gDBgPHAMmuRVMBbA333S+d+oEuryqqxbk8xLfY1Ui2dvsFip+9SmR099j08NPejmZ8lV5vfroOmATsAL4GZhmjGnnWioVuPbvh/fec7b1MlSflnCPM4qr+vFkglKOWk6jCkpeRwplgCtwCkMkUF1ExBjts6vyJvMTa7X3xnLR0aMkXXUty4Mq6GQ1H3awfkMOXno5pf5cTuVPZ7Hjjs62I6kCkNeRwq/AF8aYFkAjoAqwyLVUKiBJWhpRUycAJydJKR8mQkLnngBETX5H11ooJPJaFJoZYyYCGGNSjDEPAY+7F0sFogrffkGJXTtIrl6DpKY32o6j8mBP85akXlCJsM0bKLv4R9txVAHIbfJaNIAxJuH0x4wxP4kj0p1oKtBUm+xMhtreSddM8BemaFES77oPgGqT37GcRhWE3H4zXxORWSJyj4jEisgFIlJNRG4QkedxDiHVKYCcys9F/Lmc0st/50R4BLtad7AdR52DHe07kV48hPILv6Pk5g224yiX5TZPoT3wNFAbGAP8BMwFugPrgRuMMd+4HVL5v8xJUDva3U16aKjlNOpcnChdll2t2gM6ma0wyHUMb4xZAwwFPgXWAluAJcBMY0yqu/FUQNi5k4pfzsUEBWUdilD+JbMfUpV5H8O+fZbTKDfl9cDuJJzDRCOBUZ7tD9wKpQLM2LEEpaWxp9ktpFaJsp1G5UPyhbVJuvp6glNSYLyOFgJZXotCbWNMd2PM956vnjiHlJQ6u2xrJmROhlL+KeEe5/JURo1yFkZSASmvRWG5iFyZeUNErkDnKai8mDQJkpI4GFufg5fF2U6jzkNSk+tJrlELEhNhxgzbcZRL8loUrgAWi8hWEdkK/AJcJyJ/isgfrqVT/i09HYYNA2Dbfb11zQR/FxTEti69nO1XX9XJbAEqr20uWriaQgWm2bNh0yaoUYO9zW61nUZ5we7b2lL3rddg5Ur45hto3tx2JOVleRopGGO2ne3L7ZDKDxnjfJoEeOwxTJG8fv5QviyjeAj07evcyPzvqwKKTitV7vjxR1iyBMqXh65dbadR3tSrF4SHw7ffwtKlttMoL9OioNyR+SnyoYegZEm7WZR3lS59cnEkHS0EHC0Kyvv++AO++MIpBr17206j3NC3LxQtCjNnOueNVMDQoqC877XXnO/du0O5cnazKHdERsLdd0NGBrzxhu00yovE39bJiYuLM/Hx8bZjKHJe6jFk53aubuFMaVn8xS+kVq1W0LGUy5rVrehsrF4Nl1wCISGQkKBLq/o4EVlqjMl1spCrIwURaSEi60Vko4j8a/0FEXlURNaIyB8i8q2IVHczj3JftQ/eJSg9nb9btNaCEOhiY6FlS0hNhdGjbadRXuJaURCRYJzOqjcDdYE7RaTuabstB+KMMfWAmYCetfJjRQ7sp8rMKQBsu1fPJRQKAwY430ePhuRku1mUV7g5UmgMbDTGbDbGHAc+BFpn38HTRylzRfBfcdZ/Vn4q8qNJFEk5StLV13OkziW246iCcM01cOWVTufUCRNsp1Fe4GZRqApsz3Y70XPfmXQDvnAxj3JRUGoKUVOd7plbu/WxnEYVGBEYONDZfv11bZQXANwsCjk1usnxrLaIdALigNfO8HhPEYkXkfi9e/d6MaLylsrzZlA86R8OxdZj/xVNbMdRBalVK7joIudkszbK83tuFoVEIHvz/Ehg5+k7iUgzYBDQyhhzLKcXMsa8Y4yJM8bEVdArHHxPejrV33sLgK33PaCN7wqboCDo39/Z1kZ5fs/NorAEqCUiMSJSDOgIzMu+g4hcDryNUxD2uJhFueiCBfMpmbCFo1HVtfFdYdWpE1SqdLJRnvJbrhUFY0wa0Af4CmcZz4+NMatFZIiItPLs9hoQBswQkRUiMu8ML6d8lTFUnzgGgIQuvbTxXWEVoo3yAoWrv8HGmPnA/NPuG5xtu5mbP1+5r8ySxZRatYLjZcqys00H23GUTb16wQsvnGyU17Ch7UQqH7TNhTovmaOE7Xd1I6OENr4r1LRRXkDQoqDyLWz9Gsov/I70EiVIvOte23GUL3j4YW2U5+e0KKh8q/7+WAB2tL2bE6XLWk6jfII2yvN7WhRU/iQkUHH+HDKCg0m4p6ftNMqX9OvnfJ84EXRekd/RS0VU/gwfTlBaGrtu/a82vlOnymyU99lnbH76JTY/NDDXp2R1XlXW6UhBnbsdO2DcOAC23feA5TDKJ3laX1Sb8i5FD+yzHEadCy0K6twNHQqpqfzdvCVHLo61nUb5omuugebNKZJ8hOrjta22P9GioM7N5s0wfjwEBbGpzwDbaZQve+EFAKKmv0exPbsth1F5pUVBnZtnn4W0NLjnHo7WvMh2GuXL4uLY0+xWglNTiHl7hO00Ko+0KKi8W70apkxxrkN/5hnbaZQf2PRgf4wIVWdOISRxm+04Kg+0KKi8GzzY6YDZsydER9tOo/xA8oUXs/u2tgSlpVFj7Ou246g8EONnbW7j4uJMfHy87RiFz5Il0LgxlCjhzFStXJkFa/62nUr5gRLbt3FVyyZIRga/zPkhx8OOekmq+0RkqTEmLrf9dKSg8uapp5zvDz4IlSvbzaL8SkpUdXa2vQvJyKDmaO2J5Ou0KKjc/fgjfP01REScXKhdqXOw5X+PkF48hIpff0b46pW246iz0KKgzs4YGDTI2X7sMShXzm4e5ZeOVaxM4p1dAag58mW7YdRZaZsL9S/ZzxWU+2kBly9axPEyZVl0892k63kElU9buz9I1Y8nU/7n7ym99FcONLzSdiSVAx0pqDPLyKDmyFcA5xc6PTTMciDlz06UKUdCF2e9hZpvvqxrOfsoLQrqjC745jMi1v5J6gWVSOzY1XYcFQC2denF8VJlKLP0V8ou+sF2HJUDLQoqR5KWRs1RzpUiW3o9QkZICcuJVCBID49gWzenieKFI1/S0YIP0qKgclTp05mEbtnI0ajq7Lz9TttxVADZftd9HCt/ARGr/6DCgvm5P0EVKC0K6l/k+DFqjB0GwOYH+mOKFbOcSAWSjBIl2dLrEcBzJVJ6uuVEKjstCupfqs6cSomdiRypeRG7b7nddhwVgHa0vZuUqlGEbd5A5c9m2Y6jstGioE6VnEzMuOEAbHrocQgOthxIBSJTrBibezvLdsaMHQbHj1tOpDJpUVCnGj2a4kl7ORhbn7033mw7jQpgu25rR3KNWpRMTIAJE2zHUR5aFNRJBw7AK868hE19nwARy4FUQAsOPrlQ0/PPQ3Ky3TwK0KKgsnviCdi/n/2NrmLf1dfZTqMKgT033cqh2HqwaxcMGWI7jkKLgsq0eDGMGwdFirBu0Is6SlAFIyiIdU+97Pz/9vrrsFKb5dmmRUE5J/l69nS2BwwguVYdu3lUoXKoXgN44AHn0tSePfUSVcu0KCgYNsxZavPCC0+um6BUQXrhBahaFX7/Hd56y3aaQk2LQmG3cePJY7njxjkrqylV0CIiYNQoZ/vJJ2HHDrt5CjEtCoWZMdCrFxw7Bp07w4032k6kCrPbb4fWreHwYWeFP2WFFoXCbMoU+PZbZ+Gc13VRdeUDRo2CsDCYPRvmzrWdplDSolBY/fMPPPqosz1sGFSoYDePUgBRUc75BYA+fZxRgypQWhQKq/79ncJw/fXQpYvtNEqd9MAD0KgRJCbqhQ8WiPGzfuZxcXEmPj7edgz/9v33cMMNUKwY/PEH1K59ysMLdMlNZVnY2lU07vB/SEYGS6bP59Cll2c91qxuRYvJ/JeILDXGxOW2n6sjBRFpISLrRWSjiDyew+PXisgyEUkTkXZuZlEeqanwP2dJRAYN+ldBUMoXHKlzCQn39ESMoc6z/ZG0NNuRCo0ibr2wiAQDY4CbgERgiYjMM8asybZbAtAV6OdWDuXI/PRfY9Qr1NiwgeQatfi1VVeMjgqUj9rcux8Vv/qU8HWriJr8Dgn39rYdqVBwc6TQGNhojNlsjDkOfAi0zr6DMWarMeYPIMPFHMojdON6osePBmDts69hihW3nEipM8soGcq6wU6DxppjXiNkR4LlRIWDm0WhKrA92+1Ez33KhowMLn5uAEFpJ9jR9m4ONLzSdiKlcpXU9EZ239ya4JQULn7+cV3TuQC4WRRy6qiWr/+iItJTROJFJH7v3r3nGatwqjJrGmWW/caxcuXZ8NjTtuMolWd/DXyeExGlKL/wOyp+qXMX3OZmUUgEorLdjgR25ueFjDHvGGPijDFxFfR6+nO3eze13ngecH7B0kqVthxIqbw7XuECNj7iXJp60UtPw/79lhMFNjeLwhKglojEiEgxoCMwz8Wfp3JiDNx/P0UPHeSfJv/h71va2E6k1Dnb0e5uDlzemOJJe+Ghh/QwkotcKwrGmDSgD/AVsBb42BizWkSGiEgrABFpJCKJQHvgbRGpbxPcAAAOI0lEQVRZ7VaeQuvFF2HOHNLCwlk3+GVdJ0H5p6Ag1jw3jPSQEk57ltGjbScKWDp5LZB9/jncdhsAy8d8QNJ1N1kOpNT5qfjFHC7t1wuCg2HBAmdGvsoTn5i8pixavx7uussZZj//vBYEFRD+vrkNDBjgLMTTvj0k6GWq3qZFIRAdOgRt2jjf27Z1+tMrFShefBGaN3d6d91+O6Sk2E4UULQoBJqMDLjnHli3DmJj4f339TyCCizBwTB9OtSoAcuWOUt4+tlhcF+mRSHQPP+804e+dGmYM8fpTa9UoClb1vn/OzTUOfE8YoTtRAFDi0IgmTcPnn3WGRlMn+6suaxUoLr0Upg0ydnu399ZMEqdNy0KgWLtWujUydl+6SVo0cJuHqUKQuY5s/R06NABtmyxncjvaVEIBAcPOieWDx+GO+5wrs5QqrAYMgRuuQWSkpwTz0eP2k7k17Qo+LuMDLj7bvjrL2c4PXGinlhWhUtwMEydCrVqwcqV0K2bnng+D1oU/N0zzziT1MqUOXniTanCJvuFFR9+6Kw7rvLFtUV2VAH45BMYOhQTFMTyV8exLzUUdNEcVVjVrQuTJzuHkB5/HOrXd+YzqHOiIwV/9fPPznwEYMOjT7Hv6ussB1LKB7RpA4MHO4dVO3QAbYlzznSk4IMW5PJpv9yP31DvkR4EH0tl123tSOh6fwElU8oPPPMM/PknzJ4N//mPc6n2f/5jO5Xf0JGCn6n06UzqP9iV4GOpJLbrxOoX3tQTy0plFxQEH33k9P46csS5PHv2bNup/IaOFPxI1JTx1H7JWWxkS/cH2fTwk1oQVKGT20g6U7PJk52Zz6NHQ7t28O67cN99BfKzAZrVrXheP8sWLQr+wBhqjH6VGuOGA/BX/2f0kJFSuQkKgpEjoXx5Z6Z/t27OXIb+/W0n82l6+MjXZWRQe+gT1Bg3HBMUxOqhw7UgKJVXIs45hlGjnNsDBsDAgTqP4Sx0pODD5PhxYp98iEpfzCG9WHFWDRvH3htvth1LKb9wyqGeG9pT8ZVgYgc9RNCrr7JjUyLrBr+KKeL8CfTXQz1u0JGCjwo6mkz9Pl2o9MUc0kLDWPH2NC0ISp2Hv1v+l5Wj3ic9pARVZ03j0sd6EnQs1XYsn6NFwQcVObCfBj06UH7R9xwvW46l73/C/sZNbMdSyu8lXduMZe9+xImIUlywYD6X3d+J4OQjtmP5FC0KvmbHDuK63k7pFfGkVK5K/OR5HK5bz3YqpQLGwQaNWfr+bI6Vv4Cyv/1Mg3vbwp49tmP5DC0KvsIYZ4p+vXqEbVjHkRq1iJ/yKUeja9pOplTAOVK7LvGT53E0qjqlVq90mknOmKEnoNGi4Bu2bXNa/95zD+zbR9LV17N08lyOVapiO5lSASulWjTxk+exP+5KZ6Rwxx3w3//Czp22o1mlRcGmjAwYMwYuuQS+/NLpdPr++yx/ZzonSpe1nU6pgHe8QkWWvvcJvPUWhIc7nVbr1oUJEwrtqEGLgi3r1sG110KfPs5U/HbtYM0a6NJFZykrVZCCgqBXL1i9Gm691Vm0qnt3aNYMNm+2na7AaVEoaCdOwIsvOm19Fy2CSpVg1izneGalSrbTKVV4RUXBp586C/aUKwfffeeM4ocPd5b7LCS0KBSkZcugUSMYNAiOH3f6sKxZ4xzHVEpZs2DN387X2j0suOxGfpzzI7tvuR1SUuDRRzl4eSN+mfuD7ZgFQotCQdizx+m30rixs1xgdDR8841z3LJMGdvplFKnOVG2PKtee4sVoz8gtWJlSv25nCvaNafGyFcoemCf7XiuEuNnJ1Pi4uJMvL8snLF8Obz5Jkyf7owMRKBvXxg69KzLZp5LJ0allLuCDx+i1htDifz4AwDSQ0qwq2VbtnfqRnKtOuf9+gXVYkNElhpj4nLbT0cK3paW5pwjuPZaaNAAJk1yziPcdhv88otzfFLXUVbKb6SHR7DumVeJ/2AO/zT5D8GpKUTOnMJVbf7D5d3aU/77r50rCQOENsQ7T5mf6oscPEDVWVOJnDaRErt2AJAWGsbO/97F9rvuI6VatDbdUsqPHWh4JSveuZKSmzcQNXUCVeZ+RLlfF1Lu14UcjYpm+93d2Hl7R9LDwm1HPS96+Og8Lf50IdWmjqfyvBkEp6QAcLRaDNs7dWNnm46kh4Zl7ZvXoqCHj5TyfUUOHaTKJ9OImjaREju2A54Pgrd3ZPtd3UipHpOn1/G1w0daFM7V8ePw66/w1Vfw9denLAyedPV1JHTqTlLTG51rn0+jRUGpAJSeToXvv6LalHcps+QXAIwIB+s1YN/V15F09fUcqtcgq0336bQonKcCLwrGwMaNTgH4+mvn2uUjJ7sqpoeUYFer9my/uxvJF9YuuFxKKZ8TtnYV1aaOp+Lnswk+fizr/rSwcPZd2ZSkJteTdPV1pEZWz3pMi8J5KpCicOAAfP/9ydHAli2nPl63Lvzf/0Hz5nxXoTYZJUq6m0cp5VeCk5MpHb+Ycot+oNziHwndsvGUx5Or18gaRVx2TxuIiHA9kxaFHJx+WCb4yGFCN28gbOM6QjeuJ3TTX4Rt/IuQ3TtOfWLZsnDTTdC8ufMVGXnG11RKqdOF7Eig7OKfKLfoe8r+upCihw+dukP16hAb63zVret8r1MHwsJyfsF88ImiICItgDeBYGC8Mebl0x4vDnwANASSgA7GmK1ne818F4Xly1nz+Q9n/+PvkVG0GEFXNHZGA//3f86lpcHBOe6rRUEpdS4kLY2IVSsou/gHyi36gdKrVzqXrefk9GLRoAHUy9/6KtaLgogEA38BNwGJwBLgTmPMmmz79AbqGWN6iUhH4HZjTIezvW6+i8J118FPP51yV0bRYiTHXEjyhbU5cmFtkmteRPKFtUmJrH7Gk0JKKeVNkpZGiYQthHk+sIZuWk/YxvWU3LKJoLTTikXLlk5/pvz8nDwWBTf/8jUGNhpjNnsCfQi0BtZk26c18KxneyYwWkTEuFGpbr6Z3WHl9I+/UsqnmCJFOFqjFkdr1Drl/pyKRaWbbnA9j5t/EasC27PdTgSuONM+xpg0ETkIlAP+8Xqaxx9nVSs91KOU8g85FYtKBXClkptFIadFAU4fAeRlH0SkJ9DTc/OIiKw/hxzlcaPI+IZAfm+g78+fBfJ7A/98f9Vz38XdopAIRGW7HQmcvs5d5j6JIlIEKAX8qwWhMeYd4J38hBCR+LwcR/NHgfzeQN+fPwvk9waB/f7cbIi3BKglIjEiUgzoCMw7bZ95QBfPdjvgO1fOJyillMoT10YKnnMEfYCvcC5JnWiMWS0iQ4B4Y8w8YAIwWUQ24owQOrqVRymlVO5cvfTGGDMfmH/afYOzbacC7d3MQD4PO/mJQH5voO/PnwXye4MAfn9+N6NZKaWUe3SRHaWUUlkCtiiIyEQR2SMiq2xn8TYRiRKR70VkrYisFpG+tjN5k4iEiMjvIrLS8/6es53J20QkWESWi8hntrN4m4hsFZE/RWSFiPjO4ideICKlRWSmiKzz/P5dZTuTtwXs4SMRuRY4AnxgjLnEdh5vEpHKQGVjzDIRCQeWAm2ytxDxZyIiQKgx5oiIFAV+BvoaY361HM1rRORRIA6IMMa0tJ3Hm0RkKxBnjPG36/hzJSKTgIXGmPGeqypLGmMO2M7lTQE7UjDG/EQOcx4CgTFmlzFmmWf7MLAWZ3Z4QDCOzEUrinq+AubTi4hEArcC421nUXknIhHAtThXTWKMOR5oBQECuCgUFiISDVwO/GY3iXd5Dq+sAPYA3xhjAun9jQAGAIGz2vupDPC1iCz1dCMIFDWAvcB7nkN/40Uk1HYob9Oi4MdEJAyYBTxsjDmU2/7+xBiTboy5DGcmfGMRCYhDgCLSEthjjFlqO4uLmhhjGgA3Aw94DuUGgiJAA+AtY8zlQDLwuN1I3qdFwU95jrXPAqYaYz6xncctnuH5D0ALy1G8pQnQynPc/UPgBhGZYjeSdxljdnq+7wFm43RMDgSJQGK2UetMnCIRULQo+CHPidgJwFpjzBu283ibiFQQkdKe7RJAM2Cd3VTeYYx5whgTaYyJxpnB/50xppPlWF4jIqGeix/wHFppDgTEFYDGmN3AdhHJXIz9Rk5dCiAgBOxiAiIyHbgeKC8iicAzxpgJdlN5TROgM/Cn57g7wJOeGeSBoDIwybNQUxDwsTEm4C7dDFAVgdnO5xaKANOMMV/ajeRVDwJTPVcebQbutZzH6wL2klSllFLnTg8fKaWUyqJFQSmlVBYtCkoppbJoUVBKKZVFi4JSSqksWhSUOk8iMiK3WbsiskBEyhRUJqXyS4uCUudBRMoCV3oaMJ7NZKB3AURS6rxoUVAqj0RkkIis93zqny4i/YB2wJeex0t5Hq/tuT1dRHp4nj4PuNNOcqXyTouCUnkgIg1x2lJcDvwXaOR5qAnOehYYYw4CfYD3RaQjUMYY867nsf1AcREpV9DZlToXAdvmQikvawrMNsYcBRCReZ77K+O0UwbAGPONiLQHxgD1T3uNPUAVIMn9uErlj44UlMq7nHrCpAAhmTdEJAio47m/7Gn7hnjuV8pnaVFQKm9+Am4XkRKeLqC3ee5fC1yYbb9HPPfdCUz0tDjP7GxbCdhaYImVygctCkrlgWf504+AFTjrWCz0PPQ5TjdeROQioDvwmDFmIU4hecqzX0PgV2NMWgHGVuqcaZdUpfJBRJ4FjhhjhonIz0DLs63XKyJvAvOMMd8WVEal8kNHCkqdv8eAarnss0oLgvIHOlJQSimVRUcKSimlsmhRUEoplUWLglJKqSxaFJRSSmXRoqCUUiqLFgWllFJZ/h/v5cOy+6qVZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.random.normal(mean_x, std_x, N)\n",
    "y = np.random.normal(mean_y, std_y, N)\n",
    "\n",
    "q_of_x_y = x+y\n",
    "\n",
    "pred_mean_q = mean_x+mean_y\n",
    "pred_std_q = np.sqrt(std_x**2+std_y**2)\n",
    "\n",
    "counts, bins, patches = plt.hist(q_of_x_y, \n",
    "                                 bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), \n",
    "                                 density=True, alpha=0.3)\n",
    "\n",
    "plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)\n",
    "plt.legend(('pred','hist'))\n",
    "plt.xlabel('q(x)')\n",
    "plt.ylabel('p(x)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### same thing with an interactive widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "700be3b4a566432aba48aa19c1f36233",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_x_plus_y(mean_x, std_x, mean_y, std_y, N)>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plot_x_plus_y(mean_x, std_x, mean_y, std_y, N):\n",
    "    x = np.random.normal(mean_x, std_x, N)\n",
    "    y = np.random.normal(mean_y, std_y, N)\n",
    "\n",
    "    q_of_x_y = x+y\n",
    "\n",
    "    pred_mean_q = mean_x+mean_y\n",
    "    pred_std_q = np.sqrt(std_x**2+std_y**2)\n",
    "\n",
    "    counts, bins, patches = plt.hist(q_of_x_y, \n",
    "                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), \n",
    "                                     density=True, alpha=0.3)\n",
    "\n",
    "    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)\n",
    "    plt.legend(('pred','hist'))\n",
    "    plt.xlabel('q(x)')\n",
    "    plt.ylabel('p(x)')\n",
    "\n",
    "\n",
    "# now make the interactive widget\n",
    "interact(plot_x_plus_y,\n",
    "         mean_x=(0.,3.,.1), std_x=(.0, 2., .1), \n",
    "         mean_y=(0.,3.,.1), std_y=(.0, 2., .1),\n",
    "         N=(0,10000,1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single variable example: division\n",
    "\n",
    "As a warm up, let's consider $q(x) = 1/x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'q(x)')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX6x/HPmZmUSQ/pPYQSeq+CooCoiFhWXHtlwbZi2113fyq6tl131cVdd+2uK2JDxYKAFAFFBBKQGgKhhPQC6T2Z8/vjxgjSAmRyJ5nn/XrNayaZydzvpTxzcu65z1Vaa4QQQnR+FrMDCCGEaB9S8IUQwk1IwRdCCDchBV8IIdyEFHwhhHATUvCFEMJNSMEXQgg3IQVfCCHchBR8IYRwEzazAxwuNDRUJyYmmh1DCCE6jNTU1GKtdVhrXutSBT8xMZGUlBSzYwghRIehlMps7WtlSkcIIdyEFHwhhHATUvCFEMJNSMEXQgg3IQVfCCHchBR8IYRwE1LwhRDCTbjUOvzT4dAOUgtS0VozImqE2XGEEMJldfgR/tLMpdy65FbmbJpjdhQhhHBpHb7gnx1zNnabnS1FWzhQfsDsOEII4bI6fMH38fDh/ITzAfhy75cmpxFCCNfV4Qs+wMVJFwNGwddam5xGCCFcU6co+CMjRxJmDyOrIostxVvMjiOEEC6pUxR8q8XK5K6TAfhizxcmpxFCCNfUKQo+wJRuUwBYsn8JDU0NJqcRQgjX02kKfnJwMt2DulNaV8qa3DVmxxFCCJfTaQq+UoopScYoX1brCCHE0TpNwYefV+t8c+AbKuorTE4jhBCupVMV/EjfSIZHDqfeUc+yzGVmxxFCCJfSqQo+wCVJlwAyrSOEEL/U6Qr+xISJeFo82ZC/gfyqfLPjCCGEy+h0Bd/f059z485Fo1m4d6HZcYQQwmV0uoIPHLFaR1otCCGEoVMW/LExYwnyCiKjNIP0knSz4wghhEvolAXfw+rBBYkXAPDlHjl4K4QQ0EkLPhw5rSOtFoQQohMX/IFhA+ke1J2DtQdZmrnU7DhCCGE6pxd8pZRVKbVJKdWucytKKa7pdQ0A83bOa89NCyGES2qPEf4sIK0dtnOUKUlT8PfwZ3PRZrYf3G5GBCGEcBlOLfhKqVjgYuB1Z27neHw8fLisx2UAzEuTUb4Qwr05e4T/D+D3gMPJ2zmua5KvQaFYvG8xh2oPmRVDCCFM57SCr5SaAhRqrVNP8roZSqkUpVRKUVFRm+eIC4jj7NizqXfU88nuT9r8/YUQoqNw5gh/DDBVKbUfeB8Yr5Sa+8sXaa1f1VoP01oPCwsLc0qQnw7efpD+AY2ORqdsQwghXJ3TCr7W+o9a61itdSJwNbBCa329s7Z3ImdFn0VCQAL5VfmszFppRgQhhDBdp12HfziLssgSTSGE22uXgq+1Xqm1ntIe2zqeqd2mYrfZ2ZC/gV0lu8yMIoQQpnCLET4YbZOndpsKwHs73zM5jRBCtD+3KfgA1/a6FoCFexdSVldmchohhGhfblXwk4KSGBU1iprGGhZkLDA7jhBCtCu3Kvjw8xLN93a+J0s0hRBuxe0K/rjYccT5x5FTmcOifYvMjiOEEO3G7Qq+1WLlN/1/A8CrW16lydFkciIhhGgfblfwAaZ0m0KMXwz7y/ezZP8Ss+MIIUS7cMuC72HxaBnlv7LlFRnlCyHcglsWfDBOxIryjWJv2V6WHpArYgkhOj+3LfgeVg+m958OwCubX8GhTevgLIQQ7cJtCz7AZd0vI8IngozSDJYfWG52HCGEcCq3LvieVs+WUf7Lm1+WUb4QolNz64IPcHmPywm3h7OrZJe0ThZCdGpuX/C9rF7c2v9WwBjla61NTiSEEM7h9gUf4Fc9fkWoPZS0Q2mszl5tdhwhhHAKKfiAt82bW/reAsB/Nv9HRvlCiE5JCn6zacnTCPEOYfvB7SzJlLNvhRCdjxT8ZnabnbsH3w3ACykvUNdUZ3IiIYRoW1LwD3N598vpGdyT3Kpc3tnxjtlxhBCiTUnBP4zVYuX3w38PwGtbXqO4ptjkREII0Xak4P/CyKiRnBt3LtWN1fxz0z/NjiOEEG1GCv4xPDjsQWwWG5/u/pS0g2lmxxFCiDYhBf8YEgISuKbXNWg0z254VpZpCiE6BSn4xzFzwEyCvIJIKUhhxYEVZscRQogzJgX/OAK9Arlz0J0APJf6HPVN9SYnEkKIMyMF/wSm9ZxGt8BuZFVkMS9tntlxhBDijEjBPwGbxcaDwx8EjEshyjJNIURHJgX/JMbGjOWc2HOobKjk6XVPmx1HCCFOmxT8Vnh45MP42HxYmrmUZZnLzI4jhBCnRQp+K0T5RXHf0PsAeGrdU5TVlZmcSAghTp0U/Fa6KvkqhoQPobimmOdSnjM7jhBCnDIp+K1kURYeO+sxPC2efJrxKWtz15odSQghTokU/FPQNbArdwy6A4DH1z5OdUO1yYmEEKL1pOCfopv63kSvLr3IqcyR5mpCiA5FCv4p8rB48PhZj2NVVt5Ne5fNRZvNjiSEEK0iBf809Anpw019b0Kjmb1mtrRdEEJ0CFLwT9MdA+8gISCBPWV7eCH1BbPjCCHESTmt4CulvJVS65VSm5VS25VSjztrW2bwtnnzzNhnsCkbc9PmsjJrpdmRhBDihJw5wq8DxmutBwKDgAuVUqOcuL121z+sP7OGzALgkTWPkF+Vb3IiIYQ4PqcVfG2obP7So/nW6a4kcmPfGxkTPYbSulL++O0faXI0mR1JCCGOyalz+Eopq1LqR6AQWKq1XufM7ZnBoiw8NfYpQu2hpBSk8OrWV82OJIQQx+TUgq+1btJaDwJigRFKqX6/fI1SaoZSKkUplVJUVOTMOE4TYg/h6bFPo1C8vPllUvJTzI4khBBHaZdVOlrrUmAlcOExnntVaz1Maz0sLCysPeI4xejo0dzW/zYc2sFD3z5EaW2p2ZGEEOIIzlylE6aUCmp+bAcmAjudtT1XcOegOxkYNpCC6gIe+f4Rufi5EMKlOHOEHwV8o5TaAmzAmMP/0onbM52HxYO/nvNX/D39WZm1kje3vWl2JCGEaOHMVTpbtNaDtdYDtNb9tNZ/dta2XEmMXwxPjnkSgDkb57Aqa5XJiYQQwiBn2jrB+Pjx3D3objSaP3z7BzJKMsyOJIQQUvCdZcaAGVyYeCFVDVX8dsVv5SCuEMJ0UvCdRCnFn8f8md5depNdmc2Dqx6kwdFgdiwhhBuTgu9EdpudF8e/SIh3COvy1/Hs+mfNjiSEcGNS8J0s0jeSf5z3DzwsHryf/j4fpn9odiQhhJuymR3AHQwKH8Ts0bN5eM3DPLPuGeID4hkV1an6yAknWbaj4Jjfn9gnop2TiM5ARvjt5NLul3Jz35tp1I3MWjGL7cXbzY4khHAzUvDb0X1D72Ny18lUN1Zzx7I72Fe2z+xIQgg3IgW/HVmUhSfHPMmYmDGU1JUwY+kM6aEvhGg3UvDbmYfVg+fHPc/AsIHkV+Uzc+lMWaPv7srzIDsVqopB+i8JJ5KDtibw8fDhpQkvcfPim8kozeCu5Xfx2qTX8PHwMTuaOEWnfFC1ogDyfoTcTZDbfF952G953oEQ0gNCukNod8LqoiiOHIe2ejohvXA3UvBNEugVyMsTX+bGRTeypXgL9628j3+N/xceVg+zowlnyNoAyx6DzO+Ofs4rEILioCQTassgJ8W4AQOBWnsk+3vNIDdpGg6rV7vGFp2LcqUWvsOGDdMpKe518ZDM8kxuXHQjh2oPcV7ceTw37jkp+h3ISUf4xbth+eOQ9oXxtYcvxAyB6EEQNQiiB0NwV7BYjOmcykI4mAEHd8PBDCq3LcavfDcAtfZwMnvNIKfrVYwfkNAeuyc6AKVUqtZ6WKteKwXffGkH05j+9XTK68sZFzuO5899Hk/5Fb5DOG7BjwNW/QVS3wbdBDY7jL4Txswypm1a+/7b8wjPWUrX7f/CvywdgDrvMLzG3QfDp4NNRvzu7lQKvhy0dQG9Q3rz+qTXCfIKYlX2Ku755h5qG2vNjiVOh9bEp78BLw6GlDcBDUNuhHs2woRHT6nYA6AsFMZewLpJn7F5zEuUB/XBq7YIlvwJ/jvFOCYgRCvJCN+FpB9K5zdf/4aSuhJGR41mzvg52G12s2OJEzh8hK+a6umd+ijR+z8xvpF8MUycDWHJbfL+AGhNaO4Kem18HO+afGrtEWwZ8xLlXQa0vETOwnUvMsLvoJK7JPPmBW/SxbsLa/PWcvfyu6luqDY7lmgFW30Zg7+dTvT+T2iy2uHXc+GaeWdU7I9JKYpjJrD+/E8oDR2Cd00BQ1dcS2TmZ227HdEptXqEr5QKBqKBGmC/1trR1mHcfYT/k72le7nt69sorilmaMRQXprwEr4evmbHEsewbEcB9soDDPr2N/hW7KPOO4wfx77CyLETjvv6tqKa6kne9ASxez8AYH/ybWT0f5CJ/aLbbBvC9bXZCF8pFaiU+pNSaivwA/AK8CGQqZT6SCl13pnHFb+UFJTEWxe8RbhPOKkFqdy25DYO1hw0O5Y4hsDiVIYvn4ZvxT4qApNZP+EjKrr0a5dta6snO4c9QdqQx3EoG4npbzDouxlQU9Iu2xcdz8mmdOYDWcDZWutkrfVYrfUwrXUc8BfgUqXUbU5P6YYSAxP57wX/JcYvhu0Ht3PDohs4UH7A7FjicNs+ZujKG/GsK6E48mxSxr9HnW/7j65zul/DxnPfpt4rmND8b2Hur6Cust1zCNd3woKvtT5fa/2O1vqoc/+11qla63u11m84L557iwuIY+7kufQJ6UNWRRbXf3U9W4u2mh1LAKQvho+nY3E0kNXtWjaPfYUmDz/T4pSGDWf9xE+o8Y2FnFT48AZorDctj3BNrTpo+8tRvFLKqpSa7ZxI4nCh9lDeuuCtloZrt319G6uyVpkdy71lp8L8W0A72Nf7dtKHzEZbzD9pvdY3ho3nvAm+YbBnBXw6ExxtfqhNdGCtXaUzQSn1lVIqSinVD2M+39+JucRhfDx8+Of4f3Jpt0upaazhnm/uYf6u+WbH6pSW7Sg45q3FwT0wbxo0VMOg69jT7z5QyrzAv1DjnwjXzQdPf9j+CSz6vTRkEy1aVfC11tcCbwNbga+Ae7XWDzozmDiSh8WDJ8Y8wcwBM3FoB4+vfZw5G+fgaPvFUuJ4KouM+fHqg9BtAlwyx6WKfYvoQcaSUKsnbHgNVsm1lIWhtVM6PYBZwMfAfuAGpZS0dmxnSinuHnw3j4x6BIuy8PrW17lnxT1U1FeYHa3TW7FlH2VvXg4l+ygP7ss3A/7GsvRDZsc6vq7nwJVvgrLAyqdhw+tmJxIuoLVTOl8Aj2itZwLjgN3ABqelEid0VfJVvDzxZQI8A1iVvYprF14rV89yIuVopP/aewk8tJUa31h+HPuqqQdoW633JTDlH8bjhQ/C9gXm5hGma23BH6G1Xg6gDc8BlzkvljiZ0dGjef/i9+ke1J395fu5duG1rM5ebXaszkdrkjc+TljeSuo9g9h09uvU28PMTtV6Q2+CCbMBDQvuhKJdZicSJjrZiVdjAbTW5b98Tmu9WykV0HwQV5ggLiCOdye/y8T4iVQ2VHL38rt5fevruFJ/pI4uav8nxO79gCarF5vHvkJ1QJLZkU7d2Pug/zRoqIIPb4T6KrMTCZOcbIT/K6XU90qpR5VSFyulRiilzlFK3aqUegf4EpDuXiby8fDhuXOf465Bd6HRzNk4hwdWPUB5/VGf0eIU+ZbtptfGxwHYOfRxykIHm5zoNCllTO2E9oSiNFj4gKzccVMnXDystb6vuYfOlcA0IBKjl04a8LLWeo3zI4qTsSgLtw+8neTgZP743R9ZmrmUHQd38Ldz/kb/sP5mx2sXp3ypwZOwNNbQf+0srE215CZcRl7iFWcSz3xefnDV/+C18bD5PUg4y2jbLNzKSefwtdYlQACwBVgKfAcUAyOVUvcrpe53bkTRWufFn8eHUz6kd5fe5FTmcOOiG3l7+9uydPM0JG96Ar/yDKr8k0gf0knOMQzvDVNeMB5/9TvIl7O23U1rTw8cCgwHPgMUcAmwGqPPjnAh8QHxzJ08lxdSX2Bu2lz+nvJ31uev58kxTxLsHWx2vA4hMvMzYvbNp8nqxZbRc2g6jU6lbdkVs00NvBoy18DG/xnz+TNWgXeA2alEO2ntKp1QYIjW+kGt9QMYHwCxWuvHtdaPOy+eOB2eVk/+MOIPzDlvDgGeAazOXs2VX1zJhnxZSXsyPuV76ZVqjOjTBz1MVVAb97N3BRc9CxH94dBe+Pxumc93I60t+PHA4Z2Y6oHENk8j2tT4+PHMv2Q+g8MHU1hdyG1LbuOv6/8ql088ngZj3t7WWE1+3MXkJl1ldqLTctL2EB52uOpto/3Cjs9g/avmhRXtqrUF/x1gvVLqseamaeswWi0IFxflF8WbF7zJ7QNvx6IszE2by7QvprGlaIvZ0VzPkj/hX5ZOtV8CacOecM22CW0lpBtc+i/j8dePyPp8N9HaXjpPAbcAJUApcIvW+hlnBhNtx2axcdegu3h38rskBSaxv3w/Nyy6gRc3vkh9k7TQBWDXEkh5E4fFg62j53SMM2nPVN/LYNB10FQHC26HpkazEwkna/U1bbXWG7XWc5pvm5wZSjhH39C+fHjJh9zS9xa01ry29TWuWXgNOw/tNDuauWpK4PN7AMjofz8VwX1MDtSOLnwGAmKMHvrfzzE7jXAyuYi5m/GyenH/sPt5+6K3ifOPY1fJLq7+8mqeT3nefS+YvviPUJkPcSM50ONms9O0L+9AmPpP4/E3z0DBdnPzCKdyWsFXSsUppb5RSqUppbYrpWY5a1vi1A0OH8z8S+ZzXe/rcGgHb21/iys+v8L9+vGkLzZORLJ5w6X/BovV7EROc9yDud0nwNBbwNEAn94OTQ1mRxVO4swRfiPwgNa6NzAKuEsp5Ua/K7s+Hw8fHhrxEPMunkevLr3IqczhruV38cDKByiqLjI7nvPVlMAXzeOQ8Y9AaHdz85hp0hMQFA/5W2D1381OI5zEaQVfa52ntd7Y/LgCox1DjLO2J05fv9B+vHfxezw47EHsNjtfZ37N1AVTeW/nezQ6OvGBvMOmchh1h9lpzOXlb/yGA/Dt3yH3R3PzCKdolzl8pVQiMBhjOecvn5uhlEpRSqUUFbnBqNJF2Sw2bup7EwsuXcC42HFUNlTy9LqnuerLq1iXd9RfW8eXvshtpnJarevZMPJ2cDTCgjugsc7sRKKNOb3gK6X8MK6Ude9x2iy/qrUeprUeFhbWgfqMd1LRftH8c/w/+ce5/yDGL4bdJbuZ/vV07v3mXrIqOkknjZoS+OJe47G7T+X80oTZ0KUbFO6AlX8xO41oY04t+EopD4xi/67W+hNnbku0HaUUExIm8NllnzFryCzsNjvLDyzn0gWXMmfjHKoaOng/9ZapnFEylfNLnj5w2X8ABWvmQJ6coNeZOHOVjgLeANK01s87azvCebysXkzvP50vL/+Sqd2m0uBo4PWtrzPl0yl8sPMDGhwdcDXH7qWHTeW8JFM5xxI/EkbOBN0EX9wDjiazE4k24swR/hjgBmC8UurH5ttkJ25POEm4TzhPjX2Kdye/y4DQARTXFPPkuie5bMFlLN63uOO0X66vgi+bu3mf9yeZyjmR8Q8bJ2TlboL1r5mdRrQRZ67S+U5rrbTWA7TWg5pvXzlre8L5BoQNYO7kuTw37jkSAxI5UHGA363+HVd/eTXf535vdryTW/kMlB2AyP4w6i6z07g2L3+4+Dnj8fI/Q2knOX7j5uRMW3FKlFJMSpzEJ5d+wqOjHyXMHkbaoTRmLp3J9CXT2VTool038rbA2n+DssAlc8Da2ktBuLHki6DPpca1cL96UNoodwJS8MVp8bB4MK3nNBZesZBZQ2bh7+HPuvx13LjoRqZ/PZ3UglSzI/7M0TwXrZtgxAyIGWp2oo7jomfBKxB2LTZaKYsOTQq+OCN2m53p/aez6FeLmDlgJn4efqzLW8fNi29m+pLppOSnmB3RmIPO3WTMSY9/2Ow0HYt/JJz/mPF40e+hptTUOOLMKO1Cv6YNGzZMp6S4QIEQp62sroy5aXOZu2MulQ2VAAyPHM70ftMZHT0a5aQe88e7pKBXdR6jF0/G1ljF5jH/pihmInD8i5u77KUJTdDyZ+RwwFsXQdYPRs+dS/5hbjBxBKVUqtZ6WGteKyN80aYCvQK5a9BdLP7VYu4YeAf+Hv5syN/AzGUzuerLq/hq71ft165Ba3pt/DO2xioKY85vKfbiFFmaj3tYPCD1Lchca3YicZqk4AunCPQK5M5Bd7L4ysXMGjKLEO8Qdh7ayR++/QNTPp3CvLR51DTWODVDWM7XhOUup9HmS/rgR5y6rU4vvBec3byk9YtZ0Hj0hXNOemlFYTop+MKpAjwDmN5/OkuuXMKjox8lISCBnMocnln/DJPmT+LFjS9SWF3Y5tu1NlSSvOkJwLioSZ1PZJtvw+2MvR9CukNxOnz/otlpxGmQOXzRrpocTXyT9Q1vbnuTrcVbAbApG+cnns91va9jYNjA03rfX44ke256kvjd/6OsywA2jP9AzqhtI8EFaxm66ibjTOU710KXpJbnjjeaP97xEtE2ZA5fuCyrxcrEhIm8O/ld/nfR/5iUMAmNZtG+RVz/1fVcu/BaFu5dSMMZXITD/9A24jLm4lBW42LkUuzbTEnEaBjwa2isha9+J2vzOxgp+MIUSikGhw/muXOfY9EVi7i1360EegWytXgrD337EBPnT+SF1BdOvUOno4neqY+itIOsHjdSGdTbOTvgziY9aVwaMWMZ7FhgdhpxCqTgC9NF+UVx39D7WHrlUmaPnk2P4B4cqj3Em9veZPInk5m5dCbLM5e3anVP7J55BJRso9Yeyd6+v22H9G7ILxwmPmY8XvQQ1B7V9Vy4KCn4wmXYbXau7HklH1/yMe9c9A5Tu03Fy+rF97nfc+/Ke7lg/gXM2TiHzPLMY/68Z00B3bcajVnTBz9Mk4dfe8Z3L0NuhtjhRpvpFU+anUa0khR84XKUUgwKH8RTY59i+bTl/H7470kMSKSwprClPfNNi27i092fHtGbv+ePz2BrrKIo6jyKYs43cQ/cgMUCU14AZYUNr0HORrMTiVaQVTqiQ9Bak1qQyoKMBXyd+XXLGn67zc75CefTpzyIa9Y/ibbaWXvhV9T6yuWTneWIVTdL/g/W/guiBrFszHvHPEAuq3Sc61RW6UjBFx1OdUM1S/YvYUHGAjYW/jyyDG9s5Cx7P7r3vI9oe5LT2ji4uyMKeF0lvDQCynPYOfgRsnvccOLXizYnBV+4jQPlB/hs+e9YeHALOR4/tzyOsicyJGQCQ7qcRxcvOemqLR1VwNO+hA+uo9Hmy9oLFx11kpsUfOeSdfjCbcTX1fLbrctZlJ3LnyLvYEzYJfjaAsir2c/C7Dd4Ysv1vLDjbr7J/4iSOjnN3yl6T4Fko0Fdzx+fNjuNOAG5CoTouLSGL+8DRwO5Xa8kLO5XXAlcFn8nO8tT2HhwBdtL13KgaicHqnbyedYrJPj2ZlCXcQwMPptgLxl5tpmLnqVxz0oisheTl/sNxdHnmZ1IHIMUfNFxbZoLmd+BTyi7B/y+5ds2iwf9gkbTL2g09U21pJWt58dDq9hRto7MqjQyq9L4LOtlYn16MiB4DAOCxxJhTzBxRzqBoDj29p1Fz83PkLzxcUrCRtDk4Wt2KvELMocvOpSf+rV41B7krMUX4lFfxraRfyc/YepJf7auqYYdZev48dAqdpZtoN5R2/JcuHcc/YPG0C/4LOJ9e2FRMtt5qpSjkeHLpxFQsp3Mnreye9BDgMzhO5sctBWd1k8Fv++6B4nK/JyDEWPZdM4bcIorcuoddewqS2Vr6Rq2lXxPdVNFy3N+tiD6BI2kb9BokgOG4mW1t+k+dGb+h7YxYvmVaBQbJn5MRXAfKfhOJgVfdGgn66HeJX8NQ1bfQpPVix8uWEiNX/wZba9JN7G3YitbS75jW+laSup/3r5VedDDfyB9gkbRO3AEod7RZ7Qtd9Bz01PE736bsi792TD+w+M2r5MPgrZxKgVf5vBFh2JprKXXxtkA7O1z9xkXewCrstIjYBA9AgZxefxd5NXsY0fpD2wv/YHMqjR2lqews9wYiIR6xdA7cDi9A0fQzX8AnlbvM95+Z7On3yzCsxcTeGgrsXvmHXNtvjCHFHzRoXTd8RI+lQeoDOzJgeRb2/z9lVJE+yQR7ZPExOhrqWgoIa1sPWml60kvT6W4LodvC3P4tnABHsqTJP/+JAcMpWfgUKLsXWXuH2jy8CN9yKMMXHMX3bc+T1HM+XIBGhchUzrC5RxvSse3NJ2RSy9H6SZSxr9PWejgds3VpJs4ULmTneUbSCtdT1b1riOe97MF0TNgMD0DhtIzYAjBXuHtms/VDPjuDsJzl1MYcz5bxrzU6p+TqZ5TI1M6ovPRDnqnPopFN5LV7dp2L/ZgTP109e9LV/++XBRzMxUNJewu30R6eSq7yjZS2lDExkPfsPHQN4Ax/dM9YCA9/AfTI2AQ/h7B7Z7ZTOlDHqVL4VrCc5YSlr2EotgLzI7k9qTgiw4hds97BB3cRJ13OBn9HzA7DgD+HsEMCRnPkJDxaK0prM0yin/5RvZUbKa4Lofiohx+KPoKgAjveLoHDKKb/wCS/PoT6Bli8h44V51PFHv630/ypifptfFxSsJG0ugVZHYstyYFX7g876psum/5GwDpQx6hydPf5ERHU0oRYY8nwh7PORGX06SbyK7aTUbFj2SU/8jeym0U1B6goPYAawo/ByDMK4Yk/wF08+9PN/8BBHtGdLqGb1ndryci6yuCijfSc/Mz7BjxV7MjuTUp+MK1aU3vlEewNVZTEHsRhR1kWsCqrCT49SLBrxcToq6m0dHAgap09lRsYW/FFvZVbqeoLoeiuhzWFS8CINAjhES/vnT160dX/77E2LthtXTw/6LKwo5hTzPy66lE7/+UgriLORh1jtmp3FYH/9ckOrvofR8RUrCGes8gdg551Ow4p81m8SDJvx9J/v2Aa2nSTeRUZ7RbNuoeAAAWi0lEQVR8AOyt2EZZw0E2l6xmc8lqADwt3sT7JpPo14cEvz4k+PbqkMcBqgOS2Nv3Hnps/Tu9Uh/lhwu+lKuRmUQKvnBZXtX59Nz8F8A4ANjg3XnmvK3KSrxvMvG+yZwXOQ2HdlBUm8W+yu3GrWIbRXU5ZFRsJqNic8vPhXhFkeDbm0S/PsT79iLGJwmbxdPEPWmdA8m3EpG9mICSbXTf8nfShz5mdiS3JAVfuCat6Z3yMLaGSgqjJ1AQd7HZiZzKoixE2BOIsCcwKmwyABUNJeyv3EFmZRr7q3aQVbWLg3V5HKzLY+OhFQBYlY1oexLxfr2aP0B6Ee4di0Ud++xWs2iLjR3Dn2bE0iuI2zOPgrjJlIaPMDuW25GCL1xSVOanhOavpsEjgJ1DHz/lXjmdgb9HMP2Dx9A/eAxgnAeQX7OP/ZVp7K/cQVZVOoW1WWRV7yKrehdrmn/O0+JNjE934nx7EufTkzjfHoS5wIdAZVAv9veeSdKOl+iT8n/8MOlzHDbpU9SepOALl+NZU0jPTcaFNHYN/hP1dvc+geknVmUlxqc7MT7dGRN+CQC1TVVkVe1u7vmfzoGqnZTWF7Gvchv7Kre1/KyXxU60TxIxPt2JbX6PSHsiNotHu+7Dvt53EJ79NX7lu0na/iIZA//Qrtt3d1LwhWvRml6ps/FoKKc48hzyEi43O5FL87b6tvQB+kllQylZ1bvIrtpNVpUx+jc+BIzjAz+xKhsR3gnE+HRr/jAw7n1tgU7Lq62e7Bj+NMNX/JqEXW9RFHM+ZaFDnLY9cSRprSBcy5aP4JPpNHr4sfaCr6QHSxupbCgluzqDnOZbdnUGxbU5aI7+/x/oEUq0TxJR9q4t9+HecW3620D3zX8jMf01qn3jWDfps1at2pGWC8cmrRVEx1R6ABYaZ9HuGviQFPs25OcRRK/AYfQK/Lku1DZVk1u9l9yaPcZ99R7yavZT1lBMWVkxaWXrW15rUVbCveOIsnclyp5IpD2RSHsCIV5Rp3VsYE+/WYQUfId/aRo9f3yatOFyLdz24LSCr5R6E5gCFGqt+zlrO6KTcDTBJzOhrozC6Ankdp1mdqJOz9vqc9i5AQaHdnCwLpec6r3k1ewjr3ofeTX7OFiXS37NfvJr9rPpsPfwUJ6E2+ONDwDveGOlkXc8Id7RWE/wQaCtnmwb+Rwjll1OzL75FEeNk1477cCZI/z/Av8C/ufEbYjO4rvn4cD34BdB2vCn3HJVjiuwKAth3rGEeccyiJ/PiK1rqqGgJtP4EGgu/Pk1mZQ1FLdMEx3OqmyEeccS4W20mwj3jmu5/XQFsarA7mQM+B3Jm56kd8ojlIUMot4u0zbO5LSCr7VerZRKdNb7i04kOxW+ecZ4fNl/aKjvYm4ecRQvq91Y6+/X64jv1zRWkl+TSX7tfgpqjF5BBTWZlNQXtvxGQMmR7xXoEUq43Sj+Yf4xjI4eRv/CTfRa/xBbznkD5JoCTiNz+MJcdZXw8W2gm2DUXdB9ApzkEofCddhtfi0tow9X11RDYW0WBTUHKKw9QEFtFoW1WRTV5hjHCBqK2V1uTA4t8ALiYvDQB4jYOI2ggL6EesUQ5m3cQr2iCfQMM2HvOh/TC75SagYwAyA+/swvVyc6mMV/gJJ9ENEPJnTcXjniSF5Wu3Hil2/PI77v0E0cqiugoPZAywdAUW02B6syKHFUku0oI7v0+6Pez0N58p+MeGL9Y4nzj2u5xfvHE+UXhUc7n0/QUZle8LXWrwKvgrEs0+Q4oj1tXwCb5oLNG371OnjI9WE7O4uyEuodTah3NH0ZdcRziRv+QGPOl+wMTGRtj19T3FBAUW02xbW5VDSWkFGaQUZpxjHe00KUbxSxfrHE+jff/GKJ8Yshxj+GYK/gTtd2+nSZXvCFmyrLgS9mGY8nPQnhvc3NI0yXNegRRham0KtoN0MD95E+ZHbLc7VNVSRF15JVkUVWRRbZFdktj/Or8smpzCGnMod1+euOel+7zW4U/1/covyiiPGLIcAzwG0+EJy5LPM94FwgVCmVDczWWr/hrO2JDqSx3pi3ry2FHpNg+HSzEwkX0OThx7bRLzBsxTXEZbxLWchg8hOmAsYZxX1CkugT0ueon6tvqievKo/simzjVpnd8oGQU5lDZUPlcX87APCx+RDtF020XzRRvlHGY99oIn0jifaLJtQe2mkuTu/MVTrXOOu9RQf39cNwYC34R8Ol/5YlmKJFeZcBpA9+hN6pj9I75WEqApOpCko+4c94Wj1JCEggISDhmM+X1ZWRU5lDbmUuOZU5ZFdkk1eV1/K96sbqE34g2Cw2InwiiPSNJNI3kijfKCJ9Ilu+jvSN7DC/JciUjmhfm9+H9a+A1RN+/Q74yeoLcaScpF8TeHAT0fs/ZcD3v2X9xI/P6LKWgV6BBHoFHvO3A6015fXlLcU/ryrPuFXmtTw+VHuoZcroeOw2OxE+EcbNN6LlcbhPOBG+xn0X7y6m/6YgBV+0n7zNP8/bX/Qsy8rjZAmmOJpS7BzyGP4lO/AvS6fvhofYcta/nLQpdcIPBIDaxlryq/LJr84nvyqfvKo8CqoKWu7zq/Opaqhif/l+9pfvP+62bBYb4fZwwn3CCfMJa/lA+Olx35C++Hj4OGU/WzI49d2F+En1IfjgemishcHXw9CbIa3Q7FTCRTlsdraMeYkRSy8nPGcpCelvQN//MyWLt82bxMBEEgMTj/uayvpKCqoLjA+GqnwKqwspqC5ouRVWF1JWV0ZuVS65VbnHfI/5l8wnucuJp6/OlBR84XyOJvh4utEcLXowTH5O5u3FSdX4xbN9xLMMWnMH3bf+HQafA13PNjvWMfl5+uHn6Ue3oG7HfU1NYw1F1UUUVBdQVF1EYXUhhTWFFFYXUlRdRISP89tKSMEXzvfN07BnOfiEwFXvyHp70WrFMRPY12smXXe+AvNvgRmrIDDG7FinxW6zEx8QT3yAeSeYdo61RsJ1pX0J3/7d6I9y5ZsQFGd2ItHB7Ol3LwfDR0NVEcy7CmrLzY7UYUnBF86TtcGYygGY+BgknWteFtFxWaxsHT0HQnpAwTb48AbjXA5xyqTgC+co3m2MxhprYND1cNY9ZicSHVijVxBcPx98w2DvSmO1lwtdra+jkIIv2l5FPrxzBdQcMs6kveQfcpBWnLngRLj2Q/Dwgc3zYOUzZifqcKTgi7ZVWwZzr4SyAxAzFKb9F6zSyVC0kZghcOVbxjGhVX+FjXJ9pVMhBV+0ncY6Y619wVbo0s0YjXn6mp1KdBLLdhQYt6bBpA15DADHF/fC7mXmButApOCLtuFwwKe3w77V4BcBN3wCvqFmpxKdVE63q9nXayYW3QQf3QQ5qWZH6hCk4Iszp7VxIZPtn4CnP1z3kTHfKoQT7el/H3nxU6G+Ev53OWSnmB3J5cmJV+LMOByw8H5IfQuHxYMfR7/IoZJIKPm5R87EPnJhauEEysKOEc9gcdQTkb2Yxv9OZdM5r1MWOrTlJfJv70gywhenr6kRFtwBqW+BzZvNY/7NoYgxZqcSbkRbPNg26nny4y7G1ljF4NW3EVS0wexYLksKvjg9jfXw8a2w5X3w8IXrPuJg1DizUwk3pC02to/8G3nxU7E1VjN49XSCC38wO5ZLkoIvTl1DrbEaZ8dn4BUAN3wKXc8xO5VwY9piY/uIv5KbeAXWphoGfTuDLgVrzI7lcmQOXxzXsmP0qrc0VjP+x3uNsx3twUaxjx7c/uGE+CWLlR3Dn0YrKzH7PmLgtzMhygt6XWx2MpchI3zRap41hQxZdbNR7H3D4eavpNgL16IspA17gqxu12J11MP718G3z0sbhmZS8EWrBBzczIhlVxB08EcIiIVbFkHEsa8QJISplIX0IbPJ6H8/oGH54/DJDGMq0s3JlI44qah9H9M79VEsjgZKQocRfPN74Bfe6p8/1tSQEE6lFPt73073PkONYr/1Qzi0B66eB/6RZqczjYzwxXEpRwM9Nz5B3w1/xOJoIKv7dWwc999TKvZCmKr3FLjtawiMN87GffU8yN1kdirTSMEXx1ZVzJBVtxCf8Q4Oiwc7hj1F+pDZaKun2cmEaLVlOwpYdiiMVeM+oCR0KFTk0vTGhbDpXbec15eCL46262t4+WyCi9ZT5x1Oynnvkps0zexUQpy2Bu8QNo57m5yuV2JtqoXP7jQO6FYWmh2tXSntQp9yw4YN0ykp0g/DWY43l95y+nlNCSz+I2x+D4DS0CFsGf0i9XaZwhGdhNZEZS4gedMT2BoqqfcMYufQP1MYdyHQdq0YTvp/rQ0ppVK11sNa81oZ4QvDzoXw0kij2Nu8YdJTpJz7rhR70bkoRV7i5ay9YCEHI87Cs76UAWvvoe8PD2CrKzU7ndNJwXdzHnWHYP5t8P61UFkA8aPh9jVw1t1gsZodTwinqPOJYtM5b7FzyGyarHaiDnzB6CVTYMfnnXpuX5ZluilLUx2xGfNITPsP1JeCzW5caHzEDLDIOEC4AaXI7n4dByPG0nfDHwgq3mhcID1uFEx6AuJGmJ2wzUnBdzPK0Uhk5md02/4i3tV5xjcTz4apL0KXJHPDCWGCGv8EUs59l9i979Nr50uQ9QO8cT70uRQmzIaQbmZHbDMylHMXWhOW/TUjv76Evhv+iHd1HhWByWwa+wrc9IUUe+HeLFayu18H92yCsx80fuPd8Rm8NAK++j1UFpmdsE3ICL+za6yD7Qtg3X8Y2HzCSbVvLHv73Ut+/BTjYtBKmRxSCBfhHQATHoFht8LKp431+utfgdT/Qv9pMOp2iOxvdsrTJgW/E1q2owCvqlxi97xHzL6P8Kw7BECdVwj7+txFTtJVR5xAJa0PhPiFwBi49CUYdSeseBLSF8GPc41bwlij8CdP7nALG6TgdyZNjbBvFQPW/Juw3BUo7QCgIjCZ7O7XkZcwFYfNx+SQQnQgEX3hmvfg4B5Y/ypsmguZ3xm3oHgYejP0vbzDTInKiVcdXWMd7F0FaZ9D+ldQfRAAh7JRGHsBWd2vM67xKdM2QpzUSU+Mqi2HH9+Fda9Ayb6fvx810Cj8fS6DLl1d9sQrKfgdUfUh2P8teT/MJyxvBbaGypanqvwSyU+4lJykq6i3h5kYUohOzNFESP5qIrMWEpazHFtjVctTZcH9KIo5n5LwUZR36Ye2eLQ8Z3bBlymdjqAiHzK/h8w1xn3hDgCifno6MJnC2AsojJ1EVUAPGc0L4WwWKwejz+Ng9HlYmuoIyf+W8KxFhOWuILBkG4El2wBotPlQGjqUkvBRlISPhKbzwGpe2ZURvivRGsqyoGA7FGwz7vM2w6G9R7ysyeJJecggiqPGURgziRr/BJMCCyEOZ2msJSR/NSEFawgu/AHfin1HvsDDxzguENkfIgcYt/De4Hn6x9ZcZoSvlLoQmANYgde11n9x5vY6BK2NDn0l+6E007gv2W8U9YIdUFd21I802nwoCxlCSdhwSsOGUdZloLQpFsIFOWzeFMVOoih2EgCeNQUEF66nS+EPBBetw6fyAGRvMG4/URYI7QlXvAZRA5yaz2kFXyllBV4CzgeygQ1Kqc+11jucsT1nHyQ57vv3CjOKdE0p1JZBbanxuKbEKOxVhUaPmsrm+4oCaKw57nbqvYKpCOxFZVAylS33PY+YBxRCdAz19ggKEi6hIOESADzqSvAr3Yl/6Q78S9PwL92JT/keLEU7wdf5x9ycOcIfAWRorfcCKKXeBy4F2rbg15ZB9gZC8g6BdqDQxnLE5sc4/EA7wNEIjibQTca9oxGaGsDR0Hzf/HVTvbHypanOuG+shcZ6BpeWYW2sxtpUbdw31mBtrIbGauAUp8XswRCcaNyCEiA4kY0VgVQG9KTeO1Tm4IXopBq8gimJGE1JxOiW71kaaxkfeqhdLr3ozIIfA2Qd9nU2MLLNt1KcAXN/xeA2f+MjhZzoSa8A8A4E7yCwBxmP7UHgF2HcfMN+fuwXbpzN9wuH5OQnIdySw+YN0c6uYAZnFvxjDVOPGgorpWYAM5q/rFRKpZ/gPUOB4jbI1sbKMT7PnMJF99mp3HGfwT332x33Gdp2v1u9asOZBT8biDvs61gg95cv0lq/CrzamjdUSqW09mh0ZyH77D7ccb/dcZ/BvP12ZrfMDUAPpVRXpZQncDXwuRO3J4QQ4gScNsLXWjcqpe4GlmAsy3xTa73dWdsTQghxYk5dh6+1/gr4qg3fslVTP52M7LP7cMf9dsd9BpP226XOtBVCCOE8csUrIYRwEy5X8JVSFyql0pVSGUqph47x/P1KqR1KqS1KqeVKqU7RSOZk+33Y665USmmlVIdf2dCafVZKXdX8971dKTWvvTM6Qyv+jccrpb5RSm1q/nc+2YycbUUp9aZSqlApte04zyul1IvNfx5blFJD2jujM7Riv69r3t8tSqnvlVIDnR5Ka+0yN4yDu3uAJMAT2Az0+cVrzgN8mh/fAXxgdu722O/m1/kDq4EfgGFm526Hv+sewCYguPnrcLNzt9N+vwrc0fy4D7Df7NxnuM/nAEOAbcd5fjKwCOPcnVHAOrMzt9N+n3XYv+2L2mO/XW2E39KOQWtdD/zUjqGF1vobrXV185c/YKzv7+hOut/NngCeBWrbM5yTtGaffwO8pLUuAdBaF7ZzRmdozX5r4KfTsQM5xvkrHYnWejVw6AQvuRT4nzb8AAQppaJO8PoO4WT7rbX+/qd/27RTLXO1gn+sdgwxJ3j9bRgjg47upPutlBoMxGmtv2zPYE7Umr/rnkBPpdQapdQPzd1XO7rW7PdjwPVKqWyMVW6/bZ9opjnV//edUbvUMle7AEqr2jEAKKWuB4YB45yaqH2ccL+VUhbgBeDm9grUDlrzd23DmNY5F2P0861Sqp/WutTJ2ZypNft9DfBfrfVzSqnRwDvN++1wfjxTtPr/fWeklDoPo+CPdfa2XG2E36p2DEqpicD/AVO11nXtlM2ZTrbf/kA/YKVSaj/GPOfnHfzAbWv+rrOBz7TWDVrrfUA6xgdAR9aa/b4N+BBAa70W8MbovdJZter/fWeklBoAvA5cqrU+6OztuVrBP2k7huapjVcwin1nmNOFk+y31rpMax2qtU7UWidizPdN1Vp35MuDtab1xgKMg/QopUIxpnj20rG1Zr8PABMAlFK9MQp+UbumbF+fAzc2r9YZBZRprfPMDuVsSql44BPgBq31rvbYpktN6ejjtGNQSv0ZSNFafw78DfADPlJG3/gDWuuppoVuA63c706llfu8BJiklNoBNAG/a49RkDO1cr8fAF5TSt2HMbVxs25eytERKaXew5iWC20+LjEb8ADQWr+McZxiMpABVAO3mJO0bbVivx/F6Lz+7+Za1qid3FBNzrQVQgg34WpTOkIIIZxECr4QQrgJKfhCCOEmpOALIYSbkIIvhBBuQgq+EEK4CSn4QgjhJqTgC3EcSqnhzb3KvZVSvs09+fuZnUuI0yUnXglxAkqpJzFaG9iBbK31MyZHEuK0ScEX4gSa+91swLgGwVla6yaTIwlx2mRKR4gT64LRu8kfY6QvRIclI3whTkAp9TnGVam6AlFa67tNjiTEaXOpbplCuBKl1I0YHQznKaWswPdKqfFa6xVmZxPidMgIXwgh3ITM4QshhJuQgi+EEG5CCr4QQrgJKfhCCOEmpOALIYSbkIIvhBBuQgq+EEK4CSn4QgjhJv4flCWi09XFmk0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "counts, bins, patches = plt.hist(x, bins=50, density=True, alpha=0.3)\n",
    "gaus_x = scs.norm.pdf(bins, mean_x,std_x)\n",
    "\n",
    "q_for_plot = 1./bins\n",
    "\n",
    "plt.plot(bins, gaus_x, lw=2)\n",
    "plt.plot(bins, q_for_plot, lw=2)\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('q(x)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'p(x)')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAELCAYAAAA2mZrgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8VPWd//HXhyRckoBRiEAJknAHK6AGrKJdbKmlrouP1kuh1Vatsqvr2t9269Z2f7UtvWx3bX+2du1atD60Wy+raJV6WVu3WusVguKFAHLHyC1yh3BL8vn9cWbiOOYySebMmcm8n49HHpnJ+WbmzYnmk+/lfI+5OyIiIgC9og4gIiLZQ0VBRERaqCiIiEgLFQUREWmhoiAiIi1UFEREpEVoRcHM7jSz7Wb2VhvHjzGz35vZ62a23MwuDyuLiIikJsyewl3ArHaO/z1Q6+6TgRnAT82sd4h5RESkA6EVBXd/DtjZXhOgv5kZUBpr2xhWHhER6VhhhO/9H8AiYDPQH/i8uzdHmEdEJO9FWRQ+DSwDPgGMAv5oZn9x973JDc1sHjAPoKSk5NTx48dnNKiISK5bunTpe+5e3lG7KIvC5cCPPdh8aY2ZrQfGA4uTG7r7AmABQHV1tdfU1GQ0qIhIrjOzjam0i3JJ6ibgkwBmNhgYB6yLMI+ISN4LradgZvcRrCoaZGZ1wHeAIgB3vw34PnCXmb0JGPANd38vrDwiItKx0IqCu8/t4Phm4Jyw3l9ERDovyjkFEZGMOXr0KHV1dRw6dCjqKKHq27cvFRUVFBUVden7VRREJC/U1dXRv39/KisrCS6P6nncnR07dlBXV0dVVVWXXkN7H4lIXjh06BADBw7ssQUBwMwYOHBgt3pDKgoikjd6ckGI6+6/UUVBRERaqCiIiOSo0tLStL+mJpoz5OnabR94XrSjninXXMoxoyvh4YchD7q1ItKxpqYmCgoKInt/9RQiUHDgAFOuvoRj3loGjzwCf/hD1JFEJAM2bNjA+PHj+fKXv8ykSZO48MILaWhooLKykvnz53PmmWfy4IMPsnbtWmbNmsWpp57KWWedxcqVKwFYv349p59+OlOnTuXb3/52KBlVFDLMjh7lpH+8kmOWv05zYWwd8Q9/GG0okXxjFs5HClatWsW8efN44403GDBgAL/85S+B4PqC559/njlz5jBv3jx+8YtfsHTpUn7yk59wzTXXAPDVr36Vq6++miVLljBkyJBQTo2KQiY1NzPxxn9k0AvPcOS4gSy573EoK4O//CX4EJEeb/jw4UyfPh2ASy65hOeffx6Az3/+8wDs37+fF198kYsuuogpU6bwt3/7t2zZsgWAF154gblzg80iLr300lDyaU4hg0bf/EOGLlpIY79ilv3yt+ybOAmuuw7mz4cf/QiefDLqiCL5wT2yt05eMhp/XlJSAkBzczNlZWUsW7Yspe9PN/UUMmT4bxZQeeetNBcW8ubNd7D3pJODA9ddByUl8D//A0uXRhtSREK3adMmXnrpJQDuu+8+zjzzzA8cHzBgAFVVVTz44INAcJXy66+/DsD06dO5//77AbjnnntCyaeikAn338+4f7sRgNrv38yOsz7x/rGBA+Hqq4PHP/pRBOFEJJMmTJjA3XffzaRJk9i5cydXx///T3DPPffw61//msmTJ3PiiSfy6KOPAvDzn/+cW2+9lalTp7Jnz55Q8plH2I3qipy7yc7WrTBiBBw5wuqv/V82fuXaDxyeOXEwbNkCVVVw+DAsXw4TJ0YUVqTnWrFiBRMmTIg0w4YNGzjvvPN46623Qn2f1v6tZrbU3as7+l71FNLs6dptH/h46zcPw5Ej7PjYWWy84u9b/6ahQ+GKK4LHP/5x5sKKiCQJrSiY2Z1mtt3M2iyJZjbDzJaZ2XIz+3NYWaJ07OIXAdg5fUb7S9auvx4KCuDee2GdbkAn0hNVVlaG3kvorjB7CncBs9o6aGZlwC+B2e5+InBRiFkic+ySoCjsmnpG+w2rquCLX4SmJrjppgwkExH5sNCKgrs/B+xsp8kXgIfdfVOs/fawskSlz9bNFG9aT2NJKfsmnNTxN9xwQ9CbuPNO2Lw5/IAiIkminFMYCxxrZs+a2VIz+1KEWUJxbE2w7Gz3qafhhSlcEjJhAnzuc3DkCPz0pyGnExH5sCiLQiFwKvDXwKeBb5vZ2NYamtk8M6sxs5r6+vpMZuyW+HxCh0NHib71reDzr34FPfy2gSKSfaK8orkOeM/dDwAHzOw5YDLwdnJDd18ALIBgSWpGU3ZDyvMJiU45JViSWlsLr78Op50WUjqR/Ja8c3F3zZw4uMM2bS1JvfHGG/n4xz/OzJkzW/2+Rx55hLFjxzIxA8vVo+wpPAqcZWaFZlYMnAasiDBPWnV6PiHRtGnB58WL0x9MRLLO/Pnz2ywIEBSF2trajGQJc0nqfcBLwDgzqzOzr5jZ35nZ3wG4+wrgf4A3gMXAHe6e3Wu1OqHT8wmJ4kXhlVfSnEpEotbU1MRVV13FiSeeyDnnnMPBgwe57LLLWLhwIQA33HADEydOZNKkSXz961/nxRdfZNGiRVx//fVMmTKFtWvXhpovtOEjd5+bQpubgB65/rJL8wlx8SEj9RREepzVq1dz3333cfvtt3PxxRfz0EMPtRzbuXMnv/vd71i5ciVmxu7duykrK2P27Nmcd955XHjhhaHn0xXNIenSfELcSSdBnz6wejXsbG9Vr4jkmqqqKqZMmQLAqaeeyoYNG1qODRgwgL59+3LllVfy8MMPU1xcnPF8Kgoh6LNtS9fnEwCKioIJZ4AlS9IbTkQi1adPn5bHBQUFNDY2tjwvLCxk8eLFXHDBBTzyyCPMmtXm9b+hUVEIQdmSbswnxGkISSTv7N+/nz179nDuuefys5/9rOWeCv3792ffvn0ZyaCb7ISgW0NHcVqBJBKqVJaQZtq+ffs4//zzOXToEO7OzTffDMCcOXO46qqruOWWW1i4cCGjRo0KLYOKQgjSXhTcU77/q4hkr+QN8b7+9a9/qM3iVv4QnD59eu4vSc1XfbZtoWTjuq7PJ8SNHBncgGf7dti4MX0BRUTaoaKQZi3zCadM6/p8AgQ9Aw0hiUiGafgozTo7dNTWpfYzJw4OisKTTwZF4eKL05ZRJF+5e+g3vo9ad++mqZ5CmqVlPiEuvgJJVzaLdFvfvn3ZsWNHt39pZjN3Z8eOHfTt27fLr6GeQjq9+24wn1Bcwr6Jk7r/elOnBp+XLoXGRujOcJRInquoqKCuro5c2mm5K/r27UtFRUWXv1+/ZdLpz8EdRbt1fUKiQYOCCed162D5cpg8ufuvKZKnioqKqKqqijpG1tPwUTo9+yyQpqGjOF3EJiIZpKKQTmEUBe2YKiIZpKKQLps3w+rV6ZtPiNOyVBHJIBWFdIltXLdn8qnpmU+IO/nkYIJ5+XLYvz99rysi0goVhXRZvhyA/WMnpPd1+/WDSZOguRlefTW9ry0ikiTMO6/daWbbzazdu6mZ2VQzazKz8O8eEabYfib7x4xP/2trXkFEMiTMnsJdQLubgZtZAfBvwFMh5siMWFE4MDqEoqAVSCKSIWHejvM5M6vsoNk/AA8BU8PK0V3tbkMRd/QorFwJwP5RY9MfQpPNIpIhkc0pmNkw4LPAbSm0nWdmNWZWk5VXI65eHRSGqiqai0vS//rjxkH//rBpE2zdmv7XFxGJiXKi+WfAN9y9qaOG7r7A3avdvbq8vDwD0Topvj/6Rz8azusXFLy/5YV6CyISoii3uagG7o/tWDgIONfMGt39kQgzpSxxWGnks68wElg/JMRL6KdNgz/9KSgKs2eH9z4iktciKwru3vIb1MzuAh7LlYKQrHT1KgD2hzHJHKd5BRHJgNCKgpndB8wABplZHfAdoAjA3TucR8glJWuCSeYDYSxHjUssCs3N0EuXmIhI+oW5+mhuJ9peFlaOsPU6dJDiTetpLijgQFV4N9Nm2DD4yEeC7TTWrYPRo8N7LxHJW/pzs5tK1q/Bmps5OGIk3rtPuG8Wn8hesSLc9xGRvKWi0E0la+LzCePCf7MJsS00amvDfy8RyUsqCt1Uujp20VqY8wlxEycGn1UURCQkKgrd1DLJPCoDPQUVBREJmYpCN2W0pxAfPlqxAnrwzcdFJDq6R3M3FOzfR7/NdTQX9ebgCem9cK3NPZeOPx62b4d33oETTkjre4qIqKfQDSVr3wbgwMgx6b2xTns0hCQiIVJR6IbS+MqjTAwdxSUOIYmIpJmKQjfE5xMOZGI5apx6CiISIhWFboivPAp1z6NkKgoiEiIVhW54f+VRBnsKWoEkIiFSUeiiol076PPedhr7FXPoI8Mz98ZDhkBZGezaBdtaX6EkItJVKgpdFN/e4sDocZndsdRMQ0giEhoVhS4qTSwKmRYvClqBJCJppqLQRSWrI5hkjtPGeCISktCKgpndaWbbzeytNo5/0czeiH28aGaTw8oShkiuUYjT8JGIhCTMnsJdwKx2jq8H/srdJwHfBxaEmCW93DO751EyFQURCUloRcHdnwN2tnP8RXffFXv6MlARVpZ0612/jaK9uzk6oIwj5YMzH2D4cCgpCfZA2rEj8+8vIj1WtswpfAV4MuoQqWrpJYweF6wGyjQzbXchIqGIvCiY2dkEReEb7bSZZ2Y1ZlZTX1+fuXBtaFl5lMmL1pJpCElEQhBpUTCzScAdwPnu3uY4iLsvcPdqd68uLy/PXMA2RLK9RTL1FEQkBJEVBTM7AXgYuNTd344qR1dEOskcp56CiIQgtJsAmNl9wAxgkJnVAd8BigDc/TbgRmAg8EsLxuUb3b06rDxp09z8wauZo6KiICIhCK0ouPvcDo5fCVwZ1vuHpe/mdyg82MDhgeUcPXZgdEGqqqBPH6irg717YcCA6LKISI8R+URzrildnQW9BICCAhgXy7ByZbRZRKTHUFHopJL1q4HYctSoaQhJRNJMRaGTStYFRaGhanTESdAeSCKSdioKnVS8fg0AB0aOiTgJ2i1VRNJORaEz3ClZl4VFQT0FEUkTFYVO6L3jPYr27qaxtD9HBh0fdRwYPTqYcF6/Hg4ejDqNiPQAKgqdUBybZD4wckw0ex4l690bxowJ7tW8alXUaUSkB1BR6IT4JPOBbJhkjtMQkoikkYpCJ2TVfEKcViCJSBqpKHRCfOVRQzYVBa1AEpE0UlHohJJ1wb59Gj4SkZ5KRSFFvRoO0G/LuzQXFnGwYkTUcd43Lnajn9Wr4ciRqNOISI5TUUhRycZ1ADScUIkXFUWcJkG/fsHmeE1NQWEQEekGFYUUFa+NbW8xcmzESVqhG+6ISJqoKKSopOUahSyaT4hTURCRNAmtKJjZnWa23czeauO4mdktZrbGzN4ws1PCypIOWXmNQpyKgoikSZg9hbuAWe0c/wwwJvYxD/jPELN0W0k2bYSXTEVBRNIktKLg7s8BO9tpcj7wGw+8DJSZ2dCw8nRLUxPFG2ITzdncU1i1Cpqbo80iIjktyjmFYcA7Cc/rYl/LOv3qNtHr6BEODfkITSWlUcf5sLIyGDIk2BRv48ao04hIDouyKLS2o5y32tBsnpnVmFlNfX19yLE+rGWSORt7CXEaQhKRNIiyKNQBwxOeVwCbW2vo7gvcvdrdq8vLyzMSLlHL9hYqCiLSw0VZFBYBX4qtQvoYsMfdt0SYp00tK4+ycZI5TkVBRNKgMKwXNrP7gBnAIDOrA74DFAG4+23AE8C5wBqgAbg8rCzd1bI7apWKgoj0bJ0qCmZWAhxy96aO2rr73A6OO/D3nXn/SLi/31MYlSNFwT07bgIkIjmn3aJgZr2AOcAXganAYaCPmdUT/KW/wN179IY7RTuz6xacT9dua/XrMycMhQEDYNcu2L4dBg/OcDIR6Qk6mlN4BhgFfBMY4u7D3f144CzgZeDHZnZJyBkj1XLRWtXo7P7r20z3VhCRbuto+Gimux9N/qK77wQeAh4ysyzaMjT9cmKSOW7CBHj55eDeCjNmRJ1GRHJQuz2FeEEws5nJx8zsy4lteqqcmGQmGFZafVwFAO+8sJSna7e1OdQkItKWVJek3mhm/2lmJWY22Mx+D/xNmMGyRXGsp9CQzZPMMfHeTDyziEhnpVoU/gpYCywDngfudfcLQ0uVRXLiauaYA7F7PZSoKIhIF6VaFI4FTiMoDIeBEWbZPOuaJg0N9NtcR3NhYXbdgrMNB4cNp6l3H/pu30rBvr1RxxGRHJRqUXgZeNLdZxEsTf0I8EJoqbLFqlUANJxQlV234GxLQQENlaOA91dNiYh0RqpFYaa73wng7gfd/TrghvBiZYmVKwFoyIWVRzHxC+xK1r4dcRIRyUXtFgUzqwRw903Jx9z9udi+RRXhRMsCsaKQE8tRY+JZNa8gIl3R0XUKN8Wuan4UWArUA32B0cDZwCcJ9jSqCzNkZGIXgeXCJHOcioKIdEe7RcHdLzKziQTbXFwBDAEOAisItrn4obsfCj1lVOLDRzlVFLQCSUS6rsM5BXevBX4A/J6gGKwHlgALe3RBaGqCt4Nx+VwaPmqoHIn36kW/uo30OtxzfzwiEo5UJ5rvBiYAtwC/iD3+TVihssKGDXD4MIcGD83OW3C2wXv34WDFCKy5mX4b10cdR0RyTKpbZ49z98kJz58xs9fDCJQ1cnCSOe7AyDEUb1pPybq3CW5pISKSmlR7Cq/F7o4GgJmdRgrXKZjZLDNbZWZrzOxDS1jN7AQze8bMXjOzN8zs3NSjh2z5ciC35hPiNNksIl2Vak/hNIJbZ8aXpp4ArDCzNwnulzMp+RvMrAC4FfgUweqkJWa2KDZHEfd/gQfc/T9jE9pPAJVd+6ek2ZtvArB/zISIg3SeioKIdFWqRWFWF157GrDG3dcBmNn9wPlAYlFwYEDs8THA5i68TzjiRWFsDhaFUbEVSGtVFESkc1IqCu6+sQuvPQx4J+F5HUGPI9F3gT+Y2T8AJcCHtuiOxNGjLdco7B8zPuIwnRe/rqJ4w9pgFVVBQcSJRCRXpDqn0BWtbZjnSc/nAne5ewVwLvBfsYvlPvhCZvPMrMbMaurr60OImmT1ajhyBCorc2rlUVxT/wEcOn4IBUcOB6uoRERSFGZRqAOGJzyv4MPDQ18BHgBw95cIrpYelPxC7r7A3avdvbq8vDykuAliQ0ecdFL47xWSllVTujWniHRCmEVhCTDGzKrMrDcwB1iU1GYTwVYZmNkEgqKQga5AB3pAUWhQURCRLgitKLh7I3At8BTBldAPuPtyM5tvZrNjzf4JuCp2zcN9wGXunjzElHk9oCiopyAiXZHq6qMucfcnCJaZJn7txoTHtcD0MDN0SU8qCrW17TcUEUkQ5vBRbtq3D9avh6IiGDs26jRdFt8YjxUrIAs6XyKSG1QUksWuZGbChKAw5Kgjg8o5OqAM9u6FLVuijiMiOUJFIVkPGDoCwEzzCiLSaSoKyXpKUSDhauzXe/behSKSPioKyXpQUdg3/sTgwbJl0QYRkZyhopDI/f2iMOlDe/zlnH3jPxo8UFEQkRSpKCTauhV27ICyMhg2LOo03bZ/zHgwC+YUDukubCLSMRWFRIlDR9ba1k25pbm4JFhW29io6xVEJCUqCol60HxCi5NPDj5rCElEUqCikKgnFoUpU4LPKgoikgIVhUQqCiKS50Ld+yiX/O8b7zLjreUUAM8UHk9T7baoI6VHYlFoboZe+jtARNqm3xAx/Tatp+DIYQ4OHUZT/wEdf0OuGDwYhgwJ9nTSDXdEpAMqCjGlq2O338zBezJ3KN5beO21aHOISNZTUYgpfTt+T+YeXBQ0ryAiHQi1KJjZLDNbZWZrzOyGNtpcbGa1ZrbczO4NM097SlevBFQURCS/hTbRbGYFwK3Apwju17zEzBbFbqwTbzMG+CYw3d13mdnxYeXpSEtPoYcNHz1du43i/sM5AzhU8yrPxybQZ04cHG0wEclKYfYUpgFr3H2dux8B7gfOT2pzFXCru+8CcPftIeZp24ED9KvbSHNhIQ2VoyKJEKaG4ZU09ium79bNFO3aEXUcEcliYRaFYcA7Cc/rYl9LNBYYa2YvmNnLZjYrxDxtW74cc6ehajTeu3ckEUJVUMD+cRMBKF2l7S5EpG1hFoXWNg9Kvi9kITAGmAHMBe4ws7IPvZDZPDOrMbOa+vr6tAeNX7TWI+cTYvaNC7bR7r/yrYiTiEg2C7Mo1AHDE55XAJtbafOoux919/XAKoIi8QHuvsDdq929ury8PP1J40Whh80nJNof20a7/8rlEScRkWwWZlFYAowxsyoz6w3MARYltXkEOBvAzAYRDCetCzFT61p6CuMz/taZEr/hjnoKItKe0IqCuzcC1wJPASuAB9x9uZnNN7PZsWZPATvMrBZ4Brje3TM/E5oPPYUx4/FevShet5peh3VvBRFpXah7H7n7E8ATSV+7MeGxA1+LfURj2zaor6extD+HhlZEFiNszf2KOVA5itJ1qylZvRJOHhF1JBHJQrqiOd5LGD2+R9xYpz37JgS7v2peQUTaoqKQB0NHcfs1ryAiHVBReOUV4P2/onuyfVqBJCIdUFF44QUAdp8yNeIg4dvXcgHb8uDeCiIiSfK7KLzzDtTVQVkZB0aOjTpN6I4OLOfQ8UMobDgA6zK/8ldEsl9+F4UXXww+n3563tyRLD6voB1TRaQ1+fGbsC2xoSPOOCPaHBkUn1dQURCR1uR3UYj3FKZPjzZHBu1TT0FE2pG/ReHAgeAXY0EBTO35k8xx+8bHVlmpKIhIK/K3KCxeDE1NMHkylJZGnSZjDg4fQWNxCbz7LoSx46yI5LT8LQp5OHQEQK9e7B+nISQRaV3+FoU8nGSO23vipODBSy9FG0REsk5+FoXm5vd/IeZbTwHYVX168OCZZ6INIiJZJz+LwsqVsHs3VFTA8OEdt+9hdk89Pdj876WX4ODBqOOISBbJz6KQx0NHAEfLjgsm2A8f1hCSiHxAfhaFfJ1kTnT22cFnDSGJSIJQi4KZzTKzVWa2xsxuaKfdhWbmZlYdZp4W8aKQpz0FQEVBRFoVWlEwswLgVuAzwERgrplNbKVdf+A64JWwsnxAfT28/TYUFwdDKPnq4x8P9ntavDi4kE9EhHB7CtOANe6+zt2PAPcD57fS7vvAvwOZuXFwfAx92jQoKsrIW2alY46BU06Bo0ffn2MRkbwXZlEYBryT8Lwu9rUWZnYyMNzdH2vvhcxsnpnVmFlNfXevwtXQ0fs0hCQiScIsCq3d8NhbDpr1Am4G/qmjF3L3Be5e7e7V5eXl3UuV5yuPPkBFQUSShFkU6oDEiwAqgM0Jz/sDHwWeNbMNwMeARaFONh85AkuWBI9PPz20t8kZZ54ZbAhYUwP79kWdRkSyQJhFYQkwxsyqzKw3MAdYFD/o7nvcfZC7V7p7JfAyMNvda0JL9Oqrwdr8CRPguONCe5uc0b9/sENsUxP85S9RpxGRLBBaUXD3RuBa4ClgBfCAuy83s/lmNjus922Xrk9o8XTtNp6u3cb6k6YBsGHhYzxduy3iVCIStcIwX9zdnwCeSPrajW20nRFmFkCTzK3YNW06VbffwnGLtQJJRPLpimZ3TTK3YvfJU2kuLKL/irco3Lsn6jgiErH8KQobNsDWrTBwIIwdG3WarNHcr5g9k0/Fmpspq3k56jgiErH8KQqJQ0fW2mrZ/LVrWtBz0hCSiORPUdDQUZt2TQsm3o9VURDJe6FONGeVs8+GPXvgk5+MOknW2TP5VJp696H/quXw3nswaFDUkUQkIvnTU7joIrjnnmBdvnxAc5++7JkSu2bwz3+ONoyIRCp/ioK0Kz6EpC0vRPKbioIAsFNFQURQUZCYvSedTFO/flBbC9t0ZbNIvlJREAC8d292nxxseaHegkj+UlGQFjumx7bSvv/+aIOISGRUFKTF1vMugMJCeOwx2LIl6jgiEgEVBWlxZFA5zJ4dbKV9991RxxGRCKgoyAddeWXw+Y47gk0ERSSvqCjIB51zDlRUwNq1upBNJA+Fus2Fmc0Cfg4UAHe4+4+Tjn8NuBJoBOqBK9x9Y1h5dBOZFBQUwBVXwPz5QW9hxoyoE4lIBoXWUzCzAuBW4DPARGCumU1MavYaUO3uk4CFwL+HlUc64fLLg51kFy6EXbuiTiMiGRTm8NE0YI27r3P3I8D9wPmJDdz9GXdviD19GagIMY+kqrISPvWp4H7Wv/1t1GlEJIPCLArDgHcSntfFvtaWrwBPhphHOuOqq4LPt9+uCWeRPBJmUWjtTjat/nYxs0uAauCmNo7PM7MaM6upr69PY0Rp0+zZwRbab74JNTVRpxGRDAlzorkOGJ7wvALYnNzIzGYC/wL8lbsfbu2F3H0BsACgurpaf7aGKHEyfsx5FzLirtuou+kXrPxuq/UagJkTB2cimohkQJg9hSXAGDOrMrPewBxgUWIDMzsZ+BUw2923h5hFuuDdz30BgCGPP0zBgQMRpxGRTAitKLh7I3At8BSwAnjA3Zeb2Xwzmx1rdhNQCjxoZsvMbFEbLycRaBg1lt0nT6Ow4QDH/0E/GpF8EOp1Cu7+BPBE0tduTHg8M8z3l+5794IvUPbaYoYtvJctn50bdRwRCZmuaJZ2bfv039BYUkrZsiWUrFkVdRwRCZmKgrSrubiErX/9OQAqF/w84jQiEjYVBenQxsuvoalvP4Y+/jDlTz/R8TeISM5SUZAOHTyhkjVf+xcAJnzveop2vhdxIhEJS6gTzdJzvDP3CsqffpLjFr/A+Pnf4M2b7wj2R2pHWxsQ6roGkeylnoKkplcvan9wM43FJQz+4+MMfuKRqBOJSAhUFCRlh4adwNvfmA/A+B98k9712opcpKfR8JF0yuYLvsDxTz/OoL/8iQnf+Sdev/W/dJ8KkR5EPQXpHDNWfO+nHB1wDOV/fpqhj9wfdSIRSSMVBem0w4OHsupbPwRg3I9vpM/muogTiUi6qChIl2w97wK2f/IzFO7fR/Vln2XAm69FHUlE0kBzCtI1Zqz47k302bqZY5a/TvWl5/P2P3+XurmXd7hUtbO0tFUkc8xz7K5a1dXVXtPFm75oQjT97Mhhxt70PYbfeycAWz/7TqC1AAAHA0lEQVRzPiu+91OaSko7/Vpt/ZLv7M9NxULkw8xsqbtXd9ROw0fSLd67D6v+5Ue8+ZPbaCwuYciTjzLt4k9Tuqo26mgi0gXqKUjaFG9Yy0n/eCX9315BU5++bLzsarb+zQU0VI2OOhqgHoTkt6zoKZjZLDNbZWZrzOyGVo73MbP/jh1/xcwqw8wj4WqoHMWSex/n3c/NpeDwIUb+6mbOOO9Mpl18DifcdRt9tm2JOqKIdCC0noKZFQBvA58iuF/zEmCuu9cmtLkGmOTuf2dmc4DPuvvn23td9RRyw7EvP8/Q3z/I8X98nMID+wFwM3ZNPYNdp51Jw4gqGkaMpGHEyC7NP3SFegqSz1LtKYRZFE4Hvuvun449/yaAu/9rQpunYm1eMrNCYCtQ7u2EUlHILb0OHWTQc//LkMd/x8DnnqbgyOEPtTk8sJyGylEcLj+eppL+NJaW0ljaP/a4P019+uJFhXhhEV5QQHNhEV5YiBcUgBluvYIVT7164WYfXv2U8NzbWhnVytdPGzmwzX/XK+t2tPr19r4nm+R6/p6uzZ/PSSfA+PFdes1Ui0KYS1KHAe8kPK8DTmurjbs3mtkeYCCgvZl7iOa+/dh+znlsP+c8CvbtpfzZP1C6qpbiTesp3riOfps20GdHPX121EcdtVOS/0PONbmev6dr8+dz2mnw8suhvneYRaG1P8mSewCptMHM5gHzYk/3m1kU94UcRG4Uq1zImQsZQTnTKRcyQrbnfOWVeK+2KzlHpNIozKJQBwxPeF4BbG6jTV1s+OgYYGfyC7n7AmBBSDlTYmY1qXS9opYLOXMhIyhnOuVCRlBOCHf10RJgjJlVmVlvYA6wKKnNIuDLsccXAn9qbz5BRETCFVpPITZHcC3wFFAA3Onuy81sPlDj7ouAXwP/ZWZrCHoIc8LKIyIiHQt17yN3fwJ4IulrNyY8PgRcFGaGNIp0+KoTciFnLmQE5UynXMgIypl7VzSLiEh4tPeRiIi0UFFIYGZ3mtl2M3urjeNmZrfEtuV4w8xOyXTGWI6Ocs4wsz1mtiz2cWNr7ULOONzMnjGzFWa23My+2kqbyM9nijmz4Xz2NbPFZvZ6LOf3WmkT6bYxKWa8zMzqE87llZnMmJSlwMxeM7PHWjmWFVvwdJAxnHPp7vqIfQAfB04B3mrj+LnAkwTXV3wMeCVLc84AHov4XA4FTok97k+w5cnEbDufKebMhvNpQGnscRHwCvCxpDbXALfFHs8B/jsLM14G/EeU5zIhy9eAe1v72UZ9LlPMGMq5VE8hgbs/RyvXSSQ4H/iNB14GysxsaGbSvS+FnJFz9y3u/mrs8T5gBcEV7IkiP58p5oxc7Bztjz0tin0kTwieD9wde7wQ+KRZmu941I4UM2YFM6sA/hq4o40mkZ5LSCljKFQUOqe1rTuy7hdIzOmxbvyTZnZilEFiXe+TCf5yTJRV57OdnJAF5zM2lLAM2A780d3bPJ/u3gjEt43JpowAF8SGCxea2fBWjmfCz4B/BprbOB75uaTjjBDCuVRR6JyUtuXIAq8CI9x9MvAL4JGogphZKfAQ8H/cfW/y4Va+JZLz2UHOrDif7t7k7lMIdgeYZmYfTWoS+flMIePvgUp3nwQ8zft/jWeMmZ0HbHf3pe01a+VrGTuXKWYM5VyqKHROKlt3RM7d98a78R5cK1JkZoMyncPMigh+0d7j7g+30iQrzmdHObPlfCbk2Q08C8xKOtRyPq2dbWMyoa2M7r7D3eNb5d4OnJrhaADTgdlmtgG4H/iEmf02qU3U57LDjGGdSxWFzlkEfCm2auZjwB53z7o7x5jZkPj4p5lNI/g5t74Xb3gZjOCK9RXu/v/aaBb5+UwlZ5acz3IzK4s97gfMBFYmNYt025hUMibNGc0mmMPJKHf/prtXuHslwSTyn9z9kqRmkZ7LVDKGdS5DvaI515jZfQQrTQaZWR3wHYLJMtz9NoKrs88F1gANwOVZmvNC4GozawQOAnMy+R90zHTgUuDN2BgzwLeAExJyZsP5TCVnNpzPocDdFty8qhfwgLs/Ztm1bUwqGa8zs9lAYyzjZRnO2KYsO5etysS51BXNIiLSQsNHIiLSQkVBRERaqCiIiEgLFQUREWmhoiAiIi1UFEREpIWKgoiItFBREOkmM5sa25Ssr5mVxO4lkLznj0hO0MVrImlgZj8A+gL9gDp3/9eII4l0iYqCSBqYWW9gCXAIOMPdmyKOJNIlGj4SSY/jgFKCu7f1jTiLSJeppyCSBma2iGCL4ypgqLtfG3EkkS7RLqki3WRmXwIa3f3e2A6hL5rZJ9z9T1FnE+ks9RRERKSF5hRERKSFioKIiLRQURARkRYqCiIi0kJFQUREWqgoiIhICxUFERFpoaIgIiIt/j/lCjShQVlSlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('x')\n",
    "\n",
    "q_of_x = 1./x\n",
    "\n",
    "pred_mean_q = 1./mean_x\n",
    "pred_std_q = np.sqrt((std_x/mean_x)**2)/mean_x\n",
    "\n",
    "counts, bins, patches = plt.hist(q_of_x, bins=50, density=True, alpha=0.3)\n",
    "\n",
    "plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)\n",
    "plt.legend(('pred','hist'))\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('p(x)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now let's do the same thing with an interactive widget!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_1_over_x(mean_x, std_x, N):\n",
    "    x = np.random.normal(mean_x, std_x, N)\n",
    "\n",
    "    q_of_x = 1./x\n",
    "\n",
    "    pred_mean_q = 1./mean_x\n",
    "    pred_std_q = np.sqrt((std_x/mean_x)**2)/mean_x\n",
    "\n",
    "    counts, bins, patches = plt.hist(q_of_x, \n",
    "                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), \n",
    "                                     density=True, alpha=0.3)\n",
    "\n",
    "    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)\n",
    "    plt.legend(('pred','hist'))\n",
    "    plt.xlabel('q(x)')\n",
    "    plt.ylabel('p(x)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "09cc77921c8f4f4287866c77cc69629a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_1_over_x(mean_x, std_x, N)>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now make the interactive widget\n",
    "interact(plot_1_over_x,mean_x=(0.,3.,.1), std_x=(.0, 2., .1), N=(0,10000,1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check propagation of errors under division "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b158140c73a74d55947cdef181af2f8c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, description='mean_x', max=3.0), FloatSlider(value=1.0, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_x_over_y(mean_x, std_x, mean_y, std_y, N)>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def plot_x_over_y(mean_x, std_x, mean_y, std_y, N):\n",
    "    x = np.random.normal(mean_x, std_x, N)\n",
    "    y = np.random.normal(mean_y, std_y, N)\n",
    "\n",
    "    q_of_x_y = x/y\n",
    "\n",
    "    pred_mean_q = mean_x/mean_y\n",
    "    pred_std_q = np.sqrt((std_x/mean_x)**2+(std_y/mean_y)**2)*mean_x/mean_y\n",
    "\n",
    "    counts, bins, patches = plt.hist(q_of_x_y, \n",
    "                                     bins=np.linspace(pred_mean_q-3*pred_std_q,pred_mean_q+3*pred_std_q,30), \n",
    "                                     density=True, alpha=0.3)\n",
    "\n",
    "\n",
    "    plt.plot(bins, scs.norm.pdf(bins, pred_mean_q, pred_std_q), c='r', lw=2)\n",
    "    plt.legend(('pred','hist'))\n",
    "    plt.xlabel('q(x)')\n",
    "    plt.ylabel('p(x)')\n",
    "\n",
    "\n",
    "\n",
    "interact(plot_x_over_y,mean_x=(0.,3.,.1), std_x=(.0, 2., .1), mean_y=(0.,3.,.1), std_y=(.0, 2., .1),N=(0,100000,1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}