notebooks/normal_errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstrate How a Normal Distribution Can Occur Naturally\n",
    "*This notebook shows how error terms that break down as sums of non-normal random variables can be normally distributed.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import load_boston\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import root_scalar\n",
    "from numpy.polynomial.hermite import hermgauss\n",
    "from collections import defaultdict\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decorator for Vector Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vectorizable(f):\n",
    "    def f_(x):\n",
    "        if hasattr(x, '__iter__'):\n",
    "            return np.array([f(xi) for xi in x])\n",
    "        return f(x)\n",
    "    return f_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a Non-Normal Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll make a cubic spline that gradually transitions between two different linear functions with different slopes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_spline(fpa, fpb, delta):\n",
    "    mid = delta / 2\n",
    "    fpmid = (fpb + fpa) / 2\n",
    "    b1 = fpa\n",
    "    d1 = (fpmid - b1)/(3*mid**2)\n",
    "    \n",
    "    def q1(x):\n",
    "        return b1*x + d1*x**3\n",
    "    def q1p(x):\n",
    "        return b1 + 3*d1*x**2\n",
    "\n",
    "    qmid = q1(mid)\n",
    "    \n",
    "    b2 = fpb\n",
    "    d2 = (fpmid - b2) / (3*mid**2)\n",
    "    a2 = qmid + b2*mid + d2*mid**3\n",
    "    \n",
    "    def q2(x):\n",
    "        return a2 + b2*(x - delta) + d2*(x - delta)**3\n",
    "    def q2p(x):\n",
    "        return b2 + 3*d2*(x - delta)**2\n",
    "    \n",
    "    def q(x):\n",
    "        if x < mid:\n",
    "            return q1(x)\n",
    "        return q2(x)\n",
    "    def qp(x):\n",
    "        if x < mid:\n",
    "            return q1p(x)\n",
    "        return q2p(x)\n",
    "    return q, qp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, we'll use the cubic splines to construct a function that transitions between multiple different linear functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_warp(descriptor, delta=0.5):\n",
    "    slope = descriptor[0]\n",
    "    rest = descriptor[1:]\n",
    "    if not rest:\n",
    "        return lambda x: x*slope, lambda x: slope\n",
    "    z, rslope = rest[0]\n",
    "    r, rp = make_warp([rslope] + rest[1:], delta)\n",
    "    c = (rslope - slope) / (2 * delta)\n",
    "    intercept = slope*z + c*delta - r(z)\n",
    "    q, qp = compute_spline(slope, rslope, delta)\n",
    "    @vectorizable\n",
    "    def f(x):\n",
    "        if x < z - delta:\n",
    "            return slope * x\n",
    "        if x > z:\n",
    "            return r(x) + intercept\n",
    "        return q(x - (z-delta))+slope*(z - delta)\n",
    "    @vectorizable\n",
    "    def fp(x):\n",
    "        if x < z - delta:\n",
    "            return slope\n",
    "        if x > z:\n",
    "            return rp(x)\n",
    "        return qp(x - (z-delta))\n",
    "    return f, fp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a Non-Normal Base Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, fp = make_warp([0.25, (-1.5, 1), (-0.25, 5)], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5b3H8c+PrCyBEAhrQEAWZZEtiVirdd/aul2luJKIolVbvbet1a72am/tXrVVi0oSFMW9UGvrQrXWtkrCvktAFBBI2EJYsj/3jzloCAlJyMycmcz3/XrxysyZk5xvzhy+efLMk8Scc4iISPTp4HcAERE5NipwEZEopQIXEYlSKnARkSilAhcRiVLx4TxYz5493aBBg8J5SBGRqLdw4cIdzrn0htvDWuCDBg2iqKgonIcUEYl6ZvZxY9s1hSIiEqVU4CIiUUoFLiISpVTgIiJRSgUuIhKlVOAiIlFKBS4iEqVU4CIiIbRjXyU/+fNKKqprg/6xVeAiIiFyoKqGaQVFPLvgE9aX7gv6xw/rT2KKiMSK2jrHN59dwrLNe3js2omM6tct6MfQCFxEJMicc/zkzyt5a/V2fvyVkZw/qk9IjqMCFxEJshnvbmDWfz7mptMGk3Pq4JAdRwUuIhJEf176KT/76xq+PKYv91x4YkiPpQIXEQmSDzbs5FvPLyVrUHd+PXksHTpYSI+nAhcRCYLiknJumlVERlpHHr8+k+SEuJAfUwUuItJGJeUVTJ1ZSGJ8Bwpys0ntlBiW42oZoYhIG+yvrOGG/EJ27a/iuZsnMSCtU9iOrRG4iMgxqqmt4/ZnFrHq0738/urxnJSRGtbjawQuInIMnHP8cO4K3l5byk8vG83ZJ/YOewaNwEVEjsEj76zn2QWb+PoZx3PNycf5kkEFLiLSSq8s3swvX1/LJeP68Z3zRviWo0VTKGa2ESgHaoEa51ymmaUBzwGDgI3AZOfc7tDEFBGJDP8u3sFdLy5j0pA0fnHFSSFf6300rRmBn+mcG+ecy/Tu3w3Md84NA+Z790VE2q2128q5+amFDOrRmT9el0lSfOjXeh9NW6ZQLgEKvNsFwKVtjyMiEpm2lVWQk7eAjolx5N+QTbeOCX5HanGBO+ANM1toZtO9bb2dc1u929uARl+CNbPpZlZkZkWlpaVtjCsiEn7lFdXk5hey92A1eblZ9E/t6HckoOXLCL/onNtiZr2AN81sTf0HnXPOzFxj7+icmwHMAMjMzGx0HxGRSFVdW8etsxfx4fZyZuZkheT3eh+rFo3AnXNbvLclwCtANrDdzPoCeG9LQhVSRMQPzjnueXk5/1y3g59dNoYvDU/3O9Jhmi1wM+tsZimHbgPnASuAecBUb7epwNxQhRQR8cOD89fx4sLNfPPsYUzOGuB3nCO0ZAqlN/CKmR3a/xnn3N/MrBB43symAR8Dk0MXU0QkvJ4v2sTv3lrHf03I4L/PGeZ3nEY1W+DOuQ3A2Ea27wTODkUoERE/vfthKd97eTlfHNqTn10+Bm8AG3H0k5giIvWs/LSMrz+9kKG9uvDotRNIjI/cmozcZCIiYfbpnoPckF9I144J5Odmk5Ls/1rvo1GBi4gAZQeryclbwIHKWvJys+jTLdnvSM3Sr5MVkZhXVVPHLU8t5KMd+8nPzeaEPl39jtQiKnARiWnOOb770jL+s2Env5k8llOH9vQ7UotpCkVEYtqv3ljLK4u38K1zh3P5hAy/47SKClxEYtYzH3zCH95ez5SsAdx+1lC/47SaClxEYtLba0r44dwVfGl4OvdfOjpi13ofjQpcRGLO8s1l3PbMIk7ok8IfrplAfFx0VmF0phYROUabdh0gN7+Q7p0SycvJoktS9K7liN7kIiKttOdAFTl5C6iqqeXZm06mV9fIX+t9NCpwEYkJFdW1TJ+1kE27DjJrWjbDeqf4HanNVOAi0u7V1Tm+/cJSFmzcxUNXjWfSkB5+RwoKzYGLSLv387+t4dVlW7n7whO4eGw/v+MEjQpcRNq1Wf/ZyB/f3cC1kwZy8+lD/I4TVCpwEWm33li5jXvnreScE3tx71dHReVa76NRgYtIu7T4k918c85ixvTvxkNXjY/atd5H0/4+IxGJeR/v3M+NBUWkpyTxxNQsOiW2z/UaKnARaVd27a8iJ6+QWufIz80mPSXJ70gh0z6/LIlITKqoruXGgkK27DnIMzeezPHpXfyOFFIagYtIu1Bb57hzzhIWb9rDg18bR+agNL8jhZwKXETahfv/soq/rdzG9y86kQvH9PU7TliowEUk6j353kfk/WsjuacO4sbT2tda76NRgYtIVPvr8q3c/5dVXDCqDz/48ki/44SVClxEolbRxl3c8dwSxg9I5XdTxhHXoX39oE5zWlzgZhZnZovN7FXv/mAz+8DMis3sOTNLDF1MEZHDbSjdx42ziuif2pEnpmaRnBDnd6Swa80I/A5gdb37Pwd+65wbCuwGpgUzmIhIU3bsqyQnr5A4M/Jzs0jrHJvjxxYVuJllAF8GnvDuG3AW8KK3SwFwaSgCiojUd6Cqhmn5hZSUV/DE1EyO69HZ70i+aekI/HfAXUCdd78HsMc5V+Pd3wz0D3I2EZHD1NY5vvnsYpZvKeOhKeMZP7C735F81WyBm9lXgBLn3MJjOYCZTTezIjMrKi0tPZYPISKCc457563krdUl3HvxKM4b1cfvSL5ryQj8VOBiM9sIzCEwdfIgkGpmh34UPwPY0tg7O+dmOOcynXOZ6enpQYgsIrHoj+9u4Kn3P2b66UO4/pRBfseJCM0WuHPuHudchnNuEDAF+Ltz7hrgbeAKb7epwNyQpRSRmDZv6ac88Nc1fOWkvtx9wQl+x4kYbVkH/l3gf8ysmMCc+JPBiSQi8rn3N+zk288vJXtQGr+6ciwdYmyt99G06rcROufeAd7xbm8AsoMfSUQkYN32cqbPKmJAWkdmXD8xJtd6H41+ElNEIlLJ3gpy8gpJjI8jPzeb1E6xudb7aFTgIhJx9lfWcENBIbv2V5GXk8WAtE5+R4pIKnARiSg1tXXc9swiVn26lz9cM54xGd38jhSx9Bd5RCRiOOf4wZ9W8M7aUv7vsjGcdUJvvyNFNI3ARSRi/OHtYuYUbuK2M4/n6pMH+h0n4qnARSQivLxoM79640MuHdePb583wu84UUEFLiK++1fxDu56cRmnDOnBL64YS+D35UlzVOAi4qs12/Zyy1MLGZLemceum0hivGqppXSmRMQ3W8sOkptXSKekOPJys+nWMcHvSFFFBS4iviivqCY3r5C9B6uZmZNF/9SOfkeKOlpGKCJhV11bx62zF1Fcso+ZOVmM6qe13sdCBS4iYeWc4+6XlvPPdTv4xRUncfpw/ZrpY6UpFBEJq9+9tY6XFm3mjrOHMTlzgN9xopoKXETC5vnCTTw4fx1XTszgznOG+R0n6qnARSQs/vFhKfe8spzThvXk/y4fo7XeQaACF5GQW/lpGbc+vZDhvVN45JoJJMSpeoJBZ1FEQmrLnsBa764dE8jLySIlWWu9g0UFLiIhU3awmty8BRysqiU/N5s+3ZL9jtSuaBmhiIREZU0tNz9VxEc79lOQm82IPil+R2p3VOAiEnTOOb774jLe37CL335tLF8Y2tPvSO2SplBEJOh++fpa/rTkU75z/gguG5/hd5x2SwUuIkE1+4OPeeSd9VyVPYBbzzje7zjtmgpcRILm72u288M/reDMEencd8lorfUOMRW4iATFss17uG32Ykb268rvr55AvNZ6h5zOsIi02aZdB7ghv5C0zonMzMmic5LWR4SDzrKItMmeA1VMzVtAda1jzvQseqVorXe4NDsCN7NkM1tgZkvNbKWZ/cTbPtjMPjCzYjN7zswSQx9XRCJJRXUtN80qYvOug8y4biJDe2mtdzi1ZAqlEjjLOTcWGAdcYGaTgJ8Dv3XODQV2A9NCF1NEIk1dneNbLyylcONufjV5LCcP6eF3pJjTbIG7gH3e3QTvnwPOAl70thcAl4YkoYhEpAf+toa/LNvKPReewMVj+/kdJya16EVMM4szsyVACfAmsB7Y45yr8XbZDPRv4n2nm1mRmRWVlpYGI7OI+Kzg3xuZ8e4Grj/lOKafPsTvODGrRQXunKt1zo0DMoBs4ISWHsA5N8M5l+mcy0xP159OEol2b6zcxr1/Xsk5J/bmx18dpbXePmrVMkLn3B7gbeAUINXMDq1iyQC2BDmbiESYxZ/s5ptzFnNSRioPXzWeuA4qbz+1ZBVKupmlerc7AucCqwkU+RXeblOBuaEKKSL+27hjP9MKiuiVksyTUzPpmBjnd6SY15J14H2BAjOLI1D4zzvnXjWzVcAcM7sfWAw8GcKcIuKjXfuryMlbgHOO/NwsenZJ8juS0IICd84tA8Y3sn0DgflwEWnHKqprubGgkK1lFTxz08kMSe/idyTx6CcxRaRJtXWOO+YsZvGmPTxy9QQmHpfmdySpR78LRUQa5ZzjvldX8frK7fzgyyO5cExfvyNJAypwEWnUk+99RP6/N3LDqYOZ9sXBfseRRqjAReQIry3fyk9fW82Fo/vwgy+f6HccaYIKXEQOU7hxF3c+t4QJA7vz26+No4PWekcsFbiIfGZ96T5umlVE/9SOPH59JskJWusdyVTgIgJAaXklOXkLiDMjPzeLtM76DdGRTssIRYQDVTVMKyiktLySOdNP4bgenf2OJC2gEbhIjKupreMbzyxmxZYyHr5qAuMGpPodSVpII3CRGOac494/r2T+mhLuu2QU547s7XckaQWNwEVi2GP/2MDT73/CzV8awnWnDPI7jrSSClwkRs1dsoWf/20NXx3bj++e3+Jf8S8RRAUuEoPe37CT77ywjOzBafzqypO01jtKqcBFYsy67eVMn1XEwB6dePy6TJLitdY7WqnARWJIyd4KcvIKSUqIIy8ni26dEvyOJG2gAheJEfsqa8jNL2T3gSrycrIYkNbJ70jSRlpGKBIDqmvruG32ItZsK+eJqZmM7t/N70gSBBqBi7Rzzjl++KcV/OPDUu6/dDRnjujldyQJEhW4SDv3+78XM6dwE7efOZSrsgf6HUeCSAUu0o69tHAzv37zQy4f359vnTfc7zgSZCpwkXbqX8U7+O5Ly/jC8T144L9OwkxrvdsbFbhIO7Rm215ueWohx6d34bHrJpIYr//q7ZGeVZF2ZmvZQXJmFtIpKY683Cy6Jmutd3ulZYQi7cjeimpy8wrZV1nD8zefQr/Ujn5HkhDSCFyknaiqqePWpxdRXLKPR6+dwMh+Xf2OJCGmEbhIO+Cc4+6Xl/Fe8Q5+ecVJnDYs3e9IEgbNjsDNbICZvW1mq8xspZnd4W1PM7M3zWyd97Z76OOKSGN++9Y6Xl60hTvPGcaVmQP8jiNh0pIplBrgW865kcAk4DYzGwncDcx3zg0D5nv3RSTMniv8hIfmr2NyZgZ3nD3M7zgSRs0WuHNuq3NukXe7HFgN9AcuAQq83QqAS0MVUkQa987aEr73ygpOH57OTy8bo7XeMaZVL2Ka2SBgPPAB0Ns5t9V7aBvQ6B/TM7PpZlZkZkWlpaVtiCoi9a3YUsZtsxcxoncKj1wzgYQ4rUmINS1+xs2sC/AScKdzbm/9x5xzDnCNvZ9zboZzLtM5l5merhdWRIJhy56D3JBfSLeOCeTlZtElSesRYlGLCtzMEgiU92zn3Mve5u1m1td7vC9QEpqIIlJf2YFqcmYu4GB1Lfk3ZNO7a7LfkcQnLVmFYsCTwGrn3G/qPTQPmOrdngrMDX48EamvsqaWm58uYuPO/fzxuokM753idyTxUUu+7zoVuA5YbmZLvG3fAx4AnjezacDHwOTQRBQRgLo6x10vLuP9Dbv43dfG8YXje/odSXzWbIE7594Dmnpp++zgxhGRpvzyjbXMXfIp3zl/BJeO7+93HIkAetlaJAo8/f7HPPrOeq4+eSC3nnG833EkQqjARSLc/NXb+dHcFZx1Qi/+9+JRWustn1GBi0SwpZv2cPszixnVrxsPXzWeeK31lnp0NYhEqE27DjCtoJAeXRJ5MieTzlrrLQ3oihCJQLv3VzE1bwHVtY4507PplaK13nIkFbhIhKmoruWmWUVs3n2Qp6edzNBeXfyOJBFKUygiEaSuzvGt55dS9PFufjN5LNmD0/yOJBFMBS4SQX7219X8ZflWvnfRCXzlpH5+x5EIpwIXiRD5//qIx//5EVNPOY6bThvidxyJAipwkQjw+spt/OTVVZw7sjc/+qrWekvLqMBFfLbok91889nFjM1I5aEp44nroPKWllGBi/ho44793FhQRJ9uyTw5NZOOiXF+R5IoogIX8cnOfZXk5C3AOUd+bjY9uiT5HUmijNaBi/jgYFUtN84qYmtZBc/cNInBPTv7HUmikApcJMxq6xx3zFnMkk17ePSaiUw8rrvfkSRKaQpFJIycc9z36ireWLWdH31lJBeM7uN3JIliKnCRMHryvY/I//dGpn1xMLmnDvY7jkQ5FbhImPxl2Vbu/8tqLhrTh+9fdKLfcaQdUIGLhEHhxl389/NLyDyuO7+ZPI4OWustQaACFwmx9aX7uLGgiIzUjjx+fSbJCVrrLcGhAhcJodLywFrvhDgjPzeb7p0T/Y4k7YiWEYqEyIGqGqYVFLKjvIo50ycxsEcnvyNJO6MRuEgI1NTW8Y1nFrNiSxkPXzWesQNS/Y4k7ZBG4CJB5pzjx/NWMn9NCfddOppzRvb2O5K0UxqBiwTZo/9Yz+wPPuGWLx3PdZOO8zuOtGPNFriZzTSzEjNbUW9bmpm9aWbrvLf6WWARYO6SLfzib2u5eGw/7jp/hN9xpJ1ryQg8H7igwba7gfnOuWHAfO++SEz7z/qdfPuFpZw8OI1fXnmS1npLyDVb4M65d4FdDTZfAhR4twuAS4OcSySqfLi9nOlPFTGoR2dmXJdJUrzWekvoHesceG/n3Fbv9jZAr9JIzNq+t4KcmQtITogjLzeLbp0S/I4kMaLNL2I65xzgmnrczKabWZGZFZWWlrb1cCIRZV9lDbl5hew5WE1eThYZ3bXWW8LnWAt8u5n1BfDeljS1o3NuhnMu0zmXmZ6efoyHE4k81bV13Dp7EWu3l/PINRMY3b+b35Ekxhxrgc8Dpnq3pwJzgxNHJDo45/jBKyt498NSfnrpaM4Y0cvvSBKDWrKM8FngP8AIM9tsZtOAB4BzzWwdcI53XyRmPPz3Yp4r2sQ3zhrKlOyBfseRGNXsT2I6565q4qGzg5xFJCq8uHAzv3nzQy6f0J//OXe433EkhuknMUVa4b11O7j7pWWcOrQHD1x+EmZa6y3+UYGLtNDqrXu55emFDO3VhUevnUhivP77iL90BYq0wNayg+TmFdIlKZ683Cy6Jmutt/hPv41QpBl7K6rJzStkX2UNL9xyCn27dfQ7kgigEbjIUVXV1PH1pxdSXLKPR6+dwIl9u/odSeQzGoGLNME5x90vL+NfxTv51ZVjOW2YfhBNIotG4CJN+O2bH/Lyoi38z7nDuWJiht9xRI6gAhdpxJwFn/DQ34v5WuYAvnHWUL/jiDRKBS7SwNtrS/j+n1Zw+vB07r9stNZ6S8RSgYvUs2JLGbfNXsSI3ik8cs0EEuL0X0Qil65OEc/m3QfIzS+ke6dE8nKz6JKk1/glsukKFQHKDlSTk1dIRXUts288md5dk/2OJNIsjcAl5lXW1DL9qSI+2XmAGddlMrx3it+RRFpEI3CJaXV1ju+8sIwPPtrFg1PGccrxPfyOJNJiGoFLTPvF62uZt/RT7rpgBJeM6+93HJFWUYFLzHrq/Y957B/ruXbSQL7+peP9jiPSaipwiUlvrdrOj+eu4OwTenHvV0dprbdEJRW4xJylm/bwjWcXM7p/Nx6+ejzxWustUUpXrsSUj3bsJze/kJ4piTw5NYtOiXodX6KXClxiRkl5BdfP/ACAgtxs0lOSfE4k0jYqcIkJ5RXV5MwsZEd5FXk5WQxJ7+J3JJE2U4FLu7evsoYb8gv5cHs5j147gbEDUv2OJBIUmgCUdm1fZQ05MxeweNMeHpwyjjNG9PI7kkjQqMCl3dpWVsGNswpZvbWch68az0Vj+vodSSSoVODSLv1n/U7ufG4x+ypqePz6iZx1Qm+/I4kEnQpc2pVtZRU8OH8dzy74hEE9OpGfm60/RCztVpsK3MwuAB4E4oAnnHMPBCWVSCtU1tTy/oZdvLJoM6+t2EZdnSP31EF85/wRWuct7doxX91mFgf8ATgX2AwUmtk859yqYIUTaczeimoWfbybwo27KPxoN0s276Gqpo6uyfFMzszg5tOPZ0BaJ79jioRcW4Yn2UCxc24DgJnNAS4BVOASVNv3VlC0MVDYCz7axZpte6lzENfBGN2vK9dPOo6Th/Tg9OE9SYqP8zuuSNi0pcD7A5vq3d8MnNxwJzObDkwHGDhwYBsOJ+1dbZ3j0z0HWV+6jxVbyli6uYxlm/ewfW8lAB0T4hg/MJVvnDWM7MFpjB+YqikSiWkhv/qdczOAGQCZmZku1McTf9XVOcorayivqGZfZQ3lFYHb5RU17K2oYV+9+5+/rWHn/ko27TpIVW3dZx9rSHpnThnSgzEZqUw8rjuj+nXVHxkWqactBb4FGFDvfoa3TaJUdW3dEcVav2wPFfLeJraXV9Swr7Km2ePEdTBSkuMD/5ISSEmOZ1ivFM4Z2ZvBPTozuGdnTuzXla7JCWH4rEWiV1sKvBAYZmaDCRT3FODqoKSSVnHOUVlTx16vVPdVNBz5Hjka/ux2ve0V1XXNHisxvgNdk+NJSU74rIR7pSTT5VAhJyd4j8fTJenzfQ5t75IcT8eEOP3+bZEgOOYCd87VmNntwOsElhHOdM6tDFqyGOGcY39VbWAke8TotuHI9/Cph/qlXF3b/OxU58Q4r2gDxdqtUyIZaZ1ISfq8ZA97m3R4UXdJjteLhCIRpE1z4M6514DXgpQl6tTWOa90GxvhVgfmfCubmpL4vJzrmuneDgZdGpRp767JDG0w0q0/Mq6/f9fkBLokxxPXQaNekfYkZl/Cr6qpO2wku7eRkj1Uvnsbbvdu76+qbfY4CXHWoFTjGZDW6fNiPWL0e+TtzomachCRI0VdgTvnqKiua2KE28Rqh8rPi/fQFEVlTfPzvckJHY6YTujTNfmwku2SFCjiw7Z55ds1OYGk+A4qXxEJiago8O+9spz31u34rJBrmptzgHoj20Cxdu+UyMC0Toe9yJbSyAj40HRDSnK8lqyJSESLigLvn9qRCQNTj3yR7Ygph8Dbzoma7xWR9i8qCvy2M4f6HUFEJOJojkBEJEqpwEVEopQKXEQkSqnARUSilApcRCRKqcBFRKKUClxEJEqpwEVEopQ5F74/kmNmpcDHx/juPYEdQYwTLMrVOsrVOsrVOu0113HOufSGG8Na4G1hZkXOuUy/czSkXK2jXK2jXK0Ta7k0hSIiEqVU4CIiUSqaCnyG3wGaoFyto1yto1ytE1O5omYOXEREDhdNI3AREalHBS4iEqUiqsDN7EozW2lmdWaW2eCxe8ys2MzWmtn5Tbz/YDP7wNvvOTNLDEHG58xsifdvo5ktaWK/jWa23NuvKNg5GjnevWa2pV62i5rY7wLvHBab2d1hyPVLM1tjZsvM7BUzS21iv7Ccr+Y+fzNL8p7jYu9aGhSqLPWOOcDM3jazVd71f0cj+5xhZmX1nt8fhTqXd9yjPi8W8JB3vpaZ2YQwZBpR7zwsMbO9ZnZng33Ccr7MbKaZlZjZinrb0szsTTNb573t3sT7TvX2WWdmU48pgHMuYv4BJwIjgHeAzHrbRwJLgSRgMLAeiGvk/Z8Hpni3HwO+HuK8vwZ+1MRjG4GeYTx39wLfbmafOO/cDQESvXM6MsS5zgPivds/B37u1/lqyecP3Ao85t2eAjwXhueuLzDBu50CfNhIrjOAV8N1PbX0eQEuAv4KGDAJ+CDM+eKAbQR+0CXs5ws4HZgArKi37RfA3d7tuxu75oE0YIP3trt3u3trjx9RI3Dn3Grn3NpGHroEmOOcq3TOfQQUA9n1d7DAn34/C3jR21QAXBqqrN7xJgPPhuoYIZANFDvnNjjnqoA5BM5tyDjn3nDO1Xh33wcyQnm8ZrTk87+EwLUDgWvpbO+5Dhnn3Fbn3CLvdjmwGugfymMG0SXALBfwPpBqZn3DePyzgfXOuWP9Ce82cc69C+xqsLn+NdRUD50PvOmc2+Wc2w28CVzQ2uNHVIEfRX9gU737mznyAu8B7KlXFo3tE0ynAdudc+uaeNwBb5jZQjObHsIc9d3ufRs7s4lv21pyHkPpBgKjtcaE43y15PP/bB/vWiojcG2FhTdlMx74oJGHTzGzpWb2VzMbFaZIzT0vfl9TU2h6EOXH+QLo7Zzb6t3eBvRuZJ+gnLew/1FjM3sL6NPIQ993zs0Nd57GtDDjVRx99P1F59wWM+sFvGlma7yv1iHJBTwK3EfgP9x9BKZ3bmjL8YKR69D5MrPvAzXA7CY+TNDPV7Qxsy7AS8Cdzrm9DR5eRGCaYJ/3+safgGFhiBWxz4v3GtfFwD2NPOzX+TqMc86ZWcjWaoe9wJ1z5xzDu20BBtS7n+Ftq28ngW/f4r2RU2P7BCWjmcUDlwMTj/IxtnhvS8zsFQLfvrfpwm/puTOzx4FXG3moJecx6LnMLAf4CnC28yYAG/kYQT9fjWjJ539on83e89yNwLUVUmaWQKC8ZzvnXm74eP1Cd869ZmaPmFlP51xIf3FTC56XkFxTLXQhsMg5t73hA36dL892M+vrnNvqTSeVNLLPFgLz9IdkEHjtr1WiZQplHjDFWyEwmMBX0gX1d/CK4W3gCm/TVCBUI/pzgDXOuc2NPWhmnc0s5dBtAi/krWhs32BpMO94WRPHKwSGWWC1TiKBbz/nhTjXBcBdwMXOuQNN7BOu89WSz38egWsHAtfS35v6ohMs3hz7k8Bq59xvmtinz6G5eDPLJvB/N6RfWFr4vMwDrvdWo0wCyupNH4Rak98F+3G+6ql/DTXVQ68D55lZd2+68zxvW+uE+lXaVr6iexmBuaBKYDvwer3Hvk9gBcFa4MJ6218D+sXyQqMAAADpSURBVHm3hxAo9mLgBSApRDnzgVsabOsHvFYvx1Lv30oCUwmhPndPAcuBZd4F1LdhLu/+RQRWOawPU65iAnN9S7x/jzXMFc7z1djnD/wvgS8wAMnetVPsXUtDwnCOvkhg6mtZvfN0EXDLoesMuN07N0sJvBj8hTDkavR5aZDLgD9453M59VaPhThbZwKF3K3etrCfLwJfQLYC1V53TSPwmsl8YB3wFpDm7ZsJPFHvfW/wrrNiIPdYjq8fpRcRiVLRMoUiIiINqMBFRKKUClxEJEqpwEVEopQKXEQkSqnARUSilApcRCRK/T9306qRbZZo5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVPklEQVR4nO3de3Bcd3nG8efV3ZYlS7bkW+xESUOcC4UkaAIlAdIEkpBmoOnQTjLTlgJTl3aYgWk7DBlmGCj9BzplGDotwYVMKFACBdKmGQIxJXG4JWAHJ9hO7DiJg+3YkWxL1mWlvUhv/9ijZCOvpF17z54953w/Mxrv7jm7++rs8aPfvudm7i4AQONqiroAAMDiCGoAaHAENQA0OIIaABocQQ0ADa4ljBft6+vzgYGBMF4aABJp586dx929v9y0UIJ6YGBAO3bsCOOlASCRzOyFhabR+gCABkdQA0CDI6gBoMER1ADQ4AhqAGhwFe31YWYHJY1LmpFUcPfBMIsCALyimt3zft/dj4dWCQCgrFD2owYawc4XTmr7/uMSp/JFnSxvb9EH3/Y7NX/dSoPaJT1oZi7pS+6+df4MZrZF0hZJOvfcc2tXIXAGDp3M6Patjyk3MyuzqKtBWvStaI80qK9x9yNmtkbSNjN72t0fKZ0hCO+tkjQ4OMgQBpH6+qMvyOX6xR3Xaf3KZVGXA5yVivb6cPcjwb9Dku6VdFWYRQFn66F9Q3rj+asJaSTCkkFtZp1m1jV3W9INknaHXRhwpobGprX/pQm99aK+qEsBaqKS1sdaSfdasdHXIuk/3f0HoVYFnIU9L45Jki7f1BtxJUBtLBnU7v6cpNfXoRagJvYeLQb1xeu7Iq4EqA2OTETi7D06pk2rlqm7ozXqUoCaIKiROPuPjWvz2u6oywBqhqBGori7Do1kNLB6edSlADVDUCNRhieyms7P6lyCGglCUCNRDp3MSJI2rSKokRwENRLl0MkpSdKmXoIayUFQI1F+G4yoN/ZyRCKSg6BGohw6mdHa7nZ1tDZHXQpQMwQ1EuXY2LTWcX4PJAxBjUQZGstqTVd71GUANUVQI1GGxqcJaiQOQY3EyBVmNZLJa213R9SlADVFUCMxhieyksSIGolDUCMxXhqbliSt6SaokSwENRJjaGxuRE3rA8lCUCMxhseDETWtDyQMQY3EGBrPqsmk1SsIaiQLQY3EGB7PalVnu5qbLOpSgJoiqJEYJydzWt3ZFnUZQM0R1EiM0UxePcu5/BaSh6BGYpzM5LSKETUSiKBGYoxmcuolqJFABDUSYXbWNZLJq5fWBxKIoEYijE8XNDPr6l3OiBrJQ1AjEUYyOUmiR41EIqiRCCeDoGZEjSQiqJEIo3NBzYgaCURQIxFOTuYlSasYUSOBCGokwshkcUTd08leH0geghqJMJLJqaXJ1NXeEnUpQM0R1EiEkUxOPcvbZMYJmZA8BDUSgfN8IMkIaiTC2HReK5cR1EimioPazJrN7Ndmdn+YBQFnYmyqoO4O+tNIpmpG1B+W9FRYhQBnY3w6r64ORtRIpoqC2sw2SvoDSV8OtxzgzIxNF9S9jBE1kqnSEfXnJX1U0uxCM5jZFjPbYWY7hoeHa1IcUAl319hUXt2MqJFQSwa1md0iacjddy42n7tvdfdBdx/s7++vWYHAUqbzsyrMurrZmIiEqmREfbWkd5nZQUn3SLrOzL4ealVAFcami4ePd7ExEQm1ZFC7+x3uvtHdByTdJunH7v6noVcGVGhsqhjUtD6QVOxHjdibG1HT+kBSVfVd0d0flvRwKJUAZ2hsuiBJ7EeNxGJEjdiba32wHzWSiqBG7L08omY/aiQUQY3YY2Miko6gRuyNTxfU1tKkjtbmqEsBQkFQI/bGpvNsSESiEdSIPQ4fR9IR1Ii98emCutiHGglGUCP2aH0g6QhqxN7YVJ7zfCDRCGrE3mR2Riu4+jgSjKBG7E1mC+okqJFgBDVizd01mSuos42gRnIR1Ii1qfyMZl2MqJFoBDVibTI7I0la0c5RiUgughqxNpktnpCJETWSjKBGrE0Q1EgBghqxNjeiZvc8JBlBjVibzDGiRvIR1Ii1CTYmIgUIasQaGxORBgQ1Ym0uqJdzwAsSjKBGrL2810cbrQ8kF0GNWJvMFtTR2qSWZlZlJBdrN2JtgjPnIQUIasQaZ85DGhDUiLXJLGfOQ/IR1Ii1yVyB1gcSj6BGrE1mZ9TJwS5IOIIasUaPGmlAUCPWJuhRIwUIasQaI2qkAUGN2JqddU3mZjghExJvyaA2sw4z+6WZPWFme8zsU/UoDFhKJl88cx4jaiRdJWt4VtJ17j5hZq2SfmpmD7j7oyHXBiyKM+chLZZcw93dJU0Ed1uDHw+zKKASE1zdBSlRUY/azJrNbJekIUnb3P2xMvNsMbMdZrZjeHi41nUCp8lkaX0gHSoKanefcffLJW2UdJWZvbbMPFvdfdDdB/v7+2tdJ3AaTnGKtKhqrw93H5X0kKSbwikHqFyG6yUiJSrZ66PfzHqC28skvUPS02EXBizl5RE1u+ch4SoZiqyX9FUza1Yx2L/t7veHWxawtEl61EiJSvb6eFLSFXWoBagKrQ+kBUcmIrbmWh/LW2l9INkIasRWJjfD9RKRCqzhiC3OnIe0IKgRWxnOnIeUIKgRWxPZGS3nYBekAEGN2JrMcr1EpANBjdjK5Gh9IB0IasTWRLbAUYlIBYIasZXJzbDXB1KBoEZsTbDXB1KCoEYsuXtwYVtaH0g+ghqxlC3MatY5zwfSgaBGLL1y0QCCGslHUCOWuAwX0oSgRixxGS6kCUGNWJrkXNRIEYIasTSZJaiRHgQ1YumVy3DR+kDyEdSIpfHpvCSpq6M14kqA8BHUiKXx6WLro6uD1geSj6BGLI1P52UmrWA/aqQAQY1YGpsuaEVbi5qaLOpSgNAR1Iil8ekCbQ+kBkGNWBqfzrMhEalBUCOWGFEjTQhqxNJ4Nk9QIzUIasRScURN6wPpQFAjlmh9IE0IasSOu7MxEalCUCN2soVZ5WecETVSg6BG7IwF5/noJqiREgQ1YueV83zQ+kA6ENSIHU7IhLRZMqjNbJOZPWRme81sj5l9uB6FAQsZm+IUp0iXSoYkBUl/5+6Pm1mXpJ1mts3d94ZcG1AWI2qkzZIjanc/6u6PB7fHJT0l6ZywCwMW8spFAwhqpENVPWozG5B0haTHykzbYmY7zGzH8PBwbaoDypgbUXcvo/WBdKg4qM1shaTvSvqIu4/Nn+7uW9190N0H+/v7a1kj8CojmZyam0xdXNgWKVFRUJtZq4oh/Q13/164JQGLG53Kq2dZq8y4aADSoZK9PkzSVyQ95e6fC78kYHGnMnmtXE7bA+lRyYj6akl/Juk6M9sV/Nwccl3Agkancupd3hZ1GUDdLNnkc/efSuI7JhrGyGRe61d2RF0GUDccmYjYOTVF6wPpQlAjdkYzOfUso/WB9CCoESu5wqwmczPqZUSNFCGoESujUzlJUg9BjRQhqBErpzLFw8dXstcHUoSgRqyMBmfOo/WBNCGoESujwYiajYlIE4IasTKSoUeN9CGoESsjk8Wg7u1kRI30IKgRK8cnsupobVJnW3PUpQB1Q1AjVk5M5LS6s50z5yFVCGrEyvBEVn1d7VGXAdQVQY1YOT6RUx/9aaQMQY1YOTGRVd8KRtRIF4IasTE76zoxmVNfFyNqpAtBjdg4NZXXzKxrdScjaqQLQY3YOD6RlSQ2JiJ1CGrExvGJ4sEubExE2hDUiI2h8WlJ0ppuRtRIF4IasfHiaDGo161cFnElQH0R1IiNY6em1NXRohXtS16TGUgUghqx8eKpaW1gNI0UIqgRG8dOTWvdyo6oywDqjqBGbBw9NaUNPQQ10oegRixkCzM6PpHTum5aH0gfghqxcHhkSpK0sZegRvoQ1IiFg8cnJUkDfZ0RVwLUH0GNWHg+COoLCGqkEEGNWHj++KRWLmvlWolIJYIasfDc8CRtD6QWQY2G5+7ae3RMl67viroUIBIENRre4ZEpnZrK67INK6MuBYgEQY2Gt+vQqCTpd88hqJFOSwa1md1lZkNmtrseBQHz/eSZYXV3tOiyDd1RlwJEopIR9d2Sbgq5DkBSsR89Pp3X8YmsXhyd0t4Xx/TA7mO6dvMatTTzBRDptOT5It39ETMbCL8UpJm764vbn9WXtj+nU1P5V01rMul9Vw9EUxjQAGp2Yl8z2yJpiySde+65tXpZpMSn/nev7v75Qb39kjV64/mr1d7apLbmJrW1NOnCNSv0uo09UZcIRKZmQe3uWyVtlaTBwUGv1esi+X6w+5ju/vlBve/qAX3ilktlZlGXBDQUmn6IVLYwo0/fv1eXru/WHe+8hJAGyiCoEal7fnlIR0andMfNF6uthdURKKeS3fO+KekXkjab2WEz+0D4ZSENZmddd/3seQ2e16trLuyLuhygYVWy18ft9SgE6fPzZ0/ohRMZ/e07LqLlASyC75qIzDd/9Vut6mzTTa9dF3UpQEMjqBGJ/MysHtk3rBsvW6v2luaoywEaGkGNSDz+wojGswW97aI1UZcCNDyCGpHYvn9YLU2mN1+4OupSgIZHUCMS2/cP68rzetXd0Rp1KUDDI6hRd0Pj09rz4piu3dwfdSlALBDUqLtH9h+XJL3tIoIaqARBjbrbvn9Y/V3tunQ955cGKkFQo65mZl0/eWZYb31NPwe5ABUiqFFXTxwe1WgmT38aqAJBjbravm9YTSbO7QFUgaBGXT28f1iv39Sj3s62qEsBYoOgRt2cmMjqycOj7O0BVImgRt08vG9Y7tL1F6+NuhQgVghq1M2Pnx7Smq52XbaB3fKAahDUqItcYVaP7B/WdRevUVMTu+UB1SCoURc/O3Bc49mCrr+EtgdQLYIadfG9Xx9Rz/JWNiQCZ4CgRujGp/N6cM8x3fK69VzAFjgD/K9B6L6787CyhVm95w2boi4FiCWCGqGamXXd9bODesN5vbp8U0/U5QCxRFAjVN/ZeUi/PZnRX77lgqhLAWKLoEZohsez+uwP9mnwvF7deBl7ewBniqBGKE5O5vRXX9uhyVxB/3jrazmlKXAWWqIuoNQt//ITTednX77v7q+a7gveOe3u4s+V5Kc93xefPv8Fwn6/JZ4/f46ln19lfUvMv9Tyn8rPqNlMX7j9cl28jiMRgbPRUEF9Yf8K5Wfm/Ze3he/OH6XNH7PNH8SdPn3x55/+3vPmX/L1q3z+EgXU/P1Oe/7io95qXr+1uUm3XnGONq/rWvQ1ASytoYL687ddEXUJANBw6FEDQIMjqAGgwRHUANDgCGoAaHAENQA0OIIaABocQQ0ADY6gBoAGZ/MPFa7Ji5oNS3rhDJ/eJ+l4DcupFeqqDnVVh7qqk8S6znP3spdACiWoz4aZ7XD3wajrmI+6qkNd1aGu6qStLlofANDgCGoAaHCNGNRboy5gAdRVHeqqDnVVJ1V1NVyPGgDwao04ogYAlCCoAaDBRRLUZvbHZrbHzGbNbHDetDvM7ICZ7TOzGxd4/vlm9lgw37fMrC2EGr9lZruCn4NmtmuB+Q6a2W+C+XbUuo4y7/dJMztSUtvNC8x3U7AMD5jZx+pQ1z+Z2dNm9qSZ3WtmPQvMV5fltdTvb2btwWd8IFiXBsKqpeQ9N5nZQ2a2N1j/P1xmnmvN7FTJ5/uJsOsK3nfRz8WKvhAsryfN7Mo61LS5ZDnsMrMxM/vIvHnqsrzM7C4zGzKz3SWPrTKzbWb2TPBv7wLPfW8wzzNm9t4zKsDd6/4j6RJJmyU9LGmw5PFLJT0hqV3S+ZKeldRc5vnflnRbcPtOSX8dcr3/LOkTC0w7KKmvjsvuk5L+fol5moNld4GktmCZXhpyXTdIagluf0bSZ6JaXpX8/pL+RtKdwe3bJH2rDp/deklXBre7JO0vU9e1ku6v1/pU6eci6WZJD6h4Bbc3SXqszvU1Szqm4kEhdV9ekt4q6UpJu0se+6ykjwW3P1ZunZe0StJzwb+9we3eat8/khG1uz/l7vvKTHq3pHvcPevuz0s6IOmq0hmseGG/6yR9J3joq5L+MKxag/f7E0nfDOs9QnCVpAPu/py75yTdo+KyDY27P+juheDuo5I2hvl+S6jk93+3iuuOVFyXrreQL5Xu7kfd/fHg9rikpySdE+Z71tC7Jf2HFz0qqcfM1tfx/a+X9Ky7n+kRz2fF3R+RdHLew6Xr0EI5dKOkbe5+0t1HJG2TdFO1799oPepzJB0quX9Yp6/IqyWNloRCuXlq6S2SXnL3ZxaY7pIeNLOdZrYlxDpKfSj4+nnXAl+3KlmOYXq/iqOvcuqxvCr5/V+eJ1iXTqm4btVF0Gq5QtJjZSb/npk9YWYPmNlldSppqc8l6nXqNi08WIpieUnSWnc/Gtw+JmltmXlqstxCu7itmf1I0roykz7u7v8T1vtWo8Iab9fio+lr3P2Ima2RtM3Mng7++oZSl6QvSvq0iv+xPq1iW+b9Z/N+tahrbnmZ2cclFSR9Y4GXqfnyihszWyHpu5I+4u5j8yY/ruLX+4lg+8N/S3pNHcpq2M8l2Ab1Lkl3lJkc1fJ6FXd3MwttX+fQgtrd334GTzsiaVPJ/Y3BY6VOqPi1qyUYCZWbpyY1mlmLpD+S9IZFXuNI8O+Qmd2r4tfus1rBK112Zvbvku4vM6mS5VjzuszsLyTdIul6Dxp0ZV6j5surjEp+/7l5Dgef80oV161QmVmriiH9DXf/3vzppcHt7t83s38zsz53D/UERBV8LqGsUxV6p6TH3f2l+ROiWl6Bl8xsvbsfDdpAQ2XmOaJiH33ORhW3zVWl0Vof90m6Ldgif76Kfxl/WTpDEAAPSXpP8NB7JYU1Qn+7pKfd/XC5iWbWaWZdc7dV3KC2u9y8tTKvL3jrAu/3K0mvseLeMW0qfm28L+S6bpL0UUnvcvfMAvPUa3lV8vvfp+K6IxXXpR8v9MelVoIe+FckPeXun1tgnnVzvXIzu0rF/6Oh/gGp8HO5T9KfB3t/vEnSqZKv/WFb8FttFMurROk6tFAO/VDSDWbWG7Qpbwgeq07YW0sX2IJ6q4q9mqyklyT9sGTax1XcYr9P0jtLHv++pA3B7QtUDPADkv5LUntIdd4t6YPzHtsg6fsldTwR/OxRsQUQ9rL7mqTfSHoyWFHWz68ruH+zinsVPFunug6o2IvbFfzcOb+uei6vcr+/pH9Q8Q+JJHUE686BYF26oA7L6BoVW1ZPliynmyV9cG49k/ShYNk8oeJG2TfXoa6yn8u8ukzSvwbL8zcq2Vsr5No6VQzelSWP1X15qfiH4qikfJBdH1Bxm8b/SXpG0o8krQrmHZT05ZLnvj9Yzw5Iet+ZvD+HkANAg2u01gcAYB6CGgAaHEENAA2OoAaABkdQA0CDI6gBoMER1ADQ4P4fV8ze2t5Q8SMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3Bc53nf8e+DvQAgcSEJgOIVJGVBomlXtmRUcmtLli1faNcRHSeOqfRix05Yt1WTjCdN5LijeuSZTO1M3LRTpY58mSaZOLLixCnr0lEcW/KttkpYkSVRFEmY4v22AEncsbg9/eMcQCtoQSzIs4vFi99nBsPdPQd7Hh4sf3zxnPecY+6OiIgsfTWLXYCIiCRDgS4iEggFuohIIBToIiKBUKCLiARCgS4iEoiSAt3MdprZITPrNrP7iyz/L2b2dPx12MwuJ1+qiIhcic03D93MUsBh4B3AKWA/cK+7Pz/H+v8euMXdP5JwrSIicgXpEta5Deh296MAZvYIsAsoGujAvcB/mu9NW1tbfevWrSWWKSIiAD/5yU963L2t2LJSAn0jcLLg+Sng9mIrmtkWYBvwnfnedOvWrXR1dZWweRERmWZmx+dalvRB0d3A19x9co5C9phZl5l15XK5hDctIrK8lRLop4HNBc83xa8Vsxv4i7neyN0fdvdOd+9sayv6G4OIiFylUgJ9P9BhZtvMLEsU2ntnr2Rm24HVwI+SLVFEREoxb6C7+wRwH/AYcBB41N0PmNmDZnZPwaq7gUdcl28UEVkUpRwUxd33AftmvfbArOefSq4sERFZKJ0pKiISCAW6iEggFOgStLGJKb7y5An6hscXuxSRslOgS9D+5unT/O7Xn+UPv314sUsRKTsFugTtpyej68R1Xxhc5EpEyk+BLkE7Egf5iYvDi1yJSPkp0CVoJ+Mg7x0cW+RKRMpPgS7BcveZIB/MTzA6XvQSQyLBUKBLsPpHJhibnOLG6xoA6B3SKF3CpkCXYOUG8wBsbVkJwFB+YjHLESk7BboEq2c60FujQB8YVaBL2BToEqzpQN/SsgLQCF3Cp0CXYPUMqOUiy4sCXYLVNxIF+MZV9UA000UkZAp0CdbA6Dgrsima6zOAAl3Cp0CXYA2MTtBQm2ZlbXTZf7VcJHQKdAnWYH6Cxro02XQNqRpjRCcWSeAU6BKs/tFxGuqidktduobR8alFrkikvBToEqzB/ARNdVG7pS6T0qn/EjwFugRrYDRqucB0oGuELmFToEuwBuODogB1mRpGJzRCl7CVFOhmttPMDplZt5ndP8c6v2Rmz5vZATP7SrJliizcwOg4jdM99EyKvFouErj0fCuYWQp4CHgHcArYb2Z73f35gnU6gE8Ab3L3S2a2tlwFi5RicsoZGpssGKGr5SLhK2WEfhvQ7e5H3X0MeATYNWudXwMecvdLAO5+IdkyRRZm+iSil3roNTooKsErJdA3AicLnp+KXyt0I3Cjmf3QzH5sZjuLvZGZ7TGzLjPryuVyV1exSAleEejplHroErykDoqmgQ7gLuBe4Atmtmr2Su7+sLt3untnW1tbQpsWeaWB0XGAl/XQ1XKR0JUS6KeBzQXPN8WvFToF7HX3cXd/EThMFPAii2L62ufTPfRatVxkGSgl0PcDHWa2zcyywG5g76x1/oZodI6ZtRK1YI4mWKfIggyOzu6h68QiCd+8ge7uE8B9wGPAQeBRdz9gZg+a2T3xao8BvWb2PPA48B/cvbdcRYvMp392yyWtlouEb95piwDuvg/YN+u1BwoeO/Dx+Etk0WmWiyxHOlNUgjRQpOUyMeVMTvliliVSVgp0CdLg6ASpGqM+kwIgm44+6mMTartIuBToEqSB0XEaatOYGQDZlAJdwqdAlyAN5F+6MBdAJh6h5yfVR5dwKdAlSIWXzgWojUfo45PqoUu4FOgSpMFZga4euiwHCnQJ0kD+pUvnggJdlgcFugRpdstFB0VlOVCgS5AK71YELx0UHdNBUQmYAl2CFI3QC1ou8Qg9rxG6BEyBLsHJT0wyNjlV9KCoZrlIyBToEpzZp/0D1OqgqCwDCnQJzuCsa6EDZHRQVJYBBboE56URepFpizooKgFToEtwBvLT10LXiUWyvCjQJTizbz8HBfPQdVBUAqZAl+BM99CbdKaoLDMKdAnOQHz7uQadKSrLjAJdglO05aIRuiwDCnQJzmB+gtp0zUyIA6RqjFSNaZaLBK2kQDeznWZ2yMy6zez+Iss/bGY5M3s6/vrV5EsVKU3/rNP+p2VTNTpTVIKWnm8FM0sBDwHvAE4B+81sr7s/P2vVr7r7fWWoUWRBBkbHaap75Uc7kzK1XCRopYzQbwO63f2ou48BjwC7yluWyNUbzE+87IDotGw6pYtzSdBKCfSNwMmC56fi12b7BTN7xsy+ZmabE6lO5CrMvhb6tNp0jUboErSkDor+b2Cru98MfAv4k2IrmdkeM+sys65cLpfQpkVebmB0nMbaIj30dA1jkwp0CVcpgX4aKBxxb4pfm+Huve6ej59+EXhDsTdy94fdvdPdO9va2q6mXpF5DY7O0XJJ1TCuEboErJRA3w90mNk2M8sCu4G9hSuY2fqCp/cAB5MrUWRh5mq5ZNKmEboEbd5ZLu4+YWb3AY8BKeDL7n7AzB4Eutx9L/DrZnYPMAFcBD5cxppF5jQ15QyOzT1tUT10Cdm8gQ7g7vuAfbNee6Dg8SeATyRbmsjCDY1N4A6NtcVmuSjQJWw6U1SCUuxuRdOy6ZRaLhI0BboEpdjNLaZldWKRBE6BLkEZzL/ySovTsukaxjVCl4Ap0CUo/VdouWRSmocuYVOgS1AGZm5uoXnosvwo0CUogzPXQn9lDz2jM0UlcAp0CUrfSNRDb6ovPkLXQVEJmQJdgnJ5ZIxsuob6TOoVy3QtFwmdAl2C0jc8zqr6DGb2imWZlOkGFxI0BboE5fLwOKtWvLJ/DpBNpZiccianFOoSJgW6BOXS8Bir6rNFl2XS0ahdc9ElVAp0CUrfyJVG6NHHXXctklAp0CUoV2y5pKOPu0boEioFugTl8sgYq1YUb7lMj9A1dVFCpUCXYIyOTzI6PkVzffEReialEbqETYEuwbg8HJ1UNF/LRSN0CZUCXYJxcWgMgNVztFymR+g6uUhCpUCXYPQMRvcpb2usLbq8ViN0CZwCXYJxYSAK9LVzBPpLPXSdWCRhUqBLMHJxoLc2FA909dAldCUFupntNLNDZtZtZvdfYb1fMDM3s87kShQpTW4gz8psipVFbhAN0bVcQLNcJFzzBrqZpYCHgHcDO4B7zWxHkfUagd8Anky6SJFS5AbzrG2qm3P59AhdZ4pKqEoZod8GdLv7UXcfAx4BdhVZ79PAZ4DRBOsTKdn5vtE5D4jCSycWaYQuoSol0DcCJwuen4pfm2FmtwKb3f3/JFibSFHuzgvn+umNZ7VMO3FxmM2rV8z5feqhS+iKNxsXwMxqgM8BHy5h3T3AHoD29vZr3bQsU7+37yBf+P6L1Bj8ypu28bvveTXjk1Oc6x9lS8vcga4zRSV0pQT6aWBzwfNN8WvTGoHXAk/ENxVYB+w1s3vcvavwjdz9YeBhgM7OTs0dkwX7WW6QL/7gRd5783oa6zJ86QcvMpSf4KNv3gbA5jX1c37vzAhdgS6BKiXQ9wMdZraNKMh3A788vdDd+4DW6edm9gTwW7PDXCQJf/3UKQx44Od2sLaxjtUrMvzREz/jfH906Gb7uqY5vzeji3NJ4Obtobv7BHAf8BhwEHjU3Q+Y2YNmdk+5CxQp9M3nzvGmG1pZ2xjNZvn4O26kc8tqHj+UA6BjbcOc31urEboErqQeurvvA/bNeu2BOda969rLEnml3ECeo7khPtj5UgcwnarhD3e/nvc99EN+/paNpFNzj1FmeugT6vZJmK75oKhIpew/dhGA27atednrm1av4P/ef/dMj3wuqRojVWOMTU6WrUaRxaRT/2XJ2H/sIvWZFK/d2PyKZfOF+bRMynQtFwmWAl2WjINn+9m+vnGmdXI1sqkaHRSVYCnQZUmITiYaYPu6xmt6n2y6RgdFJVgKdFkSLgzkuTw8zk3XXWOga4QuAVOgy5LwwrkBALavn3ueeSky6RqdKSrBUqDLknDoXD/AtbdcNEKXgCnQZUl4sWeYNSuzrJrjfqGlyqQ0QpdwKdBlSTjeO0T7mrkvvFWqbLpG10OXYCnQZUk43jt8xSspliqrEboETIEuVW9sYoqzfSNsSWiErh66hEqBLlXv1KVhphzaW1Ze83vpTFEJmQJdqt7xi8MAybRcNEKXgCnQpeqd6I0DPYGWi2a5SMgU6FL1Tlwcpj6TuuINoEulWS4SMgW6VL0zl0fYuLqe+BaH10SzXCRkCnSpemf6RlnfXJfIe+niXBIyBbpUvXN9I6xrSibQM6kaxtVykUAp0KWqjU9OcWEgz/pV9Ym8n0boEjIFulS1CwN53Ems5RLNcnHcNRddwqNAl6p2rm8EgHUJBXptfKs6jdIlRCUFupntNLNDZtZtZvcXWf4xM3vWzJ42sx+Y2Y7kS5Xl6GzfKAAbmpNpuWRS0UwZnS0qIZo30M0sBTwEvBvYAdxbJLC/4u7/yN1fD3wW+FzilcqydPZyFOhJjdCz8f1IdbaohKiUEfptQLe7H3X3MeARYFfhCu7eX/B0JaDhjyTibN8oK7IpmurSibxfJm65aC66hKiUfyUbgZMFz08Bt89eycz+HfBxIAu8LZHqZNk71z/C+ua6RE4qAo3QJWyJHRR194fc/VXA7wD/sdg6ZrbHzLrMrCuXyyW1aQnYmcujrE+ofw7RtEXQQVEJUymBfhrYXPB8U/zaXB4B3ldsgbs/7O6d7t7Z1tZWepWybJ3rG02sfw4aoUvYSgn0/UCHmW0zsyywG9hbuIKZdRQ8/WfAkeRKlOVqasrJDea5runaL8o1LZNSD13CNW8P3d0nzOw+4DEgBXzZ3Q+Y2YNAl7vvBe4zs7cD48Al4EPlLFqWh0vDY0xOOW0NyQX6TMtFI3QJUElTB9x9H7Bv1msPFDz+jYTrEqFncAyAtsbkWi7TI3T10CVEOlNUqlZuIA9Aa0M2sffUCF1CpkCXqpUbjE4qSuLGFtOyMz10nSoh4VGgS9WaHqEnGug6sUgCpkCXqpUbyFOXqaGhNpmzROGla7mo5SIhUqBL1coN5GlrrE3sLFHQiUUSNgW6VK3cYD7RKYugE4skbAp0qVo9A2O0Jh3o6qFLwBToUrVyg/lED4hCwTx0jdAlQAp0qUrjk1NcHBpLPNA1QpeQKdClKvXOnCWabKCnazTLRcKlQJeqNDMHPeEeupmRTdcwphOLJEAKdKlK02eJtiY8QodopotG6BIiBbpUpXKN0CHqo6uHLiFSoEtVKsdp/9MyKdMIXYKkQJeq1DM4RmNtmrpMKvH31ghdQqVAl6o0fdp/OWRSNeQV6BIgBbpUpdxgPvGzRKdlUzWMq+UiAVKgS1XqGczT2pjcjS0KRdMWFegSHgW6VKWegeQvzDUtm1IPXcKkQJeqMzo+Sf/oRNlaLhnNQ5dAlRToZrbTzA6ZWbeZ3V9k+cfN7Hkze8bMvm1mW5IvVZaL3qHotP9ynFQE6ExRCda8gW5mKeAh4N3ADuBeM9sxa7V/ADrd/Wbga8Bnky5Ulo+eMp5UBBqhS7hKGaHfBnS7+1F3HwMeAXYVruDuj7v7cPz0x8CmZMuU5aRnMAr0co3QazUPXQJVSqBvBE4WPD8VvzaXjwLfvJaiZHmbPku0taE8s1x0pqiEKrm77wJm9i+ATuAtcyzfA+wBaG9vT3LTEpCZEXq5ZrlohC6BKmWEfhrYXPB8U/zay5jZ24FPAve4e77YG7n7w+7e6e6dbW1tV1OvLAM9g2M01pXntH9QD13CVUqg7wc6zGybmWWB3cDewhXM7Bbgj4nC/ELyZcpykivjHHTQiUUSrnkD3d0ngPuAx4CDwKPufsDMHjSze+LVfh9oAP7SzJ42s71zvJ3IvM72jbCuua5s76/roUuoSuqhu/s+YN+s1x4oePz2hOuSZexc3yhvfFVL2d5fPXQJlc4UlaoyOeWcH8izvowj9EyqhimHCYW6BEaBLlWlZzDP5JSzvrm+bNvIpqOP/bjOFpXAKNClqpy5PAJQ9hE6oD66BEeBLlXlXF90c+iyHhSNR+ia6SKhUaBLVTkbB3pZWy4pAxToEh4FulSVc/2jZNM1rF6RKds2ZkboarlIYBToUlXOXB5hfXMdZla2bdSlozNQ8xOTZduGyGJQoEtVOX15hA1lbLcA1GWjQB8ZU6BLWBToUlVO9A6zpWVFWbdRH18jZmRcgS5hUaBL1RjMT9A7NEZ7hQJ9VIEugVGgS9U4eTG6R0r7mjIH+kzLRQdFJSwKdKkax3ujQN+yZmVZtzM9Qh8emyjrdkQqTYEuVaNSI/Q6tVwkUAp0qRrHLw7RVJemuYxz0KGg5aJAl8Ao0KVqHOsZZmtredstAHXxiUXqoUtoFOhSNQ6fH6BjbWPZt5NO1ZBN1WiELsFRoEtV6Bse58JAno7rGiqyvfpsSj10CY4CXarC4QsDANxYqUDPpHSmqARHgS5V4cj5QYCKtFwgGqGr5SKhUaBLVTh0rp8V2RQbV5X3Oi7T6jIKdAmPAl2qwjOn+3jthmZqasp3lcVC9Zka9dAlOCUFupntNLNDZtZtZvcXWX6nmT1lZhNm9ovJlykhG5uY4sCZfl63ubli26zPqocu4Zk30M0sBTwEvBvYAdxrZjtmrXYC+DDwlaQLlPC9cK6fsYkpXr95dcW2Wa+WiwQoXcI6twHd7n4UwMweAXYBz0+v4O7H4mU6U0MWrOvYJYCKjtDVQ5cQldJy2QicLHh+Kn5twcxsj5l1mVlXLpe7mreQAP2gu4etLSvYtLq813AppGmLEqKKHhR194fdvdPdO9va2iq5aalS+YlJfny0lzs6Kvt5WFmbZiivqy1KWEoJ9NPA5oLnm+LXRK7Z4y/kGB6b5G2vXlvR7TbWpRnMT+DuFd2uSDmVEuj7gQ4z22ZmWWA3sLe8Zcly8WjXSdoaa7njhtaKbrehNs2Uw7DaLhKQeQPd3SeA+4DHgIPAo+5+wMweNLN7AMzsH5vZKeADwB+b2YFyFi1Ln7vzt8+d5TsvXOBfvnEL6VRlT4lorIsu0TswqraLhKOUWS64+z5g36zXHih4vJ+oFSMyp/7RcX5wpIcnDl3giUM5LgzkeVXbSj765m0Vr6WxLvroD+bHgbqKb1+kHEoKdJGr4e48f7af7x7O8cShHD85fonJKaexLs0dHa3cddNa3vWadaysrfzHsCEO9H6N0CUgCnRJ1KWhMb7f3cN3D+X43pEcuYE8ADvWN/Gv77yeu25ay63tqyreYpmtKQ50tVwkJAp0uSaTU87TJy/z3cM5vns4xzOnLuMOq1ZkuKOjjbfc2MadHa2sbaqutkZDbdRDH1SgS0AU6LJg5/tHZwL8B0d66BsZp8bg9ZtX8Zt338idN7Zy86ZVpCp0oa2r0TgzQh9f5EpEkqNAl3lFo/BLfOeFCzz+Qo7nz/YDsLaxlnfuuI633NTGm29oZdWK7CJXWrqGmYOiGqFLOBToUtTFoTG+ezgK8O8dyXF5eJxUjfGGLav5nZ3bueumNrava8SsekfhV9KQ1UFRCY8CXQCYmnIOnOnn8UMXePzQBZ4+GfXCWxuy3L39Ot62fS1v7miluT6z2KUmoqbGaKxN0z+ilouEQ4G+jE3PC3/8hQs8cTiakWIGN2+KeuFv3d5W0ZtOVNqahiyXhscWuwyRxCjQlxF359D5Ab57KMfjhy7QdewSE1NOU12at9y0lrfe1MadN7bR2lC72KVWRMvKLL2DCnQJhwI9cLmBPD/s7uF7R3J8/0jPzLzwV69vYs+d1/PW7Wu5ZfPizwtfDC0NtZy8OLzYZYgkRoEemNHxSbqOXeL7R3J870gPB+MZKWtWZnnTDa3c0RF9rW+uzM2Yq1lrQ5anT15e7DJEEqNAX+ImJqd47kw/Tx7t5Yc/6+XJo73kJ6bIpIzOLWv47Z03cWdHGzvWNwXbC79aLStruTg0xtSUa99IEBToS8zE5BTPnu7jyRcv8uOjvXQduzQzl/qGtQ388u3t3NnRxu3Xr2FFVj/eK2lpyDI55fSNjLN65dKZQy8yF/2Lr3IjY5M8c+oyPzlxiSePXqTr2EWG4mt4d6xt4H23bOCN17dw+7YW2hqXx8HMpLTEB39zg3kFugRBgV5F3J1Tl0Z46sQlnjp+iadOXObg2X4mpqK76nSsbeD9t27ijde3cNu2NQrwa7RxVXQc4fSlEW68rnGRqxG5dgr0RZQbyPPcmT4OnO7jmVN9PHXiMj2D0SyU+kyK121uZs+d13Nr+2puaV81M6KUZLSviW5Kfbx3aJErEUmGAr0C3J0zfaM8d7qPA2f6OXC6j+fO9HG+Pz+zztaWFdzR0cqt7au4pX0129c1LsuphJXU2pClPpPixMWRxS5FJBEK9AS5O2f7Rjl8foDuC4McPj/AkQuDdJ8fZCA+cFlj0cHLN72qlddsbOa1G5p49YYmmurCOKV+KTEz2tes4MRFjdAlDAr0BXJ3Lg6NcfziMMd7hzjeO8zx3mFe7Bmi+8Lgy67e19qQ5Ya1Dfz8rRvpuK6R12xo4tXrmqjPphbxbyCFOq5r4Knjlxa7DJFEKNBnmZpyeofGONs3wpnLo5zrG+Fs3yinLo1wrHeIE73DM6NtADNY31THlpaVvD8O7hvXNtBxXSNrNHOi6r1+8yq+8cxZLvSPVt1NOEQWqqRAN7OdwH8FUsAX3f0/z1peC/wp8AagF/igux9LttRrMzo+ycWhMXoG8/QOxn8OjdE7mCc3kOdM3yhn+0Y435dnbHLqZd+bTdewcVU97WtW0LllNe0tK9nasoItLSvYtHoFdRmNuJeqW7esBuBHR3vZ9fqNi1yNyLWZN9DNLAU8BLwDOAXsN7O97v58wWofBS65+w1mthv4DPDBchQ8lJ8gN5Cnf3ScvpFx+kcm6B8dp38kfj760mt9I+NcGhqjd3DsZaPqQnWZGlobatnQXM+t7atZ11zHhuZ61jfXsb65nvWr6mhZmV2y1/2WK3vdplWsa6rj6/9wWoEuS14pI/TbgG53PwpgZo8Au4DCQN8FfCp+/DXgv5uZubsnWCsAf/KjY3z2bw8VXZaqMZrrMzTVpWmqz9BUl2HT6hW0NmRpbailZWWWloZaWhqytK6M/lyMO85L9UjVGP/89nb+4FuH+er+E3zgDZt1GQBZskpJs43AyYLnp4Db51rH3SfMrA9oAXqSKLLQ27avZV1THU11GZpXRKHdVJ+mqS7DimxKI2lZsF+783q+393D7/zVs3xq7/OsWZmlNl2TeLDrkynTfv3uDn7udRsSf9+KDk/NbA+wB6C9vf2q3mP7uia2r2tKsixZ5uoyKf78V2/nsQPneOr4ZfpGxslPTJLk75dO4r+syhJWrjt/lRLop4HNBc83xa8VW+eUmaWBZqKDoy/j7g8DDwN0dnbqEy5VI5Oq4b03b+C9Nyc/ahKplFJORdwPdJjZNjPLAruBvbPW2Qt8KH78i8B3ytE/FxGRuc07Qo974vcBjxFNW/yyux8wsweBLnffC3wJ+DMz6wYuEoW+iIhUUEk9dHffB+yb9doDBY9HgQ8kW5qIiCyErv4kIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIW6zp4maWA45f5be3UobLCiRAdS2M6lq4aq1NdS3MtdS1xd3bii1YtEC/FmbW5e6di13HbKprYVTXwlVrbaprYcpVl1ouIiKBUKCLiARiqQb6w4tdwBxU18KoroWr1tpU18KUpa4l2UMXEZFXWqojdBERmaVqA93MPmBmB8xsysw6Zy37hJl1m9khM3vXHN+/zcyejNf7anzp36Rr/KqZPR1/HTOzp+dY75iZPRuv15V0HUW29ykzO11Q23vmWG9nvA+7zez+CtT1+2b2gpk9Y2ZfN7NVc6xXkf0139/fzGrjn3F3/FnaWq5aCra52cweN7Pn48//bxRZ5y4z6yv4+T5Q7L3KUNsVfy4W+W/x/nrGzG6tQE03FeyHp82s38x+c9Y6FdtfZvZlM7tgZs8VvLbGzL5lZkfiP1fP8b0fitc5YmYfKrbOvNy9Kr+AVwM3AU8AnQWv7wB+CtQC24CfAaki3/8osDt+/Hng35S53j8AHphj2TGgtYL77lPAb82zTired9cD2Xif7ihzXe8E0vHjzwCfWaz9VcrfH/i3wOfjx7uBr1bgZ7ceuDV+3AgcLlLXXcA3KvV5KvXnArwH+CbR3fbeCDxZ4fpSwDmiedqLsr+AO4FbgecKXvsscH/8+P5in3tgDXA0/nN1/Hj1QrdftSN0dz/o7sXuBr0LeMTd8+7+ItBNdCPrGRbdWPRtRDesBvgT4H3lqjXe3i8Bf1GubZTBzM2/3X0MmL75d9m4+9+5+0T89MdEd79aLKX8/XcRfXYg+izdbWW+aa27n3X3p+LHA8BBonv2LgW7gD/1yI+BVWa2voLbvxv4mbtf7QmL18zdv0d0T4hChZ+jubLoXcC33P2iu18CvgXsXOj2qzbQr6DYTatnf+BbgMsF4VFsnSTdAZx39yNzLHfg78zsJ/F9VSvhvvjX3i/P8SteKfuxnD5CNJorphL7q5S//8tufg5M3/y8IuIWzy3Ak0UW/xMz+6mZfdPMXlOhkub7uSz2Z2o3cw+qFmN/TbvO3c/Gj88B1xVZJ5F9V9GbRM9mZn8PrCuy6JPu/r8qXU8xJdZ4L1cenb/Z3U+b2VrgW2b2Qvw/eVnqAv4H8Gmif4CfJmoHfeRatpdEXdP7y8w+CUwAfz7H2yS+v5YaM2sA/gr4TXfvn7X4KaK2wmB8fORvgI4KlFW1P5f4GNk9wCeKLF6s/fUK7u5mVraphYsa6O7+9qv4tlJuWt1L9OteOh5ZFVsnkRotuin2+4E3XOE9Tsd/XjCzrxP9un9N/xBK3Xdm9gXgG0UWlbIfE6/LzD4MvBe42+PmYZH3SHx/FZHYzc+TZmYZojD/c3f/69nLCwPe3feZ2R+ZWau7l/WaJSX8XMrymbWipe0AAAHOSURBVCrRu4Gn3P387AWLtb8KnDez9e5+Nm5BXSiyzmmiXv+0TUTHDxdkKbZc9gK74xkI24j+p/1/hSvEQfE40Q2rIbqBdblG/G8HXnD3U8UWmtlKM2ucfkx0YPC5YusmZVbf8ufn2F4pN/9Ouq6dwG8D97j78BzrVGp/VeXNz+Me/ZeAg+7+uTnWWTfdyzez24j+HZf1P5oSfy57gX8Vz3Z5I9BX0Gootzl/S16M/TVL4edorix6DHinma2OW6TvjF9bmEoc+b2aL6IgOgXkgfPAYwXLPkk0Q+EQ8O6C1/cBG+LH1xMFfTfwl0Btmer8n8DHZr22AdhXUMdP468DRK2Hcu+7PwOeBZ6JP0zrZ9cVP38P0SyKn1Worm6iPuHT8dfnZ9dVyf1V7O8PPEj0Hw5AXfzZ6Y4/S9dXYB+9mahV9kzBfnoP8LHpzxlwX7xvfkp0cPmfVqCuoj+XWXUZ8FC8P5+lYHZamWtbSRTQzQWvLcr+IvpP5SwwHufXR4mOu3wbOAL8PbAmXrcT+GLB934k/qx1A79yNdvXmaIiIoFYii0XEREpQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigfj/xn+a9cSvNMMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-10, 10, .0001)\n",
    "plt.plot(x, f(x))\n",
    "plt.show()\n",
    "plt.plot(x, fp(x))\n",
    "plt.show()\n",
    "plt.plot(x, norm.pdf(f(x))*fp(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Invert Warping Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert(f, z):\n",
    "    result = root_scalar(lambda x: f(x) - z, x0=z, x1=z+1)\n",
    "    return result.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Moments of Base Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_moment(f, p=1):\n",
    "    zx, wx = hermgauss(50)\n",
    "    c = np.sqrt(2)\n",
    "    f_inv = lambda z: invert(f, c*z)**p\n",
    "    result = sum([f_inv(z)*w for z, w in zip(zx, wx)])\n",
    "    result *= c / np.sqrt(2*np.pi)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_stddev(f):\n",
    "    mean = compute_moment(f)\n",
    "    moment2 = compute_moment(f, 2)\n",
    "    return moment2 - mean**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 \t -2.226769149373234 \t 3.90479258921163\n"
     ]
    }
   ],
   "source": [
    "f_volume = compute_moment(f, 0)\n",
    "f_mean = compute_moment(f, 1)\n",
    "f_stddev = compute_stddev(f)\n",
    "print(f_volume, \"\\t\", f_mean, \"\\t\", f_stddev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f33d3c93828>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3Qc53nf8e+zd+LCC0hQpMSrJIoSbSu2DCt2HNuqLduS04i5OaF6mtqxUx2nYS4nTVL5uEf1kf+pnZO0Tas2UWKfOGkSWXabhImZKI5j+ZLaCimZkk1SJCFKFEGRxIIggAXIXWCxT/+YWXC1XAALYmZx0e9zDg52Z17svBgsfnjwzjsz5u6IiMjSl1joDoiISDQU6CIiy4QCXURkmVCgi4gsEwp0EZFlQoEuIrJMNBXoZnaPmR0zs14ze7DB+v9iZofCj+NmNhR9V0VEZCY22zx0M0sCx4H3An3AAeB+dz8yTftfAt7k7h+JuK8iIjKDVBNt7gR63f0kgJk9BuwGGgY6cD/wn2Z70XXr1vm2bdua7KaIiAA8/fTTA+7e3WhdM4F+A3C65nkf8IONGprZVmA78I+zvei2bds4ePBgE5sXEZEqMzs13bqoD4ruAb7k7pPTdOQBMztoZgfz+XzEmxYReW1rJtDPAJtrnm8KlzWyB/jz6V7I3R919x537+nubvgfg4iIXKNmAv0AsMPMtptZhiC099U3MrNbgTXAt6PtooiINGPWQHf3MrAXeAI4Cjzu7ofN7GEzu6+m6R7gMdflG0VEFkQzB0Vx9/3A/rplD9U9/2R03RIRkbnSmaIiIsuEAl1EZJlQoIu0wFipzJ8+dYriRMMZvSKRaGoMXUTm5/e/cZLf/eoJKhXnZ9+2baG7I8uUKnSRFjh8Zjj4/MrIAvdEljMFukgLHDtfAODFgbEF7oksZwp0kZhNVpxXhi4DcGFsfIF7I8uZAl0kZhcvjVMJT7cbVKBLjBToIjHLF0oA7FjfwcVL40xWdDK1xEOBLhKzgdEg0G/sbscdxsbLC9wjWa4U6CIxqwb6trXtQDAnXSQOCnSRmFWHXLYq0CVmCnSRmF0YGyeTSrBhVRaA0ZLOFpV4KNBFYjZyuczKXJr2THBitip0iYsCXSRmo6UynbkU7dnU1HOROCjQRWJWKE7QmUvRkVWFLvFSoIvErFAs05FN0ZZNAgp0iY8CXSRmo8VgyCWXDgK9OFFZ4B7JcqVAF4lZoThBRzZNLhUEeqmsWS4SDwW6SMwK4UHRdNJImCp0iY8CXSRGlYozWiqzMpfCzMilk7prkcSmqUA3s3vM7JiZ9ZrZg9O0+WkzO2Jmh83sz6LtpsjSNDZexh06csEMl1w6SVFDLhKTWW9BZ2ZJ4BHgvUAfcMDM9rn7kZo2O4CPA29394tmtj6uDossJdU55525NAC5VEJDLhKbZir0O4Fedz/p7uPAY8Duujb/FnjE3S8CuHt/tN0UWZoKxSDQq3PQsxpykRg1E+g3AKdrnveFy2rdAtxiZv9kZt8xs3savZCZPWBmB83sYD6fv7Yeiywh1UDvDIdcsqrQJUZRHRRNATuAu4D7gT8ws9X1jdz9UXfvcfee7u7uiDYtsngVihPAlUDPpZOatiixaSbQzwCba55vCpfV6gP2ufuEu78IHCcIeJHXtCsVejiGnk5QUoUuMWkm0A8AO8xsu5llgD3Avro2f0lQnWNm6wiGYE5G2E+RJenKQVHNcpH4zRro7l4G9gJPAEeBx939sJk9bGb3hc2eAC6Y2RHga8BvuPuFuDotslRUh1yqB0VzKR0UlfjMOm0RwN33A/vrlj1U89iBXws/RCQ0WixjxtS10LNpHRSV+OhMUZEYjRTLdGRSJBIGqEKXeCnQRWJUvblFVS6dUKBLbBToIjEKbm6RnnqeTScZn9SQi8RDgS4So0KxPHUdF4B00piY9AXskSxnCnSRGNUPuWSSSSYrzmRFoS7RU6CLxKh6+7mqTCr4lRsva9hFoqdAF4lRoVh+1Rh6OhnMdtE4usRBgS4So+Cg6JUKPasKXWKkQBeJyXi5QqlcobPRkIsqdImBAl0kJvXXcQGNoUu8FOgiMZm6jsurxtCDX7kJVegSAwW6SEzqb24BkEmqQpf4KNBFYjIV6A3G0EsKdImBAl0kJlfuVnRlyCWjIReJkQJdJCbVg6IdOigqLaJAF4lJwzF0BbrESIEuEpNG0xY1y0XipEAXiclIcYJMMkE2lZxaphOLJE4KdJGY1F86F64cFNUsF4mDAl0kJqPFV186F3QtF4lXU4FuZveY2TEz6zWzBxus/7CZ5c3sUPjx89F3VWRpqb8wF2gMXeKVmq2BmSWBR4D3An3AATPb5+5H6pp+wd33xtBHkSWpUCzTmU2/aplmuUicmqnQ7wR63f2ku48DjwG74+2WyNJXaDDkokCXODUT6DcAp2ue94XL6v2kmT1nZl8ys82R9E5kCau/QTRAKhHc4EJDLhKHqA6K/jWwzd1vB74CfL5RIzN7wMwOmtnBfD4f0aZFFqdGFbqZkUklKCnQJQbNBPoZoLbi3hQum+LuF9y9FD79Q+DNjV7I3R919x537+nu7r6W/oosCZWKMzpeZmXu6sNU2WRCQy4Si2YC/QCww8y2m1kG2APsq21gZhtrnt4HHI2uiyJLz9h4GXeumocOkE4lNOQisZh1lou7l81sL/AEkAQ+5+6Hzexh4KC77wN+2czuA8rAIPDhGPsssuhduY5L+qp1GVXoEpNZAx3A3fcD++uWPVTz+OPAx6PtmsjS1ejCXFWZlAJd4qEzRUVi0Oha6FWZVELXcpFYKNBFYlBocKXFqnQywXjZW90leQ1QoIvEoDrk0miWiyp0iYsCXSQGMw25ZJMJJjSGLjFQoIvEoFqhd2QbTVs0VegSCwW6SAwKxQmSCaMtk7xqXTqpeegSDwW6SAwKxTId2RRmdtU6zUOXuCjQRWIwfHmC1W1Xj59DcKaohlwkDgp0kRgMXZpg9YrGgZ7VkIvERIEuEoOhS+Osbss0XJfWkIvERIEuEoOhGYZcMqkEE5M6sUiip0AXicHFsXHWqEKXFlOgi0RssuKMFMusmmYMXWeKSlwU6CIRG74cnCW6Zrohl6QxXq7grmEXiZYCXSRiQ5fGAWY8KApQrijQJVoKdJGIDYUV+kwHRQGNo0vkFOgiEWu2QtdcdImaAl0kYgOFINDXtjcO9KkKXYEuEVOgi0Ssv1AEoLsz23B9JqkhF4mHAl0kYv2FEqtWpMmlr77SIlyp0HVykUStqUA3s3vM7JiZ9ZrZgzO0+0kzczPria6LIktL/0hp2uocroyhq0KXqM0a6GaWBB4B7gV2Afeb2a4G7TqBXwGeirqTIktJf6HI+hkC/UqFrkCXaDVTod8J9Lr7SXcfBx4Ddjdo9yng00Axwv6JLDn9hdkq9OAa6SVV6BKxZgL9BuB0zfO+cNkUM7sD2OzuX46wbyKLWqXifP/M8NSZoRBU3WeHi2xe0zbt12U0bVFicvUND+fIzBLA7wAfbqLtA8ADAFu2bJnvpkUWjLvzS499ly8/d5ZMMsGvvncHv/Cumzhz8TKTFWfL2hkCXScWSUyaCfQzwOaa55vCZVWdwOuBJ8PbbW0A9pnZfe5+sPaF3P1R4FGAnp4eHeKXJevbL1zgy8+d5WffupULYyU+83fHmCg7b9yyGoCtXdMHuk4skrg0E+gHgB1mtp0gyPcA/6q60t2HgXXV52b2JPDr9WEuspx88ek+VuZSfOJHbiObSvDvv/gs//Wrx3nXLd0A3HJd57Rfq4OiEpdZx9DdvQzsBZ4AjgKPu/thM3vYzO6Lu4Mii81kxfmHI+e59/UbyaWTmBmf2v16tna18eSxPBtW5lgzzVmicKVC10FRiVpTY+juvh/YX7fsoWna3jX/boksXkdeGaFQKvNDN6+dWtaeTfHf77+DPY9+m59/x/YZvz6rE4skJvM+KCryWvOdkxcAeOuNa1+1/A2bVvHdh943NaQyHZ1YJHHRqf8ic/T0qYts6WrjupW5q9bNFua1bTSGLlFToIvM0fPnRnjd9Suv+eurJxapQpeoKdBF5mCsVObU4CVu3TCfQNflcyUeCnSROTh+voA73Lpx+mmJs9HlcyUuCnSROTh6tgDAro3XXqEnEkYqYRpDl8gp0EXm4Pj5Am2ZJDesXjGv18mkEgp0iZwCXWQOTl0YY+vadhIJm9frpJMJDblI5BToInNw6sIlts1w4a1mZVIJxnVikURMgS7SpMmKc/ripRmvpNisjCp0iYECXaRJrwxdZmLS2ba2fd6vpTF0iYMCXaRJLw9eAmBrBBV6Ommq0CVyCnSRJp26UA30+Vfo6aQqdImeAl2kSacujJFJJdjY4BoucxUcFFWgS7QU6CJNOn3xEptWr5j3lEXQtEWJhwJdpEmvDBW5fp4nFFVldVBUYqBAF2nSueEiG1bNf7gFwgpdgS4RU6CLNKE8WaG/UOT6iAI9k0wwUdaJRRItBbpIE84XSlQcNqyKZsglrYOiEgMFukgTzg1fBmDj6qiGXDQPXaKnQBdpwitDRQCuj6hCz6pClxg0Fehmdo+ZHTOzXjN7sMH6j5nZ98zskJl9y8x2Rd9VkYVzbjgI9CgPimqWi0Rt1kA3syTwCHAvsAu4v0Fg/5m7v8Hd3wh8BvidyHsqsoBeGb5MeybJylwqktfTxbkkDs1U6HcCve5+0t3HgceA3bUN3H2k5mk7oMP3sqycHSqycfUKzOZ/UhEEB0VVoUvUmik3bgBO1zzvA36wvpGZ/SLwa0AGeHckvRNZJM6OFNkY0XALhNMWJx13j+yPhEhkB0Xd/RF3vwn4D8B/bNTGzB4ws4NmdjCfz0e1aZHYnR26HG2gp8IbRatKlwg1E+hngM01zzeFy6bzGPBjjVa4+6Pu3uPuPd3d3c33UmQBTVacgdESGyK4KFdVJhn86k3orkUSoWYC/QCww8y2m1kG2APsq21gZjtqnv4IcCK6LoosrMGxcSoO3Z3ZyF4znQyGWXRgVKI06xi6u5fNbC/wBJAEPufuh83sYeCgu+8D9prZ3cAEcBH4UJydFmml/kIwZTHSQE9VK3QFukSnqTlY7r4f2F+37KGax78Scb9EFo18oQREG+jVIRdV6BIlnSkqMoupQO/QQVFZ3BToIrPIjwaBvq4zE9lrqkKXOCjQRWaRL5ToyKZoy0RzligEp/6DxtAlWgp0kVnkC6VIx8/hypCLAl2ipEAXmUW+UKK7I9pAr1boJQ25SIQU6CKzyI/GWaHrxCKJjgJdZBaxDLnooKjEQIEuMoPixCSFYjnyQE+ngjNFNYYuUVKgi8wgjpOKQBW6xEOBLjKD6hz0yCv0pE4skugp0EVmcOUs0WgDPZtShS7RU6CLzKAa6OtjqtA1hi5RUqCLzCBfKGEGXe3RnfYPOrFI4qFAF5lBfrREV1uGVDLaX5W0DopKDBToIjMYKJRYF/H4OdTc4EInFkmEFOgiM8iPliK9ymKVmZFJJlShS6QU6CIzGBiN/jouVemkaQxdIqVAF5nBQGE8liEXCA6MqkKXKCnQRaYxVipzeWKSdRFPWaxKJxOq0CVSCnSRaQxU71SkCl2WiKYC3czuMbNjZtZrZg82WP9rZnbEzJ4zs6+a2dbouyrSWlcCPfqDohBcz0Wn/kuUZg10M0sCjwD3AruA+81sV12z7wI97n478CXgM1F3VKTV8oVxIN4KXUMuEqVmKvQ7gV53P+nu48BjwO7aBu7+NXe/FD79DrAp2m6KtF5cF+aqSmvaokSsmUC/AThd87wvXDadjwJ/O59OiSwGA+F1XKI+7b8qqNB1YpFEJ7rbmANm9q+BHuBd06x/AHgAYMuWLVFuWiRyA6MlutozU6fpRy2dNFXoEqlm3qlngM01zzeFy17FzO4GPgHc5+6lRi/k7o+6e4+793R3d19Lf0VaZmC0FNsBUQiHXDSGLhFqJtAPADvMbLuZZYA9wL7aBmb2JuD3CcK8P/puirTe+ZES163Mxfb6WU1blIjNGujuXgb2Ak8AR4HH3f2wmT1sZveFzX4L6AC+aGaHzGzfNC8nsmScGy6yIcZA14lFErWmxtDdfT+wv27ZQzWP7464XyILqjxZob9QZMOq+AJd0xYlajpTVKSB/GiJihNroGvaokRNgS7SwLnhIgAbY67QdT10iZICXaSBaqBvWLkitm0E10OfjO315bVHgS7SwNkWVOjZVIKShlwkQgp0kQbOjRTJpBKsbkvHto1cOkmpXKFS0bCLREOBLtLAueEiG1flMLPYtrEikwRQlS6RUaCLNHB2+HKsc9ABcqng1+/yhMbRJRoKdJEGXh68xOautli3Ua3QFegSFQW6SJ3ixCTnR0psiTnQc+nk1PZEoqBAF6nTdzG4tH+rAv3yuAJdoqFAF6nz8mAQ6LEPuahCl4gp0EXqvHyhNRV6dQy9OKFZLhINBbpInZcHL7MinYz1WugAuZQOikq0FOgidYIZLitinYMOsCKjaYsSLQW6SJ2T+VFu6u6IfTua5SJRU6CL1CiVJ3npwhg71ivQZelRoIvUeHFgjIrDzdd1xr6tFZq2KBFToIvUOHF+FKDFFbpmuUg0FOgiNU70j5IwuLG7PfZtJRNGJpnQQVGJjAJdpMbRsyNsW9dONpxSGLdcOqExdImMAl2kxnN9Q9x+w6qWbW9FJqlAl8g0Fehmdo+ZHTOzXjN7sMH6d5rZM2ZWNrOfir6bIvE7N1zk/EiJ2zetbtk2c+mkhlwkMrMGupklgUeAe4FdwP1mtquu2cvAh4E/i7qDIq3ybN8QAD+wuXWBviKd1CwXiUyqiTZ3Ar3ufhLAzB4DdgNHqg3c/aVwnQ7Xy5J18KVBMskEr7t+Zcu2uSKT5JICXSLSzJDLDcDpmud94bI5M7MHzOygmR3M5/PX8hIisfnmiQHevHXN1HTCVujIphgtlVu2PVneWnpQ1N0fdfced+/p7u5u5aZFZtQ/UuT5cwXeccu6lm63M6dAl+g0E+hngM01zzeFy0SWjS9/7ywA77n1upZutz2TYrSoQJdoNBPoB4AdZrbdzDLAHmBfvN0SaR135/GDfezauJKdG+I/5b9Whyp0idCsge7uZWAv8ARwFHjc3Q+b2cNmdh+Amb3FzPqADwK/b2aH4+y0SFTKkxX+6P+9xNGzI/zc27e1fPud4Rh6peIt37YsP83McsHd9wP765Y9VPP4AMFQjMii118o8vVjeZ48luebJ/KMFMu8ZdsafvxN13Ssf146csGv4Nh4mc5cuuXbl+WlqUAXWcrKkxUOnR7ia8f6efJYnsOvjACwvjPL+1+3gX9x63refet6UsnWnzjdkQ1CfLSkQJf5U6DLsjRVhR/P883jQRWeTBh3bFnNb7x/J3ft7GbXxpWx35VoNtUKfbRYhtZdcUCWKQW6LAvVKvzJY3mePN7P988EVXh3Z5b3vW4Dd+3s5h03d7OqbXFVwZ3ZMNB1YFQioECXJStfKPH143mePNbPN08MMHx5goTBHVvW8Bvv38m7bgmq8ERiYavwmUxV6Ap0iYACXZaMyYpz6PRFnjyW52vHXl2Fv3fXdYu2Cp9JR7ZmyEVknhTosqjlCyW+cTwI8Poq/Nffdwt37Vy/6KvwmVQDfaQ4scA9keVAgS6LSm0V/uSxPN87Mwws7Sp8Jl3tGQAGxxToMn8KdFlwA6OlKzNSTuQZurS8qvCZtGWSZFMJBsdKC90VWQYU6NJyQRU+xNeP9fPk8TzP9QVV+LqOLHfftvyq8JmYGWvbM6rQJRIKdGmJ8yNFvnVigK8fz/ON11gVPpuujowqdImEAl1iUShO8NTJQb7VO8A/9Q5won8UCKrw99waVuE71rG6LbPAPV14Xe1ZBsfGF7obsgwo0CUSE+GJPd86EQT4d08PMVlxcukEb9nWxU+9eRNvv3nda7YKn8na9gwn86ML3Q1ZBhTock2KE5M81zfMgZcGeerFQZ5+aZCx8UkSBm/YtJqPvetG3n7zOu7Y0to7AC1FXe0ZVegSCQW6NGW0VOaZUxf55xcH+ecXBznUN8R4ObiF7M7rOvmJO4IK/G03rn1NHMyM0vrOLJfGJxkpTrBSF+iSeVCgy1UqFefkwBiHTg/x7OkhDp0e4sjZESYrTjJhvP6GVXzobVu5c/taerauYU27xsHnY9OaNgBOD17iddfrCl1y7RToQr5QmgruQ6eHeLZviEJ4KnpHNsXtm1bxi3fdxFu2d3HHljW0Z/W2idLmrhUA9F28rECXedFv5mvIZMV5cWCUw6+McPRsgSNnRzh6doR8IZgyl0wYO6/r5Ed/4HreuHk1b9y8mpu6O0jqIGasNtdU6CLzoUBfhiYrzpmLl+nNF3ihf4ze/lGePzfC8+cKlMJx73TS2LG+k3fu6Oa2jZ3cvmk1b7hhFSsyOoDZaqvb0nRkU5y6oECX+VGgL1HuzsDoOKcvXuL04CVeyI/xQv8oL+RHOTkwNnXAEoJpcTs3dPKzb93KbRtXsuv6ldzU3UEm1fo79MjVzIzbNnZy+JXhhe6KLHEK9EWqUnEGxkr0j5Q4M3SZ04OX6LsYfA5C/DKXJyan2icMtnS1cVN3B++8pZubuzu4aX07N3V36OSdJeD2Tav53985RXmysiC3wpPloalAN7N7gP8GJIE/dPf/XLc+C/wx8GbgAvAz7v5StF1dHi6NlxkcG+fi2ASDl8YZKJQ4XyhyfrjI+ZES50aK9I8U6S+UKNfdCb49k2RzVxtb17bzwzd3s6VrBZu72tjc1caWrjbN917C7tiyhs9+60WeeXmIO7d3LXR3ZImaNdDNLAk8ArwX6AMOmNk+dz9S0+yjwEV3v9nM9gCfBn4mjg4vtErFGR0vM1osUyiWGS1NMFJ89fNC+HhwbJyLl8bDAB9n8NI4xYlKw9ftzKW4bmWODStz3HjTWjaszHFd+LFxVY7NXW2saUsv+D0wJR7v2tlNNpXgL757RoEu16yZCv1OoNfdTwKY2WPAbqA20HcDnwwffwn4H2Zm7v7qEjMCo6Uyw5cnKE9WmJh0JiYrlCediUrwuTxZYaLiTJQrlCtBm+rnqbbh1xYnJqc+Lk9MUpyoTD0uTVQolie5PD4Zfg7WNXOrMLNgul9Xe4au9gwbVua4beNKutozrGnL0NWeDj+H61flaMto9Ou1rCOb4ifu2MQXD57m7tvW8+5b1+uPt8xZMylyA3C65nkf8IPTtXH3spkNA2uBgSg6WetPvn2KT//d85G9XjaVIJdOsiKdJJcOHlefr23PkFsdPM6GyzpyKVbmUnRkU3Tm0nTkUnTmUnTWPG/PJPXLKHP2m+/fyTOnLvLRzx+kI5ti1Yo02XSC6d5J073H9M5b/H75PTv40R+4PvLXbWlZaGYPAA8AbNmy5Zpe466d3XS1p0klEqSSRiaZIJUMHqfDZemkkU4mSCUSpJMWrE8YmVTwOZVMkEkmyKYSulCULBpr2jP81d6389fPvsL3zwxTKJWnppleZZr/fX26FbKorFoRzyUemgn0M8DmmuebwmWN2vSZWQpYRXBw9FXc/VHgUYCenp5reufdtnElt21ceS1fKrLo5dJJPtizmQ/2bJ69sUidZuZHHQB2mNl2M8sAe4B9dW32AR8KH/8U8I9xjJ+LiMj0Zq3QwzHxvcATBNMWP+fuh83sYeCgu+8DPgv8iZn1AoMEoS8iIi3U1Bi6u+8H9tcte6jmcRH4YLRdExGRudApaSIiy4QCXURkmVCgi4gsEwp0EZFlQoEuIrJM2EJNFzezPHDqGr98HTFcViAC6tfcqF9zt1j7pn7NzXz6tdXduxutWLBAnw8zO+juPQvdj3rq19yoX3O3WPumfs1NXP3SkIuIyDKhQBcRWSaWaqA/utAdmIb6NTfq19wt1r6pX3MTS7+W5Bi6iIhcbalW6CIiUmfRBrqZfdDMDptZxcx66tZ93Mx6zeyYmb1/mq/fbmZPhe2+EF76N+o+fsHMDoUfL5nZoWnavWRm3wvbHYy6Hw2290kzO1PTtw9M0+6ecB/2mtmDLejXb5nZ82b2nJn9hZmtnqZdS/bXbN+/mWXDn3Fv+F7aFldfara52cy+ZmZHwvf/rzRoc5eZDdf8fB9q9Fox9G3Gn4sFfjfcX8+Z2R0t6NPOmv1wyMxGzOxX69q0bH+Z2efMrN/Mvl+zrMvMvmJmJ8LPa6b52g+FbU6Y2YcatZmVuy/KD+A2YCfwJNBTs3wX8CyQBbYDLwDJBl//OLAnfPx7wC/E3N/fBh6aZt1LwLoW7rtPAr8+S5tkuO9uBDLhPt0Vc7/eB6TCx58GPr1Q+6uZ7x/4d8DvhY/3AF9owc9uI3BH+LgTON6gX3cBf9Oq91OzPxfgA8DfEtwF763AUy3uXxI4RzBPe0H2F/BO4A7g+zXLPgM8GD5+sNH7HugCToaf14SP18x1+4u2Qnf3o+5+rMGq3cBj7l5y9xeBXoIbWU+x4GaL7ya4YTXA54Efi6uv4fZ+GvjzuLYRg6mbf7v7OFC9+Xds3P3v3b16l+3vENz9aqE08/3vJnjvQPBeeo/FfLNYdz/r7s+EjwvAUYJ79i4Fu4E/9sB3gNVmtrGF238P8IK7X+sJi/Pm7t8guCdErdr30XRZ9H7gK+4+6O4Xga8A98x1+4s20GfQ6KbV9W/4tcBQTXg0ahOldwDn3f3ENOsd+Hszezq8r2or7A3/7f3cNP/iNbMf4/QRgmqukVbsr2a+/1fd/Byo3vy8JcIhnjcBTzVY/TYze9bM/tbMXteiLs32c1no99Qepi+qFmJ/VV3n7mfDx+eA6xq0iWTftfQm0fXM7B+ADQ1WfcLd/6rV/WmkyT7ez8zV+Q+7+xkzWw98xcyeD/+Sx9Iv4H8BnyL4BfwUwXDQR+azvSj6Vd1fZvYJoAz86TQvE/n+WmrMrAP4P8CvuvtI3epnCIYVRsPjI38J7GhBtxbtzyU8RnYf8PEGqxdqf13F3d3MYptauKCB7u53X8OXNXPT6gsE/+6lwsqqUZtI+mjBTbF/AnjzDK9xJvzcb2Z/QfDv/rx+EZrdd2b2B8DfNFjVzH6MvF9m9mHgXwLv8XDwsMFrRIZlEr0AAAH5SURBVL6/Gojs5udRM7M0QZj/qbv/3/r1tQHv7vvN7H+a2Tp3j/WaJU38XGJ5TzXpXuAZdz9fv2Kh9leN82a20d3PhkNQ/Q3anCEY66/aRHD8cE6W4pDLPmBPOANhO8Ff2n+ubRAGxdcIblgNwQ2s46r47waed/e+RivNrN3MOquPCQ4Mfr9R26jUjVv++DTba+bm31H36x7gN4H73P3SNG1atb8W5c3PwzH6zwJH3f13pmmzoTqWb2Z3Evwex/qHpsmfyz7g34SzXd4KDNcMNcRt2v+SF2J/1al9H02XRU8A7zOzNeEQ6fvCZXPTiiO/1/JBEER9QAk4DzxRs+4TBDMUjgH31izfD1wfPr6RIOh7gS8C2Zj6+UfAx+qWXQ/sr+nHs+HHYYKhh7j33Z8A3wOeC99MG+v7FT7/AMEsihda1K9egnHCQ+HH79X3q5X7q9H3DzxM8AcHIBe+d3rD99KNLdhHP0wwVPZczX76APCx6vsM2Bvum2cJDi7/UAv61fDnUtcvAx4J9+f3qJmdFnPf2gkCelXNsgXZXwR/VM4CE2F+fZTguMtXgRPAPwBdYdse4A9rvvYj4XutF/i5a9m+zhQVEVkmluKQi4iINKBAFxFZJhToIiLLhAJdRGSZUKCLiCwTCnQRkWVCgS4iskwo0EVElon/DyNQAVWkxrAUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-10, 10, .0001)\n",
    "plt.plot(x, norm.pdf(f(x+f_mean))*fp(x+f_mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample From Base Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample(f, n):\n",
    "    mean = compute_moment(f)\n",
    "    zs = np.random.normal(size=n)\n",
    "    result = np.array([invert(f, z) for z in zs])\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the Histogram for a Normal Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_normal_histogram(stddev, bin_size):\n",
    "    a = -3*stddev\n",
    "    b = 3*stddev\n",
    "    a = np.floor(a / bin_size)*bin_size\n",
    "    b = np.ceil(b / bin_size)*bin_size\n",
    "    num_bins = int(np.rint((b - a) / bin_size))\n",
    "    bins = np.array(range(num_bins+1))*bin_size + a\n",
    "    weights = np.zeros(num_bins)\n",
    "    for i in range(num_bins):\n",
    "        t = a + bin_size*i\n",
    "        weights[i] = norm.cdf(t + bin_size, scale=stddev) - norm.cdf(t, scale=stddev)\n",
    "    return bins, weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_normal_histogram(stddev, bin_size):\n",
    "    bins, weights = compute_normal_histogram(stddev, bin_size)\n",
    "    plt.hist(bins[:-1], bins, weights=weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAATiklEQVR4nO3df4xl513f8feHdWxIAIM2+0eyu2JW8oI6RhWB0TYtCCEW8BqQl0qOWFOoUS0sVLsJTRFdFzUIV0ZYoJhfDmKVDXKoxXq1SaUVNnGKHNQi1WuPExeyu9lqaod63bQZjONgUOyu8+0f9zgeJnd3zszcmfvjeb+kle495zlnn2fuuedznuf8uKkqJEnt+ZpxV0CSNB4GgCQ1ygCQpEYZAJLUKANAkhp11bgrsB5vfetba25ubtzVkKSp8dRTT/1VVe0aNm+qAmBubo7FxcVxV0OSpkaSv7zcPIeAJKlRBoAkNcoAkKRGGQCS1CgDQJIa1SsAkhxKciHJUpKjQ+Zfk+Shbv6ZJHPd9J1JPpHk5SS/s2qZ70ryF90yv5Uko2iQJKmfNQMgyQ7gfuBGYB64Jcn8qmK3AS9W1XXAfcC93fQvAf8e+Pkhq/5d4GeA/d2/QxtpgCRpY/r0AA4AS1X1TFW9CpwADq8qcxh4oHt9CjiYJFX1t1X1ZwyC4CuSvA34xqp6vAbPo/4w8GObaYgkaX36BMBu4LkV7y9204aWqapLwEvAzjXWeXGNdQKQ5PYki0kWl5eXe1RXktTHxN8JXFXHgGMACwsL/nqNxm7u6MO9yn32V39ki2sibU6fHsDzwN4V7/d004aWSXIVcC3wwhrr3LPGOqWJ03fnv96y0jj06QE8CexPso/BTvoI8BOrypwGbgX+G3Az8Fhd4bcmq+pzSb6Y5J3AGeCfA7+9gfpLW24zO/KVy9oj0KRZswfQjenfCTwKnAdOVtXZJHcnuakrdhzYmWQJeC/wlUtFk3wWeD/w00kurriC6F8CHwSWgP8J/PFomiSNziiP4u0RaNJkmn4UfmFhoXwaqLbLVu2w7QloOyV5qqoWhs2b+JPA0nbb6iP119dvEGjcfBSEJDXKAJBW2M5xes8JaNwMAKkzjh2yIaBxMgAkqVEGgJo3d/ThsR6J2wvQuBgAatqk7HzHHUJqkwEgSY0yANSsSTzinsQ6aXYZAJLUKANATZrkI+1JrptmiwEgSY3yWUBqyrQcXfu8IG0HewCS1CgDQM2YlqP/laaxzpoeBoAkNcoAUBOm+Uh6muuuyWYASFKjDADNvFk4gp6FNmjyGACS1CjvA9DMmrWjZu8N0KjZA5CkRhkAktQoA0AzadaGf1aa5bZpexkAktQoA0Azp4Uj5BbaqK1nAEhSowwASWqUAaCZ0tLQSEtt1dYwACSpUd4JrJnQ6tGwdwdrM+wBSFKjDABJalSvAEhyKMmFJEtJjg6Zf02Sh7r5Z5LMrZh3Vzf9QpIbVkz/10nOJvl0kj9M8rWjaJDa0+rwz0r+DbQRawZAkh3A/cCNwDxwS5L5VcVuA16squuA+4B7u2XngSPA9cAh4ANJdiTZDbwbWKiqbwd2dOUkSdukTw/gALBUVc9U1avACeDwqjKHgQe616eAg0nSTT9RVa9U1bPAUrc+GJyA/rokVwFvBv735poiSVqPPgGwG3huxfuL3bShZarqEvASsPNyy1bV88CvA/8L+BzwUlV9fNh/nuT2JItJFpeXl3tUVy1x6OMN/i20XmM5CZzkmxn0DvYBbwfekuQnh5WtqmNVtVBVC7t27drOakrSTOsTAM8De1e839NNG1qmG9K5FnjhCsv+APBsVS1X1f8DPgr8k400QG2aO/qwR7xD+DfRevQJgCeB/Un2Jbmawcna06vKnAZu7V7fDDxWVdVNP9JdJbQP2A88wWDo551J3tydKzgInN98cyRJfa15J3BVXUpyJ/Aog6t1PlRVZ5PcDSxW1WngOPAHSZaAv6a7oqcrdxI4B1wC7qiq14AzSU4Bn+ymfwo4NvrmSZIuJ4MD9emwsLBQi4uL466Gxsxhjn58PIQAkjxVVQvD5nknsCQ1ygCQpEYZAJoqDv/0599KazEAJKlRBoAkNcoA0NRwSGP9/JvpSgwASWqUASBJjfI3gTXxHMbYHH83WJdjD0CSGmUASFKjDABNNId/Rse/pVYzACSpUQaAJDXKANDEcshi9PybaiUDQJIaZQBIUqMMAE0khyq2jn9bvc4AkKRGGQCS1CifBaSJ4vDE9vD5QAJ7AJLULANAkhplAGhiOPyz/fybt80AkKRGGQCS1CgDQBPBoYjx8W/fLgNAkhplAEhSowwAjZ1DEOPnZ9AmA0CSGmUASFKjegVAkkNJLiRZSnJ0yPxrkjzUzT+TZG7FvLu66ReS3LBi+jclOZXkM0nOJ/nHo2iQJKmfNR8Gl2QHcD/wg8BF4Mkkp6vq3IpitwEvVtV1SY4A9wI/nmQeOAJcD7wd+JMk31pVrwG/CXysqm5OcjXw5pG2TBPPcefJ4gPi2tOnB3AAWKqqZ6rqVeAEcHhVmcPAA93rU8DBJOmmn6iqV6rqWWAJOJDkWuB7geMAVfVqVX1h882RJPXVJwB2A8+teH+xmza0TFVdAl4Cdl5h2X3AMvD7ST6V5INJ3jLsP09ye5LFJIvLy8s9qitJ6mNcJ4GvAr4T+N2qegfwt8BXnVsAqKpjVbVQVQu7du3azjpqCzn8M7n8bNrRJwCeB/aueL+nmza0TJKrgGuBF66w7EXgYlWd6aafYhAIkqRt0icAngT2J9nXnaw9ApxeVeY0cGv3+mbgsaqqbvqR7iqhfcB+4Imq+j/Ac0m+rVvmIHAOSdK2WTMAujH9O4FHgfPAyao6m+TuJDd1xY4DO5MsAe+lG86pqrPASQY7948Bd3RXAAH8K+DBJH8OfAfwK6NrliaZQwyTz8+oDb1+E7iqHgEeWTXtfStefwl412WWvQe4Z8j0p4GF9VRWkjQ63gksSY3q1QOQRsFhhenijWGzzx6AJDXKAJCkRhkA2hYO/0wvP7vZZQBIUqMMAElqlAGgLecQwvTzM5xNBoAkNcoAkKRGGQDaUg4dzA4/y9ljAEhSowwASWqUzwLSlnC4YDb5fKDZYg9AkhplAEhSowwAjZzDP7PPz3g2GACS1CgDQJIaZQBopBwaaIef9fQzACSpUQaAJDXKANDIOCTQHj/z6WYASFKjDABJapTPAtKmOQzQNp8PNL3sAUhSowwASWqUAaBNcfhHr3NbmD4GgCQ1ygCQpEYZANowu/xazW1iuhgAktSoXgGQ5FCSC0mWkhwdMv+aJA91888kmVsx765u+oUkN6xabkeSTyX5o802RJK0PmsGQJIdwP3AjcA8cEuS+VXFbgNerKrrgPuAe7tl54EjwPXAIeAD3fpe9x7g/GYboe01d/Rhu/q6LLeN6dGnB3AAWKqqZ6rqVeAEcHhVmcPAA93rU8DBJOmmn6iqV6rqWWCpWx9J9gA/Anxw882QJK1XnwDYDTy34v3FbtrQMlV1CXgJ2LnGsr8B/ALw5Sv950luT7KYZHF5eblHdSVJfYzlJHCSHwU+X1VPrVW2qo5V1UJVLezatWsbaqcrsXuvPhwmnA59AuB5YO+K93u6aUPLJLkKuBZ44QrLfjdwU5LPMhhS+v4k/3ED9ZckbVCfAHgS2J9kX5KrGZzUPb2qzGng1u71zcBjVVXd9CPdVUL7gP3AE1V1V1Xtqaq5bn2PVdVPjqA9kqSe1gyAbkz/TuBRBlfsnKyqs0nuTnJTV+w4sDPJEvBe4Gi37FngJHAO+BhwR1W9NvpmaDvYpdd6uc1Mtl6/B1BVjwCPrJr2vhWvvwS86zLL3gPcc4V1/ynwp33qIUkaHe8EVi8eyWmj3HYmlwEgSY0yACSpUf4msK7I7rtGwd8Nnkz2ACSpUQaAJDXKANBlOfyjUXObmiwGgCQ1ygDQUB6paau4bU0OA0CSGmUASFKjDAB9Fbvo2mpuY5PBAJCkRnknsL7CozJtJ+8OHj97AJLUKANAkhplAAhw+Efj47Y3PgaAJDXKAJBHYBo7t8HxMAAkqVEGQOM88tKkcFvcfgaAJDXKAJCkRnkncKPsbmsSeXfw9rIHIEmNMgAa5NG/Jp3b6PYwACSpUQZAYzyy0rRwW916BoAkNcoAaMTc0Yc9otLUcZvdWgaAJDXKAGiAR1GaZvZet06vAEhyKMmFJEtJjg6Zf02Sh7r5Z5LMrZh3Vzf9QpIbuml7k3wiybkkZ5O8Z1QNkiT1s2YAJNkB3A/cCMwDtySZX1XsNuDFqroOuA+4t1t2HjgCXA8cAj7Qre8S8G+qah54J3DHkHVqBDxy0qxwWx69Pj2AA8BSVT1TVa8CJ4DDq8ocBh7oXp8CDiZJN/1EVb1SVc8CS8CBqvpcVX0SoKr+BjgP7N58cyRJffUJgN3AcyveX+Srd9ZfKVNVl4CXgJ19lu2Gi94BnBn2nye5PcliksXl5eUe1dXrPGLSrHGbHq2xngRO8vXAR4Cfq6ovDitTVceqaqGqFnbt2rW9FZSkGdbnaaDPA3tXvN/TTRtW5mKSq4BrgReutGySNzHY+T9YVR/dUO01lEdJmmU+MXR0+vQAngT2J9mX5GoGJ3VPrypzGri1e30z8FhVVTf9SHeV0D5gP/BEd37gOHC+qt4/ioZIktZnzQDoxvTvBB5lcLL2ZFWdTXJ3kpu6YseBnUmWgPcCR7tlzwIngXPAx4A7quo14LuBnwK+P8nT3b8fHnHbmuTRv1rhtr55GRyoT4eFhYVaXFwcdzUmll8ItcihoCtL8lRVLQyb553AktQoA2BGePSvVrntb5wBIEmN8kfhp5xHP5KXhm6UPYAp5s5f+vv8TqyPASBJjTIAppRHOtJwfjf6MwCmkBu4dGV+R/rxJPAUcaOW+vPE8NrsAUhSowyAKeHRv7QxfncuzwCYAm7A0ub4HRrOcwATzI1WGh3PCXw1ewCS1CgDYEJ59C9tDb9bbzAAJpAbqLS1/I4NGAATxg1T2h5+1zwJPDHcGKXt1/qJYXsAE8CdvzRerX4HDYAxa3XDkyZNi99FA2CMWtzgpEnW2nfScwBj0NpGJk2Tls4LGADbyB2/ND1aCAKHgLbB3NGH3flLU2qWv7sGwBab5Y1HasWsHsQ5BLRFZnFjkVo3a8NCBsAIudOX2rDyuz7NYeAQ0Ii485faNM3ffXsAmzDNH7yk0ZnWHoEBsE7u9CVdyTSFgQHQgzt9SRsx6WFgAAzhDl/SqK3er0xCIDQfAO7sJY3DsH3PdodCrwBIcgj4TWAH8MGq+tVV868BPgx8F/AC8ONV9dlu3l3AbcBrwLur6tE+6xw1d/SSJt3l9lNbFQxrXgaaZAdwP3AjMA/ckmR+VbHbgBer6jrgPuDebtl54AhwPXAI+ECSHT3XKUnaQn3uAzgALFXVM1X1KnACOLyqzGHgge71KeBgknTTT1TVK1X1LLDUra/POiVJW6jPENBu4LkV7y8C/+hyZarqUpKXgJ3d9MdXLbu7e73WOgFIcjtwe/f25SQXetR5mLcCf7XBZSfNrLRlVtoBtmUSzUo7yL2basu3XG7GxJ8ErqpjwLHNrifJYlUtjKBKYzcrbZmVdoBtmUSz0g7Yurb0GQJ6Hti74v2ebtrQMkmuAq5lcDL4csv2WackaQv1CYAngf1J9iW5msFJ3dOrypwGbu1e3ww8VlXVTT+S5Jok+4D9wBM91ylJ2kJrDgF1Y/p3Ao8yuGTzQ1V1NsndwGJVnQaOA3+QZAn4awY7dLpyJ4FzwCXgjqp6DWDYOkffvL9n08NIE2RW2jIr7QDbMolmpR2wRW3J4EBdktQaHwctSY0yACSpUU0FQJL/kOTPkzyd5ONJ3j7uOm1Ekl9L8pmuLf8pyTeNu04bleRdSc4m+XKSqbtkL8mhJBeSLCU5Ou76bEaSDyX5fJJPj7sum5Fkb5JPJDnXbVvvGXedNirJ1yZ5Isl/79ryyyNdf0vnAJJ8Y1V9sXv9bmC+qn52zNVatyQ/xOBKq0tJ7gWoqn875mptSJJ/AHwZ+D3g56tqccxV6q17pMn/AH6Qwc2MTwK3VNW5sVZsg5J8L/Ay8OGq+vZx12ejkrwNeFtVfTLJNwBPAT82jZ9L90SFt1TVy0neBPwZ8J6qenyNRXtpqgfw+s6/8xZgKtOvqj5eVZe6t48zuI9iKlXV+ara6N3d4zZTjzSpqv/C4Cq+qVZVn6uqT3av/wY4zxtPIJgqNfBy9/ZN3b+R7beaCgCAJPckeQ74Z8D7xl2fEfgXwB+PuxKNGvaYlKnc0cyqJHPAO4Az463JxnUP0Hwa+Dzwn6tqZG2ZuQBI8idJPj3k32GAqvrFqtoLPAjcOd7aXt5a7ejK/CKD+yseHF9N19anLdKoJfl64CPAz63q/U+Vqnqtqr6DQU//QJKRDc9N/LOA1quqfqBn0QeBR4Bf2sLqbNha7Ujy08CPAgdrwk/krOMzmTY+0mRCdePlHwEerKqPjrs+o1BVX0jyCQaP1h/JifqZ6wFcSZL9K94eBj4zrrpsRvdjOr8A3FRVfzfu+jTMR5pMoO7E6XHgfFW9f9z12Ywku16/yi/J1zG44GBk+63WrgL6CPBtDK46+UvgZ6tq6o7YukduXMPggXsAj0/j1UwASf4p8NvALuALwNNVdcN4a9Vfkh8GfoM3Hmlyz5irtGFJ/hD4PgaPUf6/wC9V1fGxVmoDknwP8F+Bv2DwXQf4d1X1yPhqtTFJ/iGD31rZweCA/WRV3T2y9bcUAJKkNzQ1BCRJeoMBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhr1/wFwaxRmlSsj0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_normal_histogram(1.0, 0.025)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot a Histogram With Normal Overlay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_histogram(xs, bin_size, normal_stddev=None, plot=plt):\n",
    "    a = np.min(xs)\n",
    "    b = np.max(xs)\n",
    "    a = np.floor(a / bin_size)*bin_size\n",
    "    b = np.ceil(b / bin_size)*bin_size\n",
    "    num_bins = int(np.rint((b - a) / bin_size))\n",
    "    alpha = 1.0\n",
    "    if normal_stddev:\n",
    "        alpha = 0.5\n",
    "        bins, weights = compute_normal_histogram(normal_stddev, bin_size)\n",
    "        plot.hist(bins[:-1], bins, weights=weights, alpha=alpha, density=True)\n",
    "    plot.hist(xs, bins=num_bins, range=(a, b), alpha=alpha, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Least Squares Errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_target(X, b, f, m):\n",
    "    n = X.shape[0]\n",
    "    mean = compute_moment(f)\n",
    "    stddev = compute_stddev(f)\n",
    "    e = np.zeros(n)\n",
    "    for _ in range(m):\n",
    "        ei = sample(f, n)\n",
    "        e += (ei - mean)\n",
    "    e /= (np.sqrt(m) * stddev)\n",
    "    return np.dot(X, b) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_least_squares_errors(X_train, X_test, b, f, m):\n",
    "    y_train = compute_target(X_train, b, f, m)\n",
    "    y_test = compute_target(X_test, b, f, m)\n",
    "    model = LinearRegression(fit_intercept=False)\n",
    "    model.fit(X_train, y_train)\n",
    "    y_train_pred = model.predict(X_train)\n",
    "    rss = np.sum((y_train - y_train_pred)**2)\n",
    "    noise_var = rss / (X_train.shape[0] - X_train.shape[1])\n",
    "    y_test_pred = model.predict(X_test)\n",
    "    return (y_test - y_test_pred), np.sqrt(noise_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Some Error Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20000\n",
    "b = np.array([-2, 3])\n",
    "X_train = np.random.uniform(-1, 1, size=(N, len(b)))\n",
    "X_test = np.random.uniform(-1, 1, size=(N, len(b)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_error_terms = [1, 2, 5, 25]\n",
    "error_data = []\n",
    "for m in num_error_terms:\n",
    "    error_data.append(compute_least_squares_errors(X_train, X_test, b, f, m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAeFklEQVR4nO3df7DddX3n8edrg+Ds4ghL0qrhR9BC19Ra2N6BMMw4Uo0Ga6Fb2DHY7sJKN93dIl2UOgRm0GHbCZW1sVV2NSUZtKMga+t6nQVp0DjoetPNjSIlQbKZWxtCC8SAVISVCbz3j/s94ZuTc+49P77n+/18v9/XY+ZMzvfHuedzL+fN6/P5nM/5HkUEZmZmqfknVTfAzMysFweUmZklyQFlZmZJckCZmVmSHFBmZpYkB5SZmSXJAWVmZklyQNlRJL1W0rSkv5cUklZU3SazOpL0q5K+JelHkh6XdJukV1XdrrpwQFkvLwFfBS6puiFmNfdq4A+A1wFvBJYDt1TaohpxQCVO0g8k/b6kByX9RNJmST8r6R5JP5Z0n6QTi3zOiHgiIv4bsKPIn2tWpYpq6fMR8dWIeC4ingb+DDi/yOdoMgdUPVwCrAbOBH4NuAe4HljG/H/Dq3s9SNKp2dRCv9t7y/oFzBJRdS29Bdg1/q/RDsdU3QAbyCci4gkASd8EnoyI72bbXwLe1utBEbEPOKG0Vpqlr7JakrQauBw4d5yf0yYeQdXDE7n7z/fYPr7c5pjVViW1JGkV8Hng0ojYM4nnaCIHVINl0xLPLnD7zarbaFYH49SSpLOBaeB9EfG18lpdf57ia7BsWmKkHqGkVwJLss3jJL0yIv5fYY0zq5FRa0nSm5hfEfv+iPhK4Q1rOI+grJ/ngWez+9/Pts1sOB9kfgHG5txoy4skBiR/YaGZmaXIIygzM0uSA8rMzJLkgDIzsyQ5oMzMLEnJLTNfunRprFixoupmmC1q586dP4yIZVW3ox/XktVFv1pKLqBWrFjB7Oxs1c0wW5Skv6u6DQtxLVld9KslT/GZmVmSHFBmZpYkB5SZmSXJAWXV27ah6haYNUPDaskBZWZmSXJAmZlZkhxQZmaWJAeUpaFhc+dmNj4HlJmZJckBZVYxSVskPSnpoT7H3yrpGUkPZLcby26j1UADZyGSu9SRWQvdDnwS+OwC53wzIt5dTnPM0uARlFnFIuJ+4Kmq22E11sDREzigzOriPEnfk3SPpF/od5KkdZJmJc0eOHCgzPaZFc4BZZa+7wCnRcQvAZ8A/me/EyNiU0RMRcTUsmXJfhOI2UAcUGaJi4h/jIhns/t3A6+QtLTiZplNnAPKqtXQufMiSXqNJGX3z2G+bg9W2yqzyfMqPrOKSboDeCuwVNJ+4MPAKwAi4lPApcB/lHQIeB5YGxFRUXPNSuOAMqtYRFy2yPFPMr8M3axVPMVnZmZJckCZmVmSHFBmZk3SoIVHDiirToMKycyKN/GAknSKpG2SdkvaJen3Jv2cZmat0eCOXhmr+A4BH4yI70h6FbBT0taI2F3Cc5uZWU1NfAQVEf8QEd/J7v8YeBhYPunnNTOzeiv1PShJK4Czgb/u2u8LXJqZDavf9F5Dpv1KCyhJxwN/AfzniPjH/DFf4NLMzLqVElCSXsF8OH0uIv6yjOc0M7N6K2MVn4DNwMMR8ceTfj4zM2uGMkZQ5wP/BvgVSQ9kt3eV8LxmZlZjE19mHhHfAjTp5zEzsy7bNsAF66tuxch8JQkzsyZqwEo+B5SZmSXJAWVmZklyQFk1GjD9YFapFtSQA8rMzJLkgDKrmKQtkp6U9FCf45L0p5L2SnpQ0r8su41mVXBAWXlaMCUxotuBNQscvxA4I7utA/57CW0yq5wDysrhcOorIu4HnlrglIuBz8a87cAJkl5bTuvMquOAsnI5qEaxHHg0t72fPl9Z428GsCZxQJk1iL8ZoCWG6ejVuFPogDJL32PAKbntk7N9Zo3mgLLy1bhHV5Fp4N9mq/lWAc9ExD9U3SizSZv4xWLNbGGS7gDeCiyVtB/4MPAKgIj4FHA38C5gL/Ac8O+qaalZuRxQZhWLiMsWOR7A75bUHLNkeIrPzMyS5IAyM7MkOaDMzCxJDigzM0uSA8omz8vKzWwEDigzswRt3LqHjVv3VN2MSnmZuZlZInoFUn7fNavPbNWMhEdQZmYJGGS01LYRlQPKzKxCg07lrdq3CYCZuYOTblIyHFBmZhUZdUTUlpByQJmZ1URnFAXtCKmJB5SkLZKelPTQpJ/LEtSS760xG9Ywo6d8MOU1PaTKGEHdDqwp4XnMzGqhyMUOTQ6piQdURNwPPDXp5zEzq4O2rcQbRxLvQUlaJ2lW0uyBAweqbo6Z2URMKpwWHUXVdPo8iYCKiE0RMRURU8uWLau6OVaUmhaFWQr6ve/UTxOn+pIIKDOzpht19DRMUDUtpBxQZmYJGXbk1GRlLDO/A5gBfl7SfklXTvo5zepG0hpJj0jaK+m6HsevkHRA0gPZ7beraKeNpozRU0eTRlETv1hsRFw26ecwqzNJS4BbgdXAfmCHpOmI2N116hci4qrSG2hjGSWcxh1Fzcwd5LzXnzTWz0iBr2ZuVr1zgL0RMQcg6U7gYqA7oKxmivgwbps5oKwQvQpx1b75qYYm9OQmbDnwaG57P3Buj/MukfQWYA9wTUQ82n2CpHXAOoBTTz11Ak21QQ0aTqv2bWL7qesKf/4mjKIcUDaWQYqwMyde92Kp2FeAOyLip5J+B/gM8CvdJ0XEJmATwNTUVJTbRLNiOaBsJKPMqzuo+noMOCW3fXK277CIyL/zfRvw0RLaZRPWmdab1PRe3UdRXmZuQynia6ibtMqoIDuAMySdLulYYC0wnT9B0mtzmxcBD5fYPhtSSpczqnO9OaCsEnUumqJFxCHgKuBe5oPnrojYJekmSRdlp10taZek7wFXA1dU01pbzDDvPZWlrvXmKT4byLA9wkGKz1N+L4uIu4G7u/bdmLu/HlhfdrusWJ0FEV6xNxiPoGxRk56uqGvvzqxbSlN7TeCAsgWVVXAOKau7hWolP2KqavRUx/B0QJmZjWmQcEphWq9uIeWAsp7GWa03aiHOzB2sXQGZ1U2daswBZUep+gVc9fObDSP/es2Pljq3FNWlxhxQdoRUXriptMNsIXV+ndah7Q4oM7MCpTpqqiMHlAHFXCGiCPniTqVNZr10XpspLYJoGgeUFRICkyxOh5SlpimvydR/DwdUyxX9Ap1UUKVeSNY+KXy2qQgp15YDqsXqEk4dKReStUP3tHOdgykv1el0B1RLFTWtV/ZS2hSLyNph49Y9jQmkflKrL18stoVGufBr5xs/UyjQjVv3cM3qM6tuhrXIzOZrIaEaaAsHVIssFkyLFV5Khdn5XRxUNlHbNhy+TmRKr/9JSqkD6Cm+lugVTsMUXKrFmdqUhDVILpyaYJgaTuU9KY+gGm7cUVMdeDRlhdi24YjNJoXTqKoeTTmgGqr7Q4RtUHUxWU1t2wAXvPxdkA6mI1VZVw6optm2gY2HLgHqG075RRnD8mjKhtIZNTVsOq9o+ZmYMmurlICStAb4E2AJcFtE3FzG8zZSp6AuWH/E/ZnN1x4+ZRX1DKYiVVVQo1qsRiQdB3wW+GXgIPCeiPhB2e2sta4pvI42BdM4nb+OMkdUEw8oSUuAW4HVwH5gh6TpiNg96eeupc50Q59iyp93uLDmrl343JZLPawGrJErgacj4uckrQX+CHhP+a1NUL7T1tlepIbaFEqT0P3e9qTqShExkR98+Amk84CPRMQ7s+31ABHR89UzNTUVs7OzE21TKbrmtY/YP6S2FtO4Pb1hjFJgknZGxNS4zz1IjUi6NztnRtIxwOPAsliggBtXS72Cx0E0kjJqa5ia6ldLZUzxLQcezW3vB87NnyBpHdD5iz0r6ZES2jUpS4Efzt+9vtKGVCj3NxjHx8b/EQP6wGgPO62gp1+0RvLnRMQhSc8AJ9H1d25QLXW9hq7v+rd7f2MVVEvdJl9bQ9ZUz1pKYpFERGyCZrxxImm2iF51nflvUJ2m1JJfQ/Pa/nco44O6jwGn5LZPzvaZ2bxBauTwOdkU36uZXyxh1lhlBNQO4AxJp0s6FlgLTJfwvGZ1MUiNTAOXZ/cvBb6+0PtPZk0w8Sm+bL78KuBe5pfQbomIXZN+3grVfnqlAP4bDKFfjUi6CZiNiGlgM/DnkvYCTzEfYk3m19C8Vv8dJr6Kz8zMbBS+WKyZmSXJAWVmZklyQE2ApFskfV/Sg5K+JOmEqttUFklrJD0iaa+k66puj9WX68h15ICajK3AmyLizcAeoMclJZond8meC4GVwGWSVlbbKqsx11HL68gBNQER8VcRcSjb3M7851ra4Bxgb0TMRcQLwJ3AxRW3yWrKdeQ6ckBN3vuAe6puREl6XbJneUVtsWZxHbVQEpc6qiNJ9wGv6XHohoj4cnbODcAh4HNlts2sLlxHthAH1Igi4u0LHZd0BfBu4G0t+sS/L2tlQ3Ed9eQ6yniKbwKyL5/7EHBRRDxXdXtK5MtaWWFcR64jX0liArLL0RzHyxfz3B4R/6HCJpVG0ruAj/PyJXv+sOImWU25jlxHDigzM0uSp/jMzCxJDigzM0uSA8rMzJLkgDIzsyQ5oMzMLEkOKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7MkOaDMzCxJDig7iqS3SnpJ0rO52+VVt8ssdZJ+VdK3JP1I0uOSbpP0qtzx2yW90FVbS6psc8ocUNbP30fE8bnbZ6pukFkNvBr4A+B1wBuB5cAtXed8tKu2Xiy7kXXhgEqcpB9I+n1JD0r6iaTNkn5W0j2SfizpPkknVt1Os9RUUTsR8fmI+GpEPBcRTwN/Bpxf5HO0iQOqHi4BVgNnAr8G3ANcDyxj/r/h1b0eJOnUbKqh3+29Czznz0h6QtLfStoo6Z8V/DuZlaGK2sl7C7Cra99/kvSUpJ2SLhnpt2qJY6pugA3kExHxBICkbwJPRsR3s+0vAW/r9aCI2AecMMLzfR84K/v3NOAzwB8DvzPCzzKrUtm1c5ik1cDlwLm53X8KfBB4BngH8AVJj0fE/x7nuZrKI6h6eCJ3//ke28cX+WQR8XhE7I6IlyLib4EPMd8TNaubUmunQ9Iq4PPApRGxp7M/Ir4TEQcj4lBE3A18DviNSbShCRxQDZZNUzy7wO03B/xRgV8r1iLj1I6ks4Fp4H0R8bVFnioAFdn2JvEUX4Nl0xRD9xAlXQDMAfuAk4GbgS8X2zqzdI1RO28Cvgq8PyK+0uP4pdnx54C3A7/F/Htj1oN7xdbL2cC3gZ9k//4Nfd5MNrMjfJD5BRibc6Ot/CKJ3wMeA37E/PLzfx8R3yi/mfWgiKi6DWZmZkfxCMrMzJLkgDIzsyQ5oMzMLEkOKDMzS1Jyy8yXLl0aK1asqLoZZovauXPnDyNiWdXt6Me1ZHXRr5aSC6gVK1YwOztbdTPMFiXp76puw0JcS1YX/WrJU3xmZpYkB5SZmSXJAWVmZklyQNlkbNtQdQvM6sU1cxQHlE1ed+G5EM1sAA4oMzNLkgPKijPMyMijKDNbhAPKJschZGZjcEDZZC0WUg4xs5d16sF1ATigbFzbNrx862yb2eBcM305oKx4Ljiz0bl+DkvuWnyWro1b9xyxfc0xf1HcD9+2AS5YX9zPM0uZX+8DGWgEJWmNpEck7ZV0XY/jGyU9kN32SPpR7tiLuWPTRTbeytMdTgAzcweZmTs43A/q1Tt0j9HsaPm6aGmNLDqCkrQEuBVYDewHdkiajojdnXMi4prc+e8Hzs79iOcj4qzimmxl6xVOeTNzBznv9ScN/4NbVnSS1gB/AiwBbouIm7uObwQuyDb/KfAzEXFCduxF4G+yY/si4qJyWm2TMDN3kO2H5utq1b75+6v2jVhHDTbICOocYG9EzEXEC8CdwMULnH8ZcEcRjbPqLRZOHUOPpFom19G7EFgJXCZpZf6ciLgmIs7KOnSfAP4yd/j5zjGHU711amrVvk1HHes5K9GyjlzeIAG1HHg0t70/23cUSacBpwNfz+1+paRZSdsl/Xqfx63Lzpk9cODAgE23SesXTqv2bepbXNaXO3q2YE3Z0YpexbcW+GJEvJjbd1pETAHvBT4u6Q3dD4qITRExFRFTy5Yl+wWlhgtpDBPv6GWPdWevSVo8eoLBAuox4JTc9snZvl7W0tXri4jHsn/ngG9w5PtTlqhBp/a6eRRViJE6euDOXsqGqameddTCsBokoHYAZ0g6XdKxzBfPUavxJP0L4ERgJrfvREnHZfeXAucDu7sfa2npVUjDjJwcUj25o9dio3T4XEcDBFREHAKuAu4FHgbuiohdkm6SlH+zdi1wZ0REbt8bgVlJ3wO2ATfnV/9ZvXh6byzu6LXUqLMR0COkWjaKGuiDuhFxN3B3174bu7Y/0uNx3wZ+cYz2WcnGKaa8kZeeN1REHJLU6egtAbZ0OnrAbER0wqpfR+/Tkl5ivlPpjl5DuNO3MF/qyA4rKpw6juj9DdLzy1/Tr4Ei4u6IODMi3hARf5jtuzEXTkTERyLiuq7HfTsifjEifin7d3PZbbfRFFFTbX4/ygFl5WlJUZktxKOmwTmgDJjc5zP8Rq+11UIfyB1WW+vIAWVmVrBJdPjaGFIOKCv8vSczsyI4oMzMCtRraq+o951m5g626r1cB1TLDTJ68vtQZuPz4ojhOaCsLxeU2XCKXBjRz+EOXwtGUg4oM7MJcSdvPA6oFitzcYSn+cxsWA4oM7Oaacs0nwOqpbpHT91TEZOYmvAoysyG4YCy9DS8V2jNVMYCibZxQLVQlR/M9SjKbEIa2LFzQNlh7vmZjaaK0VMbOnsOKDtCGQU2UGE1sDdoNhENrhUHlJnZGDxlPjkDBZSkNZIekbRX0nU9jl8h6YCkB7Lbb+eOXS7p/2a3y4tsvJmZ5TRsNLVoQElaAtwKXAisBC6TtLLHqV+IiLOy223ZY/858GHgXOAc4MOSTiys9Ta0Xr29MpaYd2t6z68fd/ZsIhoWTB2DjKDOAfZGxFxEvADcCVw84M9/J7A1Ip6KiKeBrcCa0ZpqVm/u7DVPvsNX1SKjJnf2Bgmo5cCjue392b5ul0h6UNIXJZ0yzGMlrZM0K2n2wIEDAzbdhpXa9z41ubD6cGfPbAhFLZL4CrAiIt7MfOF8ZpgHR8SmiJiKiKlly5YV1CRbjJeVl86dvQZJqcN3RGevQdN9gwTUY8Apue2Ts32HRcTBiPhptnkb8MuDPtbMjuDOno1v24ZGBNUgAbUDOEPS6ZKOBdYC0/kTJL02t3kR8HB2/17gHZJOzObL35Hts5Kl1NvLa9k0nzt7ZkNYNKAi4hBwFfPB8jBwV0TsknSTpIuy066WtEvS94CrgSuyxz4F/BfmQ24HcFO2z6yN3NlriO4rR6QwXd7Ezt4xg5wUEXcDd3ftuzF3fz2wvs9jtwBbxmijjWmhpeWpFNZ5rz+p6mZMXEQcktTp7C0BtnQ6e8BsREwz39m7CDgEPEWusyep09kDd/asBQYKKDMrhjt79bfYV9VUqWmdPV/qyMzMkuSAarhUF0d0a+L8uTVPHeqpSbXkgGqZVfs2JTUlYVY3taqfmi81d0A1WMpz5b00qednzZd6PTWBA8rMbAApXHdvUE25soQDqqHyn9NIvZjyPIoysw4HlJnZAOrU0YNmdPYcUC1Rt+ICaj01YWbjc0A1UB2Wwi6kCT0/a5buSxvVTk07ew4oMzNLkgOqYeq00mghHkVZSupcS3XmgDIzW8DM5murbsLI6t7Rc0A1SN3fezIzy3NAmZn10YROX51HUQ6ohmhCIXWbmTtY29VHZsmpYS05oMzMemhSp6+uo6iBAkrSGkmPSNor6boexz8gabekByV9TdJpuWMvSnogu013P9bG12vlXlNWHdW1sHpxHVnlajaKWjSgJC0BbgUuBFYCl0la2XXad4GpiHgz8EXgo7ljz0fEWdntooLabQtoSjg1ieuoXmr/wdwe6tjZG2QEdQ6wNyLmIuIF4E7g4vwJEbEtIp7LNrcDJxfbTOunSdMQ/TTkd3Qd1UyTwqmuBgmo5cCjue392b5+rgTuyW2/UtKspO2Sfr3XAySty86ZPXDgwABNsg4XUW1MvI7AtVSEhnSIeqrbwqNCF0lI+i1gCrglt/u0iJgC3gt8XNIbuh8XEZsiYioippYtW1ZkkxqtKVeNGEST/6fRbdQ6AtfS2Gr0P+82GCSgHgNOyW2fnO07gqS3AzcAF0XETzv7I+Kx7N854BvA2WO016yuXEc1UMf3aYZVp99xkIDaAZwh6XRJxwJrgSNWEUk6G/g080X1ZG7/iZKOy+4vBc4HdhfVeGuXmo+iXEeJ84xEeo5Z7ISIOCTpKuBeYAmwJSJ2SboJmI2IaeanIo4H/ockgH3ZSqM3Ap+W9BLzYXhzRLiwCtb0YmoC15HZ8BQRVbfhCFNTUzE7O1t1M9K1bQNcsL6Ry2AXsv3UdYfvX7P6zApb8jJJO7P3hZLkWhpcvp62n7quNXV13pX/teomAP1ryVeSqKGNW/e0poB6qcv0hNVHm+spZQ4oM7OcNoVV6l8l4oCqmZm5g427nNEgun9Xj6KsKKn/T3rSUq4lB1QdZJ/NSPmFVAX/PczGl3JH1wFVIym/kMxqpUenz/WVHgdUTbR9GqIfj6JsVG1fbJSXah05oFKVu+RKnT75XYVUi8vS5Zo60qp9m5LsBDugUpQLJ//PdzD+O9nAfL292nBApcbFYzYZWW11Rk+e3jtaah09B1TCZjZf6yIaQmrFZYnYtsErYYeQ0t/IAZUSv+80tpSKy9Limlpcp0OcSh05oBIzM3fQhTSmVIrL0uFpvcGlFFIOqKrlph8cTMVJobgsAds2uK5GkEpIOaASkeISz9QM2/uturisAp1pcnf6CjOz+drKaskBNSn51XjdK/M6oyYX0cRt3LrHQdVS7vSNJ4XpUAdUkRZbIp4LpQ6H03BGLRqHVMN1rdRzXRWvis7eot+oa0PKvlDw8P38/oyLpxr54krlSw9tBH1q7PBCCKrv+TdFrw5hp47KqKGBRlCS1kh6RNJeSdf1OH6cpC9kx/9a0orcsfXZ/kckvbO4pldokJFSD16hV4wiph6qmvpzLRWrU1Ouq/J1amiSdbToCErSEuBWYDWwH9ghaToidudOuxJ4OiJ+TtJa4I+A90haCawFfgF4HXCfpDMj4sWif5FC5Htmvfb3Gx314IKph+7immSvsFW1NIp+9dc5lpmZOwhzfn+pTJ1O4fZT17Fq36Yj/oXJzU4oIhY+QToP+EhEvDPbXg8QERty59ybnTMj6RjgcWAZcF3+3Px5/Z5vamoqZmdnx/qlgIXDBl4OnAI5lMrXKZBJWKzQJO2MiKlBf15ta2lYvWpsoVrsPrcH11b6etXioGHVr5YGeQ9qOfBobns/cG6/cyLikKRngJOy/du7Hru8R+PWAZ3f7llJjwzQrgFcP+yxpcAPi3nuwrhNC/pYfqPQdn1g8VNOG/JH1riWRnG4xpbC9Yv8d1moVguX0OsXaEx7PnbUngFqqKNnLSWxSCIiNkH172xKmh2mR1wGt2lwqbarTKnUUl5q/13cnoWl1J5BFkk8BpyS2z4529fznGxa4tXAwQEfa9YWriWzIQwSUDuAMySdLulY5t+one46Zxq4PLt/KfD1mH9zaxpYm61MOh04A/g/xTTdrHZcS2ZDWHSKL5sHvwq4F1gCbImIXZJuAmYjYhrYDPy5pL3AU8wXHtl5dwG7gUPA7ya+6iipqZGM2zS4VNsFtK6W8lL77+L2LCyZ9iy6is/MzKwKvtSRmZklyQFlZmZJckDlSLpF0vclPSjpS5JOqLpNAJL+taRdkl6SVOnyz8Uu1VMFSVskPSnpoarbYr2lVlup1FRK9ZRiHTmgjrQVeFNEvBnYA/S57krpHgJ+A7i/ykbkLtVzIbASuCy7BE/VbgfWVN0IW1BqtVV5TSVYT7eTWB05oHIi4q8i4lC2uZ35z5pULiIejogKrwhw2DnA3oiYi4gXgDuBiytuExFxP/Mr3ixRqdVWIjWVVD2lWEcOqP7eB9xTdSMS0+tSPUddbsdsEa6tea6nRSRxqaMySboPeE2PQzdExJezc25g/rMmn0upXWYpS622XFP117qAioi3L3Rc0hXAu4G3RYkfElusXYnw5Xasr9RqqwY15XpahKf4ciStAT4EXBQRz1XdngQNcqkes6O4tnpyPS3CAXWkTwKvArZKekDSp6puEICkfyVpP3Ae8L+y7wIqXfYmd+dSPQ8Dd0XErirakifpDmAG+HlJ+yVdWXWb7ChJ1VYKNZVaPaVYR77UkZmZJckjKDMzS5IDyszMkuSAMjOzJDmgzMwsSQ4oMzNLkgPKzMyS5IAyM7Mk/X9qYxzWt551IwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, plots = plt.subplots(nrows=2, ncols=2)\n",
    "index = 0\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        plot = plots[i, j]\n",
    "        m = num_error_terms[index]\n",
    "        errs, noise_stddev = error_data[index]\n",
    "        plot_histogram(errs, 0.025, noise_stddev, plot)\n",
    "        plot.set_title('m = %d' % m)\n",
    "        index += 1\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}