lessons/Leave-one-out resampling.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:5ad418f924d11d61a3a1b0213375e0856af9e34b2e38607d288908944c14785b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "from numba import jit\n",
      "import statsmodels.api as sm\n",
      "from statsmodels.formula.api import ols\n",
      "from pandas import DataFrame"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%precision 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "u'%.4f'"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Jackknife estimate of parameters\n",
      "\n",
      "This shows the leave-one-out calculation idiom for Python. Unlike R, a -k index to an array does not delete the kth entry, but returns the kth entry from the end, so we need another way to efficiently drop one scalar or vector. This can be done using Boolean indexing as shown in the examples below, and is efficient since the operations are on *views* of the origianl array rather thna *copies*. Note also that "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def jackknife(x, func):\n",
      "    \"\"\"Jackknife estimate of the estimator func\"\"\"\n",
      "    n = len(x)\n",
      "    idx = np.arange(n)\n",
      "    return np.sum(func(x[idx!=i]) for i in range(n))/float(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Jackknife estimate of standard deviation\n",
      "x = np.random.normal(0, 2, 100)\n",
      "jackknife(x, np.std)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "1.8943"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def jackknife_var(x, func):\n",
      "    \"\"\"Jackknife estiamte of the variance of the estimator func.\"\"\"\n",
      "    n = len(x)\n",
      "    idx = np.arange(n)\n",
      "    j_est = jackknife(x, func)\n",
      "    return (n-1)/(n + 0.0) * np.sum((func(x[idx!=i]) - j_est)**2.0 for i in range(n))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# estimate of the variance of an estimator \n",
      "jackknife_var(x, np.std)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "0.0169"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### For fun, try optimizing jackknife"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit jackknife(x, np.std)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numba import jit\n",
      "@jit\n",
      "def jackknife_numba(x, func):\n",
      "    \"\"\"Jackknife estimate of the estimator func\"\"\"\n",
      "    n = len(x)\n",
      "    idx = np.arange(n)\n",
      "    return np.sum(func(x[idx!=i]) for i in range(n))/float(n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Leave one out cross validation (LOOCV)\n",
      "\n",
      "LOOCV also uses the same idiom, and a simple example of LOOCV for model selection is illustrated."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a, b, c = 1, 2, 3\n",
      "x = np.linspace(0, 5, 10)\n",
      "y = a*x**2 + b*x + c + np.random.normal(0, 1, len(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(12,4))\n",
      "for deg in range(1, 5):\n",
      "    plt.subplot(1, 4, deg)\n",
      "    beta = np.polyfit(x, y, deg)\n",
      "    plt.plot(x, y, 'r-o')\n",
      "    plt.plot(x, np.polyval(beta, x), 'b:')\n",
      "    plt.title('Degree = %d' % deg)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAAEKCAYAAAAfA+2IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FFX7xvFvSIDQexMLGjv6CooKoiEqCSh20Z8VRaVY\nAogFJIAoIAqIFBUVxQK+CoIgEBRiCZH+goBUgQAiSO8tIWV+f5wlJCFlk+zuzO7en+vKRTK7O/sA\nN8OzM2fOAREREREREREREREREREREREREREREREREREREREREREREREREREREREREREJVluA48Bh\n4AAwD+gEhNhYk7fVBaYB24FM4Fx7y5Fi2ELw5bYNMBfz+90BjAEq2lqRFMcWgi+7NwN/Yn6/+4HZ\nwOW2ViRFtYXgy212YzH9wgV2FyIlsxm4xfV9JeBOYBPmL9gbSnlpv0VRG+gMNEVNr78Kxtw+DMQA\n4UBVYCYw2taKpDiCMbu1gfqu70sD7wAL7StHiiEYc3vKjUAikIGaXr+XPcinXIv5y23o+rksMBT4\nG9iJ+Y82PNvzXwX+BbYBz5Dz09AXrufPBI663ussYDKwG/OPJjbbvkKAnsBGYC8wAahWot9h/sJQ\n0+uvgjm3p9yLOXsm/iXYs1sWGARM8eJ7iOcFa27DgD+AK9GZ3oCQV5DBhLaT6/v3gKmYs0sVMUMD\n3nI91hpzqfUyoBwwnjODfBBo5vq5HLAU6I0J0/lAMuYMFkBXYD4m7KWBj4D/5lP7uZjLLPl9PVTI\n711Nr/8K5tyeMryA9xDnCtbsnnptBubDWvUCnivOE6y5fcX1+wI1vQEhvyAvAF7DfJo6Ss6/6GaY\nT11gLm0MzPZYBGcG+Ytsj1+P+UeS3WucvkSyNlc99YCTeOdSh5pe/xXMuQWIxoyNvNBL+xfvCfbs\nVgPGAT94af/iHcGY23OADZjhHBAgTW+Y3QU41NmY/1RrAuUxn7hOCeF0sOoBi7M9ti3XfizMDWOn\nnIf5ZHYg27ZQICnb41Mw4TolHaiD+ZQoUpBgyG1T4GvgfsylPQkMwZBdXHW87Np3ZcyNUeK/Ajm3\nw4E3gSOcvmHP72/cU9N7pmsxYZsL7ANOYO60zStIOzCfhk45J4/nWNm+34r5xHhxPu+9FWiP+fRY\nmHOB1QU83hH4xo39SGAIhtw2xpwhexL4zY33Ev8QDNnNrjSmWUl147niXIGe21uA5sDgbNsWAF2A\nb914X3GgzcCtru8rA3dgzh59ke05wzGDxGu5fq7P6XE1rTED0y/FfMr7kjMvWfTPtq9SmE+Cr2LG\n7IQCVwBNXI93w/xnfmrIQS3grmL/7vIXjhlzlIn5RxVe8NPFYYIxt1cAu4AHPLxf8a1gzO69mONs\nKdf+J2LGe4r/CMbc1sTMPFIbcwY5E7gO9Qt+bTOn5947iJl771lynsIvixmLkwwcAtYAL2R7vCfm\nU9w2zFRgmZyenuZzzOWB7OphBpzvwFwWmc/psTkhwIvAOldNG4EBJfst5inT9ZWR7VfxH8GY27GY\ny3dHsn2t9PB7iPcFY3ZfwIztPAr8A3yM92c3Ec8KxtzmFlRTloUCy4Dprp+rAwnAesxE21Vtqstp\nLsP8x+ykOfaCnbJbOOXWmZTdwim7zqPcFk65tYm7f+BdMZ9aTo056YkJ8cXAL66fg9W9mE941TCT\njk8j5+BysZeymzfl1vmU3bwpu86m3OZNufUTZwM/Y5ZSPPXJbR1mjAeYZW3X2VCXU/yIudyxDzOR\ndJ2Cny4+pOzmT7l1NmU3f8qucym3+VNu/cR3mLumW3A6xNmn0QjJ9bOIUyi74q+UXfFHyq04WmHD\nG+7ALIG3jPznZ7PIOdWGiBMou+KvlF3xR8qtOF5h8/TegJkG43bMNBWVMavJ7MJcptiJucNwd14v\njoiIsJKTkz1WrASlZIq38laxs6vciocou+KvipNd9Qtit+Iec/OU/XLFYKCH6/uewNv5vMbytddf\nf13vGUDviWfOChQ1uz7/fVpW8PydBst7Bkt2g+Xv0473tOt9PZBd9Qt6T5+/pzu5Lep0Gad2+DYQ\njZmC5JYCQiziFMqu+CtlV/yRciuOU5RliOe4vsBMlNzS8+WIeIWyK/5K2RV/pNyKIwXcxMhRUVF6\nzwB6z2ASLH+nwfKewSJY/j7typCy6z3BkqNgeU935HeHpae4hlmIFE9ISAh4P6e5KbdSYsqu+Ctl\nV/yRO7kNuDO9IiIiIiK5qekVERERkYCnpldEREREAp6aXhEREREJeGp6RURERCTgqekVERERkYCn\npldEREREAp6aXhEREREJeGp6RURERCTgqekVERERkYCnpldEREREAp6aXnGkpPh4erdqZXcZIkWm\n7Iq/UnbFHxUltyFersWyLMvLbyGBJik+nllduzIwOflUQL2d09yUWykWZVf8lbIr/qioudWZXnGc\n2SNHMjA52e4yRIpM2RV/peyKPypqbtX0iuOEpabaXYJIsSi74q+UXfFHRc2tml5xnPSyZe0uQaRY\nlF3xV8qu+KOi5lZNrzhOTJcuxNWta3cZIkUW06ULz1W4kkRa2F2KSJHEdO7McyEXMon77S5FxG0x\nXbrQq159HuIbt57vTtMbDiwClgNrgEGu7f2AbcAy11frohYrkpfINm1oVaMGff7zn5LsRrkVn4ts\n3JgrMmox5JK7SrIbZVd8LnLPHq697BrGXXF9SXaj7IpPRbZpQ6s6tSndYJZbz3f37szywHEgDJgL\nvAzcChwBhhXwOt2NKUW3aBH83//Bxo2ElC4Nxb+LWLkV3+rZE06cgBEjCAkJAWVX/EFGBlx+OXzy\nCbRooeyK//jzT2jdGjZvJiQ8HArJbZibuz3u+rUMEAoccP3s6ylNJBgMGQLdu0OYu/HMl3IrvnP4\nMIwZA0uXemJvyq74zrRpZFSpTmhkpCf2puyKz2x/41PqxnYj1M2xve6O6S2FuVyxC/gNWO3aHgus\nAD4DqhatVJE8bNgAc+bA0097Ym/KrfjM+oHfcUd4AjRo4IndKbviG5bFujcmcMuxaRDikb5U2RXf\n2LqVN+Ov5qcLnnP7JUVNeBVgFtATM15nj2t7f6AekLtT0eUKKZpnn4VateDNNwFKepntFOVWvOvk\nSTLPjyD5w1lcdPflgLIrfmLuXGjfnj2/r6NW3VBA2RU/0b07VkgpGDqUkBD3clvU68eHgHigCZCY\nbfunwPS8XtCvX7+s76OiooiKiiriW0rQ2LWLxHHjSHz2WciWGw9QbsW7vvmGpHp1SFw20dym4znK\nrnjX4MEk3nkniR/19/SelV3xngMHSBwzhsQOHeCNfm6/zJ1PcjWBdOAgUA7zye0NzCWLna7nvAhc\nCzyS67X65Cbu69MH9u6F0aOzNpXgjINyK76Rmcn8C9vRdPQTlGoVnbVZ2RWns1avYULzEbT9Zzhh\nlcplbVd2xemSOozjQPJ+7v61a9Y2T53prQd8iRmnUwoYB/wCfAU0AixgM9CpOIWLAHD0KHz0ESxY\n4Kk9KrfiEyemzuKtfR2ZeONNlPfMLpVd8Ykjb3/A75c8wwPlyxX+ZPcou+J9KSlUmDKe9Hc+KvJL\nvX03pT65iXtGjoSkJJg0KcdmD40tKyrlVtwXFQWdOsHDD+fYrOyKo/37L1xxhbl5uEaNHA8pu+Jo\nY8bA1KkQH59jszfG9Ip4Xno6DBsGEybYXYlI0SxaBFu2wAMP2F2JSNGMGAGPP35GwyviZFZ6BulD\nhlN6zIfFer2WIRb7ffedmebp+hKtBCTic10e3ceGx97wxJzSIj6TefAwNw+7g+0Pv2x3KSJFsmLk\nHFrv+gKKOae0ml6xl2XB4MHwyit2VyJSNBs2cNeezzi3e1u7KxEpklKffsLHraZw1vXn2F2KiPss\ni0bfxTF55L/FnlNapyfEXj//DGlpcNttdlciUjTDhtGyy+VQvYLdlYi47+RJGD6ci6dP1xpp4l/m\nzYM9e6j62B3F3oWaXrHXqbO8pXTRQfzHzpV7qPHNJEqvX134k0Uc5H8DZnHpxddQqXFju0sRKZJJ\n3eZyZ9dXKRsaWux9qNMQ+yxbBuvWnXHXu4jTDe+8lgnXDIbate0uRcR9mZl8/cFBNj7wmt2ViBRJ\nyvJ1zF5dn5B2j5doP5qyTOzzyCNw9dXwcv43U2jqHHGco0fh/POx5s0n5OKL8n2asiuOEx8PvXvD\nH38UOCZS2RXHefppc8N7nz75PkVTlolzbdkCs2aZBSnykJwMhw/7tiQRt4wdCy1aFNjwijjS4MHw\n6qvFvglIxBb//gtTppg5pUtIwxvEHu+9Bx06QOXKeT68fj0sXerjmkQKMTcxne/6r9NsI+J33nl+\nK9PWXqQ5pcXvPN16O3+0es0jc0rrTK/43r59MG4crFqV71NOTebQoYOPahJxQ+WFszlZr4zmlBa/\n0zb5HcJfaKY5pcW/HD5M97+7EvHdNx7Znc70iu99+CHcey+cddYZDy1fbqbuFXEcy+I/E+K4ZVC0\n3ZWIFM2GDUQsnUj9lx6yuxKRovnkExq2aUD4Jed5ZHdqesW3TpyADz7I8+a1tDTo0QP27rWhLpEC\nWBakzPxVc0qL39m9G/YM+BiefRYqaE5p8R/bN5/k0LDPPDqcTE2v+NaXX8J118Fll53xUOnS5t62\nWrVsqEukAEuXwl3tqmpOafE7v3x/iNETq8MLL9hdikiR/PD6H3xZ6QXw4JzSGtwjvpORAUOHwhdf\n5NickgKpqVClij1liRSmSak/mFL2IXhYi1GIf3l4+1B4cq/mlBb/kpnJc388A++/59HdqukV35ky\nxRx4mzfPsfn7782ZtHfftakukcIMHUqF7p2gTBm7KxFx39GjZlrIBQvsrkSkaH780Vz+bdnSo7vV\ndTrxDcvKd47IRx6Bd96xqS6RAmRmwphBe0n76Rfo2NHuckTctnIl9P+/ldCiBVx4od3liLjtr79g\ncOw/XplTWmd6xTeSkuDQIbjrrqxNqalQtqz5XrPoiBMdPQrJk5YR9syT+c4pLeJEtaulcePSkfCD\n5pQW/1Jx/R9ceWwhPPCpx/etVkN8Y/BgM2OD6yagjAxo1gymT4f69W2uTSQfldP28famB2G6xvKK\nf6mT9B11Lt2hOaXF79T/8i3qx0V65WyYt9ci1FraYhahiI6GzZshPDxr8/79UL16wS/VGvBiF8uC\nkAH9zZLZn31W5Ncru2KXXTst6tx2NQwcCLffXuTXK7til2N/JlPh1qbmuFvEKfbcyW1hY3rDgUXA\ncmANMMi1vTqQAKwHZgNVi1SZBJehQyE2NqvhPXVcK6zhLSFlV4rNsuDmFhlsHjEtzzmlvUzZlWL7\n+2+4++ZDWCdtmVNa2ZViO3gQGkdWJK3Dc16bU7qwpjcFuBloBPzH9f2NQE9MgC8GfnH9LJIlKT6e\n3q1a0a9ZM3p//TVJEREA/PCDWYDCB5RdKZak+Hj6tG5Fo6238mnI3yRt2uTrEpRdKZak+HjGdGxF\nqx0X0CcEkmbO9HUJyq4U2al+YfjtN3Dv0QtZ0PBiu0sCoDzwP6AhsA6o49pe1/VzXiwJPnNmzLB6\nRURYljlhZllg9YqIsObMmGEdP25Zmze7vy/AE9e7ippdr/3ZiLMVlN2iUnbFl5Rd8Ue+zq07U5aV\nwlyq2AX8Bqx2hXeX6/Fd2cIswuyRIxmYnJxj28DkZBJGjaJcOWjQwGelKLtSJLNHjqRV8lkcyHb1\n9VR2fUzZlSKZOXw0J5KfIy3b/enKrjjdqX5hLO2zjrvezK07t8ZlYi5VVAFmYS5XZFdgd92vX7+s\n76OiooiKiipqjeJnwlJTc/xsAZ34mAoHfy70tYmJiSQmJnqqlGJnV7kNTmGpqczkdmqwj2oczNoe\nmpJS6GuVXbFTaEoajVhOadJzbVd2xbnCUlPJJIQtNKA0aVnbbchtnvoAL2MuTdR1bauHhjdINnEx\nMTkuVVhgzSLa6hV9W5H3hWcusxU1u174UxF/kFd2LbB6t2pV5H0pu+JLyq74I1/ntrDhDTU5fZdl\nOSAaWAZMA55wbX8CmFrYG0nwiOnShbgqVXJsS4zYRKuuz/uyDGVXisSyIKZzZ+JyzQ3ZKyKC6NhY\nX5ai7EqR7NwJMTExxIWG5tiu7IrTtejcjbiqOSfz8GZuC5uH70rgS0xzXAoYBwzBTD8yETgX2AI8\nCNmuBZ7mar4lqGzaRNJVV/Hfy9vy78HLuer8X4iOjSWyTZsi76oE80WWJLvKbRB65hl4tNz3hC54\ni4SaNQlNSSEjPFzZFUc7ehSaNbNYXCWG/119GQnr1yu74hcsC6696iSvJDdm1bW1CAWv51aLU4jn\nPfooXHIJyY/2ZeVKuOee4u9Kk6SLr2xfd4RqkVdSPuEHuOqqEu9P2RVfyZg+k9AeL8Off3pkFStl\nV3zlSLvnqVS/MgwaVPiTC6GmV3xv2TKzAtCGDVCxYol3p4Ov+Mzrr8OmTTBunEd2p+yKT2RkQOPG\n8OabJTvDkI2yKz6xdi20aAHr10PVkq9X4okV2USK5rXXmNTmcw5nlrzhFfGFhQth+/I98P770L+/\n3eWIuCU93Sx0eXTsRHOC4e677S5JxG2jR8PO7oPh1Vc90vC6q+TXQURO+fVXrPUb+OOqaG44CpUr\n212QSOEWLoSUOf+lfrt2Pp1EWqQkLAuaNUmjQt/XYPw4CPH1iVmR4rEsOL5mC+VXLoIpo3363hre\nIJ5hWXD99dC9Ozz0kMd2q8ts4nUbN0LTprBuHdSs6bHdKrvidcOHwy+/wPTpHt2tsiteZVlmWEP7\n9ubLQ9zJrc70ikcc+3oqGw5fQKMHH7S7FBG3WJbr5FifPtCtm0cbXhFv2roVzq1yyNz883Phi/6I\nOMWRI1ApaSbs2wft2vn8/dX0SsmlpbHytf8y8bqhNCqlYeLiH3r1gqsqbeKhOXPg00/tLkfELWlp\ncO+98PPN71PtttvgyivtLknEbbffZvHxrk+4/N1BkGteaV/Q8AYpuY8+gsmTISHB47vWZTbxlgMH\ngHvvpdr/xcCzz3p8/8queIv17w5CrrzCzJZz7rke37+yK96SMmYc4V9+DL//7vFx6JqyTLzun7+O\nc3bUhYTMmA7XXOPx/evgK16TkADPPQdr1kDp0h7fvbIrXvPcc1CuHLz7rld2r+yKV6SkwCWXwH//\nC82be3z3mrJMvO7Zu7fz538e80rDK+INy5fD+nWZ0LMnDBzolYZXxNMyM02vu2fRJpg40YzPEfET\nI0bAlgHjoVEjrzS87tKYXim+vXuZtucGSs1YYHclIm5bswbCl87n4lKloG1bu8sRcVvLllDj3V5m\nlpwaNewuR8RtFUodp+pHb8OcH2ytQ8MbpMgsC44fhwp9ukNqKnzwgdfeS5fZxONOnoTLL4ePP4Zb\nb/Xa2yi74nFLlphFKNavhwoVvPY2yq54XFwc7NgBY8d67S00ZZl4xYIF8M7rx/nhjy9h9Wq7yxFx\nS9YUZZ9+ChERXm14RTxp40aIuMAipEcP6NvXqw2viCft2wc1Tu4wN7wvX253OTrTK8Vz9JGOVLyw\nrlnv3Yt0xkE8ZfhwsFJSeXFEA5g5Exo39ur7KbviCZmZcMstML7DHM5+syOsWuX1cejKrnjK7bfD\nG+GDuDZiPwwZ4tX30uwN4lGZmVCqFPDnnxAdDRs2eH2tYR18xVOOH4ej/YZSe9sf5u5hL1N2xVOs\njExCmlxjLhH7YBy6siuekrbqL0pHNTdDcqpX9+p7afYG8agHH4TFizF3Dffq5fWGV8STyh/bQ+2x\nb8OAAXaXIlIkIRMnQJkycP/9dpciUiSl3+gNL7/s9YbXXTrTK27buhXO3pREqfZPwLp1ULas199T\nZxykpDZsgIMH4drxXc3lilGjfPK+yq6UhGWZKcp6dj/Jea0vM2PRb77ZJ++t7EpJDRsGt9X/k8u6\n32YOwuXLe/09dSObeNS551jwfz2gf3+fNLwinrBlC2xduodrx4+HtWvtLkfEbffcA/VnjoGLL/ZZ\nwyviCeeeY1F7RB/o188nDa+7dKZXCrRqFYwcCZ98AkyZAm+8AX/84Rrc63064yAe8dhjcOGF5gDs\nI8qulNiRI3DRRTBrFlx1lc/eVtmVEvvpJ+jWzTQRYb45v6oxvVJiF10ETz8NpKebcbyDBvms4RUp\nCcsyXyxfDj//DC+9ZHdJIm5ZuxYyMjDXiFu29GnDK1ISO3aYGy/p2RPeestnDa+73OlezgF+A1YD\nq4Auru39gG3AMtdXay/UJzYrWxauvx744guoWxda+81fs3Ib5L77ztw/wWuvmbveK1WyuyR3KbtB\nLi4ONizYay6z9e9vdzlFoewGuZdegjl9f4HwcLj3XrvLOYM7ly/qur6WAxWBpcA9wIPAEWBYAa/V\n5Qo/NWwYtGoFDRti5nq6+GKYPNnVAftOCS6zKbdBLj0d9kydR71XHzc3XpYp49P3V3alRGJjzVW1\nESN8/tbKrhRX5olUQi67lJAvv4AWLXz63p66kW2n6wvgKLAWqH/qPYpbnDjbBRdA7dquH0aNgqZN\nfd7wlpByG+TCQi3qDX4RBg70ecNbQspusEtOhm++8ccbL5XdIFdqzMfQ8HKfN7zuKurgzAZAY2Ch\n6+dYYAXwGVDVc2WJ3aqXjmfEY63o17w5vfv0ISk62u6SSqIBym3QmPj5rzzapCf9rriC3uvWkVSx\not0llUQDlN2gkBQfz3/OSiD2ukfpfeONJLVqBbVq2V1WSTRA2Q0KSfHxRF30BZ0aPUXvl18m6bbb\n7C4pX0UZYVwRmAR0xXyCGw2cWoO2P/Au8HTuF/XLdrd0VFQUUVFRxatUvO7QIZg0CS6sE8/sbl0Z\nmJyc9VjckCFw9tlEtmnj1RoSExNJTEz05C6V2yCSFB/PzD5juWl7LTqzBoC4F1+EUqWUXXGspPh4\nZnXtyic7atJkxxLCyCBuwQKIj/d6bkHZleI7ld0+yedxNX9QjTTihg+H88/3x2NultLALKBbPo83\nAFbmsd0S//HPP5Y1eLBlxcXEnLrxPcdX71atfF4TUJJBXsptkFF2lV1/5KTcWpayK+5zUnbdya07\nwxtCMJcj1gDDs22vl+37e/MJsfiRs8+GV16BsNTUPB8PTUnxcUUlotwGmZQUKJWSludjyq442eZ/\nr+Ikpc/Y7me5BWU36OzYfSHphJ6x3anZdWd4Q3PgMeBPzFQjAL2Ah4FGmM56M9DJGwWK9y1aBJdc\nAlVdo6zS81ltLSM83IdVlZhyG2Q+/RQStz6BmTEpJ2VXnCojA1btb8kxxlCGgzkf86/cgrIbVCwL\nFu5sy26mcBY7cjzm1Oy60/TOJe8zwj96uBaxSXw8ZGZCs2bm55izziKuTBkGnjyZ9ZxeERG0jo21\nqcJiUW6DzPPPpHLxByOJq1KFgYcOZW1XdsXJQkNhVNxGhsYeZmDm6e1+mFtQdoNKSAiMvutbPvh8\nLwOzXWRzcna1DLHkNHEivPQSSf37k/Dtt4SmpJARHk50bKxPbqjITcthilsyM81SwydPkvTEEyR8\n8IGyK46Xng5he3ZAs2YktW1LwqpVtucWlF0pXEYGhMZPg86dSRowgISJE23Prju5VdMbpCwLNmww\na05kmTfPrKAyezY0amRbbdnp4CuF+fprqBM/lpabx8Cvv0K5cnaXBCi7UrCjR+H6azNYEn4T5e67\nDfr0sbukLMquFOTIEWh+9XEWH7iY8B+nwLXX2l0S4F5uizpPrwSITZugSxfT/AKmA77/fvjqK8c0\nvCLuOH/lNM6aOxGmTXNMwytSmIrlMphzbjvKNb4Ueve2uxwRt1Xat4WEw00J/+wDxzS87tKZ3iBm\nWWZMDnv3mgG9r7wCHTvaXVYOOuMgBZo9Gx5/HH7/PddlC/spu5IvyzLLDP/1F8ycCaXPnLnBTsqu\n5OvgQbjhBujc2Zw5cxCd6ZUzHDp0+uxuSAhmnqe774a2bR3X8IoUZOOMdaQ+0t6sqOKwhlckP8eP\nw4A288n8bY7JrsMaXpH87N95koFXT8aKjnFcw+suNb1B5tVXYepU1w+ZmfDEE3DOOTBwoK11iRTJ\n9u28//A8Ep8ZDzfdZHc1Im5LnTSdGvNnUOrHeKhSxe5yRNxjWWR07U79qscIGfau3dUUm4Y3BJnU\nVChTxnWWt2dPmDsXfv4ZHDqnni6zyRmOHIHISHjwQXjtNburyZeyK2dYtAjuuANmzYKrr7a7mnwp\nu3KGN94w85smJkL58nZXkyd3cuvOPL0SADIyzHyQWetOfPwxfP89zJ/v2IZX5Azp6fB//2dunujZ\n0+5qRNxmJW/i3zs7U//zzx3d8IrkdvyT8ez/bBZn/2+6Yxted2l4QxD45x9o3tyMZgDgxx/h9dfN\nDRQ1a9pam4jbXDf/PLLsFf7s9IHrcoWIH9i/n5W3duOZWlPNmV4Rf/Hrr8zr8QMDmsVDnTp2V1Ni\nGt4QJHbuhLp1geXLIToafvjB3IHpcLrMJlmGDIHx41kzZh4XX12RMIdfp1J2BTBjymJi4NpryRw8\nlFJ+cKpJ2RUA1qyBm2+GCROwWkQ5/jyDFqeQnLZtM1OTDRsGDzxgdzVu0cFXAPjuO+jeHRYsgLPP\ntrsatyi7gmWZlQJTU81ql/7Q8aLsCuZMWbNm0L+/ybAf0JRlQa5vX5g+3fXD4cPQpo2ZZsRPGl4R\nwIw7f+45JnVJ4mhV/2h4RQDo25d5KyrS6/xv/KbhFeHYMbjzTt5v+CEfHvaPhtddDr9AKCXx5JNQ\nrRqQlmbudG/eHF5+2e6yRNy3cSPcdx+ZX3zF3ITzaZlud0Eibho7Fr75hitnL6TcQc3FK34iIwMe\neQSuuIJH323N8RN2F+RZGt4QoJLi45k9ciRhKSmkJycTU68ekQsW4PiBkLnoMltwycptairpISHE\n/PUXkf36+eXCKcpucMmR3aNHidm4kcjFi/1y4RRlN7jk6Bc2byamWjUi//c/M7+pH9GUZUFo/354\n6vG/uWxddwZtWp+1Pa5MGZg1i8g2bWysTiR/SfHxzOralYHJyVnb4qpWJbXmOUTbWJdIYfLK7gu1\nm7D55x2w7/wpAAAgAElEQVQ84YdNrwSPvLL78Mmb2fblPB7pcLONlXmHBhkFmIoVIWz7FzkaXoCB\nmzeTMGqUPUWJuGH2yJE5DrwAvQ+m8HC7KzhyxKaiRNyQV3Zv312Lr0eutakiEffkld1bd5Xm188n\n21SRd6npDTBlysAVVX/L87HQlBQfVyPivrDU1DO2lSOFjo07UqmSDQWJuCmv7N7Oj9xQ91sbqhFx\nX17ZfYbPOLvMKhuq8T41vQEiPh6WLDHfp2ctu5ZThlZeEwfLL7ehFTTOT5wte3Yzsv23qmOuOF32\n7O6nWtb3gZpdNb0BImu1NSDmmmuIK53zbuFeERFEx8b6uCoR98V06UJcxYoAHKM8LzGUnhdcrNyK\n48U89RRxoaEAvE1PRvGCjrniF2KaNSPOdYP7/UxmFQ0DOrvu3J15DvAVUBuwgE+AkUB1YAJwHrAF\neBA4mOu1uhvT1w4ehCuuIKlLFxJ+/ZXQlBQywsOJjo31y5vYSnAXsXLrb2bNIunxx0m46ioyUkJZ\nc+gWug9q6Je5BWU3qHTqRNI//5CQmYl1PJPU0pW4s/tTyq6y62xHj5p+4emnSZg3D46nQfnSAd0v\nuBPquq6v5UBFYClwD9Ae2AsMBnoA1YCeuV6rEHtRZiZMnQr33svp5QE7dzarAH38sa21eUoJDr7K\nrT85dgyuuAI++ghatbK7Go9QdoNEUpKZ13T1aqhSxe5qPELZDRIvvgj79sFXX9ldiUd4akW2nZgA\nAxwF1gL1gbuAL13bv8QEW3zo0CGYPdusPQHA77+bJdjeecfWuhxCufUnfftC8+Ycv6kVK1faXYzt\nlF1/kZJi5pAeNYr3x1VhzRq7C7Kdsusv/vc/+OYbDr7+HnfcYVbKDgZFHdPbAGgMLALqALtc23e5\nfhYfqlbNnBgrUwaT2I4dYeRIqFrV7tKcpgHKrXMtWQLjx8N777FyZcBcpPCUBii7zvXWW3D55XDv\nvdSpAzVq2F2QozRA2XWmtDR45hkYOpQqF9Sgb1/I5z7igFOUxSkqApOBrkDuWTMt19cZ+vXrl/V9\nVFQUUVFRRSpQzvT331CpElSvnm3jW2/BJZfAfffZVpcnJCYmkpiY6MldKrdOlpYGHTrAkCFQqxbX\n14Lrr7e7qOJRdoPM6tUwejSsWAHAAw/YXE8JKLtB5t13oV49ePRRQkLguuvsLqh4ipNbd8fslAZm\nAD8Cw13b1gFRmMsZ9YDfgEtzvU5jdLxgxAjT9D71lGvD6tUQFQXLl0P9+naW5nElXA5TuXW6wYPh\n559J+WEWZcNDTo9NDwDKbgDLzIQbb4R27Zh3ZWduuAFl9zRl18k2boSmTVk/YRnL953Dgw/aXZDn\neGpMbwjwGbCG0wEGmAY84fr+CWBq0UuU4ujaNVvDm5lphjW88UbANbwlpNw6XXKyaXo/+ojX+4Uw\ndqzdBTmGsut0o0dDaCgp7ToyeDBaMfA0ZdfJLAs6dYLXXiO19jlkZNhdkO+580nuRiAJ+JPTlyRe\nAxYDE4Fz0RQkPrF7N9SunWvj6NFmPOTvv0OpwJt2uQRnHJRbJ7MsiIkxX6+8Qmqq2RRI86EruwFq\n2zZo3NjM2nDZZXZX4xXKboD64gt4/31YuBDCijK61T94asqyklCIPeTwYbjpJpPVcuVcG7dvh0aN\nIDERGja0szyvKeFltuJSbr3tq69g+HBYvDggD76g7AYky4J77oGrryalx+sB9SEtO2U3AO3eDVde\nyV8f/ca5rS8/3UcEEE8NbxAHqFwZli4lZ1BfeAGefTZgG14JUHv2wCuvYH0yho7PhbF1q90Fibhp\n8mTYsIEl0a9xjybdEn/SrRu0a8eYeZeTlGR3MfbRmV6HS0uD0NA8Ri58/z306mVuXgvU0w3ojENA\neuwxqFsXhg5l1iy49dbAPNmr7AaYAwfMAioTJsCNN3LoUMCsRXEGZTfA/PgjPP88rFoF5cvbXY3X\naHhDAHjnHdMQvPRSto2upYb55hsz5iGA6eAbYGbNMlcnVq6EChXsrsarlN0A07GjOQMxerTdlXid\nshtAXEsNn/zwU8rc3tLuarxKTW8ASE2F9PRc/UGALTVcEB18A4hrqeGT73/CpEPRPPxwYE3zlJuy\nG0BcSw2P77ue0ErlefhhuwvyLmU3gLz4Iim7D3PN8s+YPz9wr06Ae7kNwIuKgcGyTENQtmyulVJO\nLTW8erVttYkUS9++cOON7Ls6mhXDCfjGQQJEtqWGm15Z/vSy7yJO51pqOHzVKuaGBnbD6y6d6XWg\n48fNTE4//mgWociSmmpmaxgwAO6/37b6fElnHALEkiVwxx1mWEOtWnZX4xPKboDo29eMhfz+e7sr\n8RllNwCkpUGTJvDKK+Y+iiCg2Rv8VPnyMGZMroYXYNCggFhqWIJMejp06MCxAe9xuGxwNLwSIFxL\nDc9qO4YDB+wuRqQIhg1jeYXm9PjzUbsrcRQNb3CoM+Y8X7MGPvjAzNYQyAMhJfC89x7UqsXU8IdY\n3h+GDLG7IBE3ZGZChw7Qvz9zVtUg4jqoVs3uokTcsHEjDBnCBb8u4a4j6hey0/AGB/n2WzhxAtq3\nz/VAZqaZpeHRR+G552ypzS66zObnNm2C664zi1BccAEZGeYG+GCg7Pq5Dz4wB+U5cwJytcuCKLt+\nzLIgOhpuuy3XtE+BTzey+ZmmTc143jOcmqWhc2ef1iNSIq513q1XexBywQVA8DS84ue2bYN+/dj7\nwzxqBlnDK37uyy9J+ud8arfqyqV21+JA+tfsIA0awOWX59q4fbu5keKTT4LubIP4uXHjYN8+Oqx7\nKahXABI/Y1nw/PMc6dCdmztdzIkTdhck4qbdu6FHD/5+LI5d+3ROMy/6U7FZRoaZjKF795w3riXF\nxzN75EjCli4lvUoVYrZsIVLLDYvDZeX2yBHSlywhZvBget9Tinr17K5MpGBZ2d22jfS//ybmySf5\n4w0oXdruykQKlpXdZctIL1+emKtXE9migd1lOZKaXptZFtSsCeXKnd6WFB/PrK5dGZicbDbs20dc\n164ARLZpY0OVIoU7I7dA3Pvv0+qii2jQQLkV58ozu6+8AmXK6JgrjpY9u1s4jwb8rX6hALpebrOw\nMLMkdli2jx+zR47McfAFGJicTMKoUT6uTsR92XObQSnepA+9kncot+J42bO7kOt5lg91zBW/cCq7\nGZTiMcbzL/WU3QLoTK9N1q2DAwegWbMzHwvL8242CE1J8XJVIsUXlpqa9X0GoVTjAOGkKLfieNmz\n25hlVOIIoGOuOF+Ya9B5KJn8zk1ZUxcou3nTmV6bbN9uptI7Q2Ym6Rs25PmajPBw7xYlUgLp2b4v\nQxqxvE8omcqtOF76wYNZ35flJA1ZA+iYKw5nWaRv3MhJzMDz7HN1Kbt5U9Nrk1tvhccfz+OB/v2J\nqVaNONcUT6f0ioggOjbWN8WJFNWJE8T8+y8vVr2Q/9Eka7NyK463fDkxmzZxd/UX+IPGWZuVXXG8\nIUOoTQuuqjA1x2ZlN39anMLHFiww8/Hmuajad9/Byy/DokUkLV1KwqhRhKakkBEeTnRsbFAOStck\n6X7AsszCKZbFqKtimfzZBqLqfx7UuQVl1y/s3AnXXw9DhjDwzwv4d86H1ArdpOwqu843bRo8+yzW\nwkXEz1vDoi+GqV9wI7dqen0oPR0efNCsNVGrVq4Hly6F1q0hIQEaNbKlPifSwdcPvPUWTJ1qVq7K\nPg1JkFN2HS4lBW6+GWJi4I037K7GUZRdh1u50lwunj7dfGgTwL3cujO8YSywC1iZbVs/YBuwzPXV\nulgVBpmwMPj++zwa3h074J57TDeshteTlF1vmzIFRo/mr6HTscLV8HqIcuttrtUCM+qfy8TLXke9\nlscou962Zw/cdRdPXrqQJaFqeIvKnab3c84MqQUMAxq7vn7ycF0B5fhx2Ls3nwdPnDANb6dOcN99\nPq0rCCi73rR8OXTsiPX9FF58uw6bNtldUMBQbr1tyBBYuZID745l3oJSZGTYXVDAUHa96eRJuP9+\neOQRXv3oAq66yu6C/I87U5b9DjTIY7uvL334rRkzYPFiGDo01wOWBc88AxdcAHFxttQW4JRdb9m1\ny3xYe/99Qq5tQnx8PuPUpTiUW2+aPh1GjICFC6l5TgVGjLC7oICi7HqLZcGzz0KNGtC/P5drGoJi\nKckfWyywAvgMqOqZcgLTgw/C4MF5PPD227B+PYwdq47Bt5TdkkhNhfvuI7Pdkxxq/X+A4usjym1J\nrVoFTz3F1g9ncKDiOXZXE0yU3ZIaMYKt87fx8jkTsELU8RZXcRenGA286fq+P/Au8HReT+zXr1/W\n91FRUURFRRXzLf3PsWNQoYL5vlTujP7wA3z4ISxapJt/sklMTCQxMdGbb+FWdoM5twWyLOjYEerV\n4+cb+vLV8zB+vN1FOYOXs6tjbkm5xkLy3nt8v7kxlffBU0/ZXZQzKLsO99NPMHgw1X9eSNSmMjrJ\n4FKc3Lr7R9cAmA5cWcTHgvZuzMxMaNLEXEmrXz/XgytWQMuWMHMmXHutLfX5Cw/cRdyAomc3aHNb\nqKFD4euvYe5cqFCB1FQoW9buopyphNltgI65nnPyJERHww03wKBBdlfjeMqug6xdCy1amJuGmze3\nuxpH89TsDXmpl+37e8l5p6ZgzuwmJeXR8O7eDXffDaNGqeG1h7JbXDNmwHvvkf79tKxLGGp4fUa5\nLS7Lgueeg6pV2dd9oN3VBCNlt7j27YM772TqQ9/yV001vJ7gzvCGb4AWQE3gH+B1IApohLkrczPQ\nyUv1+SXLMmMcK1bM9YBrLCTt2sFDD9lSW5BRdj1l1Spo3x6mT6d1h3MYMQIaNrS7qICl3HrSyJGw\neDHHZs+jxc2lWLgwj2OzeIqy6ylpaeaGoHvv5WDDW0hPL/wlUjgtTuFhU6fCb79x5h3BlmUGkB0+\nbFZeO2OQr+RFk6Q7wN69cN11ZgL/xx9n1y6oXVs3rxVG2XWAn34yH9YWLIAGDTh5EsqUsbso51N2\nHeC55+Dvv83Ka6GhdlfjF9zJbXFvZJN8tGmTz6iFYcPMvKZz56rhFf9xal7IBx+Exx8HoE4dm2sS\ncce6ddCuHZmTvifkvAaEoIZX/MSHH5KZmMSvgxbRUg2vR6n78rDSpfMYxztzJrz7rpmx4dR0DiJO\nZ1nw/PNQpQoL73qLTrooKf5i/3648054+20+Xn0j/fvbXZCIm375Bd58k92fxzPu+wpaOMXDNLzB\nAyzLnAQbMAAaNICk+HhmjxxJWGoq6WlpxKxaReRPP0GzZnaX6nd0mc23cmR31y5iTpwgcuVKUstU\nYuNGjeMtCmXXt7Kye+IE6atWEdOiBZFTppCaaqaPrF7d7gr9h7LrOzmOuenppl+YOhU0XVuRaXiD\nj4SEmIVSzjnHBHhW164MTE7Oejyudm3Yv59IG2sUKUye2T3vPEhKIrJNGzW84lh5ZnflSoiPJ7JN\nG80yIo6UV2671mhIyr4TxNhYVyDT8AYPad7cjDWfPXJkjgADDNy9m4RRo2yqTMQ9eWWXvzvwzaAp\n9hQk4qbc2c2gFGuShzBlqFZOEefK65hbe9/djHz9L5sqCnxqektgxYozZ2kIS03N87mhKSk+qEik\n+PLK7o3MpVrIdhuqEXFf7uyGkklP3qZy5g6bKhIpXF7H3F68xTU1ptpQTXBQ01sCNWrAJZfk3Jae\nz3W0jPBwH1QkUnx5Zfc2fiK0QvCNsxP/kld2r2cxmeV03BXnyp7bY5QHzIBU5dZ71PSWwNlnQ+vW\nObfFdOlCXOXKObb1ioggOjbWh5WJFF1Mp07EhYWxg7qMpjOg7Ip/iLn1VuJCQ1nCNcqu+I2YLl2I\nq1KFI1TkWv7Hccopt16m2RuK4auvzHy8NWrk8eCKFSS1aEHCNdcQmpFBRng40bGxRLZp4/M6A4Hu\nIvahN98kafZsJoedx+Ydjbjq/F+U3RJQdn0kIwOaNCGpZUu+W7ybgwdq0uCs1cpuCSi7PrJ2LUnX\nX09CkyZkpJUltIKl3JaAO7lV01tElgVvvQUdOphVqc54sEULePRRNKmpZ+jg6yN//w3XXANLl8J5\n59ldTUBQdn3k44/h669hzhwtE+ghyq4PWBYnou8ivM2thLzYze5qAoI7udXwhiIKCYG4uDwaXoBv\nv4WjR+GZZ3xel0hJpHd/ldevmMyR6mp4xY/s38/xPoPofu4kjp9Qwyt+5Icf6Lb0cabUf8HuSoKK\nml43HT5slm/P19Gj8OqrMGqU1skW//Lrr4QsXUK9+5qh+y3Fr/TtS9n72tAopjblytldjIibTpyA\n7t1596ta3NNWyyX4koY3uGnJEpg8GQYNyucJcXHmEvF4zQvpSbrM5mXp6dCoEbz5Jtx3n93VBBRl\n18v+/BNatoS1a/O5wUKKS9n1rtR+gyi70tVUiMdoTK+vbNwITZuag/BZZ9ldTUDRwde7Rtw3h2bb\nvuO6RaM0HtLDlF0vsizeuuBTbnu4Ko3fesDuagKOsus987/fyTsPLeOH9ZdBgwZ2lxNQNKbXAzZs\ngMzMQp7UvTu88ooaXvEve/bQ8JeRnD24ixpe8S8TJ3JdqSWc201XJ8S/3DChK192X6GG1yY601uI\ntm2hb1/4z3/yecKPP0LXrrByJVrg3fN0xsGLOnaEChXgvffsriQgKbtecuwYXHaZmbHhppvsriYg\nKbued/IklJmfCE8+CWvWQPny9hYUgNzJrUZQF+K77wo4CXbyJHTrBsOHq+EVv5GUBAsnbePV6dNh\n3Tq7yxFx26+/wrJBSbx0441qeMVvZGRAs6YW00/056x331XDayM1vXmwLDh0CKpWLeSq74gRcNFF\ncPvtPqtNpKQuvTiTsr/1g4EDoUoVu8sRcdvlFf6m4uJhsPpzu0sRcVtoKCT836dUn41uGLaZhjfk\nYeFC0w9Mn17Ak/7914x5WLDANL7iFbrM5gVffQXvv2+CXkrD+r1F2fWCe+6B66+H116zu5KApux6\nRnq6aXhD9u6Bhg3ht9/Mr+IVnrqRbSywC1iZbVt1IAFYD8wGqhavRGdq2tQMayhQz55mWTY1vE4W\ndNnNz969Znx62r7DJrujRqnhdS7lNpv9++Hhh+HkjNmwapW5cVicStnN5p13zKGW3r3hkUfU8DqA\nO5/kbgKOAl8BV7q2DQb2un7tAVQDeubxWr/65JaR4ea6EvPnw4MPmvGQFSt6va5gVsIzDsXNrl/l\n1h2WBXPnwk3TX4U9e+BzXR72thJkN2iOue6wLPh1Vhq3drsShg6FO+6wu6SAp+x6xvHjYC1bToX7\nW5t+oWrQ9Pu28NSZ3t+BA7m23QV86fr+S+CeohbnRPffD4sWFfKkjAyIjYXBg9XwOl/QZLcwISFw\nU+2/TLP79tt2lyMFU26zCQmBW1ePhAsugDZt7C5HCqbsYoY1AJQvZ1Hh1efNeEk1vI5Q3OubdTCX\nMHD9Wscz5djr44/huusKedLYsVCunLneJv4oILObn0cfNVeEsSwztd5rr0GdgP4tB6qgyi3Ac8/B\n8uXAzp1mKczhwzWftH8Kquxu3gyRkeaQy9dfm1me2re3uyxx8cTsDZbrK0/9+vXL+j4qKoqoqCgP\nvKV3FNoLHDgAffqYuXl18PWKxMREEhMTffV2+WbXn3JbkJ49zZSmTJ9ulsl+4QW7SwpYPsxuwBxz\nC/LUU67sduoJTz8NF19sd0kBS9n1nPPPhylTIOToEejRwyw1rPsnvKI4uXW3c2sATOf0GJ11QBSw\nE6gH/AZcmsfrHD9GZ9UqM9h83Dg3ntylC6SlwejRXq9LDA/cRdyAomfX8bktyKnSsz6XpaSYGyg+\n+giio22rK9iUMLsNCNBjbpEsXGjGna1bB5Uq2V1N0FB2i+6Me4J69IBdu+CLL+wqKeh4cxniacAT\nru+fAKYWcz+2u/RSePllN564ciV8+y307+/1msSrAia7+fngAzOFdJZhw8z0emp4/VnA5xbMIfad\nd1w/ZGaa+yfeeUcNr38LiuzeeadrOA7A+vXw2WdmWI44ijuf5L4BWgA1MeNx+gI/ABOBc4EtwIPA\nwTxe69ef3LJYFtxyi5nz6fnn7a4mqJTwjENxs+vXud23z/xaowbwzz/QuDEsXmxuBBKfKUF2g/aY\nu38/HDzoiupnn5mvefM0nMzHlN2i27ED6taFECxzw+Utt7h5Rk08xZ3cBu3iFAMHmlw2auTGk7/7\nDgYMgKVLIUyL2PmSJkl3z/HjcOwY1KqV64GHHzZzSb/5pi11BTNl1z2WBUeP5jqZe/CguQw3cyZc\nfbVttQUrZdc9Bw6YSZxKl862ccYM0+z++SeUKWNbbcHIm8Mb/N5118G557rxxOPHTYBHjlTDK471\nzTfw4Ye5NiYlmbNkPfOaElPEGWbMMBOL5NCvH9x9txpecbQ33oBJk7JtSE2Fbt3MTCNqeB0paM/0\nuq1vXzM+59tv7a4kKOmMg3ssy3xl3SScng7XXANxcWYhFfE5Zdc9lmXutSxXzrVh9Wq4+Wbz6xmX\nLsQXlF33pKWZc2FZo28GDTI3X/7wg611BSt3chtUpy7374fx4829EQUNEUuKj2f2yJGEHTxI+h9/\nEPPJJ0T6rkwRt6xbB1u3QkyMyXNISLbsbtpE+r59xJQvr+yK4xw+bO4Nbt7c5LZcuWzZXbyY9Lp1\niVm8mEgtRiEO888/ZurdiAgzrCErt4cPk75kCTGjR+uY62BB1fSePJnrU1kekuLjmdW1KwOTk7O2\nxQ0cCLVr6wAsjnLwIOzeffrnPLPbrRuEhCi74igbNkB8vGl6IY/sHjxInGvMg7IrTjJnjrl/IiIi\nn2Pu229DvXrKrUNpeEMuvVu1YsDs2Wds79OqFf1/+smGioKbLrO5T9l1FmXXfcqusyi77lFunUU3\nsrkkJcGePe49NyyfJ4ampHiwIpHimT3brCR8hsxMwjZtyvM1yq44wbp1ZkHLMxw4QFjWBKc5Kbvi\nBOvWmZuFcwvbuzfP5yu3zhUUTe+8ebBlSyFPSkmBHj1IX706z4czwsM9XpdIUTVtCk88kWvj5s3Q\nsiXp2cc6ZKPsihOcdRY0a5Zr44wZcOWVpJctm+drlF1xgpCQ0ytdAmaWhj59SF+5Ms/nK7fOFRRN\n72uvwbXXFvCEhQvNBP7JycSMHUtcRESOh3tFRBAdG+vdIkXyYVnmJkyAypXN9KWAWbHqgw9MuG+/\nnZhx45RdcZy0NPNr5cpw++2ujQcOQLt2Zq6y8eOJGT1a2RXHuuQSeOQR1w9Ll0KTJrBiBTGffabc\n+pmAvZHNsszc0FddVcCTTpyA11+Hr76CUaPggQfMXZdVq9Jn1ChCU1LICA+ndWysBqWLbRISzKwj\nX32VbeOmTfDUU+aMw9y5cOmlJruhocquOMbq1Wa2nF9+yXYD8fTp0Lkz3H+/OUhXqJB1t7uyK06x\ndi289x588olrQ2oq9O8PY8aYpd0feYTIkBCoXl259SMBeyPb5s3w4oswZUo+szUsXAhPPglXXmnO\nltWu7esSxQ26ocI4edI113lmplmFol8/cwmjWzcIDbW7PMmDsmvs2eOabnf/fnNmd/58GDsWWrSw\nuzTJh7JrjrlLl7qG5CxZAu3bm/WxP/oI6tWzuzzJg5YhzsuJE2bBiXHjss7uinMF68H3wAFznI2O\nzrYxORmeftocjT//3FxzE8cK1uxmZpoTuDmWeJ82DZ59Ftq2hbfeggoVbKtPChes2QUzHCdrWeHU\nVLOE+6efZp3dLXDOU7FVUM7esHcvZGTk8+CCBWbs7tatZmZ0NbziUNu3mxNigOkiRo2C66+HO++E\n339XwyuOtXWrWZ7VsoB9++Cxx6B7d3P7+4gRanjFsXbvNofZ9HTMWYdrrjFjdFasgEcfVcMbAAJu\nTG/fvnBWtaWkLOlFWGoq6WXLEtOpE5Hz58PXX5vmoW1bu8sUydOcGfEkjBqZld2kMW2JHD/eHIXn\nzVOzK46UtSqVK7cvdulCyLR0c3b3wQdN06BmVxwod3b7delCWN958NlnMHw4PPSQmt0AEnBNb9vW\nM/n5xS68tSnbCim//QZNmxL5559ay10c660+S5g6ajeLD52e7Dzu55+hfXsiP/5YY3fFkZLi4/nv\ns8Mp/89tDOMlAOLmz4eKFYmcOBFuusnmCkXylhQfz09dutJ2U2WuZhlg+oWq11xD5IoVULeuzRWK\npwXM8IbUVPPrrx+MyNHwAgxMSyOhfHk1vOJoB+e/ybRDOVeeGJiZScK2bWp4xbFmjxzJ8H9+53Zm\nZm0bePQoCVdcoYZXHG32yJH02LSH1xhECmau6IFpaSRUqaKGN0AFxJnebdvgnjvSWBw3jbBVq/J8\njlZIESea9oPF+SFbuHLHbMqvXEhdzlwRUNkVJ9q4JpUTi1cRtn494aTSkl9yPB56aoJeEafZu5fM\nOb8Ttno1VTjMLFrneFjH3MDlyKY39xibmC5dzpz37t9/Yc4cSEri7DlzSNh+jFJfXEF6PuPGtEKK\n+EKh2c3MhFWrsrKbNrsSqRX2QnQ10mvXznO9bGVXvM2tY+7x42aqxzlzYM4cli08h5NnNSA9nzvu\nlVvxBbeyu2MHJCWZrzlzWLC5Lu9W6sellcrluU9lN3A5rulNio9nVteuDEzONiY3ORl27SIyNDQr\ntBw4gHVTJCEtIqFDB6pddRWEhhITH09crtf3ioigtVZIES/LN7sbNhCZkWGy+/vvbK/akLNuuZSQ\nu+/m/ncj4dxzAZRdsUW+uT12jMjKlbM+oLFiBYcbNqPSzU0I6dmTB264ASpXJim+qXIrtsg3u7t3\nExkWdrpf2LvXDLWJjIT27bnuykaM3BfGxmXDld0gU9JbErcAh4EMIA24LtfjRZ53r3erVgyYPfuM\n7X3KlKH/PfeY0LZoAZdfTs9epbjuOrjvvpzPTYqPJyHbCinRWiHFb3lxvsgt5J/dYs0XmW92K1Sg\n/4BG3lUAAAu4SURBVBNPQGQk1k2R3PpYPT7+GC666Mx9KLuBw0vZ3YKvjrmlStH/ppvM8TYyEpo1\n4+6Hy/Pqq9C8ec7nKreBxe+zW7o0/e+6y2S3RQu44goSfilFSAi0bJnzucpu4HAntyU902sBUcD+\nEu4nS9ipO9JyCW3aFCZMyLEtNhYqVTrzuZFt2ii0UhjfZbdJE7PqH+Zf488/Q6l8biFVdqUQvstt\n8+aQmJhj24QJkNeVX+VW3OC77DZrBpMm5dgWHp73cVfZDS6emL3Bo58G00+cyHN7Rrkzx97Urw+V\nK3vy3SXIeC67lkX6v//m+dDR0JrccYdZSA3yb3hF3OTZY24+M4NklC/P8eNmmt1jx8w2DXWUEvJs\ndvNpek/1C8nJroUmMKMbcl+hkOBT0v9+LeBnYAnQocTVzJpFzNq1xOWaKqRXRATRrjE2+/bBE0+c\nDrJIMXkuu5mZ8PLLxJw8Sdx55+V4qFdEBHd1b0/fvlCmTIneRQQ8fczds4eYv/8mLtfZg1PH3PLl\noXNnKF++xO8k4tnsJiYSs2ZNgf1Cnz7wxx8lficJICX91FUP2AHUAhKAWOD3bI+7P0ZnwgTo0gW+\n/56kgwfzHWOTkQG//AIxMSWsXPyCF8f0FpRd93Obng7PPAPr18OMGSQtWEDCqFGkHA5ne2pDOr95\ngy6dBSkvZddzx9ytW82BtG1bkpo2JeH99wlNSeFIqVrUiXyFV/vlHnIpwcLx2Z06FTp2hAkTSDp+\nPEe/0PKFWFrc0ca1Qy2mFkx8MaZ3h+vXPcAUzMD07CGmX79+Wd9HRUURFRV15l5Gj4aBAyEhAf7z\nHyIh30YhNFQNbyBLTEwkMdc4Qi8pMLtu5fbECbNE5cmTJrsVKmSND1u7FiZPhkj1u0HDR9n1zDF3\n7Vpo1QpefBFefNEcc++4A4CNG2HcOC9ULo7lV9n9/HPo1Qt+/BGuuSZHv3DiBNx6K/wUaYY+quEN\nbMXJbUkiUR4IBY4AFYDZwBuuX08p+JObZcGAAfDFF6ZpuOCCfJ86eTLUqQM33liCisXveOmMQ2HZ\nLfyMw6FDcPfdUK8efPkllClDcjJUqwbVq3u4WvFLXshuyY+5AP/7H9x5JwweDO3aAbBmDdSoYY6x\nIo7N7tCh8P77MGsWXHJJnk9Zvx4uvtgTJYu/cSe3JRnTWwfzKW05sAiYQc4AFywz05xlmDQJ5s4t\nsOEFqFIFKlYsQbUip5Usu7t3w803Q8OGMH581mDdTz+FpUu9UK2IUbLcgpk6pE0bGDMmq+EFcxhe\nudKDlYrkVLLsWhb07Aljx8Lvv+doeBcvhrffPv1UNbxSEG+f/M/7k1taGjz1FGzZAtOnQ9WqXi5D\n/JUXx/QWJP8zDlu2mPE1Dz/MyV79WLQ4hJtu8mlt4iccl93Jk+HZZ2HSJLZHRDJ7NrRv79vixD84\nKrsZGeZuyhUrYOZMMqrVZMUKuPpq8/DOnWZIjq4Ci7fP9BbP8eNw771w4IC5RFFAw5uZacaWZWb6\nsD6R/Kxebea9eeEFeOMN9u0P4aOPzEkIEUcbM8ZMbD5rFkRGUrYs7PfYbKkiXpKSYubM27LF3MFe\nsyYnTsArr5yewaluXTW84j7fNr0HD5qbJ6pVgylTCp0H58gRWL7cfNATsdXChXDLLSx8Zgzb7usC\nmOG8X3+tmyXEwSzLXPsdNIjYW9eQXLkxADVrwksv2VybSEGOHDFDcUqV4sm6P7Fyi1mJqmJF0/+G\nlfQ2fAlKvmt6d+6EqChzTeLLL6F06UJfUqUKvPuuW08V8Z7Zs82NP2PHsrBKazZtsrsgETdYFpkv\nv2rGnc+dywMdqlK7tt1FiRRu7fwDrG72DFx4IXz7LS/3CM1z2XaRovLNmN7NmyE62qwq0bt3oafG\nkpOhQgVz2UKCm91jyw5/8T2zus7kgfgndQ1NisTW7KanM63VB8SvOo+P10ZqShEpEluzu3Ur/206\nkjJNr6bt5Id1KU3c5ot5egvVu3lzYtatI3LAAHMThRsSEqBcOdMji9ihd6tWxEREcNWUX1lw1++0\nbV7L5/8DiBRH7+hoYo4c4daKdWix7L9QvYLdJYm4pfeNNxKzfj2P9OgBLz1idzkSgLx/pheIq1OH\nVp99ppWppMjsOuOwlksYF5ZMq9GjiXzmGR+/vQQCu7K7l+oMq5BKq/HjibznHh+/vQQC2870AnG1\na9Nq7Fj1C1Jkjpm9YeCuXSSMGlXo8/75xwfFiLjhKBUZmJ5OwqRJdpciUiQ12M/AY8dI+Ogju0sR\nKbKBu3e71S+IFIfPbmQLTUkp8PEjR8xMZseP+6ggkQI0wawyUVhuRZxK2RV/peyKt/hs0o+M8PAC\nH69UyaysUsr3MweL5Kuw3Io4lbIr/krZFW/xSYvZKyKC6NjYPB9LSfn/9u4nNK4qCsD4F5MWrRFE\nXIi2MC61CApiRasVLbYV0bW7unAlVRBC6aJrl11U6EZcSKQuFEURXAgGRbAiNlStihED/q2CEhCL\nJDQu7kTT2CZjJ+e8d/O+HzwmKaEnw3wZ3uTd3Pl3H15PeNUmq3UrtZntqla2q0jhv+k9vGcPew8c\nuOii9KNHyzuuHTwY/Z1Ig1urW6mtbFe1sl1Fy9mndxULCzA/X7Yok1Zqep9e6VLZrmplu6pRK/bp\nvZjFxbLn9NiYbycoSZKkWI2soj17FnbsKDs2SJIkSdEaW94wOwu9XvB0Vc/LbKqV7apWtqsaDdJt\n42t6pdX45Kta2a5qZbuqUWvekW3J5CQcO5Y5UZIkSUr+Q7bdu13HK0mSpHwub1CreZlNtbJd1cp2\nVaOM5Q17gS+Br4ELvr3EwgJMTMDc3JCTpPWzZrdSS9muamW7atwwJ72jwHOUkG8GHgNuWvlFIyOw\nfTuMjw8x6X+YmprKGeTMWg3UbRO68ph2ZWaAVrbblcezqYZsN05XOurKzEEMc9J7BzADzALzwMvA\noyu/aHQU9u8vtxm68uB2ZWaAgbptQlce067MDNDKdrvyeHrSOxTbdWYrDHPSewPw3bLPv+//m9Rm\ndqta2a5qZbtqhWFOel1xrhrZrWplu6qV7ap6dwJvL/v8EP9dnD5Did3D41KPGdaX3XpkHbbrUeth\nux41Huvd7XnGgG+AHrAZmKYFC9OlNditamW7qpXtakPYB3xFObs+1PD3Ig3KblUr21WtbFeSJEmS\natXEJtQvAGeAT5PmAWwD3gU+Bz4DnkqYeTlwgnJ56DTwbMLMJaPASeDNpHmzwKn+zI+SZma3a7fx\nsrsF241iu7Fm2fjdgu1m6EK7QLmjM5S1O5vIW7tzD3AbuRFfB9za/3iccukm475u6d+OAR8COxNm\nAjwDvAS8kTTvW+CapFnQTLt2Gy+7W7DdKLYbqwvdgu1maGW7w74N8YU0tQn1+8DvCXOW+5nyQwrw\nB/AFcH3C3D/7t5spTxq/JczcCjwEPE/ue7JnzmqiXbuN1VS3JM+z3Vhdanejdwu2G6217Uac9HZ1\nE+oe5ZXjiYRZl1F+eM5QLpecTph5BJgAziXMWrIIvAN8DDyRMK+L7faw2wi2G6+H7a43u83Rw3bX\n20DtRpz0Lgb8n203DrwCPE15BRftHOUyyVbgXuC+4HkPA79Q1spkvmq7m/LEsA94knJJKlLX2rXb\nOLYby3Zj2G08240xULsRJ70/UBZsL9lGefW2UW0CXgUmgdeTZ88BbwG3B8+5C3iEsmbmOHA/8GLw\nTICf+re/Aq9RLoVF6lK7dhvLduPYbhy7jWW7cbLb/UeTm1D3yF2YPkJ5MI8kzrwWuLr/8RXAe8AD\nifN3kfPXmFuAq/ofXwl8ADwYPLOpdnvYbbSsbsF2I9lunC51C7abYSO3e54mNqE+DvwI/EVZI/R4\nwsydlEsH05Rf5Z+kbL8S6Rbgk/7MU5R1M5l2kfPXmDdS7uM0ZXuXrI6y27XbHFndgu1Gst04XekW\nbDfLRm9XkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJbfI3yVKl4tZGLm8AAAAA\nSUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x11877ad50>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def loocv(x, y, fit, pred, deg):\n",
      "    \"\"\"LOOCV RSS for fitting a polynomial model.\"\"\"\n",
      "    n = len(x)\n",
      "    idx = np.arange(n)\n",
      "    rss = np.sum([(y - pred(fit(x[idx!=i], y[idx!=i], deg), x))**2.0 for i in range(n)])\n",
      "    return rss"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# RSS does not detect overfitting and selects the most complex model\n",
      "for deg in range(1, 5):\n",
      "    print 'Degree = %d, RSS=%.2f' % (deg, np.sum((y - np.polyval(np.polyfit(x, y, deg), x))**2.0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Degree = 1, RSS=38.85\n",
        "Degree = 2, RSS=4.19\n",
        "Degree = 3, RSS=4.12\n",
        "Degree = 4, RSS=3.92\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# LOOCV selects the correct model\n",
      "for deg in range(1, 5):\n",
      "    print 'Degree = %d, RSS=%.2f' % (deg, loocv(x, y, np.polyfit, np.polyval, deg))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Degree = 1, RSS=408.47\n",
        "Degree = 2, RSS=43.66\n",
        "Degree = 3, RSS=48.67\n",
        "Degree = 4, RSS=94.94\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Calculation of Cook's distance"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "D_i = \\frac{\\sum_{j=1}^n (\\hat Y_j - \\hat Y_{j(i)})^2}{p \\text{MSE}}\n",
      "$$\n",
      "\n",
      "The calculation of Cook's distance involves the fitting of $n$ regression models, so we want to do this as efficiently as possible."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def cook_dist(X, y, model):\n",
      "    \"\"\"Vectorized version of Cook's distance.\"\"\"\n",
      "    n = len(X)\n",
      "    fitted = model(y, X).fit()\n",
      "    yhat = fitted.predict(X)\n",
      "    p = len(fitted.params)\n",
      "    mse = np.sum((yhat - y)**2.0)/n\n",
      "    denom = p*mse\n",
      "    idx = np.arange(n)\n",
      "    return np.array([np.sum((yhat - model(y[idx!=i], X[idx!=i]).fit().predict(X))**2.0) for i in range(n)])/denom"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# create data set with outliers\n",
      "nobs = 100\n",
      "X = np.random.random((nobs, 2))\n",
      "X = sm.add_constant(X)\n",
      "beta = [1, .1, .5]\n",
      "e = np.random.random(nobs)\n",
      "y = np.dot(X, beta) + e\n",
      "y[[7, 29, 78]] *= 3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# use Cook's distance to identify outliers\n",
      "model = sm.OLS\n",
      "d = cook_dist(X, y, model)\n",
      "plt.stem(d);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGFxJREFUeJzt3X+QHOVhp/FnTytWGAyWOTsKQik5C3fG5TgnOcjCxGhs\ni13FUJHrnDJQwaaMz6XKFVold3WR9YNoUoSq43JOzAqCdTnF8blIlDofdmGtza4oZyQqhwE5WNiJ\nxEmLVZEEtomwHYuyZC3M/fH27vbMzkz37szsLPs+n6qpnenut/udt7u/0/12zw5IkiRJkiRJkiRJ\nkiRJkiRpjlsHHAGOApsbTHcNMAZ8ZAZlJUlzxALgGLAcWAh8G7i6znTfAPYyGfx5y0qSZtG/yhi/\nihDex4HzwB5gfY3pNgJfAl6aQVlJ0izKCv6lwInU65PJsOpp1gMPJq/L0ygrSZplWcFfzhgP8Fng\n08m0Xckjb1lJ0izrzhh/CliWer2McOSe9m5CNw7AvwZ+g9C1k6csvb295dHR0WlUWZIEjAJXtmPG\n3cnMlwMXkH2B9vPAv59m2bKCHTt2dLoKc4ZtMcm2mGRbTKKJXpWsI/4x4E5gmHCXzm7gMLAhGb9r\nBmUlSR2UFfwAX08eafUC/xM5ykqSOijr4q5mUaFQ6HQV5gzbYpJtMcm2aI2u7EnaLumukiTl1dXV\nBTPMcI/4JSkyBr8kRcbgl6TIGPySFBmDX5Iik+c+fklqu6GhAwwOjnDuXDc9PWMMDPRx443Xd7pa\n85LBL6njhoYOsGnTMKOj90wMGx3dBmD4t4FdPZI6bnBwpCL0AUZH72Hnzn0dqtH8ZvBL6rhz52p3\nPpw9u2CWaxIHg19Sx/X0jNUcvmjRq7NckzgY/JI6bmCgj97ebRXDenu3snHjDR2q0fzm/+qRNCcM\nDR1g5859DA/fTX//XWzceIMXdhto5n/1GPyS5pSuLjASsvlP2iRJuRn8khQZg1+SImPwS1Jk8gT/\nOuAIcBTYXGP8euAQ8AzwLeADqXHHgWeTcU81U1FJUmtkXRFeADwHrAVOAU8DtwKHU9NcBLySPP8V\n4MvAlcnr7wHvBl5usAzv6pE0wbt68mnnXT2rgGOEI/fzwB7CEX7aK6nnFwP/XF2/mVRMktQeWcG/\nFDiRen0yGVbtw4SzgK8DA6nhZeAx4CDwqZlXU5LUKln/ljnvCddXksf7gC8C/zYZfh3wIvAWYB/h\nWsHj06+mJKlVsoL/FLAs9XoZ4ai/nseTeV4GnCaEPsBLhL7/VdQI/mKxOPG8UChQKBQyqiVJcSmV\nSpRKpZbMK6v/vZtwcfeDwAuEO3OqL+72As8Tzg5WAv87GfYGwsXhnxIuAI8Af5j8TfPirqQJXtzN\np5mLu1lH/GPAncAwIcR3E0J/QzJ+F/AR4OOEi79ngFuScUuAh1PLeYipoS9JmmVz4Y4bj/glTfCI\nPx//SZskKTeDX5IiY/BLUmQMfkmKjMEvSZEx+CUpMga/JEXG4JekyBj8khQZg1+SImPwS1JkDH5J\niozBL0mRMfglKTIGvyRFxuCXpMgY/JIUGYNfkiJj8EtSZPIE/zrgCHAU2Fxj/HrgEPAM8C3gA9Mo\nK0maZVk/1LsAeA5YC5wCngZuBQ6nprkIeCV5/ivAl4Erc5YFf2xdUoo/tp5PO39sfRVwDDgOnAf2\nEI7w015JPb8Y+OdplJUkzbKs4F8KnEi9PpkMq/ZhwpH814GBaZaVJM2i7ozxeU+4vpI83gd8EXj7\ndCpRLBYnnhcKBQqFwnSKS9K8VyqVKJVKLZlXVv/QaqBIuEgLsAV4Dbi3QZlRQjfPVTnL2scvaYJ9\n/Pm0s4//ICHAlwMXADcDj1RN05ta+Mrk7+mcZSVJsyyrq2cMuBMYJtyls5vQl78hGb8L+AjwccIF\n3DPALRllJUkdNKPThBazq0fSBLt68mlnV48kaZ4x+CUpMga/JEXG4JekyBj8khQZg1+SImPwS1Jk\nDH5JiozBL0mRMfglKTIGvyRFxuCXpMgY/JIUGYNfkiJj8EtSZAx+SYqMwS9JkTH4JSkyBr8kRSZP\n8K8DjgBHgc01xv82cAh4Fvg74F2pcceT4c8ATzVTUUlSa3RnjF8A3A+sBU4BTwOPAIdT0zwPXA/8\nhPAh8T+A1cm4MlAAXm5ZjSVJTck64l8FHCMcuZ8H9gDrq6Z5ghD6AE8CV1SNn9GvwEuS2iMr+JcC\nJ1KvTybD6vkk8LXU6zLwGHAQ+NRMKihJaq2srp7yNOb1fuAO4LrUsOuAF4G3APsI1woery5YLBYn\nnhcKBQqFwjQWK0nzX6lUolQqtWReWd0wq4Eioe8eYAvwGnBv1XTvAh5OpjtWZ147gDPAZ6qGl8vl\n6Xy+SJrPurrASMjW1dUFM+xKz+rqOQhcBSwHLgBuJlzcTfslQujfRmXovwF4Y/L8IqAP+M5MKilJ\nap2srp4x4E5gmHCHz27CHT0bkvG7gD8AFgMPJsPOEy4KLyF8IIwv5yFgpFUVlyTNzFy448auHkkT\n7OrJp51dPZKkecbgl6TIGPySFBmDX5IiY/BLUmQMfkmKjMEvSZEx+CUpMga/JEXG4JekyBj8khQZ\ng1+SImPwS1JkDH5JiozBL0mRMfglKTIGvyRFxuCXpMgY/JIUmTzBvw44AhwFNtcY/9vAIeBZ4O+A\nd02jrCRplmX9UO8C4DlgLXAKeBq4FTicmuZa4B+BnxCCvgiszlkW/LF1SSn+2Ho+7fyx9VXAMeA4\ncB7YA6yvmuYJQugDPAlcMY2ykqRZlhX8S4ETqdcnk2H1fBL42gzLSpJmQXfG+OmccL0fuAO4brpl\ni8XixPNCoUChUJjGYiVp/iuVSpRKpZbMK6t/aDWhz35d8noL8Bpwb9V07wIeTqY7Ns2y9vFLmmAf\nfz7t7OM/CFwFLAcuAG4GHqma5pcIoX8bk6Gft6wkaZZldfWMAXcCw4S7dHYT7srZkIzfBfwBsBh4\nMBl2nnBht15ZSVIHzeg0ocXs6pE0wa6efNrZ1SNJmmcMfkmKjMEvSZEx+CUpMga/JEXG4JekyBj8\nkhQZg1+SImPwS1JkDH5JiozBL0mRMfglKTIGvyRFxuCXpMgY/JIUGYNfkiJj8EtSZAx+SYqMwS9J\nkckT/OuAI8BRYHON8W8HngDOAv+5atxx4FngGeCpGddSktQy3RnjFwD3A2uBU8DTwCPA4dQ0p4GN\nwIdrlC8DBeDlZisqSWqNrCP+VcAxwpH7eWAPsL5qmpeAg8n4Wmb0K/CSpPbIOuJfCpxIvT4JvGca\n8y8DjwGvAruAP59W7dRWQ0MHGBwc4dy5bnp6xhgY6OPGG6/vdLUktVlW8JebnP91wIvAW4B9hGsF\nj1dPVCwWJ54XCgUKhUKTi1WWoaEDbNo0zOjoPRPDRke3ARj+0hxUKpUolUotmVdWN8xqoEi4wAuw\nBXgNuLfGtDuAM8Bn6syr3vhyudzs54umq79/OyMjf1Rj+F08+ujdHaiRFHR1gZGQraurC2bYlZ7V\nx38QuApYDlwA3Ey4uFuzHlWv3wC8MXl+EdAHfGcmlVTrnTtX+2Tv7NkFs1wTSbMtq6tnDLgTGCbc\n4bObcEfPhmT8LmAJ4W6fSwhnA5uAdwBvBR5OLechYKSFdVcTenrGag5ftOjVWa6JpNk2F+64saun\nA2r18ff2buW++9bZx6+Osqsnn2a6egz+iA0NHWDnzn0MD99Nf/9dbNx4g6GvjjP48zH41RR3NM0l\nbo/5tPPiriRpnjH4JSkyBr8kRcbgl6TIGPySFBmDX5IiY/BLUmQMfkmKjMEvSZEx+CUpMga/JEXG\n4JekyBj8khQZg1+SImPwS1JkDH5JikzWb+6+bg0NHWBwcIRz57rp6RljYKDPX5eSJPIF/zrgs4Qf\nW/+fwL1V498OfB5YAWwDPjONsm1R6/dkR0e3ARj+kqKX1dWzALifEODvAG4Frq6a5jSwEfjvMyjb\nFoODIxWhDzA6eg87d+6bjcVL0pyWFfyrgGPAceA8sAdYXzXNS8DBZPx0y7bFuXO1T2TOnl0wG4uX\npDktK/iXAidSr08mw/JopmxTenrGag5ftOjV2Vi8JM1pWX38zfzWfe6yxWJx4nmhUKBQKDSxWBgY\n6GN0dFtFd09v71Y2blzX1HwlqVNKpRKlUqkl8+rKGL8aKBL66QG2AK9R+yLtDuAMkxd385Ytl8vN\nfL7UNjR0gJ079zE8fDf9/XexceMNXtito6sL2rAKpBlxe8ynq6sLsjO8dtmM8d3Ac8AHgReApwgX\naQ/XmLYI/JTJ4M9bti3BP86NKJttpLnE7TGfZoI/q6tnDLgTGCbcpbObENwbkvG7gCXA08AlhCP6\nTYS7eM7UKStJ6qAZfVq0mEf8HWYbaS5xe8ynmSN+/2WDJEXG4JekyBj8khQZg1+SImPwS1JkDH5J\niozBL0mRMfglKTIGvyRFxuCXpMgY/JIUGYNfkiJj8EtSZAx+SYqMwS9JkTH4JSkyBr8kRcbgl6TI\n5An+dcAR4Ciwuc40g8n4Q8CK1PDjwLPAM4QfW5ckdVjWj60vAO4H1gKnCD+q/giVP5r+IeBK4Crg\nPcCDwOpkXBkoAC+3rMaSpKZkHfGvAo4RjtzPA3uA9VXT/CbwheT5k8CbgF9IjZ8LP+guSUpkBf9S\n4ETq9clkWN5pysBjwEHgUzOvpiSpVbK6eso551PvqP7XgReAtwD7CNcKHs85T0lSG2QF/ylgWer1\nMsIRfaNprkiGQQh9gJeALxO6jqYEf7FYnHheKBQoFAoZ1ZKkuJRKJUqlUkvmldX/3g08B3yQEOJP\nAbcy9eLuncnf1cBnk79vIFwc/ilwETAC/GHyN61cLuc9sZi+ri5o4+znBdtIc4nbYz5dXV0ww2uo\nWUf8Y4RQHyaE+G5C6G9Ixu8CvkYI/WPAK8AnknFLgIdTy3mIqaEvSZplc+GOG4/4O8w20lzi9phP\nM0f8fnNXkiJj8EtSZAx+SYqMwS9JkTH4JSkyBr8kRcbgl6TIGPySFBmDX5IiY/BLUmSy/lePEkND\nBxgcHOHcuW56esYYGOjjxhuv73S1pOi4LzbP4M9haOgAmzYNMzp6z8Sw0dFtAG5w84Rh8vrgvtga\n/pO2HPr7tzMy8kc1ht/Fo4/e3dzM54B2/VOsuR6m4/U7deolnn++i5/97HMT43p7t3Hfff1zqr6x\naLQ9zvd9cTra+W+ZBZw7V7uZzp5dMMs1ef2Y60dmlfXbDlSGyejoPezcedecqKsmuS+2hhd3c+jp\nGas5fNGiV2e5Jq8fg4MjFaEP42G6ryXzHxo6QH//dgqFIv392xkaOtBE/QyT1wv3xdbwiD+HgYE+\nRke3VQRZb+9WNm5c18FazW3tPDJrxdlEZf2KNacxTOYe98XWMPjrqO6fvu22pXzzm3cxPHw3/f13\nsXHjuqa6AV4v/d8zrV87j8zqn03k75qZWr9twNwKk7m+jXTC+PvfubN1+6I6o9wKe/fuL/f1bSuv\nWbOj3Ne3rbx37/5yuGo8s3n19m4th0tM4dHbu7W8d+/+Gc1vOvPvhOr3NLV++8sXXvjR8jvfuami\nbRup/R63tOQ9rlmzo2K+4481a3bknke99wjlcn//9o6ti/r16+w2Uku9fa5ZefexrOnaVb+5Anhd\n/05Z0w3Q6qDu69tWM1j6+7e3JPgbzb8Tqt9TZf32l2FmAbR37/6JNkuHabM7ZKvar1b9WnQc0rSs\n95i3DdsVfrX2uSVL7iivWPE75TVrdpRXrPjkxPPpLnd8HWTVvdG6ej18cJbLza0f2hz864AjwFFg\nc51pBpPxh4AV0yzbdOO1OqgbHVG2IhiyjlhbvbNOdweqV7+skM1z1pW1Q6bnUS88Wn02ka5fu4J/\nukF96aW3191G8obaTMOvUV3Hxy1efHNV3dIHCI0PFvJsj3nq3mhd1c6E/eXLLvvotPardp41NPvh\nRBPBn9XHvwC4H1gLnAKeBh4BDqem+RBwJXAV8B7gQWB1zrJNGe8DffLJkzXH57mQWKsfNat/Om/f\na73pps7/ADDCt771NCtX/gdefPESvv/9P5kY28xtkDO5EFrv/afr+9RTRykUivzLv5wELuDnP391\nyr3w48uByeXU7p/v5/bbH2DJkodS87gPOE+637263ul+3tWrr2BwcIQ//uNvTNTpkkve2pK+8fR6\nrJ73tddezhNPvDBlHVeXqVynB3j88Qfo7X2YhQvPTMyvcrrtqRqUCDfgjfDssye4/fYHOH36b6ra\ncOo1jqltfYDR0S4+9rHdXHPNSM26A3W3l8pxxapWGmFyXaWfVy53+fL/lWv7rlf3W265h/e+dySp\na/11OvXmggPAMKdP/w3799dfLjT+fke6TLM5kLV+2nlNJ+vm/2uBHYQjd4BPJ3//a2qazwF/C4xv\niUeAAvC2HGUByn192yo2wvTOVW9Hq1wpU+/DhvCljuHhu9m7t/aOW2uHvPDCB3jrW7v54Q/fWPWF\nnq3cd986brrpenp703cVhDK9vb/I5ZdfXKd+ldMtXHimarnjFxaLwFjVewkfCosXn2D58oU126JR\nm333u4erQmJyftdcs4yBgT5uuul6yuXGG3xl+eGkvunn+dbBk0+e5Cc/+cs680vP4/pkXOWyL7vs\nAd75zqsndqCbbrqevXvTH27p+dVfP+mdcPz9w+SXh2q3xdR5d3f/FWNjlet46vaTfl/12q/RdEUq\nPwSLVAZvWKeXXnqSX/7lbsbX/aFD3+PHP/5Carr0BezadV+w4FXOnPlSzXZ/7bUufvSjPTXqWqtO\naenlVpervT2uWVNk//5ijXmE5Yzvg+l1ld4Pjh9/OVXXxstN71dTP3yrt+nQFkuWvKnu/p03B06e\n/FHD9ZO93a6BNn0J97eAP0+9vg3YWTXNV4H3pl4/Brwb+EiOsgBl2F/u7t6Q45QxPV1WP/SWiT7b\nydOp6ukazaP2xb76p7fTq9+SJXeUV678j1Xz25E80sPytEWjNqs3v8lTy9qn1pPv/8ILN6SGb6vz\nvHH32OS8q0/B682jen61uw+g+rS+8TqdbLPKeUyePle3Rb15N1pW9XS12yZ7uv1l2F5jeHW5PNt3\no2WnyzVq90bbUqNlNXqPedZp/fKV62om+0GjNmvUFnm3s0bT1XuPebfb9vXx5wnvrwLXpV7PIPjr\n7Vx5QyK9kzQK6rw7ZOWjsl+t0SM7CNOP8WsQlWXqbQyN2iJvm9V/1N/RKi+CLl788Vzzq/8es3bI\nesPr1y/PNYmseaTXb2VbNJp383XPN131/KZ+CNYuX2+6RnXPu52Nzz+s3xUrPllesuT3UuPqLTdf\nW9Tq/65V7/of+pX1u+yyj+asT6P1nXe7yDtdo/WTp83aF/yrgUdTr7cw9SLt54BbUq+PAL+QsyzQ\nO/4GfPjw4cNH7sein9Em3cAosBy4APg2cHXVNB8CvpY8Xw18cxplJUlz0G8AzwHHCEftABuSx7j7\nk/GHgJUZZSVJkiTFIs8XvOarZYTbYP8B+C4wkAx/M7AP+H+EG6Lf1JHazb4FwDOEmwUg3nZ4E/Al\nwvdd/pHw3ZhY22ILYf/4DvBXQA/xtMVfAD8gvPdxjd77FkKOHgH6ZqmOM7KA0AW0HFhIfNcAlgD/\nLnl+MaFL7GrgvwG/nwzfzNTvPcxX/wl4iPAlP4i3Hb4A3JE87wYuJc62WA48Twh7CN8Tup142uJ9\nhP+CkA7+eu/9HYT8XEhot2PM4X+5fy2Vd/18mskvecXoK4RvOY/fFQXhw+FIx2o0e64g3Ab8fiaP\n+GNsh0sJYVctxrZ4M+FgaDHhA/CrwA3E1RbLqQz+eu+9+o7JRwk32tTVyU+FpcCJ1OuTybAYLSd8\nuj9JWLE/SIb/gMkVPZ/9KfBfgNdSw2Jsh7cBLwGfB/6e8D2Yi4izLV4GPgP8E/AC8GNCN0eMbTGu\n3nu/nJCf4zKztJPBX+7gsueSi4H/A2wCflo1bvye3fnsJuCHhP79el8/j6EdIBzZrgT+LPn7ClPP\ngmNpi17gdwkHRZcT9pPbqqaJpS1qyXrvDdulk8F/inCBc9wyKj+1YrCQEPpfJHT1QPgkX5I8/0VC\nKM5n7wV+E/ge8NfABwjtEVs7QNj+TxL+oSGEi7wrge8TX1v8GvB/gdOEf2D1MKF7OMa2GFdvn6jO\n0iuSYXV1MvgPEv6j53LCF7xuZvLCXgy6gN2EOzc+mxr+COEiFsnfrzC/bSVstG8jfAP8G8DHiK8d\nIITaCeDfJK/XEu5q+SrxtcURQj/1hYR9ZS1hX4mxLcbV2yceIew7FxD2o6uAp2a9dtMQ8xe8fp3Q\np/1tQjfHM4TbW99MuNA5329Xq2UNkx/+sbbDrxKO+A8RjnIvJd62+H0mb+f8AuEMOZa2+GvCtY2f\nEw4GPkHj976VkKNHgP5ZrakkSZIkSZIkSZIkSZIkSZIkSZIkadL/B26cKGEjjy/gAAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x118723b10>"
       ]
      }
     ],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}