SeddonStettler_HW05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ASTR 211 - HW05 - Sadie Seddon-Stettler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note: As encouraged, Cole Meldorf, Kathryn Downey, and I all worked together collaboratively on this homework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import emcee\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from codes.cosmology import d_L_vectorized\n",
    "from codes.constants import clight\n",
    "from codes.plot_utils import plot_pretty, plot_2d_dist\n",
    "plot_pretty(fontsize=12)\n",
    "\n",
    "from scipy.special import gamma\n",
    "from scipy.interpolate import RectBivariateSpline\n",
    "from scipy.optimize import minimize, differential_evolution\n",
    "\n",
    "from timeit import default_timer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read sample of 740 supernovae...\n"
     ]
    }
   ],
   "source": [
    "zCMB, mB, emB, x1, ex1, csn, ecsn = np.loadtxt('data/jla_lcparams.txt', \n",
    "                                               usecols=(1, 4, 5, 6, 7, 8, 9), unpack=True)\n",
    "\n",
    "print(\"read sample of %d supernovae...\"%(np.size(zCMB)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished in 37 sec\n"
     ]
    }
   ],
   "source": [
    "Om0min, Om0max = 0., 1.\n",
    "OmLmin, OmLmax = 0., 1.\n",
    "\n",
    "nOm0 = 13; nOmL = 9\n",
    "Om0tr = Om0min + (Om0max-Om0min) * np.linspace(0.,1.,nOm0)**2 # square-law \n",
    "OmLtr = OmLmin + (OmLmax-OmLmin) * (1.- np.linspace(0.,1.,nOmL)**1.75)\n",
    "# the values of OmLtr are in descending order, so reverse the order before use\n",
    "OmLtr = OmLtr[::-1]\n",
    "\n",
    "nSN = np.size(zCMB)\n",
    "\n",
    "dLgrid = np.zeros((nSN, nOm0, nOmL))\n",
    "H0 = 70.\n",
    "H0c = H0 / clight\n",
    "\n",
    "tstart = default_timer()\n",
    "atol = 5.e-16; rtol = 5e-16\n",
    "for i, Omd in enumerate(Om0tr):\n",
    "    for j, Omld in enumerate(OmLtr):\n",
    "        # cancel the c/H0 factor to tabulate d_L without it\n",
    "        dLgrid[:,i,j] = H0c * d_L_vectorized(zCMB, H0, Omd, Omld, atol=atol, rtol=rtol)\n",
    "\n",
    "print(\"finished in %.3g sec\"%(default_timer()-tstart))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished in 0.144 sec\n"
     ]
    }
   ],
   "source": [
    "dLz = []\n",
    "tstart = default_timer()\n",
    "for iz, zi in enumerate(zCMB):\n",
    "    #construct approximation for tilde(d)_L(z_i, Om0, OmL)\n",
    "    spl2d = RectBivariateSpline(Om0tr, OmLtr, dLgrid[iz], s=0, kx=3, ky=3)\n",
    "    dLz.append(spl2d)\n",
    "\n",
    "print(\"finished in %.3g sec\"%(default_timer()-tstart))\n",
    "dLz = np.asarray(dLz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2(x):\n",
    "    mu_model = np.zeros(nSN)\n",
    "    for i in range(nSN): \n",
    "        mu_model[i] = 5.*np.log10(dLz[i](x[0], x[1])) \n",
    "\n",
    "    delta_mu_2 = (mB - (mu_model + x[2]) + x[3]*x1 - x[4]*csn)**2\n",
    "    sig_2 = (emB**2) + (x[3]**2)*(ex1**2) + (x[4]**2)*(ecsn**2)\n",
    "    return np.sum(delta_mu_2/sig_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ln_L(x):\n",
    "    chi_2 = chi2(x)\n",
    "    sig_2 = (emB**2) + (x[3]**2)*(ex1**2) + (x[4]**2)*(ecsn**2)\n",
    "    ln = np.sum(np.log(2*np.pi*sig_2))\n",
    "    return (-1/2)*(chi_2 + ln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neg_2_lnL(x):\n",
    "    return -2*ln_L(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "completed in 21 sec\n",
      "minimum at : [ 0.25748662  0.63500693 24.07861696  0.11942862  2.57935328]\n",
      "-2lnL min = 3.3830e+02\n"
     ]
    }
   ],
   "source": [
    "tstart = default_timer()\n",
    "b = [(0, 1), (0, 1), (20, 28), (0.05, 0.3), (1., 5.)]\n",
    "res = differential_evolution(neg_2_lnL, bounds=b)\n",
    "print(\"completed in %.3g sec\"%(default_timer() - tstart))\n",
    "print(\"minimum at :\",res.x)\n",
    "\n",
    "ln_Lmin = ln_L(res.x)\n",
    "print(\"-2lnL min = %.4e\"%ln_Lmin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously, our derived value for Omega0 was 0.394, and our tM0 value was 24.1. These are both very similar to the new values (of 0.257 and 24.07 respectively). However, our new Omega0 is 0.635, quite different from our previously derived value of 0.363."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reduced chi2 = 0.9316979440494364\n"
     ]
    }
   ],
   "source": [
    "x = [0.25752049, 0.63500829, 24.07862607, 0.11942735, 2.57935899]\n",
    "red_chi2 = chi2(x) / (nSN - 3)\n",
    "print(\"reduced chi2 =\", red_chi2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this $chi^2$ value is very close to 1, this suggests that our model describes the supernova measurements well - almost to the point of overfitting (if it were below 0.8, I would start to worry about overfitting), but not there. The $chi^2$ value of our previous model was around 6, so this suggests that the 5-parameter model describes the measurements much better than the 3-parameter model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ln_prior_pdf(x, args = None):\n",
    "    px0, px1, px2, px3, px4 = (0, 0, 0, 0, 0)\n",
    "    if 0 <= x[0] <= 1:\n",
    "        px0 = 1\n",
    "    if 0 <= x[1] <= 1:\n",
    "        px1 = 1\n",
    "    if 20 <= x[2] <= 28:\n",
    "        px2 = .125\n",
    "    if 0.05 <= x[3] <= 0.3:\n",
    "        px3 = 4\n",
    "    if 1 <= x[4] <= 5:\n",
    "        px4 = .25\n",
    "    if px0*px1*px2*px3*px4 == 0:\n",
    "        return -3000\n",
    "    else:\n",
    "        return np.log(px0*px1*px2*px3*px4)\n",
    "\n",
    "def posterior_pdf(x, args = None):\n",
    "    return ln_prior_pdf(x) + ln_L(x)\n",
    "\n",
    "def neg_post_pdf(x):\n",
    "    return posterior_pdf(x)*(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vectorized_mcmc(x, nsteps=1, step=0.1, modelpdf = None, args = None):\n",
    "    \"\"\"\n",
    "    MCMC sampler implementing a simple Metropolis algorithm\n",
    "    to follow a number of chains (\"walkers\") in parallel\n",
    "    \n",
    "    This version also can sample distributions of arbitrary number of \n",
    "    dimensions ndim\n",
    "    \n",
    "    Parameters:\n",
    "    ------------\n",
    "    x0 - a real numpy array of size [nwalkers,ndim]\n",
    "        initial x value\n",
    "    nsteps - integer\n",
    "        number of MCMC steps (iterations) to take\n",
    "    step - float\n",
    "        step size controlling step proposal distribution\n",
    "    modelpdf - python function object\n",
    "        ln(pdf) where pdf is the target pdf to sample\n",
    "    args - pointer to a list\n",
    "        list of arguments to pass modelpdf\n",
    "        \n",
    "    Returns:\n",
    "    ---------\n",
    "    xchain - numpy array of size [nsteps*nwalkers, ndim]\n",
    "        coordinates of samples in the MCMC chains of size\n",
    "    pacc   - float\n",
    "        acceptance ratio of the MCMC steps (ratio of accepted to the total proposed number of steps)\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # the input array here contains initial values for multiple MCMC sequences\n",
    "    # or \"walkers\" in the MCMC jargon (because they \"walk\" the multi-d space of our target pdf)\n",
    "    nwalkers = np.shape(x)[0]\n",
    "    # make sure input is sensible\n",
    "    assert(nwalkers>0)\n",
    "    \n",
    "    # number of dimensions that we will be sampling\n",
    "    ndim = np.shape(x)[1]\n",
    "        \n",
    "    # initialize some auxiliary arrays and variables \n",
    "    chain = x;\n",
    "\n",
    "    naccept = 0; ntry = 0; nchain = 0\n",
    "    # initialize arrays that will be used to hold pdf values of samples\n",
    "    # at old and proposed sample locations\n",
    "    gxold = np.empty(nwalkers)\n",
    "    gxtry = np.empty(nwalkers)\n",
    "    for i in range(nwalkers):\n",
    "        gxold[i] = modelpdf(x[i,:], *args)\n",
    "    eps = 2.e-16\n",
    "    \n",
    "    nsample = 0 \n",
    "    while nsample < nsteps:\n",
    "        # proposal step using uniform pdf in range [-step, step]\n",
    "        xtry  = x + step*np.random.uniform(-step, step, np.shape(x)) \n",
    "        for i in range(nwalkers):\n",
    "            gxtry[i] = modelpdf(xtry[i,:], *args) \n",
    "        gx = np.copy(gxold) \n",
    "        # compare pdf values at the old and proposed sample locations  \n",
    "        gr   = gxtry - gx\n",
    "        u = np.random.uniform(eps,1.0,np.shape(x)[0])\n",
    "        # accept proposal with probability min[1.0, e^gr]\n",
    "        iacc = np.where(gr > np.log(u))\n",
    "        # update those \"walkers\" for which proposal step was accepted \n",
    "        x[iacc,:] = xtry[iacc,:]\n",
    "        gxold[iacc] = gxtry[iacc]\n",
    "        naccept += np.size(gxtry[iacc])\n",
    "        # add all walkers to the chain, regardless of whether their proposal step was accepted\n",
    "        chain = np.vstack((chain, x))\n",
    "        nsample += 1\n",
    "\n",
    "    return chain, 1.*naccept/(nsample*nwalkers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: -336.2202456224824\n",
      "     jac: array([ 1.25055521e-04, -6.82121026e-05, -6.82121026e-05, -1.70530257e-03,\n",
      "       -1.87583282e-04])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 1413\n",
      "     nit: 16\n",
      " success: True\n",
      "       x: array([ 0.25751647,  0.63500183, 24.07862665,  0.11942716,  2.57935803])\n"
     ]
    }
   ],
   "source": [
    "snbounds=[(0,1), (0,1), (20, 28), (0.05, 0.3), (1, 5)]\n",
    "\n",
    "x0 = differential_evolution(neg_post_pdf, bounds=snbounds)\n",
    "\n",
    "print(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emcee time = 98.8 sec, acceptance fraction = 0.534\n"
     ]
    }
   ],
   "source": [
    "ndim = 5\n",
    "nwalkers = 10\n",
    "nsteps = 1000\n",
    "args0 = []\n",
    "\n",
    "#differential evolution yielded [ 0.25747937,  0.63495553, 24.07862841,  0.11942547,  2.5793592 ]\n",
    "x0 = np.zeros((nwalkers, ndim))\n",
    "for i in range(nwalkers):\n",
    "    x0[i][0] = np.random.normal(0,scale=.01) + 0.25747937\n",
    "    x0[i][1] = np.random.normal(0,scale=.01) + 0.63495553\n",
    "    x0[i][2] = np.random.normal(0,scale=1) + 24.07862841\n",
    "    x0[i][3] = np.random.normal(0,scale=0.01) + 0.11942547\n",
    "    x0[i][4] = np.random.normal(0,scale=0.1) + 2.5793592\n",
    "\n",
    "tstart = default_timer()\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, posterior_pdf, args=[])\n",
    "echain, prob, state = sampler.run_mcmc(x0, nsteps)\n",
    "tend = default_timer()\n",
    "# acceptance fraction averaged over all walkers\n",
    "eacc = np.mean(sampler.acceptance_fraction)\n",
    "print(\"emcee time = %.3g sec, acceptance fraction = %.3g\"%(tend-tstart, eacc))\n",
    "x_1 = sampler.flatchain[:,0];\n",
    "x_2 = sampler.flatchain[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sadie\\Anaconda3\\lib\\site-packages\\scipy\\signal\\signaltools.py:974: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  res *= (1 - noise / lVar)\n",
      "C:\\Users\\sadie\\Anaconda3\\lib\\site-packages\\scipy\\signal\\signaltools.py:974: RuntimeWarning: invalid value encountered in multiply\n",
      "  res *= (1 - noise / lVar)\n",
      "C:\\Users\\sadie\\Documents\\GitHub\\a211w20\\codes\\plot_utils.py:829: UserWarning: Log scale: values of z <= 0 have been masked\n",
      "  norm = LogNorm(), extent = [xbins[0], xbins[-1], ybins[0], ybins[-1]])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOQAAADWCAYAAADB5S4rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dd3yURfrAv7ObZLNJSE9IgBASSOiEEqQqqKCAir1g7556tlPvvNPzmp71iuXO36HooSKWUxFRAUFQem8htBQChIT0upvNZnd+f2wSl7Cb7S3s9/PJJ7vvvu/MM/POM+WZZ2aElJIgQYL4BwpfCxAkSJCfCSpkkCB+RFAhgwTxI4IKGSSIHxFUyCBB/IigQgYJ4kd4TSGFEO8KISqEEHlWfhdCiNeFEAVCiL1CiLHeki1IEH/Bmy3kf4FZ3fw+G8hq/7sXeMsLMgUJ4ld4TSGllD8BNd3ccjnwvjSxGYgVQqR6R7ogQfwDfxpD9gWOm30/0X4tSJCzhhBfC2CGsHDNql9fVFSUNHf7S0hMJCkpyS2RuhNXHROdlc883katvvNzL3WoxesA0Wa/OYs7HTG7pt0bTp628tseGXbu2FElpXS8MOJfCnkCSDP73g84ae3mwUOGsGHLVpcjDVda7iS0GIzd/m4vHeF4iyUbjwJwxeQBndfW55V1fp46ItXidYAZOa53SDyV3nClwit56Y541CHKEmef9SeFXAr8UgjxMTABqJdSltl4xmFcVTBPxOeJgtahmADNGj1tGj0xiRGszyvrVEpz5eyQw9X8sfW8tysoR/G1fF5TSCHEYmA6kCiEOAH8AQgFkFL+H/AtMAcoADTAHfaG3bUQdJep9ma4txXXHhxJpyWaNaYuqrlS2oM7CmmH7I7ma0fcvlYU8E4r7TWFlFLOs/G7BB50JmxPZJKtMH3R8nkinR3p8HTB9weFsoU/VML+1GV1CEeMHf6Q0d6ioxXs4KYZWSzZeJR6jZ6WWi0AidGqzlbSna2fOYGggM7g6VYyYBUS7Fc0RzLQncrrahfTVW6akXXa9xBVCJERoXZ1V31RiXk7v7pLo68qlIBVSMmZmeaJQmQehzsMHt580YtWHTlDKcFk8EmMVp12zZExpSWcqfQs5a0vFMGdcbpaRgJWIS3hjnGfv/KfL7ZSsH8nYybP4MYLB50xZWFJoTqUseZUAwd2bWbwiAlnKGIHXSsLTyqGpbB7ahfXUXqUQtrCnpfeVWn9RYlf/8P9FB3Yzd1PvEKE+qrOsWJ6ShQl5U2UlB/p/H6koAZlmJL3vz0IwGf/fYX1yz/jmrueJCzqJjL7x1iMo7uWq4OuUyPOKJK58nu712AepydaaFef94/S1sPpePmuvKzJF15On/6DGDh0jM17C/J3sOzjf2MwtHHrnCGUFOwnLCSco0csLrRxCEsVlvmfPXirJbYnfn8jYFtIgXtqNX9pAW3x9j+fhX8+2znhHxlhcnObOiKVqoajnS1mVYOOj/7zZ8pLi8nIGMWShCjm3f07vl+6kCtueYRb5ww5I2xbxo2zwYpqyR5hLY2eLDMBq5CWjDq28LbyeaLQmrvEdVB8oBCj0UBi734kpkTxiydfp6GumgGDRgKQOTKH+0b+vVOJLcnpS6OKP2LLucRTZSlgFdIZXPHSsde31dMF2twlDiB/z25e/8P9GA0GLrzsFi654RaTYvbu57Y4nR0ndhdOoFcAnpL7rFJIW5gP9q0Zdzp+s/RCrCmrq0aQDqSUfPDPV1CFRzBtxg0YhJ6Xn7yVi0ddR37DXpZ8+E+KD+/l/qdeI6S9NTR3FGjW6Fm06gh3XTzYahxdFcXZuUFr+eNvCuhIS+cN2YMK6QTWXoyj1x3l1QVf8cPXC1GEhZKYnMYP375P79QBLN36HpOGzaIiNoGpM69mYGZc5zTIolVHzghnwYpDZ8xPWmux3F0ILVV0vlZSd02XuaMbGxgWDS/R1frnSCvoafR6Pe/9/Y+06LWEhaspLthD4eFd6HTNgGBT/nK0mmYysnOsTvK31Gi8K3TX+M2sze6wPLtLJnvu8dYcd8C2kB1WVk+/UHeMO23da+2+JRuP8sHrLxESHkZNZRnHiw6SNXwck2ZcSmioisHHz6GhpByFQkFErxiGjpxIxckS1u1LR4juvX27evH4WjGCmAhYhexqZbXW5fJll6hr3N2NMS1NlFeXNvDt/+YzecRsGsoraGvTU1NRTvGhPNYt/5zIXrFEhEZhNBrITh/Lri2r2L5hOWu+u5jbH3ueqJheRESEcmW7Zda8+2rJpc4T+Mpf1FnDnKPOI+5OQ8AqpDUsKaA3nJbt7fo4ek+4OopWRRuHGvYzeOQ5jD5nJp8seB6A5sY6mqljQNpwdu5dRfaI8bQ1tFC8YxdLF73JrQ/9DiklX248SlK0ivSUKJcdy/2pJfUHB3h3V/g9RiHttQp2eJT4U8Hqirmf6sPP/Ieln76J0WCg8OAeDu3bSlJKGm1tOmqrKpidexPbitYCEJvUG1VfNWHlkUycNpen7riEhroa7v3NC1x0cXc7cNqPr5XA3uknW+F0xV/8oHuMQnZgK2P9WRHNKa2oo+jAHt557XEa6qsJV6q5aeSDbGneSPHJfPStOgC0ITqqakoJU4RzcNcmHn9uIeUnD/DXX10PBklEVAyb1yxnYNYEwPqqDl8rWgf2Op53HQp425ZgreIPrvbogej1rfz+vms5XnyAc3IupkJzkqNH9rHkyIdUN5Yz+cIryd+1ntuv/iNZEYMYMugcQo638Omhd1izbDEFB3Ywdfq1rP/xMy669nYuuuY2DFrJbXOHWY3TlgHKncvQbMXpq3DcEaersgSsQnp6+0ZvM3/Rrs7PxUV7aNO3MiZpAskksef4OgCMbW0AbFz9JREhkcz/6Dedz/z+6vk0bqtl7fJFTJo+lyOHt3PhJbdw3oXzECKEkAjr6yMdxZ+mg3oaAauQntqj05vjy46xopSS776cz6G8zYSGhVN56ji11WVExIaR1TcJ3R4NU8ddwfodSwCYkHkhW4pWA3DzyF8SEhOOetogdEs0JCT15eKr7iKlXyYhEaF2K6Cn0+3PQwVHHcm7S0uwy+pm7MlsV30yO/xRO9zaokO1rFy6gHvv/welJw5xOH8rj9zwdw7pD1J18iSPznyBBWtfBKBv/2wUQoECBb1j0/hw35u8fNUiGpsbMBoMpPYbiDo8ultltDT+sZYmR7qy3sSdlmB/qiyC/Qw/YNfuvURGxZIxdDTTZ9/MC2+soa1BS1neQSp05YT0iuD83KsBuOr8B4hXJZHSO4OyuhKS+6Sj6BNFjDaCF/6xhqioeD5Z8EK38bnDS8aR591l8LAUvyc9f7qbN7a1wbazBBWyHWsLbc2vWXrhrhaCiuLjvPHHR5l30zMo6/RIvQGVSk16zEASY1MJLzOwvvoHNhSuYGDWWCrLS1CP6oder+OPv1/CE08sZGPlOtbv+4bQBgMXzLqVosO7aanVsuCzvWeky5E5WUv3W8ob83zoTkGczSdfjU+7S1NwtYcPcXfmXzF5AOvzyqgsO8Hv7p/JkBGTyE4ZCSebMDS2YKzQom/QElmlYHXZekZHnEt1RSnVFaXcMvx+ojLTmPjsZTRWVLB40fMcObwdpTKEvPz1TJp0OQ311ShCFcQlR7k9vV0VzBvKYmmFTde5SXvkcHQ44gu8qpBCiFnAa4ASeEdK+WKX3/sDC4HY9nueklJ+6w3ZPL1SwxLqyCgUSiWzc+bRWlnHkYPbaUg1EHPQSFHlASpiG1CHRKJR67h61kPsP76dqCH9oEaLlJIFi54ma/g4Xvjd9+gbW9mzbQ1rly/mvBnXo1BYXi7maCtpL57qwtkTl7OrMZxd2dIjdgwQQiiBfwEzMR2ss00IsVRKmW922zPAp1LKt4QQwzAdLzDAWzJ6m2ZNK9ExiWwoXon+kI5d279HKZTcMfox5ky4hdW6VeTvXY+xUcd3Py5Eo21AXiqhzciKnz7AYGwjPCKCBa/8hhf/9S733/wUi1Zd3Rl+x+4ClrqN7ihU9njQuMP44u6xpyuyeBpvtpDnAAVSyiKA9kN1LgfMFVIC0e2fY+jm9Ctf445pgvnPP07O2AvIjZ2ARtVKeHQEOzd+T6mijGwg42giieoUCkr2EBIaSt8BWfxrwWNcM+IuNuz6mjCVmu8+fxsAaTSesTVkdzjrKmZLoS11H13NJ2dbMk/hyXlYbyqkpQNZJ3S554/ASiHEQ0AkMMM7ojmHrZ0DLBXOO594EWOLkfT0ERzK28pNc+9lQ+WP7N37I5X1J9EZWygq2UfF8WLaDHqqtOUohJI+SZmUnTyKvlXHq2t/DYBW0wjA48/9l5BQ0w4BkRGhVvdedQRL4zbztFlLb9cwHMFfxnG+xJsKac+BrPOA/0op/yaEmAR8IIQYIaU84w1VVlYyZcI5nd/vvPtu7rrnXrcK3B2OGD46KCkp4b2//RaA2PgUzh91Jc8suYP7xj3F0GkP8dGKVxmvyqWXMoadhr1c+fgzDPjsG/bsX0PJycNIs+waOHgshYd2AlBd/nNHIjFa5fIu5N2l0ZOGEWcm5/1loYC7ZPCmPdmeA1nvAj4FkFJuAsKBREuBJSUlsWHL1s6/rsroj65c6enpXHfPb4nsFcvIpHFMzL4IKSUL89/kaMRJpoydS7HyBB9XLKJJU0f+R0uJH5hOn94DO5WxX2oWF0+7lbqaU8RExDNs9GTGnXv6Sg57uq6WxnrdTY/4mq7TDY7KaGk6y1ncFY4lvJnr24AsIUSGECIMuAHTIa3mHAMuBBBCDMWkkJXOROYPtWZX/vPZWspLC9E3t7Du0DecOHEIhVCSljGEzz96hW+2vo+xXsOvBzyDCFVScvIg9aWlXP3M8zw496+oI3oxvF8uK358n+rKUpIGZHDfU/8gTO3aOR22JsDt/bMVnrPY6grbo5zumhO1hDsrMm+eD9kmhPglsALTlMa7Usr9Qog/A9ullEuBx4G3hRCPYerO3t5+bmRA8smPBXyzeD4ZQ0ZTsH8nK/73HiFtSoxKAxl9RvBT1feEqsLp238w1fmFVOtrSJmcw/YTu8iIGkRrFEzMvgj9Z/mkx2Vw98ynKVacIGR3GLfe/xfOnTOb8IhIwL7jyJ0pkPbc19240tGw7MXby6/M47Imhzvw6jxk+5zit12uPWv2OR+Y4k2ZPMm65Z+z8oP3aGprwGBoQ6EMoU/6QJLC0ynI39l53+7V36EPbYNWOHhgC2mGVDTxkpbyet746jecP/oqGhRtHCvYy879P9ArJoHY2BR0GkF4hCkMf+wROCOTM4psqcV0V35YCsuTeR301LEDZ40VaZmDqdfVkDZwKBnZI2lt1pG/Zz0nGo4wacwljFLl8E3+YtLVGeyt2cGV4XNJixtMmDGMJbVfkBKTRrQhkoP7N5O3fhdTJ14JQGN9NX3TszvjsbSbuSM4Ykyx9nzXPHJHt7UjXHu6o55qJR1Z2eGO+P1r5O6nOOu7OGXSBL74KY83Fn4BEg7mbSI9LptrB99Oep/BbKr5iaFRI0kwxGDEyJctS6mrKqNCXUe/2EzSYwaxu3Y7B+r2ce09T1NYuBuAm+7/A5FxzrnFOYI9c5XWVoa4Ymm1lt/dvQdPjQ87MPdn7vhuS15nCLaQHsRkXEnl49X5/LT8Ux48/09UZerIOpXGc8vuZ2TueVSfKmPyiEv4re4CjOX1PHviGcIbw2mRLdz+4Ivcl/Es+w5uoG9NPIq4cCiDnIkTiUk09VW7to7ensvzdDzubHHdIYenEYFqMxmXmys3bNlq8Td/m2Ceet4cjJVaxsy6jK/+8xKh4WoStdHk6fczNno8B3UH0LVq+dOQl/nw2AIONuejDo8iQhlJyuBszu0zk7XHvsXQqGPipZdz7qwr+fg/L3Hh5TfRLyOb687LBOzvQrmaP77MX08bi7rG5UzY6hDlDillrjNxBrusHmTJxqMsWnWEmupy+oX3J/NoAv3HjSQpIpWBvUdwU9JtFOuK6KWO44FRT9F2soYCrWnv1FunPUFixgD271zHpvwVHNm3nSZDI9feeCPffDyfH5Z+xPwXngTg05+KzojbnnGhtblIW8+5C1tzobZkciRsZ/BFhdMju6z+0DK2trby+p9/hV6ro/jQHorFPjJ6D+XK6HksUv6b6BYVtBmp1VVza//rSK4M51hLMZnR2eQOPh+1UUXfjCz0bS3sObgBgNnX3019i4KLr7kDraaZiRdcCtDZQtqDvR5G1n5zt/Gku/AcXcXh6DjS0bRYG0O6kx6pkL7C3EOmpqqCTau+YmD2WKaMncuGnUt5fdXT/P6Gd2gqqyA7cwS61DBitiwjr3YX6a0xEBOKVt/C+uKVZEcMoUmlobm5nrjE3jRUV1NdXkpzrZa4xDi++mTBGfF7oyLy1Fxi1++OLrFyRi5HlbHrZ0/kd1AhzehaAzqa4YVFtZ2fjxUeQSEUzB52LWFaJbvD1nL/pX9h//4NaMPbOFFbyEbtbnqJKEpbT/B+26cojGHMHn8TG4tXseqIyYkpM2kocTGJGFraiIvvTUxiRKchx1/8OF3F1e6ls07slp631Apae9YT0x5BhTTDnabz44VHGJU5mWN5+9hQ/yOP3/0Wa9YuRlXaykVJc0jSJjI4LJvB465n99F19JZJrAhbQ1Job66f9ABDM8eTMGQQgzLHgN7Ali3fsHXtci6+eh5LNh7liskDeoQyguV5zA5sXXc2Pnt/667VtjdMRwgqpJNYKigdh6TqGrXk7fqRkFYFJ9PrqC4o5d0Pn8YQHcL1Q26lefdRInP7M7Apgi8Pf0B+1S4GRQ5GGSoJaxaE91Yz+bobMbQYECEK9I0aVq1+n3n3PAO47ghgK13edkNzZh7PG9MhlpTT1V6ULYIKaQVbmd31twVLNoEiHE1dI99//V8qy44xbMxUavYWADA+eiIRUTEknFSSkJBBY+UpPjzxDvWGOmanX8Pemu20tGhYfnIJ52deT/xBNYQqMCoE337/H5JS+jP2/BkeVUbzdHnSeGNpfGhLwfxlKsvTDghBhbSCI5mdt28fv75tDuHqSHQtzUTHJRKpjKHwxw0U1h9k9ribOHWqnEsjp2I0aJEtrejD2giPjuYe4x38oNtIlfYUs8fO47vtH5GTOJ64IbHsO7Ce9SUr0YcZ+eXTbxHZ3gLbwh3K5E5fUGthu7pCxHxqxN1zqtYWZ3u6QggqpBv4bNlaNM0NDI8bTUnYMWS5llOyirCoSJJUyWzPX809UXdgbNBiqGogNLM3ccDY0DH88fifiFLFMHBADuWlxSiFkkV5bzHIOJamxjr6jBjCb/70IuFqNWCy5Ha3vMrfDD2uTGnYmmJwNZ32trq2jDv2jjPtIegY4AZWLvuSi8ZcT01UE/WnTpEYlUJ8ZDL1dZU0yiaSI1JRRIYjEiMxGNqQzTraSir57PgHDIoaQpOunuxJE9lbvpUYVRzNjfXs2rwKZWgI8x58ulMZreGucZSr4TjzvLt9Qe3BlgNER6XW3X3W5HY1HUGFdBEpJdUVFaRHD+Tk8cPoDFqmTbmW/snZRKqj6aPqxyXxl/Hc0T+wu3Y7oTFRGKoaCBvZnwvSL6Wo6TB90rLQby7n9t+8xLQRczFKAyPGnsftjz5HdnqyTRncNa7xZDfXke5pR0HvbuGzKwuCzZWpO4XydkUBwS6ryxwrKaG+rpKClAKi1QlcnH0tS9a/TUavLOq1NWTEZRKX1IfftjxL0fF8WhT11BsbCatTUlR9kIzEIQw2ZLO/fjfTTl1D1JwEnhiTw2vvPkzhgV08cvP5vk6iR3FH6+6IwcfeLr35fY5YdANmx4Ceyv99tJzIqBiOnNqNrqaeJVXvER+WSIOihpTY/gxQpGOobESEKogfNwx1ZQzPH3iZKeUXcrSpAJqgVlFOYkI/tm1fwWcL3ubciVeRmT2aCedf5uvkuQ1nvWK6e9YT3jmuWnNdbVGDCukkSzYeZdn7H/DVJ68xfehcVmxfTEp4XyKiEojTRpLRnEKsCGds71yMtU3oFE38VPAj6QmmE6kK6g8AMCxsKHWxLWiMzaw+uYzQUBVbdy/nod/bb1V1FGsGCW9hyRGgu5bFW6v2bcnQ3TV3yRMcQzrB+rwyNq1czeeLXuXumU9TU34CvbGV45piZrVM5eZxDzO+/zRmj70ZNDqUidHExaUyrnU4Bwu3c32va4lTJ3HJ+FsoNpaQmpRJleYUyX37o9frGDpmEkPGTfDYGZjgm/GRP2PLd9baeNPd+RhUSCdZ/sV8rpr+C/qo+nH5mNuZN+GXxIcmENEvFX3ecWRLK037imjVtXC4cBfakxVkRQ3lwuFX8knjZxiMbXyz7QMevvxlkozxDBk5gdwpswkNVZG/cwO7NqzmSiecAOyxIHoLeyf5rc1Vdh3DeRNX5zWdJaiQTpC38yTHig/w+fdvIlMjEOpQdNX11Oir+fTwu4SN7I8ID8OIpEHfgFIoKdOdoqy8mAHNKTyf8xrqMNNucV/sWECf2AwO79/OVx++wWVjbiXUGEpVxQm3y+3tKQZL8XhzgbE92FtJdV2jaY3gtIcXWZ9Xxvq8MrZsW0J8WCI5w6fz8r/vBnUI4TP7ECtimTnkGtpO1BDSLx6VIozEqGQykrMIVYRSo6+luHAfAGO1g0mMSiFCG0b2pCkk9k6jsb6a7/Yv5rI7HiJ94FCHzuroirOeL55sQS1ZLjvoWlF4ozW3NWa0NLXi6QXcQaOOHXQoRlWDDoDq8lKyM8ZwUn+S6pYKfv3eNcQqY7lv7G9QHWmgwXiS8FM1tBhbiCIGTX09UcoITukqaNY1k16v5b8NHxAV0osS7SGe+vvl7TEJbn/0L7z1/GMAjBu7E3D8WABfWQidjccZ7yJfeSS50+JrCa+2kEKIWUKIQ0KIAiHEU1buuU4IkS+E2C+E+Mib8nVlfV4Zi1YdYdOGPaz/aQf1VRoAJJKfdi8ltSWZP1z8Fo/f/m/UIhzjqQZCY6KITu2Dqk8iirBQ2qJDiUxJJjoshuFRQ+itSkZTXc1TCb+mxdhCnaGOJ656jT7pWYDkm4/nAzB5xlziE207BXTFnw01thwHHBn7+nM6XcFrCml2PuRsYBgwr/0MSPN7soDfAlOklMOBR70lnyVKypsoP3KMF5+cx4tPzqO6vIzigkK2/LAMI0Z+LFyGUh1G8Y5tNIfr2VmzBZRKak8ex9jcgm5YDA01p5DNLSjUYahUaiIUairjm6kWtYgQJWmZQzhedIDhCWNRqSM5VniAUedM447Hn3fboTn+jK3xmCfHvd155XSN217ngEByDLDnfMh7gH9JKWsBpJQVXpTPItJoRKNpIEylRkrJ4n8/x6Qhs6BvOOuWf8q6vV9TF9nM5Vk3kqDpRVNpJfFp6YhIFYlVOmRiGhWtp9gr8ykuzyMqtBcFJwrRGVqICImiubQK5QUpNB4uQ6VSY2jRsXfrj+za8D03zxxmW8AAwhGHcm9hKT5XuqWBZNSxdD5k1wMpsoFsIcQGIcTm9iPQfUZkRCipgwfwwturePy590gZmE5ZSSErt3zEKO0ILp5yC/vb9jOyOZv1x1ZS2lRCr+RkFLERaAtPoqOV8ohaXi95lT3HNxKRmMB+3X5iYpNQCiVXT7ufe+94iUEVqVx0yT1IKbnstoeITUgmfeAQlmw86svkO427jTHemq5xdFMt8++u+Naa42/nQ4YAWcB0TMfVrWs/H7Ku64PeOB+yYzHwRz8IQsJDWf7xAnQ6LTdd+Vte/fJXpPXKYHS/CXxetJh6XQ3nxU9E0ScKFIJP2r7k8MEDPJv1V6LCoglJjGbi8FnUrD+Ftr6B8Kgovtj6NoaNeubN/TUjYnPoP2A4I4YP5oFHd3YenhOIYyVLi5w9teWGO7C2YNqa4cia15A7Kg1vKqQ950OeADZLKfVAsRDiECYF3dY1sI7zIb3BhpWf89G/nychqQ8P3fpPTp48Qmq/gcSGJPPtgY+Ji0jkwsRZnFLUcGL/52T3HsXladex8tQ3vFf+NhcNvppR0y7im10fUKOrpF9SFjJGSUXpUa658xE+X/QGm/KWU1KcR2XVVRQW1faY8aM9c3bedFYwj8+WsvmiMvSmQnaeDwmUYjof8sYu9yyh/RRlIUQipi5skRdltMiRvF20aJopLTlCq6YZfaqCshOFlFFIlCqGPqp+VBoqOaYtoU5TRYiIZXBpClUh1ST3TSe333moczIZfjCHbcZvKT6VT0xDAul9h3Fw7zZ++/LHHN67ndlX30d2zlhum9szxo5d5xU9PWVgL5ZaOHdUCu4Ix9/Oh1wBXCSEyAcMwJNSympn4rO2GZEt8/n6vDJKyptoqdWCkHz9yb85sHMTAEkxfVAoQ1j/7vvkRuTSmqZm76F1DB+Yy8DGPoQrwpGhGk6ENPCO9l0UoSFQ0sa6iFX8dP9DPHLb6wAIBPXaGuoLtnDV1IeJjk5kyqVXOJNMh/HW/J21LTDA8k4Avt7pwF1xuxpOjzzbwxXW55XxzbLv+OztV1CFR9BW0cTMyfPInjOT7/7+MmsLlzEhZgpFxiIqG8uYOfo6Dhfv5K4+90FlE8rEaLR1tZQpq2nT6ygLrULRO4pBsUMZMHosC398mX07fgRgeL/xVOjLeeh3b/PkgzPdnhZfYk3BLFWOgThO7o7g2R5u5n/vvEpbk5biw3spqz9GxoxzKfjpJ9YWLiM2oTdhhJKgTGRIRi7l5cWMSDuHT4sWcrj5EDqNhrDUBD6rXMwPTatZU76C9KwRpEWko01WIo2SlFjTUHr48MmINtj845LOuDvc81xxm3MFe9zJLN3fnauZNaukP9I1He6yntqLf+aKD5k6IpWYuCROVR9HSiN6YyvrvvgQo5SEKsOoqz5FmaGc45oSmiuqiW+Npl5bQ5mikg+1H7OveQ/5hzdikAYGpY4iRz2KLXtXIFva2PzVYvJ2/USTvhF1RDSbDq5kxDnTufCS23yd7E66a626m6Dv7rolQ4mlCXdvFnxr8XWVy9vL1FweQwohzpVSrnOHMN5mycajNGv0AKSnRFFS3oSUkmOF+VyXex+h4WrWF69g5+E1aE/VEKmIIm30CCa3TrYnD/oAAB5tSURBVGDR4fkkR/fFqDGwsWgNs6IuJua8LI7u2UlB62HC43rRSxVLnm4LlQcqKDqex4ghk7ni5ke58PKrCVFFExqmIj3F+sGrlnaY8/VYqyvOWCUdmYx3N44ovC8sru6ojq51Qxgex1KNWF+lobG0nsbSekrKmwAwaNu44JKbWbLvfWqryzknZxa9Q1IZP2AaD9/9GmF1gk+K3iOxfxoV2pPs5yDXxF/Ht43fUZ5/gKLWImZcfBvjJ1zGcV0JF8y4mcReKTQ0VVNZWEhLSyNDsjKIjY06bUeAji7q1BGpp/11TUNPxBPpsrfb6U+VGzjRQgohlgLFwE5ghzNheBNnXvb5c27m5PFCDh3Yg/HIZppVrSSOyubjT15iz/GNzL3hYSoPFVJ4Yic5CWMwGgX3xNxFWK9YylRV5BhHYEwOoziqF3uaDlLVWM7YrOmUaoq54IqbLcZpz7yjvxUeOLPF9pcWvDuDkr/IaAmHlUlKOVcIkQ6MxTSXmO52qdxIdxnfcSw4mNzkEqNVHCo4yjMPXEZoqIo4ZRwlbceYOekmThzIY8/xjfzu2rdYsmo+Ea1hzFCch6q1F5/Vf0qCSOACxWwSYlL5sXkN0+Ov4oeW1US0RJISnUaZ7ji/e+1Thg3q6i3oWaxNP3jKzG9tqsOeZ32FPRZhb+FQ8yGECAWQUpZIKb+UUj4jpbzEM6J5nismDyAmMaJTMSvrW/hh2efodFqefmAhrc0aJgyZydC4kXz4w9+ZN/oBlLUGDPpWtGo9n2j+R0bvYfRTpxMbk8QHRfNJDEtGF9VGhahm7jWPMP3yWyhvOM7wnKmEou7siiZGq0iMVnk8jb4cr3WN05qBxFPyWOuu2rthll976ggh3gGuEkI0Y3J52wvslVK+4SnhvEGHv+r73x7ku8Xv8PWnb2AwtLFhzf9oMDaw98gGQowh1GurWbz730zPvISjTUXcY7idG8fdiSIsBFVbCH0TB6CRGsZOm01kRAyG4lp0KSGoqyLIGXsBanUvwqLCOuMNdNc4e1dueNM9zloc/tIS24MjuXQu0FtKmQZcBXwJRHpEKhfozhO/K+ZzfUaDgb07fuDSmfdyde49fLf5QybPvpqxQ85nQPwgBILbom8huVyFzqDlYGQJv9/xCGWlBdSGNLKzfCP6XoKqomLqNhykbWIyqjA1jc21FBzZSZ/+gzy2raMvsLZe0Nb9XT87gi+U2lKcnpTDkTHkZiAOqJBSlmLyR/3WI1K5EWv+hR3K2PH/WFE+RYf3UHR4D88+tIjImFiWf/sBsamp7C/ajESyXbGPWllFsiqFXfXbeXr0C9QbGxEn9ShVSupqK2g6Uckbu/7GuNZZxMQm8dOqT2hra2Xhm08zfe7pGx93141y5/jFW2Mhd8This+rtd9daTm93bo6opDzgR+FEAuALZi6q/WeEct5LPmvdmDtxVTUaenTL5vM7NEUHd7NyW37WL5vEeEinAgRQVpYGnNUF/NR3UckRaYSnpKAprKGYwf2slysptJQQW9jH3rHphGX2Y+Xr1yJNqKNf//tQe785Usk90sne/QQ4Oe5RVu1rDsLgr902cy9eszx1QoLS657vlZSRxTyQ+D99mceAEYJIcKllAM9IpmL2Dtw37JtL//56+M0NdaQkNCX2659lujoFJIPJNGgryWkspU92jyOyINcNeVe8vI3cHHytXx2aj6lqgqUynAevuglDioOsfzLtzlUtB3xuUClimDqzOsYOf481AmRtLbBtdMHnCGPP5vgHaU739TuegPujt8cb3vauIojCnlCSvkH8wtCCM+bCe3Enr1OzF/MgcIyPpv/CutWfEp4WCQ6vYb+mcOQtS3E6FVMmXMD/cri2FK7nvLCCnLTp6MamUJ6eTYLdrzCoEFjMKJGXRnOh6tfJSI6pjPs1PRMXv3vMjS6n9dk3zQjy6JcgVRYusPSVom20mbvWMxZD6BAdKRwRCF3CyEekVK+1nFBSqnzgEx2IbBc69q61sGKzxaw8fsvAOiXOYSCgzvYs3UNoWkGtslV1JeXc83shynem0d8azSlzSVcGjKbrOwUjNLIsMsvIXWnkuWhX1G4bwnXj72GYycPERYewY2/+A2KEBWRSunxI8g9jb0teFfngO5+d+Q3VwjEys6RKqQ38AshxEkhxDIhxPNCCJ+7zVnqkljrpphbWy+YfTMGYxv9E7KY2t+09GlU2iS2H/+JkroCQkQIB47vJCfrPKJ7pzJx2EW8seARdJFG0jKHEntMgSFaQaIqmXOyL+STn94kNS2TR//8JqMnXUCgLmtzB77oJvpqdYa7sauFFEL0Bb4B1gOHgVPAKGAC8JnHpLOBvZPMFhchF+0A4Fj1EbIvmcnDYwehXXuEBGUCW4+tYVD/HI4V7ydfbuXB3Geobq3BKCSLV/8DhVFJm7INERJCaW0RfdOyGDpqMrnnzmbMxGmeSayP8OR4z9OGK1tx+OP43aZCCiEuAhYCa4BW4A6gF3CHlHKhZ8Wzjivtz5eLl/Dea0/zxLMfkNKWhF4Ioo9J4nNHk9E8hlnK22mpa2TBN3/CGCUIGZuCskhDWIgKaZQM6jeS/U27KS05iCo8gjHxE4gY0JuvF7/JrGuuR6FUBnxX1ZOYW1s9qRDOTpP4EntayOeAc6WUBR0XhBCTgLeFEPdKKTd6TDo3Yr6lYnxSKtqWJthQRnlmMwc3b2fZD+9w18RfM/TKOXCqmQijYOCg0azZt4SPN7yB4Wg9UqtnwuhLyc6dRO8jGXxQ8jK6Fg2oFOQVbiE6LgGhsD6xbG2c5Y8Fw1H8sbUJROzpaIeZKyOAlHITJm+dv3pEKjvoMOo4OmZYt+JLPlnwVwanj6NZ2cKrb/+C4w1FPHnhqwyeOBV5qBqSI6jtZ2THwbVEhEZxSeIVRKYk0WRo5NNtb/HcW7eS2ds0r5jcJ50S41GGjpvAb15+GyHEaa2jNU8Wez1cLOGp5UrOxGXvjgDm78nZtAfquNAR7GkhW4QQSVLKSvOLUsrDQogYaw95GoljLUt9lQYpJf9+/gkAonrF89bKZ7kq507W7/+W3g89ieFAFahDkDvLOdmynzb0/OLefxDdGEdmYxYbjV8D0CduABt3fsNN4x9izanl3PP438kaFA94z0e1u8l1Z3C3dwx4Z5xoHpe99/oz9ijkK8ASIcS1UsrOfVTbt2kMqCpLCEFoWDiRvaJ5+YNV3DtnFDtLN5DdbzQffPVXhsXm0NhQw7ri5Vxyzq3MOedWtq38ErUigoyUn7dm1Co1bCpdTU78OYSGqkBKrzqL++P6Q2t4S1H8OQ8cwaZCSSk/x3RIziYhxBdCiD8IIZ4HNgF/87SA3WFvdzVcqehcZpWRPZLIyFgen3c+AEerDhGblsapkkI+3vIvSiJKqW+pIbZPX2g1MubSK9hY9D2LNv4TAGVIKI11tfTpn8nRliJCQ8NQeLkr1TXdnuzK2TMscCV+R4YbjsYTiF1cu7eBFEJEYFqQPAKoB5ZJKXd4ULZuyc3Nldu3b7f428frijo/m685LClvorEsnz8/ej+zz7uNikMFrDr8BRFRMdx+xe/ZtWkF08dehTIpkmPGY4w0DAMpKSzNY+PObzCoYP+xrUTHJZKRPYzjhQXc+sCfmTDj/LPOqtrVG8fVwm+t1bflGG6tBbbXW6g7nO19eGUbSCmlRkr5rpTyV1LKP/lSGbujuwzUt7aycP6/GJI2lkGadFq1Gvr0z0LTVM/Hy14lfHBvmk6colnbwPC6gSAlxiotaW190bfpuGHkvSTF96WlqYnUXoMICQ1lyOgJZ50ydoctg401A5ej2NuyWnLpsxe/XqDsDtpPs3oN087l70gpX7Ry3zWYHA7GSyktN4Pt2HLZ6jj12Gg08vqz99Hc1MCxgv3oGzRsP/5j530Z0VlsWP4p+6OTeezm91HsrkNkx0N8Ex9/9QIzRl7NyoplNGsbCQlTsfKHhUyaZtptfOHSfO67coRDedFTxjzOYq8rnbP55C87JTiK1xTS7MDWmZgO1dkmhFgqpczvcl8v4GFMS7zcRm3VKU4UHyJWxgIwbdRcth8zKWTf3gPZdWIjRqOB2LB4aDSAlKA2LSiOT+jDBu1G2ow6Lrn+fmITeqNUhHDunFkkx0VYjbMnY89eOv5kfPJ1/Pbibwe2AvwFeBl4orvALE17dHw370Iu2XiU+ioNdZWNhIap6B2ZRmuogTj1z8eFjx42nYadS+ifnM3EWVfBtlMQp6Z1ZQFNfSWX3/k4330+n907V7F3xxoAXn7np05lnDoi1epC6EDAE04K7uiWuovujF+2FjX78wJlV7F0YOsE8xuEEGOANCnlMiFEtwppD0s2HiUxWkWzRk9UQiyapgauvPER1i9dxB++uIMpuZdzvLqAIfpBGBPOZ8C4CcS1JhKaGU+9poaF+17i2NJ8zs25jHV7viY0REVmVg6Dho0jIrLXGfHZa+DwdU3tq8LmK5xJp6/yxm8ObBVCKIB/ALfbE1hVZSVTzQ5svcPKga0dJ1mFEM7IsdNYse9TLp10Ha3DohiZeg6jxQUUa4+RM/s6ojc2s7/XMVRaNW+//1uICWXCtMtY9+PXDB05mdR+A7ny9kcZMrg30DM2quowsniqSxcoXUV/wZ8ObO2FaUplrRACIAVYKoSYa8mwk5SUhPm0h62Xbmg1cM2dT/LXJ66n79xsZg+/Afr2AinJqBtGU0U1q44sIf14BkuqX6e88Tg0Ql31KZ5+8XPSRw3F2GIgJjGCqgadQ1s4+luB9KbBw9buAZ6KI1DxmwNb2/fnSez4LoRYCzxhy8raQcdLNp+DbNboTec8AhXlJQwYNYxX317MPddcyGLgjy9+Q4SqF8dL8vly7ZucKD7C1ePuYVepyV/+F4+/gbalnlC1CoVCgSJCQbNGb3P3OHduTuXuMN0lgyNK5ex40tGFzj0Bfzuw1f7wOH3MtmDFIYDTlEVT38Cit/7Mobxt1NWc4qJL7+SCa24AYFj2RIRBsviNp9lV9PNZQYrwn7Nk07ov2bP1B4aPmcpjL7zjVLpdwR8Kn70+rp4waPk6/b4Ya/eIA1vNFbKjRVTHCv704I0kqHsza/QNbK/bzPffvIcqXI2uxXRP9uDxHD60DZU6Ap1W0xl2ZGQssfHJNDXUUl9fyW9f/ZhR54wnMVp1WnfV2hjS3406nqAndkOdVciz8sBW8z11LGVY3s7ttFTUc+OlTxLbHMlFs+5iWPZElCKEqSNN+6NGRyfwys2fMr3PbACSovoQG59MRvYo7nrqJerrK1EolCT37W9xzGjtUFVzbxVLfz0dW+m3J1/8YUsOX7wvvz65qjs62vUz1t3FqQFY8/VHpMcOZNXBz0lUJNN7dSVHju5C36ojoU8/zlHN4YLx10IjTB53CTnDz+PdLS9zx4MvEqaGoSNH8b81e5k4LIWEhITT4joblMpR3D1/6awjuTt9WX1BwCpkV+66eDAA7648jF6nR6WKZFvZd6TuLkfXoqFN38qtEx5jcGoOir69ID2GE4fyeHPj3+grU9l6/EeGjZpCRuZIhnZ2RVUcKGtlqkkfA+7lugNLHjfmeGq9o6Ph+pMjgisEtEJaKig71n3P2y8/gVKpxGgwcNUdjzH5gotY8t47rPjhS4YMyoU2I8ZWPd+ueZcWbRPV8fUMGjqG3v0HoE6Icvv8oj/PxbnaFXSmJfIH67G/ErBjyK509Pe//uhf3DjmfmYPvZ7hOVMYNnIKzRo96Zkj0bVoEJGhkB7D96vfxyANPPXyR7TqWzh6YB8bV32BrkXDolVH3C6bv2LPuM6VcZyl586m8bSjBHQLaYmhORNZvf5rquvL0LY28+hNE7jilgfYu2UTg4aMA4WCnfvXsHXrN/zulU8wGAyUHj1Cq17L0y//jxAZRnpK1GlHjJ/NuKo0rsw5ulsWR/BVKx6wCtl15/IOvv3sbe58+AUqykrI37MRTVMDjdXNaDWNFB/Zw7Hi/bS0NPPIc/MZm5sNwP/W7gFArTZ54ZhPbXjKuhdsHRzHm11/T7oTdkfAKqQ1lEol55x7KVUVpSz/8h2Sk/qzbtUnNDbUcsn1DzBq3HlkjRrKvBmnr1+0NoXhKQLJC8VSwfTWRsf2yOIpevwCZW8RHqemb+xArpz3K3ZtW8kVlz1GfX0lEy+4DHVEFOoI62sYbU36O4I7Xqg/G0A8PbVgPoZ1djok0OgRnjpw+kvraO1KypvO8Dt1ZbsNawXQ1r4vruCM14+1AuwJH1tHZbMUlr9VNK7iiqdOj2khrb3UZo2+83N6SpRH4vBkgfLHwhrIi7H9nYBWSPPjAeDn1q+juzl1xOljw5LyJqoajp52r714snvW0Uq4Oq50ZUMnT+Cp7R39sZJyFwGrkJZWO3saT4+VXF1B0dVRwlPyutt5PpAMXJ4mYBWyK9ZavI7WsqOltGc9ozm+bGnOlkLpjXQGylg1oBXyhnMzz8hk8y5qhzIuWnWENrOxZKKTY0lHWgZfF4CeMq51l+9sd8/5kyN6QCukszgypeGJAnA242j39GzLx4BXyI4XPP9r026ShlYDrY06ho9JZX1eGTNy+gIQ4kA3NYjncMdOA/a2mq5sM+IrAlYhu+7L2qpr4aP/+ws54y9gVO70zustBiM3zcjyvoBB7MIZRbD3GV8omas2h4BVyK5Ul5exafUSDPo2RuVOP2ucwrtzVvCXWv9swR0GwIBVyNomXec85BWTB5Cc2p83P9nlW6F8gDedFfzJ+OEu3FlxucNhImAVsoPEaBXr88qIS448bd+b9XllZ0x5dGCp9XQ2I515mYE+72ZrDGdJcX2VZlsK56oM7q6kAl4hO+hQxhk5fZ32bvEWgaR85vjz2M3XuKsc+d6/yknio1TccG6m1+KzZ6e0noQntnXsaXnnico8oFvIVXtKz7hm6UV3Z+DpbiGqt1tPT64acYaeYBiyZ1znT3PNfnVgqxDiV8DdQBtQCdwppSyxN/xAt6y6Yz7N3jC7o6uTek9XSkeWqHUXjjvyyd8ObN0F5EopNUKI+zGdE3m9o3HZU4jO2M/VRqH35ULhQFYIf9srxxbdGWm8YWX2qwNbpZRrzO7fDNxsb+DmraM9GeZKpvZE87859i4sDiTc8a688b796sDWLtwFfGftx6Mnyrjnulmd35947EGL50N6EntM+Y4qb6AquzOW7UBLozfwmwNbT7tRiJuBXGCatcCiomPJPe9iLrhsHpG9YrjuPM9YXF3drsLRQufLQuqtVs9bafT0eM8TFYs/HdgKgBBiBvA0ME1KqbMWWHVFGasXLeTw3u08/Of/uEVAbzoH+COOtuK2rjkarjextpOer/GbA1sBhBBjgP8As6SUFd0FJhDUaqtIVyjdJqA/FhxLBJrl01cGse4UzpIcjli5PZUOr+46J4SYA/yTnw9sfd78wFYhxCpgJNDh63ZMSjnXUlgZ2cPlpTfcw8jxU4iONZ2GY2nBck/B02dO+lN3tbsKxxeVkaUDh7qTwZVd5wJ2G8jc3Fy5ffuZp50HokI6owzusBI7i7NxB1rL3oE3FTKgPXUsEagv3R5PIVcUwR24w0fYW2lydxkwL1c9ZR7Sa7jrpfsady4LcgZnFTnQnAFs4U1Ze6RCdsXZQbm3BvOBUjht9T46fvdGejrisLcldjZ8bxOwCtl1Cw9ruPOFBYrieIKzLe2OOPq70yAWsArpDGdboXIX7liE7YkxnSvY09J7WgZL9FiF9IdJ3kDAnflkXoj9pfLzZIsW8Muv3Im1A1uD2KYn5JurLZilqQx/IGAVMoj38MTuAd7C3nWQ3f3mze7rWaWQwVUGjhPIeRWIsp8VCukPLyaQW5lAxZXpLld2DnCFs0Ih/QF3j1ecddh2NG5nCrU/OGbYM93lypI5T6UpoBWy6+SwvQWlp7QuzqTDW5t5+cIB3NcyuIOAVkhzfFmreRNPegp113p03fiqu7B6Qj47gjsr+oBVSPM1KmdbAfAU3bmjeTKP3XFKlacdD+xZK3lWn+0RCHiiK+gN45A/7Krn7nF2d3HaE5e35imDCulFbDlEOxKGNfxlgttT2LvsytvbsZz185Ddeer4g0HBVZxNQ6B23z1VkQSabSFgFbI7rG1g5KmM9+S0w9mCu9+Nv3gDOUqPUsjuXoIjO4x58mXaWmXg6fjN47EHd3Xh7Em3u2XwFJ6sVHuUQjqCN6yG/laQOvDFOLRrnjgShy+PcfA2AauQXbfm8qcX5arV0NdjY2e74F2fc6THYi3Ms42AVUhw79jN04XdH5f6uIMWg+WTpVzJT3+qXL1NQCukI/jrS/ZXuRylp6TDHqxVQu7Aq1W1EGKWEOKQEKJACPGUhd9VQohP2n/fIoQYYDUsPNvSWNsuv+PPnfF4Kh2e2gCqK55QRk/mi7vocBPs+ucKXkux2fmQs4FhwDwhxLAut90F1EopBwH/AF6yFl5lZaXbCoKll2/Ng8U80xe8Pd/leNxBd3KYVyKeNDY5mhfd4cjudV3T9/6Cd9wmR3fx2SDR2fC9WQV1ng8ppWwFOs6HNOdyYGH75/8BFwohLJ2aRVVVld0R28pEZ2u2d99x7OVbq1FdVRBLcni699D1z9G88BTz58/3aNptvbP2uJOcDd/fzofsvEdK2SaEqAcSAPu1zwJn0/img45xjifmPd2Zn852qztk6CqLvduD+iv+dj6k3WdIajQanTpEaTC7VImLiusEieoQpbfjtIQ/yOEPMoB/yDHY2Qf97XzIjntOCCFCgBigxlJgUspwTwgZJIgv8eYYsvN8SCFEGKbzIZd2uWcpcFv752uAH2SgHs8VJIgTeK2FbB8T/hJYwc/nQ+43Px8SWAB8IIQowNQy3uAt+YIE8QcC9nzIIEF6Iv4984p7nQk8KMOvhBD5Qoi9QojVQoh0d8tgjxxm910jhJBCCKcODXVVBiHEde35sV8I8ZG7ZbBHDiFEfyHEGiHErvb3MsfN8b8rhKgQQuRZ+V0IIV5vl2+vEGKsXQFLKf32D1PXthDIBMKAPcCwLvc8APxf++cbgE98IMP5QET75/vdLYO9crTf1wv4CdgM5PogL7KAXUBc+/dkH5WL+cD97Z+HAUfdLMN5wFggz8rvc4DvMM0cTAS22BOuv7eQbnUm8JQMUso1UkpN+9fNmCzI7saevAD4C/Ay0OIjGe4B/iWlrAWQUlb4SA4JRLd/juFMi75LSCl/wsoMQDuXA+9LE5uBWCFEqq1w/V0hLTkT9LV2j5SyDehwJvCmDObchalmdDc25RBCjAHSpJTLPBC/XTIA2UC2EGKDEGKzEGKWj+T4I3CzEOIE8C3wkAfk6A5Hyw3g/6s93OpM4EEZTDcKcTOQC0xzY/x2ySGEUGDy/73dA3HbJUM7IZi6rdMx9RTWCSFGSCnrvCzHPOC/Usq/CSEmYbLej5BSesuNx6ly6e8tpCPOBNhyJvCgDAghZgBPA3OllDo3xm+vHL2AEcBaIcRRTOOWpW427Nj7Pr6SUuqllMXAIUwK6k7skeMu4FMAKeUmIBwXnL6dwK5ycwbuHnC7eeAcAhQBGfw8eB/e5Z4HOd2o86kPZBiDyciQ5cu86HL/Wtxv1LEnL2YBC9s/J2LqtiX4QI7vgNvbPw9tVwbhZjkGYN2ocwmnG3W22hWmpwqQGxM9BzjcXuCfbr/2Z0wtEZhqvs+AAmArkOkDGVYBp4Dd7X9LfZEXXe51u0LamRcC+DuQD+wDbvBRuRgGbGhX1t3ARW6OfzFQBugxtYZ3Ab8AfmGWD/9ql2+fve8i6BgQJIgf4e9jyCBBziqCChkkiB8RVMggQfyIoEIGCeJHBBUySBA/IqiQQYL4EUGFDBLEjwgq5FmAEOLq9rWie4QQ24UQF/tapiCW8Xfn8iAuIoS4EdNKh8ullOVCiCxMDt/jpZTHbTwexMsEPXV6MEKISOAAMMVc+YQQHwPLMflbngJGY3KEvgm4F5Pv5Top5V1eF/osJ9hl7dncAOy00BLqMK2KGQkUSSmnYlrkvQD4DaZVI1cJIVTeFDZIUCF7OiMwOVd3JQc4AsQC/2y/pgUWSCnLpGkVvgZoFUJECiEWCiHeFkLc5BWpz2KCCtmzqce0PKmT9sW60Zi6qjvlzwt2c4At7ff0A05K03jmKuB/Usp7gLneEvxsJaiQPZtlwHVCiCQAIUQ28A5wB6buqnnrOQrY2/45x+xzP37eisL86IYgHiBoZe3BSCm3CyH+AqwWQkhMLeYvpJTrhBBzMa0fRQgRDqhl+8ZUnK6cJzAp5W6CFbjHCVpZg3RLu6X2TUy72K2XUi7ysUg9mqBCBgniRwS7IEGC+BFBhQwSxI8IKmSQIH5EUCGDBPEjggoZJIgfEVTIIEH8iKBCBgniRwQVMkgQPyKokEGC+BH/D8+PvY74qVU8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "conflevs = [0.6827, 0.9545]\n",
    "\n",
    "plot_2d_dist(x_1, x_2, xlim=[0,1], ylim =[0,1], nxbins=200, nybins=200, \n",
    "             clevs=conflevs, smooth=4, xlabel='$Om_0$', ylabel='$Om_L$', figsize=(3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing these results to our results from HW03, we can see that these results map very well to our previous ones! The values are localized in a similar area, with the shape and slope of the oval matching. This suggests that our parameters previously obtained were pretty good, and also that the MCMC sampling is working correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "AutocorrError",
     "evalue": "The chain is shorter than 50 times the integrated autocorrelation time for 5 parameter(s). Use this estimate with caution and run a longer chain!\nN/50 = 20;\ntau: [45.61647444 46.80794996 33.57586347 41.56337973 39.44313983]",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAutocorrError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-16-5c662610fd07>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\ensemble.py\u001b[0m in \u001b[0;36mget_autocorr_time\u001b[1;34m(self, **kwargs)\u001b[0m\n\u001b[0;32m    524\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    525\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 526\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbackend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    527\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    528\u001b[0m     \u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBackend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\backends\\backend.py\u001b[0m in \u001b[0;36mget_autocorr_time\u001b[1;34m(self, discard, thin, **kwargs)\u001b[0m\n\u001b[0;32m    148\u001b[0m         \"\"\"\n\u001b[0;32m    149\u001b[0m         \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_chain\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdiscard\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdiscard\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mthin\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 150\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mthin\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mautocorr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mintegrated_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    151\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    152\u001b[0m     \u001b[1;33m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\autocorr.py\u001b[0m in \u001b[0;36mintegrated_time\u001b[1;34m(x, c, tol, quiet)\u001b[0m\n\u001b[0;32m    108\u001b[0m         \u001b[0mmsg\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;34m\"N/{0} = {1:.0f};\\ntau: {2}\"\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_t\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtau_est\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    109\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mquiet\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 110\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAutocorrError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtau_est\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    111\u001b[0m         \u001b[0mlogging\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwarning\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    112\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAutocorrError\u001b[0m: The chain is shorter than 50 times the integrated autocorrelation time for 5 parameter(s). Use this estimate with caution and run a longer chain!\nN/50 = 20;\ntau: [45.61647444 46.80794996 33.57586347 41.56337973 39.44313983]"
     ]
    }
   ],
   "source": [
    "sampler.get_autocorr_time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emcee time = 189 sec, acceptance fraction = 0.547\n"
     ]
    },
    {
     "ename": "AutocorrError",
     "evalue": "The chain is shorter than 50 times the integrated autocorrelation time for 4 parameter(s). Use this estimate with caution and run a longer chain!\nN/50 = 40;\ntau: [44.58065362 46.51340912 22.86415133 40.62321482 68.37274824]",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAutocorrError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-17-8f1d2f244c68>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     22\u001b[0m \u001b[0mx_1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatchain\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     23\u001b[0m \u001b[0mx_2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatchain\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 24\u001b[1;33m \u001b[0msampler\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\ensemble.py\u001b[0m in \u001b[0;36mget_autocorr_time\u001b[1;34m(self, **kwargs)\u001b[0m\n\u001b[0;32m    524\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    525\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 526\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbackend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    527\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    528\u001b[0m     \u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mBackend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_autocorr_time\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\backends\\backend.py\u001b[0m in \u001b[0;36mget_autocorr_time\u001b[1;34m(self, discard, thin, **kwargs)\u001b[0m\n\u001b[0;32m    148\u001b[0m         \"\"\"\n\u001b[0;32m    149\u001b[0m         \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_chain\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdiscard\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdiscard\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mthin\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 150\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mthin\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mautocorr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mintegrated_time\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    151\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    152\u001b[0m     \u001b[1;33m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Anaconda3\\lib\\site-packages\\emcee\\autocorr.py\u001b[0m in \u001b[0;36mintegrated_time\u001b[1;34m(x, c, tol, quiet)\u001b[0m\n\u001b[0;32m    108\u001b[0m         \u001b[0mmsg\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;34m\"N/{0} = {1:.0f};\\ntau: {2}\"\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn_t\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtau_est\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    109\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0mquiet\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 110\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mAutocorrError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtau_est\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    111\u001b[0m         \u001b[0mlogging\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwarning\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    112\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAutocorrError\u001b[0m: The chain is shorter than 50 times the integrated autocorrelation time for 4 parameter(s). Use this estimate with caution and run a longer chain!\nN/50 = 40;\ntau: [44.58065362 46.51340912 22.86415133 40.62321482 68.37274824]"
     ]
    }
   ],
   "source": [
    "ndim = 5\n",
    "nwalkers = 10\n",
    "nsteps = 2000\n",
    "args0 = []\n",
    "\n",
    "#differential evolution yielded [ 0.25747937,  0.63495553, 24.07862841,  0.11942547,  2.5793592 ]\n",
    "x0 = np.zeros((nwalkers, ndim))\n",
    "for i in range(nwalkers):\n",
    "    x0[i][0] = np.random.normal(0,scale=.01) + 0.25747937\n",
    "    x0[i][1] = np.random.normal(0,scale=.01) + 0.63495553\n",
    "    x0[i][2] = np.random.normal(0,scale=1) + 24.07862841\n",
    "    x0[i][3] = np.random.normal(0,scale=0.01) + 0.11942547\n",
    "    x0[i][4] = np.random.normal(0,scale=0.1) + 2.5793592\n",
    "\n",
    "tstart = default_timer()\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, posterior_pdf, args=[])\n",
    "echain, prob, state = sampler.run_mcmc(x0, nsteps)\n",
    "tend = default_timer()\n",
    "# acceptance fraction averaged over all walkers\n",
    "eacc = np.mean(sampler.acceptance_fraction)\n",
    "print(\"emcee time = %.3g sec, acceptance fraction = %.3g\"%(tend-tstart, eacc))\n",
    "x_1 = sampler.flatchain[:,0];\n",
    "x_2 = sampler.flatchain[:,1]\n",
    "sampler.get_autocorr_time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emcee time = 481 sec, acceptance fraction = 0.55\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([67.19498625, 58.14679249, 36.92435422, 62.30623738, 61.44491265])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ndim = 5\n",
    "nwalkers = 10\n",
    "nsteps = 5000\n",
    "args0 = []\n",
    "\n",
    "#differential evolution yielded [ 0.25747937,  0.63495553, 24.07862841,  0.11942547,  2.5793592 ]\n",
    "x0 = np.zeros((nwalkers, ndim))\n",
    "for i in range(nwalkers):\n",
    "    x0[i][0] = np.random.normal(0,scale=.01) + 0.25747937\n",
    "    x0[i][1] = np.random.normal(0,scale=.01) + 0.63495553\n",
    "    x0[i][2] = np.random.normal(0,scale=1) + 24.07862841\n",
    "    x0[i][3] = np.random.normal(0,scale=0.01) + 0.11942547\n",
    "    x0[i][4] = np.random.normal(0,scale=0.1) + 2.5793592\n",
    "\n",
    "tstart = default_timer()\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, posterior_pdf, args=[])\n",
    "echain, prob, state = sampler.run_mcmc(x0, nsteps)\n",
    "tend = default_timer()\n",
    "# acceptance fraction averaged over all walkers\n",
    "eacc = np.mean(sampler.acceptance_fraction)\n",
    "print(\"emcee time = %.3g sec, acceptance fraction = %.3g\"%(tend-tstart, eacc))\n",
    "x_1 = sampler.flatchain[:,0];\n",
    "x_2 = sampler.flatchain[:,1]\n",
    "sampler.get_autocorr_time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the errors thrown above, we need at least somewhere in the vincinity of 5000 steps to start getting reasonable chain lengths. However, as you can also see by the runtimes, that starts to get expensive fast - so 5000 is the most I will run with (more than that and my machine really starts to suffer)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsphere_volume(r = 1., d=1, nrnd=1):\n",
    "    '''\n",
    "    Parameters:\n",
    "    -----------\n",
    "    r: real\n",
    "        radius of the sphere\n",
    "    d: integer\n",
    "        number of dimensions\n",
    "    nrnd: integer\n",
    "        number of Monte Carlo samples to use in the estimate\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    real number - estimate of the d-dimensional sphere volume\n",
    "    '''\n",
    "    k = 0.\n",
    "    \n",
    "    for i in range(0, nrnd):\n",
    "        points = np.random.uniform(-r, r, d)\n",
    "        \n",
    "        vector = sum(points**2)\n",
    "        if vector < r:\n",
    "            k += 1\n",
    "            \n",
    "    ratio = float(k)/float(nrnd)\n",
    "    volume = (2*r)**d * ratio\n",
    "    \n",
    "    return volume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsphere_volume_exact(r=1., dim=1):\n",
    "       return 2. * r**dim * np.pi**(0.5*dim) / (gamma(0.5*dim) * dim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sadie\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "D = 10\n",
    "    \n",
    "estimates = []\n",
    "exacts = []\n",
    "errors = []\n",
    "x = np.linspace(0, D, D)\n",
    "\n",
    "for d_i in range(0, D):\n",
    "    estimate = dsphere_volume(d=d_i, nrnd=1000)\n",
    "    exact = dsphere_volume_exact(dim=d_i)\n",
    "    error = np.abs(1-(estimate/exact))\n",
    "    estimates.append(estimate)\n",
    "    exacts.append(exact)\n",
    "    errors.append(error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x21178147f48>]"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAD2CAYAAAD7/h1SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAe+ElEQVR4nO3dd3jUVaLG8e8hlSTUJJTQQpcAoYVuF1QW1gK6NhBdBRELYru6tnXdRVx3V2W9llgQBSkiKmJXVgFXgYSQ0EmkhNCSEAihpJ/7R6LXdUEmZGZ+k5n38zzzkDJM3pHk9eTMOb9jrLWIiIjvqed0ABEROTEVtIiIj1JBi4j4KBW0iIiPUkGLiPioYHc+WExMjI2Pj3fnQ4qI+LXU1NR8a23siT7n1oKOj48nJSXFnQ8pIuLXjDE7T/Y5TXGIiPgolwraGNPYGLPQGLPZGLPJGDPY08FERAKdq1MczwGfWmuvMMaEAhEezCQiIrhQ0MaYhsDZwA0A1tpSoNSzsURExJUpjg5AHjDTGJNmjHnVGBN5ojvm5eWRlJT00y05OdmtYUVEAokrUxzBQF/gDmvtSmPMc8ADwCO/vGNsbKxWcYiIuIkrI+gcIMdau7L6/YVUFbaIiHjQKQvaWrsP2GWM6Vr9oQuAjR5NJSL/pbS8knmrstl/uNjpKOIlrq7iuAOYU72CYxtwo+ciicgv7cg/yp3z0sjIKaR/fBPmTxxMvXrG6VjiYS6tg7bWrrXWJllrE621l1lrD3o6mIhUeS8th5EzlrPzwDGuHdiW1TsO8tb3J918Jn7ErVu9RcR9jpSU8+j761mUtpsB8U159uretGwUzq6CYzz16WbOP6MZbZpqS4I/01ZvER+UkXOIUTOW8/7a3Uwd1oW5EwcR17g+xhieHN0TA/zhvXXoyDr/poIW8SGVlZZXlm1jzIv/prS8kvm3DGbKsM4E/Wy+uXWTCP5nxBksz8znnZQcB9OKp2mKQ8RH5BWVcM876SzbmsfF3VswfUxPGkeEnvC+Ywe2Y0n6Xp74aCPndI2lecNwL6cVb9AIWsQHfLM1jxHPLWPltgP85fIevDi270nLGaBePcP0MT0pLa/k4ffXa6rDT6mgRRxUWl7JtI83Mf71VURHhvHhHWdy3cB2GHPqJXQdYqO4e3gXvti4nyUZe72QVrxNUxwiDvn52uZxg9rx0MhuhIcE1egxbjqzPR+t28sfF29gSMdooqPCPJRWnKARtIgDfr62+aWx/Xjish41LmeA4KB6/PWKRA4Xl/H4h9rg629U0CJedKSknLvnr2Xq/HS6xzXikylncXGPFrV6zDNaNOS28zqxOH0PX27c76ak4gs0xSHiJRk5h7hzbhrZBceYOqwLt5/f6T+Wz9XG5HM78en6fTz0/jr6t29Ko/ohbnlccZZG0CIeVllpSV72A6NfOPna5toKDa7HU2MSySsq4cmPN7ntccVZGkGLeFBuUTH3LEhneWY+I3q0YProRBpFeGZ026tNYyac1YGXl21jVGIcZ3aO8cjXEe/RCFrEQ77ZmsdvnlvOqu0FTLu8Jy9c19dj5fyjqcO70D4mkgcWZXC0pNyjX0s8TwUt4mYnWtt87cC2Lq1trq3wkCCeGpNIzsHjPP3ZFo9/PfEsTXGIuNH2/KPcOTeNdbtPf21zbQ1o35TrB7dj1nc7GJXYkqT4pl79+uI+GkGLuMmiNTmMmrGc7IJjvDzu9Nc2u8P9F59BXKP63P9uBsVlFY5kkNpTQYvU0pGScqbOX8vdC9Lp3qpqbfNF3Wu3trm2osKCeXJ0T7blHeW5rzIdzSKnT1McIrWQkXOIO+amsavgGHcP78Jt57lvbXNtnd0lliv7tSZ52TZG9mxJj1aNnI4kNaQRtMhp+Pna5rLqtc13XuDetc3u8PDIBJpGhnLfwgzKKiqdjiM1pIIWqaHcomLGz1zFtI83MzyhOZ9MOZv+PvpCXKOIEP58WQ827T3MS1//4HQcqSFNcYjUwNdbcrn3nXSKisuZdnlPrhnQxivL52rjou4tGJXYkhlLM7moRwu6NG/gdCRxkUbQIi4oLa/kLx9t5IaZq72+ttkdHr+kO1Fhwdy/MIOKSl3cv67QCFrkJA4Xl/H9Dwf4NiufpVty2VVw3LG1zbUVHRXGHy/pzpR5a5n57XZuPquD05HEBSpokWplFZWkZR9iRVY+KzLzSM8ppKLSEhEaxMD2TXlsVHeGJTR3OuZpu6RXHIvX7uFvn29hWLfmxMdEOh1JTsG48yyzpKQkm5KS4rbHE/Ekay2ZuUdYkZnPiqx8Vm47wNHSCuoZSGzdmLM6x3Bmpxj6tG1CaLB/zAbuKyxm+D++ISGuIXMnDKKej606CUTGmFRrbdKJPqcRtASU3MPFVSPkrHy+zcpn/+ESANrHRDK6b2uGdophcMdov72ecotG4Tw0shsPLFrH26uyGTuondOR5FeooMWvHS0pZ+X2A6zIPMCKrDy27j8CQNPIUIZ0jOaszjEM7RRD6yYRDif1nqv6t+HDjD1M/2Qz55/RjLjG9Z2OJCehgha/Ul5RSXpOId9m5bMiM5812Qcpr7SEBddjQPumjKkeJSe0bBiwv94bY5g+OpELn1nGH95bx8wb+teZ1SiBRgUtdZq1lm35R/k2K5/lmfl8/8MBikrKMQZ6xDXi5rM6cFbnGPq1a1LnVl54UpumEdx/cVce/3Aji9bsZky/1k5HkhNQQUudk3+khG+r55BXZOazp7AYgNZN6jOqV0uGdophSMcYmkaGOpzUt40fHM+SjL38aclGzuoSQ7MG4U5Hkl9QQYvPO15awaodBT+NkjftPQxAw/BghnaKYfJ5MZzVOYa2TSP0q3oN1KtneGpMIr+ZsZzHPtjAi2P7OR1JfsGlgjbG7ACKgAqg/GRLQkTcyVrL9E83M3PFDkorKgkJMvRr14T7LurK0E4x9GzVyOcuTlTXdGoWxZQLOvP0Z1v4ZN1eRvRs6XQk+ZmajKDPs9bmeyyJyC8kL9vGy99s45JecYzu24oB7ZsSEapf+txt4tkd+HjdXh75YAODOkTTRFNDPsM/Vt+L3/l43V6e/GQzIxNb8uxVvTm3azOVs4eEBNXjr1ckcuhYKU8s2eh0HPkZVwvaAp8bY1KNMRNPdqe8vDySkpJ+uiUnJ7snpQSUNdkHmTp/LX3bNubvV/YK2OVw3tQ9rhG3ntuRRWm7+dfmXKfjSDWXtnobY+KstXuMMc2AL4A7rLXLfnk/bfWW2so+cIzLX/iWqPBgFt06hOioMKcjBYyS8gpGzVjBkZJyPp96Ng3C/XM3pa/5ta3eLo2grbV7qv/MBd4DBrgvnkiVQ8dKueGNVVRYy8wb+qucvSwsOIinrkhk3+Fipn+y2ek4ggsFbYyJNMY0+PFt4EJgvaeDSWApKa/glrdSySk4TvK4JDrERjkdKSD1bduEm4a2Z87KbL774YDTcQKeKyPo5sAKY0w6sAr4yFr7qWdjSSCx1vLAu+tYub2Ap69MZEB73zw+KlDcc2FX2kVH8MCiDI6XVjgdJ6CdsqCttdustb2qb92ttX/xRjAJHM98mcl7abu5Z3gXLu3dyuk4Aa9+aBDTRyey88Ax/v75FqfjBDQtsxNHLUzNYcZXmVzRrzW3n9/J6ThSbXDHaK4d2JbXv93OmuyDTscJWCpoccy/f8jnwUUZDO0UzbTLe2qbto95cMQZNG8Yzv0LMygp11SHE1TQ4ois3CJueSuV+OhIXriun9+cWOJPGoSHMG10T7Jyj/D80iyn4wQk/VSI1+UVlXDDzNWEBQcx88b+fnt6iT84r2szRvdpxYtf/8DGPYedjhNwVNDiVcdLK7j5zRTyj5Tw2vikgDrJpK56ZFQCjSNCuP/ddMorKp2OE1BU0OI1FZWWu+ankZFziBlX96FXm8ZORxIXNIkM5YlLe7B+92GSl29zOk5AUUGL1zz58SY+27Cfh0cmcGH3Fk7HkRoY0bMlI3q04NkvM8nKPeJ0nIChghaveOu7Hby6Yjs3DInn90PjnY4jp+HxS7tTPySI/3k3g4rKU1/DR2pPBS0et3Tzfh5bvIFh3ZrxyKgELaero5o1COfRUQmk7jzIm9/tcDpOQNAFdsWj1u8u5Pa300iIa8hzV/fRCSh13Oi+rfgwYw9PfbqZzzfsdzqOz2hUP4SXxrn/yDAVtHjMnkPHuWnWahrXD+H18f2JDNO3W11njOHJ0T159IMNFB4rczqOz6hw4bLNp0M/MeIRRcVl/P6N1RwtqWDhrYNp1lAnRvuLlo3q88r1OpbUG1TQ4nZlFZXc9nYamblHmHlDf85o0dDpSCJ1kl4kFLey1vLoBxtYtjWPv1zWg7O7xDodSaTOUkGLW728bBtzV2Uz+dyOXD2grdNxROo0FbS4zUcZe5n+yWZGJbbk3gu7Oh1HpM5TQYtbpO4sYOqCtSS1a8LfdBK3iFuooKXWdh44yoQ3U4lrFE7y9UmEhwQ5HUnEL6igpVYOHi3lxpmrsdYy88YBNI0MdTqSiN/QMjs5bT+dxH3wOHMmDKR9TKTTkUT8igpaTou1lvsXZrBqRwEzrulD/3idxC3ibprikNPyzBdb+WDtHu67qCuX9IpzOo6IX1JBS429k7KLGUuzuCqpDZPP7eh0HBG/pYKWGvk2K58HF63jzE4x/PnyHrp0qIgHqaDFZZn7i5g0O5UOsZG8MLYvIUH69hHxJP2EiUtyi4q5YeZqwkOCmHnjABqG6yRuEU9TQcspHS+tYMKsFAqOlvL6+P60alzf6UgiAUHL7ORXVVRapsxLI2N3IcnjkujZupHTkUQChkbQ8qumfbyJzzfu59FRCQxPaO50HJGAooKWk5r17x28Vn0S941D2zsdRyTgqKDlhJZu3s/jH25gWLfmPDIqwek4IgHJ5YI2xgQZY9KMMUs8GUict3V/EXdUn8Q945reOolbxCE1GUFPATZ5Koj4hoKjpdw8K4WIsGBeuT6JiFC9jiziFJcK2hjTGhgJvOrZOOKksopKJs9JZd/hYpLH9aNlIy2nE3GSqyPoZ4H7gUoPZhEHWWt5bPEGvt9WwFNjetKnbROnI4kEvFMWtDFmFJBrrU091X3z8vJISkr66ZacnOyWkOJ5b32/k7dXZjPpnI5c3qe103FEBNc2qgwFLjHG/AYIBxoaY2Zba8f+8o6xsbGkpKS4O6N42LdZ+Tz+4UaGdWvGfRfpsFcRX3HKEbS19kFrbWtrbTxwNbD0ROUsddP2/KNMnrOGjrGRPHt1H63YEPEhWgcdwA4Xl3HzrNXUM/Dq9f2JCtOKDRFfUqOfSGvt18DXHkkiXlVRabnj7TR2HjjG7JsH0jY6wulIIvILGjIFqOmfbOKbrXlMu7wngzpEOx1HRE5AUxwBaEHKLl5ZXnWNjWsHtnU6joichAo6wKTsKOCh96qOrHp4ZDen44jIr1BBB5Ccg8eYNDuV1k0i+N9r+xKsI6tEfJrmoAPE0ZJyJryZSkl5JfOuT6JRhI6sEvF1GkIFgMpKyz0L0tmy7zD/vKYPnZpFOR1JRFyggg4Az365lU837OMPv+nGuV2bOR1HRFykgvZzH6bvYcbSLH6X1JqbztSpKCJ1iQraj63LKeTed9JJateEJy7rgTHaxi1Sl6ig/VTu4WImvJlCTFQYL43rR1hwkNORRKSGtIrDDxWXVTDhrVQOF5excNIQYqLCnI4kIqdBBe1nrLU8uGgd6bsO8dLYfiTENXQ6koicJk1x+JmXvtnGe2m7uWd4Fy7u0cLpOCJSCypoP/Llxv389bPN/LZXHLef38npOCJSSypoP7FlXxFT5qXRI64Rfx2TqBUbIn5ABe0HCo6WctOs1USGBfPK9UnUD9WKDRF/oBcJ67jS8kpunZ1KblEJC24ZTItG4U5HEhE30Qi6DrPW8tjiDazcXsDTVyTSu01jpyOJiBupoOuwN7/bydxV2Uw+tyOX9m7ldBwRcTMVdB21IjOfPy3ZyLBuzbn3wq5OxxERD1BB10Hb848yeU4qnWKjePbq3tSrpxUbIv5IBV3HFB4v46ZZqwkOqser45OICtPrvCL+SgVdh5RXVHLH3DSyDxzjxev60qZphNORRMSDNPyqQ578ZDPLtubx5OieDOwQ7XQcEfEwjaDriAWrd/Haiu3cMCSeawa0dTqOiHiBCroOWL2jgIfeX8dZnWN4eGQ3p+OIiJeooH1czsFjTHorldZNInj+mr4EB+mfTCRQ6Kfdhx0tKefmWSmUVlTy6vgkGkWEOB1JRLxIBe2jKistU+evZev+Ip6/ti8dY6OcjiQiXqaC9lHPfLmVzzfu5+GRCZzTJdbpOCLiABW0D1qcvod/Ls3iqqQ23Dg03uk4IuIQFbSPWZdTyH3vpDMgvilPXNZDF94XCWCnLGhjTLgxZpUxJt0Ys8EY87g3ggWiQ8dKmTQ7lejIUF4Y25fQYP3/UySQubKTsAQ431p7xBgTAqwwxnxirf3ew9kCSmWl5a75a8krKmHBpMHERIU5HUlEHHbKIZqtcqT63ZDqm/VoqgA0Y2kmX2/J49HfJujC+yICuDgHbYwJMsasBXKBL6y1K090v7y8PJKSkn66JScnuzOr3/p6Sy7PfZXJ6L6tuG6gtnGLSBWXLpZkra0AehtjGgPvGWN6WGvX//J+sbGxpKSkuDujX9tVcIwp89bStXkD/nJZT70oKCI/qdGrUNbaQ8DXwMUeSRNgissqmDxnDZXW8vK4fjqNW0T+gyurOGKrR84YY+oDw4DNng4WCP64eAPrdhfyzO960y460uk4IuJjXJniaAnMMsYEUVXoC6y1Szwby/8tWL2Leat3cdt5HRmW0NzpOCLig05Z0NbaDKCPF7IEjPW7C3n4g/Wc2SmGu4frwFcROTHthPCyHzejxESG8tzVvQnSga8ichI68sqLfrxC3f7DxSy4ZTDR2owiIr9CI2gv+ufSLP61JY9Hf9udPm2bOB1HRHycCtpLvt6Sy7NfbWV0n1aM1WYUEXGBCtoLdhUc46751ZtRLtdmFBFxjQraw37cjFJRYXlxrDajiIjr9CKhhz3+YdVmlORx/Wgfo80oIuI6jaA9aEHKLuau2sWt53bkwu4tnI4jInWMCtpD1u8u5JH31zOkYzT3DO/idBwRqYNU0B5QeKyMW+ek0iQilBnX9CE4SP+ZRaTmNAftZlUno6Sxr7CY+bfoZBQROX0a2rnZ8/+q2ozyyKgE+moziojUggrajb7ZmsczX27lst5xjBvUzuk4IlLHqaDdJOfgMabMS6NLswZMG63NKCJSeypoN/j5ZpSXxvUjIlRT+yJSe2oSN/jTko1k5BTysjajiIgbaQRdSwtTc3h7ZTaTzunIRdqMIiJupIKuhQ17CnnovXUM7hDNvRdqM4qIuJcK+jQVHivj1tlraBIRyj+v1WYUEXE/zUGfhspKy90L1rK38DjzJmozioh4hoZ9p+GFr7P4anMuD49MoF87bUYREc9QQdfQ8sw8/v7FVi7tHcf1g7UZRUQ8RwVdA7sPHefOuVWbUZ7UZhQR8TAVtItKyiuYPDuV8grLi2P7ajOKiHicWsZFf/pwI+k5hbw0th8dYqOcjiMiAUAjaBe8m5rDnJXZ3HJOBy7uoc0oIuIdKuhT2LjnMH+o3oxy34VdnY4jIgFEBf0rCo9XnYzSOCJEJ6OIiNdpDvokKist9yxYy+6Dx5l/yyBiG2gzioh4l4aEJ/HiNz/w5aZcHh7ZjX7tmjodR0QCkAr6BFZk5vP3z7dwSa84xg+JdzqOiASoUxa0MaaNMeZfxphNxpgNxpgp3gjmlD2HjnPnvDQ6NYti+hhtRhER57gyB10O3GOtXWOMaQCkGmO+sNZu9HA2ryspr+DWOWsoLa/kxbE6GUVEnHXKEbS1dq+1dk3120XAJqCVp4M54YklG0nfdYinr0ikozajiIjDajRENMbEA32AlSf6fF5eHklJST+9P3HiRCZOnFiLeN4zb1U2s7/PZuLZHRjRs6XTcUREXC9oY0wU8C5wl7X28InuExsbS0pKiruyec13Pxzg4ffXc3aXWO6/SJtRRMQ3uLSKwxgTQlU5z7HWLvJsJO/akX+UW+ekEh8TyfM6GUVEfIgrqzgM8BqwyVr7D89H8p7C42XcNGs1AK+NT6JheIjDiURE/p8rw8WhwDjgfGPM2urbbzycy+PKKyq5/e01ZBcc46Wx/WgXHel0JBGR/3DKOWhr7QrA7xYD//mjTSzPzOepMT0Z1CHa6TgiIv8lICdc3/p+J2/8ewcTzmrPVf3bOh1HROSEAq6gV2Tm88fFGzj/jGY8MKKb03FERE4qoAp6W94RJs9JpVNsFM9d3Zugen43cyMifiRgCrrwWBk3zUohJKger45PooFWbIiIjwuIi02UVVQy+e1Udh88zpwJA2nTNMLpSCIip+T3BW2t5Y+LN/Bt1gH+dmUv+sfr2s4iUjf4/RTHm9/tZM7KbCad05Er+rV2Oo6IiMv8uqCXbc3j8Q83MDyhua6xISJ1jt8WdFZuEbfNWUPXFg159qre1NOKDRGpY/yyoA8eLeX3b6QQFlK1YiMyzO+n2kXED/ldc5WWVzJpdir7Dhczd8IgWjWu73QkEZHT4lcjaGstj7y/npXbC/jrmET6tWvidCQRkdPmVwX92ortzE/ZxR3nd+KyPn55KpeIBBC/Keilm/cz7eNNjOjRgqnDujgdR0Sk1vyioLfsK+LOuWtJiGvI33/XSys2RMQv1PmCPnCkhJtmrSYiNIhXrk8iItTvXvcUkQBVp9uspLyCSbNTySsqYcEtg2nZSCs2RMR/1NmCttby0HvrWb3jIM9f24debRo7HUlExK3q7BRH8rJtLEzN4a5hnRmVGOd0HBERt6uTBf3Fxv1M/3QzoxJbMuWCzk7HERHxiDpX0Jv2HmbKvDQSWzXib1f2whit2BAR/1SnCjqvqISbZ6XQMDyE5OuTCA8JcjqSiIjH1JkXCYvLKrjlrRQOHC1h4aQhNG8Y7nQkERGPqhMFba3lwUXrWJN9iBev60uPVo2cjiQi4nF1Yorjha9/4L203dx7YRdG9GzpdBwREa/w+YL+dP1env5sC5f1juO28zo5HUdExGt8uqDX7y5k6vx0+rRtzPQxiVqxISIBxWcLOvdwMTfPSqFJRAjJ47RiQ0QCj0++SFhcVsGEN1M4XFzGwklDiG0Q5nQkERGv87mCttZy7zvpZOwuJHlcEglxDZ2OJCLiCJ+b4pjxVRZLMvbyPxefwfCE5k7HERFxjE8V9JKMPTzz5VbG9G3NLWd3cDqOiIijTlnQxpjXjTG5xpj1ngySvusQ9yxIJ6ldE6aN7qEVGyIS8FwZQb8BXOzJEPsKi5nwZgqxDcJ4eVw/woK1YkNE5JQFba1dBhR4KsDx0gpufnM1R0vKeW18f6KjtGJDRATcvIojLy+PpKSkn96fOHEiEydO/NW/Ywx0btaAu4d3oWuLBu6MIyJSp7m1oGNjY0lJSanR3wkPCeKZq3q7M4aIiF/wqVUcIiLy/1TQIiI+ypVldnOB74CuxpgcY8xNno8lIiKnnIO21l7jjSAiIvKfNMUhIuKjfKagk5OTnY7gdXrO/i/Qni/oObuTCtpBes7+L9CeL+g5u5PPFLSIiPwnY61134MZkwfsPM2/HgPkuy1M3aDn7P8C7fmCnnNNtbPWxp7oE24taBERcR9NcYiI+CgVtIiIj1JBi4j4KBW0iIiPcrygjTEXG2O2GGOyjDEPOJ3H04wxbYwx/zLGbDLGbDDGTHE6k7cYY4KMMWnGmCVOZ/EGY0xjY8xCY8zm6n/vwU5n8jRjzNTq7+v1xpi5xphwpzO524mOATTGNDXGfGGMyaz+s4k7vpajBW2MCQL+FxgBJADXGGMSnMzkBeXAPdbabsAg4LYAeM4/mgJscjqEFz0HfGqtPQPohZ8/d2NMK+BOIMla2wMIAq52NpVHvMF/HwP4APCVtbYz8FX1+7Xm9Ah6AJBlrd1mrS0F5gGXOpzJo6y1e621a6rfLqLqh7aVs6k8zxjTGhgJvOp0Fm8wxjQEzgZeA7DWllprDzmbyiuCgfrGmGAgAtjjcB63O8kxgJcCs6rfngVc5o6v5XRBtwJ2/ez9HAKgrH5kjIkH+gArnU3iFc8C9wOVTgfxkg5AHjCzelrnVWNMpNOhPMlauxv4G5AN7AUKrbWfO5vKa5pba/dC1SAMaOaOB3W6oM0JPhYQO2eMMVHAu8Bd1trDTufxJGPMKCDXWpvqdBYvCgb6Ai9aa/sAR3HTr72+qnre9VKgPRAHRBpjxjqbqm5zuqBzgDY/e781fvgr0S8ZY0KoKuc51tpFTufxgqHAJcaYHVRNY51vjJntbCSPywFyrLU//na0kKrC9mfDgO3W2jxrbRmwCBjicCZv2W+MaQlQ/WeuOx7U6YJeDXQ2xrQ3xoRS9YLCYoczeZQxxlA1L7nJWvsPp/N4g7X2QWtta2ttPFX/xkuttX49srLW7gN2GWO6Vn/oAmCjg5G8IRsYZIyJqP4+vwA/f2H0ZxYD46vfHg984I4Hdeup3jVlrS03xtwOfEbVK76vW2s3OJnJC4YC44B1xpi11R/7g7X2YwcziWfcAcypHnxsA250OI9HWWtXGmMWAmuoWq2UBvjdtUerjwE8F4gxxuQAjwHTgQXVRwJmA1e65WvpYkkiIr7J6SkOERE5CRW0iIiPUkGLiPgoFbSIiI9SQYuI+CgVtIiIj1JBi4j4qP8DH9CSfq1OItcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, estimates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x211781a7908>]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD2CAYAAADLcgxzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUZdr/8c+VAoEAoSRIqBFCB2kRVBQVGyLiqlhWESyIbRVdd32U37Pq+qyru+raXaWooK511UXFChYUBUNXQwmhtyS0kEBCyv37IwMbMJAJTHKmfN+v17ymnJOZy9F8vXPuc5/LnHOIiEjoivK6ABEROToKchGREKcgFxEJcQpyEZEQpyAXEQlxMbX9gYmJiS4lJaW2P1ZEJKTNmzcv1zmXVNm2Wg/ylJQU0tPTa/tjRURCmpmtOdQ2HVoREQlxCnIRkRCnIBcRCXEKchGREKcgFxEJcQpyEZEQpyAXEQlxtX4euUgw2ppfxHsLNtCoXiwdkuJpn9iAJvF1vC5LxC8KcoloeYXFTPomi8nfrqJgb+kB25rUj6V9UoPyYE9qQPvE8vt2zeoTG60/ZiV4KMglIu3ZW8qU71fzz69WsnNPMef1TGbcmR2pEx1FVm4+WTkFrMzJZ2VOATOX5vBW+vr9PxsdZbRtWp/2ifF0aP7fgG+fFE+z+DqYmXf/YBKRFOQSUfaWlPFm+jqenrGC7F1FnNY5iT+c3ZkerRL275OSGM/gLgf+3M49xazKLSAr578hn5VTwKzMXPaWlO3fr1FczP5Q71BhNN+uWX3qxkTX1j+mRBgFuUSE0jLHfxZu4PEvlrNu2x6OT2nCM1f0pf+xTf36+YR6sfRu05jebRr/6n037tizP9izcvNZmV3Ad5m5vDt/w/79ogxaN6lPe9/x9w7NffdJ8SQ1rKtRvBwVBbmENeccn/68hcc+W8aK7Hy6t2zES9f04LROSQEJz+goo03T+rRpWp/TOh+4Lb+ohFX7wj2nfDS/MqeAH7K2Ulj831F8g7oxvoCP9x2TLx/RH5sYT1ysRvFSNQW5hCXnHN9m5vLIp8tYvH4n7ZPiefaKvpzbowVRUbUz+m1QN4aerRPo2TrhgNfLyhyb8gr3H6bJysknK7eAuau28f7Cjfv3M4OWCfVIbd6AsYPaMzA1sVbqltBjzrla/cC0tDSny9hKTZq3ZjuPfLqUH7K20apxPcad2ZGL+rQiJgTONNm9t8R3LL5g/6Ga9NXb2bBjD9efcix/OKezjrVHKDOb55xLq2ybRuQSNjI25fHop8uYsTSbxAZ1uP/8bvx2QNuQCr76dWLo3jKB7i3/O4rfs7eUv07PYOKsVcxakcuTl/ehc4uGHlYpwUYjcgl5q3ILePzz5UxbtJFGcTHccGoHrhmYQv064TVOmbl0C3e9s5i8whLuHtKFq09KqbXDROK9w43IFeQSsjbu2MPTM1fwVvp66kRHcc3AFG4Y1IGE+rFel1ZjcvOL+J93FjNjaTandEzk0Ut6cUyjOK/LklqgIJewsjW/iOe+WskrP6zBOceVA9px8+kdaN4wMgLNOce/5q7l/z78hbjYaB6+qCdDeiR7XZbUMB0jl7BQcTn9nuJSLu7bmnFndqR1k/pel1arzIwrB7TjhPbNuP2Nhdz46nwuTWvNved3p0Fd/UpHIv1bl6BX2XL6O87qRGrzBl6X5qkOSQ34900n8eSM5Tz31UrmrNrG45f1pm/bJl6XJrXMr0MrZrYa2AWUAiUHD++tfGXFk8BQYDdwtXNufmXvpUMr4q+9JWW8+eNanpqZSc4hltNLubmrtnHHmwvZnFfI705P5dbBqSFxuqX4L1CHVk53zuUeYtu5QEffbQDwT9+9SLWVljneX1C+nH799vLl9M9WYzl9JOp/bFM+vv0U7vvPzzw5YwVfL8/hict6k5IY73VpUgsCdWjlAmCqKx/e/2Bmjc0s2Tm3KUDvLxGgfDn9Zh77bPn+5fR/+U0PTg3Qcvpw1ygulscv683pXZrzv+8tYehTs7j//O5cktZa31+Y8zfIHfCZmTngBefchIO2twLWVXi+3vfar4I8JyeHtLT//nUwduxYxo4dW62iJbw455i1IpdHP/vvcvrnruzLkO61t5w+nAzv1ZK0dk34/VsLuevfi5m5NJuHLuqpRhlhzN8gH+ic22hmzYHPzWypc+6bCtsr+22r9OB7UlISOkYu+8xbs42/f7KMOavKl9M/MuI4LgyR5fTBrGXjevxrzAlMnJXFo58t45wntvPoJb0Y1CnJ69KkBvgV5M65jb77bDN7D+gPVAzy9UCbCs9bAxsROYRfNubx2GehvZw+2EVFGTec2oGTOyYy7o2FjHpxLtcMTOF/hnTRVRXDTJVBbmbxQJRzbpfv8dnAAwftNg34nZm9Qfkk504dH5fK5BeVMP7dJfuX0//xnM5huZw+mHRvmcCHt57Mwx8v5aXvVvNdZvn1WromN/K6NAkQf357jgHe802WxAD/cs59YmY3AjjnngemU37qYSblpx9eUzPlSigrLXOMe30BXy3P4ebTOoT9cvpgEhcbzf3Du3Nq5yT++PZiLnjmO+4a0plrBx6reYgwoCX6Umse/OgXJs5axf9d0J2rTkzxupyItTW/iLvfXcLnv2xhYGozHr2kF8kJ9bwuS6pwuPPINaMkteL1uWuZOGsVo09spxD3WLMGdZlwVT8euqgn89fsYMgTs/hosY6EhjIFudS42Stz+dP7PzGoUxJ/GtbN63KE8uu1/LZ/W6aPO4WUZvW55V/zufOtRewqLPa6NDkCCnKpUVk5+dz06nyOTYznmSv66LTCIHNsYjzv3HQStw1O5b0F6xn61CzSV2/zuiypJv1WSY3ZsXsv101JJzrKmDz6eBrFaWIzGMVGR/H7szvz1g0nAnDpC9/z2GfLKC4tq+InJVgoyKVGFJeWcdOr89mwfQ8vXNWPts0i61KzoSgtpSnTbzuFC/u05umZmYz452yycvK9Lkv8oCCXgHPO8af3f+L7rK08dFFPjk/Rxa5CRcO4WB67tBfPXtGX1Vt3c95T3/KvOWup7bPbpHoU5BJwk79dxRs/ruOW0ztwcb/WXpcjR+C845L55PZT6NuuMePfW8L1U+exNb/I67LkEBTkElBf/LKFB6dncG6PFtx5Vmevy5GjkJxQj1euHcD/nteVb5bncM4Ts/hyWbbXZUklFOQSMBmb8hj3xgJ6tEzgsUt7acVgGIiKMsac0p73bxlI0/hYrnnpR+77z08UFpd6XZpUoCCXgMjeVciYKek0iIth4qg0XTslzHRr2YhpvzuZawamMOX7NQx7+lt+2rDT67LER0EuR62wuJSxU+exrWAvk0cfT4uEyOhmH2niYqO57/zuTL22P3l7irnwue94f8EGr8sSFORylJxz/PGdxSxct4PHL+utfpoRYFCnJD69fRD92jXhzrcX8clPm70uKeIpyOWoPDljBR8s2shdQzozpEcLr8uRWtIkvg6TRx9Pz1YJ3Pb6AmatyPG6pIimIJcjNm3RRp74YgUX923NTad28LocqWXxdWN4+ZrjaZ8Uz9ip87S030MKcjki89du5w9vL6J/SlP+elEPNfeNUI3r1+GV6wbQIiGOa176UROgHlGQS7Vt2LGHsVPn0aJRHM9f1U/t2SJcUsO6vDpmAA3jYhj14lwys3d5XVLEUZBLteQXlXDdyz9SVFzK5NFpNFVndgFaNa7Hq2MGEGUwctJc1m3b7XVJEUVBLn7b16ptRXY+z1zZl47HNPS6JAki7ZMa8Mp1A9i9t4SRk+eQnVfodUkRQ0EufntoegYzlmZz//ndOLVTktflSBDqmtyIl6/tT86uIkZOnsP2gr1elxQRFOTil9fnrmXSt2rVJlXr27YJk0alsXrrbka/NFddh2qBglyqpFZtUl0npSby3BV9+XljHtdNSWfPXl2bpSb5HeRmFm1mC8zsw0q2XW1mOWa20HcbE9gyxStq1SZH6sxux/CPS3vx4+pt3PTaPPaWqONQTanOb+U4IOMw2990zvX23SYdZV0SBNSqTY7WBb1b8eBvevLVshzueHMhpWVqUFET/ApyM2sNnAcooCOEWrVJoFwxoC3jh3bhoyWbuOfdxZQpzAPO32uNPgHcBRzufLOLzWwQsBy4wzm37miLE29UbNX22CW91KpNjtrYQR3ILyzhqZmZNKgby5+GddVq4ACqckRuZsOAbOfcvMPs9gGQ4pw7DvgCmHKoHXNyckhLS9t/mzBhQrWLlpqlVm1SE+44qxNXn5TCi9+t4okvVnhdTljxZ0Q+EBhuZkOBOKCRmb3qnBu5bwfn3NYK+08E/naoN0tKSiI9Pf1I65UaplZtUlPMjHuHdaOgqIQnZ6ygYVwMY05p73VZYaHKEblz7h7nXGvnXApwOTCzYogDmFlyhafDOfykqAQptWqTmhYVZTx88XEM7dmCv3yUwetz13pdUlg44n5cZvYAkO6cmwbcZmbDgRJgG3B1YMqT2qJWbVJboqOMJy7rQ0FROuPfW0J83RiG92rpdVkhzZyr3RnktLQ0p0MrwaWwuJTLJ/zAss27ePvGE9XlR2rFnr2ljH5xLvPXbueFq/pxRtdjvC4pqJnZPOdcWmXbtLojwh3Yqq2XQlxqTb060Uy+Oo1uLRtx02vzmb0y1+uSQpaCPMId2KotueofEAmghnGxTLmmP+2a1uf6KeksWLvd65JCkoI8gqlVmwSDJvF1eHXMAJo1qMvVL/3I0s15XpcUchTkEUqt2iSYHNMojtfGDKBebDQjJ81lVW6B1yWFFAV5BFKrNglGbZrW59Ux/SlzjpGT5rBxxx6vSwoZCvIIo1ZtEsxSmzdk6rX9ydtTzMhJc8jZVeR1SSFBQR5B1KpNQkGPVgm8eM3xbNy5h1EvzmXnbjWmqIqCPIKoVZuEiuNTmjLhqjRWZudz9ctzKSgq8bqkoKYgjxBq1SahZlCnJJ76bW8WrdvB2FfSKSxWl6FDUZBHALVqk1A1pEcyj4zoxXeZW7n19QUUl6rLUGUU5GFOrdok1F3crzUPXNCdz3/Zwh/fXqTGFJXQlZHCmFq1SbgYdWIKuwpLeOTTZcTXjeEvv9Hah4oU5GGqYqu2164foFZtEvJuOT2VXYUlPP/1ShrExXD3kC4Kcx8FeZh68dtVfJ+1lUfVqk3CyP8M6Ux+UTEvfJ1Fo7hYbjk91euSgoKCPAxt2rmHJ2es4MyuzRmhVm0SRsyMB4b3oKCotPwwS51orh54rNdleU5BHoYe/CiD0jLHfed397oUkYCLijIeGXEcBUUl3P/BLzSIi434AYtOYQgz32Xm8uHiTdx8Wiptmuq4uISnmOgonr6iDyenJnLXO4v4eMkmr0vylII8jOwtKeO+aT/Ttml9bjhVTW0lvNWNiWbCqH70aduE295YwNfLc7wuyTMK8jDy0neryMzO5/7h3YiL1RUNJfzVrxPDi1cfT8fmDbnhlXR+XL3N65I8oSAPExUnOAd3Ue9DiRwJ9WKZel1/Wjaux7Uv/chPG3Z6XVKtU5CHCU1wSiRLbFCX18YMoFG9WK6aPIfM7F1el1SrFORhYLZvgvOm0zpoglMiVnJCPV4bM4DoqCiufTmd7QV7vS6p1vgd5GYWbWYLzOzDSrbVNbM3zSzTzOaYWUogi5RD21tSxr2+Cc4b1XdTIlxKYjwTRvVjc14hN782P2IuslWdEfk4IOMQ264DtjvnUoHHgb8dbWHiH01wihyob9smPHxRT77P2sqfP/jZ63JqhV9BbmatgfOASYfY5QJgiu/xO8AZposg1DhNcIpU7qK+rbnh1Pa8+sNaXvl+tdfl1Dh/R+RPAHcBh/o7pRWwDsA5VwLsBJpVtmNOTg5paWn7bxMmTKhmybKPJjhFDu2uc7owuEtz7v/gF2Zn5npdTo2qcom+mQ0Dsp1z88zstEPtVslrlV40OCkpifT0dP8rlErtm+C8/cyOmuAUqUR0lPHk5b256LnZ3PTafP5zy0BSEuO9LqtG+DMiHwgMN7PVwBvAYDN79aB91gNtAMwsBkgAIvPM/FqgCU4R/zSMi2Xy6OMxgzFT08krDM9GzlUGuXPuHudca+dcCnA5MNM5N/Kg3aYBo32PR/j2URuPGvLy7PIJzvvO1wSnSFXaNqvPc1f2ZXVuAeNeX0BpGHYYOuLzyM3sATMb7ns6GWhmZpnA74G7A1Gc/NqmnXt44ovyCc4zumqCU8QfJ3VI5P7h3flyWQ5/+2Sp1+UEXLUuY+uc+wr4yvf43gqvFwKXBLIwqZwmOEWOzMgT2rF8yy4mfJNFp2MahtWlb7WyM4RoBafI0fnTsG4MTG3G+HeXMG9N+EzjKchDhCY4RY5ebHQUz17Rl5aN47jhlXls2LHH65ICQkEeIjTBKRIYjevXYdLoNIqKy7h+Sjq795Z4XdJRU5CHgM07CzXBKRJAqc0b8tQVfVi6OY8731pEWYifyaIgDwEPTs+gpMxx7zBNcIoEyumdmzN+aFc+/mkzT85Y4XU5R0VBHuRmZ+bywaKN3HxaB9o20wSnSCBdd/KxjOjXmidnrOCjxaHb91NBHsQ0wSlSs8yMBy/sQb92Tbjz7YUh211IQR7ENMEpUvPqxkTz/Mh+NK1fh+unppO9q9DrkqpNQR6kNMEpUnuSGtZl4ug0duwu5oZX5lFYXOp1SdWiIA9SmuAUqV3dWybw+GW9WLB2B+PfXUIoXS5KQR6EZq/UBKeIF4b0SOb3Z3Xi3QUbmPBNltfl+K1a11qRmldcWsa9//mZNk3raYJTxAO3Dk5l2ZZdPPzJUlKbNwiJQ5sakQeZ/T04z++uCU4RD5gZj47oRfeWjRj3xkKWb9nldUlVUpAHEU1wigSHenWimTgqjXp1ohkzJZ3tBXu9LumwFORBRBOcIsEjOaEeE67qx+a8Qm56bR7FpYdqWew9BXmQ0ASnSPDp07YJf7u4Jz9kbeP+aT97Xc4habIzCGiCUyR4XdinNcs25/P81yvp0qIhV52Y4nVJv6IReRDQBKdIcPvjOZ05o0tz7v/gF77LzPW6nF9RkHts885CnvxiBWd00QSnSLCKjjKeuLw3HZLiufm1+azOLfC6pAMoyD324PQMitWDUyToNYyLZdKo44kyuG7Kj+QVFntd0n4Kcg/tm+C86VRNcIqEgrbN6vPclf1Ys3U3t/5rAaVB0pCiyiA3szgzm2tmi8zsZzP7cyX7XG1mOWa20HcbUzPlho+KE5w3naYJTpFQcWKHZvz5gu58vTyHhz/O8LocwL+zVoqAwc65fDOLBb41s4+dcz8ctN+bzrnfBb7E8LRvgnPy6DRNcIqEmCsHtGP55l1MnLWKTsc05JK0Np7WU+WI3JXL9z2N9d2C4++JEKUJTpHQ96dh3RiY2oz/995PpK/e5mktfh0jN7NoM1sIZAOfO+fmVLLbxWa22MzeMbND/u8pJyeHtLS0/bcJEyYcYemhSxOcIqEvJjqKZ6/oS8vGcdz46jw27NjjWS1WnWvumllj4D3gVufcTxVebwbkO+eKzOxG4FLn3ODK3iMtLc2lp6cfZdmha/bKXK6YOIdxZ3TkjrM6eV2OiBylzOx8Lnz2O1o3rc87N55IfN2aWWdpZvOcc2mVbavWWSvOuR3AV8CQg17f6pwr8j2dCPQ7gjrDXnFpGfdpglMkrKQ2b8BTV/Rh2eY87nxrEWUenMniz1krSb6ROGZWDzgTWHrQPskVng4HgmMqN8i8/N1qVmTnc98wreAUCSend27O+KFd+eTnzTwxY0Wtf74/fwMkA1PMLJry4H/LOfehmT0ApDvnpgG3mdlwoATYBlxdUwWHqvJL1C7njC7NObObJjhFws11Jx/Lss27eGrGCjod04Bhx7Wstc+uMsidc4uBPpW8fm+Fx/cA9wS2tPCiCU6R8GZm/OXCHmTlFvCHtxfRrmk8PVsn1Mpna2VnLdAKTpHIUDcmmudH9qNZfF2un5pOdl5hrXyugryGaYJTJLIkNazLxFFp7NxTzNhX5lFYXFrjn6kgr2Ga4BSJPN1aNuLxy3qzcN0O7nl3CdU5zftIKMhr0JY8TXCKRKohPVpw51mdeG/BBl74JqtGP0sdgmrQgx9pglMkkv1ucCrLtuzib58sJTWpQY0N6DQiryGzV+YyTROcIhHNzHhkRC96tExg3BsLWLZ5V418joK8BmiCU0T2qVcnmomj0qhfN4Z/z19fI5+hQys1YN8E56RRukStiECLhDim/W4gLRrF1cj7K8gDTBOcIlKZ5IR6NfbeOrQSYJrgFJHapiAPIE1wiogXFOQBoglOEfGKgjxAtIJTRLyiIA8ATXCKiJcU5AGgCU4R8ZKC/CjNydqqCU4R8ZSC/CiUlTkenJ5BckKcJjhFxDMK8qPw4ZJNLF6/kzvP7qwJThHxjIL8CBWVlPL3T5bSNbkRF/Zp5XU5IhLBFORH6JXv17B++x7GD+1CdJR5XY6IRDAF+RHYubuYp2dmMqhTEqd0TPK6HBGJcFUGuZnFmdlcM1tkZj+b2Z8r2aeumb1pZplmNsfMUmqi2GDxzJcryCss5p5zu3hdioiIXyPyImCwc64X0BsYYmYnHLTPdcB251wq8Djwt8CWGTzWbdvNlNlrGNG3NV2TG3ldjohI1UHuyuX7nsb6bgd3Er0AmOJ7/A5whpmF5YHjRz5dRlQU/P7sTl6XIiIC+HmM3MyizWwhkA187pybc9AurYB1AM65EmAn0Kyy98rJySEtLW3/bcKECUdefS1bvH4H0xZt5LqTj63RawuLiFSHX40lnHOlQG8zawy8Z2Y9nHM/VdilstH3waN2AJKSkkhPT69+pR5zzvHgRxk0i6/Djadq8Y+IBI9qnbXinNsBfAUMOWjTeqANgJnFAAnAtgDUFzRmLs1mzqptjDuzIw3jYr0uR0RkP3/OWknyjcQxs3rAmcDSg3abBoz2PR4BzHTOVToiD0UlpWU89PFS2ifG89v+bb0uR0TkAP4cWkkGpphZNOXB/5Zz7kMzewBId85NAyYDr5hZJuUj8ctrrGIPvJW+nszsfJ4f2Y/YaJ16LyLBpcogd84tBvpU8vq9FR4XApcEtrTgUFBUwj8+X05auyac013XGheR4KPhZRUmfJNFbn4R48/rSpieUSkiIU5BfhjZeYVM+CaL83om07dtE6/LERGplIL8MB7/YgUlZWX88ZzOXpciInJICvJDWLFlF2/+uJYrB7QjJTHe63JERA5JQX4ID3+8lPg6Mdx2RkevSxEROSwFeSW+X7mVGUuzufn0VJrG1/G6HBGRw1KQH6SszPHX6Rm0TIjjmoEpXpcjIlIlBflBPli8kSUbdvKHc9SHU0RCg4K8gsLiUv7+yTK6JTfiN73Vh1NEQoOCvIKp369mw449jB/alSj14RSREKEg99mxey/PzMzk1E5JnNwx0etyRET8piD3eWZmJvlFJdwzVH04RSS0KMgp78M59fs1jOjXmi4t1IdTREKLghz4+74+nGdpKb6IhJ6ID/JF63bwwaKNXH9Ke1okxHldjohItUV0kDvneHB6eR/OG9SHU0RCVEQH+RcZ2cxdtY3bz+xIg7p+9aEWEQk6ERvkJaVlPPxxBu0T47lcfThFJIRFbJC/mb6OlTkF/M+5XdSHU0RCWkQmWH5RCY9/voLjU5pwdjf14RSR0BaRQb6/D+dQ9eEUkdBXZZCbWRsz+9LMMszsZzMbV8k+p5nZTjNb6LvdWzPlHr0teYVM/CaL845Lpo/6cIpIGPDnVI0S4E7n3HwzawjMM7PPnXO/HLTfLOfcsMCXGFiPf76ckrIy7lIfThEJE1WOyJ1zm5xz832PdwEZQEhe43X5ll28lb6OkSe0o10z9eEUkfBQrWPkZpYC9AHmVLL5RDNbZGYfm1n3ANQWcA9NzyC+bgy3DVYfThEJH36vgjGzBsC/gdudc3kHbZ4PtHPO5ZvZUOB9oNK0zMnJIS0tbf/zsWPHMnbs2GoXXl2zM3P5clkOd5/bhSbqwykiYcSvIDezWMpD/DXn3LsHb68Y7M656Wb2nJklOudyD943KSmJ9PT0o6m52srKHH/9OINWjetx9UkptfrZIiI1zZ+zVgyYDGQ45/5xiH1a+PbDzPr73ndrIAs9GtMWbeSnDXn84ZxO6sMpImHHnxH5QOAqYImZLfS9Nh5oC+Ccex4YAdxkZiXAHuBy55yrgXqrrbC4lEc+XUaPVo24oFdIztGKiBxWlUHunPsWOOyqGefcM8AzgSoqkKbMLu/D+ciI49SHU0TCUliv7NxesJdnvszktM5JnJSqPpwiEp7COsifnplJQVEJ95zb1etSRERqTNgG+ZqtBbzyw2ou6deGzi0ael2OiEiNCdsgf+TTZcRERfH7szt5XYqISI0KyyBfuG4HHy7exPWnHMsxjdSHU0TCW9gFuXOOv36UQWKDOoxVH04RiQBhF+Sf/7KFuau3cfuZndSHU0QiQlgFeXFpGQ9/spT2SfFcdnwbr8sREakVYRXkb/y4jqycAu4eoj6cIhI5wibt8otKePKL5fRPacpZ6sMpIhEkbA4iv/D1SnLz9zJptPpwikhkCYsR+Za8QibOymLYccn0btPY63JERGpVWAT5Pz5bTmmZ465zunhdiohIrQv5IF+2eRdvz1vHqBNTaNusvtfliIjUupAP8oc+zqBB3RhuHZzqdSkiIp4I6SD/LjOXr5blcMvpqTSurz6cIhKZQjbIy8ocf51e3odztPpwikgEC9kgf3/hBn7emMcfz+msPpwiEtFCMsgLi0t51NeHc3ivll6XIyLiqZAM8pdnr2bjzkLGD+2qPpwiEvFCLsi3F+zl2S8zGdylOSd1UB9OEZGQC/KnZq6goKiEu8/V4h8REfAjyM2sjZl9aWYZZvazmY2rZB8zs6fMLNPMFptZ35oods3WAl79YQ2XprWh0zHqwykiAv5dNKsEuNM5N9/MGgLzzOxz59wvFfY5F+jouw0A/um7D6i/f+Lrw3mW+nCKiOxT5YjcObfJOTff93gXkAG0Omi3C4CprtwPQGMzSw5kofPXbuejJZu4flB7mqsPp4jIftU6Rm5mKUAfYM5Bm1oB6yo8X2r/OrQAAAOISURBVM+vwx6AnJwc0tLS9t8mTJjg32cDgzolccOg9tUpWUQk7Pl9PXIzawD8G7jdOZd38OZKfsRV9j5JSUmkp6f7X6FPn7ZNmHpt/2r/nIhIuPNrRG5msZSH+GvOuXcr2WU9ULFJZmtg49GXJyIiVfHnrBUDJgMZzrl/HGK3acAo39krJwA7nXObAliniIgcgj+HVgYCVwFLzGyh77XxQFsA59zzwHRgKJAJ7AauCXypIiJSmSqD3Dn3LZUfA6+4jwNuCVRRIiLiv5Bb2SkiIgdSkPvJ39MkI4W+j1/Td3IgfR+/VlPfiYLcT/qP8kD6Pn5N38mB9H38moJcREQqZeXzlLX4gWY5wJpa/dDASARyvS4iiOj7+DV9JwfS9/FrR/OdtHPOJVW2odaDXEREAkuHVkREQpyCXEQkxCnIRURCnIJcRCTEKcgPw582d5HKzKLNbIGZfeh1LV4zs8Zm9o6ZLfX9t3Ki1zV5zczu8P3O/GRmr5tZxHWDMbMXzSzbzH6q8FpTM/vczFb47psE4rMU5Ie3r81dV+AE4BYz6+ZxTcFiHOXdogSeBD5xznUBehHh34uZtQJuA9Kccz2AaOByb6vyxMvAkINeuxuY4ZzrCMzwPT9qCvLD8LPNXcQxs9bAecAkr2vxmpk1AgZRfqlnnHN7nXM7vK0qKMQA9cwsBqhPBPYncM59A2w76OULgCm+x1OA3wTisxTkfjpMm7tI9ARwF1DmdSFBoD2QA7zkO9Q0yczivS7KS865DcCjwFpgE+X9CT7ztqqgccy+Xg2+++aBeFMFuR+qaHMXUcxsGJDtnJvndS1BIgboC/zTOdcHKCBAfy6HKt9x3wuAY4GWQLyZjfS2qvCmIK+CH23uIs1AYLiZrQbeAAab2aveluSp9cB659y+v9TeoTzYI9mZwCrnXI5zrhh4FzjJ45qCxRYzSwbw3WcH4k0V5IfhZ5u7iOKcu8c519o5l0L5BNZM51zEjracc5uBdWbW2ffSGcAvHpYUDNYCJ5hZfd/v0BlE+ARwBdOA0b7Ho4H/BOJN/Wn1FskqbXPnnJvuYU0SfG4FXjOzOkAWEd7q0Dk3x8zeAeZTfubXAiDirmlrZq8DpwGJZrYeuA94GHjLzK6j/H94lwTks3TRLBGR0KZDKyIiIU5BLiIS4hTkIiIhTkEuIhLiFOQiIiFOQS4iEuIU5CIiIe7/A+NIxkINbcypAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, exacts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x21178212e88>]"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD2CAYAAADcUJy6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxcdZnv8c/Ta/a1Owl0QtbuLGwJNGHNwqZBlCCC4FxHnGGG0QsjA4rCyBUHZVBQ0LnDKLmKOuMSWdTJaDAi3RB20oGEkNCddDoJ6SSkqrOvvdVz/6jqWGm6k+p0dZ2q6u/79apX1znnd+o8VUk/feo5v/P7mbsjIiLZJyfoAEREpGcowYuIZCkleBGRLKUELyKSpZTgRUSyVF7QAbQpKirycePGBR2GiEhGWb58eYO7F3e0LW0S/Lhx46iqqgo6DBGRjGJmmzrbphKNiEiWUoIXEclSSvAiIllKCV5EJEsllODNbJ6Z1ZhZrZnddYx215qZm1l53Lq7Y/vVmNmHkxG0iIgc33F70ZhZLvAocDlQDywzs0XuvqZdu4HAF4DX49ZNA24ATgVOBv5sZmXu3pq8tyAiIh1J5Ax+JlDr7nXu3gQsBOZ30O4bwIPA4bh184GF7t7o7huA2tjriYhID0skwZcAm+OW62PrjjCzGcAYd/99V/dtEw6HKS8vP/JYsGBBAqGJiGS2p5fX88SyzcdveAISudHJOlh3ZBB5M8sBHgE+29V94xUXF+tGJxHpdR5bup6iAYV88pwxSX/tRM7g64H4I48GtsYtDwROA543s43AecCi2IXW4+0rItJr1e86yNrt+7lkyogeef1EEvwyoNTMxptZAdGLpovaNrr7Hncvcvdx7j4OeA24yt2rYu1uMLNCMxsPlAJvJP1diIhkoMqaMAAX91CCP26Jxt1bzOxWYAmQCzzu7qvN7D6gyt0XHWPf1Wb2BLAGaAFuUQ8aEZGoyuoQY4f3Y0JR/x55/YQGG3P3xcDiduu+1knbue2W7wfuP8H4RESy0uHmVl5Z38AN55yCWUeXK7tPd7KKiATg1bodHG6O9Fh5BpTgRUQCUVkdom9+LueOH9Zjx1CCFxFJMXenojrEhZOK6JOf22PHUYIXEUmx2tB+6ncd4uIpHU7ElDRK8CIiKVZZEwLg4sk9V38HJXgRkZSrqA4xZdRATh7St0ePowQvIpJCew83U7VxV4/dvRpPCV5EJIVeXNtAS8R7tHtkGyV4EZEUqqwJMbhvPjPGDOnxYynBi4ikSCTiPF8TYk5ZMXm5PZ9+leBFRFJk1ZY9NOxvSkn9HZTgRURSpqI6hBnMLuvZ/u9tlOBFRFKksibEjDFDGNa/ICXHU4IXEUmB8L5G3q7fk7LyDCjBi4ikxPNtd6+mW4I3s3lmVmNmtWZ2VwfbP2dmq8xshZm9ZGbTYuvHmdmh2PoVZvbDZL8BEZFMUFkTYuSgQqadNChlxzzuhB9mlgs8ClxOdI7VZWa2yN3XxDX7pbv/MNb+KuBhYF5s23p3n57csEVEMkdza4QX1zZw5Rkn9djkHh1J5Ax+JlDr7nXu3gQsBObHN3D3vXGL/QFPXogiIpmtauMu9jW2pLQ8A4kl+BJgc9xyfWzdUczsFjNbDzwIfCFu03gze8vMXjCzWd2KVkQkA1XWhMjPNS6aVJTS4yaS4Dv6PvGBM3R3f9TdJwJfAe6Jrd4GnOLuM4A7gF+aWYcFqHA4THl5+ZHHggULEnsHIiJprqI6xLnjh9O/MKFpsJMmkaPVA2PilkcDW4/RfiHwAwB3bwQaY8+Xx87wy4Cq9jsVFxdTVfWB1SIiGW3zzoPUhvbzqZmnpPzYiZzBLwNKzWy8mRUANwCL4huYWWnc4pXAutj64thFWsxsAlAK1CUjcBGRTNA2uUcq+7+3Oe4ZvLu3mNmtwBIgF3jc3Veb2X1AlbsvAm41s8uAZmAXcGNs99nAfWbWArQCn3P3nT3xRkRE0lFFdYjxRf0ZX9Q/5cdOqCDk7ouBxe3WfS3u+W2d7Pc08HR3AhQRyVSHmlp5df0O/urc1JdnQHeyioj0mFfWN9DYEgmkPANK8CIiPaayJkS/glxmjh8WyPGV4EVEeoC7U1kd5qJJRRTm5QYSgxK8iEgPWLt9P1t2H0r53avxlOBFRHpARXVs9MjJSvAiIlmlsibEtJMGMWpwn8BiUIIXEUmyPQebWb5pV2C9Z9oowYuIJNnSdWFaI87FU1Iz92pnlOBFRJKssjrE0H75TB8zNNA4lOBFRJIoEnGeXxtmTlkxuTmpm9yjI0rwIiJJtLJ+NzsPNAXaPbKNEryISBJVVofIMZhTFmz9HZTgRUSSqqImxFmnDGVIv4KgQ1GCFxFJltDew7yzZW9alGdACV5EJGmerwkDwUzu0REleBGRJKmoDnHS4D5MGTUw6FCABBO8mc0zsxozqzWzuzrY/jkzW2VmK8zsJTObFrft7th+NWb24WQGLyKSLppaIrxU28DcySMwC7Z7ZJvjJvjYnKqPAlcA04BPxSfwmF+6++nuPh14EHg4tu80onO4ngrMA/6jbY5WEZFsUrVxJ/sbW9KmPAOJncHPBGrdvc7dm4CFwPz4Bu6+N26xP+Cx5/OBhe7e6O4bgNrY64mIZJWK6hAFeTlcOGl40KEckcicrCXA5rjleuDc9o3M7BbgDqAAuCRu39fa7VtyQpGKiKSxipoQ500YTr+ChKa6TolEzuA7Kib5B1a4P+ruE4GvAPd0ZV+AcDhMeXn5kceCBQsSCE1EJHibdhygLnyAiycHf3NTvET+1NQDY+KWRwNbj9F+IfCDru5bXFxMVVVVAuGIiKSXtsk90qn+DomdwS8DSs1svJkVEL1ouii+gZmVxi1eCayLPV8E3GBmhWY2HigF3uh+2CIi6aOyJsyE4v6MHd4/6FCOctwzeHdvMbNbgSVALvC4u682s/uAKndfBNxqZpcBzcAu4MbYvqvN7AlgDdAC3OLurT30XkREUu5gUwuv1e3gM+eNDTqUD0joaoC7LwYWt1v3tbjntx1j3/uB+080QBGRdPZy7Q6aWiJpMzxBPN3JKiLSDRXVIQYU5nHOuGFBh/IBSvAiIifI3Xm+JsRFk4ooyEu/dJp+EYmIZIjq9/exbc/htOs900YJXkTkBLV1j5ybZv3f2yjBi4icoMrqEKeVDGLEoD5Bh9IhJXgRkROw+2ATb763i0smp2d5BpTgRUROyAtrw0SctOwe2UYJXkTkBFRWhxjWv4AzRg8JOpROKcGLiHRRa8R5YW2YuWXF5Oakx+QeHVGCFxHpohWbd7PrYHNal2dACV5EpMsqq0Pk5hizS9Oze2QbJXgRkS6qqA5x9ilDGdwvP+hQjkkJXkSkC97fc5g12/amfXkGlOBFRLrk+Zr0nNyjI0rwIiJdUFEdomRIX8pGDgg6lONSghcRSVBjSysv1TYwd3IxZunbPbJNQgnezOaZWY2Z1ZrZXR1sv8PM1pjZ22b2nJmNjdvWamYrYo9F7fcVEckUb2zYycGm1owoz0ACMzqZWS7wKHA50Um0l5nZIndfE9fsLaDc3Q+a2eeBB4HrY9sOufv0JMctIpJyldVhCvNyuGBiUdChJCSRM/iZQK2717l7E7AQmB/fwN0r3f1gbPE1YHRywxQRCV5lTYjzJw6nb0Fu0KEkJJEEXwJsjluuj63rzE3AM3HLfcysysxeM7OrO9spHA5TXl5+5LFgwYIEQhMRSY0NDQfY0HCAi9N49Mj2Epl0u6MrCd5hQ7NPA+XAnLjVp7j7VjObAFSY2Sp3X99+3+LiYqqqqhKJWUQk5dom98iU+jskdgZfD4yJWx4NbG3fyMwuA74KXOXujW3r3X1r7Gcd8DwwoxvxiogE4vmaEJNGDGDMsH5Bh5KwRBL8MqDUzMabWQFwA3BUbxgzmwE8RjS5h+LWDzWzwtjzIuBCIP7irIhI2jvQ2MLrdTsz6uwdEijRuHuLmd0KLAFygcfdfbWZ3QdUufsi4CFgAPBkrG/oe+5+FTAVeMzMIkT/mHyrXe8bEZG091JtA02tkbSde7UzidTgcffFwOJ2674W9/yyTvZ7BTi9OwGKiAStsjrEwMI8zhk3LOhQukR3soqIHIO7U1kTYlZZEfm5mZUyMytaEZEUW7NtL9v3NmZU98g2SvAiIsdQGeseOSfD6u+gBC8ickwV1SHOGD2YEQP7BB1KlynBi4h0YueBJt7avDsjyzOgBC8i0qmla8O4Z9bdq/GU4EVEOlFRHaJoQAGnlwwOOpQTogQvItKBltYIL6wNM6dsBDk56T+5R0eU4EVEOvDW5t3sOdScseUZUIIXEelQZXWI3BxjVllmTO7RESV4EZEOVFSHKB87lEF98oMO5YQpwYuItLN19yGq39+X0eUZUIIXEfmAyprMm9yjI0rwIiLtVFaHGT20L5NGDAg6lG5RghcRiXO4uZWXaxu4ePIIYvNbZCwleBGROK9v2Mmh5taML89AggnezOaZWY2Z1ZrZXR1sv8PM1pjZ22b2nJmNjdt2o5mtiz1uTGbwIiLJVlkdok9+DudPHB50KN123ARvZrnAo8AVwDTgU2Y2rV2zt4Bydz8DeAp4MLbvMOBe4FxgJnCvmQ1NXvgiIsnj7lRUh7hgYhF98nODDqfbEjmDnwnUunuduzcBC4H58Q3cvdLdD8YWXwNGx55/GHjW3Xe6+y7gWWBeckIXEUmuuoYDvLfzIBdn4NjvHUkkwZcAm+OW62PrOnMT8ExX9w2Hw5SXlx95LFiwIIHQRESSp21yj4uzoP4OiU263dFlZO+wodmngXJgTlf3LS4upqqqKoFwRER6RkV1iLKRAxg9tF/QoSRFImfw9cCYuOXRwNb2jczsMuCrwFXu3tiVfUVEgrbvcDPLNu7MmrN3SCzBLwNKzWy8mRUANwCL4huY2QzgMaLJPRS3aQnwITMbGru4+qHYOhGRtPJybQPNrZ6xszd15LglGndvMbNbiSbmXOBxd19tZvcBVe6+CHgIGAA8Gbsx4D13v8rdd5rZN4j+kQC4z9139sg7ERHphorqEAP75HH22Ozp6JdIDR53Xwwsbrfua3HPLzvGvo8Dj59ogCIiPS0ScSprwswuKyY/N3vu/8yedyIicoLWbNtLeF8jl2RReQaU4EVEqKgOYQZzsqT/exsleBHp9SqqQ5wxeghFAwqDDiWplOBFpFfbsb+RlfW7s648A0rwItLLvbA2jHvmT+7RESV4EenVKqpDFA0o5NSTBwUdStIpwYtIr9XSGmHp2jAXTy4mJyezJ/foiBK8iPRayzftYu/hlqwsz4ASvIj0YpU1YfJzjYtKi4IOpUcowYtIr1VZHeKcccMY2Cc/6FB6hBK8iPRKW3Yfomb7vqwaXKw9JXgR6ZUqsmxyj44owYtIr/R8dYhThvVjYnH/oEPpMUrwItLrHG5u5eX1DVwyZQSxIc6zkhK8iPQ6r9bt4HBzhLlZNrhYe0rwItLrVFaH6Jufy3kThgcdSo9KKMGb2TwzqzGzWjO7q4Pts83sTTNrMbNr221rNbMVscei9vuKiKSSu1NRHeLCScPpk58bdDg96rgzOplZLvAocDnRSbSXmdkid18T1+w94LPAlzp4iUPuPj0JsYqIdNv68H7qdx3i83MnBh1Kj0tkyr6ZQK271wGY2UJgPnAkwbv7xti2SA/EKCKSNG3dI+dmcf/3NomUaEqAzXHL9bF1iepjZlVm9pqZXd1Zo3A4THl5+ZHHggULunAIEZHEVFSHmDJqICVD+gYdSo9L5Ay+oz5E3oVjnOLuW81sAlBhZqvcfX37RsXFxVRVVXXhZUVEumbv4WaqNu7i72dPCDqUlEjkDL4eGBO3PBrYmugB3H1r7Gcd8DwwowvxiYgkzUvrGmiJeNaOHtleIgl+GVBqZuPNrAC4AUioN4yZDTWzwtjzIuBC4mr3IiKpVFEdYnDffGaMGRJ0KClx3ATv7i3ArcAS4F3gCXdfbWb3mdlVAGZ2jpnVA9cBj5nZ6tjuU4EqM1sJVALfatf7RkRO0K4DTTS2tAYdRsaIRJzna0LMLismL7d33AKUSA0ed18MLG637mtxz5cRLd203+8V4PRuxigi7exvbGHe95dSMqQvT37uAnKzcDaiZFu1ZQ8N+5u4ZEp2370ar3f8GRPJMj94vpbtext5873d/PiluqDDyQiVNSHMYE5Z76i/gxK8SMbZsvsQP3pxA1dPP5nLp43kO39aS21of9BhpbUDjS08tbyeGWOGMKx/QdDhpIwSvEiGeeiP1QDcOW8K93/8NPrm53LnUytpjXSl93Lv8s0/vMuW3Ye464qpQYeSUkrwIhlk5ebd/G7FVv5u1nhKhvRlxMA+/MtVp/KWSjWdqqjezq/eeI+bZ09g5vhhQYeTUkrwIhnC3fnmH9ZQNKCAz8+ddGT9fJVqOrXzQBNffmoVU0YN5I7Ly4IOJ+WU4EUyxJLV77Ns4y7uuHwyAwr/0gHOzLj/6mip5ssq1Rzh7nz1t6vYc6iJR66fTmFedo8c2REleJEM0NQS4YFnqpk8ciCfLP9Aj2RGDIqWatSr5i9++9YWnnnnfe64fDJTTxoUdDiBUIIXyQD/+epGNu04yD9fObXTm3RUqvmLrbsPce9/r+accUO5uZeMO9MRJXiRNLfrQBP/9tw65pQVM6es85t0VKqJikScLz25kog7371ueq++CUwJXiTN/VvFOvY3tvDVK4/fxS++VPP4SxtSEF36+ekrG3ll/Q7+z0enccrwfkGHEygleJE0Vhfez3+9uokbZp5C2ciBCe0zf/rJXDZ1JN/5Uw3rw72rVFMb2se3/1jNpVNGcP05Y46/Q5ZTghdJY996pprCvBxuvyzxLn5mxr9+/DT65Ody55O9p1TT3Brh9l+vpF9BLg984nTMem9ppo0SvEiaeq1uB39as53/ffEkigcWdmnf3liq+b8VtazasocHrjmdEQP7BB1OWlCCF0lDkUj0pqaSIX256aLxJ/QavalUs2Lzbh6trOWas0qYd9pJQYeTNpTgRdLQb9/awjtb9vLleZPpk39iN+j0llLNoaZW7vj1CkYN6sPXrzo16HDSihK8SJo51NTKQ0tqOHP0YD52xsndeq0Rg/rw9aumZXWp5oFn3qWu4QAPXXcGg/rkBx1OWkkowZvZPDOrMbNaM7urg+2zzexNM2sxs2vbbbvRzNbFHjcmK3CRbPX/Xqzj/b2Hueej08hJQh/uq6eXZG2pZunaMP/56ib+9sLxXDCxKOhw0s5xE7yZ5QKPAlcA04BPmdm0ds3eAz4L/LLdvsOAe4FzgZnAvWY2tPthi2Sn0N7D/PCF9Vxx2ijOGZeckQ+ztVSz+2ATdz61kkkjBvDleZODDictJXIGPxOodfc6d28CFgLz4xu4+0Z3fxuItNv3w8Cz7r7T3XcBzwLzkhC3SFb67p/W0twa4a4rpiT1dbOxVPO1/17Njv1NfO/66Sd8nSLbJZLgS4DNccv1sXWJSHjfcDhMeXn5kceCBQsSPIRIdlizdS9PLN/MjeePY+zw/kl//Wwq1fzPyq0sWrmV2y4t5bSSwUGHk7YSmXS7oyJgot/xEt63uLiYqqqqBF9WJLu4O/cvXsPgvvn84yWlPXKMtlLN5Y8s5c4nV2bsZN3v7znMPb97h+ljhvD5uRODDietJXIGXw/E3/M7Gtia4Ot3Z1+RXqOyJsTLtTu47dJSBvfruZ4g8aWan7yceaUad+fLT79NY0srD3/yzE5H1pSoRD6dZUCpmY03swLgBmBRgq+/BPiQmQ2NXVz9UGydiMQ0t0a4/w/vMqGoP58+b2yPHy9aqhnBQ0syr1Tz89c2sXRtmK9+ZCoTigcEHU7aO26Cd/cW4Faiifld4Al3X21m95nZVQBmdo6Z1QPXAY+Z2erYvjuBbxD9I7EMuC+2TkRiFr7xHuvDB7jriinkp+CMNFqqOT3jetXUhfdz/+J3mV1WnJI/hNnA3NPjH7e8vNxVg5feZu/hZuY+9DxlIwfwq78/L6UDZP32rXpu//VK7rlyKn83K70nxWhpjXDtD19lQ8MB/nT7bEYO0lgzbcxsubuXd7RNBSyRAD1aWcuug03cc+W0lI9+mEmlmh88v54Vm3fzzatPU3LvAiV4kYBs3nmQn7y0kWtmjA6kq19bqaYwL4cvP/V22pZqVtXv4fvPreOqM0/mY2d2b+iG3kYJXiQg3/5jNTk5cOeHg7sLc0RsgK7lm3alZa+aw82t3P7ECoYPKOC++RpIrKuU4EUCsHzTLn7/9jZunj2RUYODLTl8fEb6lmoe/GMNtaH9PHTtmQzpVxB0OBlHCV4kxdyjY70XDyzkH2YHf3HTzLg/DUs1r9Q28PjLG/jM+WOZfYzJxqVzSvAiKfaHVdt4673d3PmhyfQvTORm8p43Ms1KNXsPN/OlJ1cyoag/d19x/MnGpWNK8CIpdLi5lW89U82UUQP5xNmjgw7nKPGlmrqASzVfX7Sa7fsaefj66fQt0EBiJ0oJXiSFfvbKRup3HeKeK6el3Tgw8aWaOwMs1fzxnW385s0t3DJ3ItPHDAkkhmyhBC+SIjv2N/LvFbVcMmUEF5Wm5+QUQZdqQvsOc/dvVnF6yWD+8dKeGXStN1GCF0mR7z+3joPNrfzzR5I71nuyfXxGCZdOSX2pxt256+lVHGxq5ZHrz0zJsA3ZTp+gSArUhvbxi9ff469mnsKkEQODDueYzIx/vSb1pZpfL9tMRXWIr8ybkvafUaZQghdJgQcWV9MvP5d/uiwzyg6pLtW8t+Mg3/j9Gi6YOJzPXjCux4/XWyjBi/Swl2sbeK46xC2XTGL4gMKgw0lYqko1rRHnjidWkGPGQ9edmZSJxiVKCV6kB7VGnG/+4V1GD+2bcWem8aWanrwBasHSOqo27eJf5p9KyZC+PXKM3koJXqQHPb28nne37eUr86Zk5MTQIwf14d6PnUpVD5Vq1mzdy8PP1nDFaaP4+IxEp3qWRCnBi/SQA40tfOdPNcw4ZQgfPeOkoMM5Ydec1TOlmsaWVu54YgWD+xZw/8dPT/lwyb1BQgnezOaZWY2Z1ZrZXR1sLzSzX8e2v25m42Lrx5nZITNbEXv8MLnhi6Svx5bWEdrXGMhY78nUU6Wah59dS/X7+3jw2tMZ1l8DifWE4yZ4M8sFHgWuAKYBnzKzae2a3QTscvdJwCPAt+O2rXf36bHH55IUt0hae3/PYRYsXc9HzziJs8cODTqcbkt2qeaNDTtZsLSOT80cwyVTRiYhQulIImfwM4Fad69z9yZgITC/XZv5wM9iz58CLrVMPmUR6aaHltQQicBX5qX3TU1dEV+q2dBw4IRfZ39jC198cgVjhvbjnivbnytKMiWS4EuAzXHL9bF1HbaJTdK9Bxge2zbezN4ysxfMbFZnBwmHw5SXlx95LFiwIOE3IZJO3tmyh9+8Vc/fXDSOMcP6BR1O0hx1A1Q3Juv+xv+sYcuuQzz8yTPTZjTNbJXIp9vRmXj7f9nO2mwDTnH3HWZ2NvA7MzvV3fe2b1xcXIwm3ZZM1zbW+9B+Bdxy8aSgw0m6tlLNF59cyU9e3tDlybqfXbOdX1dt5vNzJ1I+blgPRSltEjmDrwfGxC2PBrZ21sbM8oDBwE53b3T3HQDuvhxYD5R1N2iRdPXnd0O8VreT2y8rZVCf/KDD6RHXnFXCJSdQqtmxv5G7f/M2U08axO2XKQ2kQiIJfhlQambjzawAuAFY1K7NIuDG2PNrgQp3dzMrjl2kxcwmAKVAXXJCF0kvza0RHlj8LhOL+/OpmacEHU6PMTMe6GKpxt3559+uYu+hFh65/kwK8tRDOxWO+ynHauq3AkuAd4En3H21md1nZlfFmv0YGG5mtcAdQFtXytnA22a2kujF18+5+85kvwmRdPCL1zZR13CAr145lbwsHwkxvlfNT1/ZeNz2T7+5hSWrt/PFD5UxZdSgng9QADD39Jh/sby83FWDl0y152Azc75TyaknD+LnN52b0f3eE+Xu3PSzKl5Z38Azt81mfFH/DtvV7zrIvO+9yLSTBvGrm89Lu4lOMp2ZLXf38o62ZfdphkiK/HvlOvYcauarH8nsm5q6oq1UU5DbeakmEnG+9ORK3J3vfvJMJfcUU4IX6aZNOw7w01c2ct3Zo5l2cu8qP4wc1IevHaNU8/jLG3itbif3fuzUrOoymimU4EW66VvPVJOfm8MXPzQ56FAC8YkjvWqqj+pVs3b7Ph5cUsNlU0dyXXl6TTDeWyjBi3TDso07eead9/mH2RMZOahP0OEEoqNSTVNLhNt/vYKBhXk8cI0GEguKbiMTOUGRiPPN369h1KA+/P3s8UGHE6i2Us2XnlzJT1/ZyK4DTazeupfH/vpsigdmziQn2UYJXlLC3Yk4RNzxuJ/O0esH9cnLmLO9/3l7Kyvr9/Cd686kX4F+lT5xVgmLV23jwT9W09wa4dqzR/PhU0cFHVavpv+VGaK5NcKh5lYONbVysCn681BzS9zz1nbPWzjUFOm0TWvEo0mVuGTbLgl73PaIE9cmtj7SlqTb2rTtxwdeO1EjBhYyp6yYOZOLuWhSEUP6pecwsoebW/n2M9WcVjKIazRRBRAbq+bjp3P5Iy9QNKCQez+mgcSCpgTfQyIRp/r9few40HhCCfgvbaLPm1u7dr9CjkG/gjz6FuTSryCXvvm5R54P7ZdPXk4OZpBjhln0lzPHooMKRddZbHtHbWI/j7T54HLb6+QYEPt5ZH1O7PWwo14/4s7K+j0sWf0+Ty6vJ8dg+pghzCkbwZzJxZxeMjhtutn9+KUNbN1zmO9+crrmEI0zanAfFt16EX3zcxmYpUM1ZBIl+CTauvsQL61rYOm6MC/XNrDrYHOnbeMTcN/8WBKOPR/aL5++BXn0zc/ptE30eV4n63MpyM3JmFJHey2tEVbW7+GFtWFeWBvme8+t5ZE/r2Vov3xmlRYzp6yYWWVFjBgYzEXN8L5G/qOylsunjeT8icOPv0Mv09kNT5J6SvDdcLCphdfrdrJ0XZgX1zVQG4pOZzZiYCGXTBnJRaXDKRnSL+sScE/Ly83h7LFDOXvsUO64vIydB5p4cV002S9d28CilVQhyLoAAAeZSURBVNGx7k49eVC0nFNWzFljh5KfouEBHvnzWhpbItx9RfaM9S7ZSUMVdEEk4qzZtjea0Nc2ULVpJ82tTmFeDudOGM7s0iJmlRZTNnKAkncPafs3aDu7f3PTLloizoDCPC6cNJzZsYQ/emjP3FRT8/4+rvj+Uj5z/ji+ftWpPXIMka441lAFOoM/jvf3HObF2Bn6S7UN7DzQBMDUkwbxtxeOZ1ZpMeXjhtInPzfgSHuHnBzjtJLBnFYymFsunsTew828UrsjdnYfZsnq7QBMLO5/pHZ/7vhhSfv3+dfF7zKgMI/bLi1NyuuJ9CQl+HYONbXy+oYdvLiugRfXhVm7PVp2KRpQyNxY7ffCScHVf+Vog/rkM++0Ucw7bRTuzvrwfl5Y28ALa8P8/PVNPP7yBgrzcjhvwvAjvXMmFPU/oW9Ybd8a7rlyKkM1SbRkgF5fomn7yh89Qw+zbMMumlojFOTlcO74YcyKlV2mjBqoskuGaftj3ZaY68LR2+hHD+17pHZ/waQiBiQwbVxLa4SP/NuLHG6O8OwdsynM0zc2SQ8q0bSzfe/hI2foL61rYEes7DJl1EA+e+E4ZpUWcc645H2tl2D0Lchl7uQRzJ08AoDNOw8eSfa/e2sLv3j9PfJyjPJxQ6PlnLJipp7U8R/yJ6rqWbt9Pz/4X2cpuUvG6BVn8IeaWnlj405eXButpdds3wdA0YACZpUWM6u0iIsmFTGil44l0hs1tURYvmnXkYT/7rboNMHFbTdalUX/XwzpV8D+xhbmPlTJ+KL+PPEP5+ubnKSVY53BJ5TgzWwe8H0gF/iRu3+r3fZC4D+Bs4EdwPXuvjG27W7gJqAV+IK7L+noGMlM8G03GbVdHH1j406aWqJll5njji676CYVgei3uqWxZP/iugb2HGomx+DMMUMYUJjHi+sa+N0tFzJ9zJCgQxU5SrdKNLE5VR8FLic6ufYyM1vk7mvimt0E7HL3SWZ2A/Bt4Hozm0Z0DtdTgZOBP5tZmbu3du8tfVBo72Feqm2IlV4aaNjfCMDkkQP5zHljmVVWzMxxw+hboK/X8kEjB/XhuvIxXFc+htaIs7J+Ny/URBP+S7UNfOKs0UruknESqcHPBGrdvQ7AzBYC84H4BD8f+Hrs+VPAv1v0e+x8YKG7NwIbYnO2zgReTU74sGX3IW766TKq34+WXYb3L+Ci2Bn6rNKiXjuEq5y43BzjrFOGctYpQ7n98jL2HW6mr67HSAZKJMGXAJvjluuBcztr4+4tZrYHGB5b/1q7fTscmSkcDlNe/pdvGTfffDM333zzcYMbObCQkiF9uXpGCbNKi5g6apDKLpJUGlNFMlUiCb6jbNm+cN9Zm0T2BaC4uJgTqcHn5ebw48+e0+X9RESyXSKDd9QDY+KWRwNbO2tjZnnAYGBngvuKiEgPSCTBLwNKzWy8mRUQvWi6qF2bRcCNsefXAhUe7Z6zCLjBzArNbDxQCryRnNBFRORYjluiidXUbwWWEO0m+bi7rzaz+4Aqd18E/Bj4r9hF1J1E/wgQa/cE0QuyLcAtPdGDRkREPqhX3OgkIpKtjtUPPjUDaIuISMopwXfTggULgg4h7egzOZo+jw/SZ3K0nvo8lOC7Sf9RP0ifydH0eXyQPpOjKcGLiEiXpM1FVjMLA5uCjuMEFAENQQeRZvSZHE2fxwfpMzladz6Pse5e3NGGtEnwIiKSXCrRiIhkKSV4EZEspQQvIpKllOBFRLKUEvwJMrMxZlZpZu+a2Wozuy3omNKBmeWa2Vtm9vugY0kHZjbEzJ4ys+rY/5Xzg44pSGZ2e+z35R0z+5WZ9boZeczscTMLmdk7ceuGmdmzZrYu9nNoMo6lBH/iWoAvuvtU4DzgltgUhb3dbcC7QQeRRr4P/NHdpwBn0os/GzMrAb4AlLv7aUQHL7wh2KgC8VNgXrt1dwHPuXsp8FxsuduU4E+Qu29z9zdjz/cR/cXtcLaq3sLMRgNXAj8KOpZ0YGaDgNlER1vF3ZvcfXewUQUuD+gbmzeiH71wfgh3X0p01N1484GfxZ7/DLg6GcdSgk8CMxsHzABeDzaSwH0P+DIQCTqQNDEBCAM/iZWtfmRm/YMOKijuvgX4DvAesA3Y4+5/CjaqtDHS3bdB9OQRGJGMF1WC7yYzGwA8DfyTu+8NOp6gmNlHgZC7Lw86ljSSB5wF/MDdZwAHSNJX70wUqyvPB8YDJwP9zezTwUaV3ZTgu8HM8okm91+4+2+CjidgFwJXmdlGYCFwiZn9PNiQAlcP1Lt72ze7p4gm/N7qMmCDu4fdvRn4DXBBwDGli+1mdhJA7GcoGS+qBH+CzMyI1lbfdfeHg44naO5+t7uPdvdxRC+cVbh7rz47c/f3gc1mNjm26lKis5v1Vu8B55lZv9jvz6X04ovO7cRPe3oj8N/JeNHjTtknnboQ+GtglZmtiK37Z3dfHGBMkn7+EfhFbD7jOuBvAo4nMO7+upk9BbxJtBfaW0CvGzfYzH4FzAWKzKweuBf4FvCEmd1E9A/hdUk5lgYbExHJTirRiIhkKSV4EZEspQQvIpKllOBFRLKUEryISJZSghcRyVJK8CIiWer/A8eyMQG5nr/HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, errors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsurprisingly, we see the lowest error around 1 or 2 dimensions, with still low values at 3 and 4. This makes sense, since spheres really only have physical meaning in low dimensions. Interestingly, there is another dip at just before 8. I personally like 2 dimensions and a radius of 1, because those make the most intuitive physical sense to me, so that's what I'm going to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_nrnd = 10000\n",
    "\n",
    "errors = []\n",
    "x = []\n",
    "n = 100\n",
    "\n",
    "while n < max_nrnd:\n",
    "    estimate = dsphere_volume(d=2, nrnd=n)\n",
    "    exact = dsphere_volume_exact(dim=2)\n",
    "    error = np.abs(1-(estimate/exact))\n",
    "    x.append(n)\n",
    "    errors.append(error)\n",
    "    n += 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x211785d9588>]"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD2CAYAAADbPoDqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAc90lEQVR4nO3df6xk5X3f8feHWd8LdW0wy3VEdlEWBI20BMk2t8BtomqVrfGCHK+rEHkxCrSlumrtTetaacVVupULDWuqKiYOyAYMDqASoBCSlRN3myxeVQmTLbsCGxa89WWNzWXdcvlhStyam9399o957jI73LlzZubMnDNzPi9pdGfOPHPuec6v73me55znUURgZmbVdUrRC2BmZsVyIDAzqzgHAjOzinMgMDOrOAcCM7OKW1P0AmRx1llnxYYNG4peDDOzkXLgwIFXI2KqU7qRCAQbNmxg//79RS+GmdlIkfSDLOlcNWRmVnEOBGZmFedAYGZWcQ4EZmYV50BgZlZxDgRmZhVXyUBQr9fZuXMn9Xq96EUxMyvcSDxHkKd6vc7mzZtZWlpiYmKCPXv2MDMzU/RimZkVpnIlgr1797K0tMSxY8dYWlpi7969RS+SmVmhKhcINm3axMTEBLVajYmJCTZt2lT0IpmZFapyVUMzMzPs2bOHvXv3smnTJlcLmVnlVS4QQCMYOACYmTVUrmrIzMxOlikQSNoi6ZCkeUk3rPD9pKSH0vf7JG1I0y+R9HR6fVvSP8w6TzMzG46OgUBSDbgduALYCFwtaWNLsuuBNyLifOBLwC1p+rPAdER8CNgC3CFpTcZ5mpnZEGQpEVwCzEfE4YhYAh4Etrak2Qrcm94/AmyWpIj4vxFxNE0/FYgu5mlmZkOQJRCsA15q+ryQpq2YJp343wTWAki6VNJB4Bngn6Xvs8zTzMyGIEsg0ArTImuaiNgXERcCfxeYk3RqxnmesLi4yPT09InXnXfemWGxzcwsiyy3jy4A5zR9Xg8caZNmQdIa4HTg9eYEEfG8pJ8Av5BxnidMTU2N5VCV9XrdzzOYWeGyBIIngQsknQu8DGwDPt2SZhdwHVAHrgIej4hIv3kpIo5K+jng54EXgR9nmOdYc59HZlYWHauGUp3+dmA38DzwcEQclHSjpE+kZHcDayXNA58Hlm8H/SXg25KeBh4DPhMRr7abZ54ZKzv3eWRmZZHpyeKI+FPgT1um/bum9z8Ffm2F390P3J91nlWy3OfRconAfR6ZWVEq2cVEGbjPIzMrCweCArnPIzMrA/c1ZGZWcQ4EZmYV50BgZlZxDgRmZhXnQGBmVnEOBGZmFedAYGZWcQ4EZmYV50BgZlZxDgRmZhXnQGBmVnEOBGZmFedA0Kd6vc7OnTup1+tFL4qZWU/c+2gfPMqYmY0Dlwj64FHGzGwcOBD0YXmUsVqt5lHGzGxkuWqoD1lHGavX6x6JzMxKy4GgT51GGXM7gpmVnauGBsztCGZWdg4EA+Z2BDMrO1cNDVjWdgQzs6JkCgSStgC/C9SAr0XEF1u+nwTuAy4GXgM+FREvSvoo8EVgAlgC/nVEPJ5+sxc4G/h/aTaXR8QrfeeohDq1I5iZFaljIJBUA24HPgosAE9K2hURzzUlux54IyLOl7QNuAX4FPAq8CsRcUTSLwC7gXVNv7smIvbnlBczM+tBljaCS4D5iDgcEUvAg8DWljRbgXvT+0eAzZIUEU9FxJE0/SBwaio9mJlZSWQJBOuAl5o+L3DyVf1JaSLiKPAmsLYlza8CT0XE203Tvi7paUk7JKndAiwuLjI9PX3ideedd2ZYbDMzyyJLG8FKJ+joJo2kC2lUF13e9P01EfGypPcBjwK/TqOd4V2mpqbYv981SGZmg5ClRLAAnNP0eT1wpF0aSWuA04HX0+f1wGPAtRHxwvIPIuLl9Pct4AEaVVBmZjZkWQLBk8AFks6VNAFsA3a1pNkFXJfeXwU8HhEh6QzgT4C5iPjL5cSS1kg6K71/D/Bx4Nn+smJmZr3oGAhSnf92Gnf8PA88HBEHJd0o6RMp2d3AWknzwOeBG9L07cD5wI7UFvC0pA8Ck8BuSd8BngZeBu7KM2NmZpaNIlqr+8tneno63EZgZtYdSQciYrpTOncxYWZWcQ4EZmYV50BQUh4L2cyGxZ3OlZDHMDCzYXKJoIQ8hoGZDZMDQQl5DAMzGyZXDZWQxzAws2FyICgpj2FgZsPiqiEzs4pzIDAzqzgHAjOzinMgMDOrOAcCM7OKcyAwM6s4BwIzs4pzIDAzqzgHAjOzinMgMDOrOAcCM7OKcyAwM6s4BwIzs4pzIDAzqzgHAjOzissUCCRtkXRI0rykG1b4flLSQ+n7fZI2pOkflXRA0jPp7y83/ebiNH1e0pclKa9MmZlZdh0DgaQacDtwBbARuFrSxpZk1wNvRMT5wJeAW9L0V4FfiYiLgOuA+5t+8xVgFrggvbb0kQ+zSqnX6+zcuZN6vV70otgYyDJC2SXAfEQcBpD0ILAVeK4pzVbgC+n9I8BtkhQRTzWlOQicKmkSOBN4f0TU0zzvAz4JfLOPvJRKvV73UJM2EPV6nc2bN7O0tMTExAR79uzxPmZ9yRII1gEvNX1eAC5tlyYijkp6E1hLo0Sw7FeBpyLibUnr0nya57mu3QIsLi4yPT194vPs7Cyzs7MdF7yok7EPVBukvXv3srS0xLFjx1haWmLv3r3ev6wvWQLBSnX30U0aSRfSqC66vIt5njA1NcX+/fs7LObJijwZ+0C1Qdq0aRMTExMn9u1NmzYVvUg24rI0Fi8A5zR9Xg8caZdG0hrgdOD19Hk98BhwbUS80JR+fYd59mWlk/GwLB+otVrNB6rlbmZmhj179nDTTTe5tGm5yFIieBK4QNK5wMvANuDTLWl20WgMrgNXAY9HREg6A/gTYC4i/nI5cUT8SNJbki4D9gHXAr/Xd26aFHnVtHyguo3ABmVmZsb7leVGEW1rZN5JJF0J3ArUgHsi4rcl3Qjsj4hdkk6lcUfQh2mUBLZFxGFJ/xaYA77XNLvLI+IVSdPA7wOn0Wgk/o1oszDT09PRbdUQuMHWzKpN0oGImO6YLksgKFqvgcDMrMqyBgI/WWxmVnEOBGZmFedAYGZWcQ4EZmYV50BgZlZxDgRmZhXnQGBmVnEOBGZmFedAYGZWcQ4EZmYV50BgZlZxDgQ2MB5O0Ww0ZOmG2qxrHqXNbHS4RGADUeTAQGbWHQcCGwiP0mY2Olw1NMLKPPCOR2kzGx0OBCNqFOrgPZyi2Whw1VCTUbrLxXXwZpYXlwiSUbjCbrZcB7+8vK6DN7NeORAkK11hlzkQuA7ezPLiQJCM4hW26+DNxktRN4A4ECS+wjazIhVZPZ2psVjSFkmHJM1LumGF7yclPZS+3ydpQ5q+VtK3JP21pNtafrM3zfPp9PpgHhnqx8zMDHNzcw4CZjZ0Rd4A0jEQSKoBtwNXABuBqyVtbEl2PfBGRJwPfAm4JU3/KbAD+M02s78mIj6UXq/0kgEzs3FQ5EOYWaqGLgHmI+IwgKQHga3Ac01ptgJfSO8fAW6TpIj4CfAXks7Pb5HNzMZPkdXTWQLBOuClps8LwKXt0kTEUUlvAmuBVzvM++uSjgGPAv8hImKlRIuLi0xPT5/4PDs7y+zsbIZFNzMbHUXdAJIlEGiFaa0n7CxpWl0TES9Leh+NQPDrwH0rJZyammL//v0dF9TMzLqXpbF4ATin6fN64Ei7NJLWAKcDr68204h4Of19C3iARhWUmZkNWZZA8CRwgaRzJU0A24BdLWl2Adel91cBj7er5oFGsJB0Vnr/HuDjwLPdLryZmfWvY9VQqvPfDuwGasA9EXFQ0o3A/ojYBdwN3C9pnkZJYNvy7yW9CLwfmJD0SeBy4AfA7hQEasCfA3flmjMzM8tEq1y4l8b09HS4jcDMrDuSDkTEdKd07n3UzKziHAjMzCrOgSBHZR/PoOzLZ2bFqHync3n19lf28QzKvnxmVpxKlwiWT447duxg8+bNfV0pl33EsLIvn5kVp9KBIM+TY5EdRmVR9uUzs+JUumooz8Foyj6eQdmXz8yKU/nnCIoaEagbo7CMZlY+WZ8jqHSJAMo/3KMbec1s0CrdRjAK3MhrZoPmQFBybuQ1s0GrfNVQ2bmR18wGzYFgBJS9HcPMRpurhqy03CWG2XC4RGCl5LulzIbHJQIrJd8tZTY8DgRWSr5bymx4XDVkpeS7pbLzk+fWLwcCKy3fLdWZ21IsD64aMhthbkuxPDgQmI0wt6VYHlw1ZDbC3JZieXAgMBtx7dpS3IhsWWWqGpK0RdIhSfOSbljh+0lJD6Xv90nakKavlfQtSX8t6baW31ws6Zn0my9LUh4ZyoufarVRlucwrDb+OgYCSTXgduAKYCNwtaSNLcmuB96IiPOBLwG3pOk/BXYAv7nCrL8CzAIXpNeWXjIwCD6IbNS5Edm6kaVEcAkwHxGHI2IJeBDY2pJmK3Bvev8IsFmSIuInEfEXNALCCZLOBt4fEfVoDJF2H/DJfjKSJx9EvXNJqhzciGzdyNJGsA54qenzAnBpuzQRcVTSm8Ba4NVV5rnQMs917RZgcXGR6el3RlubnZ1ldnY2w6L3Js+xjKskj3vaXa+dDzciWzeyBIKV6u5bBzrOkqbn9FNTUwxqzOKV+CDqzUolqW7WnR+OypcfyLOssgSCBeCcps/rgSNt0ixIWgOcDrzeYZ7rO8yzUD6IutdvSarfQGJmvckSCJ4ELpB0LvAysA34dEuaXcB1QB24Cng81f2vKCJ+JOktSZcB+4Brgd/rYfmtRPotSblKzqwYWuV8/U4i6UrgVqAG3BMRvy3pRmB/ROySdCpwP/BhGiWBbRFxOP32ReD9wATwY+DyiHhO0jTw+8BpwDeB32gXPKanp2OYVUNWHLcRmOVH0oGImO6YLksgKJoDgZlZ97IGAvc1ZGZjx7cxd8ddTJgNkau+Bs93n3XPgcCsC/2cyH2CGo5e7j6reoB2IDDLqN8TuW+PHY5u7z5zgHYbgY24LHXBedUX99v1iLt9GI7l25hvuummTCd1dynjEoGNsCxXcq1pbr31Vl577bVCnnPwE+vD080DoX5+xYHAhmQQdbCtV3L33Xffu/5Hc5q3336b7du3c/z48Z6qAHo5kbfm20+sl48DNBARpX9dfPHFYaPriSeeiNNOOy1qtVqcdtpp8cQTT+Q+34mJiZicnHzX/2hOs2bNmjjllFMCiFqtFjfffHMuy5Fl+Trl+4knnoibb7657zRmzWg89NvxHOsSgb3rqjXvq/dBNZI2X8n98Ic/5K677nrX/2hOs3btWj73uc8NrQoga757qeKqYoOmDY4DQcWtVIfefLLM44QzyDrY5ZN9vV7n3nvvXfF/NFfHXHTRRUOrAsia7ywBw3cc2SA5EFRc6wnm0Ucfzf2EM4w62Kz/Y5h19FmXKUvAcIOmDZL7Gqq4YZQIrLMs1XFVf+jJuudO5yyzQbcRmFkxHAjMuuQAaOMmayBwG4EZvivHqs1dTJjhbgasOGXoMtslgjHhao3++K4cK0JZSqIOBGOgLDvTKCtTNwMO6tVRludDKhMI8jy4ynaglmVnGnVl6AfIQb1aylISrUQgyPPgKuOButrOVLagZatzUK+WspREKxEI8jy4ynigttuZyhi0bHVFXyH6wiGbPNdTGUqilQgEeR5cRR+o7ay0M5UxaNnqirxC9IVDNuO4njIFAklbgN8FasDXIuKLLd9PAvcBFwOvAZ+KiBfTd3PA9cAx4F9ExO40/UXgrTT9aJaHHnqV58FVdFGumyuRsgYtW11RV4i+cMhmLNdTp36qaZz8XwDOAyaAbwMbW9J8Bvhqer8NeCi935jSTwLnpvnU0ncvAmdl6Svb4xE09NKvv/uwH23D3H6DGjei12UZdL57/R9lWk+dkHE8giyBYAbY3fR5DphrSbMbmEnv1wCvAmpN25LOgaBLN998c9RqtaENrFImVQxoRZxwyrCeV8t3XsvX70VVGdZTFlkDQZaqoXXAS02fF4BL26WJiKOS3gTWpul/1fLbdcuFEeC/SQrgjoi4s90CLC4uMj39Ts3R7Owss7OzGRZ9vFS1qmcc62SzKKIKIku1VC8Npd38pl2+89wPul23K/3vubm5nv53GWUJBFphWmtPde3SrPbbX4yII5I+CPyZpO9GxH9faQGmpqZwp3PFt0/kpdsTyVjWyWZQxsDfy8m429+0y3ee+0G363bc98EsgWABOKfp83rgSJs0C5LWAKcDr6/224hY/vuKpMeAS4AVA4G9owy3mvWjlxNJGU+Iw1DGwN/LCbHb37TLd577Qbfrttf/3XzRA5RqW56kU90RjWBxmEZj73Jj8YUtaT7LyY3FD6f3F3JyY/FhGo3P7wXel9K8F3gC2NJuGdxGMD56becYlTrZcddr3XpebR1F7gfd/u/mfE9MTMTk5OTQG5jJq40gGnX+22k09NaAeyLioKQb0z/ZBdwN3C9pnkZJYFv67UFJDwPPAUeBz0bEMUk/AzwmaTnQPBAR/7XHWDYUftAmH71eWY16SWhcdHMl3XzM5Hn7dlH7Qbf/u7kkdPz4caBx4V3GqiUPTJNB3l1UVD2geB2Mv6o28DdrXge1Wg1JHD16dKjrwwPT5CivhiIfHA2+uh9/4964mkVr6QnK20bgQJBBXo1U43Rw+KreVlPVBv5WrRc9ZT1WHAgyyOvujXE5OFyyqZ52gb/d9DLe8WTtORBklEd1xrgcHHmXbLotXbg0MlztAn+nC4JRrQKs4v7lQDBko3pwNMuzZNNt6SJr+ioezIPSLvCPU1XnsqqWdj14vXVtuWRz00039X2grHQy6Tf98sG8Y8cONm/eXOig4KsZ9qDlvf6/5cBfq9VOCvztpo+ybvfHXpVhwPqTZHnYoOiXHygrh0E8zNPtw0ZZ0o9C53zD7lCu3//XbtuP24N+w9guw9z25NjpnNnAiszdtptkST8KjfLDrlbp9/+1q9Ich6rOZsNoxytjlZoDQUmUvU57kDtvtyeTTulHoVF+2MFqFIJjWfQb3Dody6XcFlmKDUW/xr1qaBQGuhiFZRw1w65WGfT/63f+41D9lPU4GVaeyGtgmjK8xj0QjEKddsRoHZA2XHm0Qaz0+1G7ACnbsZw1EPiuoRIo890XzXc3zMzMMDc3V8qqljIp3R0hHeSxvP3ebdPu96vNt4zruczH8qqyRIuiX+NeIogo59X2qF2NlUHR66yfrpL7Wd5hlwiKXs+rKdOxjO8aGi1lvPuijHc3lF2R66yXO7vyWt5eG+g7dVXdbr5l3jfLeCx34kBgbZXy7oaSK3Kd9XJyzHN5uz0BZh0HeKX5et/MlwOBtTUKt2GWTZHrrJeTY5HL289VfRX2zWHeUu6BaaxvZX8GokpGaVtUtV+fLPJaNx6YxlaV1wlj1A/mUTpxZpFn/fSg100Vrup7New2EAeCCsrz5F3mRrtOVlsP4xYgWnXK37AC/Cg2rA7DsNtAHAgqKM+T9yg32rVbD+Pe1XWW/I1ygC9C3vvCsEtLDgQVlPedIqNavG+3HrKcBFtPprfeeiuvvfbaSKyDLPkb5QA/bIPskHFY+5IDwRjqdHWS98l7VIv37dZD60lw7dq17Ny5s+197G+//Tbbt2/n+PHjI9FOkuUkP8oBvlXz8QD5DyA/FqWnLE+dFf2qwpPFeSnzE5ejZPnp0DvuuKPjk61r1qyJU045pTT9y2RR5NOvgxrXol2HdcvbaWJiIiYnJ3M/Nsp8zJFnp3PAFuAQMA/csML3k8BD6ft9wIam7+bS9EPAx7LOs/nlQJBd2Tq9GnWrrc9OwcLebRAnzdXm2bz9JIWkgRwbZepWolnWQNCxakhSDbgd+CiwADwpaVdEPNeU7HrgjYg4X9I24BbgU5I2AtuAC4GfBf5c0t9Jv+k0T+uB63bztdr6bK4Su+iii8aiGmXQBlGNsto8m7dfrVZDEkePHs18bGRtBB7V6tFlWdoILgHmI+IwgKQHga1A80l7K/CF9P4R4DZJStMfjIi3ge9Lmk/zI8M8rQfjVLdbBlnX56ifCIZlEBcqnYJ18/aD7G0Eo/6MTDeyBIJ1wEtNnxeAS9uliYijkt4E1qbpf9Xy23Xpfad5nrC4uMj09DsPx83OzjI7O5th0avJJ6V8eX3mZxAXKp3m2br9sv7PsWgEzihLINAK01r7pWiXpt30lcZBaNvXxdTUFO5iwmw8DCKwDmKeVapmzRIIFoBzmj6vB460SbMgaQ1wOvB6h992mqeZWWGqVM2aJRA8CVwg6VzgZRqNv59uSbMLuA6oA1cBj0dESNoFPCDpd2g0Fl8A/A8aJYVO8zQzK1RVqgU7BoJU578d2A3UgHsi4qCkG2ncmrQLuBu4PzUGv07jxE5K9zCNRuCjwGcj4hjASvPMP3tmZtaJu6E2MxtTWbuh9uD1ZmYVN5aB4M477yx6EQrhfFeL810tg8y3A8EYcb6rxfmuFgcCMzMbmJFoLJa0CPygi5+cBbw6oMUpM+e7Wpzvaukl3z8XEVOdEo1EIDAzs8Fx1ZCZWcU5EJiZVZwDgZlZxTkQmJlV3FgFAklbJB2SNC/phqKXp1+SzpH0LUnPSzoo6V+m6WdK+jNJ30t/P5CmS9KXU/6/I+kjTfO6LqX/nqTrispTNyTVJD0l6Rvp87mS9qU8PCRpIk2fTJ/n0/cbmuYxl6YfkvSxYnKSnaQzJD0i6btpu89UYXtL+ldpH39W0h9IOnVct7ekeyS9IunZpmm5bWNJF0t6Jv3my2mQsNVlGc9yFF40Oq97ATgPmAC+DWwsern6zNPZwEfS+/cB/xPYCPxH0jjPwA3ALen9lcA3afTuehmwL00/Ezic/n4gvf9A0fnLkP/PAw8A30ifHwa2pfdfBf55ev8Z4Kvp/TbgofR+Y9oPJoFz0/5RKzpfHfJ8L/BP0/sJ4Ixx3940Bqv6PnBa03b+R+O6vYG/D3wEeLZpWm7bmEYPzzPpN98Erui4TEWvlBxX7gywu+nzHDBX9HLlnMc/pjHO8yHg7DTtbOBQen8HcHVT+kPp+6uBO5qmn5SujC8aY1TsAX4Z+EbaqV8F1rRubxq92M6k92tSOrXuA83pyvgC3p9OiGqZPtbbm3dGODwzbb9vAB8b5+0NbGgJBLls4/Tdd5umn5Su3WucqoZWGlJzXZu0IycVfz8M7AN+JiJ+BJD+fjAla7cORnHd3Ar8G+B4+rwW+HFEHE2fm/Nw0lCpQPNQqaOU7/OAReDrqUrsa5Ley5hv74h4GfhPwA+BH9HYfgcY/+3dLK9tvC69b52+qnEKBFmG1BxJkv428CjwuYj4P6slXWHaakOGlpKkjwOvRMSB5skrJI0O341Uvmlc3X4E+EpEfBj4CY1qgnbGIt+pPnwrjeqcnwXeC1yxQtJx295ZdJvXntbBOAWCLENqjhxJ76ERBP5zRPxhmvy/JZ2dvj8beCVNb7cORm3d/CLwCUkvAg/SqB66FThDjaFQ4eQ8nMifsg+VWkYLwEJE7EufH6ERGMZ9e/8D4PsRsRgRfwP8IfD3GP/t3SyvbbyQ3rdOX9U4BYITQ2qmuwu20RhCc2Sl1v67gecj4neavloeGpT094+bpl+b7jS4DHgzFTN3A5dL+kC6+ro8TSuliJiLiPURsYHGdnw8Iq4BvkVjKFR4d76X18eJoVLT9G3pLpNzeWeo1FKKiP8FvCTp59OkzTRG9xvr7U2jSugySX8r7fPL+R7r7d0il22cvntL0mVpXV7bNK/2im40ybkB5koad9a8APxW0cuTQ35+iUax7jvA0+l1JY360D3A99LfM1N6Aben/D8DTDfN658A8+n1j4vOWxfrYBPv3DV0Ho0Dex74L8Bkmn5q+jyfvj+v6fe/ldbHITLcPVH0C/gQsD9t8z+icUfI2G9v4N8D3wWeBe6ncefPWG5v4A9otIX8DY0r+Ovz3MbAdFqPLwC30XLzwUovdzpnZlZx41Q1ZGZmPXAgMDOrOAcCM7OKcyAwM6s4BwIzs4pzIDAzqzgHAjOzivv/Wzi/M1JkKD0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, errors, 'k.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the number of samples increases, there is a definite decreasing trend in the fractional error. This makes sense - since we're randomly sampling, more samples means better accuracy. Since the method depends on ratio of samples found inside to outside, more samples increases the odds that we won't be thrown off by randomly grabbing, say, a single point outside."
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}