{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Linear Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from scipy import stats\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLM: Binomial response data\n",
    "\n",
    "### Load Star98 data\n",
    "\n",
    " In this example, we use the Star98 dataset which was taken with permission\n",
    " from Jeff Gill (2000) Generalized linear models: A unified approach. Codebook\n",
    " information can be obtained by typing: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "::\n",
      "\n",
      "    Number of Observations - 303 (counties in California).\n",
      "\n",
      "    Number of Variables - 13 and 8 interaction terms.\n",
      "\n",
      "    Definition of variables names::\n",
      "\n",
      "        NABOVE   - Total number of students above the national median for the\n",
      "                   math section.\n",
      "        NBELOW   - Total number of students below the national median for the\n",
      "                   math section.\n",
      "        LOWINC   - Percentage of low income students\n",
      "        PERASIAN - Percentage of Asian student\n",
      "        PERBLACK - Percentage of black students\n",
      "        PERHISP  - Percentage of Hispanic students\n",
      "        PERMINTE - Percentage of minority teachers\n",
      "        AVYRSEXP - Sum of teachers' years in educational service divided by the\n",
      "                number of teachers.\n",
      "        AVSALK   - Total salary budget including benefits divided by the number\n",
      "                   of full-time teachers (in thousands)\n",
      "        PERSPENK - Per-pupil spending (in thousands)\n",
      "        PTRATIO  - Pupil-teacher ratio.\n",
      "        PCTAF    - Percentage of students taking UC/CSU prep courses\n",
      "        PCTCHRT  - Percentage of charter schools\n",
      "        PCTYRRND - Percentage of year-round schools\n",
      "\n",
      "        The below variables are interaction terms of the variables defined\n",
      "        above.\n",
      "\n",
      "        PERMINTE_AVYRSEXP\n",
      "        PEMINTE_AVSAL\n",
      "        AVYRSEXP_AVSAL\n",
      "        PERSPEN_PTRATIO\n",
      "        PERSPEN_PCTAF\n",
      "        PTRATIO_PCTAF\n",
      "        PERMINTE_AVTRSEXP_AVSAL\n",
      "        PERSPEN_PTRATIO_PCTAF\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sm.datasets.star98.NOTE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data and add a constant to the exogenous (independent) variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.star98.load(as_pandas=False)\n",
    "data.exog = sm.add_constant(data.exog, prepend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The dependent variable is N by 2 (Success: NABOVE, Failure: NBELOW): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[452. 355.]\n",
      " [144.  40.]\n",
      " [337. 234.]\n",
      " [395. 178.]\n",
      " [  8.  57.]]\n"
     ]
    }
   ],
   "source": [
    "print(data.endog[:5,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The independent variables include all the other variables described above, as\n",
    " well as the interaction terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.43973000e+01 2.32993000e+01 1.42352800e+01 1.14111200e+01\n",
      "  1.59183700e+01 1.47064600e+01 5.91573200e+01 4.44520700e+00\n",
      "  2.17102500e+01 5.70327600e+01 0.00000000e+00 2.22222200e+01\n",
      "  2.34102872e+02 9.41688110e+02 8.69994800e+02 9.65065600e+01\n",
      "  2.53522420e+02 1.23819550e+03 1.38488985e+04 5.50403520e+03\n",
      "  1.00000000e+00]\n",
      " [1.73650700e+01 2.93283800e+01 8.23489700e+00 9.31488400e+00\n",
      "  1.36363600e+01 1.60832400e+01 5.95039700e+01 5.26759800e+00\n",
      "  2.04427800e+01 6.46226400e+01 0.00000000e+00 0.00000000e+00\n",
      "  2.19316851e+02 8.11417560e+02 9.57016600e+02 1.07684350e+02\n",
      "  3.40406090e+02 1.32106640e+03 1.30502233e+04 6.95884680e+03\n",
      "  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(data.exog[:2,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit and summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:           ['y1', 'y2']   No. Observations:                  303\n",
      "Model:                            GLM   Df Residuals:                      282\n",
      "Model Family:                Binomial   Df Model:                           20\n",
      "Link Function:                  logit   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -2998.6\n",
      "Date:                Fri, 24 Apr 2020   Deviance:                       4078.8\n",
      "Time:                        14:17:04   Pearson chi2:                 4.05e+03\n",
      "No. Iterations:                     5                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1            -0.0168      0.000    -38.749      0.000      -0.018      -0.016\n",
      "x2             0.0099      0.001     16.505      0.000       0.009       0.011\n",
      "x3            -0.0187      0.001    -25.182      0.000      -0.020      -0.017\n",
      "x4            -0.0142      0.000    -32.818      0.000      -0.015      -0.013\n",
      "x5             0.2545      0.030      8.498      0.000       0.196       0.313\n",
      "x6             0.2407      0.057      4.212      0.000       0.129       0.353\n",
      "x7             0.0804      0.014      5.775      0.000       0.053       0.108\n",
      "x8            -1.9522      0.317     -6.162      0.000      -2.573      -1.331\n",
      "x9            -0.3341      0.061     -5.453      0.000      -0.454      -0.214\n",
      "x10           -0.1690      0.033     -5.169      0.000      -0.233      -0.105\n",
      "x11            0.0049      0.001      3.921      0.000       0.002       0.007\n",
      "x12           -0.0036      0.000    -15.878      0.000      -0.004      -0.003\n",
      "x13           -0.0141      0.002     -7.391      0.000      -0.018      -0.010\n",
      "x14           -0.0040      0.000     -8.450      0.000      -0.005      -0.003\n",
      "x15           -0.0039      0.001     -4.059      0.000      -0.006      -0.002\n",
      "x16            0.0917      0.015      6.321      0.000       0.063       0.120\n",
      "x17            0.0490      0.007      6.574      0.000       0.034       0.064\n",
      "x18            0.0080      0.001      5.362      0.000       0.005       0.011\n",
      "x19            0.0002   2.99e-05      7.428      0.000       0.000       0.000\n",
      "x20           -0.0022      0.000     -6.445      0.000      -0.003      -0.002\n",
      "const          2.9589      1.547      1.913      0.056      -0.073       5.990\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())\n",
    "res = glm_binom.fit()\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantities of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of trials: 807.0\n",
      "Parameters:  [-1.68150366e-02  9.92547661e-03 -1.87242148e-02 -1.42385609e-02\n",
      "  2.54487173e-01  2.40693664e-01  8.04086739e-02 -1.95216050e+00\n",
      " -3.34086475e-01 -1.69022168e-01  4.91670212e-03 -3.57996435e-03\n",
      " -1.40765648e-02 -4.00499176e-03 -3.90639579e-03  9.17143006e-02\n",
      "  4.89898381e-02  8.04073890e-03  2.22009503e-04 -2.24924861e-03\n",
      "  2.95887793e+00]\n",
      "T-values:  [-38.74908321  16.50473627 -25.1821894  -32.81791308   8.49827113\n",
      "   4.21247925   5.7749976   -6.16191078  -5.45321673  -5.16865445\n",
      "   3.92119964 -15.87825999  -7.39093058  -8.44963886  -4.05916246\n",
      "   6.3210987    6.57434662   5.36229044   7.42806363  -6.44513698\n",
      "   1.91301155]\n"
     ]
    }
   ],
   "source": [
    "print('Total number of trials:',  data.endog[0].sum())\n",
    "print('Parameters: ', res.params)\n",
    "print('T-values: ', res.tvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First differences: We hold all explanatory variables constant at their means and manipulate the percentage of low income households to assess its impact on the response variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = data.exog.mean(axis=0)\n",
    "means25 = means.copy()\n",
    "means25[0] = stats.scoreatpercentile(data.exog[:,0], 25)\n",
    "means75 = means.copy()\n",
    "means75[0] = lowinc_75per = stats.scoreatpercentile(data.exog[:,0], 75)\n",
    "resp_25 = res.predict(means25)\n",
    "resp_75 = res.predict(means75)\n",
    "diff = resp_75 - resp_25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interquartile first difference for the percentage of low income households in a school district is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-11.8753%\n"
     ]
    }
   ],
   "source": [
    "print(\"%2.4f%%\" % (diff*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plots\n",
    "\n",
    " We extract information that will be used to draw some interesting plots: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "nobs = res.nobs\n",
    "y = data.endog[:,0]/data.endog.sum(1)\n",
    "yhat = res.mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot yhat vs y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.graphics.api import abline_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2de3hU5bXwfyvDQAICQUCUAIIoIFQRRdFSq6IVrVWpWqmX2mNr+/kd7alWqfj1glhbaKmXHq21HrU3bb2XiqhUC9WKpYhctEHCQZFLUAQhXCcwSdb3x8yEyWTvPXsuezKZrN/z5GFm73e/e82E7PW+6yqqimEYhmGkUtbWAhiGYRjFiSkIwzAMwxFTEIZhGIYjpiAMwzAMR0xBGIZhGI6YgjAMwzAcMQVhGICIDBYRFZFOPsb+h4i8nuH8p4pITfYStpjLt6yGkQumIIx2h4h8ICL7RaRPyvHl8Qfn4LaRrMXDe3fSzwpV/YeqDk8a94GInOUxz+ki0hS/fpeI1IjI1VnIc5uIPJrt5zE6NqYgjPbKWuCyxBsROQaoaDtxWlGpqgfFf0ZnOccmVT0I6AHcAvyPiIzMn4iG4Y0pCKO98gfgqqT3XwV+nzxARHqKyO9FZIuIrBOR74tIWfxcSER+LiJbReR94DyHax8WkQ9FpFZE7hCRULbCxncEG+Ov/wAMAubEdwjf9bpWY8wGtgOtFISI9BeR50Rkm4isEZFvxI+fA/w/YHJiJ5Ot/EbHxGyYRntlEfAVETkaWA1MBj4D3JE05l6gJ3AE0Bv4K/Ah8DDwDeALwBhgD/BMyvy/AzYDRwLdgOeBDcCvcxVcVb8iIqcC16jqK+nGx5XahUAl8I7DkD8B1UB/YATwsoi8r6ovichPgCNV9cpc5TY6HraDMNoziV3E54BVQG3iRHy1Pxm4VVV3qeoHwJ3AV+JDLgXuUdUNqroNmJF0bT/gXOAGVd2jqh8DdwNfzkC2rSJSF/+5OcvP119E6oCtwDTgK6rawtEtIgOJKcZbVLVeVZcDDyV9TsPIGttBGO2ZPwCvAUNIMS8BfYDOwLqkY+uAqvjr/sR2BMnnEhwOhIEPRSRxrCxlfDr6qGpD4o2InJ7BtQk2qeqANGP6A9tUdVfSsXXA2CzuZxgtMAVhtFtUdZ2IrAU+D3w95fRWIErsYb8yfmwQB3YZHwIDk8YPSnq9AdhHykM+z+SrjPIm4GAR6Z6kJJI/p5VrNrLGTExGe+frwARV3ZN8UFUbgSeBH4tIdxE5HPgOkAj5fBL4LxEZICK9gKlJ135IzF9xp4j0EJEyERkqIqflUe7NxHwjOaGqG4A3gBkiUi4ixxL7Th5Lus/ghHPeMDLB/tMY7RpVfU9Vl7ic/hYxB/T7wOvAH4FH4uf+B5gHrACWAs+mXHsVMRPVSmLRQ08Dh+VR9BnA93P0USS4DBhMbDfxZ2Caqr4cP/dU/N9PRGRpjvcxOhhiDYMMwzAMJ2wHYRiGYThiCsIwDMNwxBSEYRiG4YgpCMMwDMORdpcH0adPHx08eHBbi2EYhtGueOutt7aqat9Mrml3CmLw4MEsWeIW1WgYhmE4ISLr0o9qiZmYDMMwDEdMQRiGYRiOmIIwDMMwHDEFYRiGYThiCsIwDMNwpN1FMRmGYXQEZi+rZda8GjbVRehfWcGUicOZNKYq/YV5xBSEYRhGkTF7WS23PvsOkWgjALV1EW59NtZttpBKwkxMhmEYRcaseTXNyiFBJNrIrHk1LlcEgykIwzCMImNTXSSj40FhCsIwDKPI6F9ZkdHxoDAFYRiGUWRMmTicinCoxbGKcIgpE4cXVA5zUhuGYRQZCUe0RTEZhmEYrZg0pqrgCiEVMzEZhmEYjtgOwjAMo4iZvayW6XOq2b43CkBlRZjbLhhVkN2FKQjDMIwiZfayWqY8vYJoozYfq4tEmfLUCiD4pDlTEIZhlDRBl6xwmh/y42CeNa+mhXJIEG1SZs2r8TXnJ7v3cefLqzO+N5iCMAyjhAm6ZIXT/FOeXgEae4jnek+vxLh0SXPRxiZ+/8913PPKavbub/Qc64Y5qQ3DKFmCLlnhNH+0UZuVQ6739EqM8zr36uotnHPPa/zo+ZUcN7CSeTecmvG9wXYQhmGUMEGXrMhknmzuOWXi8FY+CIBwmTgmza3duocfz13JK+9+zODeXXnoqrGcefQhiEjG9wZTEIZhlDD9KyuodXgw56tkhdv8bmMzJWGSShfFtKs+yn0L1vDI62vpHCpj6rkjuHr8YLp0CjnO6xdTEIZhlCxTJg5v4SOA/JascJo/HJIWPohc7+mVMNfUpDy9dCM/e6mGrbv38aUTBjDlnOEc0r08q3ulYgrCMIySJeiSFW7zB3nPBEvXb2f6c9Ws2LiDMYMqefirYxk9sDKv9xDV1iFUxczYsWN1yZIlbS2GYRhGm/DRjnp++tIq/rysln49ujD13BFcOLqKsjJvP4OIvKWqYzO5l+0gDMMw2gH10UYefn0tv1ywhoYm5bozhvKfpx9Jty7BPcZNQRiGYRQxqsq86s38+IWVbNgWYeKofnzv8yMZ1Ltr4Pc2BWEYhlGk1Hy0i+lzqnnjvU8Y1u8gHrtmHOOP7FOw+5uCMAzDKDLq9u7nrpdX8+iidXQvD3P7haO4/KRBdAoVNrfZFIRhGIYH+a7l5DVfQ2MTf1y8nrteXs3OSJQrTz6cG88aRq9unfP1cTLCFIRhGIYL+a7l5DXfId27MH3OSmo27+KUI3oz7YKRjDi0R54+SXaYgjAMw3DBq5ZTttVZneab+uzb1EebGNCrggeuPJ6Jow7NujxGPjEFYRiG4UK+azm5XVcfbeLms4dxzalHUB7OrTxGPrFqroZhGC641U/KtpaT23WH9ijn+glHFZVyAFMQhmGUCLOX1TJ+5nyGTJ3L+Jnzmb2s1tc5L6ZMHE5FykM7l7pKk08cSGrCc0U4xNRzR2Q1X9CYickwjHaPl/MXyNrRnK9aTh/vqmfWSzU89dZGupd3IiRCXSRKVUB1mvJFoLWYROQc4BdACHhIVWemnO8JPAoMIqasfq6qv/Ga02oxGYaRyviZ8x3LblfFTTpu5xZOnRCoXPsbmvjNwrXcO38N+xoa+dr4IVw/4Ui6l4cDva8TRVWLSURCwC+BzwEbgTdF5DlVXZk07DpgpaqeLyJ9gRoReUxV9wcll2EYpUc2zuR8NQ1yQlWZv+pj7pj7Lmu37uHMEYfwvfOO5oi+BwV2zyAI0sR0ErBGVd8HEJHHgQuBZAWhQHeJxXMdBGwDGgKUyTCMIiWXhLR0jYGCbBqUypqPd3P78yt5bfUWhvbtxm+vPpHThx8SyL2CJkgFUQVsSHq/ERiXMuY+4DlgE9AdmKyqTakTicg3gW8CDBo0KBBhDcNoO3JNSEvXGCjIpkEJdkSi/Pff/pffvfEBFZ1D/OALI7nqlMMJF7g8Rj4JUkE4ZXmkOjwmAsuBCcBQ4GUR+Yeq7mxxkeqDwIMQ80EEIKthGG1IrglpfpzJs+bVUFsXISTSPHfytdnS2KQ8uWQDP59Xw7a9+/nyiQO56ezh9DmoS07zFgNBKoiNwMCk9wOI7RSSuRqYqTFP+RoRWQuMABYHKJdhGEVGrglp6cxTidf5LJsBsHjtNm57rpqVH+7kpMEH87vzR/Kpqp5ZzVWMBKkg3gSOEpEhQC3wZeDylDHrgTOBf4hIP2A48H6AMhmGUYSk8yF44dc8lc+yGbV1EWa88C7Pv/0h/XuWc+9lY/jCsYcVRXmMfBKYglDVBhG5HphHLMz1EVWtFpFr4+cfAH4E/FZE3iFmkrpFVbcGJZNhGPknH9VO0/kQvO7rpFicHvz5KJsR2d/Ir197jwdefQ9V+PaZR3HtaUOp6FxcGdD5ItBEOVV9AXgh5dgDSa83AWcHKYNhGMGRr2qnmSakpd7XidQHfy67FFVl7jsfMuOFVdTWRTjv2MO49dwRDOgVfFe3tsQyqQ3DyJpMzTZOu43EPJvqIvSsCFPZNcymuoinE9npvqmkPviz2aUAVG/awfQ5K1m8dhtHH9aDuy4dzbgjenteUyqYgjAMI2syMds47TamPL0CFKJNseDEuki0ebzXbiSdWcjpwZ/pLuWT3fv4+V9X8/ib6+nVtTM/+eIxTD5xIKHUYkoljCkIw2jH5LvbWaZkYrZxWvVHG72j1t12I273BTzrG00aU5X2+4k2NvH7f67jnldWE9nfyNWfHsK3zzyKnl0LXx6jrTEFYRjtlHx3O8sGL7NNqvJye6Cnw2m34HbfGRcd0+KzeylQp3O9unXm9jnVvLdlD6ce1Ydp54/kyEO6ZyV3KWAKwjDaKfnudpYNbmYbaJ1zILTOlPWD027Ej7ko0wqv33lyOU0Kg3t35aGrxnLm0YeUXNhqppiCMIx2Sr67nSWTienKyWwzfub8VspLoZWSCJUJjU3uasPLiZzOXOSlQBOvk2lS6FHeiXk3fpYunUozbDVT2m+REMPo4OS721mCxMq7ti6CcmDl7bfJDrgrKSXmIxCgV9ew4wOoW+cQEh+XajLKBC8F6mbu2lnfYMohCVMQhtFOyXe3swTpVt7pmL2sljIX00yiB8PamefRtXOn5uilZCq7dmbtzPNYOHVCTqYyN0XpVSMp1MFNSqmYgjCMdsqkMVXMuOiY5hV5rivuBLmYrhK7j0aHRmSpyitIExk4K9CQCFt273O9xknujoz5IAyjHeMnbDNTcsk4dktgC4m0Ul653McPiXv97KVVbNpRjwAicN3pQ3l2aS0f7qhvdU1Vnu5dKtgOwjCMFuRiunJb/TeptlJkQZnIEqgq5eEQoVDMbHT2qH7Mv+l0pkwcwS3njAj03qWC7SAMw2hBphnHyWSyK8jlPumo+WgXtz9fzcI1nzCs30E8ds04xh/ZpyD3LiVE25nNbezYsbpkyZK2FsMwDAeciug5JbAFRd3e/dz18moeXbSO7uVhbjp7GJefNIhOobI2zzpva0TkLVUdm8k1toMwDCNvtNXKvKGxiT8uXs9dL69mZyTKlScfzo1nDaNXt85AcWSdt0dMQRiGkdfVdRCOcy/eWLOV6XNWUrN5F6cc0ZtpF4xkxKE9Wowphqzz9ogpCMPo4LTX1fWGbXv58dx3ean6Iwb0quCBK49n4qhDHctjBB1SW6qYgjCMDk57W13v2dfAr/7+Hg/+431CItx89jCuOfUIysPuGdBBh9SWKqYgDKODUwyraz8mLlXlL8s3MePFd9m8cx+TjuvP1HOP5tCe5Wnnz7ZZUEfHFIRhdHDaenXt2EjoqRVMn1NN3d4o/SsrmHziQP5e8zFL19dx7ICe3H/F8Zxw+MG+72FhrdlhYa6G0Y7Jh3M5qNBUv7KNnznfV6+I7uWd+MEXRnLJ8QMo60Bd3fKFhbkaRgciX87lIFbXmcjmt5HQQV06cenYgVnLZGSOKQjDaKekq7qayQM/XWhqYjdQWxchJEKjqmNrz+Rxqbg5vhPzpePDHfXMXlZrZqECYgrCMNopbk7kxGo9X2GrqbuBxMM8dV4nU5UfmTOpoNoewm9LiYyK9YlImYj0SD/SMIxsmL2slvEz5zNk6lzGz5zv2aTHy4mcST+H789+h6G3vsDgqXMZeusLfH/2Oy3Ou1VoTZ3Xa5yXzIf1SB+F5HQ/I3jSKggR+aOI9BCRbsBKoEZEpgQvmmF0LDLt5OZUDdULp9X792e/w6OL1jev4htVeXTRekb+4MXm+6YLd02cTzcuNay0sUn50+L17NrX0GpsOOTuhLbktsLhZwcxUlV3ApOAF4BBwFcClcowOiCZdnJLbhjkB6fV+5/+tcFx7N5oU7NyShfumjjvNS61mdHitds4/97XufXZdxh5WA9uPntYi8ZHsy4Z7fq5LLmtcPjxQYRFJExMQdynqlERaV+xsYbRDsgmYS3hXB4ydS5ef5RuSWFe9v+EcnJKMnOa1y0ZLVkx1NZFmPHCuzz/9of071nOvZeN4QvHHoaIcP2Eo1rNb8ltbYsfBfFr4ANgBfCaiBwO7AxSKMPoiOSSsOZ2LYAAF5/gHKWULoJoU12kRRisVxSTV7hsZH8jv37tPR549T1U4dtnHsW1pw2lorO7icyS29qerBLlRKSTqrY2HBYAS5Qzgqat+gbkkrCWLoKoqrKChVMntDqe8EG44XadX1SVue98yIwXVlFbF+G8Yw/j1nNHMKBXV8fPYMogOAJJlBORfsBPgP6qeq6IjAROAR7OTkzDKF7asrJpLivmxJgbnljueN7NTHXHpGMAeOxf60ldK+ZqzqnetIPpc1ayeO02jj6sB3ddOppxR/R2HNteK8qWOml3ECLyIvAb4HuqOlpEOgHLVPWYQgiYiu0gjCBxK/uQ60q6UOQify4r+ORr+/UoZ0ifbixa+wm9unbm5rOHM/nEgYQ8ymO09++9PZDNDsJPFFMfVX0SaAKIm5a8g50No51SDJVNc8Ep9NXvTmDSmCoWTp3A3ZOPA+DGJ5anzcWA1uG5H+2s55/vf8Jnj+rLgptO5/JxgzyVA7T/771U8aMg9ohIb4gFSYjIycCOQKUyjDbCzSHcXkIrk0NfEyGjmRTdyzQXA9wT5NZ8vJueXcO+7tvev/dSxU8U03eA54ChIrIQ6AtcEqhUhtFGFGvfgEzMP7m0/My0edDarXtco6dq6yKMnznfl8zF+r13dNIqCFVdKiKnAcOJRczVqGo0cMkMow0oxtDKQjpw/Zp6dtVHuW/BGh55fS0CjjkYwoFKrelkLsbv3fAXxXRVyqHjRQRV/X1AMhlGm5LLCjwICtkSNF0uRlOT8vTSjfzspRq27t7HJScM4Niqnsx4cVULGZ2URjqZi+17N/yZmE5Mel0OnAksBUxBGEYBKKQD94wRfR3zIs4Y0Zel67cz/blqVmzcwZhBlTz81bGMHlgJQI+KcIvVv5vZyZzO7Qs/JqZvJb8XkZ7AH/xMLiLnAL8AQsBDqjrTYczpwD1AGNiqqqf5mdswOgqZZFjnmmy2YNUWx+PPvFXLo4vW069HF+6ePJoLR1e16OqWuvp3C1v1cjpbolxLiuH7yKjcd5y9QOuiKSmISAj4JXAuMBK4LJ5klzymErgfuEBVRwFfykIewyhp/IauZhOBlIrbCj8SbeS6M4Yy/6bT+eKY9C0/Mw23zYfspUSxfB9+yn3PEZHn4j/PAzXAX3zMfRKwRlXfV9X9wOPAhSljLgeeVdX1AKr6cWbiG0bp4zd0NdNqsE64rfD7de/ClIkj6NbFX4+xTMNt3WS/6ckVvnpjlBr5+F3mAz+/7Z8nvW4A1qnqRh/XVQHJtYQ3AuNSxgwjVi3270B34BdOzm8R+SbwTYBBgwb5uLVhFA/5MBX4ceDmw1dx5bhBzPprDU1JHuaKcIhbP39083u/nycTp7ObjG7d60qdYkkc9OODeDXLuZ32oKmBDZ2AE4g5viuAf4rIIlVdnSLDg8CDECu1kaU8hlFwChmimks12Lq9+7nr5dU8umgd5eEQ4VAZOyLRVhVbZy+rZcpTK4g2HXhwT3lqBZDb5/FybCcIKnKrGMnld5lPXBWEiOzCPbxZVTVd69GNwMCk9wOATQ5jtqrqHmIZ268Bo4HVGEZABOH8S8yZWg57z74GR1PBDU8s54YnllNZEea2C0bl5aGXTbJZQ2MTf1y8nrteXs3OSJQrTz6cG88aRq9unR3H3/ZcdbNySBBtUm57rrq5L3U2361Xz4lkOkoUVLEkDroqCFXtnuPcbwJHicgQoBb4MjGfQzJ/Ae6LFwDsTMwEdXeO9zUMV4JY0afOmWwWSUddJJqXFXjy9X4f0G+s2cr0OSup2byLU47ozbQLRjLiUO91X13EOUe2LhLN6btNlb3MpU9FRym9USyJg/48ToCIHEIsDwKAhGPZDVVtEJHrgXnEwlwfUdVqEbk2fv4BVX1XRF4C3iZWDPAhVf13Fp/DMHyRTdJZulWxWy0iv0SblOlzqh3vn8mKPN3Y2ctqmT6nmu17DzzkD+7WmQeuPJ6Jow5FxDsyKR25JvQl+yzcemN0pNIbxZA46CeT+gLgTqA/8DFwOPAuMCrdtar6ArE+1snHHkh5PwuY5V9kw8ieTJ1/flbF+TB7bN8bW4GnPtCnPL2CaGOSvf9p591GOjlnL6vl5qeW09DU8r479+6nPtrkWzn06hpuoWCSj+fTsVosK+iOjp88iB8BJwOrVXUIMYfywkClMoyAyLRqqJ9ww3Rmj15dw1T5MI2khjBOn1PdrBwSRBtju41M5FSNXZOqHAAatPV9vZh2/ijCoZbKJBwSpp0/Ku8VWRPlx9fOPI+FUyeYcmgD/CiIqKp+ApSJSJmqLgCOC1guwwiETBO4/KyKneZMnnva+aNYOHUC90w+rtXD1eteTiv1xPHU3AA3OWvrIlz8qzdc50qMcWL2slrGz5zf4l6TxlQx65LRLfIbZl0ymkljqnLqRWEUJ358EHUichDwGvCYiHxMLB/CMPJKIUoLZGq68BNumDxnahRT8tyJf298cnmr9p6pc6YjObvWS06A//14t+dcAo7mLS+TlVVk7Rj4aTnaDYgQ221cAfQEHovvKgqOtRwtTdyckpk0uylmub4/+x3+9K8NNKoiEnsoJ0eLCvDpoQfzwSeR5ofrtj37iEQd7EIpVFVWMLh3BQvf29bq3MjDurPqw12kmyW1tadbLaXKijDLp52dViaj+Aiq5eg3gf6q2qCqv1PV/24r5WCULsVSWiCVXDu0QUw5PLpofXPYpmpL5QCxHcHC97a1qL3T0KS+/kBr6yKOygFiu4f0Kqa1icrNZJUIZzU6Bn5MTD2AeSKyjVg9padVdXOwYhkdjWIpLeBEruGGf/rXhvSDHIg2Kr26hunauZNnbkC6OfyQat7yMlllm81cDNVJjcxIu0BR1enxSqvXEQt1fVVEXglcMqNDUco9iTN9qCdTtzfaHMlz56WjPZ3c2eLkSPZyLGejtIulOqmRGZmU+/4Y+Aj4BDgkGHGMjkopR8CEckhAa6Ug81yJzM1kNmlMFb26hv3J5INiNSEa3vhJlPu/wGSgL/A08A1VXRm0YEbHopARMIU2dZx8RC9HH0EZePoHEgoyuc5TPhFo4ZhOZdr5o/KWzVzMJkTDHT8+iMOBG1R1edDCGB2bQpQWKGR11cT9Fn+w3fHc5ScPao5sSiUkwoyLjgHwVcQuG9LtBPKptIulOqmRGX7KfU8thCCGUQhyrReUjFsF18G9K3jjvW1prUELVm3hzktHe4bRnjLjbzkrh3BIQGlRhVWI9ZlOR76UdrFUJzUyw3exPsMoBfJl6vCq4OrXFLSpLuKYZBeJNvKzl1axbP12PtxRn5FcyQg0r/qXrNvGY4vWNystJdZneuzhBxckksiS6NonpiCMDkW+TB25VnBNvmfiIZmscDbtqOd3/1xHmbTOmUgQ8gh7TU18mzWvptWOptANeIqhOqmRGZlEMRlGu8fNrLJtz76MQi5zda6Gy6SFecVN4fQo70S4rHUUVDgk3HnpaO6ZfJyv6C9zEhvZ4KogRGSXiOx0+ymkkIaRLxas2uJ4PBJtyiguP1fn6qwvjW5eTUcbm1zNUjsiDRxU3nqjH23U5tW/n0zvUs4zMYIjbUc5EbmdWP7DH4iZNa8Acu02ZxhtgteKOTkuP52t3G+LTDcS9+nVrTO3O5TvTtC/siLt6j+10c6seTXc+MTyFrKbk9jIBj8mpomqer+q7lLVnar6K+DioAUzjCBIt2JOhL2my/hNXrnDgWS4qsoKxg89mHSpcbV1Eb7z5HK++shiGpuUaz4zhPJOLf8cEw9wv6t/r2zlfNSUMjoefpzUjSJyBbE6TApcBuQ/KNswcsRPAly6lX8iiiiZSLSRG55Yzg1PxFKBOoeEbl06Ubc3Sv/KCu6ZfFyrUtmJqqyVXcPsizay16Eqa5PGfAzzbvwsXTqF+FRVT1f5nWTeu7+h+eE/e1ktNz25opXTOtkRbU5iI1P8KIjLgV/Ef5RYN7nLgxTKMDLFbwJc4nVqb2aIrdj9mIz2Nyr749c6tfZMlsOrUQ/ArvoGunQKNV/v1WfhtueqqYscmG/73ii3PvsOS9Zt45m3al0jmjJxRFtBPSMZP8X6PlDVC1W1j6r2VdVJqvpBAWQzDN9kUutn0pgqlv3wbO6ZfFwrk4uf1qCppPouMvFL+HUSTxpTRbcurddzkWgjf/rXBs97lon4cr7ns6CeUzc6o/3hpxbTMOBXQD9V/ZSIHAtcoKp3BC6dYfgkmzBOtxV7Ns7nxH0yWa1n6iR2mztdtdhGVV/lRPKVZV7ociZGcPhxUv8PcCsQBVDVt4EvBymU0bHJZvWZrzDOVGeuXxL38Xu/RK2lTB6YbnP7qRbrp3JqvnIlrHJr6eBHQXRV1cUpx6wntREI2Zo5Mi0X7qWEJo2pau7BcOXJg3zJfcaIvqgqQw/plnZsRTjEnZeOzng17fYZLxs3sNVxJ7we9LOX1VLmomgyVbKWlFc6+HFSbxWRocQr0YvIJcCHgUpllBx+nZ+ZmDlS57z4hCoWrNqS9h7pTCCp844fenDa4nt/rd7M4rXbWL15d4vjTr2ms3X8etUzGnv4wc3H3TrPuT3oE9+H0zXZ5EpY5dbSwY+CuA54EBghIrXAWmLJcobhi0xs0n5Xn05zPrpoPV3DZfSsCLOpLtJs0kgNQfUKB00taldbF2Hbnv3cHQ9lHTJ1rqOi+HjXPrbs3tfquAKL3t9Ok2peooKc/CapCu2MEX155q1a30lxbo71bMxgYJVbSwk/JqZ1qnoWsYZBI1T1M6q6LmC5jBIiE5u0X1+C20Ntb7SJukjU0TzltVImPj5ZOTjJ6iZft84h3HzFjaqBtdl0Msk981YtF59Q5Tspzk0pN6lmvdOxpLzSwM8OYq2IvAQ8AcwPWB6jBMnEJu139enXnp1snpo+p9ozOsmrOmriflMmDue7T7/N/sYDiW+dyiAcKsNP/qhXVFA2OQhOnykSbWTBqi0tqrkmfC5OcwdhErKkvNLAzw5iOPAKMVPTWhG5T0Q+E6xYRimRSYRRrsXnnNhUF2H2stq0SWte4aL9KyvYsG0vL/37I/Y3NhGKV1itrOiEIJrhObsAACAASURBVC0S2PzIk0o2znmvz5R8j3Rzl3I/cCM3/HSUiwBPAk+KSC9iGdWvAunDJowOg9fqN1ObtJ/VZybF8vpXVuQcYjmyfw/OvOtVQiLcfPYwrjn1CMrDIcbPnE9dxDmoT8DRX+Gk3LLJQfD6TMn3SGfiS5xP7opnGdQG+OwHISKnicj9wFKgHLg0UKmMdkW6FWoQNunEnL26hj3HJRRRriGWL6/cDKqxzOXFG7hj7krGz5zv2T1OoXmnkUCIfT+pobXZhIZ6nUtWvm7jkgsTQmwHlfi+TDkY4C+Tei2wnNguYoqq7glcKqNd4Wf1G4RNOjFn8u4ltThel3iFVDc7uxNuvoj9jQfaij66aL2vuZqaYivy2rpIix1FaiRXNn4At2sqK8Itvmu3cW6FCQvZZc4objx3ECISAn6jql9U1T+ZcjCcaOvEqOTEtmnnj0KTcqDrIrGCdmeM6OsrmQzSl67IhMRMlRVh1+io2ctq2bu/tZkqnR/AzXdw2wWjfI3LR3E/o7TxVBCq2gicUSBZjHZKMXUrc9vNLFi1pYWZq7IijEMnz0CorYu4OrETO4lUZ3NlRTitGc6v6c5tnFthQktoMxL4CXN9Q0TuIxbm2ryDUNWlgUlltCvykRiVbZnp1OvczEib6iKtMpEru4bZFYnSkL8NQ8Y4mXkAunXp5Ovz+zXd+S1MaNFLRjJ+FMSn4//ennRMgQkOY40OiFcJCD9kW/3T6TqvyCGnXg3Jm4je3TrzyZ79vmTOB179Jwph5sn192aUPqJ5tLcWgrFjx+qSJUvaWgwjj4y5/a+O8fxVlRUtkr1ScYsiSlUSFeEQMy46hlnzahzH9ywPs+QHZxEOlXnOKRLrApctvbqG6dq5U4uHsZtMTmPtwW3kgoi8papjM7kmbZiriPQTkYdF5MX4+5Ei8vVshTSMZPwmeznhZk5ScLTLu43fWR+NZ0K7O3TvnnycaymNVMJlQjjU0sFREQ4x7fxRzT2mE7WizhjRl3CKM6RMYHd9Q14a9xhGLvgxMf0W+A3wvfj71cT8EQ8HJJNRJBSi/aTfZC8n3MJRQyKtdh6L124jHBKijd5VTr3MLpms9p3mAFqZxJ5YvIHUbtVNGquDlIyFnxptgR8F0UdVnxSRWwFUtUFEfLXbEpFziGVeh4CHVHWmy7gTgUXAZFV92p/oRpBk6hfIVpl45SZ4OUtnL3Pvwdyo2lx36JAeXejfs4JlG+qo7Bpmz76GFkoi4ZT1I/8ZI/q2KuaX2Bl49ZJOMH7m/FY+h2gGNisLPzUKjZ9M6j0i0psD/SBOBnakuyieQ/FL4FxgJHCZiIx0GfdTYF4GchsBk0kF1myb/MxeVuvatS012cvpfm4kspUV2LxzH8s21DFxVD/+OfVMZl0yupX5CUgr/+xltTzzVm0L5SDAxSf4TwDM9QFv4adGofGzg/gO8BwwVEQWEiv7fYmP604C1qjq+wAi8jhwIbAyZdy3gGeAE/0KbQSPV3mG8TPnt1hhZ9vLeNa8GseII4FWyV6p13nVYHKa89+1O6noHHIM93Ra2afK73RPBRas2uIoQ2JHUlsX8awS6xcLPzXaAj/F+paKyGnEqroKUKOqfkpXVgEbkt5vBMYlDxCRKuCLxEJmXRWEiHwT+CbAoEH+WkAaueGVU5Bqbso2k9rtvOId3prNStzrGq/Pme56r8qsCYWSTjmEyoRGD1NTto17DCNX/EQxfQmoUNVqYBLwhIgc72NuJ+tB6l/BPcAt8YxtV1T1QVUdq6pj+/bt6+PWRq44RfMk46eJTjqTiNt5twxfv/Nmek3IpRdz8vFMPmO6HU4q3bt08vzM2TbuMYxc8eOD+IGq7or3gJgI/A74lY/rNgIDk94PADaljBkLPC4iHxAzW90vIpN8zG0EzKQxVVx8QpWrjwBaNtHJpp9ALtdlQro5vZzdyff0K2umO5wdkSgLp06w0hdG0eFHQSSWQucBv1LVvwCdfVz3JnCUiAwRkc7Al4n5MppR1SGqOlhVBwNPA/+pqrN9S28EyoJVWxzt+QkSD65sy3nncl23zv4K7/kxz7g9mEMiWZUsz/SBXhkvWW6Ne4xiw4+TulZEfg2cBfxURLrgQ7HEw2GvJxadFAIeUdVqEbk2fv6BHOTuMOSrRlE2OQxeK+HUB1cm5bxzkW3t1j1c/8el7Nnvz4Tjxzzj1nyoUbWFr8XvZ8ykmRHQnIBnpS+MYsOPgrgUOAf4uarWichhwBQ/k6vqC8ALKcccFYOq/oefOTsS+axR5Oe6VLz6CGTrNM1Wtl31Ue5bsIaH/rHW0aFb5lIGw89qPnHfm55c0crclE2CWvKD3k8U046kSq/Wy9koJvxEMe2N+wjOjSe+LVTVvwYumZFT+Ggm17mt6N2qtHoph3S7g0xla2pSnlm6kZ/Nq2HLrn10DglO6/Ie5WH2NTRlXZl00pgqbnxiueO5bKKm3MJpM20KZBhtiZ8oph8Sc0z3BvoAvxGR7wctmJF9I55sQjKdksQy9RH4SZjLRLal67fzxfsXMuXptxnQq4LvnDWsuatbKjsi0Zzbmgbd18J8DEZ7w4+J6TJgjKrWA4jITGK9qe8IUrBipBC1iZLJpg1lptelW9FnYvJwm+uGJ5Yza15Nc6G6dLJ9tKOen760ij8vq6Vfjy5cOW4Q81d9zF2vrHa9t8bvn8vvJF1fi1x//+ZjMNobfhTEB0A5UB9/3wV4LyiBipV82fUzIdtGPJlcl892oemS0W599h0uPqGKZ96qdZStPtrIw6+v5ZcL1tDQqFx3xlAG9erKbXNW+nL45vo78XqA5+v3bz4Goz3hqiBE5F5iC7N9QLWIvBx//zng9cKIVzxk6w/IhuSVas+KMOXhMur2Rn2vODNZqWa7S3HCK/saWrb+TJbt5rOHUR4O8bm7X2XDtggTR/Xje58fyaDeXR3LYHiR6+/E7QFeyN+/YRQLXjuIRFeet4A/Jx3/e2DSFDH5XGl7kbpSrYtEm/sRBLFSzWaXkolTO5VE68+EbDUf7eL256tZuOYThvU7iMeuGcf4I/u0GJ8pXkoqWwr1+zeMYsJVQajq7wBEpBw4ktju4b2EL6Kjkc+VtheFXqlmahf3Y2px65sAB76vur37uevl1Ty6aB3dy8PcfuEoLj9pEJ1CZa3GZ/rAdyudkQuF+v0bRjHh2nJURDoBPwG+BqwjFvE0gHjzIJ8F+/JOW7UcTX0wQvqQz2wYMnWua4XTtTPPy9t9ssUtVLOyIky3Lp2azWLRxqZWyWwV4RB3TPoUe/c3cOfLq9kZifLpoX1Y8/FuNu+sd1ROs5fVcuMTyz0zup248uRBLFi1xVXpZepwztfvv9CBDoaRIJuWo14mpllAd2CIqu6K36AH8PP4z7ezFbQ9UqgIlGJfqbqZVOoiUeriCV91kdZrh8qKMJePG8SDr71PzeZdnHJEb049qg/3zl/juRuZNKaKG1zyE7x4dNH65tep82bjcM7H778tAh0MIxe8dhD/CwzTlAHxBj+rVPWoAsjXirbaQRSKfKxUgyzP4baDSEd5uIz6aBMDelXw/fOOZuKoQ/nMTxc4zlVVWdGiZajbPXt1DVMfbfLtxE7M6zZf6n3zTVvd1zAgux2EV6KcpiqH+MFGnHuyGHkg2wJ2CXLp7ubnunRlwN2ojzbRo7wTG7dH+NHz7/KX5Zt8O37dEsymnT+quSOcHxLztpXD2RzdRnvDS0GsFJGrUg+KyJXAquBEMiaNqWLh1AmsnXkeC6dOyMj8kEmr0Gyuc1JgveLVSNOxs74BOKB8KsLO//1SzWleSnPSmKq0/SNS53Uz1ymxVX46ZZotQWdqG0a+8fJBXAc8KyJfIxbqqsS6vlUQ6wJnFCF+V6mp5iQ3s5HTfKkhtE5msXR4jd22Z19zqQ+3eybjJ7w2OXTXa3yQfoFsEx8No63wCnOtBcaJyARgFLFAmhdV9W+FEs7IHD9ObidnqeBsN0y9zslHMWlMFTsjUWa+tIq9PstwexGJNjHlqRWAv4f0pDFVLFm3jccWrW/xGRKfqSrFn5IuHDeosGIrtWG0N/xUc50PzC+ALIZPvJzJUyYOZ8pTK4gm1b4Ol0mLVaqTOUmhlZJIrUPkFIHT0NjEJ3v2c+/8NUQbmzioSyd272vI+TNGmzSjh7RTc6OEcnByACcUm1tYcW1dhNnLapk+p5rte2NRWZUVYW67YFROD3QrtWG0J/x0lDOKCF/O5NQ8sZT3bmaoxAPVyTnu5qO45dl3mPHiKsYNOZh5N3yWPT6Vg59Utkyct9k6gL3s/zc9taJZOUAsfHfKUysC81EYRrHhp1ifEQDZhqKmy7SeNa+GaEpJ7GijNjubZ82rcQ1B8wq3dPNRNDYpv736ROr2RvnKw4t9hbdVuiTSpZKJ8zbb/JEpE4e7JuI5NSbKdGdjGO0Z20G0AdmGokL6lbLb+cQ93B70Xs7S789+x1We/j3Lqdsb9Zw7lX0N6ZVDqlksHdn2Wpg0pirjmG0LSzU6CraDCBinnUIu9ZbSrZS92oS6RfmkOnFT5U/OSk5GgO+eM4Lbnqv2HcHkJUeCbGz9Tm0+k8N0veaqyrDek4WlGh0FUxAB4ubYdXtA+lmZpguVdDvvdk8BzyzeO+audD2XWHk7ldZwu5dbb+Z81JpKKIFMy1k4fWfhMqGJ1mamTHc2htGeMRNTgLjtFNyqjfpZmabLtHY775ZM5nbP2roI1/9xKVt373eVpaqyIm0CXgK3MNp0cmRKNomCTt/ZrC+N5s4vjW6RBFhZEWbWl0ab/8HoMNgOIkDcdgSNqoRD0sKZnEvC1JJ121qZsZx2BX7aaaaaWsJl0iJkNoFwwMHrBy/lkM9ksWyjmdzCT4NWBlbd1ShmbAcRIJ6rYo0Vm8u03pKTg/vRRevTOry9dh6zl9Uy9Zm3He3w0SYlVNZyxyPAFScPYtKYqpxX/pnWmkpHeypnkUuwgmEUAttBBIhXSYdok9K1cyeW/fDsjOZ0MqGkEok2Mn1OdfP4xOr0jBF9Hcf/eO671Dc0uc7X1KRUVVa4JuZl068BYg7rfK+Y21M5C2tjahQ7toMIkMSq3Y1swiX9XrN9b5QpT6/w3GlMfeZtLntwEVt27/OcS4k9ePvHlcSseTXNq9xswkQTNKrmfcWcazXcQmLVXY1ix3YQAZNIXstXE6BMWnCmJsylUt/QxD/f/4RuXULs2ee+KxG8I4PcwkRDIq5RSwmCWDG3l3IWxd4cyjBsB1EAsk3imr2slvEz5zNk6tzmMtRBmEp+POmYtD0evCKD3D7fZeMGEi5LX1Sjo66Ys/1/YRiFwhREDjg9wJ3Ixuzh5sAEfPdf8ENVZUUL+Zxw2wMkHuxun++OScdwUHn6TWqZSNrvsBRpT+Ywo2Pi2nK0WCmWlqP5amLvhld7Sj/9D/yQkBdaOrP37GvwlfzWq2s4rZPdrVpqOpny2V7VQkkNI/8tRw0Psu3c5hc3P0NtXaTVyrOyItwqFDWVrp1DXHx8VavVKtBqp+I3M3p3fUPaFb+XPd0pYTCT79BPmKiFkhpG9piTOkuCjkBxc/AmHqrJjtjxM+e7PtRF4MYzh/FfZx3leH78zPlZ70T8VDZ1CzudcdExrkl2fr9DP2GiFkpqGNljCiJLgo5AcYv+aVT13S4U4M4vjeai4we4ns9VofnJUAbnLmq5Rnf5UdIWSmoY2WMKIkuCTshyCx3t1TXcKuTUjZAIFx0/wNMG76ZgenUN07Vzp+Zr9u5vaNE8J4Hf+lFOq/Vcv0M/StpCSQ0je8wHkSVBR6C4hUCqtg45dSOx2/CywbvdZ9r5o1g4dQJrZ57HwqkTmHb+qLyHZOb6HfoJE7VQUsPIHotiKmKcVv6ZlLVIhK26RUMlCvr5jfLJRzRQviOKLIrJMPyRTRSTmZjaCD8PLSfTzE9fXMWHO+t93cOr2mqyDd5v5nGuGcpu/TEScwdFe8msLjSmOI10mImpDcgm9LKhsYk//PMDdtT7C0GtrAh7VlttCxt8vkODLYQ1e+y7M/wQqIIQkXNEpEZE1ojIVIfzV4jI2/GfN0RkdJDyZIPfbOlMyPRB+caarZz336/zg79UM3pAJd+dOLzZbt+ra7hVOYuKcIjbLhgFFJcNPt8RRUHnopQy9t0ZfgjMxCQiIeCXwOeAjcCbIvKcqib3sFwLnKaq20XkXOBBYFxQMmVKUCYRvw/KDdv28uO57/JS9UcM6FXBA1cez8RRhyIi/OcZR7aQ081U4BVmWmjyHVFkIazZY9+d4YcgfRAnAWtU9X0AEXkcuBBoVhCq+kbS+EWAe8B+GxBUkpXbg1KJJa7914Qj2bA9woP/eJ+QCDefPYxrTj2CcpeCeuls7MVig893aLCFsGaPfXeGH4I0MVUBG5Leb4wfc+PrwItOJ0TkmyKyRESWbNmyJY8iehPUKsvJ7JOgti7CLc++w30L1vD5Tx3KgptP5/oJR1EeDgVi7iok+Q4NLqT5rL1/96kUk+nRKF6C3EE4FQdyjNAUkTOIKYjPOJ1X1QeJmZ8YO3ZsweJyg1plJZt93BLdyoC/LN/Emx9sb/6jTTV33fjEcpas28Ydk9ybEhUb+dzNFMp81lbRV0FSTKZHo3gJLA9CRE4BblPVifH3twKo6oyUcccCfwbOVdXV6eYtZB5EthVbk30ClV3DqMKOSNTxj3Dw1Llp5agIhygPlzlmMif6Qy9YtcXXH7qTvwLsQeGFV2XdRC6JYRQ7xZYH8SZwlIgMAWqBLwOXJw8QkUHAs8BX/CiHQpPNKitVqSQ/1JNXnp8/5jB+s3Atgnu/hQSRaKNr9rQCjy1a3zyH1+rWaSU85ekVoLHCe+mu76iYQ9foqASmIFS1QUSuB+YBIeARVa0WkWvj5x8Afgj0Bu6XWJXShkw1nB9ySQjK1CTi5NhOJhJt5Ht/foebnlpBY5MSLhNUlYYcNnKpl7o50p1kc2pLatVOW2IOXaOjEmgmtaq+ALyQcuyBpNfXANcEKUOh7cd+VpV79h94SEeblHBIqOzciR2RKGU++jhnK0cmK15bHR8g6MKMhlGslHwmdaETgrJZVUYblbq4j+KycQMzutatTZCTHJnI5jS21CJ5/GKtQY2OSsnXYiq0/TiXdqC1dRGeeauWbp1DLXYZblRVVnDGiL4881atr9Wtk2zhkLTwQbhdX4qRPJlQLLkkhlFISl5BFNp+nHiI3DF3JVt37wcgVCY0NilVHn0VEkSijVRWhKkIe5f1To6gGXv4wb58LG5Od6djfvwX5qswjNKm5BVEoe3HtXURXnl3M1t376d/z3Ju/fzRfOHYw4g74R1DZ1PZEYly9+TjmvMkUiOdUuXPZHXrNjbd9RbJYxgdj5JXEIVKCIrsb+TXr73HA6++hyp8+8yjuPa0oVR0bpmt6idJrn9lRYsHeTGUZbZIHsPoeFjDoBxRVea+8yEzXlhFbV2E4wZW8uGOCB/v3OcraS2bRLy2oD3JahhGa4otUa7kqd60g+lzVrJ47TaOPqwHF42p4qHX1/p25LancgftSVbDMPKD7SCy4JPd+7jz5dU8vng9PSvC3DxxOF8+cRCf/dkCK8lgGEZRYjuIgIk2NvH7f67jnldWE9nfyH98egjfPvMoenYNA+bINQyjtDAF4ZNXV2/h9jnVvLdlD6ce1Ydp54/kyEO6txjT0R25xeBMNwwjf5R8JnWurN26h2t+9yZffWQxDU3KQ1eN5fdfO6mVcgDnGvtCzBdRzJnH+ciQth7HhlF62A7ChV31Ue5bsIZHXl9L51AZU88dwdXjB9Olk3OjH2gdwpqcv1Csmcf5ypC2RDrDKD1MQaTQ1KQ8s3QjP5tXw5Zd+7jkhAF895zhHNK93PUaJ9OKU55DJNrITU+u4MYnlheNCSZfD3bzvxhG6WEKIoml67cz/blqVmzcwZhBlTx01VhGD6z0vMZtBe6WKZ2o1FosO4p8Pdg7uv/FMEoR80EAm3fWc+MTy7no/jf4aGc9d08ezTPXfjqtcgC47blqxxV4SNzqrLYcF1RVWb+4PcAzfbBbj2PDKD069A6iPtrIw6+v5ZcL1tDQqFx3xlD+8/Qj6dbF39cye1ktdRHnwnuNqlSEQ2mrura1CSZftaoskc4wSo8OqSBUlXnVm/nxCyvZsC3CxFH9OHHwwfxm4Qfcv+A93w83r9V/SIQZFx3T/MB0awTU1iaYfD7YrSS2YZQWHU5B1Hy0i9ufr2bhmk8Y1u8gHrtmHFt27cs4kmf2slrXYntAK2XQo6ITu+sb0vZdaAvswW4YhhMdRkHU7d3PXS+v5tFF6+heHmb6BaO4YtwgOoXKGD9zfkaRPAnHtBeVFeEWSmf73mistWhFmB3x7nFmgjEMo5gpeQXR0NjEnxav586XV7MzEuXKkw/nxrOG0atb5+YxmUbyOIWGJhMuE0RaN/yJNirdunRi+bSzs/gk+cWyng3DSEdJK4g31mxl+pyV1GzexSlH9GbaBSMZcWiPVuMyDdH0cixXVoS57YJR3PjEcsfztXURhkyd26YP5Y7ePtQwDH+UZJjrhm17ufYPb3H5Q/9iz/4GHrjyeP74jXGOygEyD9F0UxxVlRUsn3Y2k8ZUeTqf27oUhVdynGEYRoKSUhB79jXw83k1nHnXq7y6egs3nz2MV75zGud86kDLTycmjalixkXHUFVZgRB70Hs1wvGjUJzGpNJWD2XLejYMww8lYWJSVf6yfBMzXnyXzTv3Mem4/kw992gO7eleHiOVTPs6g3doaOoYt64bbfFQtqxnwzD80O4VxNsb67jtuWqWrq/j2AE9uf+K4znh8IMDv68fhZI8ZvzM+UXzUM5XcpxhGKVNu1UQH++qZ9ZLNTy9dCO9u3XhZ5ccyyXHD6CsLH2Ji7agmB7KlvVsGIYf2p2CUIVfv/oe985fw76GRr556hFcP+FIupeH21o0T4rtoWzJcYZhpKPd9aTuPnC49r7iLs4ccQjfO+9ojuh7UFuLZBiGUfR0iJ7UAvz26hM5ffghbS2KYRhGSdPuwlyPOqS7KQfDMIwC0O4UhI82C4ZhGEYeaHcKwjAMwygMpiAMwzAMR9pdFJOIbAHW5ThNH2BrHsQJgmKWDYpbvmKWDYpbvmKWDUy+XEjIdriq9s3kwnanIPKBiCzJNNyrUBSzbFDc8hWzbFDc8hWzbGDy5UIuspmJyTAMw3DEFIRhGIbhSEdVEA+2tQAeFLNsUNzyFbNsUNzyFbNsYPLlQtaydUgfhGEYhpGejrqDMAzDMNJgCsIwDMNwpGQVhIicIyI1IrJGRKY6nB8hIv8UkX0icnMRyneFiLwd/3lDREYXkWwXxuVaLiJLROQzhZLNj3xJ404UkUYRuaRYZBOR00VkR/y7Wy4iPyyUbH7kS5JxuYhUi8irxSSfiExJ+u7+Hf/9Bt8hzJ9sPUVkjoisiH93VxdCrgzk6yUif47/7S4WkU+lnVRVS+4HCAHvAUcAnYEVwMiUMYcAJwI/Bm4uQvk+DfSKvz4X+FcRyXYQB/xXxwKrium7Sxo3H3gBuKRYZANOB54v5P+3DOWrBFYCg+LvDykm+VLGnw/MLxbZgP8H/DT+ui+wDehcRPLNAqbFX48A/pZu3lLdQZwErFHV91V1P/A4cGHyAFX9WFXfBKJFKt8bqro9/nYRMKCIZNut8f9lQDdwbbndJvLF+RbwDPBxEcrWVviR73LgWVVdD7G/kyKTL5nLgD8VRDJ/sinQXUSE2CJqG9BQRPKNBP4GoKqrgMEi0s9r0lJVEFXAhqT3G+PHioVM5fs68GKgEh3Al2wi8kURWQXMBb5WINnAh3wiUgV8EXiggHKB/9/rKXEzxIsiMqowogH+5BsG9BKRv4vIWyJyVcGky+DvQkS6AucQWwQUAj+y3QccDWwC3gG+rapNhRHPl3wrgIsAROQk4HDSLDxLVUE4FQUvpnhe3/KJyBnEFMQtgUqUdEuHY61kU9U/q+oIYBLwo8ClOoAf+e4BblHVRoexQeJHtqXEauKMBu4FZgcu1QH8yNcJOAE4D5gI/EBEhgUtWJxM/m7PBxaq6rYA5UnGj2wTgeVAf+A44D4R6RG0YHH8yDeTmPJfTmyHvYw0O5x211HOJxuBgUnvBxDT6sWCL/lE5FjgIeBcVf2kmGRLoKqvichQEemjqoUoVuZHvrHA47GdPn2Az4tIg6oG/TBOK5uq7kx6/YKI3F9k391GYKuq7gH2iMhrwGhgdZHIl+DLFM68BP5kuxqYGTe/rhGRtcRs/YuLQb74/72rAeJmsLXxH3cK4UAp9A8xxfc+MIQDDptRLmNvo/BO6rTyAYOANcCni1C2IzngpD4eqE28Lwb5Usb/lsI5qf18d4cmfXcnAeuL6bsjZiL5W3xsV+DfwKeKRb74uJ7E7PvdCiFXBt/dr4Db4q/7xf8u+hSRfJXEnebAN4Dfp5u3JHcQqtogItcD84h59x9R1WoRuTZ+/gERORRYAvQAmkTkBmJe/52uExdQPuCHQG/g/vhKuEELUC3Sp2wXA1eJSBSIAJM1/r+uSORrE3zKdgnwf0Wkgdh39+Vi+u5U9V0ReQl4G2gCHlLVfxeLfPGhXwT+qrFdTkHwKduPgN+KyDvETD63aGF2hn7lOxr4vYg0EotU+3q6ea3UhmEYhuFIqTqpDcMwjBwxBWEYhmE4YgrCMAzDcMQUhGEYhuGIKQjDMAzDEVMQRskRr/C5POlnsIi8ET83WEQuTxp7nIh8Pot7/F1Ecg47ztc8hhEEpiCMUiSiqscl/Xygqp+OnxtMrCBdguOAjBWEYXQETEEYHQIR2R1/ORM4Nb6zuAW4HZgcfz9ZRLqJyCMi8qaILBORC+PXV4jI4/Fa+k8AFQ73OFdEnkx6f7qIzIm//pXEemdUi8j0NDIiIpeIyG/jr/uKyDNxmd4UkfHx46cl7ZKWiUj3fHxXEQMU4AAAAiBJREFUhpGgJDOpjQ5PRbwgGcBaVf1i0rmpxEqrfAFARDYDY1X1+vj7nxDrMfA1EakEFovIK8D/Afaq6rHxGllLHe77MvBrEekWz/KdDDwRP/c9Vd0mIiHgbyJyrKq+7fPz/AK4W1VfF5FBxLJljwZuBq5T1YUichBQ73M+w/CFKQijFImo6nFZXns2cIEc6DJYTqwu1meB/wZQ1bdFpNXDPV7u4CXgfBF5mlhF1O/GT18qIt8k9jd3GLHa/H4VxFnAyHjJFYAe8d3CQuAuEXmMWA+HjZl9VMPwxhSEYbREgItVtabFwdjD2U9dmieA64gVk3tTVXeJyBBiq/0TVXV73HRU7nBt8vzJ58uAU1Q1kjJ+pojMJeZDWSQiZ2msEYxh5AXzQRgdjV1Ad4/384BvxcshIyJj4sdfA66IH/sUsVarTvydWIXbb3DAvNQD2APskFgHr3Ndrt0sIkeLSBmxgnQJ/gpcn3gjIsfF/x2qqu+o6k+JFZ4c4TKvYWSFKQijo/E20CCxjm43AguImW+Wi8hkYhU5w8DbIvJvDjRD+hVwUNy09F1cavxrrEnR88SUwPPxYyuINWepBh4hZhpyYmr8mvnAh0nH/wsYG3eQrwSujR+/QUT+LSIriFWGLVTXQaODYNVcDcMwDEdsB2EYhmE4YgrCMAzDcMQUhGEYhuGIKQjDMAzDEVMQhmEYhiOmIAzDMAxHTEEYhmEYjvx/h3U9Vh+rVTEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(yhat, y)\n",
    "line_fit = sm.OLS(y, sm.add_constant(yhat, prepend=True)).fit()\n",
    "abline_plot(model_results=line_fit, ax=ax)\n",
    "\n",
    "\n",
    "ax.set_title('Model Fit Plot')\n",
    "ax.set_ylabel('Observed values')\n",
    "ax.set_xlabel('Fitted values');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot yhat vs. Pearson residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Fitted values')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2df5wdZXXwv2c3l2QTfmyQaGFLDFANikgiq6KxFqySggUjqFSxWmtLW7WWVFND61uCYokipW/1rYrVaivFoNAVjK9BTcBXLGrSbAjRpCq/FyrBZBGSJdkf5/1jZjazd+fHM3Pnzp1793w/n3yyd+78OPPcmXOe5zznnEdUFcMwDMPIQlerBTAMwzDaDzMehmEYRmbMeBiGYRiZMeNhGIZhZMaMh2EYhpEZMx6GYRhGZsx4GC1FRC4WkdsSvr9dRP6ogOucKSIPN3qeKiAia0TkS62Ww4VOandjKmY8DGdE5H4RGRGRp0Tkf0TkCyJyeCPnVNXrVfXsomTMi4ioiOzz7+2XIvIdEbmo1XK1AyLyByIy7rfdr0RkUER+N8d5viAiVzZDRqN4zHgYWTlPVQ8HlgBLgctaLE+RnObf22LgC8AnReTy1orUNvyn33a9wOeAG0Xk6BbLZDQRMx5GLlT1f4ANeEYEABE5Q0S+LyLDIrJNRM4MffcHInKviDwpIveJyMWh7d8L7fcaEdkpIk+IyCcBCX03xV0jIov8EcMs//M7ROQn/jXuFZE/yXlvj6vqvwF/BlwmIs/wz3+UiHxORB4VkSERuVJEukP3caeIfMKXfaeI/HZI1rRjvyciHxeRvX77nBM69gQRucO/r28Bx4TlTWn320Xkw75sT4rIbSJyTOj7V4SOfUhE/sDfPtuX50ER+YWIfFpEehzabgL4PNADnFj/vYg8z5dpWER2iMj5/vZLgIuBv/JHMLemXctoLWY8jFyIyK8D5wA/8z/3AeuBK4GjgfcDN4nIAhGZB/wjcI6qHgG8HBiMOOcxwE3AB/EU5M+BZRnEegz4XeBI4B3AtSLyolw36PE1YBbwEv/zF4Ex4DfwRl1nA+H5mJcC9/qyXw7cHOp9uxy7yz/2Y8DnRCQwnP8ObPG/+zDw9uCgpHYPnfsteO3xTOAwfx9EZCHwf4FPAAvwOgLB7/JR4Ln+tt8A+oC/TWsw35D/EfAU8NO672rArcBtvix/DlwvIotV9TrgeuBjqnq4qp6Xdi2jtZjxMLIyICJPAg/hKevArfNW4Buq+g1VnVDVbwGbgXP97yeAF4hIj6o+qqo7Is59LvBjVf2qqo4C/wD8j6tgqrpeVX+uHnfgKanfzHWX3vlGgceBo0XkWXjG8lJV3aeqjwHXAr8XOuQx4B9UdVRV1+EZg9c6HvuAqn5WVcfxDM2xwLN8Bf9i4H+p6gFV/S6eAg5Ia3eAf1HV/1bVEeBGDo0WLwa+rao3+DL/UlUHfaP1x8BKVd2jqk8Cf1cnbz1niMgw3u/1ZuD1qvpE/T7A4cBaVT2oqhuBr/v7G23GrFYLYLQdK1T12yLyW3g94mOAYeDZwBtFJNxjrAGbVHWfP/n8frwe9Z3A+1R1Z925j8MzSgCoqorIQzjiu3oux+sxdwFzge2Z7/DQ+Wp4PfI9ePdXAx49NCCgKywvMKRTK40+gHdPLsdOGklV3e/vdzhe++5V1X115z3e/zu23aPODez3z4t/jp9H3PoCvLbbEpJXgO6IfQPuUtVXJHwP/u/ru7YCHsAb1RhthhkPIxeqeoeIfAH4OLACTxH+m6r+ccz+G4ANvt/8SuCzTB8VPMohpYjfAz4+9P0+PKUW8GuhfWfjubzeBnxNVUdFZIDQnEkOXofnavohnrvnAHCMqo7F7N8nIhIyIAuBW/DaJu3YOB4F5ovIvJABWQgE10hs9xQe4pBLLszjwAhwiqoO5ThvHI8Ax4tIV8iALAT+2//bSny3Eea2MhrhH4DXiMgS4EvAeSKyXES6RWSOeDH+vy4izxKR8/25jwN4/vDxiPOtB04RkQt83/l7CRkIPH/8K0VkoYgcxdRIr8OA2cBuYMwfheQKARaRo8Wb0P8/wEd9d86jeG6wa0TkSBHpEpGT/BFYwDOB94pITUTeCDwPz6XkcmwkqvoAnhvqChE5TEReAYRHGbHt7nCr1wOvFpE3icgsEXmGiCzxFftn8eaMnum3SZ+ILHc4ZxI/wOsA/JXfRmf69/Jl//tfEDHJblQTMx5GblR1N/CveP74h/B66n+Np8AfAlbhPWNdwPvwep57gN8C3hVxvseBNwJrgV8CzwHuDH3/LWAdcDfeBPLXQ989iWdsbgT24k0S35LxlraJyFN4QQB/hOfzD08Svw3PSP3Yv8ZX8eYmAn7gy/w48BHgDar6S8djk3gL3oT6Hjy33L8GX6S0eyKq+iDe3Mj7/HMPAqf5X38Arx3uEpFfAd/GC2HOjaoeBM7Hm/95HPgn4G0h9+XngOf7kVgDjVzLaD5ii0EZRuP4Ia5/5OD3N4yOwEYehmEYRmbMeBiGYRiZMbeVYRiGkRkbeRiGYRiZ6Yg8j2OOOUYXLVrUajEMwzDaii1btjyuqgvS95xORxiPRYsWsXnz5laLYRiG0VaIyAN5jzW3lWEYhpEZMx6GYRhGZsx4GIZhGJkx42EYhmFkxoyHYRiGkZmOiLYy2oOBrUNcvWEXjwyPcFxvD6uWL2bFUlvKwTDaETMeRikMbB3ispu3MzLqVWIfGh7hspu9dZrMgBhG+2FuK6MUrt6wa9JwBIyMjnP1hl0tksgwjEYw42GUwiPDI5m2G4ZRbcx4GKVwXG9Ppu2GYVQbMx5GKaxavpieWveUbT21blYtb2hxOsMwWkRLjYeIfF5EHhORe0Lb1ojIkIgM+v/ObaWMRjGsWNrHVRecSl9vDwL09fZw1QWn2mS5YbQprY62+gLwSUJrMvtcq6ofL18co5msWNpnxsIwOoSWjjxU9bvAnlbKYBiGYWSnqnMe7xGRu3231vyoHUTkEhHZLCKbd+/eXbZ8hmEYM5oqGo9PAScBS4BHgWuidlLV61S1X1X7FyzItZaJYRiGkZPKGQ9V/YWqjqvqBPBZ4CWtlskwDMOYSuWMh4gcG/r4euCeuH0NwzCM1tDSaCsRuQE4EzhGRB4GLgfOFJElgAL3A3/SMgENwzCMSFpqPFT1zRGbP1e6IIZhGEYmKue2MgzDMKqPGQ/DMAwjM2Y8DMMwjMyY8TAMwzAyY8bDMAzDyIwZD8MwDCMzZjwMwzCMzJjxMAzDMDJjxsMwDMPIjBkPwzAMIzNmPAzDMIzMmPEwDMMwMmPGwzAMw8iMGQ/DMAwjM2Y8DMMwjMyY8TAMwzAy09LFoAzDMJIY2DrE1Rt28cjwCMf19rBq+WJWLO1rtVgGZjwMw6goA1uHuOzm7YyMjgMwNDzCZTdvBzADUgHMbWUYRiW5esOuScMRMDI6ztUbdrVIIiOMGQ/DMCrJI8MjmbYb5WJuqzbG/MFGJ3Ncbw9DEYbiuN6eFkhj1GMjjzYl8AcPDY+gHPIHD2wdarVohlEIq5YvpqfWPWVbT62bVcsXt0giI4wZjzbF/MFGp7NiaR9XXXAqfb09CNDX28NVF5xqo+uKYG6rNsX8wcZMYMXSPjMWFcVGHm1KnN/X/MGGYZSBGY82xfzBhmG0EnNbtSnBUN6irQzDaAVmPNoY8wcbhtEqzG1lGIZhZMZGHoYRgSVgGkYyZjwMow4ryGcY6bTUbSUinxeRx0TkntC2o0XkWyLyU///+a2U0Zh5WAKmYaTT6jmPLwC/U7dtNfAdVX0O8B3/s2GUhiVgGkY6LXVbqep3RWRR3ebXAWf6f38RuB34QGlClYj51auJFeQzjHRaPfKI4lmq+iiA//8zo3YSkUtEZLOIbN69e3epAhaBFTasLpaAaRjpVNF4OKGq16lqv6r2L1iwoNXiZMb86tXFCvIZRjpVjLb6hYgcq6qPisixwGOtFqgZmF+92lgCZmOYS7bzqeLI4xbg7f7fbwe+1kJZmoYVNjQ6FXPJzgxaHap7A/CfwGIReVhE3gmsBV4jIj8FXuN/7jjMr250KuaSnRm0OtrqzTFf/XapgrQAK2xodCrmkp0ZVHHOY8ZgfnWjE7FQ55lBFec8DMNoY8wlOzOwkYdhGIViLtmZgRkPo22w8M/2wVyynU+q20pE/kJEjhSPz4nIf4nI2WUIZxgBFv5pGNXCZc7jD1X1V8DZwALgHXRo+KxRXSz80zCqhYvxEP//c4F/UdVtoW2GUQoW/mkY1cLFeGwRkdvwjMcGETkCmGiuWIYxFcvIN4xq4WI83om3psaLVXU/cBie68owSsPCPw2jWsRGW4nIi+o2nShi3iqjNVj4p2FUi6RQ3WsSvlPgVQXLYhiJWPinYVSHWOOhqmeVKYhhGIbRPjglCYrIC4DnA3OCbar6r80SyjCM9iAucdMSOjufVOMhIpfjrSn+fOAbwDnA9wAzHoYxgwkSN4P8myBxc/MDe7hpy9C07YAZkA7CJdrqDXgl0v9HVd8BnAbMbqpUhmFUnrjEzRt+8JAldM4AXNxWI6o6ISJjInIk3rKwJzZZLqMDMVdGZxGXoDmumml/oz1xMR6bRaQX+CywBXgK+GFTpTI6jjgXB5gro12JW7ejWyTSgFhCZ2eR6rZS1Xep6rCqfhpvWdi3++4rw3DGalN1HnGJm29+6fGW0DkDcJkwf2XUNlX9bnNEMjoRq03VeSQlbvY/+2hzUXY4Lm6rVaG/5wAvwXNfWZJgDmaq39+WJm0vXJ/TuMRNS+jsfFKNh6qeF/4sIscDH2uaRB3MTPb7r1q+eMq9Q+e5MjqlYzCTn1PDnTxrmD8MvKBoQWYCM9nvv2JpH1ddcCp9vT0I0Nfbw1UXnNoxyqiTFquayc+p4Y7LnMcn8GpZgWdslgDbmilUpzLT/f6d7MpIUrjtds8z/Tk13HAK1Q39PQbcoKp3Nkmejsb8/p1LJylce04NF1xCdb8Y+ne9GY782JoUnUsnLVZlz6nhQtJ6Hts55K6ahqq+sCkSdTBVXpOiUyZ7W0WVAwKy/rZVfk6N6iAaU0pARJ7t//lu//9/8/+/GNivqh9qsmzO9Pf36+bNm9N3NCKpj64BT/F10oR2GVTRANtvayQhIltUtT/XsXHGI3TyO1V1Wdq2VmLGozGWrd0Y6ePu6+3hztWdm85TtrJvhXGZqb+t4UYjxsNlwnyeiLxCVb/nX+zlwLw8FzOqSSdN9rpSdi5Dq3InZuJva5SDS57HO4H/IyL3i8j9wD8Bf9hUqYxSadZk78DWIZat3cgJq9ezbO3GSuU8lJ3L0KrciU6ayDeqhUuG+RbgNL8cu6jqE80XC3xD9SQwDozlHVp1Cs10eTRjsrfqWcpl98iLuF6eZ+CskxfwpbsejNxuGI2QFG31VlX9koj8Zd12AFT175ssG8BZqvp4CdepNM1WxM2Irql60lzZuQyNXi/vM7Bp5+5M22cCVQxsaEeSRh7BvMYRZQhixBOniK+4dUdhD33R2d9V97WXHVrb6PWuuHVHLmNc9d+hbKo+Im4nYo2Hqn7G//+K8sSZKgJwm4go8BlVvS78pYhcAlwCsHDhwhaIVx5xL/re/aMMbB1KfOhb1cuqQpZy0r0H/6+5ZQfDI6MAzKnlKfXmRiOju4GtQ+zdPxr5XZoRqMLvUCWqPiJuJ1LfFhH5mIgcKSI1EfmOiDwuIm8tQbZlqvoi4Bzg3fXriqjqdarar6r9CxZ0tv826UVPmnBtZbG+Vmcpu977gbGJyb/37h9tavusWNrHnatfxX1rX8udq1/lrKySfuM0I9Dq36Fq2EisOFy6Wmer6q+A38WrqPtcpq7x0RRU9RH//8eA/8BbR2RGkvSiJz30rayO2uoqui733i7VY5N+4zQj4Po7VDkyrkgs+qw4XPI8av7/5+IVRdwTTJo3CxGZB3Sp6pP+32cDlcloL5sVS/umuFfCJD30re5ltbKKrsu9N7t9wm6zo3pqiMDw/tHM7sM411NvT62Q9q36PECRrtcql5FpN1xGHreKyE6gH/iOiCwAnm6uWDwL+J6IbAN+CKxX1W82+ZqVZs35p2R2P8zkXpbLvTezferdZsMjo+zdP5rLfRjnelpz/imZ5Yi6dpVHYEW7Xls9Iu4kXKrqrgZeBvSr6iiwH3hdM4VS1XtV9TT/3ymq+pFmXq8dyPPQz2R/t8u9N7N9ohRymCzKuRGF52IYsozAynZvNcOw5Z17MqbishjUXLziiAvxopuOAxYDX2+uaEY9Wd1A7V4dNXBXDA2P0C3CuCp9jvfgcu/NbB8X11cW91heF6CLYUiLyAr/DsKhUttluLda7Xo14nGZ8/gXYAvwcv/zw8BXMONRaer9xNdetKRtjAZM98OP+wU8sygsF4XbrHmZOIVcv0+zcQnVTZoHqP8d6suoNjvM1UKNq4vLnMdJqvoxYBRAVUeA5s6YGw3RCetpJ7l9quKPh3g3TpRLLExZ7kMX11ySWyzN/QbNHQXMZNdr1XEZeRwUkR78ToeInAQcaKpURkM0OxGqjMTDNIXUbLeFyz26RCkVEW3VCK6uubgRmEs7N3MU0O6u107GxXhcDnwTOF5ErgeWAX/QTKGMxmimn7issM40t08zFZbrPaYZ6WaHKrsa8Xo5gtGSizJO+x3KGAW0MuTbiMcl2upbwAV4BuMGvJDd+5orltEIzQxBLSusM8nt02yFFXePl64bnOKaapaRdoloyuuazHpc1O8Q+KxnepjrTEmsjCNx5CEiLwP6gO+q6noReSHwj8BvAseXIJ+Rg2YmQjVTYdb3oq+64NTc0VaNkHQvQ8MjrPrKNqA5k7lFjXriyHpcmW6jIt2hzXatVj2xsgySSrJfjVeSZBD4gIh8HXgX8HfYYlCFU+TD3owXPpAvbtHi+rDOLNeNexGvuuDUliyVmuaqGZ1Q1tyygzXnn1K4kXZV7nmNeJ7jynAbFamMy1DsVmAxeeTxWmCpqj4tIvOBR4AXqupPyxGtPShC6TfjYY964fPKWi9fPXFhna73UbUXMWrkVs/wyGghRrr+N4kzWvXKPe+op6qhr0U+A1nPlee9sPyTZOMxoqpPA6jqXhHZZYZjKkUp/aKVZ9TLAOSWNSlcM+xGWrZ2Y0esORE2Cmm5Gq69ctffJJyEF6Zeued1TVa1tlORz0DWjPk870VVjXCZJBmPk0TkltDnReHPqnp+88RqD4pS+kW+OFEvw8p1g8ypdTEyOjFl35HRcf7yxkEuXTcIQE+tizm17mnhpHFyCExxK+W9j9jCf3NrzlFBeUhb72PF0j6Wfui2yLU05s+tTduWdJ0oBTV7Vte050dhmgGJUu55Rz1VDX0tUhlnOVfed7iqRrhMkoxHff2qa5opSDtSlNIv8sWJehkUphmOgImQlhoZnZjcL9wDc5Uv731EvYi1buGpp8cmFXfRfmvXHufl553Cqq9uY3T8UEPVuoXLz0svShgQp6DiRnOKN6LLGoLrSvi4wICuXDfYUkNSpDLOcq6873BVjXCZJK0keEeZgrQjRSn9Il+cIl09QQ/MVb689xH1Iu47MDatBP3I6Djvu3HblGPysuYWt2Vdi1ASWX+Tvt6e1ECBqs615ZWrSGWc5VyNvMMzPf/EJUnQiKEopV/ki+NSUykLjwyPOMkXKI2R0fFcYbX1L+IJq9dH7jeuWoiCi1obBYhsuzQlkaYw436T+XNrPD06kfn5yaP0o2RsxlxbI8aoSGXsei5zP+XHjEcDFN1byhodFfX9quWLWbluMHLSda4/7xEXbhtF0ANLehmjihgGL2BeZZBkBBuNxEpKaBRIXRc+TJrCHNg6xP6DY9OO66l1T7q+sj4/eaKJomSMc5vlHb02YozKKHkThbmf8mPGo0GaNXR1UUpxuREXn7GQ6+96cNqk6+xaF/tj5j6iqHWJUw+sGaG2aeGyjbjnko5VyCR3WsZ91D309tRYc/4pUybnA1xKh7iG86bJGIwQ68kbMZR3/qDVCXcz3f2Ul9TyJCLyXBH5rIjcJiIbg39lCDeTSVNKSd9fueJUrr1oybQqqcMRUUNJHD5nltNL1YxQ2xVLvUqv3TFLHjcSEpl2bBa5k+49LsR53uzodnUpHTKwdSi2pHXcfcXJGIwQwzTisslbFqfKKxka8biMPL4CfBr4LJBcm9kojDSFnPZ9VG/KJW8hjKuxyTLpmMU9EWwv2iedNqrJInfSvWc1qi4juLgsf/HvK+5+omTsC819tHJ98Krl+WSlVS63VuNiPMZU9VNNl8SYQppCzhMlctbJCyLdWXNqXZG5DK69+7jz1iuNPO6JRn3SSS/2FbfumHbfWeWOUpji75PVLeSiROP2UeLbMEmpFz1JDdl+q4GtQ3QV7D4rk1a73FqJi/G4VUTeBfwHoXU8VHVP06QyUntxWXt5A1uHuGnL0BQFL8CFp/fR/+yjc/fuk84bNfLJMzeSV8GlvdjB3NGaW3ZMRl/NqU335CbJHYTURi3TGqUQa93x80guHYKkyK04ypwUzvJbBb9PVDu1S8RT1UrrlImL8Xi7//+q0DYFTixeHCMg7YXPqhDikgc37dzNlStOzXQu1/PWU5Z7Irzmdj1RL/aBsUNBBHv3j07rOaZNUAcKc9najeluwYRQN5cOwarli6clLQI89fRYYpRYFSeF4+aEukXaptR7u7vcGiHVeKjqCWUIYkwn7YXPohDyzJEkkaSg464X12tWYNnajYX0htOKONbLltZzDCaoo3R+vdwuCmN0Qhsqf75iad+UkZLLeatKXHtNqLbNfczkGlepxkNEasCfAa/0N90OfEZVs4XuGC2lyIfcRUFHnTdporooX3FSEcco2dKMX1IZ+uD4tDIuUcdA/HxM2v0/EZPg2G693U5QvDM5yTA1VBf4FHA68E/+v9P9bUYbEbUiXN6HPE1Bx503CL/ti1EORYRnpinQsGwuYa8uCjmQ+6yTF8SeL0y3SO6VAMOyuW6vKlmfySqu3LdiaR8Xnt43GVLeLRI539eJuMx5vFhVTwt93igi25olkNEcipw0TVKoaSVJgp71CavXR/boG+k9J0XuRMmWNKrYf3CME1avTzxfmKHhEf69LuIsjnHVhiZay+7tuoaiZg1ZzfJMVjWqKQgYCZ6RcVVu2jJE/7OP7ngD4mI8xkXkJFX9OYCInIjle5ROEbHkRU2aJuUNuK78V7TL4oMD26eFCwf01LojJ2CTDFUQwutiOAJcc/e7JN5dNjQ8wgmr16dW04XyloZ1Udp5lXta2ZvgHqOMeBWimizaKplVwCYRuRcvCvPZwDuaKpUxhVb2uuLqZzXa8y2y9zywdSjWcHQJsZE7WYpIdoswocpRPTX2HRybFu2UhYmUQwM31qqvbOOKW3dMW18FyoueclWOzS6yGGfEWz3PM5OjrVLnPFT1O8BzgPf6/xar6qZmC2YcolXlG+L88sDk3EW4/EkWJRGe/wifA8js105yPyUp6iifexwTqty39rUMXn42V7/htPQDCmB0Qtm7fzTznEiRuCrHopWoS+ADtH6ep1Pmn/LgEm31RuCbqnq3iHwQuFxErlTV/2q+eAa0rneTlhyXtUcZNYoJu7nyjrDS2uGvb757ynXPOnkBm3bu5pHhEY7qqTGn1jXZu49aRwSmKoMVS/ucS710CagWUyq/CHdIVvdnsxcCi8Pl2a5CVJNFWyXzv1T1SRF5BbAc+CIWbVUqze7dhKNYllxxG0s/dBsnrF6fuXpr2jXSoovyjrDS2mH/6MSU637prgcnPw+PjPL06ATXXrSEO1e/ijXnnzJtNCJ4JVjCuI5aJnzDsWr54tgoMyC2AGQ9jQYUZI3wirvPfQfGphxXZDQfxP+m3SK5R7vNIG4E3Wq5ysDFeARv82uBT6nq14DDmieSUU/RL2aYeoUyPDI66SqJo6jlcesNQ94RVhb3UxQjo+OsuWUHy9ZuZOW6QepTAhW4acvQFGUZpTTeesbCSQMRNgWBkj7r5AXUuqcbiVqX8OaXHu90D8f19uQOWc1jnIP7rC9/MjwyOsXw1Idhd4tMnjuPqy3umb/mTadx39rX5hr5NosVS/u4c/WrKidXs3GZMB8Skc8ArwY+KiKzcTM6DSEivwP8b6Ab+GdVXdvsa1aVIqNr6t0We/cdiF3fPArXle7qZY0zAEPDI5NZ2nldH0mFDl0ZHhmddFdFtUd9OfykBbred+O2yMigTTt3c/UbTpsiZ3htj/5nHz157qiJ+Z5aN2edvCB38EQj63VfvWHXtLatd6NFVUHOG9xRZkSZkQ/RlFBEEZkL/A6wXVV/KiLHAqeq6m1NE0qkG/hv4DXAw8CPgDer6o+j9u/v79fNmzc3S5yOwSUzPA4Bpxc4KmQ2qXJveJ8LT+/jpi1D0/zHWdwASz90W24D4kJPrTtWvrT2FeC+ta91vtYHB7Zzww8eYlyVbvFGJ5t27s4dJh1Xe8vl2Li8nPp7auQaRvmIyBZV7c91bJLxEJEu4G5VfUFe4fIgIi8D1qjqcv/zZQCqelXU/kcccYSefvrpJUrYnmx9cJgDY/lSdGbP6ub4o3s45vDZsfs8/tQBfvbYU5HfzerqYkKViYTnLbjGQ3tGODA27nTNeu57fB+/+NXTzvtnQRA0QoXOntXN0oW9qe0rCCc9c96U+3n8qQOR9/v4Uwe4d/e+xPaq54wTn5H4fVzbPOvIOZxwzLzEY+PuLbj3gLvu/WVu+apG3G/TSdxxxx25jUei20pVJ0Rkm4gsVNUH84mXiz7godDnh4GXhncQkUuASwBmz+6sH7RZ5DUcwbH37t4HEPsC3f/4/tjjxyYm+I1nHj75MsZd45jDZzf0gsYtYBUo/tmzvFFQXH2oOLr8PI8ogvtJa19Fp7RhvYEIt/FDe0YyGY7Zs9LnS+LaxmXRr+OP7plmzLpEOP7oqS7F2bO6Y9th64PDhSjgMpR60m/TaQYkLy5zHscCO0Tkh8C+YKOqnt80qYgsETTlTVLV64DrwHNb3X777U0Up/qE5xl659ZQ9QrohV1NTiXD8cJL4/Ijjurt4fYY98Oi1etjzxl2WyS5NuLO7cLA1iEuXTcY+V3YveLSDszppZIAABs8SURBVLUu4fA5syZDeM86ecGkCylObtf2PSq0/zMj9j+qt4d9fgCDC10Cf/+mJamuvSTX0+0O7jSXMN+BrUOs+so2RmMeoPFaN5c2EI0UuAbnh1yDjZ4ziqTfppFntGqIY5RfFC7G44rcZ8/Pw8Dxoc+/DjzSAjkaoqzlKet97WGff9qqd1G85aULYzO2HxkeyXVf9WtSFB0bH7RBHOFJ96QJ4qi5HddFi1zb12Up4Sx5IUmrCIbJGpCQlpcTS4I+ajRXpaxyIDM5c9wVl/U87ihDkDp+BDxHRE4AhoDfA97SAjlyU2ZJkbRs3KhV7+LqBYG3kFOcojmqpzbtvi5dN8iaW3Yw77Bu9h2cLse8w7qnrUkRlqMIw5rWBvsPHlosKWttrqRFiy483YtEWrlukON6e7jw9L7YSe0Al6WEVy1fzMp1g06jjyjvVp6yMvWj16eeHpscQbg+v1dv2JVauqURBZxUE6xIOqFcfLNJDbkVkTNE5Eci8pSIHBSRcRH5VTOFUtUx4D3ABuAnwI2quqOZ1yyaMkuKuLyM4UWfgpj0OJ/6I8MjsXH2IkQq0uGRUQ6OTdDdNbXbWesWPvL6U6ftX3RsfFobBCsEDmwdypw3E3fuoIJqOOnupi3e+Xt7opeFFZgyUomTY8XSPi4+Y6FzifcwecrK1B+zd//oNNeTy/Pr8iw2ooDjkildkyxdaWZuVafg4rb6JF7P/ytAP/A2vFpXTUVVvwF8o9nXaRZlDntdXBxRL2zv3FpkWOtxvT2xo4OVMXMK4NVi6ql1cVCZDC+96MXHOxmGqJ5y1PXjzuXSBiOj47zvxm1c86bTuOqCUxs+d5AIV3+NNbfsYN/BschzXXzGQuelhK9cceqU3I+5MSO7N7/0kIc3Kc8kqayMay2ptOc37XdoVAHHFUjMUv3YBcszScfFeKCqPxORblUdB/5FRL7fZLnanjKHvWm+9qjyGgNbh3jq6ekKrtYtky93VOXWtJpO4QS7pLUNwsaiPiFuaHiEVV/dBoqz28TVzTOuymU3b+eqC051zjuIc/fEtXdUbSyA+XNrk+vFB2RZanhg6xB/uW5wWun3L931IJt27uaskxdMWVuiniTF79qpSXt+o9oqWMY3ba0XF/oSXI5FU1bl4nbFJVN8v4gcBgyKyMdEZCWQHBRulDrsrS+VMbc29WeNKq9x9YZdkREx8w6blfjCZC0FEuXqiCqJUu8nHx3XTG6TFUv7nKOTsroP4+oXZVVYLiGxSVy9YVfsmiFDwyNcf9eDmZcGdvkuwOX5jWqray9awv0FuSfNnVQdXEYev49nZN4DrMSLgrqwmUJ1AmUv2BO+zv6DY+yvK7ExMjrOFbfumNwvTtEG+Q9p62v/9c13T7tGHPW9WlcXicu5wsT1SrOeB9wjjaJGJHHZ9C4KOimSLU3mJOOZpmBXLV/Mqq9um2LEu7uEI2bPmhbynUYze+zmTqoOLtFWD4hID3CsqrYibLdtaeZLFCiZoeGRSbcAJEed7N0/mlq6Iyi8lxQpFtQ62p/T1dHIvE+SAnYNlU07j2uk3IqlfWx+YM+UEiIXnu7VqMoTipx23bxl3btF3Eq81FmfLpisu1UlzJ1UDVyirc4DBoFv+p+XiMgtzRbMiCfs9oHkHmcWAgXXSAXcuHOGaWTep37uBg6VlF+5bpDZs7qYP7c2pdJtVjeHa6Rc3PrVkG+xrLTr5qkeLMA1bzrN6dr1bsLRCW36gmNG++LitloDvAS4HUBVB0VkUdMkmsG4Jt814vaJoj4xLi6iKmwwXHrBUROkA1uH2B8RiRTO6E4yhpt27p7yub63PjwySk+tm2svOpRxHY5YcnFzJFUAXnLFbZNunH0HxiKV/aXrBnNNDqdF6IVdNvUjzjhcEwibHR1YVsKsUR4uxmNMVZ9oJI3dSCfKZbFy3SCbH9gzLULH5YXu7akxb/YshoZHEsuNRCXGuUSKuUR41Z83rupsuCw5JJcPCe497LarJ6pUeBZFlWQYg0iqNMOZJyk0S4Se4CVsingT8XEJn3H5Jo1cOytlJswa5eESbXWPiLwF6BaR54jIJwAL1S2YqNGEAtff9eC0xXTSXuieWjdrzj+FVcsXU+uWWMMR575xiWgJomrikrOiZIwbMc2bPTXCa9XyxbHJceE5mSQFHjawWRdPanRxqYCsUV1R1w2HWUdFqQWrIF7zptOodU1vtX0Hx1LvN+7aRUUxlZkwa5SHi/H4c+AU4ADw78ATwKXNFGomEjeaUJj2ksUpGZjqX08qFZE0iZq2tGZ4juHInlnTlFaU0hnYOuS8rG1cdnXSnEw9gfHKs/RqcP9FEB4ppRmwFUv7uPD0vin3HQ6zTqvrdPic6Y6E0XG3eYtmLqdqdaI6k1i3lYjMAf4U+A1gO/Ayv2yIUTADW4di3Q4w9SULLxAUEOdfT3o5J1Sdk9PqZa0vwljrFnp7arEhnVmKFgbUZ1e7zMkEhI1X3kJ6gfFttGaSS/RamE07d0+bxwjkTZqLGdg6FJtH4qqkmxXFZHWiOpOkOY8vAqPA/wPOAZ6HjTgKJ6lia0Dwkn1wYDtfumv6sipnnbwg8qVP8t3nfXGjlPHouDJv9iwGLz/b+ZiAJNdInDJLuq9efx5g5brBROUfpVDrJ3WDjO28wQm1LkmNXnM1+GmVdpOy64tS0nknvZtRRdloPUluq+er6ltV9TPAG4BXliTTjCLNBRN+yW74wUOR+8RtD+Y86gmUWh7yuCCSvsvjGonzz7/1jIUcGJtgrx+xFUQkRVGvUKPcWzdtGeLC0/smXTnz59boqbl4ej0On+PN5WRpszhFHyjruLmYOMNRXzU3y9xPmDzuv4BmusSM1pE08pgcA6vqmEVbNYe0jOlw7y5rUbjguCtu3TGZHFgf2ZSVPC6IpBLoeeSoD1kNChRGLdakMC2kNarXGzc62LRz92TUWJr7rZ69+0cZ2DqU2mb1db5q3TJlripcaReIXfCqnvDz02jEU6PraFhiX+eRZDxOC5VeF6DH/yyAquqRTZduBhCnWIIQy8D9smr5Yrpj5kWSylEX/dLmcUHEFcuLSviLIsldEj5vnBENivIluVtcRgd58msuu3k7F57eN839VesS9h8cY9Hq9VOM2/DIKLUuYf7c2uQqhmF5Xedi6kOlG1X+Nult1BNrPFS18VhFI5UoxVrrEvYdHJuSU3DZzds548T53PnzPdPOES7J3Wzy1BYKyniEVycMoojCFXfjyrLH9ZhdlXncIk9hXEZUSYqyfrQQEIxewiXggyrCwWiw/qjRCWXuYbPY+rdT55DiStLE3U+YRpW/TXob9TiVZDeaR5Qy3h9SLAEjo+Pc/8sR3nrGwim1lN780uOnJRG6kmUCtNEM4aQooji3ysp1g8ypdU0p8x4+zkXxuU7MnnXygshghD37DjitQLhq+eJYd9IjwyNTRoDL1m6MLdsePiZMffuE3XEubrlGlX+rJ70tQ716mPGoAPWupRNWr4/cb2h4hE07dzvVKkojiw+8iAzhpDDTE1avjwxVVphmOMLnO6qnFqmEu0WYUM2kZOrLngSMjE4krgEfno+IcyflKQxZf0xcEun8uTUuP++UVMUa5zocGh5h2dqNTqPHQI6yFbhlqFcTMx4VJCkks6gXJ4sPvFF/OSTfk5J9JbjA9VNPrUu4+o3ZjWuSQo9bA75egbr2zvOsthcnXzBCTXPLJdXFcn2mWjXpXcTzZxSPe9yhURpp5TGKKO2QxQdexGRpUSU/4NBa6lFzDIfN6uLqDbsyh6OmuW+i1oCvX9woyBAPAhiCEu1Ro4CoEGrwRhJRYaxJ8l26bnDyXpPCcQPZ+3p7Yl2IVcQm66uJGY8KEo6Lj6PRFycpn6CRfeOoV6x5CcqqxGVT7zs4nisXIc24uS7kFFWivf76K5b2Me+w6EH/3JiVHNPmFoaGR1j1lW2s+uq21PvPqowbyQ8pgiKeP6N4zHhUlHAvMYpGX5wshfCKKJpXr1jzEpRVcb1/1x51YNwiagvSU+vmrJMXpCrQLAUAn4iZMG+kUzA6odNGY1HXz6KMG0kOLApberaamPGoOM16cbJk/RaRIewSVhtkiSeNuHrnevkvWdxgLgo5MG71FYh7e2qTeRpF9uhdFHi4x/++G7el3kMc9dfP8kxVoSKuZahXE5swrzjNjHLJMgHa6GRpkgKvX4wKYMkVt0VGUgUDl6h22XdgLPIYl1FKUrn4TTt3O03YZgmHTZtcr48wamTEVn/9LM9UVeYbLEO9epjxaAM64cVJypGIihSKc+uEt9e3y8DWIVZ9ZduU5VRd63gVUbMrziAELq8oRR2nwPNks9e6BOoCCeJGFK7PlCUHGnGY8TBKYdXyxZkUe26lVT9n4Tg/n3Y919pUvXNrzJ7VNVmeftEzeqZk1teHxcYpcJeefa1bmHfYrCml8KHYUWqrkwON6mLGow2oWnZtbnkyKPY8Sitq8atgMaQ0+dKuF/dd1PomwRrqEF0q3SVHIc6YuSRAFvls5HWb1hd7DJbLrcLzaxSDGY8KkbW2UytewLzZvlkVex6l1Yh/Pu56wbaR0fHJwpTharXL1m5MnFCOm6lIkynOmLViojhuhBTXiah/RsLzUK1+fo3iMONREeKU8uxZXZXKrs2b7ZtHsWed62nUPx81h1I/aV1fHj2vwQqixsLXqlfE4WKKVeuxJ3Ui0uZrLDu8M7BQ3YoQp5TjCui1Krs2r7IsI9Gr6LBmlzDVpPtKurdw8FRcLgUQm83eapLaxuXZtOzw9seMR0XI+jK1KtrlqJ5apu0BZSR6FZ0P4GIok+4r6d6GR0an5HC0OpciK2nL5abRCdFarc68bzXmtqoIcS6X+XNrPD06UZlol7jqImlVR4I1PcLl5KPqPjVKkWHNLm6wtLmZNbfsiBw9BhVtIT6Ho8q986S2iZqvCdMJ0VpW6beCIw8RWSMiQyIy6P87t9UylUFcD/by806pVHZtXE2puO0BrnWfmk2W3qLraCmqWGJwnTi3o0vKX5V750ltUz8C7O2pMX9urRLPb1FUIfO+1VR15HGtqn681UKUSVoPtiovW95J6SqU1c7aW2wkTDWp5+1Cnt55mYt7uTyvVXlmm0FVMu9bSVWNx4ykHV64vEljSYtBxWVfF00eA5bnN8mTHQ75FrEKKHtxr2Dfqj+vzcIy7yvotvJ5j4jcLSKfF5H5UTuIyCUisllENu/eHb0KnFE8eSel416qwPdfRsXWRnuLri4vl/PVTxH11Lq55k2n5Y6syuJGMZdL41il3xaNPETk28CvRXz1N8CngA/juYU/DFwD/GH9jqp6HXAdQH9/f2N1vo1M5Olxxi2Dmif7Oi+N9Baz9NZdVgq88PQ+Nu3cXdiIq+zFvWY6zSxY2i60xHio6qtd9hORzwJfb7I4RsEk+dPD2+MUbLOUWCN1mrK4vJIMZV+TlEwWw2gul2KYyW47qOCch4gcq6qP+h9fD9zTSnmMbKT10MMv27K1G0tVYo30FrP01lvRK81S0deKHRpFINrgym5FIyL/BizB66jdD/xJyJhE0t/fr5s3by5BOiONOIMQVXo9Kiqp2T30rASjqLhRUlxJ+ahzNNuQfHBg+5Q8mjNOnM9/PfhEZH0smNkuF8NDRLaoan+eYys38lDV32+1DEZ+8vbQh4ZHpsyBFJF01ajSTgu5demtl5VMFpVH8/2f74mdU6pauROj/ahqtJXRpmStYRUk2PX19sQqujx8cGA7K9cNNhTJlRRy6xplVlZkU9R14nwKQ8MjM66UhlE8ZjyMQskbwljk5PnA1qEpCzAFZFXacdcWcO65lxXZlPV8zQyJNmYGZjyMQsmTBzKwdSh2Xag8k+dXb9iVex0Nl2tnkamMasJJ54trV8vrMBrFjIdROFG1npK44tYdkcpeIFcEUJKByKK0i0gEKyuZLO46F5+xMPYYy+swGqFyE+ZG+5O1xtLemKKKSr5J5bg8hqzGqIiQ27LCdpOus2nn7oZDoqu2FLLReioXqpsHC9WtDlERSknLp8aF9oJbGKyrDAJcfMZCrlxxaubzFU3Zijjrb1L08UZ1aSRU19xWRqFkjS5Kcp3kde1Ezbtce9GSyhiOqFUDmzl53egiWVYLy4jC3FZGoWSNLopzMfX21Brq1Yaz2YOe/sp1gy13ubSqNH0jpTSsFpYRhY08jELJGl0UN9G75vxTCpGnFT39JNpREZcVMWa0F2Y8jELJGl1U9Lrj9VTN5dKOitjKjxtRmNvKKJQ80UXNrE5atZ5+OxYltPLjRhRmPIzCqVKp6izlx8uIgmpXRVyl39SoBmY8jI7GtadfVgHD4HymiI12x+Y8jI7GdU6lanMjhlF1bORhdDwuPf2qzY0YRtWxkYdh0J5RUIbRSsx4GAadGY46sHWIZWs3csLq9Sxbu9FKsBuFYm4rw6B9o6DiKDMAwJiZmPEwDJ9OioJqVRkUY+ZgbivD6EAsAMBoNmY8DKMDsQAAo9mY8TCMDqQTAwCMamFzHsaMplNXyGtlAECntqkxFTMexoyl0yOSWhEA0OltahzC3FZGIp2cK2AlSYrH2nTmYCMPI5ZO70VaRFLxWJvOHGzkYcTS6b1Ii0gqHmvTmYMZDyOWTu9FWkRS8VibzhzMbWXEkmUhpXak00qSVAFr05mDqGqrZWiY/v5+3bx5c6vF6Djq5zzA60UWuca4YRitQ0S2qGp/nmNt5GHEYr1IwzDiaInxEJE3AmuA5wEvUdXNoe8uA94JjAPvVdUNrZDR8OikYoGGYRRHq0Ye9wAXAJ8JbxSR5wO/B5wCHAd8W0Seq6rj009hGIZhtIqWRFup6k9UNSre83XAl1X1gKreB/wMeEm50hmGYRhpVC1Utw94KPT5YX/bNETkEhHZLCKbd+/eXYpwhmEYhkfT3FYi8m3g1yK++htV/VrcYRHbIsPBVPU64Drwoq1yCWkYhmHkomnGQ1VfneOwh4HjQ59/HXikGIkMwzCMoqhaqO4twL+LyN/jTZg/B/hh2kFbtmx5SkQ6o2ZG4xwDPN5qISqCtcUhrC0OYW1xiNyp/60K1X098AlgAbBeRAZVdbmq7hCRG4EfA2PAux0jrXblTXTpNERks7WFh7XFIawtDmFtcQgRyZ1d3RLjoar/AfxHzHcfAT5SrkSGYRhGFqoWbWUYhmG0AZ1iPK5rtQAVwtriENYWh7C2OIS1xSFyt0VHFEY0DMMwyqVTRh6GYRhGiZjxMAzDMDLTVsZDRH5HRHaJyM9EZHXE9yIi/+h/f7eIvKgVcpaBQ1tc7LfB3SLyfRE5rRVylkFaW4T2e7GIjIvIG8qUr0xc2kJEzhSRQRHZISJ3lC1jWTi8I0eJyK0iss1vi3e0Qs5mIyKfF5HHROSemO/z6U1VbYt/QDfwc+BE4DBgG/D8un3OBf4vXpmTM4AftFruFrbFy4H5/t/nzOS2CO23EfgG8IZWy93C56IXL49qof/5ma2Wu4Vt8dfAR/2/FwB7gMNaLXsT2uKVwIuAe2K+z6U322nk8RLgZ6p6r6oeBL6MV4U3zOuAf1WPu4BeETm2bEFLILUtVPX7qrrX/3gXXqmXTsTluQD4c+Am4LEyhSsZl7Z4C3Czqj4IoKqd2h4ubaHAESIiwOF4xmOsXDGbj6p+F+/e4silN9vJeLhU3HWuytvmZL3Pd+L1LDqR1LYQkT7g9cCnS5SrFbg8F88F5ovI7SKyRUTeVpp05eLSFp/EW5DuEWA78BeqOlGOeJUil96sWm2rJFwq7jpX5W1znO9TRM7CMx6vaKpErcOlLf4B+ICqjnudzI7FpS1mAacDvw30AP8pInep6n83W7iScWmL5cAg8CrgJOBbIvL/VPVXzRauYuTSm+1kPFwq7s6UqrxO9ykiLwT+GThHVX9Zkmxl49IW/cCXfcNxDHCuiIyp6kA5IpaG6zvyuKruA/aJyHeB04BOMx4ubfEOYK16jv+fich9wMk4FGPtMHLpzXZyW/0IeI6InCAih+EtV3tL3T63AG/zowfOAJ5Q1UfLFrQEUttCRBYCNwO/34G9yjCpbaGqJ6jqIlVdBHwVeFcHGg5we0e+BvymiMwSkbnAS4GflCxnGbi0xYN4IzBE5Fl4FWbvLVXKapBLb7bNyENVx0TkPcAGvEiKz6tXhfdP/e8/jRdJcy7e8rX78XoWHYdjW/wt8Azgn/we95h2YCVRx7aYEbi0har+RES+CdwNTAD/rKqRIZztjONz8WHgCyKyHc918wFV7bhS7SJyA3AmcIyIPAxcDtSgMb1p5UkMwzCMzLST28owDMOoCGY8DMMwjMyY8TAMwzAyY8bDMAzDyIwZD8MwDCMzZjyMGYFfTXcw9G+RiHzf/26RiLwltO8SETk3xzVuF5GGw6GLOo9hNBMzHsZMYURVl4T+3a+qL/e/W4RXMDBgCV7cu2EYMZjxMGYsIvKU/+davKzrQRH5APAh4CL/80UiMs9fE+FHIrJVRF7nH98jIl/210BYh1crqv4a54jIjaHPZ4rIrf7fnxKRzf5aElekyIiIvEFEvuD/vUBEbvJl+pGILPO3/1ZodLVVRI4ooq0Mo562yTA3jAbpEZFB/+/7VPX1oe9WA+9X1d8FEJFfAP2q+h7/898BG1X1D0WkF/ihiHwb+BNgv6q+0K8j9l8R1/0W8BkRmefXk7oIWOd/9zequkdEuoHviMgLVfVux/v538C1qvo9vxTNBrwKse8H3q2qd4rI4cDTjuczjEyY8TBmCiOquiTnsWcD54vI+/3Pc4CFeIvs/COAqt4tItMUv18m45vAeSLyVeC1wF/5X79JRC7Bew+PBZ6PVzbEhVcDzw9VCT7SH2XcCfy9iFyPt27Hw9lu1TDcMONhGOkIcKGq7pqy0VPcLvV91gHvxluQ50eq+qSInIA3Snixqu713VFzIo4Nnz/8fRfwMlUdqdt/rYisx5uzuUtEXq2qOx1kNIxM2JyHYcCTwBEJnzcAfy6+tRCRpf727wIX+9teALww5vy34y0D+sccclkdCewDnvArup4Tc+wvROR5ItKFt6BVwG3Ae4IPIrLE//8kVd2uqh8FNuOVGDeMwjHjYRieq2hMRLaJyEpgE55LaFBELsKrvloD7haRe/zPAJ8CDvfdVX9FzDoQqjoOfB3PQHzd37YN2ArsAD6P526KYrV/zEYgXCb7vUC/P1n/Y+BP/e2Xisg9IrINGKFzV5A0WoxV1TUMwzAyYyMPwzAMIzNmPAzDMIzMmPEwDMMwMmPGwzAMw8iMGQ/DMAwjM2Y8DMMwjMyY8TAMwzAy8/8Bu1gG89x8lXoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(yhat, res.resid_pearson)\n",
    "ax.hlines(0, 0, 1)\n",
    "ax.set_xlim(0, 1)\n",
    "ax.set_title('Residual Dependence Plot')\n",
    "ax.set_ylabel('Pearson Residuals')\n",
    "ax.set_xlabel('Fitted values')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of standardized deviance residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAYdklEQVR4nO3de5hcdX3H8ffHEOUulywYCbiKYAVawuMasaBFLjUCGrBSwUqTioZWUbBUDd4K9RYfFWirYoPwEOSiqCAIXojBNKIIbjRgYlBUIreQLJeYRCoa8u0f5zdwmMzsnN2d2dlf9vN6nnn23M/3nD3zmd+cc2ZGEYGZmeXnGd0uwMzMhscBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWKQd4Imm5pMO6XUc3STpe0r2SNkg6qNv11Eg6W9JlbVzeJZI+mrpfIemX7Vp2aR0h6YUVp23L9kn6gqQPjXQ5Y9Fgz09Jh0m6r03rWSTpre1Y1mgYFwEuaaWkI+uGzZJ0c60/IvaPiEUtltObnphbdajUbvs0cFpEbB8RP6s6UzkQcxMRP4iIF3W7jnaIiH+OiI90u45OqPL8HI/GRYDnYgy8MDwPWN7lGtpK0oRu12CFMXB8b3Ec4Em5lS5pmqR+SeskrZZ0bppscfq7Np1meLmkZ0j6oKTfSVoj6VJJzy4t9x/TuIclfahuPWdL+pqkyyStA2aldd8iaa2kVZI+K+mZpeWFpLdLukvSekkfkbR3mmedpKvK09dtY8NaJT1L0gZgAnC7pN80mFeSzkvz/V7SHZIOkDQb+AfgvWmffDNNP0fSb1KNv5B0fGlZsyTdLOnTkh6VdLek15TGP1/S/6Z5FwCT6mr5qqQHUx2LJe1fGneJpAskfUvSH4BXSTpI0k/T8r4CbF2a/sm335LemLah9nhc0qI07lmp3nvSMfEFSduUlvOe9P96QNJbGu3/IWzfwZJ+lI6B25VOHUg6UVJ/3bTvlnRdadtrp4Z2lnS9pIG0j6+XNKU036J07Pww1XGjpEml8YeWarhX0qwq+6Gutllp+edJegQ4e7D5JU1Kda6V9IikH0h6RhpXft5sk7b1UUm/AF5at96nnb4ayn6pW84L0//p95IeSsfO2BIRW/wDWAkcWTdsFnBzo2mAW4CTU/f2wMGpuxcIYKvSfG8Bfg28IE17NfClNG4/YANwKPBMilMUfy6t5+zUfxzFi+k2wEuAg4Gt0vpWAGeU1hfAdcCOwP7A48DCtP5nA78AZjbZD01rLS37hU3mfTWwBNgJEPBiYHIadwnw0brpTwCem7brjcAfStPPStv9NooXjX8BHgBU2v/nAs8CXgmsBy6r244d0vjzgaWlcZcAvwcOSeveEfgd8G5gIvCGtO6PpukPA+5rsL07pn1/auo/P+33XdK6vwl8Io2bDqwGDgC2A65osS+bbh+wB/AwcHSq/6jU3wNsm6bdp7SsnwAn1v8fgF2Bv0vz7AB8FfhGab5FwG+AfSmOu0XA3DRur7Sek9I+2xWY2mo/NNjOWcBG4J0Ux/M2LfbjJ4AvpHVOBF5ROiZW8tTzZi7wg7SMPYFl5f9h/b4fxn55a+q+EvhA+j9sDRza7SzbbB93u4BR2cjin78BWFt6PEbzAF8MnANMqltOL5sH+ELg7aX+F1EExFbAh4ErS+O2Bf7E0wN8cYvazwCuqTs4Dyn1LwHeV+r/DHB+k2U1rbW07GahczjwK4oXl2fUjXvyCTLIdiwFZqTuWcCv6/ZLAM+hCI+NwHal8VdQCvC65e6U5n12qZZLS+NfSenFIQ37EYMEeHrCXg9ckPpF8QK0d2malwN3p+6LSeGX+vdtti9bbR/wPkovqmnYd0kvysBlwIdT9z4UQbttq/8DMBV4tNS/CPhgqf/twHdS91nlY640zaD7ocH0s4B7qs4P/AdwbZP9tpKnnje/BaaXxs2mYoBX3C+1AL8UmAdMGezY7uZjPJ1COS4idqo9KA7YZk6heBLeKeknko4dZNrnUrTwan5HEd67p3H31kZExGMUramye8s9kvZNb+seVHFa5ePUvcWmaO3V/F+D/u2HUeugIuIm4LPA54DVkuZJ2rHZ9CpOHS1Nb4fXUrROy9vxYGnZj6XO7VONj0bEH+rqrC13gqS5Kk7PrKN4YlO37PI+fS5wf6RnZP3ymvgYRevsXam/1vpdUtqe76ThtXWU1znY8gfdPorrECfU1pPWdSgwOY2/gqJlDPAmitbjY9SRtK2k/1FxumwdRaNkJz39msCDpe7HeOq42ZOidV6v1X5opLxfWs3/KYp3iDdK+q2kOU2WOZT9/TQV90vNeyledG5TcRfMoKfGumE8BXhlEXFXRJwE7AZ8EviapO0oXtnrPUDxpKuptbBWA6uA8nnHbSjewj1tdXX9FwB3UrxN3hF4P8VB1A6D1dpSRPxXRLyE4tTNvsB7aqPK00l6HnAhcBqwa3rBXEa17VgF7Jz2d7nOmjcBM4AjKU4Z9dZWWy61bnl7SCqPLy/vaSSdSBGQb4iIP6fBD1G8MO5fagQ8OyJqgbeKIvRaLr/C9t1L0QLfqfTYLiLmpvE3ApMkTU11XtFkPWdSvMN6WTqOXlnbxEFqK9ewd4PhrfZDI+X/xaDzR8T6iDgzIl4AvBb4V0lHNFhmq/39GMULRc1zSt2V90tEPBgRb4uI5wKnAp9XxVtDR4sDvAFJb5bUExGbKE63ADwBDACbKM4h11wJvFvFhantKVrMX4mIjcDXgNdK+msVFxbPofUTaAdgHbBB0l9QnB9ul8FqHZSkl0p6maSJFG+D/0ixT6B4ASjvk9qL3UCa958oWuAtRcTvgH7gHEnPlHQoxZO5ZgeK8/4PUzxJP95ikbdQvEi9S9JWkl4PTGuyjQcB/03xbm2gVNMmihek8yTtlqbdQ9Kr0yRXUVyA3k/StsC/j2D7LqM4Zl6d3m1sreJC65Q0f+24+hTFOeAFTVa1A0VYrpW0y2A1NXA5cKSkv0/7bFdJUyvsh0G1ml/SsenCoSieA0/w1DFWdhVwVrogOYXiHHvZUuBNaf9NB/6mNK7yfpF0QukC56MUx3SjerrGAd7YdGC5ijsz/pPiItEf01vVjwE/TG8BD6Y4//klirdid1ME2zsBImJ56v4yRathPbCGIoCa+TeKVuZ6ioO9nVe+m9ZawY6pnkcp3rI+THFRFuAiYL+0T74REb+gOBd/C0W4/yXwwyHU+SbgZcAjFE+wS0vjLk3rv5/igu2PB1tQRPwJeD3F+dhHKS6oXt1k8hnAzsDNeupOlG+nce+jeHv/4/TW+3sULTki4tsUF+duStPcNNzti4h7Ux3vp3gBvJfinU75uXoFxTuQrw7y4ns+xUXDhyj20Xda1PSkiLiH4iLqmanGpcCBaXTT/VDRYPPvk/o3UBw7n4/G936fQ3EM3E3xjuRLdeNPp3hRXEtxh9Q3SuOGsl9eCtyacuA64PSIuLvSVo6S2hVeGwWp1buW4vTImDoQzCw/boF3mKTXpgsn21G0WH/OUxfezMyGzQHeeTMoLh4+QPEW8cTw2x4zawOfQjEzy5Rb4GZmmRrVL5eZNGlS9Pb2juYqzcyyt2TJkociYrMPTI1qgPf29tLf3996QjMze5Kkhp829SkUM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMjeonMW386Z1zw5CmXzn3mA5VYrblcQvczCxTlQM8/b7czyRdn/p3kbRA0l3p786dK9PMzOoNpQV+OrCi1D8HWBgR+wALU7+ZmY2SSgGefpn5GOCLpcEzgPmpez5wXHtLMzOzwVRtgZ8PvBfYVBq2e0SsAkh/d2s0o6TZkvol9Q8MDIyoWDMze0rLAJd0LLAmIpYMZwURMS8i+iKir6dns+8jNzOzYapyG+EhwOskHQ1sDewo6TJgtaTJEbFK0mRgTScLNTOzp2vZAo+IsyJiSkT0AicCN0XEm4HrgJlpspnAtR2r0szMNjOS+8DnAkdJugs4KvWbmdkoGdInMSNiEbAodT8MHNH+kszMrAp/EtPMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8tUlR813lrSbZJul7Rc0jlp+NmS7pe0ND2O7ny5ZmZWU+UXeR4HDo+IDZImAjdL+nYad15EfLpz5ZmZWTMtAzwiAtiQeiemR3SyKDMza63SOXBJEyQtBdYACyLi1jTqNEl3SLpY0s5N5p0tqV9S/8DAQJvKNjOzSgEeEU9ExFRgCjBN0gHABcDewFRgFfCZJvPOi4i+iOjr6elpU9lmZjaku1AiYi3Fr9JPj4jVKdg3ARcC0zpQn5mZNVHlLpQeSTul7m2AI4E7JU0uTXY8sKwzJZqZWSNV7kKZDMyXNIEi8K+KiOslfUnSVIoLmiuBUztXppmZ1atyF8odwEENhp/ckYrMzKwSfxLTzCxTDnAzs0w5wM3MMuUANzPLVJW7UMye1Dvnhm6XYGaJW+BmZplygJuZZcoBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOBmZpnyJzHHMX+q0ixvboGbmWWqyk+qbS3pNkm3S1ou6Zw0fBdJCyTdlf42/FV6MzPrjCot8MeBwyPiQIpfoJ8u6WBgDrAwIvYBFqZ+MzMbJS0DPAobUu/E9AhgBjA/DZ8PHNeRCs3MrKFK58AlTZC0FFgDLIiIW4HdI2IVQPq7W5N5Z0vql9Q/MDDQrrrNzMa9SgEeEU9ExFRgCjBN0gFVVxAR8yKiLyL6enp6hlunmZnVGdJdKBGxFlgETAdWS5oMkP6uaXt1ZmbWVJW7UHok7ZS6twGOBO4ErgNmpslmAtd2qkgzM9tclQ/yTAbmS5pAEfhXRcT1km4BrpJ0CnAPcEIH6zQzszotAzwi7gAOajD8YeCIThRlZmat+ZOYZmaZcoCbmWXKAW5mlikHuJlZpvx1sjamDPUrblfOPaZDlZiNfW6Bm5llygFuZpYpB7iZWaYc4GZmmfJFzC2If+PSbHxxC9zMLFMOcDOzTDnAzcwy5QA3M8uUL2Ja1vzJTRvP3AI3M8tUlZ9U21PS9yWtkLRc0ulp+NmS7pe0ND2O7ny5ZmZWU+UUykbgzIj4qaQdgCWSFqRx50XEpztXnpmZNVPlJ9VWAatS93pJK4A9Ol2YmZkNbkjnwCX1Uvw+5q1p0GmS7pB0saSd21ybmZkNonKAS9oe+DpwRkSsAy4A9gamUrTQP9NkvtmS+iX1DwwMtKFkMzODigEuaSJFeF8eEVcDRMTqiHgiIjYBFwLTGs0bEfMioi8i+np6etpVt5nZuFflLhQBFwErIuLc0vDJpcmOB5a1vzwzM2umyl0ohwAnAz+XtDQNez9wkqSpQAArgVM7UqGZmTVU5S6UmwE1GPWt9pdjZmZV+ZOYZmaZcoCbmWXKAW5mlikHuJlZphzgZmaZ8veB27ji7w+3LYlb4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpapKr+Juaek70taIWm5pNPT8F0kLZB0V/q7c+fLNTOzmiot8I3AmRHxYuBg4B2S9gPmAAsjYh9gYeo3M7NR0jLAI2JVRPw0da8HVgB7ADOA+Wmy+cBxnSrSzMw2N6Svk5XUCxwE3ArsHhGroAh5Sbs1mWc2MBtgr732Gkmt485Qv/rUzMaXyhcxJW0PfB04IyLWVZ0vIuZFRF9E9PX09AynRjMza6BSgEuaSBHel0fE1WnwakmT0/jJwJrOlGhmZo1UuQtFwEXAiog4tzTqOmBm6p4JXNv+8szMrJkq58APAU4Gfi5paRr2fmAucJWkU4B7gBM6U6KZmTXSMsAj4mZATUYf0d5yzMysKn8S08wsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsU0P6Olmz8WY4X+m7cu4xHajEbHNugZuZZcoBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWqSo/qXaxpDWSlpWGnS3pfklL0+PozpZpZmb1qrTALwGmNxh+XkRMTY9vtbcsMzNrpWWAR8Ri4JFRqMXMzIZgJOfAT5N0RzrFsnOziSTNltQvqX9gYGAEqzMzs7LhBvgFwN7AVGAV8JlmE0bEvIjoi4i+np6eYa7OzMzqDSvAI2J1RDwREZuAC4Fp7S3LzMxaGVaAS5pc6j0eWNZsWjMz64yW30Yo6UrgMGCSpPuAfwcOkzQVCGAlcGoHazQzswZaBnhEnNRg8EUdqMXMzIbAn8Q0M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMtv07Wmuudc0O3S7AxaKjHxcq5x3SoEtvSuQVuZpaplgGefnV+jaRlpWG7SFog6a70t+mv0puZWWdUaYFfAkyvGzYHWBgR+wALU7+ZmY2ilgEeEYuBR+oGzwDmp+75wHFtrsvMzFoY7jnw3SNiFUD6u1uzCSXNltQvqX9gYGCYqzMzs3odv4gZEfMioi8i+np6ejq9OjOzcWO4Ab5a0mSA9HdN+0oyM7Mqhhvg1wEzU/dM4Nr2lGNmZlVVuY3wSuAW4EWS7pN0CjAXOErSXcBRqd/MzEZRy09iRsRJTUYd0eZazMxsCPxJTDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMtf9BhMJJWAuuBJ4CNEdHXjqLMzKy1EQV48qqIeKgNyzEzsyHwKRQzs0yNNMADuFHSEkmz21GQmZlVM9JTKIdExAOSdgMWSLozIhaXJ0jBPhtgr732GuHqzMysZkQt8Ih4IP1dA1wDTGswzbyI6IuIvp6enpGszszMSoYd4JK2k7RDrRv4W2BZuwozM7PBjeQUyu7ANZJqy7kiIr7TlqrMzKylYQd4RPwWOLCNtZiZ2RC04z7wLUbvnBu6XYKNQ0M97lbOPaZDlVhufB+4mVmmHOBmZplygJuZZcoBbmaWqWwuYvpCj9nw+Lmz5XIL3MwsUw5wM7NMOcDNzDLlADczy1Q2FzGHyp+qtC3VWDu2h1OPL5S2h1vgZmaZcoCbmWXKAW5mlikHuJlZprbYi5hmNjyjcZG0058OHWvLH846qnAL3MwsUyMKcEnTJf1S0q8lzWlXUWZm1tpIftR4AvA54DXAfsBJkvZrV2FmZja4kbTApwG/jojfRsSfgC8DM9pTlpmZtTKSi5h7APeW+u8DXlY/kaTZwOzUu0HSL0ewzk6aBDzU7SKGyDWPjtxqzq1eaFGzPtnZlQ9z+UPazyPchuc1GjiSAFeDYbHZgIh5wLwRrGdUSOqPiL5u1zEUrnl05FZzbvWCax6ukZxCuQ/Ys9Q/BXhgZOWYmVlVIwnwnwD7SHq+pGcCJwLXtacsMzNrZdinUCJio6TTgO8CE4CLI2J52yobfWP+NE8Drnl05FZzbvWCax4WRWx22trMzDLgT2KamWXKAW5mlikHeImkT0m6U9Idkq6RtFO3a2pF0gmSlkvaJGnM3oaV29cuSLpY0hpJy7pdS1WS9pT0fUkr0jFxerdrGoykrSXdJun2VO853a6pKkkTJP1M0vXdrMMB/nQLgAMi4q+AXwFndbmeKpYBrwcWd7uQZjL92oVLgOndLmKINgJnRsSLgYOBd4zx/fw4cHhEHAhMBaZLOrjLNVV1OrCi20U4wEsi4saI2Jh6f0xxb/uYFhErImKsfrq1JruvXYiIxcAj3a5jKCJiVUT8NHWvpwiYPbpbVXNR2JB6J6bHmL+rQtIU4Bjgi92uxQHe3FuAb3e7iC1Eo69dGLPBsiWQ1AscBNza3UoGl05FLAXWAAsiYkzXm5wPvBfY1O1Cxt0POkj6HvCcBqM+EBHXpmk+QPF29PLRrK2ZKjWPcZW+dsHaQ9L2wNeBMyJiXbfrGUxEPAFMTdebrpF0QESM2esOko4F1kTEEkmHdbuecRfgEXHkYOMlzQSOBY6IMXKTfKuaM+CvXRglkiZShPflEXF1t+upKiLWSlpEcd1hzAY4cAjwOklHA1sDO0q6LCLe3I1ifAqlRNJ04H3A6yLisW7XswXx1y6MAkkCLgJWRMS53a6nFUk9tTu9JG0DHAnc2d2qBhcRZ0XElIjopTiOb+pWeIMDvN5ngR2ABZKWSvpCtwtqRdLxku4DXg7cIOm73a6pXrowXPvahRXAVWP9axckXQncArxI0n2STul2TRUcApwMHJ6O36WppThWTQa+L+kOihf5BRHR1dvycuOP0puZZcotcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8vU/wNHlagK/NEd0wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "resid = res.resid_deviance.copy()\n",
    "resid_std = stats.zscore(resid)\n",
    "ax.hist(resid_std, bins=25)\n",
    "ax.set_title('Histogram of standardized deviance residuals');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QQ Plot of Deviance Residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3gU5drH8e9NaCIoElA5agARGx6xYC/HfmxHRUVFFASPCBHBglLynqNHDSAgig2IFAVWsGBFUbH3AooKoghKELFAEEEBE5Pn/WMmuGQ3ySRkdjfJ73NduXZndmb2Xsve+5S5H3POISIiEq1OsgMQEZHUo+QgIiIxlBxERCSGkoOIiMRQchARkRh1kx1AVWjevLlr3bp1ssMQEalW5s2bt9o51yLeazUiObRu3Zq5c+cmOwwRkWrFzHJLe03dSiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIxlBxERGqoSARat4Y6dbzHSCT4uTViKquIiGwpEoFevWDDBm87N9fbBujatfzz1XIQEamBsrL+SgzFNmzw9geR1ORgZpPM7GczWxC172Yz+97M5vt/pyczRhGR6mj58ortLynZLYcHgVPj7L/TOXeA//d8gmMSEan2MjIqtr+kpCYH59ybwJpkxiAiUhNlZ0OjRlvua9TI2x9EslsOpelrZp/53U47xDvAzHqZ2Vwzm7tq1apExyciktK6doWcHGjVCsy8x5ycYIPRAJbsNaTNrDUwyzm3n7+9E7AacMCtQEvnXM+yrtGxY0enwnsiIhVjZvOccx3jvZZyLQfn3E/OuULnXBHwAHBosmMSEaltUi45mFnLqM1OwILSjhURkXAk9SY4M5sOHAc0N7MVwE3AcWZ2AF630jLgyqQFKCJSSyU1OTjnusTZPTHhgYiIyBZSrltJRESST8lBRERiKDmIiEgMJQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOIiIpbmvWgq4srSEtIpLCtnYt6MpSy0FEJIVt7VrQlaXkICKSwrZ2LejKUnIQEUlhW7sWdGUpOYiIpLCtXQu6spQcRERS2NauBV1Zmq0kIpLiunYNPxmUpJaDiIjEUHIQEZEYSg4iIhJDyUFERGIoOYiISAwlBxERiaHkICIiMZQcREQkhpKDiIjEUHIQEUmiZCzkE4TKZ4iIJEmyFvIJQi0HEZEkSdZCPkEoOYiIhKy0rqNkLeQThLqVRERCVFbXUUaGt11S2Av5BJHUloOZTTKzn81sQdS+ZmY2x8y+9h93SGaMIiJbo6yuo2Qt5BNEsruVHgROLbFvEPCKc64d8Iq/LSJSLZXVdZSshXyCSGpycM69Cawpsfts4CH/+UPAOQkNSkSkCpW3BnTXrrBsGRQVeY+pkBgg+S2HeHZyzv0A4D/uGO8gM+tlZnPNbO6qVasSGqCISGlKDj6ffnoSu46cg9mz4ccfK3xqKiaHQJxzOc65js65ji1atEh2OCIimwefc3O97+XcXHjoIejePcFdR0VFMHMmdOzoZaecnApfIhVnK/1kZi2dcz+YWUvg52QHJCISRGmDz88/73UZha6gAKZPh2HD4MsvoV07mDSpUpkoFVsOzwDd/efdgaeTGIuISGBJu29h0yYYNw723NNrptSvDzNmwKJF0KOHt11ByZ7KOh14D9jLzFaY2eXAcOBkM/saONnfFhFJeeUNPle533+H0aNh992hTx/YaSd45hmYPx8uvBDS0ip96aR2KznnupTy0okJDUREpApkZ295wxuENPi8di3ccw+MGQN5eXDCCTB1qvdoViVvkYrdSiIi1VLo9y38/DMMHuw1Rf77XzjiCHjvPXjlFTjxxCpLDJCaA9IiItVW164hzET67jsYORIeeAD++AM6d4YhQ6BDhyp+o78oOYiIpKolS2D4cJgyxZsbe+mlMGiQN/AcMnUriYhshVAW61mwAC6+GPbaC6ZN8wYylizxpqUmIDGAWg4iIhUSiXj3MyxfDs2awfr1kJ/vvbbVi/V8+CEMHQpPPw2NG8OAAXDttbDzzlUWf1BqOYiIBFTyDui8vL8SQ7EKL9bjHLz+Opx8Mhx2GLz5Jtx8s/cmt9+elMQAajmIiAQW7w7oeALd9FZc9yg7G95917tHYcQI6N0bmjTZ6li3lpKDiEhAQe90LvOmt8JCeOIJr/to/nzv4HvvhZ49YZttqiTOqqBuJRGRgILc6VzqTW8FBV4Vvv32gwsu8Jogkyd7A81XXZVSiQECJAcza2tmDfznx5lZPzNrGn5oIiKpJd7KbfXqQXp6GTe9bdoEY8d6s4wuuwwaNIBHHoEvvvC269VL4CcILkjLYSZQaGZ7ABOBNsDDoUYlIpJCiqerXnqp9wM/OhlMngyrV8dZrOe332DUKGjTBjIzvYHlZ5+FTz7xWg5bUfcoEYKMORQ55/40s07AXc65e8zsk7ADExFJlrKmq+blea2HqVNLma76yy9/1T1as8Yra/Hww3DccVVa3iJsQVoOBWbWBa989ix/X2q2g0REtlKlp6v+9JN393KrVnDTTXDUUfD++/Dyy3D88dUqMUCwlkMPoDeQ7Zz71szaANPCDUtEJDkqPF11+XKv7tGECV4WueACrzje/vuHGmfYyk0OzrkvzGwgkOFvf4vWWBCRGirodNVjW34Nl/t1jwC6dfNaDu3ahRdcAgWZrfQvYD7wgr99gJk9E3ZgIiLJUN501f34nEfTuvDqD3t7Ywm9e8PSpTBxYo1JDBBszOFm4FBgLYBzbj7ejCURkWovunBe8+bezKOS6tWDk7f7gKc5i8/Zn3PqzaLOjTd405PuuSfEpd6SJ8iYw5/OuV9ty8EUF1I8IiIJUzz4XDzGkJdX8gjH2du9ztjdsmm58BVv6lL//1Hv6qthhx0SHW5CBUkOC8zsYiDNzNoB/YB3ww1LRCR8pQ8+O87gOYYwlCPXvQd5O3v3LFx5pVcttRYI0q10NdAe+AOYDqwDrgkzKBGRsEUi3nTVaHUopDOP8gkHMot/8TdWksn98O23cP31tSYxQLDZShuALP9PRKTaK+5OKlaXAroSYRDD2Zuv+JK96M6DPMzF7NKqHjRMXqzJUmpyMLNnKWNswTl3VigRiYiEKBKB7t294qgN2UhPJnEjI2jFcj7hADrzKE9wLkWklV5ErxYoq+UwKmFRiIgkQHGLYZvC9fRmHNdzBzvzE+9wJH0Yy2xOIz3dcGugVYaXGCq1olsNUGpycM69kchARETCNmLQGgZsuIf+jKEZvzCHk7iIGbzBPwCjVStvdqqU3a30qHPuAjP7nDjdS8656n1vuIjUGjPv+5FVQ0bz9rqxNOE3nuJshjKEjzh08zG1uQspnrK6lfr7j2cmIhARkSqXm8vrZ47k9AUTqU8+j3AhwxjMAv6+xWFpaXHWYajlSp3K6pz7wX+a6ZzLjf4DMhMTnohIJSxezNJje1DQeg+OXJBDhK7sxVd05eGYxNCokbdAmxLDloLc53BynH2nVXUgIiKVFYl4pS862Kc8YhdStNfetHzrEe4nk7Ys5QomsJQ9Ys6Lu3KbAGWPOfTBayHsbmafRb3UBHgn7MBERMoTiUD//rBH3vtMJpt/MYt1NOF2BnIn17KKHUs9V4PPZStrzOFhYDYwDBgUtX+9c25NqFEBZrYMWA8U4tV36hj2e4pI9ZHZx7F43Ks8wlBO5FXyaMZ/uIV76ctayq57ZKbB5/KUNZX1V+BXoIuZpQE7+cc3NrPGzrmAVc+3yvHOuTg1EkWktopMczyfOYur12dzOB+wkpZcxx3k0IvfKb+8hZlXZVtdSWUrt3yGmfXFK9v9E1Dk73aAprKKSOIUFjLxtMfpOGcoXfmMb2lNb8byIJfxR8D6Funp3tLOSgzlC1KV9RpgL+dcTDHbkDngJTNzwHjnXE6C319EUkF+Pu/1ncaOE4dzedHXLGJvuvEQ0+nCnwGXs1dSqLggyeE7vO6lRDvKObfSzHYE5pjZl865N4tfNLNeQC+AjBq40IZIrbdxI0ycyJrBIzjit+/4mAM5j8d5kk64ciZa1qkDRUXeoHNtLoGxNYIkh2+A183sObyy3QA450aHFpV3/ZX+489m9iTeanRvRr2eA+QAdOzYUYsPiVRzxTOP8vPW0ZtxXMdoduYnvuAoshnPC5wKWKnnF48l3H9/4mKuyYIkh+X+X33/L3Rmti1Qxzm33n9+CnBLIt5bRBIrEvHW0Gnwex79uJt+3M0OrOVFTuECsniLY8u9hrqNql6Q9Rz+l4hAStgJeNJfmrQu8LBz7oUkxCEiIcrMhCfH/sBNjKYPY2nM7zzJOQxlCHM5pNzz1VoIT5DZSi2AG/FWg9s8JcA5d0JYQTnnvgE6hHV9EUm+/7s0l/bTRjCaidSjgBlcxDAGs5D9Ap2v1kK4gnQrRYBH8Arw9Qa6A6vCDEpEaq5nRn7FuiHDuOnPCA7jQS5jBDfGLW9Rmj591FoIW5DkkO6cm2hm/f01Ht4wM631ICIVM38+884fyplLH2cTDbmPqxjFAL5n1wpdRokhMYIkhwL/8QczOwNYCRX8tykitde773rzSZ9/nj3YjuEM4i6uKbPuUTzqRkqsIMnhNjPbHrgeuAfYDrg21KhEpFqLTHPMzHyFq9dnczyvs5p07uQ27uMqfqVpmec2bgzjxikJJFuQ2Uqz/Ke/AseHG46IVGtFRbw+YBbt7srmCfch3/M3rmU0OfRiA9uWe7q6jFJHkNlKk4m/TGjPUCISkeqnsJC3+z3KDuOHcVzh53xDG3oxnofoTj4NAl1CiSG1BOlWmhX1vCHQCW/cQURqu/x8pv1zKoe9PpyjWcIX7MMlTGUGF1EY6OtF3UipKki30szobTObDrwcWkQikvJmTN7I/L4TyNwwkkv4jnkcxLnM5CnOKbfuUTS1FlJXsNS+pXaAKt2J1Ebr1jHhoPv519I7uYifeYuj6UUOL/JPyqp7FI8SQ2oLMuawHm/MwfzHH4GBIcclIinksXF5LO03hisL7uHfrOUF/kk2WbzNMRW+lqakVg9BupWaJCIQEUktkQjcdMVKem8cTW/G0Zjfmcm5DGMw86jYqr1qJVQ/ZSYHM9sG6Ars6++aCzzunMsPOzARSbzistmN85YxkNtZyCTSKGQ6XRjGYBZt/ioIpk4dr+KqEkP1U+rIkZn9HVgEHAMsA3KBfwLvmFlTM7stIRGKSKgyM73qpmZw6yVfckded5awBz2ZxINcxp4sphtTAyeGhg1h2jRwDgoLlRiqq7JaDncDVzjn5kTvNLOTgAXAwjADE5HwFK+h8Pvv3vYBfMIQhnIeM9lEQ+7hakYxgJXsEviaGkuoWcpKDi1LJgYA59zLZlaAd7+DiFQjkQj07An5fsfwkbxDFtmczmx+ZTuGMoQx9Gc1LQJf88QT4WVNbq9xypqQXMfMYm5tNLOGQIFzbkN4YYlIVYtE4NJLIT/fcRJzeI3jeIejOYSPGEI2GSznP9wWODE0bux1Hykx1ExlJYcpwEwza128w3/+KDA1zKBEpGpFItDtkiLOck/xAYcxh1PYgyVcw520ZhnDGMI6ti/3OsUJwTlYv15dSDVZqcnBOXcb8ALwppmtNrPVwBvAHOfcrYkKUEQqLxKBRvX/5LlLHmY+HXiKTqSTxxXk0JaljOGacgvimXlTUZUQapcyp7I65+4F7jWzJv72+oREJSJbJTMTJozNpxtT+Izh7MFSFrIvXZnGI1wYqO6RBphrt0DlM5QURKqPg/fZwNFfPsBSRrEbK/iIjpzDkzzDWWXWPVIBPIlWmdpKIpJiIhG4odevdN9wP7O5kx1ZxRscy+VMZA4nU1bdI802kniUHESquXOPXc2Bb43hC+6hKb8ym1PJJot3OLrcc1XWQkoTpPBeI7wlQjOcc1eYWTtgr6gV4kQkCQZ3X0mLKaOYyni2YSNPcC5DGcInHBTofCUGKUuQlsNkYB5whL+9AniMLRcBEpEEyMyE2WO/ZSC3czOTSaOQh7mY4QyqUHmLCRM0tiBlC7IqR1vn3AigAMA5t5GKFm4XkUqLRLzB4n3tCw4f242vaUcPJjOZHuzJYrozJVBiSE/37lHYuFGJQcoXpOWQ71dndQBm1hb4I9SoRGRzqYv2+R/zENl04kk2sg1j6M9orgtc92jffWGhKqFJBQVJDjfh3Qy3m5lFgKOAy8IMSqS2i0Rg3CVv8xTZnMYLrGV7ssliDP3Jo3mga5hB794aV5DKCbLYzxwz+xg4HK87qb9zbnXokYnUQpl9HEvHvUQW2bzFW/xMCwYzlPvJDFTeAtRSkKpRanIws5JTHn7wHzPMLMM593F4YYnUHpEI9O5VxEkbnmYIQzmEuXzHrvRjDBP4NxtpFPhamoEkVaWslsMdZbzmgBOqOBaRWqdv7z/5ZfwjvMcw9mMhS2jLv3mAKXSjgPqBrqHuIwlDqcnBOXd8IgMpycxOBcYAacAE59zwZMYjUlUyM2Hi2D/ozkMM5Hba8g0LaM/FRHiUCwLVPQK1EiRc5U5lNbOGZnadmT1hZjPN7Bp/TYfQmFkacB9wGt761V3MrGKL14qkmEgEtkv7nfpj72IpbcnhSvJI52yeYn8+YzoXl5sYopfgVGKQMAX5iTIFWA/c4293wVvPoXNYQQGHAkucc98AmNkM4GzgixDfU6RKRS/FuR2/chX3sZQ7acFqXucf9GAyL3MSQW8bUg0kSaQgyWEv51yHqO3XzOzTsALy7QJ8F7W9Ajgs+gAz6wX0AsjIyAg5HJHgopfibM4qBnMXfbmX7VnHc5zOUIbwLkcFvp7GFCQZgtwh/YmZHV68YWaHAe+EF5L3NnH2uS02nMtxznV0znVs0SL4erciVa34DmYz7++SS6B5/veM5lqW0ZrBDOMlTuFAPuZMngucGIoX2CkqUmKQxAvScjgM6GZmy/3tDGCRmX0OOOfc/iHEtQLYLWp7V2BlCO8jUmnRLYRiu7OUGxnBZTxIGoVM4xJuZyBfsk/g62qgWVJBkORwauhRxPoIaGdmbYDvgYuAi5MQh0hcJ50Er7zy1/a+LGQww+jCdAqox0QuZwQ3kkvrwNfUYjuSSoLcIZ1rZjvg/ZKvG7U/tJvgnHN/mllf4EW8qayTnHO651OSLhKB7t2hsNDbPoh5ZJHNuTzJb2zLnVzLHVzPj7QMdD0NMkuqCrKew614tZSW8le/f+g3wTnnngeeD/M9RILIzISxY7fcdwxvMoShnMqL/EJTbuE/jKE/a0gPdE21EiTVBelWugCvbHd+uUeK1CAlWwng+CcvkkU2x/A2P7EjAxnOWPqwnu3KvJbWUJDqJkhyWAA0BX4OORaRlBE9pmAUcQ5PMYShdGQey9mNq7mbiVxebt0jtRCkugqSHIbhTWddQNQ6Ds65s0KLSiTBom9YK5bGn3RhOoMZxr4s4mv2oCcTmcYlpdY9UgtBaoogyeEh4Hbgc6Ao3HBEEq/kzKP6/MFlPMhAbmd3vuUz/s5FTOcxOlNEWtxr1K0LDz6opCA1R5DksNo5d3fokYgkWMkxhUb8Ti9yGMAodmElH3Ao13AXszgTV8b9oppxJDVRkOQwz8yGAc+wZbeS1nOQaicz0xsDcFH322/PWvpyL9dwF83J41WOpxtTeJUTKKvukcYTpCYLkhwO9B8Pj9qn9RykWomdeeTVPbqWO7mK+9iedcziDLLJ4n2OKPU6aiVIbRHkJrikrusgsjXi3aOwCysYwCh6kUNDNvE45zOUIXzKAaVeR2MKUtsEWlXEzM4A2gOb13Fwzt0SVlAiWyteUmjLEgZyO915CMMxjUsYziAWs1ep11FFVKmtgtwhPQ5oBBwPTADOBz4MOS6RCos3HRWgPQsYzDAuYgYF1OMBrmAkN5Ra90itBJFgLYcjnXP7m9lnzrn/mdkdwBNhByYSVLxWAkBHPmIIQ+nEU6ynMXdwPXdybZl1jzSmIOIJsp7DRv9xg5n9DSgA2oQXkkgwkQg0aFAyMTiO5Q1e5BQ+4lD+wRvczE20IpeBjCg1MTRu7C2/qcQg4gnScphlZk2BkcDHeDOVHgg1KpFylLxxDRyn8gJZZHM07/AjO3EjtzOWPvxGk7jXUCtBpHTlthycc7c659Y652YCrYC9nXP/DT80kViZmd4gcXTdo/N4nHkczGxOJ4Pl9OUe2vAtI7kxbmJQK0GkfKUmBzM7xMx2jtruBjwK3GpmzRIRnAhsuQxncRdSXQq4lCkspD2P05nG/EYPJrEHS7iPvmximy2uUbeulxCcg/XrNdgsUp6yWg7jgXwAMzsWGA5MAX4FcsIPTWqreGsyF89AasAmrmQci9mTKXQnn/pcyAz2YREP0iOmIJ6Zt+xmQYESgkhFlDXmkOacW+M/vxDI8buWZprZ/PBDk9qotJlH2/IbVzKe67mDv/ED73MY/bibWZxJvBIXWodZZOuUmRzMrK5z7k/gRKBXwPNEAivt3oRiTfmFvtxLf8bQnDxe4QQuYRqvcTxKCiLhKetLfjrwhpmtxpvO+haAme2B17UkUmnlJYUW/Ly57tF2rOdZziSbLD7YosTXX7SOgkjVKjU5OOeyzewVoCXwknOb61jWAa5ORHBSM5XWdQSwK99xAyO5ggdowB88ygUMYzCf0SHu8UoKIuEos3vIOfd+nH2LwwtHaqp4pbKjtWUJgxhON6ZgOKZyKcMZxNfsGfd4lcsWCZfGDiR0ZbUU9uNzBjOMC3mEAuqRQy9GcgPLaRVzbJ06XleUxhREwqfkIKGKROInhkP4kCyyOZtnWE9jRjGAO7mWn9h5i+PS0uChh9RCEEm0ILWVRALLzPR+4Uffo/AXxz94nZc4mQ85jGN4i5u4mVbkMojbYxJDeroSg0iyqOUgVSISgZ49IT8/3quO03meIQzlKN7lB3ZmACMZz5VblLfQNFSR1KHkIJUSiUD//pCXV/oxdSjkXJ5gCEM5kPnkkkEm9zGJnvzx17pRGlwWSUHqVpIKy8z0uotKSwx1KaAbD7GQ9jzGBTRiA5cxmT1YwlgyNyeG9HSv3pFqHYmkHrUcJJAgLYUGbKInk7iREbQml/l0oDOP8gTnUkTa5uPq14dJk5QQRFKZkoOUKxKBHj284nXxbMtv9GYc13MHLfmRdzmCq7iP5zmdkiUu1IUkUj2kXLeSmd1sZt+b2Xz/7/Rkx1QbRSLQvPlfM47iJYYdWMN/uIVcWjGKG1hIe47nVY7iHZ7nDKITg7qQRKqXVG053OmcG5XsIGqbIF1HADvyE9cxmkzupwm/8TRnMZQhfMhhm49RC0GkekvV5CAJVtZdzMV2Yzk3MJJ/M4H65G+ue/Q5+28+Jj0dxoxRUhCp7lKuW8nX18w+M7NJZrZDsoOpyYq7j8pKDO1YzER6spS29GYcD3Mxe/MlFzOdlen7b15hzTlYvVqJQaQmSErLwcxehhK3w3qygLHArYDzH+8Aesa5Ri/8NSYyMjJCi7UmK68Y3t/5jCEMpTOPkU99xtKHUQzgOzKoUwemTVEiEKmpzJX2zZACzKw1MMs5t19Zx3Xs2NHNnTs3ITHVFJEIXHpp/MRwKB+QRTZn8SzraML9ZHIn1/IzOwGaiipSU5jZPOdcx3ivpVy3kpm1jNrsBCxIViw1VSQC3buXTAyO43mVOZzEBxzOUbzDf7iFVuQymOGbE0N6uhKDSG2QigPSI8zsALxupWXAlckNp2aJ7UpynMFzZJHNEbzPD+zM9YxiPFdijRtrxpFILZVyycE5d2myY6gpIhHIyoLcXO9+heiWQh0KOY+ZDGEoB/Apy2hFH+5nMj1onN6Q8ZpxJFKrpVy3klRcJAKtW3sJoG5d77FOHe/mtdxc75jixFCXAi5jMl+wL49yIQ3ZRHcepB1fM9760LNPQ804EpHUazlIcPFuWiss9B5LDjQ3ZOPmuketWM4nHMD5PMaTdKKINNLSYKrWThARn5JDNVNWV1E8jVm/ue7RzvzEOxxJH8Yym9MoLm9hpkV1RGRLSg7VRLxWQlmJYQfW0I+76cfdNOMXXuJkLiSLNzmW6JpHZtC7txKDiGxJyaEaiESgVy/YsKH8Y3fiR65jNH0YSxN+4ynOZihD+IhDY45VqQsRKY2SQ4op7jZavhyaNfP2lVcIDyCD3M11j+pRwCNcyDAGs4C/xxyrpCAi5VFySAGljSMESQp78hWDGM4lTMNhTKE7wxnIsrQ9KCyEVq0gO1uJQEQqRskhyUp2GQWtZrI/n26ue7SJhkxrkknT2wbw73678e/wwhWRWkLJIQmiu47q1Plr+mkQh/MeWWRzJs+x3pqw6MyBtJ9wLT123DG8gEWk1lFySLCSLYVgicFxAq+SRTYn8Bpr6qTz6bm30uGBvrRv2jTMcEWkltId0gmWlRVs1pHHcSbP8h5H8AonsY99ybyL76DZr8vo8Nj/gRKDiIREySHBli8v/5g0CrmQGXyedgDPchY78jNDmo3j9YnfcHDkOm8NThGRECk5hKy47lGdOt5j8fTUktLSoD75DEifxJqd92EGXdhvzwKYMoXdCxYzNO9KuvRomMjQRaQW05hDiEqOL+TmQr163mI5+fl/Hddsm4280Hkih7w2Ar77Dg48EO59HDp18rKKiEiC6ZsnRPHGFwoKoEkT7/6D7VjHsKa3s6Jeaw6ZcjVkZMDzz8O8eXDeeUoMIpI0ajmEqNTxhbw8ll11N9x9N6xdC6ec4mWSY49NaHwiIqXRT9MqVN74ws78wAhuINdawS23wPHHw4cfwosvKjGISEpRy6GKlDW+0DJ/GTcygp5Moh4FLD+iC21yBkP79skNWkSkFEoOVSTe+MLuBV9yU4PhnE8Eh/F448vY9n8DOfu6tskJUkQkICWHKhI9vtCB+WSRzXnMZNMfDanX/yoYMICLd901eQGKiFSAxhwqoeTYQiTiTTQ6gneZxRnM50BO4SWGMZijd82Fu+4CJQYRqUbUcqig2LEFx/TLX+Gl9Gz25HVWk04Wt3EfV1HQqCk5w5Mbr4hIZajlUEHFYwtGEWfxNO9zOLP+OJntflrMvK6jOWa3XIZZFk1bNSUnR+soiEj1pORQQrwuo2grcgu5iOl8Sgee5hxasIpejKd14TccPO1aFpOBQicAAAlrSURBVC3flqIiWLZMiUFEqi91K0WJNx21Vy/vedfO+TB1Kl/XHU6bP5ewkH25hKnM4CIKqUurVsmLW0Skqik5RIk3HdVt2MDXV0+AQSNhxQq2a3MwF33/BI/mn43zG16NGnlLcYqI1BQ1tlupvO6heKKnozZhHQMZzjJac/Mv/aFNG3jhBdKXfsS/JnUio1UdzLwaSRpbEJGaxlzQRYtTWMeOHd3cuXM3b5fsHgLv1315X+KtW8Nvuavpx91czT3swFpmcyqTdhrCYz8eE94HEBFJAjOb55zrGO+1GtlyiNc9tGGDt79UK1cye9/rWUZr/sutvMoJHMxczm80m3PuUGIQkdqlRiaH0qqhxt3/7bfQpw+0acM+L41h1VGdOKnlQjrbTPJaHawuIxGplZKSHMyss5ktNLMiM+tY4rXBZrbEzL4ys39W5voZGQH2L1oE3btDu3YwaRJcdhksXkybt6fy8sp9NR1VRGq1ZLUcFgDnAm9G7zSzfYGLgPbAqcD9ZpZW0YtnZ3tjDNE2zyj65BM4/3yvIurjj0O/fvDNNzB+POy+eyU/johIzZKUqazOuUUAZlbypbOBGc65P4BvzWwJcCjwXkWuX/xrPyvL60rKyICc7u9wSiQbZs+G7beHIUOgf39o0WJrP46ISI2TamMOuwDfRW2v8PfFMLNeZjbXzOauWrUq5vWuXWHZt46iF15iWat/cMotR8NHH8HQod7dbbfdpsQgIlKK0FoOZvYysHOcl7Kcc0+XdlqcfXHn2jrncoAc8KayxhyQlwenneYlhF128SqjXnFFbH+TiIjECC05OOdOqsRpK4DdorZ3BVZWKoBmzbwb1664Arp1gwYNKnUZEZHaKNXKZzwDPGxmo4G/Ae2ADyt1JTN45JEqDE1EpPZI1lTWTma2AjgCeM7MXgRwzi0EHgW+AF4ArnLOFSYjRhGR2ixZs5WeBJ4s5bVsQGXsRESSKNVmK4mISApQchARkRhKDiIiEkPJQUREYig5iIhIDCUHERGJUSNWgjOzVUBusuPYSs2B1ckOIkFq02eF2vV59Vmrl1bOubhF5mpEcqgJzGxuacv11TS16bNC7fq8+qw1h7qVREQkhpKDiIjEUHJIHTnJDiCBatNnhdr1efVZawiNOYiISAy1HEREJIaSg4iIxFBySCFmNtLMvjSzz8zsSTNrmuyYwmJmnc1soZkVmVmNnA5oZqea2VdmtsTMBiU7njCZ2SQz+9nMFiQ7lrCZ2W5m9pqZLfL/G+6f7JjCoOSQWuYA+znn9gcWA4OTHE+YFgDnAm8mO5AwmFkacB9wGrAv0MXM9k1uVKF6EDg12UEkyJ/A9c65fYDDgatq4r9bJYcU4px7yTn3p7/5Pt4a2jWSc26Rc+6rZMcRokOBJc65b5xz+cAM4OwkxxQa59ybwJpkx5EIzrkfnHMf+8/XA4uAXZIbVdVTckhdPYHZyQ5CKm0X4Luo7RXUwC+Q2s7MWgMHAh8kN5Kql5RlQmszM3sZ2DnOS1nOuaf9Y7Lwmq6RRMZW1YJ81hrM4uzTvPEaxMwaAzOBa5xz65IdT1VTckgw59xJZb1uZt2BM4ETXTW/CaW8z1rDrQB2i9reFViZpFikiplZPbzEEHHOPZHseMKgbqUUYmanAgOBs5xzG5Idj2yVj4B2ZtbGzOoDFwHPJDkmqQJmZsBEYJFzbnSy4wmLkkNquRdoAswxs/lmNi7ZAYXFzDqZ2QrgCOA5M3sx2TFVJX9iQV/gRbwBy0edcwuTG1V4zGw68B6wl5mtMLPLkx1TiI4CLgVO8P8/nW9mpyc7qKqm8hkiIhJDLQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOkjLMLD1qauCPZva9/3ytmX2R4FjOiS6mZma3mFmFb+ozs9alVSo1s/Zm9qqZLTazpWb2PzOr8v8ny/osZvZ6Ta2KK1tHyUFShnMuzzl3gHPuAGAccKf//ACgqKrfz8zKqhBwDl411eLY/uuce7kK33sbvJvihjvn9gT+jlesL4zyz6F+FqmZlBykukgzswf8+vkv+V+umFlbM3vBzOaZ2Vtmtre/v5WZveKvjfGKmWX4+x80s9Fm9hpwe7zzzexI4CxgpN9yaeufd75/jUPM7F0z+9TMPjSzJn4L4S0z+9j/O7Kcz3Mx8I5z7iUA/474vsAN/nvcbGYDig82swV+kTfM7Ck/3oVm1ivqmN/MLNuP630z26m8zxLNzE4xs/f8+B/zawdhZsPN7Av/n+WoCv+bk2pJyUGqi3bAfc659sBa4Dx/fw5wtXPuYGAAcL+//15gir82RgS4O+paewInOeeuj3e+c+5dvF/1N/gtmaXFJ/qlMB4B+jvnOgAnARuBn4GTnXMHAReWeL942gPzonf477ONlb/IU08/3o5APzNL9/dvC7zvx/UmcEVZnyWamTUH/s//53IQMBe4zsyaAZ2A9v4/y9vKiU1qCBXek+riW+fcfP/5PKC1/8v2SOAxr9wNAA38xyPwFhMCmAqMiLrWY865wnLOL81ewA/OuY8Aiqtxmtm2wL1mdgBQiJeAymLEr9Iar5prSf3MrJP/fDe8xJkH5AOz/P3zgJMDXKvY4XhdT+/4/yzq45XDWAdsAiaY2XNR15caTslBqos/op4XAtvgtXzX+uMS5Yn+Iv7df6zI+cVK+1K/FvgJ6OBfd1M511kIHLvFhc12B1Y759aa2Z9s2bJv6B9zHF5r5Qjn3AYze734NaAgqpJvIRX7/9uAOc65LjEvmB0KnIhXPLAvcEIFrivVlLqVpNryf7V/a2adwauWaWYd/JffxfsyA+gKvF3B89fjFUEs6Uvgb2Z2iH9OE39ge3u8FkURXlG2tHLCjwBHR80a2gavK+om//VlwEH+awcBbfz92wO/+Ilhb7xf/OUp7bNEex84ysz28N+zkZnt6beutnfOPQ9cgzc5QGoBJQep7roCl5vZp3i/xouX4uwH9DCzz/C+rEubBVTa+TOAG8zsEzNrW3ywv+TnhcA9/jlz8H653w90N7P38bqUfqcMzrmNeAPFWWa2GFiNN0BdvMDTTKCZmc0H+uCtKQ7wAlDX/1y34n2plyfuZykRzyrgMmC6f+33gb3xksosf98beC0kqQVUlVUkBZjZOcBo4HjnXG6y4xFRchARkRjqVhIRkRhKDiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIx/h/Ny7kFWvsEHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3gU5drH8e9NaCIoElA5agARGx6xYC/HfmxHRUVFFASPCBHBglLynqNHDSAgig2IFAVWsGBFUbH3AooKoghKELFAEEEBE5Pn/WMmuGQ3ySRkdjfJ73NduXZndmb2Xsve+5S5H3POISIiEq1OsgMQEZHUo+QgIiIxlBxERCSGkoOIiMRQchARkRh1kx1AVWjevLlr3bp1ssMQEalW5s2bt9o51yLeazUiObRu3Zq5c+cmOwwRkWrFzHJLe03dSiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIxlBxERGqoSARat4Y6dbzHSCT4uTViKquIiGwpEoFevWDDBm87N9fbBujatfzz1XIQEamBsrL+SgzFNmzw9geR1ORgZpPM7GczWxC172Yz+97M5vt/pyczRhGR6mj58ortLynZLYcHgVPj7L/TOXeA//d8gmMSEan2MjIqtr+kpCYH59ybwJpkxiAiUhNlZ0OjRlvua9TI2x9EslsOpelrZp/53U47xDvAzHqZ2Vwzm7tq1apExyciktK6doWcHGjVCsy8x5ycYIPRAJbsNaTNrDUwyzm3n7+9E7AacMCtQEvnXM+yrtGxY0enwnsiIhVjZvOccx3jvZZyLQfn3E/OuULnXBHwAHBosmMSEaltUi45mFnLqM1OwILSjhURkXAk9SY4M5sOHAc0N7MVwE3AcWZ2AF630jLgyqQFKCJSSyU1OTjnusTZPTHhgYiIyBZSrltJRESST8lBRERiKDmIiEgMJQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOIiIpbmvWgq4srSEtIpLCtnYt6MpSy0FEJIVt7VrQlaXkICKSwrZ2LejKUnIQEUlhW7sWdGUpOYiIpLCtXQu6spQcRERS2NauBV1Zmq0kIpLiunYNPxmUpJaDiIjEUHIQEZEYSg4iIhJDyUFERGIoOYiISAwlBxERiaHkICIiMZQcREQkhpKDiIjEUHIQEUmiZCzkE4TKZ4iIJEmyFvIJQi0HEZEkSdZCPkEoOYiIhKy0rqNkLeQThLqVRERCVFbXUUaGt11S2Av5BJHUloOZTTKzn81sQdS+ZmY2x8y+9h93SGaMIiJbo6yuo2Qt5BNEsruVHgROLbFvEPCKc64d8Iq/LSJSLZXVdZSshXyCSGpycM69Cawpsfts4CH/+UPAOQkNSkSkCpW3BnTXrrBsGRQVeY+pkBgg+S2HeHZyzv0A4D/uGO8gM+tlZnPNbO6qVasSGqCISGlKDj6ffnoSu46cg9mz4ccfK3xqKiaHQJxzOc65js65ji1atEh2OCIimwefc3O97+XcXHjoIejePcFdR0VFMHMmdOzoZaecnApfIhVnK/1kZi2dcz+YWUvg52QHJCISRGmDz88/73UZha6gAKZPh2HD4MsvoV07mDSpUpkoFVsOzwDd/efdgaeTGIuISGBJu29h0yYYNw723NNrptSvDzNmwKJF0KOHt11ByZ7KOh14D9jLzFaY2eXAcOBkM/saONnfFhFJeeUNPle533+H0aNh992hTx/YaSd45hmYPx8uvBDS0ip96aR2KznnupTy0okJDUREpApkZ295wxuENPi8di3ccw+MGQN5eXDCCTB1qvdoViVvkYrdSiIi1VLo9y38/DMMHuw1Rf77XzjiCHjvPXjlFTjxxCpLDJCaA9IiItVW164hzET67jsYORIeeAD++AM6d4YhQ6BDhyp+o78oOYiIpKolS2D4cJgyxZsbe+mlMGiQN/AcMnUriYhshVAW61mwAC6+GPbaC6ZN8wYylizxpqUmIDGAWg4iIhUSiXj3MyxfDs2awfr1kJ/vvbbVi/V8+CEMHQpPPw2NG8OAAXDttbDzzlUWf1BqOYiIBFTyDui8vL8SQ7EKL9bjHLz+Opx8Mhx2GLz5Jtx8s/cmt9+elMQAajmIiAQW7w7oeALd9FZc9yg7G95917tHYcQI6N0bmjTZ6li3lpKDiEhAQe90LvOmt8JCeOIJr/to/nzv4HvvhZ49YZttqiTOqqBuJRGRgILc6VzqTW8FBV4Vvv32gwsu8Jogkyd7A81XXZVSiQECJAcza2tmDfznx5lZPzNrGn5oIiKpJd7KbfXqQXp6GTe9bdoEY8d6s4wuuwwaNIBHHoEvvvC269VL4CcILkjLYSZQaGZ7ABOBNsDDoUYlIpJCiqerXnqp9wM/OhlMngyrV8dZrOe332DUKGjTBjIzvYHlZ5+FTz7xWg5bUfcoEYKMORQ55/40s07AXc65e8zsk7ADExFJlrKmq+blea2HqVNLma76yy9/1T1as8Yra/Hww3DccVVa3iJsQVoOBWbWBa989ix/X2q2g0REtlKlp6v+9JN393KrVnDTTXDUUfD++/Dyy3D88dUqMUCwlkMPoDeQ7Zz71szaANPCDUtEJDkqPF11+XKv7tGECV4WueACrzje/vuHGmfYyk0OzrkvzGwgkOFvf4vWWBCRGirodNVjW34Nl/t1jwC6dfNaDu3ahRdcAgWZrfQvYD7wgr99gJk9E3ZgIiLJUN501f34nEfTuvDqD3t7Ywm9e8PSpTBxYo1JDBBszOFm4FBgLYBzbj7ejCURkWovunBe8+bezKOS6tWDk7f7gKc5i8/Zn3PqzaLOjTd405PuuSfEpd6SJ8iYw5/OuV9ty8EUF1I8IiIJUzz4XDzGkJdX8gjH2du9ztjdsmm58BVv6lL//1Hv6qthhx0SHW5CBUkOC8zsYiDNzNoB/YB3ww1LRCR8pQ8+O87gOYYwlCPXvQd5O3v3LFx5pVcttRYI0q10NdAe+AOYDqwDrgkzKBGRsEUi3nTVaHUopDOP8gkHMot/8TdWksn98O23cP31tSYxQLDZShuALP9PRKTaK+5OKlaXAroSYRDD2Zuv+JK96M6DPMzF7NKqHjRMXqzJUmpyMLNnKWNswTl3VigRiYiEKBKB7t294qgN2UhPJnEjI2jFcj7hADrzKE9wLkWklV5ErxYoq+UwKmFRiIgkQHGLYZvC9fRmHNdzBzvzE+9wJH0Yy2xOIz3dcGugVYaXGCq1olsNUGpycM69kchARETCNmLQGgZsuIf+jKEZvzCHk7iIGbzBPwCjVStvdqqU3a30qHPuAjP7nDjdS8656n1vuIjUGjPv+5FVQ0bz9rqxNOE3nuJshjKEjzh08zG1uQspnrK6lfr7j2cmIhARkSqXm8vrZ47k9AUTqU8+j3AhwxjMAv6+xWFpaXHWYajlSp3K6pz7wX+a6ZzLjf4DMhMTnohIJSxezNJje1DQeg+OXJBDhK7sxVd05eGYxNCokbdAmxLDloLc53BynH2nVXUgIiKVFYl4pS862Kc8YhdStNfetHzrEe4nk7Ys5QomsJQ9Ys6Lu3KbAGWPOfTBayHsbmafRb3UBHgn7MBERMoTiUD//rBH3vtMJpt/MYt1NOF2BnIn17KKHUs9V4PPZStrzOFhYDYwDBgUtX+9c25NqFEBZrYMWA8U4tV36hj2e4pI9ZHZx7F43Ks8wlBO5FXyaMZ/uIV76ctayq57ZKbB5/KUNZX1V+BXoIuZpQE7+cc3NrPGzrmAVc+3yvHOuTg1EkWktopMczyfOYur12dzOB+wkpZcxx3k0IvfKb+8hZlXZVtdSWUrt3yGmfXFK9v9E1Dk73aAprKKSOIUFjLxtMfpOGcoXfmMb2lNb8byIJfxR8D6Funp3tLOSgzlC1KV9RpgL+dcTDHbkDngJTNzwHjnXE6C319EUkF+Pu/1ncaOE4dzedHXLGJvuvEQ0+nCnwGXs1dSqLggyeE7vO6lRDvKObfSzHYE5pjZl865N4tfNLNeQC+AjBq40IZIrbdxI0ycyJrBIzjit+/4mAM5j8d5kk64ciZa1qkDRUXeoHNtLoGxNYIkh2+A183sObyy3QA450aHFpV3/ZX+489m9iTeanRvRr2eA+QAdOzYUYsPiVRzxTOP8vPW0ZtxXMdoduYnvuAoshnPC5wKWKnnF48l3H9/4mKuyYIkh+X+X33/L3Rmti1Qxzm33n9+CnBLIt5bRBIrEvHW0Gnwex79uJt+3M0OrOVFTuECsniLY8u9hrqNql6Q9Rz+l4hAStgJeNJfmrQu8LBz7oUkxCEiIcrMhCfH/sBNjKYPY2nM7zzJOQxlCHM5pNzz1VoIT5DZSi2AG/FWg9s8JcA5d0JYQTnnvgE6hHV9EUm+/7s0l/bTRjCaidSjgBlcxDAGs5D9Ap2v1kK4gnQrRYBH8Arw9Qa6A6vCDEpEaq5nRn7FuiHDuOnPCA7jQS5jBDfGLW9Rmj591FoIW5DkkO6cm2hm/f01Ht4wM631ICIVM38+884fyplLH2cTDbmPqxjFAL5n1wpdRokhMYIkhwL/8QczOwNYCRX8tykitde773rzSZ9/nj3YjuEM4i6uKbPuUTzqRkqsIMnhNjPbHrgeuAfYDrg21KhEpFqLTHPMzHyFq9dnczyvs5p07uQ27uMqfqVpmec2bgzjxikJJFuQ2Uqz/Ke/AseHG46IVGtFRbw+YBbt7srmCfch3/M3rmU0OfRiA9uWe7q6jFJHkNlKk4m/TGjPUCISkeqnsJC3+z3KDuOHcVzh53xDG3oxnofoTj4NAl1CiSG1BOlWmhX1vCHQCW/cQURqu/x8pv1zKoe9PpyjWcIX7MMlTGUGF1EY6OtF3UipKki30szobTObDrwcWkQikvJmTN7I/L4TyNwwkkv4jnkcxLnM5CnOKbfuUTS1FlJXsNS+pXaAKt2J1Ebr1jHhoPv519I7uYifeYuj6UUOL/JPyqp7FI8SQ2oLMuawHm/MwfzHH4GBIcclIinksXF5LO03hisL7uHfrOUF/kk2WbzNMRW+lqakVg9BupWaJCIQEUktkQjcdMVKem8cTW/G0Zjfmcm5DGMw86jYqr1qJVQ/ZSYHM9sG6Ars6++aCzzunMsPOzARSbzistmN85YxkNtZyCTSKGQ6XRjGYBZt/ioIpk4dr+KqEkP1U+rIkZn9HVgEHAMsA3KBfwLvmFlTM7stIRGKSKgyM73qpmZw6yVfckded5awBz2ZxINcxp4sphtTAyeGhg1h2jRwDgoLlRiqq7JaDncDVzjn5kTvNLOTgAXAwjADE5HwFK+h8Pvv3vYBfMIQhnIeM9lEQ+7hakYxgJXsEviaGkuoWcpKDi1LJgYA59zLZlaAd7+DiFQjkQj07An5fsfwkbxDFtmczmx+ZTuGMoQx9Gc1LQJf88QT4WVNbq9xypqQXMfMYm5tNLOGQIFzbkN4YYlIVYtE4NJLIT/fcRJzeI3jeIejOYSPGEI2GSznP9wWODE0bux1Hykx1ExlJYcpwEwza128w3/+KDA1zKBEpGpFItDtkiLOck/xAYcxh1PYgyVcw520ZhnDGMI6ti/3OsUJwTlYv15dSDVZqcnBOXcb8ALwppmtNrPVwBvAHOfcrYkKUEQqLxKBRvX/5LlLHmY+HXiKTqSTxxXk0JaljOGacgvimXlTUZUQapcyp7I65+4F7jWzJv72+oREJSJbJTMTJozNpxtT+Izh7MFSFrIvXZnGI1wYqO6RBphrt0DlM5QURKqPg/fZwNFfPsBSRrEbK/iIjpzDkzzDWWXWPVIBPIlWmdpKIpJiIhG4odevdN9wP7O5kx1ZxRscy+VMZA4nU1bdI802kniUHESquXOPXc2Bb43hC+6hKb8ym1PJJot3OLrcc1XWQkoTpPBeI7wlQjOcc1eYWTtgr6gV4kQkCQZ3X0mLKaOYyni2YSNPcC5DGcInHBTofCUGKUuQlsNkYB5whL+9AniMLRcBEpEEyMyE2WO/ZSC3czOTSaOQh7mY4QyqUHmLCRM0tiBlC7IqR1vn3AigAMA5t5GKFm4XkUqLRLzB4n3tCw4f242vaUcPJjOZHuzJYrozJVBiSE/37lHYuFGJQcoXpOWQ71dndQBm1hb4I9SoRGRzqYv2+R/zENl04kk2sg1j6M9orgtc92jffWGhKqFJBQVJDjfh3Qy3m5lFgKOAy8IMSqS2i0Rg3CVv8xTZnMYLrGV7ssliDP3Jo3mga5hB794aV5DKCbLYzxwz+xg4HK87qb9zbnXokYnUQpl9HEvHvUQW2bzFW/xMCwYzlPvJDFTeAtRSkKpRanIws5JTHn7wHzPMLMM593F4YYnUHpEI9O5VxEkbnmYIQzmEuXzHrvRjDBP4NxtpFPhamoEkVaWslsMdZbzmgBOqOBaRWqdv7z/5ZfwjvMcw9mMhS2jLv3mAKXSjgPqBrqHuIwlDqcnBOXd8IgMpycxOBcYAacAE59zwZMYjUlUyM2Hi2D/ozkMM5Hba8g0LaM/FRHiUCwLVPQK1EiRc5U5lNbOGZnadmT1hZjPN7Bp/TYfQmFkacB9wGt761V3MrGKL14qkmEgEtkv7nfpj72IpbcnhSvJI52yeYn8+YzoXl5sYopfgVGKQMAX5iTIFWA/c4293wVvPoXNYQQGHAkucc98AmNkM4GzgixDfU6RKRS/FuR2/chX3sZQ7acFqXucf9GAyL3MSQW8bUg0kSaQgyWEv51yHqO3XzOzTsALy7QJ8F7W9Ajgs+gAz6wX0AsjIyAg5HJHgopfibM4qBnMXfbmX7VnHc5zOUIbwLkcFvp7GFCQZgtwh/YmZHV68YWaHAe+EF5L3NnH2uS02nMtxznV0znVs0SL4erciVa34DmYz7++SS6B5/veM5lqW0ZrBDOMlTuFAPuZMngucGIoX2CkqUmKQxAvScjgM6GZmy/3tDGCRmX0OOOfc/iHEtQLYLWp7V2BlCO8jUmnRLYRiu7OUGxnBZTxIGoVM4xJuZyBfsk/g62qgWVJBkORwauhRxPoIaGdmbYDvgYuAi5MQh0hcJ50Er7zy1/a+LGQww+jCdAqox0QuZwQ3kkvrwNfUYjuSSoLcIZ1rZjvg/ZKvG7U/tJvgnHN/mllf4EW8qayTnHO651OSLhKB7t2hsNDbPoh5ZJHNuTzJb2zLnVzLHVzPj7QMdD0NMkuqCrKew614tZSW8le/f+g3wTnnngeeD/M9RILIzISxY7fcdwxvMoShnMqL/EJTbuE/jKE/a0gPdE21EiTVBelWugCvbHd+uUeK1CAlWwng+CcvkkU2x/A2P7EjAxnOWPqwnu3KvJbWUJDqJkhyWAA0BX4OORaRlBE9pmAUcQ5PMYShdGQey9mNq7mbiVxebt0jtRCkugqSHIbhTWddQNQ6Ds65s0KLSiTBom9YK5bGn3RhOoMZxr4s4mv2oCcTmcYlpdY9UgtBaoogyeEh4Hbgc6Ao3HBEEq/kzKP6/MFlPMhAbmd3vuUz/s5FTOcxOlNEWtxr1K0LDz6opCA1R5DksNo5d3fokYgkWMkxhUb8Ti9yGMAodmElH3Ao13AXszgTV8b9oppxJDVRkOQwz8yGAc+wZbeS1nOQaicz0xsDcFH322/PWvpyL9dwF83J41WOpxtTeJUTKKvukcYTpCYLkhwO9B8Pj9qn9RykWomdeeTVPbqWO7mK+9iedcziDLLJ4n2OKPU6aiVIbRHkJrikrusgsjXi3aOwCysYwCh6kUNDNvE45zOUIXzKAaVeR2MKUtsEWlXEzM4A2gOb13Fwzt0SVlAiWyteUmjLEgZyO915CMMxjUsYziAWs1ep11FFVKmtgtwhPQ5oBBwPTADOBz4MOS6RCos3HRWgPQsYzDAuYgYF1OMBrmAkN5Ra90itBJFgLYcjnXP7m9lnzrn/mdkdwBNhByYSVLxWAkBHPmIIQ+nEU6ynMXdwPXdybZl1jzSmIOIJsp7DRv9xg5n9DSgA2oQXkkgwkQg0aFAyMTiO5Q1e5BQ+4lD+wRvczE20IpeBjCg1MTRu7C2/qcQg4gnScphlZk2BkcDHeDOVHgg1KpFylLxxDRyn8gJZZHM07/AjO3EjtzOWPvxGk7jXUCtBpHTlthycc7c659Y652YCrYC9nXP/DT80kViZmd4gcXTdo/N4nHkczGxOJ4Pl9OUe2vAtI7kxbmJQK0GkfKUmBzM7xMx2jtruBjwK3GpmzRIRnAhsuQxncRdSXQq4lCkspD2P05nG/EYPJrEHS7iPvmximy2uUbeulxCcg/XrNdgsUp6yWg7jgXwAMzsWGA5MAX4FcsIPTWqreGsyF89AasAmrmQci9mTKXQnn/pcyAz2YREP0iOmIJ6Zt+xmQYESgkhFlDXmkOacW+M/vxDI8buWZprZ/PBDk9qotJlH2/IbVzKe67mDv/ED73MY/bibWZxJvBIXWodZZOuUmRzMrK5z7k/gRKBXwPNEAivt3oRiTfmFvtxLf8bQnDxe4QQuYRqvcTxKCiLhKetLfjrwhpmtxpvO+haAme2B17UkUmnlJYUW/Ly57tF2rOdZziSbLD7YosTXX7SOgkjVKjU5OOeyzewVoCXwknOb61jWAa5ORHBSM5XWdQSwK99xAyO5ggdowB88ygUMYzCf0SHu8UoKIuEos3vIOfd+nH2LwwtHaqp4pbKjtWUJgxhON6ZgOKZyKcMZxNfsGfd4lcsWCZfGDiR0ZbUU9uNzBjOMC3mEAuqRQy9GcgPLaRVzbJ06XleUxhREwqfkIKGKROInhkP4kCyyOZtnWE9jRjGAO7mWn9h5i+PS0uChh9RCEEm0ILWVRALLzPR+4Uffo/AXxz94nZc4mQ85jGN4i5u4mVbkMojbYxJDeroSg0iyqOUgVSISgZ49IT8/3quO03meIQzlKN7lB3ZmACMZz5VblLfQNFSR1KHkIJUSiUD//pCXV/oxdSjkXJ5gCEM5kPnkkkEm9zGJnvzx17pRGlwWSUHqVpIKy8z0uotKSwx1KaAbD7GQ9jzGBTRiA5cxmT1YwlgyNyeG9HSv3pFqHYmkHrUcJJAgLYUGbKInk7iREbQml/l0oDOP8gTnUkTa5uPq14dJk5QQRFKZkoOUKxKBHj284nXxbMtv9GYc13MHLfmRdzmCq7iP5zmdkiUu1IUkUj2kXLeSmd1sZt+b2Xz/7/Rkx1QbRSLQvPlfM47iJYYdWMN/uIVcWjGKG1hIe47nVY7iHZ7nDKITg7qQRKqXVG053OmcG5XsIGqbIF1HADvyE9cxmkzupwm/8TRnMZQhfMhhm49RC0GkekvV5CAJVtZdzMV2Yzk3MJJ/M4H65G+ue/Q5+28+Jj0dxoxRUhCp7lKuW8nX18w+M7NJZrZDsoOpyYq7j8pKDO1YzER6spS29GYcD3Mxe/MlFzOdlen7b15hzTlYvVqJQaQmSErLwcxehhK3w3qygLHArYDzH+8Aesa5Ri/8NSYyMjJCi7UmK68Y3t/5jCEMpTOPkU99xtKHUQzgOzKoUwemTVEiEKmpzJX2zZACzKw1MMs5t19Zx3Xs2NHNnTs3ITHVFJEIXHpp/MRwKB+QRTZn8SzraML9ZHIn1/IzOwGaiipSU5jZPOdcx3ivpVy3kpm1jNrsBCxIViw1VSQC3buXTAyO43mVOZzEBxzOUbzDf7iFVuQymOGbE0N6uhKDSG2QigPSI8zsALxupWXAlckNp2aJ7UpynMFzZJHNEbzPD+zM9YxiPFdijRtrxpFILZVyycE5d2myY6gpIhHIyoLcXO9+heiWQh0KOY+ZDGEoB/Apy2hFH+5nMj1onN6Q8ZpxJFKrpVy3klRcJAKtW3sJoG5d77FOHe/mtdxc75jixFCXAi5jMl+wL49yIQ3ZRHcepB1fM9760LNPQ804EpHUazlIcPFuWiss9B5LDjQ3ZOPmuketWM4nHMD5PMaTdKKINNLSYKrWThARn5JDNVNWV1E8jVm/ue7RzvzEOxxJH8Yym9MoLm9hpkV1RGRLSg7VRLxWQlmJYQfW0I+76cfdNOMXXuJkLiSLNzmW6JpHZtC7txKDiGxJyaEaiESgVy/YsKH8Y3fiR65jNH0YSxN+4ynOZihD+IhDY45VqQsRKY2SQ4op7jZavhyaNfP2lVcIDyCD3M11j+pRwCNcyDAGs4C/xxyrpCAi5VFySAGljSMESQp78hWDGM4lTMNhTKE7wxnIsrQ9KCyEVq0gO1uJQEQqRskhyUp2GQWtZrI/n26ue7SJhkxrkknT2wbw73678e/wwhWRWkLJIQmiu47q1Plr+mkQh/MeWWRzJs+x3pqw6MyBtJ9wLT123DG8gEWk1lFySLCSLYVgicFxAq+SRTYn8Bpr6qTz6bm30uGBvrRv2jTMcEWkltId0gmWlRVs1pHHcSbP8h5H8AonsY99ybyL76DZr8vo8Nj/gRKDiIREySHBli8v/5g0CrmQGXyedgDPchY78jNDmo3j9YnfcHDkOm8NThGRECk5hKy47lGdOt5j8fTUktLSoD75DEifxJqd92EGXdhvzwKYMoXdCxYzNO9KuvRomMjQRaQW05hDiEqOL+TmQr163mI5+fl/Hddsm4280Hkih7w2Ar77Dg48EO59HDp18rKKiEiC6ZsnRPHGFwoKoEkT7/6D7VjHsKa3s6Jeaw6ZcjVkZMDzz8O8eXDeeUoMIpI0ajmEqNTxhbw8ll11N9x9N6xdC6ec4mWSY49NaHwiIqXRT9MqVN74ws78wAhuINdawS23wPHHw4cfwosvKjGISEpRy6GKlDW+0DJ/GTcygp5Moh4FLD+iC21yBkP79skNWkSkFEoOVSTe+MLuBV9yU4PhnE8Eh/F448vY9n8DOfu6tskJUkQkICWHKhI9vtCB+WSRzXnMZNMfDanX/yoYMICLd901eQGKiFSAxhwqoeTYQiTiTTQ6gneZxRnM50BO4SWGMZijd82Fu+4CJQYRqUbUcqig2LEFx/TLX+Gl9Gz25HVWk04Wt3EfV1HQqCk5w5Mbr4hIZajlUEHFYwtGEWfxNO9zOLP+OJntflrMvK6jOWa3XIZZFk1bNSUnR+soiEj1pORQQrwuo2grcgu5iOl8Sgee5hxasIpejKd14TccPO1aFpOBQicAAAlrSURBVC3flqIiWLZMiUFEqi91K0WJNx21Vy/vedfO+TB1Kl/XHU6bP5ewkH25hKnM4CIKqUurVsmLW0Skqik5RIk3HdVt2MDXV0+AQSNhxQq2a3MwF33/BI/mn43zG16NGnlLcYqI1BQ1tlupvO6heKKnozZhHQMZzjJac/Mv/aFNG3jhBdKXfsS/JnUio1UdzLwaSRpbEJGaxlzQRYtTWMeOHd3cuXM3b5fsHgLv1315X+KtW8Nvuavpx91czT3swFpmcyqTdhrCYz8eE94HEBFJAjOb55zrGO+1GtlyiNc9tGGDt79UK1cye9/rWUZr/sutvMoJHMxczm80m3PuUGIQkdqlRiaH0qqhxt3/7bfQpw+0acM+L41h1VGdOKnlQjrbTPJaHawuIxGplZKSHMyss5ktNLMiM+tY4rXBZrbEzL4ys39W5voZGQH2L1oE3btDu3YwaRJcdhksXkybt6fy8sp9NR1VRGq1ZLUcFgDnAm9G7zSzfYGLgPbAqcD9ZpZW0YtnZ3tjDNE2zyj65BM4/3yvIurjj0O/fvDNNzB+POy+eyU/johIzZKUqazOuUUAZlbypbOBGc65P4BvzWwJcCjwXkWuX/xrPyvL60rKyICc7u9wSiQbZs+G7beHIUOgf39o0WJrP46ISI2TamMOuwDfRW2v8PfFMLNeZjbXzOauWrUq5vWuXWHZt46iF15iWat/cMotR8NHH8HQod7dbbfdpsQgIlKK0FoOZvYysHOcl7Kcc0+XdlqcfXHn2jrncoAc8KayxhyQlwenneYlhF128SqjXnFFbH+TiIjECC05OOdOqsRpK4DdorZ3BVZWKoBmzbwb1664Arp1gwYNKnUZEZHaKNXKZzwDPGxmo4G/Ae2ADyt1JTN45JEqDE1EpPZI1lTWTma2AjgCeM7MXgRwzi0EHgW+AF4ArnLOFSYjRhGR2ixZs5WeBJ4s5bVsQGXsRESSKNVmK4mISApQchARkRhKDiIiEkPJQUREYig5iIhIDCUHERGJUSNWgjOzVUBusuPYSs2B1ckOIkFq02eF2vV59Vmrl1bOubhF5mpEcqgJzGxuacv11TS16bNC7fq8+qw1h7qVREQkhpKDiIjEUHJIHTnJDiCBatNnhdr1efVZawiNOYiISAy1HEREJIaSg4iIxFBySCFmNtLMvjSzz8zsSTNrmuyYwmJmnc1soZkVmVmNnA5oZqea2VdmtsTMBiU7njCZ2SQz+9nMFiQ7lrCZ2W5m9pqZLfL/G+6f7JjCoOSQWuYA+znn9gcWA4OTHE+YFgDnAm8mO5AwmFkacB9wGrAv0MXM9k1uVKF6EDg12UEkyJ/A9c65fYDDgatq4r9bJYcU4px7yTn3p7/5Pt4a2jWSc26Rc+6rZMcRokOBJc65b5xz+cAM4OwkxxQa59ybwJpkx5EIzrkfnHMf+8/XA4uAXZIbVdVTckhdPYHZyQ5CKm0X4Luo7RXUwC+Q2s7MWgMHAh8kN5Kql5RlQmszM3sZ2DnOS1nOuaf9Y7Lwmq6RRMZW1YJ81hrM4uzTvPEaxMwaAzOBa5xz65IdT1VTckgw59xJZb1uZt2BM4ETXTW/CaW8z1rDrQB2i9reFViZpFikiplZPbzEEHHOPZHseMKgbqUUYmanAgOBs5xzG5Idj2yVj4B2ZtbGzOoDFwHPJDkmqQJmZsBEYJFzbnSy4wmLkkNquRdoAswxs/lmNi7ZAYXFzDqZ2QrgCOA5M3sx2TFVJX9iQV/gRbwBy0edcwuTG1V4zGw68B6wl5mtMLPLkx1TiI4CLgVO8P8/nW9mpyc7qKqm8hkiIhJDLQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOkjLMLD1qauCPZva9/3ytmX2R4FjOiS6mZma3mFmFb+ozs9alVSo1s/Zm9qqZLTazpWb2PzOr8v8ny/osZvZ6Ta2KK1tHyUFShnMuzzl3gHPuAGAccKf//ACgqKrfz8zKqhBwDl411eLY/uuce7kK33sbvJvihjvn9gT+jlesL4zyz6F+FqmZlBykukgzswf8+vkv+V+umFlbM3vBzOaZ2Vtmtre/v5WZveKvjfGKmWX4+x80s9Fm9hpwe7zzzexI4CxgpN9yaeufd75/jUPM7F0z+9TMPjSzJn4L4S0z+9j/O7Kcz3Mx8I5z7iUA/474vsAN/nvcbGYDig82swV+kTfM7Ck/3oVm1ivqmN/MLNuP630z26m8zxLNzE4xs/f8+B/zawdhZsPN7Av/n+WoCv+bk2pJyUGqi3bAfc659sBa4Dx/fw5wtXPuYGAAcL+//15gir82RgS4O+paewInOeeuj3e+c+5dvF/1N/gtmaXFJ/qlMB4B+jvnOgAnARuBn4GTnXMHAReWeL942gPzonf477ONlb/IU08/3o5APzNL9/dvC7zvx/UmcEVZnyWamTUH/s//53IQMBe4zsyaAZ2A9v4/y9vKiU1qCBXek+riW+fcfP/5PKC1/8v2SOAxr9wNAA38xyPwFhMCmAqMiLrWY865wnLOL81ewA/OuY8Aiqtxmtm2wL1mdgBQiJeAymLEr9Iar5prSf3MrJP/fDe8xJkH5AOz/P3zgJMDXKvY4XhdT+/4/yzq45XDWAdsAiaY2XNR15caTslBqos/op4XAtvgtXzX+uMS5Yn+Iv7df6zI+cVK+1K/FvgJ6OBfd1M511kIHLvFhc12B1Y759aa2Z9s2bJv6B9zHF5r5Qjn3AYze734NaAgqpJvIRX7/9uAOc65LjEvmB0KnIhXPLAvcEIFrivVlLqVpNryf7V/a2adwauWaWYd/JffxfsyA+gKvF3B89fjFUEs6Uvgb2Z2iH9OE39ge3u8FkURXlG2tHLCjwBHR80a2gavK+om//VlwEH+awcBbfz92wO/+Ilhb7xf/OUp7bNEex84ysz28N+zkZnt6beutnfOPQ9cgzc5QGoBJQep7roCl5vZp3i/xouX4uwH9DCzz/C+rEubBVTa+TOAG8zsEzNrW3ywv+TnhcA9/jlz8H653w90N7P38bqUfqcMzrmNeAPFWWa2GFiNN0BdvMDTTKCZmc0H+uCtKQ7wAlDX/1y34n2plyfuZykRzyrgMmC6f+33gb3xksosf98beC0kqQVUlVUkBZjZOcBo4HjnXG6y4xFRchARkRjqVhIRkRhKDiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIx/h/Ny7kFWvsEHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from statsmodels import graphics\n",
    "graphics.gofplots.qqplot(resid, line='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLM: Gamma for proportional count response\n",
    "\n",
    "### Load Scottish Parliament Voting data\n",
    "\n",
    " In the example above, we printed the ``NOTE`` attribute to learn about the\n",
    " Star98 dataset. statsmodels datasets ships with other useful information. For\n",
    " example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "This data is based on the example in Gill and describes the proportion of\n",
      "voters who voted Yes to grant the Scottish Parliament taxation powers.\n",
      "The data are divided into 32 council districts.  This example's explanatory\n",
      "variables include the amount of council tax collected in pounds sterling as\n",
      "of April 1997 per two adults before adjustments, the female percentage of\n",
      "total claims for unemployment benefits as of January, 1998, the standardized\n",
      "mortality rate (UK is 100), the percentage of labor force participation,\n",
      "regional GDP, the percentage of children aged 5 to 15, and an interaction term\n",
      "between female unemployment and the council tax.\n",
      "\n",
      "The original source files and variable information are included in\n",
      "/scotland/src/\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sm.datasets.scotland.DESCRLONG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Load the data and add a constant to the exogenous variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[7.12000e+02 2.10000e+01 1.05000e+02 8.24000e+01 1.35660e+04 1.23000e+01\n",
      "  1.49520e+04 1.00000e+00]\n",
      " [6.43000e+02 2.65000e+01 9.70000e+01 8.02000e+01 1.35660e+04 1.53000e+01\n",
      "  1.70395e+04 1.00000e+00]\n",
      " [6.79000e+02 2.83000e+01 1.13000e+02 8.63000e+01 9.61100e+03 1.39000e+01\n",
      "  1.92157e+04 1.00000e+00]\n",
      " [8.01000e+02 2.71000e+01 1.09000e+02 8.04000e+01 9.48300e+03 1.36000e+01\n",
      "  2.17071e+04 1.00000e+00]\n",
      " [7.53000e+02 2.20000e+01 1.15000e+02 6.47000e+01 9.26500e+03 1.46000e+01\n",
      "  1.65660e+04 1.00000e+00]]\n",
      "[60.3 52.3 53.4 57.  68.7]\n"
     ]
    }
   ],
   "source": [
    "data2 = sm.datasets.scotland.load()\n",
    "data2.exog = sm.add_constant(data2.exog, prepend=False)\n",
    "print(data2.exog[:5,:])\n",
    "print(data2.endog[:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Fit and summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                   32\n",
      "Model:                            GLM   Df Residuals:                       24\n",
      "Model Family:                   Gamma   Df Model:                            7\n",
      "Link Function:          inverse_power   Scale:                       0.0035843\n",
      "Method:                          IRLS   Log-Likelihood:                -83.017\n",
      "Date:                Fri, 24 Apr 2020   Deviance:                     0.087389\n",
      "Time:                        14:17:04   Pearson chi2:                   0.0860\n",
      "No. Iterations:                     6                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1          4.962e-05   1.62e-05      3.060      0.002    1.78e-05    8.14e-05\n",
      "x2             0.0020      0.001      3.824      0.000       0.001       0.003\n",
      "x3         -7.181e-05   2.71e-05     -2.648      0.008      -0.000   -1.87e-05\n",
      "x4             0.0001   4.06e-05      2.757      0.006    3.23e-05       0.000\n",
      "x5         -1.468e-07   1.24e-07     -1.187      0.235   -3.89e-07    9.56e-08\n",
      "x6            -0.0005      0.000     -2.159      0.031      -0.001   -4.78e-05\n",
      "x7         -2.427e-06   7.46e-07     -3.253      0.001   -3.89e-06   -9.65e-07\n",
      "const         -0.0178      0.011     -1.548      0.122      -0.040       0.005\n",
      "==============================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/genmod/generalized_linear_model.py:274: DomainWarning: The inverse_power link function does not respect the domain of the Gamma family.\n",
      "  warnings.warn((\"The %s link function does not respect the domain \"\n"
     ]
    }
   ],
   "source": [
    "glm_gamma = sm.GLM(data2.endog, data2.exog, family=sm.families.Gamma())\n",
    "glm_results = glm_gamma.fit()\n",
    "print(glm_results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLM: Gaussian distribution with a noncanonical link\n",
    "\n",
    "### Artificial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "nobs2 = 100\n",
    "x = np.arange(nobs2)\n",
    "np.random.seed(54321)\n",
    "X = np.column_stack((x,x**2))\n",
    "X = sm.add_constant(X, prepend=False)\n",
    "lny = np.exp(-(.03*x + .0001*x**2 - 1.0)) + .001 * np.random.rand(nobs2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit and summary (artificial data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                  100\n",
      "Model:                            GLM   Df Residuals:                       97\n",
      "Model Family:                Gaussian   Df Model:                            2\n",
      "Link Function:                    log   Scale:                      1.0531e-07\n",
      "Method:                          IRLS   Log-Likelihood:                 662.92\n",
      "Date:                Fri, 24 Apr 2020   Deviance:                   1.0215e-05\n",
      "Time:                        14:17:04   Pearson chi2:                 1.02e-05\n",
      "No. Iterations:                     7                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1            -0.0300    5.6e-06  -5361.316      0.000      -0.030      -0.030\n",
      "x2         -9.939e-05   1.05e-07   -951.091      0.000   -9.96e-05   -9.92e-05\n",
      "const          1.0003   5.39e-05   1.86e+04      0.000       1.000       1.000\n",
      "==============================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-21-254e2844f190>:1: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n",
      "Use an instance of a link class instead.\n",
      "  gauss_log = sm.GLM(lny, X, family=sm.families.Gaussian(sm.families.links.log))\n"
     ]
    }
   ],
   "source": [
    "gauss_log = sm.GLM(lny, X, family=sm.families.Gaussian(sm.families.links.log))\n",
    "gauss_log_results = gauss_log.fit()\n",
    "print(gauss_log_results.summary())"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
