{
 "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:                Mon, 24 Feb 2020   Deviance:                       4078.8\n",
      "Time:                        22:49:06   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": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZgU1dW43zNNA8M6IKAwOoALIAQBQVyIcYmKiguuuCR+Wf2SmHyfJiHi74uKxsQxaDSJGmISTTRGcQuKqLjgFhUVBFQ2kZ0BRYQBZWZglvP7o7uHnp6q6qrepnvmvM8zD91Vt26dqmHuufecc88RVcUwDMMwilpaAMMwDCM/MIVgGIZhAKYQDMMwjCimEAzDMAzAFIJhGIYRxRSCYRiGAZhCMAxEZICIqIi089H2WyLynxTucamIPJ+ahM36Ol5ENmaiL8OIxxSCUVCIyFoR2SMivRKOL4oO6gNaRrImiuXLuJ/FAKr6oKqeEtdWReRgj76+JSL10T52Rp/vjBRk+ruI3JTaExltDVMIRiGyBrg49kVEhgPFLSdOM0pUtUv0Z0Qa/bylql2AEuBvwCMi0jMzIhpGc0whGIXIA8Blcd//C7g/voGIdBeR+0XkMxFZJyK/FJGi6LmQiNwqIltFZDUwweHav4nIZhGpEJGbRCSUjsDxpiYReS16eHF0BTDJ61pVbQDuJaL0DnTo+1AReUVEKkVkiYicFT1+OXAp8IvofWal8wxG68cUglGIzAO6RQfCEDAJ+GdCmz8C3YkMoMcRUSDfjp77PnAGMAoYA5yfcO0/gDrg4GibU4DvZUp4Vf1a9OOI6Cpihlf7qG/je8CXwMqEc2FgFvA80Af4CfCgiAxW1XuAB4HfRu9zZqaewWidmEIwCpXYKuFkYDlQETsRpySuUdUvVHUtcBvwzWiTC4E7VHWDqm4Dbo67dl/gNOBKVd2lqluA24GLAsi2NTpbrxSRn6f8hHCUiFQCnxAxkZ2jqjsS2wBdgHJV3aOqc4GniTOpGYZfkkZVGEae8gDwGjCQBHMR0AtoD6yLO7YOKI1+7gdsSDgXoz8QBjaLSOxYUUL7ZPRS1boA7d2Yp6pfTdKmH7AhalaKEf+shuEbUwhGQaKq60RkDXA68N2E01uBWiKD+9LosTL2riI2AwfEtS+L+7wB2E3mBvVsswk4QESK4pRCGfBR9LOlMzZ8YyYjo5D5LnCiqu6KP6iq9cAjwK9FpKuI9Ad+yl4/wyPA/4jI/iLSA5gSd+1mIvb420Skm4gUichBInJchmX/FAcHcQq8Dewi4jgOi8jxwJnAwxm+j9EGMIVgFCyqukpV57uc/gmRgXI18B/gX0QidQD+AswBFgPvAU8kXHsZEZPTUmA78BjQN6PCw1TgH1E/w4WpdqKqe4CziPg9tgJ3A5ep6vJok78BQ6P3mZmmzEYrR6xAjmEYhgG2QjAMwzCimEIwDMMwAFMIhmEYRhRTCIZhGAZQgPsQevXqpQMGDGhpMQzDMAqKBQsWbFXV3l5tCk4hDBgwgPnz3SINDcMwDCdEZF2yNmYyMgzDMABTCIZhGEYUUwiGYRgGYArBMAzDiGIKwTAMwwBMIRiGYRhRTCEYhmEYgCkEwzAMI0rBbUwzDMNojcxcWMG0OSvYVFlNv5JiJo8fzMRRua2EagrBMAyjhZm5sIJrnviA6tp6ACoqq7nmiQ8AcqoUzGRkGIbRwkybs6JRGcSorq1n2pwVOZXDFIJhGEYLs6myOtDxbGEKwTAMo4XpV1Ic6Hi2MIVgGIbRwkweP5jicKjJseJwiMnjB+dUDnMqG4ZhtDAxx7FFGRmGYRhMHFWacwWQiJmMDMMwWjlf1NT6amcrBMMwjDxi5sIKbpi1hO1VkUG8pDjM1LOGpbR6qG9QHluwwXf4qikEwzCMPGHmwgomP7aY2nptPFZZXcvkRxcDwTapvbNmGzfMWsKSTTsZ3b8HC3xcYwrBMIxWR7bTQDj1D+k7hafNWdFEGcSobVCmzVnhq78N26oof3Y5sz/YTN/uHfnDxaM487C+FP0o+f1NIRiG0arIdhoIp/4nP7YYNDJwp3NPr41oyTap7dpdx/RXV/Hn11ZTJHDlSYfw3187iOL2Ic/r4jGFYBhGq8IrDUQmFIJT/06z+lTu2a+kmAqXgd9tk1pDgzJzUQW3PLecT3fu5uyR/bj61CEpbWozhWAYRqsi22kggvQT9J6Txw9u5kMACBeJ4ya199Zv54ZZS1m8oZLD9u/O3Zcezuj+PQPdMx5TCIZhtCrcZtmZSgPhNYtP956x1USyKKPNO6q55dnlzFy0iT5dO3DrBSM4d1QpRUUS6H6JmEIwDKNVMXn84CY2fshsGgin/sMhaeJDSOeeXhvUamrruee11fzplVXUq3LFCQfxo+MPpnOHzAzlphAMw2hVZDsNhFv/2bynqvL0+5spf3Y5FZXVnD58P6457VAO6NkpI/3HENXmzpB8ZsyYMTp//vyWFsMwDCMnfLBxBzc+vYR3127n0L7duP7MoRx14D6B+xGRBao6xquNrRAMwzDykC1f1HDrnBU8umAjPTu15+Zzh3PhmAMIpekn8MIUgmEYRh6xu66ee/+zlrte/pjddfV8/9gD+fGJB9OtYzjr9zaFYBiGkQeoKnOWfMpvnlnG+m1VnHRoH/5vwlAG9uqcMxlMIRiGYSSQ6dQXyfpbtnknN85aylurP2fQvl144LtjOfaQ3pl4lECYQjAMw4gj06kvvPo79pBe/O6Fj3jonfV0Kw5z49nDuGRsGe1CLVOZwBSCYRhGHJlOfeHW3/VPLqEBpWpPPZcdPYArTzqEkk7t05I9XUwhGIZhxJHp1Bdu1+2oqeVrg3pz7YRDOWTfrin1nWmsYpphGEYcbukmUk194XbdPp3b849vH5E3ygBMIRiGUeDMXFjBuPK5DJwym3Hlc5m5sMLXOTcmjx9Mcbhpyuh0Ul/86PiDmu0d6NiuiGvPGIpI9vYUpIKZjAzDKFi8HLZASs7hTKW+qKtv4F/vrOd3L3xEgyqd2oeo2lNPaRYK9mSKrKauEJFTgd8DIeCvqlqecL478E+gjIhyulVV7/Pq01JXGIYRY1z5XMfMo6VRM43buTemnJhVuf6zcis3Pr2Ejz79kqMP3IfrzhzKoX27ZfWeyWjR1BUiEgLuAk4GNgLvishTqro0rtkVwFJVPVNEegMrRORBVd2TLbkMw2g9pOIAzlRdBCfWbt3FTbOX8eKyTynr2Ynp3xjN+GH75p1pyI1smozGAh+r6moAEXkYOBuIVwgKdJXI2+oCbAPqsiiTYRh5RjqbwJLVPshmXYR4dtbUcufcj7nvjTW0DxXxi1MH851xA+kY9l++Mh/IpkIoBTbEfd8IHJnQ5k7gKWAT0BWYpKoNiR2JyOXA5QBlZWVZEdYwjNyT7iawZLUPslkXAaC+QXl0/gZufX4Fn+/aw/mH78/k8YPp061jxu6RS7KpEJzWSIkOi/HAIuBE4CDgBRF5XVV3NrlI9R7gHoj4ELIgq2EYLUC6m8D8OICnzVlBRWU1IZHGvuOvTZW3V3/ODbOWsnTzTsb078F93xrL8P27p9VnS5NNhbAROCDu+/5EVgLxfBso14hn+2MRWQMMAd7JolyGYeQJ6WwCSzQ13T5pZLNBPvY9k6koNmyr4uZnl/HMB5/Qr3tH/njxKM44rG/B+Am8yKZCeBc4REQGAhXARcAlCW3WA18HXheRfYHBwOosymQYRh6Rav3jIKamTKWi2LW7jrtf+Zi/vL6GIoGrThrE5V87kOL2heUn8CJrCkFV60Tkx8AcImGn96rqEhH5QfT8dOBXwN9F5AMiJqarVXVrtmQyDCOzpJsVNJX6xzMXVvCzRxZTnxAy7zbIp5uKoqFB+ffCCm55bjlbvtjNxJH9uPq0IfTtnnnndEuT1Y1pqvoM8EzCselxnzcBp2RTBsMwskMmsoIG3QQWu2eiMojhNMinugoBWLBuOzc+vZTFGyoZcUAJ0785msPLeiS9rlCxncqGYaREEFOM00oi1semymq6F4cp6RRmU2W1p9PX6Z7xOA3yqaxCNu+opvzZ5Ty5aBN9unbgtgtGcM6oUoqyWL4yHzCFYBhGSvg1xTitJCY/thgUahsiM/3K6trG9l4rDS8zj9sgH2QVUr2nnnteW82fXv2YBoUfn3AwPzz+IDp3aBtDZdt4SsNopWS6slcQ/JpinGb1tfXe0eNuKw23e4ZEuPnc4UwcVer6TuLPXTVjEdPmrGg8p6rMen8z5c8sY9OOGiYM78uU04ZwQM9Ofl9Hq8AUgmEUKJmu7BUUL1NM/KCc6sYhp9WA2z3jlUHQZHfrt1Xx2kefMX/ddob27cbtk0Zy5IH7pCh1YWMKwTAKlExX9orhd9XhZoqB5juEU8HJH5DM/OP1TmKfE8/97oWP6NWlPeXnDueCMQc0S1XdljCFYBgFSqYre0HwVUfMFBPPuPK5aSsDL6ev0z1jpPpOfnbyYC4aa2lxrECOYRQoma7sBcln2H7wGnyFiL3fiSKJnC8tKW40AQXF65307e6eX+jOlz8OfK/WiCkEwyhQMl3ZC9JPJTGufK6rz6C0pJg15RNocNlDoApryifwxpQTUzZ5ub2TS8aW0ckjUiibKbELCTMZGUaBkqnKXvFkKpVEIvGKKp2NYslIfCf7duvIgb07c9sLK+heHKY4XER1bbOEyllJiV2ImEIwjALGy56eCqls4gLvDWOJJSNTvYdfJo4q5fThfbn/rbX8/qWVvL1mG/91zACu/PogXl6xJespsQsZUwiGYTSS6qrDzeQi0KxcZTZWNjFUlbnLt/Dr2ctYvXUXxw3qzbVnHMrBfbpm/d6tAVMIhmE0IZVVR1AzUKZXNgArP/2CX81exmsffcaBvTtz37eO4IQhfZi5sIL/uvfdJgog2zWVCxVTCIZhpE22zUBeVFbt4Y4XV/LAvHV0ah/i2jOG8s2j+tO+XVGLb94rNEwhGIaRNi1hiqmrb+DBt9dz+4sfsbO6lovHlvHTkwexT5cOjW2ytXmvtWIKwTDaOJnKh5QNM5Abr6/8jBtnLWXlli855qB9uPaMoRzat1uzdtnYvNeaMYVgGG2YQjOprNm6i1/PXsqLy7ZQ1rMTf/7maE4Zuq9r+cpshri2RkwhGEYbpqVNKn5XJztravnjSyv5+5traR8qYsppQ/j2uAF0aOddvrIlfRuFiCkEw2jDtKRJxc/qpL5BeWT+Bm6ds4JtVXu4YPT+/Hz8YPp0dU9DEY+FmQbDFIJhFCiZsP23pEnFbXXys0cWc9WMRezTpT3t2xWxqbKGIwb04O9njGX4/t0D3yeXvo1CxxSCYRQgmbL9Z8uk4kdZOSkioLFe8tYv9yDAZUf354azhrn6CYzMYQrBMAqQZFlJ/a4ckplU4gf2kk5hVGFHda1jv7G2FZXVCDQmuXNTViGRxsHfDQUenLeew8t62Cw/BwRSCCJSBHRR1Z1ZkscwDB+42fhjg2+QlYObSSVxFbK9yr3ucWLbxGHeyVGdTBnEt8vnyKfWRNL01yLyLxHpJiKdgaXAChGZnH3RDKPtEUshPXDKbMaVz2XmwgrHdl42/nTrGcTwSliX2G+yttBciZUG8FOk+gxGMPysEIaq6k4RuRR4BrgaWABMy6pkhtHGCOIXcLL9e+G2ovjlzA946O0NTWbrseykfiKNYm38tI1XYpsqq9mvW0dXP4LXvYzs4adATlhEwsBE4ElVraX5itAwjDQJUq1s4qhSbj53OKUlxZ5VyGI4rSh+OfMD/jlvfTPTTUwRlXQKJ5U51m+yqKSYo7p6Tz13vPgRJ972Ch9u2sEpQ/elb/eOjZXSvnFUmeuz2Gay7ONnhfBnYC2wGHhNRPoD5kMwjAwTdE9AvO1/4JTZrv26RQ099PYG12uqa+vp0K6I4nDIV9EbpxVLzLFcWlLMz08ZhAh8/bZX2LSjhgnD+zLltCEc0LNTs37H9O9pm8laiKQKQVX/APwh7tA6ETkheyIZRtsknT0BbtcCnDfa2WmczKm7o7qW2yeN9BVl5BWt9P7GSm6YtZQF67YzrF83bp80kiMP3Mf1vraZrOVIqhBEZF/gN0A/VT1NRIYCRwN/y7ZwhtFSZCrhWxDS2RMwefxgrpqxyNGW+/LyzxyvSRb22a+kONCmrsS2W3bW8PNHF/PYgo306tKeW84bzvmjDyBU1NQk5PauTQHkHj8+hL8Dc4B+0e8fAVdmSyDDaGlizt2KymqUvTZ1t4ifTJHoFygtKebmc4f7Ghgnjip1dey5mZwuPvIA1/7SMdHU1NZz18sfc8Ktr/Dkogr++7gDefnnxzPpiDJHZdAS79pwxo8PoZeqPiIi1wCoap2I+AttMIwCpCUTvqUzMy4NaHK6aeJwgGZRRiGRJs5sL3niZ/d9u3dk/LD9eHH5p2zYVs3JQ/fl/04/lAG9Orte39LJ9Yym+Fkh7BKRfYhGFonIUcCOrEplGC1IoebQnzx+MMXhptk/k830b5o4nFU3n87a8gncMWkkxeFQo3JINltPnN1v2lHDfW+upa5eefB7R/KXy8Z4KgMo3HfdWvGjEH4KPAUcJCJvAPcDP8mqVIbRgrjNqPM97DEdkxMEC3t1aw+R6KJxB/fydc9CfdetFT9RRu+JyHHAYCK/6xXRvQiG0SrJtxz6QRzc6ZicgszW99Q1uEY1bdpRw7jyub7kzbd33dbxE2V0WcKhw0UEVb0/SzIZRouST2GPuaxo5ifsVVWZu3wLN81e5tqPsDeTaTJ58+ldGyCaJBZZRP4Y97Uj8HXgPVU9P5uCuTFmzBidP39+S9zaMHLOuPK5joN0aUkxb0w5MaP3mrmwgsmPLqa2Ye+YEC4Spl0wgomjSln56Rfc+PRSXl+5lQN7d+akQ/flgbfWOW5Gy4W8RjBEZIGqjvFq48dk1MRfICLdgQfSlM0wDB/k3OmamDVCYNfuOq5/8kP++fZ6OrcPcd0ZQ/nm0f0Jh4oY2rdbk9m9qxnJnMQFQSr1EKqAQ/w0FJFTgd8DIeCvqlru0OZ44A4gDGxV1eNSkMkwWiVBdi+nu5lu2pwV1NY3nd/X1iu/fPJDBLjkyDJ+evJgenZu33g+0WfhtqLxchK3xCbAfCQf3oMfH8Is9q4Ci4ChwCM+rgsBdwEnAxuBd0XkKVVdGtemBLgbOFVV14tIn+CPYBitF79O10z4Gtxm8arw7JXHMmS/bhmTN5Nytwby5T34WSHcGve5Dlinqht9XDcW+FhVVwOIyMPA2URqKsS4BHhCVdcDqOoWX1IbRgGRzszPr9M1Exu8XFcj3Ts2KoNkzxLUSZysrnJbWTHkywY9Pz6EV1PsuxSIT6e4ETgyoc0gIum1XwG6Ar93il4SkcuBywHKyspSFMcwck8mZn5+QknT9TXsqK7lkD5dmimE4nCIX5w6BGjudK6orGbyo4sbZQwibzL5EjfHJd6jtZEvG/RcN6aJyBcistPh5wsR8ZP+2impeWIAQjtgNDABGA9cKyKDml2keo+qjlHVMb179/Zxa8PID9xmflfOWORZES0oqW7wqm9Q/vX2ek689RVeXfkZRw3syX7dOjpubJv61JImEUgAtQ3K1KeWAP6rvQWRD9pGtbR82aDnukJQ1a5p9r0RiM+etT+wyaHNVlXdRSRFxmvACCIJ9Awja2TagRdfYD6WRdQtt1CMispqrpqxiPnrtjXmFUqVVDZ4vbXqc258einLNu9k7ICe/OPMoXyltLtr+8pq5/2oldW1Ka+E/FZ+a+1RSvmyQc93lFHU4dsx9j1m9/fgXeAQERkIVAAXEfEZxPMkcKeItAPaEzEp3e5XJsNIhUw78BL7izd3uMXlx1DgwXnrGdO/Z1oKKYjtfv3nVfzmmWU8t+QTSkuKufOSUUwY3hdJUnXNi1Rt4IlyF7mk5G7tqSzyZYOenyijs4DbiKS/3gL0B5YBw7yui2ZF/TGR1Nkh4F5VXSIiP4ien66qy0TkOeB9oIFIaOqH6TyQYSQj6OCVbDXhVWBecd+sFd/mhllLmt175sIKpj61pHFm3qNTmOvPHJaSjJG+PqSyuq7x2Olf2Y/fTRpJx4SEeG706BRme1XzVUKPTuG0bODxPodE5QptJ5VFPtSA8JPc7lfAUcBHqjqQyE7lN/x0rqrPqOogVT1IVX8dPTZdVafHtZmmqkNV9SuqekcKz2AYgQgyePnJ159s0IuVkfRie1Vtkz5jDtx4M832qlomP7a4mW0+mYxPvLeRnz2yqIkyAHhh6ac89+EnnnLFc/2ZwwiHmq4iwiHh+jOHZcwGnm6CPiM9/CiEWlX9HCgSkSJVfRkYmWW5DCNrBBm8/GQATTboxdI23DFppGOkRfy94j8nOnAhslEs0SHtJeP8tduY/Nj71DssUWob1NFZ6+YcnjiqlGnnj2gyWE87P5LWIpXU225MHFXKG1NOZE35BN6YcqIpgxzix4dQKSJdgNeAB0VkC5H9CIZRkARx4PlZTXg5RuP7nTiqlPnrtvHPec7ut/g+k6064v0ebm0rKqs5f/pbSfuJJ5l/xc2skS82cCM9/CS36wxUE1lNXAp0Bx6MrhpyjiW3a93kavu+3/v4TS7nFmUU6/eXMz9oVpkskU7hInp07uDpXHWSo2pPnaNtH5L7LwS4fdLIxmd3e96S4jCLrj8lqTxG/uInuZ0fhXAV8KjP3clZxxRC68XNodiSNuRMyPTLmR+4rgrygXjlNnDKbFcFckec4jAKDz8KwY8PoRswR0ReF5ErRGTfzIhnGE0JWrErF2TCyfnQ2xuSN3IhWSRoKI1Q0RjxJicvf0hr3xxm+EtdcQNwg4gcBkwCXhWRjap6UtalM9oU+bJ9P5F0wwH9mH5cUVhbPsFxpZJ231HilcDk8YO5csYix3bp/B7yIZOnkRw/K4QYW4BPgM8By0pqZJx82b6fadKZxceefeKoUs4bXeoZpZQKic70iaNK6dEp7ClLUPyE7hr5QVKFICI/jCafewnoBXxfVQ/LtmBG2yOToYvJSCXvTqocdWCPlK4rDoc4YUjvRjn/NW+9p4M4KG7mr+vPHJbR30M+mgINZ/yEnfYHrlRV53WkYWSIXIUu5jL3/MyFFbyzdrvjueJwEdW1DY7nSkuKOWFIbx5fUNEoZyaVgYBrSctM/x7y1RRoNMePD2FKLgQxDMjN9v1M556Pt4+XdArzZU0tLuN8EyIpI8Q1gmlc+dykSd9SpUiEmQsrXJ83k7+HIFXfjJYliA/BMFoFmZqxzlxYwcgbnufKGYsa7ePbq/wpA4DKqtomEUwlxWE6hou4asYijrn5Jc9MqX4QInmGwkXNPQ/1qjmz4+fSFGikhykEo82RCed1zOzklhLarxyxNA23TxrJ7roGtlfVosCmHTUp9wsRk9Oa8gksvO4Upl0wwtGxnSs7vuUnKhx8p782jNbCCUN6O24U27Zrt6cZJR6vDKd+CBdJkxlykP6KBLoXh6msqo2aqOqa5D1yihy6KguhpEHIh0yeRnJcFYKIfIGHH0tVk1fcNow85OXlnzker65t8O1cTncgnXbBiMZ7bN+1J5B5qEGhU/t2LLwukkrCT4y/2fENPyStmCYiNxLZf/AAEbPkpUTqHxtGQeI1mFfX1nPDrCUpD7B+iJlvausb+Oe8ddzx4srAfcQ/Q2I9gWlzVjQrUJ8vFbmM/MaPyWi8qh4Z9/1PIvI28NssyWQYWSXZYL69qrYxWZxbSKpbhtPO7UPsqav3dCzXq/KLx97n5meX8enO3Xz14F6MO2gf/jD342YDdsdwkWPiOqeZvZ9wWtstbHjhJ7ndm8BdwMNETEgXA1eo6jHZF685ltzO8MKP+cQtDYRf2oeEzh3asb2q1jGzaWIYamXUUZxIqEiY/o3RnHRoHyQaBpooO+Aoa3z1tPhMq04kZmY12iZ+ktv5WSFcAvw++qNEqqUl1kY2jBbH74az2OcbZi1xTRvtxZ56ZU/0unrVRtNLbHCOl8Gr//oG5eShe3NFejle40tpxvq95okPmL9uW5PNa0748XdYriEDfISdqupaVT1bVXupam9Vnaiqa3Mgm2EEIkiKhImjSll43SncMWlkk3DIkmLnPD5exN8jSLRQsrKa8bJ27tB87lZdW89Db29Ier8Sl9xEMTKZayiXKUGMzOMnl9EgEXlJRD6Mfj9MRH6ZfdEMIxipbDhLLNc49azmeXyC3DtI9FEQh65bv36ynX5ZU+c5MGcq15AlsSt8/GxM+wtwDVALoKrvAxdlUyjDSGWmGWTDmVfd4PhNVH6J3aO7zxVGSXE4kEnG7dn8ZFJ1q50cw833EDS01pLYFT5+FEInVX0n4ZjVVDayRqozTb8pEpL1H79q+MZRZb5kPmFIbx5+Zz07a5L7JIrDIaaeNcxXvzHcnu3iIw/wtaJxG9xnLqxwVXxB9yhYErvCx49TeauIHER0k5qInA9szqpURqvEr+MySPK5xD7PG13Ky8s/87yHn/7j++3cPsSuPd52+qcXb+Zfb6+nwcGC07l9iJJO7dNy2HqFjY7p37PxuFstZrfBfdqcFY4RUEIwk1bsHrb5rbDxoxCuAO4BhohIBbCGyOY0w/BNkJTTfmeaTn0+OG89xxzUs7F9zFwRP9AnM5H8cuYHPBhXe2DXnvrGLKRXzVjkOIB65TTataeeqj3pR+84RSElKsTElNngvQHN7V0rwVOB2+a3wsePQlinqieJSGegSFW/yLZQRusjyKzf70zTqU8F3li1rfF7vOIBmnx26n/mwoomyiBRVjfZBO96BfGmKchM3QUnhfj4ggpfq6QYbs/jNwIqHtv8Vvj4UQhrROQ5YAYwN8vyGK2UIPZlvzNNv7bp2GBetafONURTiPgB3EwosfuVnzec//v3h9Ql2Ib8Fq/xa/ryM5C6KdmXl3/WZCNazIHu1HemZ/WWxK6w8eNUHgy8SMR0tEZE7hSRr2ZXLKO1ESQCyG+65CC26YrKas9NYgo8vsDdnASReP7bX1hJXYM2OnJTqZfsZvoK6kR3kzX+uB8HuqWmNmL4qZhWDTwCPCIiPYjsWH4VCB6sbbR63Ga6QWeifmaak8cPdrXpJxJycRYS5hEAACAASURBVLbGk2yD1/aqWg7o0Y5Hf3A0FdurPdNFQCS9hJMSUmBc+dzGd5NqBTe3Z4pXUslCQeN/V7dPGmmKoI3jqx6CiBwHTAJOA94FLsymUEZhkuvkahNHlTJ/3TZHm388xeFQxkpRbthezbfufYc9dQ1NahA4saOqlnBIqK1v3i7+3aQarumm4OKPu/URu38u6kobhUNShSAia4BFRFYJk1V1V9alMgqSZDPdbNiXb5o4vEnYZSzSZvb7mxtn5x3aFdGhXVFa1c3iSRaCGqMBaKhX15l87N10Lw47ypbMJFbqwyHs5jQOiaS0KjFaN54+BBEJAfep6jmq+pApA8OLltqYlJh+Ykz/ntTE5Z+urK5l1546x9rCucDLVFVRWc2uPc33eSZWVHPCz0Y8tzZuMtkmsraNp0JQ1XrghBzJYhQ46dYqzlRiNKeVSm290qVjuybO028cVZaSUziThMTZpNSlY7ukM3U/DmG3Nm5hpbaJrG3jx4fwpojcSSTstHGFoKrvZU0qoyBJJ4QxyMY1p2vjTUZujt7KqloWXndKY/sH562nW3E7dlbX+Q4bzSRevo1Kn2m5/Zjh3NrYJjIjET8KIVYI58a4YwpYxQ2jCek4jlONtJm5sILJjy5udPB6Rf3ENp7FD4Q7qiPmmg7tithd1+ArGilVenQK06l9uybvxi1SKdszddtEZjiRtGJavmEV01onA6bMdjwuwJryCa7XjbzheV/O4ljqid8+t5xNO2qanY9VFXOqplYcDnHe6FIeenuDb2WRuHM5dn+gmQN8xjsbmkQshYuESWMP8L3b2DD84Kdimp96CPuKyN9E5Nno96Ei8t1MCWnkP9kuepJOxk0vZZBoNz+0bzdHZQB7naluNvebJg6nwacyKA6HuPSosmZ9AM02ic14ZwOJ5ZcbgBnvbLC6AkbO8WMy+jtwH/B/0e8fEfEn/C3ZhSJyKpGNbCHgr6pa7tLuCGAeMElVH/Mhk5Ejgtj2Uy3DOPWpJSll3PQ7QNY3KI8t2Mibq7YiAk7jekmnsGt6h/g2ThvNnExBTs89rnxuc2e3w16G+gYl0bNgIaFGLvCjEHqp6iMicg2AqtaJSNJA7GjI6l3AycBG4F0ReUpVlzq0uwWYE1h6I+v4te2n6hSeubDCdZbvlXEzdj8vYrb5T3bW8MnOGr56cC9O+8p+3DR7WZNnCoeEL2vqGgd7J9lnLqzgyxqH8NCQNBa7T0a6IZ0WEmpkGz8KYZeI7MPeeghHATt8XDcW+FhVV0evexg4G1ia0O4nwOPAEX6FNnKH107XmQsrmjgnU3EKe1XT8sq4GaR2cYw1W3dx6VH96dyhXZOVzK7ddc2UUqLs0+ascJzNd27vHB4aWy1VVFZnzFFtIaFGtvGjEH4KPAUcJCJvAL2B831cVwpsiPu+ETgyvoGIlALnEIlYclUIInI5cDlAWZm/ClZGZvAK44yfRae6Kc3rvJe5KJXZcryfIH4Qd3Noxz+32/12OKxuEldLXsqgCAi5pLeIx0JCjVyQ1Kkc3W9wHJHw0/8GhkXrKifDyU+Y+L/+DuDq6AY4LxnuUdUxqjqmd+/ePm5tZAqnna4x4pOkpbopze18j07eNYdTmS0HrUscfzzI8wVZvXTvFGba+SOS1h+wDKRGLvATZXQBUKyqS4CJwAwROdxH3xuBA+K+7w9sSmgzBnhYRNYSWXXcLSIT/Qhu5IaJo0o5b7T7QBSbOfutZ5yI23XXn+ldczjobNlLFj9J4oI8X5DVy/aq2sbUG25KobSk2JSBkRP81EO4VlW/iNZAGA/8A/iTj+veBQ4RkYEi0h64iIjpqRFVHaiqA1R1APAY8CNVnRnoCYys8/Lyz1zPxWbIqebVT+e6zu39ZWAPiXj26TYQh0RSqhsQZPUi7I2WSlWpGkam8ONDiK19JwB/UtUnRWRqsoui0Ug/JhI9FALuVdUlIvKD6PnpKcrcZkk1rDPV62L4tfMHyWaajkzrPt/Fj/+10H/WUVXPvp1SbkBkhRDvJ/H7fG79OaHQJBss2O5ho+XwoxAqROTPwEnALSLSAX8rC1T1GeCZhGOOikBVv+Wnz7ZKOmGd6ea9d3MslxR72/kzLdMXNbXc9fIq/vL6auodIn6KBJxKFCSbscfu+bNHFjczH6US/x8/sPuJMopXuFaC0mhJ/AzsFxKZ5Z+qqpVAT2ByVqUympGs8lUmrnPbkexmyph6lrOdP9nO5qDP0tCgPPLuBk649VWmv7oKtyzW3TqGUza5TBxV6roTOZWIpphfYG35BFbdfDpryydYhlEj7/ETZVQFrAVOE5GfAH1V9flsC2Y0JdNhnUHq+gaxn/upDxzkWd5du42z7voPv3j8fcp6FnPVSYe4hmjuqK5Nqz5wuum7k2E+AiPf8VMx7TrgAuCJ6KH7RORRVb0pq5LlOena5YPiZrbxE9bp57pMVTtz6+fKGYuYNmcFk8cP9iVTRWU1Nz+zjKff30zf7h355lH9eWnZp9z+4krXexdFw0TfmJJaIt5k6bvT/Z2bj8DId/z4EC4GRqlqDYCIlAPvAW1WIWTCLh+UVGsN+L0uU9XOvNrH3tN5o0t5fEGFo0xVe+qY/soq/vzaakTgf79+CP1KOjL1qaVJnbSJTuCgeA3Ymfqdm4/AyGf8KIS1QEcgliayA7AqWwIVAqmmaUiF+Flp9+IwHcNFVFbV+p5d+p2VprICcZoxe+1shsh7enn5Z9x87vAm1/7slEEoyom3vsonO2s4c0Q/ppw2hNKSYsekcF793zBrScq/B7cBO5e/c8NoKVwVgoj8kUhU3G5giYi8EP1+MvCf3IiXn+SqdnDirLSyupbicIjbJ43M+Kw06ArEbcbsNPtPZFNldROZFm2o5IZZS1i4vpLhpd2585JRjBnQs0n7IGyvqm2SZykTtFS9aMPIJV4rhFgVmgXAv+OOv5I1aQqEVO35QcnlrDSofdtNtvjZv9tKIfaePtlRw2+fW84TCyvo3bUD084/jPMO35+ihDCiZKsON/ky+Y5y9Ts3jJYkacU0EekIHExkdbAq5ktoKfKhYppbVa1M55sZOGW2a50ArypiucBNNojI16+kmAH7FPPmqm3NKofdcNYwPt1Zw92vrKK+QfnaoF4s2bSTT3bUOCqimQsruHLGosAyfuOoMteqY0EdxJn6nec6GMEwYvipmOZlMmoH/Ab4DrCOSIjq/iJyH/B/quqvCngrJFfRIvk8K/WatcfCTZ3OjxnQg9+/tJKKympOHbYfYwb04LbnP/J01k4cVcoNs5Y4FqdJLFUZzz/nrW/8HN8vENhBnInfeUsEIxhGEFxXCCJyO9AVuEpVv4ge6wbcClSr6v/mTMo48mGFkCvSnZVmM9WFk2xB6NW5Pb88Y6iraSlW49jrfrFax8l8Fon9Ar7umWnGlc9tkfsaBqS5QgDOAAZpnMZQ1Z0i8kNgOdAiCqEtkc6sNNupLhJlC1r+ZeuuPZ4KJdFZ6/UuxvTv6duk5OUEzraD2BzTRr7jtUL4SFUHBT2XbdrSCiEd/M5GE1cDVXvqHE0zyWaxbvdLhlv+odg9/SpAv/f3WiGERGhQzZoJ0FYIRkviZ4XglbpiqYhc5tDpN4isEIw8xs9s1CnNhJMy8OovhlchHS/clAE4p73wun845JLkKEosjNZN1npV13QbmcBSVxj5jpfJ6ArgCRH5DpHQUyVS5rKYSNlLI4/x45AOUtkrdp2bfyFWRvP3L61kd11DZh4C/2G2E0eVMvWpJc1qI8dwWm3EnqPIIRtpfLK9TAUPWOoKI99xVQiqWgEcKSInAsOIBHQ8q6ov5Uo4wxsv5+8JQ3o3ibKJccKQvSVI/dquY7NYN//Cl7vrWP7JTv719nq6dgxzzWlDmDpraQaeMJicTvWNIfIfN9EkE78xbqBHTeWrZixq9I9kIirIUlcY+UzS1BWqOheYmwNZjAAkc/66VTmLP+5V56Bzh3bNFI1TConq2nquffJDikT45lH9ufKkQfTo3J6/vL7G06ZfWlLMrt11rjP6ePyG2aYaplvSKexqKku0aFm6CqM14yeXkZFlUgkPTbaL2cuHELtfRWV1szj+WJ0Dp/u77jtQePbKYxm0b9cmfbshROzpN8xa4vmMMXn82thTTQCYZG9mMywqyGitmEJoYVIND03mNHabLXcvDje5n7J3c5dXVI+Xg7UIGpWBn70JiTI06y8aeRQkyghSt9G7mZrcyIeNgYaRDUwh5BCnlUCq+YqSmUfcZssiNLtfTBl4hT7e8qx7YFlD9NmcSlAmIuDp+E03/DJmo4+966vi6jC4vc8guZJiqxvDaI34qo1spI9bJTG3gSiVMM9484hblbPKgGGldfUNPDBvHZt3uqewKonO+P0oA68WmTLF+KnaFo/TuwwXSbMwVgEuParM/AdGq8VWCDnCbSXghp/C8PPXbeOhtzdQr0pIhPNGRwaqceVzXU0mbvb9+LBSp/DNdkVCncOmAQHHVYcTyUz1mTLFBF11uZmanI5lUhlYojsj3zCFkCOCzH79OEJnLqzg8QUVjbPyelVmvLOBGe9uaKw57OSP8HK8zlxYweRHF1PrMPDXNyjhkDSpZxybMT/oEN4alExu0EolRYRbOGi2BmhLdGfkI2YyyhHJZr8hkUCF4Z1mwbUN2qwAfayC2LjyuQycMptpc1ZweFl3QtH6w7GVxcRRpdzy3HJHZQCR2X3n9u2amKBunzSSmyYOT3tm36NTOKOpw93kySdnsNcqxjBaClMIOSJZaocGVdaUT+CNKSf6GhiDrDi2V9U2sae/sWpbk5XFY/M3MuXx99m8w7vUxY7q2sYymZsqq5k2Z0VkVZHmzL6mNnM7m6EwUkRYojsjHzGFkCNiTt7YzDyRoLPXTM52a+oaePjdDbQPef93iIWLJjprIeJYTpVMz4zdHOr5ZIophFWM0fYwhZBDJo4q5bYLRwSevc5cWNFo8hlXPjcjs3InbjlvOOEi9wRxO2pqXc0cU88a5vhc3ziqzFfSu0zPjCeOKuWNKScGWnXlkkJYxRhtD3MqZwi/ESNBN0+5OR9vPnc4PTxSLgSltKSYcw7fHxFxTRLnFlW6qbI6ab0Cr8RzAEUiDJwyu81E21iiOyMfSVpTOd/Ix3oI2ayx7JVD3yliKBViskLTAcpvrqEencIsvO4UzzZB6iUEfXd+K7zZ4Gu0ZdKth2D4JJsRI26DaEV0Vh5vKy8pDhPyMPkAhEPC+GH7NrOvA838A36UAcCXNXVJawcEMQkFeXd+NqEF3ahmGG0VMxllgGxGjIQccvXHjkPT+Plx5XM9B/FLxpZx08SvUOSgNJwymfqltkFTTrXhht9352cTWqrpQQyjrWEKIQOkmnbZD27pIGLH400hXsa/W84dzqSxZa7n01VeflJtOJnVOrQrclRift+dH2VsIZ6G4Q8zGWWAbEaMlLoMjKUlxc1MIW6ERJg0tswxWimG2wDco1O4iXmpRyfn8FI/qTacQkHdopP8vjs/4ZsW4mkY/rAVQgbIZsSIV6oJvyUw61WTpkpwu8/1ZzatjeDmQPczgHtVC0v13fmpgZBqnQTDaGtYlFEB4BYhM3DK7KQJ42DvKsMtWimWctpvJE66ETuZjvixKKNg2Ltom/iJMjKFkAek8ge6u66eI256kZ01dUn7v2PSyCa1geMRYE35hNQET4FMh+ja4BaMbIZIG/mNhZ0WAEFDIlWVOUs+4eTfvcbOmjqSRJlSUhxm4qjSvLGjZzJE18JJg2NJ9QwvsqoQRORUEVkhIh+LyBSH85eKyPvRnzdFZEQ25UkXL6dsqgT5A13+yU4u/evb/PcDC+jQroj7vzOW3104stEklKgbYvWRIX9SJWQy4scGt+BYxJXhRdacyiISAu4CTgY2Au+KyFOqujSu2RrgOFXdLiKnAfcAR2ZLpnTIVv56P3+gn3+5m9+98BEPvbOebsVhbjx7GJeMLaNdNBld7P5e5pN8SZWQyRBdG9yCk80QaaPwyWaU0VjgY1VdDSAiDwNnA40KQVXfjGs/D9g/i/KkRbY2N7n9gSpwzM0vMXZgT15avoWqPfVcdvQArjzpEF5Z8RnHTXul2cDuFcUD3lE+uSKTET82uAXHIq4ML7KpEEqBDXHfN+I9+/8u8KzTCRG5HLgcoKzMfXNVNsnWbNQrH9GmHTXMXLSJdkVCfYPywtJPqWto4PEFFQVbaSuTK5VcDm6txXmdLytFIz/JpkJwcnc6hjSJyAlEFMJXnc6r6j1EzEmMGTOmRcKiUpmNJg4iJwzpzcvLP3P8Q3SrdQw01jKuqKzmwXnrm73E6tp6pj61pLEfv1lUc11DOEamViq5GtxaW7nLfFgpGvlJ1sJOReRoYKqqjo9+vwZAVW9OaHcY8G/gNFX9KFm/LRV2GjRcz6l9IvHX76iqZcSNz6clY7hImpTAdJPPSbZwSEDxdX1bwyvjbGwPh2HkO37CTrO5QngXOEREBgIVwEXAJfENRKQMeAL4ph9lkClSWf4HnY362UVcXVvPb59bzhc1tfzuhfQfP7EespuPw7Eec33ziYElgItgzmujrZA1haCqdSLyY2AOEALuVdUlIvKD6PnpwHXAPsDdEsneWZdMg6VLOsv/IEttv4PFph01XPtkxNzTtWM7avbUuxa6TwUnOYIMZDbomfPaaDtkdR+Cqj6jqoNU9SBV/XX02PSoMkBVv6eqPVR1ZPQnq8oAche7nspg8UVNHUhkM5mAayK5dOUIIpsNevmzh8Mwsk2b26mcq+W/0yDih9p6pbK6ln4lxVx/5jDf14VD0qwestug5SSb3+uzsTkv33HL1NrWTWlG66PNZTvN1fI/MXqopDhMTV09NbUNAOzXrSOf7KxxvT5myvJTN7k0YJSQmz8k2fWtLdomCBaZY7QF2lxyu1wn93p79efc+PRSlmzayZj+PbjuzKEctn8J4K/OcElxmN11Da4O6lxGuli0jWEULpbczoFcLf83bKviigffY9I989i+aw9/uHgUj/7g6EZlAP7MSjuqa7n53OGO/oRc27Et2sYwWjdtzmQE2V3+79pdx59eWcU9r69GVenasR2bdtRwy7PLaWjQJvf1symtX0lxo7wtvVvWom0Mo3XTJhVCNmhoUGYuquCW55bz6c7djC7rwYebdkQih3C3t8cP9snSMLS0Hdvy4BhG66bNmYyywXvrt3POn97kp48sZr9uHXn8h8fwyc4adtc1NGnnFd5aCJEshSCjYRipYyuENNi8o5pbnl3OzEWb6NO1A7ddMIJzRpVSVCQp2dtbegXgh0KQ0TCM1DCFkALVe+q557XVTH91FfWqXHHCQfzo+IPp3GHv62yr9vaW9nMYhpE6phACoKo8/f5myp9dTkVlNacP349rTjuUA3p2atbWyd4uRHwJ48rn5uVAme5g3pb3KRhGa8AUgk8+2LiDG59ewrtrtzO0bzduu3AERx24j2v7xAgiYW/u73wcKDMxmGeriJBhGLnBFEIStnxRw7TnVvDYexvp2ak9N587nAvHHEDIpbq90yzbKay0uraenz2ymKtmLMoL00omBnPbp2AYhY0pBBd219Vz73/Wcufcleypb+D7xx7Ij088mG4d3RPOuc2y3XYZ16s2aQctt2LIxGDeVv0mhtFasLDTBFSV5z78hJN/9xq3PLecow/qxfNXHcf/O/1QT2UAMPWpJY6z7JA4ryYS22U642oQ3AbtIIO5ZQU1jMLGVghxLNu8kxtnLeWt1Z8zaN8uPPDdsRx7SG9f185cWEFltXMSunpVisOhpAVzWtK0kolNZ1av1zAKG1MIwOdf7ua2Fz7i4XfW0604zHmHl/LWqs+57G/v+B7UvGb3IRFuPnd440BZJNJoLoqnJU0rmRrMbZ+CYRQubVoh7Klr4P631vL7l1ZStaeey44ewCF9unDT7GWBom1mLqzwzFqaOPh3K27HlzV1zeoXt7RpxQZzw2jbtEmFoKq8vGILNz29jNVbd/G1Qb257oxDObhPV8aVzw0UbRNzJHtRUhxuYo7ZXlVLOCSUFIfZES2GY6YVwzBamjanED7e8gU3Pr2M1z76jAN7debeb43hhMF9iNZ0Dhxt4xSuGU+4SBDBsah95w7tWHT9KSk+SWawncWGYcRoMwqhsmoPd7y4kgfmraNT+xC/nHAolx09gPbtmgZaBQ2d9HIElxSHmXrWMK6ascjxfEVlNQOnzG6xgdh2FhuGEU+rDzutq4/4CY6/9RXuf2stFx1xAK/8/Hi+d+yBzZQBBA+ddFMUpSXFLLr+FCaOKvV0Fit7B+Jc1yf22oxmGEbbo1UrhNdXfsbpf3id655cwqH7dWP2/xzLr88Zzj5dOrheEzTFsx8F4qcyWksMxLaz2DCMeFqlyWjN1l38evYyXlz2KWU9OzH9G6MZP2zfRj9BMoJE2/gJ10xs41bFOtcDse0sNgwjnlalEHbW1HLn3I+57401tA8VcfWpQ/jOVwfQoZ337Dxd/CiQ+DZuxepzPRBbBTTDMOJpFQqhvkF5ZP4Gbp2zgm1Ve7hg9P78fPxg+nTt2NKiOZIvA7HtLDYMI56CVwjzVn/OjbOWsnTzTsb078HfzxzL8P27t7RYnuTTQGyb0QzDiFGwCmHDtipufnYZz3zwCaUlxfzx4lGccVhf336ClsYGYsMw8o2CUwgNqkybs5y/vL6GkAg/PXkQl3/tQDomieIxDMMwvCk4hbDiky+46+VVTBzZj6tPG0Lf7hYRYxiGkQkKTiGEQ0U88aNjOLysR0uLYhiG0aoouI1pB/fpYsrAMAwjCxScQjAMwzCygykEwzAMAzCFYBiGYUQxhWAYhmEAIOpQ2zefEZHPgHUZ7LIXsDWD/WWSfJUtX+UCky0V8lUuMNlSwU2u/qra2+vCglMImUZE5qvqmJaWw4l8lS1f5QKTLRXyVS4w2VIhHbnMZGQYhmEAphAMwzCMKKYQ4J6WFsCDfJUtX+UCky0V8lUuMNlSIWW52rwPwTAMw4hgKwTDMAwDMIVgGIZhRGkzCkFEThWRFSLysYhMcTg/RETeEpHdIvLzPJLrUhF5P/rzpoiMyCPZzo7KtUhE5ovIV/NFtrh2R4hIvYicnw9yicjxIrIj+s4Wich1uZDLj2xx8i0SkSUi8mq+yCYik+Pe2YfR32nPPJCru4jMEpHF0Xf27WzLFEC2HiLy7+jf6Dsi8pWknapqq/8BQsAq4ECgPbAYGJrQpg9wBPBr4Od5JNcxQI/o59OAt/NIti7s9UMdBizPF9ni2s0FngHOzwe5gOOBp3PxnlKQrQRYCpRFv/fJF9kS2p8JzM0HuYD/B9wS/dwb2Aa0zxPZpgHXRz8PAV5K1m9bWSGMBT5W1dWqugd4GDg7voGqblHVd4HaPJPrTVXdHv06D9g/j2T7UqP/24DOQK4iFJLKFuUnwOPAljyTqyXwI9slwBOquh4ifxN5JFs8FwMP5YlcCnSVSO3eLkQUQl2eyDYUeAlAVZcDA0RkX69O24pCKAU2xH3fGD3W0gSV67vAs1mVaC++ZBORc0RkOTAb+E6+yCYipcA5wPQcyeRLrihHR00Mz4rIsNyI5ku2QUAPEXlFRBaIyGV5JBsAItIJOJWIos8Hue4EDgU2AR8A/6uqDXki22LgXAARGQv0J8mEsq0oBHE4lg/xtr7lEpETiCiEq7MqUdwtHY41k01V/62qQ4CJwK+yLlUEP7LdAVytqvU5kCeGH7neI5JTZgTwR2Bm1qWK4Ee2dsBoYAIwHrhWRAZlWzCC/X2eCbyhqtuyKE8MP3KNBxYB/YCRwJ0i0i3bguFPtnIiCn4RkdXyQpKsXgquhGaKbAQOiPu+PxGN3tL4kktEDgP+Cpymqp/nk2wxVPU1ETlIRHqparYTfvmRbQzwcGQlTy/gdBGpU9VsDsBJ5VLVnXGfnxGRu/PonW0EtqrqLmCXiLwGjAA+ygPZYlxEbsxF4E+ubwPlUdPpxyKyhoi9/p2Wli36f+3bAFGT1projzvZdn7kww8RxbcaGMheB8wwl7ZTyZ1TOalcQBnwMXBMvr0z4GD2OpUPBypi31tatoT2fyc3TmU/72y/uHc2FlifL++MiOnjpWjbTsCHwFfyQbZou+5EbPSdsy1TgHf2J2Bq9PO+0b+BXnkiWwlRBzfwfeD+ZP22iRWCqtaJyI+BOUS88/eq6hIR+UH0/HQR2Q+YD3QDGkTkSiJe+52uHedALuA6YB/g7uhst05zkGHRp2znAZeJSC1QDUzS6P++PJAt5/iU63zghyJSR+SdXZQv70xVl4nIc8D7QAPwV1X9MB9kizY9B3heIyuYrONTrl8BfxeRD4iYca7W7K/2/Mp2KHC/iNQTiR77brJ+LXWFYRiGAbQdp7JhGIaRBFMIhmEYBmAKwTAMw4hiCsEwDMMATCEYhmEYUUwhGK2KaBbMRXE/A0RkjIj8IXr+eBE5Jq79RBEZmsJ9vsyQvBnpxzAyQZvYh2C0KapVdWTCsbVE9phAJNvol8Cb0e8TgaeJxGkbRpvGVghGqye6KnhaRAYAPwCuiq4ejgPOAqZFvx8U/XkumtztdREZEu1joETqZbwrIo45m0TkFhH5Udz3qSLyMxHpIiIvich7IvKBiDTL5BmTMe77nSLyrejn0SLyalSmOSLSN3r8f0RkaTTf/cMZe2FGm8VWCEZroziazAtgjaqeEzuhqmtFZDrwpareCiAiTxGpT/BY9PtLwA9UdaWIHAncDZwI/B74k6reLyJXuNz7YSJJ9e6Ofr+QSGbOGuAcVd0pIr2AeSLylJ8dyiISJpIE72xV/UxEJhGp2fEdYAowUFV3i0iJ3xdkGG6YQjBaG04mI1+ISBciBYkejaYJAegQ/XcckVQdAA8AtyRer6oLRaSPiPQjUixlu6qujw7qvxGRrxFJCVFKJO/NJz7EGgx8BXgh65lahAAAAV9JREFUKlMI2Bw99z7woIjMJHdZU41WjCkEw9hLEVDpoVD85Hl5jEi+ov2IrBgALiWiIEaraq2IrAU6JlxXR1MTbuy8AEtU9WiHe00AvkbE7HWtiAxT1VwUZzFaKeZDMNoaXwBdnb5HExmuEZELIJIyWPbWsH6DSOpliAzwbjwcbXc+EeUAkSydW6LK4AQihUoSWQcMFZEOItId+Hr0+Aqgt4gcHZUpLCLDRKQIOEBVXwZ+QSSzZRdfb8AwXDCFYLQ1ZgHnRJ3IxxIZwCeLyEIROYjIYP9dEVkMLGFvWcL/Ba4QkXeJDPCOqOoSIgqmQlVjpp0HgTEiMj/a/3KH6zYAjxA1AxEpZoJGyiOeD9wSlWkREbNWCPhnNMvmQuB2Va1M9aUYBli2U8MwDCOKrRAMwzAMwBSCYRiGEcUUgmEYhgGYQjAMwzCimEIwDMMwAFMIhmEYRhRTCIZhGAYA/x9n/lF3RxdVugAAAABJRU5ErkJggg==\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO29fZwdZXnw/712c0g2gGyQoLISEtSGBxpJTJRYnlagalALBrBQi61vfaifah9BTA0WS/DlIf6ipU8ff76gtT6KYlDoCmINVUBblNbETYhR0iogsFEJkhVIFthkr+ePmdnMzs7LPXNmzplzzvX9fPaz55w5M3PNfWbu676v+3oRVcUwDMMw8tDXbgEMwzCMzsOUh2EYhpEbUx6GYRhGbkx5GIZhGLkx5WEYhmHkxpSHYRiGkRtTHkbbEZEdInJawrbTROShks5zh4j8WRnHajcioiLy/HbL4UI3tbtxEFMehjMicr+IjIvIEyLySxH5nIgc1uxxVfUkVb2jBBELIyLrRGRCRB73//5TRD4mIs9pp1ydQuTe+JWI/GPee0NEFvpKcVZVchrlYcrDyMtZqnoYsBRYBlzWZnnKZKOqHg4cCZwDPBvYYgrEmeDeeBHwYuDyNstjVIgpD6MQqvpLYBOeEgFARGaLyEdE5AF/9PlJERnwtx0lIl8XkTEReVRE/lVE+vxt94vIy/3XA/6MZo+I/BivEyJ0jmnmGv+7H/Rfz/PPsdvf/+si8twC1zahqjuAC4DdwKWh8/2BiGz1r+N7IvLC0Lb7ReQyEfmxf/5/FJE5OfZ9t4jcLSK/EZGNkX3XiMgvRGSXiLwl0iZp7X6aiDwkIpeKyMP+Md4c2ndARD4qIj/3z/tvoX1X+nKOici2JNNiTPuNAv8M/HZ0m4j0icjl/vkeFpHPi8gR/ubv+v/H/BnMS13OZ7QHUx5GIfxO+VXAT0Mffxj4LTyF8nxgCPgbf9ulwEPAfOBZwHuBuNw4VwDP8/9WAW/MIVYf8I/AccACYBz4WI79p6GqB4CvAb8LICIvAj4L/DnwTOBTwE0iMju024W+3M/Da4vLc+x7PnAmsAh4IfAmf98zgXcDrwBeALw8Impau4M3gzrC//ytwP8vIvP8bR8BlgO/gzfj+itgUkSGgFuAD/qfvxu4QUTmZ7WbiBwLvBoYidn8Jv/vdOB44DAO/ka/5/8fVNXDVPX7Wecy2oiq2p/9Of0B9wNPAI/jdfzfxnvQAQTYCzwv9P2XAvf5r9+P1xE/P+G4L/df3wucGdp2EfBQ6L2GjwF8DvhggrxLgT2h93cAf5bw3XXAtTGfvw34L//1J4APRLbvBF4Wuo63hba9GvhZjn3fENr2/wGf9F9/Flgf2vZbQTs4tPtpeEp0Vmj7w8BKPGU7Dpwcc93vAb4Q+WwT8MaMe2MM+DnwcWAg2u7+PfMXof0WAxPALGChf12z4s5hf/X6s4UpIy+rVfVbIvIy4EvAUXgdxnxgLt4aQfBdAfr91xvwOuhb/e3XqOr6mOMfAzwYev9zV8FEZC5wNd7oPRhZHy4i/erNIoowBDzqvz4OeKOI/GVo+yG+zAFR2YNtLvv+MvR6X2jbMcCWyHEDstod4Nequj9y7MPwfrs5wM+YyXHAH4rIWaHPGsDtMd8NWK2q30rZDt61hOX/OZ7ieFbGfkbNMLOVUQhV/Q7eqP8j/keP4I1iT1LVQf/vCPUWUFHVx1X1UlU9HjgLeJeI/H7MoX8BHBt6vyCyfR9eZxnw7NDrS/FGsqeo6jM4aAYRCuCvyZwF/Kv/0YPAh0LXN6iqc1X1utBuUdl35dg3ibQ2SW33DB4BnsQzsUV5EG/mEZb30ASFn4ddeIopYAGwH/gV8WZMo6aY8jCa4e+AV4jIUlWdBD4NXC0iRwOIyJCIrPJf/4GIPF+84fFjwAH/L8r1wGX+4vdzgb+MbN8K/LGI9PtrAS8LbTscryMdE5Ej8dZPciMiDRH5b8B1eMrpb/1NnwbeJiKniMehIvIaETk8tPvbReS5/vnfC2zMsW8S1wNvEpET/dnV1HVltXsa/r6fBf5WRI7x2/Sl/jrMtcBZIrLK/3yOv/ie2wEhwnXAJSKySDxX3v+F5+W2H885YRJvLcSoOaY8jMKo6m7g88D7/I/eg7eAfpeIPAZ8C28mAN5C77fw7OLfBz6u8bEdV+KZMu4DbgW+ENn+TrzZwBje4vRwaNvfAQN4I+q7gG/mvKQLRCSw298E/BpYrqq7/OvdDPwPvAXePf61vilyjC/5ct/r/30wx76xqOo/+9d2m7/fbZGvpLV7Fu8GtgM/wDPPfRjoU9UHgdfiKcDdeDORNTTfZ3wW7zf9Lt5v/CT+AEFV9wEfAu70PbxWNnkuo0JE1WaKhlEGInI/3sJwlt3fMDoem3kYhmEYuTHlYRiGYeTGzFaGYRhGbmzmYRiGYeSmK4IEjzrqKF24cGG7xTAMw+gotmzZ8oiqZqaciaMrlMfChQvZvHlzu8UwDMPoKETEOYNDFDNbGYZhGLkx5WEYhmHkxpSHYRiGkRtTHoZhGEZuTHkYhmEYuekKbyujMxgeGWXDpp3sGhvnmMEB1qxazOplQ+0WyzCMApjyMFrC8Mgol924nfEJLwv76Ng4l924HcAUiGF0IGa2MlrChk07pxRHwPjEATZs2tkmiQzDaAZTHkZL2DU2nutzwzDqjSkPoyUcMziQ63PDMOqNKQ+jJaxZtZiBRv+0zwYa/axZ5VrwzjCMOtFW5SEinxWRh0XkR6HP1onIqIhs9f9e3U4ZjXJYvWyIq85dwtDgAAIMDQ5w1blLbLHcMDqUdntbfQ6vpvPnI59fraofab04RpWsXjZkysIwuoS2zjxU9bvAo+2UwTAMw8hPXdc83iEid/tmrXlxXxCRi0Rks4hs3r17d6vlMwzD6GnqqDw+ATwPWAr8Avho3JdU9RpVXaGqK+bPL1TLxDAMwyhI7ZSHqv5KVQ+o6iTwaeAl7ZbJMAzDmE7tlIeIPCf09hzgR0nfNQzDMNpDW72tROQ64DTgKBF5CLgCOE1ElgIK3A/8edsENAzDMGJpq/JQ1dfHfPwPLRfEMAzDyEXtzFaGYRhG/THlYRiGYeTGlIdhGIaRG1MehmEYRm5MeRiGYRi5MeVhGIZh5MaUh2EYhpEbUx6GYRhGbkx5GIZhGLkx5WEYhmHkxpSHYRiGkRtTHoZhGEZuTHkYhmEYuTHlYRiGYeTGlIdhGIaRG1MehmEYRm7aWgzKMAwjjeGRUTZs2smusXGOGRxgzarFrF421G6xDEx5GIZRU4ZHRrnsxu2MTxwAYHRsnMtu3A5gCqQGmNnKMIxasmHTzinFETA+cYANm3a2SSIjjCkPwzBqya6x8VyfG63FzFYdjNmDjW7mmMEBRmMUxTGDA22QxohiM48OJbAHj46Noxy0Bw+PjLZbNMMohTWrFjPQ6J/22UCjnzWrFrdJIiOMKY8OxezBRrezetkQV527hKHBAQQYGhzgqnOX2Oy6JpjZqkMxe7DRC6xeNmTKoqbYzKNDSbL7mj3YMIxWYMqjQzF7sGEY7cTMVh1KMJU3byvDMNqBKY8OxuzBhmG0CzNbGYZhGLmxmYdhxGABmIaRjikPw4hgCfkMI5u2mq1E5LMi8rCI/Cj02ZEi8i8i8l/+/3ntlNHoPSwA0zCyafeax+eAMyOfrQW+raovAL7tvzeMlmEBmIaRTVvNVqr6XRFZGPn4tcBp/uv/C9wBvKdlQrUQs6vXE0vIZxjZtHvmEcezVPUXAP7/o+O+JCIXichmEdm8e/fulgpYBpbYsL5YAKZhZFNH5eGEql6jqitUdcX8+fPbLU5uzK5eXywhn2FkU0dvq1+JyHNU9Rci8hzg4XYLVAVmV683FoDZHGaS7X7qOPO4CXij//qNwNfaKEtlWGJDo1sxk2xv0G5X3euA7wOLReQhEXkrsB54hYj8F/AK/33XYXZ1o1sxk2xv0G5vq9cnbPr9lgrSBiyxodGtmEm2N6jjmkfPYHZ1oxsxV+feoI5rHoZhdDBmku0NbOZhGEapmEm2NzDlYXQM5v7ZOZhJtvvJNFuJyB+KyOH+68tF5EYReVH1ohnGQcz90zDqhcuax/tU9XER+e/AKrx8U5+oVizDmI65fxpGvXBRHsET+xrgE6r6NeCQ6kQyjJmY+6dh1AsX5TEqIp8Czge+ISKzHfczjNKwiHzDqBcuSuB8YBNwpqqOAUcCayqVyjAimPunYdSLRG8rETky9PaO0GdPAZurFcswpmPun4ZRL9JcdbcACkjof4ACx1col2HMwNw/DaM+JCoPVV3USkEMwzCMzsEpSFBE5gEvAOYEn6nqd6sSyjCMzsYCOrufTOUhIn8GvBN4LrAVWImXRv2MakUzDKPuxCkJgMtu3D4VlxMEdAKmQLoIF2+rdwIvBn6uqqcDy4DOKxpuGEapJEX9r7tphwV09gAuZqsnVfVJEUFEZqvqPSJi/pFGbsyU0V0kRf1HPwuwgM7uwkV5PCQig8Aw8C8isgfYVa1YRrcRjFLNlNE95FUGFtDZXWSarVT1HFUdU9V1wPuAfwBWVy2Y0V1YbqruI0kZzJvbsIDOHsAlq+6C4A+4D2/R/NmVS2Z0FZabqvtIivq/4qyTuOrcJQwNDiDA0OAAV527xGaYXYaL2eoWDgYJzgEWATuBkyqUq2vpVbu/lSbtLFzu06yo/164r3uZTOWhqkvC7/1aHn9emURdTC/b/desWjzt2qH7TBndMjDIc59a1H/vkjs7rqr+EM9118hJL9v9Vy8b6mpTRjcVq+rl+9RwxyVI8F2ht33Ai7A4j0L0ut2/m0epaR1up11zr9+nhhsuM4/DQ3+z8dZAXlulUN2K1aToXrqpw7X71HDBxVX3ytDfh1T1i6r6ZCuE6zasJkX30k0drt2nhgtp9TxuxvOyikVVz65Eoi6mzjUpumWxt13U2SEg729b5/vUqA+iGq8fRORl/stz8eI6rvXfvx64X1XfW714bqxYsUI3b7b6VEWJeteA1/F104J2K6ijArbf1khDRLao6opC+yYpj9DBv6uqv5f1WTsx5dEcp66/LTYGY2hwgDvXdm/y5FZ39u1QLr362xpuNKM8XIIE54vI8ap6r3+yRcD8Iicz6kk3Lfa60uqYm3bF+PTib2u0Bhdvq0uAO0TkDhG5A7gduLhSqYyWUtVi7/DIKKeuv41Fa2/h1PW31SrmodWxDO2KneimhXyjXrhEmH9TRF4AnOB/dI+qPlWtWCAi9wOPAweA/UWnVt1ClSaPKhZ76x5N3+oReRnnK3IPnH7CfK6964HYzw2jGdK8rc5Q1dtE5NzIpueJCKp6Y8WyAZyuqo+04Dy1puqOuArvmroHzbU611az5yt6D9x+T3w8b9LnvUAdHRs6kbSZx8uA24CzYrYp0ArlYZDcEV95847Sbvqyo7/rbmtvtWtts+e78ubk6nxpv1vdf4dWU/cZcSeRqDxU9Qr//5tbJ850EYBbRUSBT6nqNeGNInIRcBHAggUL2iBe60h60Pfsm2B4ZDT1pm/XKKsOWXTTrj34v+6mHYyNTwAwp5E71ZszzczuhkdG2bNvInZblhKow+9QJ+o+I+4kXOp5vFNEniEenxGRH4rIK1sg26mq+iLgVcDbRWSaa7CqXqOqK1R1xfz53W2/TXvQ0xZc25msr91Ryq7X/tT+yanXe/ZNVNo+q5cNcefaM7hv/Wu4c+0Zzp1V2m+cpQTa/TvUDZuJlYfLUOstqvoY8ErgaODNwPpKpQJUdZf//2Hgn4CXVH3OupL2oKfd9O3MjtruLLou194p2WPTfuMsJeD6O9TZM65MzPusPFziPMT//2rgH1V1m4hI2g7NIiKHAn2q+rj/+pXA+6s8Z51ZvWxomnklTNpN3+5RVjuz6Lpce9XtEzabHTHQQATG9k3kNh8mmZ4GBxqltG/d1wHKNL3WOY1Mp+Ey89giIrfiKY9NInI4MJmxT7M8C/g3EdkG/Adwi6p+s+Jz1pp1Z5+U2/zQy6Msl2uvsn2iZrOx8Qn27JsoZD5MMj2tOzu7mKeL+a7OM7CyTa/tnhF3Ey4zj7cCS4F7VXWfiDwTz3RVGX40+8lVnqPTKLLg2sujLJdrr7J94jrkMHkWaZtZbHdZIM4zA2u1A0YVC9zdXFemlbgoDwVOBP4Az3R0KF4tc6PF5L3pOz07atBRjY6N0y/CAVWGHK/B5dqrbB8X01ce81jRDs9FMWR5ZIV/B+Fgqu1WmLfabXo1knFRHh/HM1Odgac8HgduwErR1proCPHqC5Z2jNKAmXb4A34CzzwdlkuHW9UoNKlDjn6nalxcddNmYNHfIZpGtWo3V3M1ri8uax6nqOrbgScBVHUPcEilUhlN0Q31tNPMPnWxx0Oyl1LcOkWYVpkPXVx109YBssxvUO0swFyN64vLzGNCRPrxBx0iMp/qF8yNJqg6EKoVdu+sDqlqs4XLNbp4KZXhbdUMrqa5pBmYSztXOQvodNNrN+OiPP4eL87iaBH5EPA64PJKpTKaoko7cavcOrPMPlV2WK7XmKWkq16YdVXiUTmC2ZJLZ5z1O7RiFmAL3PXEpYb5F4G/Aq4CfgGsVtWvVC2YUZwqXVBb5daZZvapusNKusaLN26dZpqqSkm7BOwVNU3m3S/udwiCvHrdzbVXAiuTcJl5oKr3APcAiMigiPy1qn6oUsmMwlTpglplhxkdRV917pLC3lbNkHYto2PjrPnKNqCaxdyyZj1J5N2vlWajMs2hVZtW6x5Y2QrSUrIfC7wPOAYYBr4EfAD4E+C6lkjXQ5R5s1fxwAfyJRUtjrp15jlv0oN41blL2lIqNctUMzGprLtpB+vOPql0Je3auRdV4kX2a4XZqMzOuBUduyVYTJ95fB74Dp5b7pnAXcAO4IWq+ssWyNYRlNHpV3Gzxz3wRWWNyhclya3T9Trq9iDGzdyijI1PlKKko79JktKKdu5FZz11dX0t8x7Ie6wiz4XFn6QrjyNVdZ3/epOI/Ap4cSuqCHYKZXX6ZXeecQ8DUFjWNHfNsBnp1PW3dUXNibBSyIrVcB2Vu/4m4SC8MNHOvahpsq5ZB8q8B/JGzBd5LuqqhFtJ6pqHiMzj4PrYL4G5fqJCVPXRimWrPWV1+mU+OHEPwyUbtzKn0cf4xHQP6/GJA7zr+q1cvHErAAONPuY0+me4kybJITDNrFT0OhIT/81tOHsFFSGr3sfqZUMse/+tsbU05s1t5DpPXAc1e1bfjPtHYYYCievci8566ur6WmZnnOdYRZ/huirhVpKmPI4AtnBQeQD80P+vwPFVCdUplNXpl/ngxD0MCjMUR8BkqJcan5ic+l54BOYqX9HriHsQG/3CE0/un+q4y7Zbu444rzjrJNZ8dRsTBw42VKNfuOKs7KSEAUkdVNJsTvFmdHldcF0J7xco0Es2bm2rIimzM85zrKLPcF2VcCtJqyS4sIVydCRldfplPjhlmnqCEZirfEWvI+5B3PvU/hkp6McnDnDp9dum7VOUdTfFl3WNlvYto5PI+5sMDQ5kOgrUda2tqFxldsZ5jtXMM9zr8SdOrrpGPGV1+mU+OC45lfKwa2zcSb6g0xifOFDIrTb6IC5ae0vs9w6oltLBxdVGgfjSvlmdRFaHmfSbzJvb4MmJydz3T5FOP07GKtbamlFGZXbGrscy81NxRDXJ+bJzWLFihW7evLkt567Snzzr2EmLsJds3Bq76DrXX/fI84u7joLjHsBmAshOXX9bqhJ0kasVx8669uGRUa68eceMdZPgO5B/0JAkf5LcSTImmc0EuG/9a1JlKEOuqIztMgG189ztRkS2qOqKIvvazKNJqpq6Zo3i0mIjLly5gC/e9cCMRdfZjT72Jax9xNHoE6cRWBWutlnuss2Y58rMm5UVcR93DYMDDdadfdK0xfkAl9Qhru68WTIGM8QoRT2Giq4ftDvgrtfNT0VxyaqLiPSLyDEisiD4q1qwXierU0rb/sHVS7j6gqUzsqSOxXgNpXHYnFlOD1UVrrarl3mZXvsTKh434xLpEgvhStq1J7k4Hzo7vl1dUocMj4ySVAM6Se4kGQ+olpqxtmhanDpXMjSSyZx5iMhfAlcAv+JgNl0FXlihXD1PVoectT1uNOUStxDGVdnkWXTMYyIIPi/bJp02q0k6dpLcadeeV6m6zOCSovzFv644kmQcCq19tLM+eN3ifPLSq2YvF7PVO4HFqvrrqoUxDpLVIRfxEjn9hPmx5qw5jb7YWAbXEXjScaOdRhHzRLPOBGkPtmverDS54zpM8b+T1yzk0okmfUdJbsO0Tr3sRWrI91sNj4zSV7L5rJW02+TWTlyUx4PAb6oWxJhO1igu7yhveGSUG7aMTuvgBThv+RArjjuy8Og+7bhxM58iayNFO7isBzsa6xCYmoLtLnIHC8FxZVrjOsRGf/I6ksuAIM1zK4lWxiTk+a2C3yeunTrF46luqXVaiYvyuBe4Q0RuAaZSk6jq31YmlZH5wOftEJKCB2+/ZzcfXF3M6yfruFFaZZ4I19yOEn2wXUaOWQvUQYeZ5cUFxOce8XEZEKxZtXhG0CLAE0/un+FiHKaOi8JJa0L9Ih2T6r3TTW7N4KI8HvD/DsHKz7aUrAc+T4dQZI0kjbQOOul8SaNmxXPzLGM0nJXEMSpb1sgxWKCO6/Ojcrt0GBOT2lT689XLhlh3044ZcSppx60rSe01qdox19HLOa4ylYeqXgkgIod7b/WJyqUySqfMm9ylg447btpCdVm2Ypea22HZspRfWhr6YP+sNC5x+0DyekzW9f8mIcCx00a73dDx9nKQYaarroj8toiMAD8CdojIFhFxT+xj1IK4inBFb/KsDjrpuIH77VBC51CGe2ZWBxqWzcXt1aVDDuQ+/YT5iccL0y9SuBJgWDbXz+tK3nuyjpX7Vi8b4rzlQ1Mu5f0iset93YhLnMc1wLtU9ThVPQ64FPh0tWIZZRPuuMOxH0Vu8rQONeu4q5cNcefaMxI72WZGz4HnjqtsabOKfU/vZ9HaW1KPF2Z0bJwvRTzOkjig2lRsQ5kDARdcO+28nXuee7IZZVslgcNIsOh/QJUbtoy2Xa5W4LLmcaiq3h68UdU7grTsRusow5e8rEXTtLgB17QeZZssLh/ePsNdOCApVUqaogpcl+M8gZJwjd3vk2Rz2ejYOIvW3pKZTRdaVxrWxRW1qMtq2j0Zvufj3Hnr4NVk3lbp3Csi7wO+4L9/A3BfdSIZUdrpSx6ntMqw85ZpKx4eGU1UHH1C4mg2TxLJfhEmVTlioMHep/fP8HbKw2TGrsHIes1XtnHlzTtm1FeB1nlPuXaOVSdZTFLi7V7n6WVvKxez1VuA+cCNwD8BRwFvrlIoYzrtSt+QZCoAmjaBJZksgNx27TTzU1pHHWf+SWJSlfvWv4atV7ySDa872WmfZpmYVPbsm2irmca1cyy7E3VxfID2r/N0y/pTEVy8rfYA/xO8HFd4ZqzHqhbMOEi7RjdZwXF5R5Rxs5iwmavoDCurHd57493Tznv6CfO5/Z7d7Bob54iBBnMafVOj+7g6IjC9M1i9bCh3qpd5cxuxUfx5KMMcktf8WXUhsCRc7u06eDWZt1UKIvIlEXmGv86xA9gpImuqF80IqHp0E17oXHrlrSx7/60sWntL7uytWefIWvAsOsPKaod9E5PTznvtXQ9MvR8bn+DJiUmuvmApd649g3VnnzRjNiJ4KVjC5Jm1ADw5McngQHIUeFICyCjNOhTkXXROus69T+2ftl/Zi/hJv2m/SNMOH2VSpiNKp+Gy5nGiqj4mIhcC3wDeg1eedkOlkhlTVDm6iY72k4okhSmrPG50JF10hpWVvj2L8YkDrLtpx9SIfE5j+phKgRu2jLLiuCNTI/yDGU1SZPucRh+NfpmxXtLoEy54ybHcsGXUKUalqPNEkXWJ4PNoTZKx8YnYVC/hfGFhxZ+3M0265+vYMdcxer8VuCiPhog0gNXAx1R1QkQqryAlImcC/xvoBz6jquurPmddKdO7Jtrx7Nn7VGJ98zhcK91FZU1SAKNj41NR2kVNH0kdXB7GxiemFGdce0TT4acV6Lp449b4c+yb4OoLlk6TM1zbY8VxR04dO25hfqDRz+knzC/sPNFMve4Nm3bOaNuo4onLglzUuaOVHmVGMTIrCfop2dcC24DXAAuAa1X1dysTyltb+U/gFcBDwA+A16vqj+O+385Kgp2ES2R4EgJOD3Ccy2xa5t7wd85bPjRj9J13tLns/bc2vbaQRrQCX7RqYFr75q1+ePnwdq779wc5oEq/CK8/5djEmY3LsZup9Ldo7S2JqeDDVQebOYfRepqpJJiqPESkD3idql4f+kyAflXdX+SETkKJvBRYp6qr/PeXAajqVXHfP/zww3X58uVVidM1jDwwxlP7i5l2Zs/q59gjBzjqsNmJ33nkiaf46cPx2Wtm9fUxqcpkyv0WnOPBR8d5av8Bp3NGue+RvfzqsSedv58HQdCYLnT2rH6WLRjMbN9Z/X0sfObcadfzyBNPxV7vI088xb2796a2V5SVxz8zdXtS2zzrGXNYdFR66FbStQXXHnDXvcmVG7LkqxtJv0038Z3vfKeaMrSqOiki7wCuD32mQGWKw2cILxV8wEPAKeEviMhFwEUAs2d31w9aFUUVR7Dvvbv3AiQ+QPc/si9x//2Tkzz/6MOmHsakcxx12OymHtCkAlZBxz97ljcLSsoPlUSfH+cRR3A9We27/8DktDaMKohwGz/46HguxTF7VvbifVLbuBT9OvbIgRnKrE+EY4+cblKcPas/sR1GHhgrpQNuRaee9tt0mwIpisuax7+IyLuBjcDe4ENVfbQyqYjNXjHtSVLVa/BSp7BixQq94447KhSn/oTXGQbnNlD1EuiFTU1OKcPxAuuS4iOOGBzgjgTzw8K1tyQeM2y2SDNtJB3bhbT1hrB5xaUdGn3CYXNmTbnwnn7C/CkTUpLcru17ROj7R8d8/4jBAfb6HlEu9An87flLM017aaanO0KmpyRcFuqHR0ZZ85VtTCTcQAca/VzcxKJ3YBqcFzINNnvMONJ+m2bu0bohjl5+cbgoj7f4/0vrxKUAABqVSURBVN8e+kyB4wufNZuHgGND758L7KrwfJXQqvKUUVt72OafVfUujj8+ZUFixPausfFC1xWtSVG291jQBkmEF93TFojj1nZcixa5tq9LKeE80e9pVQTD5HVIyIrLSSSlP2o2VqVV6UB6OXLcFZcgwUWtECTCD4AXiMgiYBT4I+CP2yBHYVqZUiQrGjeu6l1SviDwCjkldTRHDDRmXNfFG7ey7qYdHHpIP3ufninHoYf0z6hJEZajDMWa1Qb7nj5YLClvbq60okXnLfc8kS7ZuJVjBgc4b/lQ4qJ2gEsp4TWrFnPJxq1Os48461aRtDLR2esTT+6fmkG43r8bNu3MTN3STAeclhOsTLohXXzVuKQnCdKyny8ifxr8VSmUvxj/DmAT8BPgelXdUeU5y6aVKUVcHsZw0ac7157Bfetfk2hT3zU2nhj0JUJsRzo2PsHT+yfp75s+7Gz0Cx86Z8mM74flKBKtHidzGnv2TUwFxOUNaEs6dpBBNRx0d8MW7/hJAYEC02YqSXKsXjbEhSsXOKd4D1MkrUx0nz37JmaYnlzuX5d7sZkOOCmY0jXI0pVWZy7uRDJnHiJyBXAacCJekOCrgH8DPl+lYKr6Df98HUkrp70uJo64B3YwIWXGMYMDibODSxLWFMDLxTTQ6ONpZcq99IIXH+ukGOJGynHnTzqWSxuMTxzg0uu38dHzT+aqc5c0fewgEC56jnU37WDv0/E+JReuXOBcSviDq5dMi/2YmzCze/0pBy28wyOjXHr9tsQMtEmK2jWXVNb9m/U7NNsBJyVIzJP92AWLM8nGZc3jdcDJwIiqvllEngV8plqxOp9WTnuzbO1x6TWGR0Z54smZHVyjX6Ye7rjI2aycTuEAu2BkHo7MDp8/KSBudGycNV/dBoqz2cTVzHNAlctu3M5V5y5xjjtIMvcktXdSlP68uY2pevEBeUoND4+M8q6NW2ekfr/2rge4/Z7dnH7C/Gm1JaKkdfyug5qs+zeurYIyvkMldMBDKSbHsunVyHFXXMxW46o6CewXkWcAD1PtYnlX0MppbzS/ztyE9BrRPFJxHjGHHjIr9YHJm9MpztQRNZGMjU/MsJNPHNBcZpPVy4acvZPymg+T8hfl7bBcXGLT2LBpZ2LNkNGxcb541wO5SwO7bAtwuX/j2urqC5Zyf0nmSTMn1QeXmcdmERnEqx64BXgC+I9KpeoCWl2wJ3yefU/vZ18kxcb4xAGuvPlg/qakjjaIf8iqr/3eG++ecY4koqNaVxOJy7HCJI1K8x4H3D2N4mYkSdH0Lh10midblsxpyjOrgz39hPlce9cDMz6f2+hjfGIy1/1b5YjdzEn1wcXb6i/8l58UkW8Cz1DVu6sVqzuo8iEKOpnRsfEpswCke53s2TeRmbojSLyX5ikW5DraV9DU0cy6T1oHnCdBYtpxXD3lVi8bYvPPH52WQuS85V6OqiKuyFnnzeO+G6ZfJDPFy+337I79fN6hs/lxzeIazJxUD1xSsouIvEFE/kZV7wfGROQl1YtmJBE2+0D6iDMPQQfn4inmqgDiOs1m1n2iazdwMKX8JRu3MntWH/PmNqZMJm9YuSC3mcPVUy6pfjUUK5aVdd68JkPw1hs+ev7JlSVNdCVvfXOj/riYrT6OV575DOD9wOPADcCLK5SrJ3ENvmvG7BNHNDAuyaMq3JG4jILjFkiHR0bZF+OJFI7oTlOG0RFyXEr5gUY/V19wMOI67LHkYuZIywC89MpbpyL39z61P7azv3jj1kKLw1kdeDTteXjGmURVAYR5aGcZZaM6XJTHKar6IhEZAa+yoIgcUrFcPUfcA3bJxq1s/vmjMzx0XEaDgwMNDp09i9Gx8dR0I3GBcS4diYuHV/S4SVlnw2nJIT19SHDtYbNdlLhU4Xk6qTTFGHhSZSnOIh1kng5c8AI2RbyF+KSAz7QCVGGqrBnTqqhwo7W4eFtN+CnSFUBE5kOi04dRkLgHTIEv3vXAjCl+1mhwoNHPurNPYs2qxTT6JVFxJHUOLh4tgVdNUnBWnIxJM6ZDZ0/38FqzanFicFx4TSatAw8r2LwmkyLmoTjyenXFnTfsZh3npRZUQfzo+SfT6JvZanuf3u9kIqqyIp6l+uhOXGYefw/8E3C0iHwIL+7j8kql6kGSHiSFGSM0V1/6U9fflpgqIm0RNcujJS2NBcQrpeGRUeeytsFCdFxdkKQ1mSiB8ipiMgk+T0qymIfoTCnNdBZ33eEqhkkj+CDw8bA5s2Y4REwcUOcRflUL0Zbqoztx8bb6oohsAX4fr49arao/qVyyHmJ4ZDTR7ADTO9dwgaCAJPt62shuUtU5OC0qazQJY6NfGBxozMjiG90nibhOJBpd7bImEyAcrFCYtC6R1aEGHmXN5kxy8V4Lc/s9u2esYwTypqVJSTMhtnuEX6VJzGgficpDROYAbwOeD2wHPlVlAaheJS1ja0DQuV4+vD3WF//0E+bHdoRptvuio7640e/EAeXQ2bPYesUrnfcJSOtEkhRY2nW5ui3HdajR2UEQsV3UOaHRJ5nea64KPyvTbtHgwDwUzRJtsRndSdqax/8FVuApjlcBH2mJRD1Glgkm3Lle9+8Pxn4n6fNgzSNK0KkVoYj9Om1bEbt60prMvLkNZ7flaIcal0zwhi2jnLd8aGodYN7cBgMNp1yiABw2x1vLydNmSR190OHmXYuJZs0t6i6blGzR9RhlJ8I02k+a2epEVV0CICL/gEWVV0JWxHR4hJY3KVyw35U375iyhUc9m/JSxH6dlgK9iBxRl9UgQaHrDCFutpM0O7j9nt1TXmNZ5rcoe/ZNMDwymtlm0TxfjX6ZtlYVzrQLxCY+jCN8/zTrLmseU0aUNOUxtfKmqvubqThlJJPUsQQulpds3MqGTTtZs2ox/QnrImnpqMteBC1iv05a4I8L+IsjzVziElE+b26DuYfMSjWZuMwOisTXXHbjds5bPjTD/NXoE/Y9vZ+Fa2+ZZmobG5+g0SfMm9uYqmIYltf1uqOu0s12/uYxZURJUx4ni8hj/msBBvz3glfK/BmVS9cDxHWsjT5h79P7p8UUXHbjdlYeP487fzaz+m84JXfVFLFfZ3kRxXlxhdOyJ42YXTrzgUY/V5yVPdNymVGldZTR2UJAMHsJp4APsggHs8HoXhOTytxDZjHyN9PXkKIzlKQcWlG502TPk03XPKaMMInKQ1Wbd3Q3MonrjPeFOpaA8YkD3P/rcd6wcsG0XEqvP+XYGUGEruRZAG22pG6aF1GSWeWSjVuZ4yfmi9svbznZNJISAz669ymnCoRrVi1OdO3dNTY+bQZ46vrbEtO2h/cJkxRJ/4aVC2bMapJSwjTT+bfbY6pVJZ0Nd0RLLqLSDlasWKGbN29utxilsWjtLYkLv2XURID4aO+BRn/sAnae7yaRdk0Cqa7KSfscMdCI7YSTysmmkRbVHlwrxGfRDdoh6RhRedLaImmfpGPPm9vgirNOyuxY437DvHU22tWBl3H/GfGIyBZVXVFkX5cgQaPFpLlklpUXKI8NvIzF0rRrUvJXggtMP1GKepKlzWKSasBHO1DX0XmRantJ8gUz1CxlmZYXy/Wealc2W1usryfufodGy8hyySyjFnoeG3gZi6VJbsNFCGqpx60xHDKrjw2bduZ2R80y38TVgI+6nK5eNsR5y4emHBiCFO3RDi6tLebNbcSOqNPku3jj1qlrTXPHDWQfGhxINCHWEVusryemPGpIOM9QEs0+OGnxBM18N5USLKRBWpWkqnx7nz5QKBYhS2G7FnKKS9EePf/qZUMcekj8pH9uQiXHrNnU6Ng4a76yjTVf3ZZ5/Xk743anUy/t/jNKxZRHTQmPEuNo9sHJU86zjNKfSWVv83LAT6viev2uI+pg1hCTW5CBRj+nnzA/swN1rQMCBys2RmlmUDAxqTNmY3Hnz9MZNxscWAZWeraemPKoOVU9OHmyqJaRcdWlUwy8h9JmXIFJKE+0tcu5g1lDVL8NDjSm4jTKHNG7dODhEf+l12/LvIYkoufPc0/lUYhVUWXGX6M4tmBec6rMC5RnAbTZxdKkReJ+ESZVZ1zXwrW3xB4nMAnFtcvep/bHel+5zFLS0sXffs9upwXbPO6wWYvrUQ+jvA4FaefPc0/VZb2hXYv1RjKmPDqAbnhwkjrLpBHkUEpMRUC0XYZHRlnzlW3TzGOu3ldl5OxKusbA5BXXUSd14EWi2Rt9AhFHgqQZhes9ZcGBRhKmPIyWEESZhwMc4zyRAgoHpUXXLBwdvLI6SdfcVINzG8ye1TeVnn7hMwemRdZH3WKTrt9lZN/oFw49ZNa0VPhQ7iy13cGBRn0x5dEB1C26tog8SZ5I4fQkYYqY6zZs2jljwdi1GFJWJ5m0La6+SVBDHbzcZGmR9UnkNfOFKfPeKGo2jaZSCcrl1uH+NcrBlEeNyJvbqR0PYNHsrEUCvfKa65qxzyd1kmHZg8SU0WqNaQvKSSsVWTLlNfNVSdLvkDSIiEulEtDu+9coD1MeNSGpU549q69W0bVFo31bsfDarH0+bg0lumgdTY9e9LoG5zamvY/riMPJFOs2Yk8bRGSt11h0eHdgyqMmJHXKdSstWrSzbMXCa9n2eRdFWWStBCDsPJXUEV917pLcObpaRVrbuNyb3RAdXjdzcquxOI+akPdhape3yxEDjVyfB7Qi0KvseAAXRZl2XWnXNjY+MS2Go92xFHnJKpebRad7a9UheLLdmPKoCUkPk1f6tD7RtUl1p7JqhbnmfWqWtNxTeXEJ5EtTWKuXDU0V9YoiMNXxJMVw1Hl03ky53G7w1qpD8GS7qZ3ZSkTWAf8D2O1/9F5V/Ub7JGoNSSaXK846CagmSLAISTmlkj4PyOttVRV5TA2uZrC4BeXgPEl1O1xC/uo8Ok9rm6jzQTd6W9UleLKd1E55+Fytqh9ptxCtJMslsi4PW9G1izqk1c7rKdaMm6pLedw0iozOW1ncy+V+rcs9WwUWPFnDYlD+zOOJPMqj24pB1ZmihXmyCly1YlblWqypqvNk4RLDkUSri3v1Ot3Sht1YDOodIvKnwGbgUlXdE/2CiFwEXASwYMGCFovXuxQdjSeN1ALbP1QfA9CsqcF1tO5yvHAxJmi+42l1ca9ep8qcc51CW5SHiHwLeHbMpr8GPgF8AO/Z+gDwUeAt0S+q6jXANeDNPCoT1phBEZNEnI082oFCtZ1YM6aGPCYvl0qB5y0f4vZ7dpfW8bS6uJfR/aa5LNqiPFT15S7fE5FPA1+vWByjZNJG6OHPkzrYqjqxZuJA8ozW0xRlWTXoo+RRjGavN8qgdmYrEXmOqv7Cf3sO8KN2ymPkI2uEHu40k9YGqurEmjE15Bmtt8OkkSejryU7NMqgjgvmXwCW4g3U7gf+PKRMYrEF8/qQZ1E6btGx6hF6XoJZVNIsyWWxvVWRyJcPb5+WtXjl8fP44QO/iV3Uhd621xseXbVgrqp/0m4ZjOIUHaGPjo1PWwMpY/G82U47y+XWZbReNJFkXuLiaL73s0cT15SaDaA0DIswN0olT31smF6rPamjK8Llw9u5ZOPWptJHpCX4c0190qpI5LjzJNkURsfGeyqNhlENpjyMUimaw6rMxfPhkdFpBZgC8nbaSecWcB65t8qzKe/xei0Pk1E+pjyMUimSnHB4ZDSx4F+RxfMNm3YWrqPhcu48MpVxjGbOk9SuvZaHySgfUx5G6eRNTnjlzTtiO3uBQh5AaQoiT6ddRibgVmQTTjvPhSuTA2gtrsNohtotmBudT94cS3sSkioqxRaV06LZ83TaZbjcxh3j9BPms2HTTi7ZuLU0T6c0WW+/Z3fTLtG9XrvCmEntXHWLYK669SFvzp+0PFBFc04luQBfuHIBH1y9JPfxyqQdOZGaPWe35HEyZtKMq66ZrYxSyetdlGY6KWraiVt3ufqCpW1XHNCeOhDNFsmy2hVGHGa2Mkolr3dRkolpcKDR1Kg2HM0emFzKNBMVpV15pZrJw2S5sIw4bOZhlEpe76Kkhd51Z59Uijx1KxfaKu+rMulEmY3qMeVhlEpe76Ky645HqZvJpVXeV2XSiTIb1WNmK6NUingoVZnaum4ml06sA9GJMhvVY95WRleTN1GjdZBGL2HeVoaRgKvJpW5rI4ZRd0x5GF2N65pK3dZGDKPu2JqH0fW4rKnUbW3EMOqOzTwMA3NHNYy8mPIwDLrTHXV4ZJRT19/GorW3cOr622z9xigVM1sZBt3njtqqCoZG72LKwzB8qow3aTVpDgDdco1GezGzlWF0IeYAYFSNKQ/D6ELMAcCoGlMehtGFdKMDgFEvbM3D6Gm6NSVJOx0AurVNjemY8jB6lm73SGqHA0C3t6lxEDNbGal0c6yApSQpH2vT3sFmHkYi3T6KNI+k8rE27R1s5mEk0u2jSPNIKh9r097BlIeRSLePIs0jqXysTXsHM1sZiRwzOBBbSKlbRpHdlpKkDlib9g5WSdBIJLrmAd4osswa44ZhtI+OqyQoIn8oIjtEZFJEVkS2XSYiPxWRnSKyqh3yGR6uhZQMw+g92mW2+hFwLvCp8IciciLwR8BJwDHAt0Tkt1T1wMxDGK2gm5IFGoZRHm2ZeajqT1Q1zmXntcCXVfUpVb0P+CnwktZKZxiGYWRRN2+rIeDB0PuH/M8MwzCMGlGZ2UpEvgU8O2bTX6vq15J2i/ksdkVfRC4CLgJYsGBBIRkNwzCMYlSmPFT15QV2ewg4NvT+ucCuhONfA1wDnrdVgXMZhmEYBamb2eom4I9EZLaILAJeAPxHm2UyDMMwIrQlzkNEzgH+DzAfGAO2quoqf9tfA28B9gMXq+o/OxzvcaA7cmY0z1HAI+0WoiZYWxzE2uIg1hYHWayqhxfZsSuCBEVkc9FAl27D2uIg1hYHsbY4iLXFQZppi7qZrQzDMIwOwJSHYRiGkZtuUR7XtFuAGmFtcRBri4NYWxzE2uIghduiK9Y8DMMwjNbSLTMPwzAMo4WY8jAMwzBy01HKQ0TO9FO1/1RE1sZsFxH5e3/73SLyonbI2Qoc2uJCvw3uFpHvicjJ7ZCzFWS1Reh7LxaRAyLyulbK10pc2kJEThORrX5ZhO+0WsZW4fCMHCEiN4vINr8t3twOOatGRD4rIg+LyI8SthfrN1W1I/6AfuBnwPHAIcA24MTId14N/DNejqyVwL+3W+42tsXvAPP816/q5bYIfe824BvA69otdxvvi0Hgx8AC//3R7Za7jW3xXuDD/uv5wKPAIe2WvYK2+D3gRcCPErYX6jc7aebxEuCnqnqvqj4NfBkvhXuY1wKfV4+7gEEReU6rBW0BmW2hqt9T1T3+27vw8oR1Iy73BcBfAjcAD7dSuBbj0hZ/DNyoqg8AqGq3todLWyhwuIgIcBie8tjfWjGrR1W/i3dtSRTqNztJebika++VlO55r/OteCOLbiSzLURkCDgH+GQL5WoHLvfFbwHzROQOEdkiIn/aMulai0tbfAz4b3jJV7cD71TVydaIVysK9ZvtqiRYBJd07c4p3TucPKnrT8dTHv+9Uonah0tb/B3wHlU94A0yuxaXtpgFLAd+HxgAvi8id6nqf1YtXItxaYtVwFbgDOB5wL+IyL+q6mNVC1czCvWbnaQ8XNK1O6d073CcrlNEXgh8BniVqv66RbK1Gpe2WAF82VccRwGvFpH9qjrcGhFbhusz8oiq7gX2ish3gZOBblMeLm3xZmC9eob/n4rIfcAJ9F4m70L9ZieZrX4AvEBEFonIIXi1zm+KfOcm4E9974GVwG9U9RetFrQFZLaFiCwAbgT+pAtHlWEy20JVF6nqQlVdCHwV+IsuVBzg9ox8DfhdEZklInOBU4CftFjOVuDSFg/gzcAQkWcBi4F7WyplPSjUb3bMzENV94vIO4BNeJ4Un1XVHSLyNn/7J/E8aV6NV/t8H97IoutwbIu/AZ4JfNwfce/XLswk6tgWPYFLW6jqT0Tkm8DdwCTwGVWNdeHsZBzviw8AnxOR7Ximm/eoatelaheR64DTgKNE5CHgCqABzfWblp7EMAzDyE0nma0MwzCMmmDKwzAMw8iNKQ/DMAwjN6Y8DMMwjNyY8jAMwzByY8rD6Hr8TLpbQ38LRWSFiPy9v/00Efmd0PdXi8iJBc7zREnylnIcw6iSjonzMIwmGFfVpZHP7gc2+69PA54Avue/Xw18HS/7rGEYMdjMw+hJ/NnG10VkIfA24BJ/VvIy4Gxgg//+ef7fN/1Egv8qIif4x1gkIt8XkR+IyAcSzvNhEfmL0Pt1InKpiBwmIt8WkR+KyHYRmZEJOJAx9P5jIvIm//VyEfmOL9OmIAuqiPxPEfmxX5fhy6U1mGFEsJmH0QsMiMhW//V9qnpOsEFV7xeRTwJPqOpHAETkJuDrqvpV//23gbep6n+JyCnAx/GS6f1v4BOq+nkReXvCub+Ml5jx4/7784EzgSeBc1T1MRE5CrhLRG5Sh6hdEWkA/wd4raruFpELgA8BbwHWAotU9SkRGXRtIMPIiykPoxeIM1s5ISKH4RXW+kooI+9s//+pwHn+6y8AH47ur6ojInK0iByDV3Boj6o+4CuA/yUiv4eXJmQIeBbwSwexFgO/jZcFFrz0G0EuoruBL4rIMNCN+buMmmDKwzDS6QPGUpSPS36frwKvA56NNxMBuBBPmSxX1QkRuR+YE9lvP9NNy8F2AXao6ktjzvUavMpxZwPvE5GTVLXrChwZ7cfWPAwDHgcOj3vv13a4T0T+EKbqPQf14O/Ey9YKnjJI4sv+916Hp0gAjgAe9hXH6cBxMfv9HDhRRGaLyBH4GWCBncB8EXmpL1NDRE4SkT7gWFW9HfgrvJKzhzm1gGHkxJSHYcDNwDn+Avnv4nX2a0RkRESeh6cY3ioi24AdHCxn+k7g7SLyAzxlEIuq7sBTRqOhVNdfBFaIyGb/+PfE7PcgcD2+KQoY8T9/Gk8RfdiXaSueaa0fuNbPEjsCXK2qY0UbxTDSsKy6hmEYRm5s5mEYhmHkxpSHYRiGkRtTHoZhGEZuTHkYhmEYuTHlYRiGYeTGlIdhGIaRG1MehmEYRm7+HxrVO0V94Pj2AAAAAElFTkSuQmCC\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZFElEQVR4nO3dfZRcdX3H8feHEEzkGbNgIAlbEVGgNdQtYkMt5UEDqIAWFRVDxYaqVPDQalDbhtPWxlMF7dGDJxZKeBBMBYSCFmOQpihCNxpiYqIgBBKIyfIQk4BgE7794/62XiYzO3d3Z3b2l/28ztmzcx/m3u+9M/OZ3/3dOzOKCMzMLD+7dLoAMzMbGge4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOCJpJWSjut0HZ0k6QxJayVtlXRUp+vpJ2mupGtbuLxzJN1dGt4q6RWtWn5a5l2SPlhx3uMkrWvBOt8r6TvDXc5o1GzbBrO/m6ynJY/FSBkTAS5pjaQTa8a96EUcEUdExF1NltMtKSTt2qZSO+1zwPkRsUdE/LjqnSRdJekf2lhXW6XtfajTdQxXRFwXEW/qdB3tsDNv23CMiQDPxSh4YzgYWNnhGlpK0rhO12CFUfD83uk4wJNyK13S0ZJ6JW2WtEHSpWm2Jen/pnTY/QZJu0j6tKRHJG2UdLWkvUvLfX+a9qSkv6lZz1xJ35B0raTNwDlp3fdI2iRpvaQvSdqttLyQ9GFJD0jaIunvJR2S7rNZ0sLy/DXbWLdWSS+RtBUYB9wv6Rd17itJl6X7/UrScklHSpoNvBf4eNon/5HmnyPpF6nGn0o6o7SscyTdLelzkp6W9LCkk0vTf0fSf6X7LgIm1dTy75J+mepYIumI0rSrJF0u6VuSngH+RNLLJN2a9s99wCE1ywtJr5R0YNqG/r9nJUVpvg9IWpVqvkPSwaVpJ0lanWr6EqB6j0Gad2Kq82lJPwX+oGb6gZJulNSX9s1HS+N/LWm/0rxHSXpC0njt2DX0RRVdYpslLZX0R6Vpc9Nz5eq0n1dK6ilNnyrpplTDk2mbmu6Hmu3oP2I9V9KjwJ1p/DGSfpCe4/er1HWZtuGhVNPDkt5bGl/etob7WzVdbqo5cpb0Z6n+LWld5w3wWH1C0mNp3p9JOqHRvB0RETv9H7AGOLFm3DnA3fXmAe4Bzk639wCOSbe7gQB2Ld3vA8CDwCvSvDcB16RphwNbgWOB3Si6KP63tJ65afh0ijfTicDrgGOAXdP6VgEXltYXwK3AXsARwPPA4rT+vYGfArMa7IeGtZaW/coG930zsBTYh+LF8hpgcpp2FfAPNfOfCRyYtutdwDOl+c9J2/3nFG8aHwIeB1Ta/5cCLwHeCGwBrq3Zjj3T9C8Ay0rTrgJ+BcxI654A3AAsBHYHjgQeq3ns6243cB1wfbp9etp3r0mPzaeBH6Rpk4DNwJ8C44GPAduADzbYl/OA/wb2A6YCK4B1adouaT//bXrOvAJ4CHhzmn4n8OelZf0z8JUGz+n3AS9L9V4E/BKYUHruPQeckh6DfwJ+mKaNA+4HLkv7bAJwbLP9UGc7u9O+vTotZyJwEPBkWu8uwElpuCvNsxk4LN1/MnBE7bY1299p266tU8euafhUijdxAX8MPAv8fpp2XOmxOAxYCxxYWs4hnc6zF+3jThcwIhtZhPNWYFPp71kaB/gS4BJgUoMnZDnAFwMfLg0fRhFOu1K8CK8vTXsp8BteHOBLmtR+IXBzaTiAGaXhpcAnSsOfB77QYFkNay0tu1GAHw/8nOLNZZeaaVdRE+B17r8MOC3dPgd4sGa/BPByYFp6Me5emv618guyZrn7pPvuXarl6tL0cWkbX10a9xmaBDjwibRvJ6bhbwPnlqbvkp5DBwPvJ4VfmiZgHY0D/CFgZml4Nr8NjdcDj9bMfzHwb+n2B4E7S+tZC7yxtF/vrrfONP1p4LWl5953S9MOB36dbr8B6KP0PC/N13A/1Jm3O+3bV9Ts12tq5rsDmEUR4JuAd/Tv99I8/79tzfY3TQK8Tp3fBC5It48rPRavBDYCJwLjB3p+d+pvLHWhnB4R+/T/AR8eYN5zgVcBqyX9j6S3DDDvgcAjpeFHKML7gDRtbf+EiHiWorVRtrY8IOlVkm5LXQSbKcJmUs19NpRu/7rO8B5DqHVAEXEn8CXgy8AGSfMl7dVofhVdR8vSYfImipZveTt+WVr2s+nmHqnGpyPimZo6+5c7TtI8Fd0zmyneeKlZdnmfdqVtLI8r74N6tZ8MXEDxnPl1Gn0w8MXS9jxFERwHsePjHDXrq3Ugjes5GDiwfz1pXZ/kt4/RN4A3SDqQ4ugkKFrz9bbjotRV8Ku0nL1p8BhQhPCE1M0wFXgkIrbVWexA+6GR8rYeDJxZs33HUhydPUNxtPYXwHpJt0t6dZ3lDXZ/v4ikkyX9UNJTaf2nsONrjIh4kKIBNRfYKOmGtN9HjbEU4JVFxAMRcRawP/BZ4BuSdqd4sdR6nOJJ2a+/BbkBWA9M6Z8gaSLFIe2LVlczfDmwGjg0IvaiePE27E8dpIFqbSoi/iUiXkfRdfMq4K/7J5XnS32iXwXOB16W3jBXUG071gP7pv1drrPfe4DTKFpFe1O0rqhZdrmePoptnNpgeS8i6TBgAfDOiCiHwlrgvHIjICImRsQPUs1TS8tQzfrqbWOjetYCD9esZ8+IOAUgIjYB3wHeSbEvrk8BVrsdf0TR2n0nsG96DH5FtcdgLTBN9U86DrQfGinXt5aiBV6+/+4RMS9t3x0RcRJF98lqiudRrWb7+xmKo7p+Ly/N+xLgRoruzAPSfvkWDfZLRHwtIo6leN0ERR6MGg7wOiS9T1JXRLxAcUgHsJ0iDF6g6Jfsdz3wMRUn3vagaDF/PbVevgG8VdIfqjixeAnNX0B7UvTvbU2tjw+1bMMGrnVAkv5A0usljad4gTxHsU+geAMo75P+N7u+dN8/o2iBNxURjwC9wCWSdpN0LPDW0ix7UvT7P0nxIv1Mk+Vtp+jrnyvppZIOpzhcr7eNewG3AJ+OiLtrJn8FuFjphKmKk79npmm3A0dIensKvY9SCo06FqZl7StpCvCXpWn3AZvTybOJ6YjjSEnlE51fo+hGeEe6Xc+eFG9cfcCukv6W4rxJFfdRhOQ8SbtLmiBpRpo20H6o4lqK18Sb07ZNUHHt9RRJB0h6W3rzfp6i23N7nWU029/LgDdKmqbigoKLS9N2ozh30gdsS0dbdS9PlHSYpONT6D9HcXRbr56OcYDXNxNYqeLKjC8C746I59Kh/j8C30+Hf8cAVwLXUPSbP0zxQP8lQESsTLdvoHhBbKHoU3t+gHX/FUXLagtF6+PrLdyuhrVWsFeq52mKQ/4nKVoxAFcAh6d98s2I+ClFX/w9FOH+u8D3B1Hneyj6gp8C/o7iJFi/q9P6H6M4YfvDCss7n6J75pcUfeT/1mC+36c4L3CpSlejAETEzRStrxtS180K4OQ07QmKk7bzKPbLoQy8vZekbXiYojV9Tf+E9IbzVmB6mv4E8K8URxv9bk3r2BAR9zdYxx0U/dU/T+t6jordDKUaXgk8StG//K40reF+qLjstRRHUJ+kCNG1FEdyu6S/iyiOFJ+iOMG4Q1dns/0dEYsoXjfLKc5j3FaatoUi8BdSPJffQ7E/63lJWscTFM+d/VPdo0b/WX8bAanVu4mie+ThTtdjZnlzC7zNJL01HbrvTtFi/Qm/PfFmZjZkDvD2O43ikPBxikO9d9c76WRmNljuQjEzy5Rb4GZmmRrRL5eZNGlSdHd3j+Qqzcyyt3Tp0icioqt2/IgGeHd3N729vSO5SjOz7Emq++lhd6GYmWXKAW5mlikHuJlZphzgZmaZcoCbmWXKAW5mlikHuJlZphzgZmaZcoCbmWVqRD+JaWNP95zbBzX/mnmntqkSs52PW+BmZpmqHODp9+t+LOm2NLyfpEWSHkj/921fmWZmVmswLfALgFWl4TnA4og4FFichs3MbIRUCvD0y9mnUvy4ar/TgAXp9gLg9NaWZmZmA6naAv8C8HHghdK4AyJiPUD6v3+9O0qaLalXUm9fX9+wijUzs99qGuCS3gJsjIilQ1lBRMyPiJ6I6Onq2uH7yM3MbIiqXEY4A3ibpFOACcBekq4FNkiaHBHrJU0GNrazUDMze7GmLfCIuDgipkREN/Bu4M6IeB9wKzArzTYLuKVtVZqZ2Q6Gcx34POAkSQ8AJ6VhMzMbIYP6JGZE3AXclW4/CZzQ+pLMzKwKfxLTzCxTDnAzs0w5wM3MMuUANzPLlAPczCxTDnAzs0w5wM3MMuUANzPLlAPczCxTDnAzs0w5wM3MMuUANzPLlAPczCxTDnAzs0w5wM3MMuUANzPLVJUfNZ4g6T5J90taKemSNH6upMckLUt/p7S/XDMz61flF3meB46PiK2SxgN3S/p2mnZZRHyufeWZmVkjTQM8IgLYmgbHp79oZ1FmZtZcpT5wSeMkLQM2Aosi4t406XxJyyVdKWnfBvedLalXUm9fX1+LyjYzs0oBHhHbI2I6MAU4WtKRwOXAIcB0YD3w+Qb3nR8RPRHR09XV1aKyzcxsUFehRMQmil+lnxkRG1KwvwB8FTi6DfWZmVkDVa5C6ZK0T7o9ETgRWC1pcmm2M4AV7SnRzMzqqXIVymRggaRxFIG/MCJuk3SNpOkUJzTXAOe1r0wzM6tV5SqU5cBRdcaf3ZaKzMysEn8S08wsUw5wM7NMOcDNzDLlADczy1SVq1DM/l/3nNs7XYKZJW6Bm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llyp/EHMP8qUqzvLkFbmaWqSo/qTZB0n2S7pe0UtIlafx+khZJeiD9r/ur9GZm1h5VWuDPA8dHxGspfoF+pqRjgDnA4og4FFichs3MbIQ0DfAobE2D49NfAKcBC9L4BcDpbanQzMzqqtQHLmmcpGXARmBRRNwLHBAR6wHS//0b3He2pF5JvX19fa2q28xszKsU4BGxPSKmA1OAoyUdWXUFETE/Inoioqerq2uodZqZWY1BXYUSEZuAu4CZwAZJkwHS/40tr87MzBqqchVKl6R90u2JwInAauBWYFaabRZwS7uKNDOzHVX5IM9kYIGkcRSBvzAibpN0D7BQ0rnAo8CZbazTzMxqNA3wiFgOHFVn/JPACe0oyszMmvMnMc3MMuUANzPLlAPczCxTDnAzs0z562RtVBnsV9yumXdqmyoxG/3cAjczy5QD3MwsUw5wM7NMOcDNzDLlk5g7Ef/GpdnY4ha4mVmmHOBmZplygJuZZcoBbmaWKZ/EtKz5k5s2lrkFbmaWqSo/qTZV0vckrZK0UtIFafxcSY9JWpb+Tml/uWZm1q9KF8o24KKI+JGkPYGlkhalaZdFxOfaV56ZmTVS5SfV1gPr0+0tklYBB7W7MDMzG9ig+sAldVP8Pua9adT5kpZLulLSvi2uzczMBlA5wCXtAdwIXBgRm4HLgUOA6RQt9M83uN9sSb2Sevv6+lpQspmZQcUAlzSeIryvi4ibACJiQ0Rsj4gXgK8CR9e7b0TMj4ieiOjp6upqVd1mZmNelatQBFwBrIqIS0vjJ5dmOwNY0fryzMyskSpXocwAzgZ+ImlZGvdJ4CxJ04EA1gDntaVCMzOrq8pVKHcDqjPpW60vx8zMqvInMc3MMuUANzPLlAPczCxTDnAzs0w5wM3MMuXvA7cxxd8fbjsTt8DNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsU1V+E3OqpO9JWiVppaQL0vj9JC2S9ED6v2/7yzUzs35VWuDbgIsi4jXAMcBHJB0OzAEWR8ShwOI0bGZmI6RpgEfE+oj4Ubq9BVgFHAScBixIsy0ATm9XkWZmtqNBfZ2spG7gKOBe4ICIWA9FyEvav8F9ZgOzAaZNmzacWsecwX71qZmNLZVPYkraA7gRuDAiNle9X0TMj4ieiOjp6uoaSo1mZlZHpQCXNJ4ivK+LiJvS6A2SJqfpk4GN7SnRzMzqqXIVioArgFURcWlp0q3ArHR7FnBL68szM7NGqvSBzwDOBn4iaVka90lgHrBQ0rnAo8CZ7SnRzMzqaRrgEXE3oAaTT2htOWZmVpU/iWlmlikHuJlZphzgZmaZcoCbmWXKAW5mlikHuJlZphzgZmaZcoCbmWXKAW5mlqlBfZ2s2VgzlK/0XTPv1DZUYrYjt8DNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy1SVn1S7UtJGSStK4+ZKekzSsvR3SnvLNDOzWlVa4FcBM+uMvywipqe/b7W2LDMza6ZpgEfEEuCpEajFzMwGYTh94OdLWp66WPZtNJOk2ZJ6JfX29fUNY3VmZlY21AC/HDgEmA6sBz7faMaImB8RPRHR09XVNcTVmZlZrSEFeERsiIjtEfEC8FXg6NaWZWZmzQwpwCVNLg2eAaxoNK+ZmbVH028jlHQ9cBwwSdI64O+A4yRNBwJYA5zXxhrNzKyOpgEeEWfVGX1FG2oxM7NB8Ccxzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8tU06+Ttca659ze6RJsFBrs82LNvFPbVInt7NwCNzPLVNMAT786v1HSitK4/SQtkvRA+t/wV+nNzKw9qrTArwJm1oybAyyOiEOBxWnYzMxGUNMAj4glwFM1o08DFqTbC4DTW1yXmZk1MdQ+8AMiYj1A+r9/oxklzZbUK6m3r69viKszM7NabT+JGRHzI6InInq6urravTozszFjqAG+QdJkgPR/Y+tKMjOzKoYa4LcCs9LtWcAtrSnHzMyqqnIZ4fXAPcBhktZJOheYB5wk6QHgpDRsZmYjqOknMSPirAaTTmhxLWZmNgj+JKaZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpappj/oMBBJa4AtwHZgW0T0tKIoMzNrblgBnvxJRDzRguWYmdkguAvFzCxTww3wAL4jaamk2a0oyMzMqhluF8qMiHhc0v7AIkmrI2JJeYYU7LMBpk2bNszVmZlZv2G1wCPi8fR/I3AzcHSdeeZHRE9E9HR1dQ1ndWZmVjLkAJe0u6Q9+28DbwJWtKowMzMb2HC6UA4AbpbUv5yvRcR/tqQqMzNrasgBHhEPAa9tYS1mZjYIrbgOfKfRPef2TpdgY9Bgn3dr5p3apkosN74O3MwsUw5wM7NMOcDNzDLlADczy1Q2JzF9osdsaPza2Xm5BW5mlikHuJlZphzgZmaZcoCbmWUqm5OYg+VPVdrOarQ9t4dSj0+UtoZb4GZmmXKAm5llygFuZpYpB7iZWaZ22pOYZjY0I3GStN2fDh1tyx/KOqpwC9zMLFPDCnBJMyX9TNKDkua0qigzM2tuOD9qPA74MnAycDhwlqTDW1WYmZkNbDgt8KOBByPioYj4DXADcFpryjIzs2aGcxLzIGBtaXgd8PramSTNBmanwa2SfjaMdbbTJOCJThcxSK55ZORWc271QpOa9dn2rnyIyx/Ufh7mNhxcb+RwAlx1xsUOIyLmA/OHsZ4RIak3Ino6XcdguOaRkVvNudULrnmohtOFsg6YWhqeAjw+vHLMzKyq4QT4/wCHSvodSbsB7wZubU1ZZmbWzJC7UCJim6TzgTuAccCVEbGyZZWNvFHfzVOHax4ZudWcW73gmodEETt0W5uZWQb8SUwzs0w5wM3MMuUAL5H0z5JWS1ou6WZJ+3S6pmYknSlppaQXJI3ay7By+9oFSVdK2ihpRadrqUrSVEnfk7QqPScu6HRNA5E0QdJ9ku5P9V7S6ZqqkjRO0o8l3dbJOhzgL7YIODIifg/4OXBxh+upYgXwdmBJpwtpJNOvXbgKmNnpIgZpG3BRRLwGOAb4yCjfz88Dx0fEa4HpwExJx3S4pqouAFZ1uggHeElEfCcitqXBH1Jc2z6qRcSqiBitn27tl93XLkTEEuCpTtcxGBGxPiJ+lG5voQiYgzpbVWNR2JoGx6e/UX9VhaQpwKnAv3a6Fgd4Yx8Avt3pInYS9b52YdQGy85AUjdwFHBvZysZWOqKWAZsBBZFxKiuN/kC8HHghU4XMuZ+0EHSd4GX15n0qYi4Jc3zKYrD0etGsrZGqtQ8ylX62gVrDUl7ADcCF0bE5k7XM5CI2A5MT+ebbpZ0ZESM2vMOkt4CbIyIpZKO63Q9Yy7AI+LEgaZLmgW8BTghRslF8s1qzoC/dmGESBpPEd7XRcRNna6nqojYJOkuivMOozbAgRnA2ySdAkwA9pJ0bUS8rxPFuAulRNJM4BPA2yLi2U7XsxPx1y6MAEkCrgBWRcSlna6nGUld/Vd6SZoInAis7mxVA4uIiyNiSkR0UzyP7+xUeIMDvNaXgD2BRZKWSfpKpwtqRtIZktYBbwBul3RHp2uqlU4M93/twipg4Wj/2gVJ1wP3AIdJWifp3E7XVMEM4Gzg+PT8XZZaiqPVZOB7kpZTvMkvioiOXpaXG3+U3swsU26Bm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWab+D4QIsf2F8nYNAAAAAElFTkSuQmCC\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZdrH8e9NQBBBkYDKqgFELOAKrtjLim1tq6IiIkpzRYgIFpSSd9VVQ1UUGxClCIwIihXBhnXtoKggLoISQCwQRFDUxOR5/zgTHDIzyUnImZkkv8915Zo5Z86cuU+UufOUcz/mnENERCRSrWQHICIiqUfJQUREoig5iIhIFCUHERGJouQgIiJRaic7gMrQpEkT16JFi2SHISJSpSxatGiDc65prNeqRXJo0aIFCxcuTHYYIiJVipnlxntN3UoiIhJFyUFERKIoOYiISBQlBxERiaLkICIiUZQcRESqqVAIWrSAWrW8x1DI/3urxVRWERHZXigEffrA1q3edm6utw3QrVvZ71fLQUSkGsrK+jMxFNu61dvvR1KTg5lNNrMfzGxJxL5bzewbM1sc/jkrmTGKiFRFq1eXb39JyW45TAXOiLH/budc+/DPvATHJCJS5WVklG9/SUlNDs65N4GNyYxBRKQ6ys6G+vW331e/vrffj2S3HOLpb2afhruddo91gJn1MbOFZrZw/fr1iY5PRCSldesGOTnQvDmYeY85Of4GowEs2WtIm1kLYK5z7pDw9p7ABsABtwPNnHO9SztHhw4dnArviYiUj5ktcs51iPVayrUcnHPfO+cKnXNFwEPAkcmOSUSkpkm55GBmzSI2OwFL4h0rIiLBSOpNcGY2EzgJaGJma4FbgJPMrD1et9Iq4KqkBSgiUkMlNTk457rG2D0p4YGIiMh2Uq5bSUREkk/JQUREoig5iIhIFCUHERGJouQgIiJRlBxERCSKkoOIiERRchARkShKDiIiKW5H1oKuKK0hLSKSwnZ0LeiKUstBRCSF7eha0BWl5CAiksJ2dC3oilJyEBFJYTu6FnRFKTmIiKSwHV0LuqKUHEREUtiOrgVdUZqtJCKS4rp1Cz4ZlKSWg4iIRFFyEBGRKEoOIiISRclBRESiKDmIiEgUJQcREYmi5CAiIlGUHEREJIqSg4iIRFFyEBFJomQs5OOHymeIiCRJshby8UMtBxGRJEnWQj5+KDmIiAQsXtdRshby8UPdSiIiASqt6ygjw9suKeiFfPxIasvBzCab2Q9mtiRiX2Mze9nMvgw/7p7MGEVEdkRpXUfJWsjHj2R3K00FziixbwiwwDnXGlgQ3hYRqZJK6zpK1kI+fiQ1OTjn3gQ2lth9HvBI+PkjwPkJDUpEpBKVtQZ0t26wahUUFXmPqZAYIPkth1j2dM59CxB+3CPWQWbWx8wWmtnC9evXJzRAEZF4Sg4+n3VWEruOnIP58+G778r91lRMDr4453Kccx2ccx2aNm2a7HBERLYNPufmet/LubnwyCPQo0eCu46KimDOHOjQwctOOTnlPkUqzlb63syaOee+NbNmwA/JDkhExI94g8/z5nldRoErKICZM2HECPjiC2jdGiZPrlAmSsWWw7NAj/DzHsAzSYxFRMS3pN238NtvMGECHHCA10zZaSd47DFYtgx69fK2yynZU1lnAu8CB5rZWjO7AhgJnGZmXwKnhbdFRFJeWYPPle6XX2DsWNhvP+jXD/bcE559FhYvhi5dIC2twqdOareSc65rnJdOSWggIiKVIDt7+xveIKDB502b4L77YNw4yMuDk0+G6dO9R7NK+YhU7FYSEamSAr9v4YcfYOhQryly881wzDHw7ruwYAGcckqlJQZIzQFpEZEqq1u3AGYirVkDY8bAQw/B779D584wbBi0a1fJH/QnJQcRkVS1YgWMHAnTpnlzYy+/HIYM8QaeA6ZuJRGRHRDIYj1LlsCll8KBB8KMGd5AxooV3rTUBCQGUMtBRKRcQiHvfobVq6FxY9iyBfLzvdd2eLGeDz6A4cPhmWegQQMYNAiuuw722qvS4vdLLQcREZ9K3gGdl/dnYihW7sV6nIPXX4fTToOjjoI334Rbb/U+ZNSopCQGUMtBRMS3WHdAx+LrprfiukfZ2fDOO949CqNHQ9++0LDhDse6o5QcRER88nunc6k3vRUWwpNPet1Hixd7B99/P/TuDTvvXClxVgZ1K4mI+OTnTue4N70VFHhV+A45BC6+2GuCTJniDTRffXVKJQbwkRzMrJWZ1Q0/P8nMBphZo+BDExFJLbFWbqtTB9LTS7np7bffYPx4b5ZRz55Qty7MmgWff+5t16mTwCvwz0/LYQ5QaGb7A5OAlsCjgUYlIpJCiqerXn659wd+ZDKYMgU2bIixWM/PP8Odd0LLlpCZ6Q0sP/ccfPyx13LYgbpHieBnzKHIOfeHmXUC7nHO3WdmHwcdmIhIspQ2XTUvz2s9TJ8eZ7rqjz/+Wfdo40avrMWjj8JJJ1VqeYug+Wk5FJhZV7zy2XPD+1KzHSQisoMqPF31+++9u5ebN4dbboHjjoP33oNXXoGOHatUYgB/LYdeQF8g2zn3tZm1BGYEG5aISHKUe7rq6tVe3aOHH/ayyMUXe8XxDj000DiDVmZycM59bmaDgYzw9tdojQURqab8Tlc9sdmXcEW47hFA9+5ey6F16+CCSyA/s5X+CSwGXghvtzezZ4MOTEQkGcqarnoInzE7rSuvfnuQN5bQty+sXAmTJlWbxAD+xhxuBY4ENgE45xbjzVgSEanyIgvnNWnizTwqqU4dOG3X93mGc/mMQzm/zlxq3XSjNz3pvvsCXOotefyMOfzhnPvJth9McQHFIyKSMMWDz8VjDHl5JY9wnLfr64zfN5tmSxd4U5cG/oc611wDu++e6HATyk9yWGJmlwJpZtYaGAC8E2xYIiLBiz/47Dib5xnGcI7d/C7k7eXds3DVVV611BrAT7fSNUBb4HdgJrAZuDbIoEREghYKedNVI9WikM7M5mMOYy7/5C+sI5MH4euv4YYbakxiAH+zlbYCWeEfEZEqr7g7qVhtCuhGiCGM5CD+xxccSA+m8iiXsnfzOlAvebEmS9zkYGbPUcrYgnPu3EAiEhEJUCgEPXp4xVHr8Su9mcxNjKY5q/mY9nRmNk9yAUWkxS+iVwOU1nK4M2FRiIgkQHGLYefCLfRlAjdwF3vxPW9zLP0Yz3zOJD3dcBuheYaXGCq0ols1EDc5OOfeSGQgIiJBGz1kI4O23sdAxtGYH3mZU7mEx3iDvwNG8+be7FQpvVtptnPuYjP7jBjdS865qn1vuIjUGHMe+I71w8by383jacjPPM15DGcYH3LktmNqchdSLKV1Kw0MP56TiEBERCpdbi6vnzOGs5ZMYifymUUXRjCUJfx1u8PS0mKsw1DDxZ3K6pz7Nvw00zmXG/kDZCYmPBGRCli+nJUn9qKgxf4cuySHEN04kP/RjUejEkP9+t4CbUoM2/Nzn8NpMfadWdmBiIhUVCjklb5oZ58wy7pQdOBBNHtrFg+SSStWciUPs5L9o94Xc+U2AUofc+iH10LYz8w+jXipIfB20IGJiJQlFIKBA2H/vPeYQjb/ZC6bacgoBnM317GePeK+V4PPpSttzOFRYD4wAhgSsX+Lc25joFEBZrYK2AIU4tV36hD0Z4pI1ZHZz7F8wqvMYjin8Cp5NObf3Mb99GcTpdc9MtPgc1lKm8r6E/AT0NXM0oA9w8c3MLMGzjmfVc93SEfnXIwaiSJSU4VmOOZlzuWaLdkczfusoxnXcxc59OEXyi5vYeZV2VZXUunKLJ9hZv3xynZ/DxSFdztAU1lFJHEKC5l05hN0eHk43fiUr2lBX8YzlZ787rO+RXq6t7SzEkPZ/FRlvRY40DkXVcw2YA54ycwcMNE5l5PgzxeRVJCfz7v9Z7DHpJFcUfQlyziI7jzCTLryh8/l7JUUys9PcliD172UaMc559aZ2R7Ay2b2hXPuzeIXzawP0AcgoxoutCFS4/36K0yaxMahoznm5zV8xGFcyBM8RSdcGRMta9WCoiJv0Lkml8DYEX6Sw1fA62b2PF7ZbgCcc2MDi8o7/7rw4w9m9hTeanRvRryeA+QAdOjQQYsPiVRxxTOP8vM205cJXM9Y9uJ7Puc4spnIC5wBWNz3F48lPPhg4mKuzvwkh9Xhn53CP4Ezs12AWs65LeHnpwO3JeKzRSSxQiFvDZ26v+QxgHsZwL3sziZe5HQuJou3OLHMc6jbqPL5Wc/hP4kIpIQ9gafCS5PWBh51zr2QhDhEJECZmfDU+G+5hbH0YzwN+IWnOJ/hDGMhR5T5frUWguNntlJT4Ca81eC2TQlwzp0cVFDOua+AdkGdX0SS7/8uz6XtjNGMZRJ1KOAxLmEEQ1nKIb7er9ZCsPx0K4WAWXgF+PoCPYD1QQYlItXXs2P+x+ZhI7jljxAOYyo9Gc1NMctbxNOvn1oLQfOTHNKdc5PMbGB4jYc3zExrPYhI+SxezKKLhnPOyif4jXo8wNXcySC+YZ9ynUaJITH8JIeC8OO3ZnY2sA7K+V9TRGqud97x5pPOm8f+7MpIhnAP15Za9ygWdSMllp/kcIeZ7QbcANwH7ApcF2hUIlKlhWY45mQu4Jot2XTkdTaQzt3cwQNczU80KvW9DRrAhAlKAsnmZ7bS3PDTn4COwYYjIlVaURGvD5pL63uyedJ9wDf8hesYSw592MouZb5dXUapw89spSnEXia0dyARiUjVU1jIfwfMZveJIzip8DO+oiV9mMgj9CCfur5OocSQWvx0K82NeF4P6IQ37iAiNV1+PjP+MZ2jXh/J8azgcw7mMqbzGJdQ6OvrRd1IqcpPt9KcyG0zmwm8ElhEIpLyHpvyK4v7P0zm1jFcxhoW8TcuYA5Pc36ZdY8iqbWQuvyl9u21BlTpTqQm2ryZh//2IP9ceTeX8ANvcTx9yOFF/kFpdY9iUWJIbX7GHLbgjTlY+PE7YHDAcYlICnl8Qh4rB4zjqoL7+BebeIF/kE0W/+WEcp9LU1KrBj/dSg0TEYiIpJZQCG65ch19fx1LXybQgF+YwwWMYCiLKN+qvWolVD2lJgcz2xnoBrQJ71oIPOGcyw86MBFJvOKy2Q3yVjGYUSxlMmkUMpOujGAoy7Z9FfhTq5ZXcVWJoeqJO3JkZn8FlgEnAKuAXOAfwNtm1sjM7khIhCISqMxMr7qpGdx+2RfcldeDFexPbyYzlZ4cwHK6M913YqhXD2bMAOegsFCJoaoqreVwL3Clc+7lyJ1mdiqwBFgaZGAiEpziNRR++cXbbs/HDGM4FzKH36jHfVzDnQxiHXv7PqfGEqqX0pJDs5KJAcA594qZFeDd7yAiVUgoBL17Q364Y/hY3iaLbM5iPj+xK8MZxjgGsoGmvs95yinwiia3VzulTUiuZWZRtzaaWT2gwDm3NbiwRKSyhUJw+eWQn+84lZd5jZN4m+M5gg8ZRjYZrObf3OE7MTRo4HUfKTFUT6Ulh2nAHDNrUbwj/Hw2MD3IoESkcoVC0P2yIs51T/M+R/Eyp7M/K7iWu2nBKkYwjM3sVuZ5ihOCc7Bli7qQqrO4ycE5dwfwAvCmmW0wsw3AG8DLzrnbExWgiFRcKAT1d/qD5y97lMW042k6kU4eV5JDK1YyjmvLLIhn5k1FVUKoWUqdyuqcux+438wahre3JCQqEdkhmZnw8Ph8ujONTxnJ/qxkKW3oxgxm0cVX3SMNMNdsvspnKCmIVB2HH7yV4794iJXcyb6s5UM6cD5P8Sznllr3SAXwJFJFaiuJSIoJheDGPj/RY+uDzOdu9mA9b3AiVzCJlzmN0uoeabaRxKLkIFLFXXDiBg57axyfcx+N+In5nEE2WbzN8WW+V2UtJB4/hffq4y0RmuGcu9LMWgMHRqwQJyJJMLTHOppOu5PpTGRnfuVJLmA4w/iYv/l6vxKDlMZPy2EKsAg4Jry9Fnic7RcBEpEEyMyE+eO/ZjCjuJUppFHIo1zKSIaUq7zFww9rbEFK52dVjlbOudFAAYBz7lfKW7hdRCosFPIGi9vY5xw9vjtf0ppeTGEKvTiA5fRgmq/EkJ7u3aPw669KDFI2Py2H/HB1VgdgZq2A3wONSkS2lbpom/8Rj5BNJ57iV3ZmHAMZy/W+6x61aQNLVQlNyslPcrgF72a4fc0sBBwH9AwyKJGaLhSCCZf9l6fJ5kxeYBO7kU0W4xhIHk18ncMM+vbVuIJUjJ/Ffl42s4+Ao/G6kwY65zYEHplIDZTZz7Fywktkkc1bvMUPNGUow3mQTF/lLUAtBakccZODmZWc8vBt+DHDzDKccx8FF5ZIzREKQd8+RZy69RmGMZwjWMga9mEA43iYf/Er9X2fSzOQpLKU1nK4q5TXHHByJcciUuP07/sHP06cxbuM4BCWsoJW/IuHmEZ3CtjJ1znUfSRBiJscnHMdExlISWZ2BjAOSAMeds6NTGY8IpUlMxMmjf+dHjzCYEbRiq9YQlsuJcRsLvZV9wjUSpBglTmV1czqmdn1Zvakmc0xs2vDazoExszSgAeAM/HWr+5qZuVbvFYkxYRCsGvaL+w0/h5W0oocriKPdM7jaQ7lU2ZyaZmJIXIJTiUGCZKfP1GmAVuA+8LbXfHWc+gcVFDAkcAK59xXAGb2GHAe8HmAnylSqSKX4tyVn7iaB1jJ3TRlA6/zd3oxhVc4Fb+3DakGkiSSn+RwoHOuXcT2a2b2SVABhe0NrInYXgscFXmAmfUB+gBkZGQEHI6If5FLcTZhPUO5h/7cz25s5nnOYjjDeIfjfJ9PYwqSDH7ukP7YzI4u3jCzo4C3gwvJ+5gY+9x2G87lOOc6OOc6NG3qf71bkcpWfAezmfdz2WXQJP8bxnIdq2jBUEbwEqdzGB9xDs/7TgzFC+wUFSkxSOL5aTkcBXQ3s9Xh7QxgmZl9Bjjn3KEBxLUW2Ddiex9gXQCfI1JhkS2EYvuxkpsYTU+mkkYhM7iMUQzmCw72fV4NNEsq8JMczgg8imgfAq3NrCXwDXAJcGkS4hCJ6dRTYcGCP7fbsJShjKArMymgDpO4gtHcRC4tfJ9Ti+1IKvFzh3Sume2O95d87Yj9gd0E55z7w8z6Ay/iTWWd7JzTPZ+SdKEQ9OgBhYXe9t9YRBbZXMBT/Mwu3M113MUNfEczX+fTILOkKj/rOdyOV0tpJX/2+wd+E5xzbh4wL8jPEPEjMxPGj99+3wm8yTCGcwYv8iONuI1/M46BbCTd1znVSpBU56db6WK8st35ZR4pUo2UbCWA4x+8SBbZnMB/+Z49GMxIxtOPLexa6rm0hoJUNX6SwxKgEfBDwLGIpIzIMQWjiPN5mmEMpwOLWM2+XMO9TOKKMuseqYUgVZWf5DACbzrrEiLWcXDOnRtYVCIJFnnDWrE0/qArMxnKCNqwjC/Zn95MYgaXxa17pBaCVBd+ksMjwCjgM6Ao2HBEEq/kzKOd+J2eTGUwo9iPr/mUv3IJM3mczhSRFvMctWvD1KlKClJ9+EkOG5xz9wYeiUiClRxTqM8v9CGHQdzJ3qzjfY7kWu5hLufgSrlfVDOOpDrykxwWmdkI4Fm271bSeg5S5WRmemMALuJ++93YRH/u51ruoQl5vEpHujONVzmZ0uoeaTxBqjM/yeGw8OPREfu0noNUKdEzj7y6R9dxN1fzALuxmbmcTTZZvMcxcc+jVoLUFH5ugkvqug4iOyLWPQp7s5ZB3EkfcqjHbzzBRQxnGJ/QPu55NKYgNY2vVUXM7GygLbBtHQfn3G1BBSWyo2IlhVasYDCj6MEjGI4ZXMZIhrCcA+OeRxVRpabyc4f0BKA+0BF4GLgI+CDguETKLdZ0VIC2LGEoI7iExyigDg9xJWO4MW7dI7USRPy1HI51zh1qZp865/5jZncBTwYdmIhfsVoJAB34kGEMpxNPs4UG3MUN3M11pdY90piCiMfPeg6/hh+3mtlfgAKgZXAhifgTCkHduiUTg+NE3uBFTudDjuTvvMGt3EJzchnM6LiJoUEDb/lNJQYRj5+Ww1wzawSMAT7Cm6n0UKBRiZSh5I1r4DiDF8gim+N5m+/Yk5sYxXj68TMNY55DrQSR+MpsOTjnbnfObXLOzQGaAwc5524OPjSRaJmZ3iBxZN2jC3mCRRzOfM4ig9X05z5a8jVjuClmYlArQaRscZODmR1hZntFbHcHZgO3m1njRAQnAtsvw1nchVSbAi5nGktpyxN0pgE/04vJ7M8KHqA/v7HzdueoXdtLCM7Bli0abBYpS2kth4lAPoCZnQiMBKYBPwE5wYcmNVWsNZmLZyDV5TeuYgLLOYBp9CCfnejCYxzMMqbSK6ognpm37GZBgRKCSHmUNuaQ5pzbGH7eBcgJdy3NMbPFwYcmNVG8mUe78DNXMZEbuIu/8C3vcRQDuJe5nEOsEhdah1lkx5SaHMystnPuD+AUoI/P94n4Fu/ehGKN+JH+3M9AxtGEPBZwMpcxg9foiJKCSHBK+5KfCbxhZhvwprO+BWBm++N1LYlUWFlJoSk/bKt7tCtbeI5zyCaL97cr8fUnraMgUrniJgfnXLaZLQCaAS85t62OZS3gmkQEJ9VTvK4jgH1Yw42M4Uoeoi6/M5uLGcFQPqVdzOOVFESCUWr3kHPuvRj7lgcXjlRXsUplR2rFCoYwku5Mw3BM53JGMoQvOSDm8SqXLRIsjR1I4EprKRzCZwxlBF2YRQF1yKEPY7iR1TSPOrZWLa8rSmMKIsFTcpBAhUKxE8MRfEAW2ZzHs2yhAXcyiLu5ju/Za7vj0tLgkUfUQhBJND+1lUR8y8z0/sKPvEfhT46/8zovcRofcBQn8Ba3cCvNyWUIo6ISQ3q6EoNIsqjlIJUiFILevSE/P9arjrOYxzCGcxzv8C17MYgxTOSq7cpbaBqqSOpQcpAKCYVg4EDIy4t/TC0KuYAnGcZwDmMxuWSQyQNMpje//7lulAaXRVKQupWk3DIzve6ieImhNgV05xGW0pbHuZj6bKUnU9ifFYwnc1tiSE/36h2p1pFI6lHLQXzx01Koy2/0ZjI3MZoW5LKYdnRmNk9yAUWkbTtup51g8mQlBJFUpuQgZQqFoFcvr3hdLLvwM32ZwA3cRTO+4x2O4WoeYB5nUbLEhbqQRKqGlOtWMrNbzewbM1sc/jkr2THVRKEQNGny54yjWIlhdzbyb24jl+bcyY0spS0deZXjeJt5nE1kYlAXkkjVkqoth7udc3cmO4iaxk/XEcAefM/1jCWTB2nIzzzDuQxnGB9w1LZj1EIQqdpSNTlIgpV2F3OxfVnNjYzhXzzMTuRvq3v0GYduOyY9HcaNU1IQqepSrlsprL+ZfWpmk81s92QHU50Vdx+Vlhhas5xJ9GYlrejLBB7lUg7iCy5lJuvSD922wppzsGGDEoNIdZCUloOZvQIlbof1ZAHjgdsBF368C+gd4xx9CK8xkZGREVis1VlZxfD+yqcMYzideZx8dmI8/biTQawhg1q1YMY0JQKR6spcvG+GFGBmLYC5zrlDSjuuQ4cObuHChQmJqboIheDyy2MnhiN5nyyyOZfn2ExDHiSTu7mOH9gT0FRUkerCzBY55zrEei3lupXMrFnEZidgSbJiqa5CIejRo2RicHTkVV7mVN7naI7jbf7NbTQnl6GM3JYY0tOVGERqglQckB5tZu3xupVWAVclN5zqJboryXE2z5NFNsfwHt+yFzdwJxO5CmvQQDOORGqolEsOzrnLkx1DdREKQVYW5OZ69ytEthRqUciFzGEYw2nPJ6yiOf14kCn0okF6PSZqxpFIjZZy3UpSfqEQtGjhJYDatb3HWrW8m9dyc71jihNDbQroyRQ+pw2z6UI9fqMHU2nNl0y0fvTuV08zjkQk9VoO4l+sm9YKC73HkgPN9fh1W92j5qzmY9pzEY/zFJ0oIo20NJiutRNEJEzJoYoprasolgZs2Vb3aC++522OpR/jmc+ZFJe3MNOiOiKyPSWHKiJWK6G0xLA7GxnAvQzgXhrzIy9xGl3I4k1OJLLmkRn07avEICLbU3KoAkIh6NMHtm4t+9g9+Y7rGUs/xtOQn3ma8xjOMD7kyKhjVepCROJRckgxxd1Gq1dD48bevrIK4QFkkLut7lEdCphFF0YwlCX8NepYJQURKYuSQwqIN47gJykcwP8YwkguYwYOYxo9GMlgVqXtT2EhNG8O2dlKBCJSPkoOSVayy8hvNZND+WRb3aPfqMeMhpk0umMQ/xqwL/8KLlwRqSGUHJIgsuuoVq0/p5/6cTTvkkU25/A8W6why84ZTNuHr6PXHnsEF7CI1DhKDglWsqXgLzE4TuZVssjmZF5jY610Prngdto91J+2jRoFGa6I1FC6QzrBsrL8zTryOM7hOd7lGBZwKgfbFyy69C4a/7SKdo//HygxiEhAlBwSbPXqso9Jo5AuPMZnae15jnPZgx8Y1ngCr0/6isND13trcIqIBEjJIWDFdY9q1fIei6enlpSWBjuRz6D0yWzc62AeoyuHHFAA06axX8FyhuddRdde9RIZuojUYBpzCFDJ8YXcXKhTx1ssJz//z+Ma7/wrL3SexBGvjYY1a+Cww+D+J6BTJy+riIgkmL55AhRrfKGgABo29O4/2JXNjGg0irV1WnDEtGsgIwPmzYNFi+DCC5UYRCRp1HIIUNzxhbw8Vl19L9x7L2zaBKef7mWSE09MaHwiIvHoT9NKVNb4wl58y2huJNeaw223QceO8MEH8OKLSgwiklLUcqgkpY0vNMtfxU2MpjeTqUMBq4/pSsucodC2bXKDFhGJQ8mhksQaX9iv4AtuqTuSiwjhMJ5o0JNd/jOY865vlZwgRUR8UnKoJJHjC+1YTBbZXMgcfvu9HnUGXg2DBnHpPvskL0ARkXLQmEMFlBxbCIW8iUbH8A5zOZvFHMbpvMQIhnL8Prlwzz2gxCAiVYhaDuUUPbbgmHnFAl5Kz+YAXmcD6WRxBw9wNQX1G5EzMrnxiohUhFoO5VQ8tmAUcS7P8B5HM/f309j1++Us6jaWE/bNZYRl0ah5I3JytI6CiFRNSg4lxOoyirQ2t5BLmMkntOMZzqcp6+nDRFoUfsXhM65j2UGUDRsAAAl0SURBVOpdKCqCVauUGESk6lK3UoRY01H79PGed+ucD9On82XtkbT8YwVLacNlTOcxLqGQ2jRvnry4RUQqm5JDhFjTUd3WrXx5zcMwZAysXcuuLQ/nkm+eZHb+ebhww6t+fW8pThGR6qLadiuV1T0US+R01IZsZjAjWUULbv1xILRsCS+8QPrKD/nn5E5kNK+FmVcjSWMLIlLdmPO7aHEK69Chg1u4cOG27ZLdQ+D9dV/Wl3iLFvBz7gYGcC/XcB+7s4n5nMHkPYfx+HcnBHcBIiJJYGaLnHMdYr1WLVsOsbqHtm719se1bh3z29zAKlpwM7fzKidzOAu5qP58zr9LiUFEapZqmRziVUONuf/rr6FfP2jZkoNfGsf64zpxarOldLY55DU/XF1GIlIjJSU5mFlnM1tqZkVm1qHEa0PNbIWZ/c/M/lGR82dk+Ni/bBn06AGtW8PkydCzJyxfTsv/TueVdW00HVVEarRktRyWABcAb0buNLM2wCVAW+AM4EEzSyvvybOzvTGGSNtmFH38MVx0kVcR9YknYMAA+OormDgR9tuvgpcjIlK9JGUqq3NuGYCZlXzpPOAx59zvwNdmtgI4Eni3POcv/ms/K8vrSsrIgJweb3N6KBvmz4fddoNhw2DgQGjadEcvR0Sk2km1MYe9gTUR22vD+6KYWR8zW2hmC9evXx/1erdusOprR9ELL7Gq+d85/bbj4cMPYfhw7+62O+5QYhARiSOwloOZvQLsFeOlLOfcM/HeFmNfzLm2zrkcIAe8qaxRB+TlwZlneglh7729yqhXXhnd3yQiIlECSw7OuVMr8La1wL4R2/sA6yoUQOPG3o1rV14J3btD3boVOo2ISE2UauUzngUeNbOxwF+A1sAHFTqTGcyaVYmhiYjUHMmaytrJzNYCxwDPm9mLAM65pcBs4HPgBeBq51xhMmIUEanJkjVb6SngqTivZQMqYycikkSpNltJRERSgJKDiIhEUXIQEZEoSg4iIhJFyUFERKIoOYiISJRqsRKcma0HcpMdxw5qAmxIdhAJUpOuFWrW9epaq5bmzrmYReaqRXKoDsxsYbzl+qqbmnStULOuV9dafahbSUREoig5iIhIFCWH1JGT7AASqCZdK9Ss69W1VhMacxARkShqOYiISBQlBxERiaLkkELMbIyZfWFmn5rZU2bWKNkxBcXMOpvZUjMrMrNqOR3QzM4ws/+Z2QozG5LseIJkZpPN7AczW5LsWIJmZvua2Wtmtiz8//DAZMcUBCWH1PIycIhz7lBgOTA0yfEEaQlwAfBmsgMJgpmlAQ8AZwJtgK5m1ia5UQVqKnBGsoNIkD+AG5xzBwNHA1dXx/+2Sg4pxDn3knPuj/Dme3hraFdLzrllzrn/JTuOAB0JrHDOfeWcywceA85LckyBcc69CWxMdhyJ4Jz71jn3Ufj5FmAZsHdyo6p8Sg6pqzcwP9lBSIXtDayJ2F5LNfwCqenMrAVwGPB+ciOpfElZJrQmM7NXgL1ivJTlnHsmfEwWXtM1lMjYKpufa63GLMY+zRuvRsysATAHuNY5tznZ8VQ2JYcEc86dWtrrZtYDOAc4xVXxm1DKutZqbi2wb8T2PsC6JMUilczM6uAlhpBz7slkxxMEdSulEDM7AxgMnOuc25rseGSHfAi0NrOWZrYTcAnwbJJjkkpgZgZMApY558YmO56gKDmklvuBhsDLZrbYzCYkO6CgmFknM1sLHAM8b2YvJjumyhSeWNAfeBFvwHK2c25pcqMKjpnNBN4FDjSztWZ2RbJjCtBxwOXAyeF/p4vN7KxkB1XZVD5DRESiqOUgIiJRlBxERCSKkoOIiERRchARkShKDiIiEkXJQVKGmaVHTA38zsy+CT/fZGafJziW9pHTE83s3IpWVjWzVWbWJMb+3cxsmpmtDP+EzGz3HYk7zufHvRYzu9XMBlX2Z0rVp+QgKcM5l+eca++caw9MAO4OP28PFFX255lZaRUC2gPbvlCdc88650ZWcgiTgK+cc62cc62AFXjVTStbIq5FqhklB6kq0szsoXD9/JfMbGcAM2tlZi+Y2SIze8vMDgrvb25mC8JrYywws4zw/qlmNtbMXgNGmdku4bUIPjSzj83svPAdzbcBXcItly5m1tPM7g+fY8/wehufhH+ODe9/OhzHUjPrU9rFmNn+wOHA7RG7bwPamdmBZnaSmc2NOP5+M+sZfn5zON4lZpYTvmMXM3vdzEaZ2QdmttzMTijrWkrEFO932Tn8WZ+YWbUssS7RlBykqmgNPOCcawtsAi4M788BrnHOHQ4MAh4M778fmBZeGyME3BtxrgOAU51zNwBZwKvOuSOAjsAYoA5wMzAr3JKZVSKWe4E3nHPtgL8BxXc+9w7H0QEYYGbppVxPG2Cxc66weEf4+cfAwWX8Lu53zh3hnDsE2BmvFlex2s65I4FrgVvC5cJLu5ZI8X6XNwP/CF/vuWXEJtWECu9JVfG1c25x+PkioEW4KuaxwOPhP54B6oYfj8FbTAhgOjA64lyPR3wpnw6cG9HvXg/IKCOWk4HusO0L/afw/gFm1in8fF+8hJYX5xxG7Cqtsaq5ltTRzG4C6gON8ZLTc+HXiovALQJa+DiX96Gl/y7fBqaa2eyI80s1p+QgVcXvEc8L8f5irgVsCo9LlCXyi/iXiOcGXFhy4SEzO6o8wZnZScCpwDHOua1m9jpeoolnKXCYmdVyzhWFz1ELOBT4CC9BRbbs64WPqYf3F30H59waM7u1xOcU/54KKd+/77i/S+dc3/Dv42xgsZm1d87FS3pSTahbSaqscA39r82sM3jVMs2sXfjld/AqoQJ0A/4b5zQvAtdE9NsfFt6/Ba8IYiwLgH7h49PMbFdgN+DHcGI4CG/5yNJiX4HXhfR/Ebv/D1jgnFsN5AJtzKyume0GnBI+pjgRbAj/tX9RaZ/j41qK44n7uzSzVs65951zNwMb2L4UuVRTSg5S1XUDrjCzT/D+Gi9einMA0MvMPsWroBlvEfjb8cYYPjWzJfw5QPwa3pfzYjPrUuI9A/G6dj7D675pC7wA1A5/3u14y7yWpTdeWe8VZrYeL6H0BXDOrQFmA5/ijZl8HN6/CXgI+Ax4Gq80eFlKu5ZI8X6XY8zss/Dv503gEx+fKVWcqrKKpAAzOxCYhzcgPC/Z8YgoOYiISBR1K4mISBQlBxERiaLkICIiUZQcREQkipKDiIhEUXIQEZEo/w+hDraJRQbL2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZdrH8e9NQBBBkYDKqgFELOAKrtjLim1tq6IiIkpzRYgIFpSSd9VVQ1UUGxClCIwIihXBhnXtoKggLoISQCwQRFDUxOR5/zgTHDIzyUnImZkkv8915Zo5Z86cuU+UufOUcz/mnENERCRSrWQHICIiqUfJQUREoig5iIhIFCUHERGJouQgIiJRaic7gMrQpEkT16JFi2SHISJSpSxatGiDc65prNeqRXJo0aIFCxcuTHYYIiJVipnlxntN3UoiIhJFyUFERKIoOYiISBQlBxERiaLkICIiUZQcRESqqVAIWrSAWrW8x1DI/3urxVRWERHZXigEffrA1q3edm6utw3QrVvZ71fLQUSkGsrK+jMxFNu61dvvR1KTg5lNNrMfzGxJxL5bzewbM1sc/jkrmTGKiFRFq1eXb39JyW45TAXOiLH/budc+/DPvATHJCJS5WVklG9/SUlNDs65N4GNyYxBRKQ6ys6G+vW331e/vrffj2S3HOLpb2afhruddo91gJn1MbOFZrZw/fr1iY5PRCSldesGOTnQvDmYeY85Of4GowEs2WtIm1kLYK5z7pDw9p7ABsABtwPNnHO9SztHhw4dnArviYiUj5ktcs51iPVayrUcnHPfO+cKnXNFwEPAkcmOSUSkpkm55GBmzSI2OwFL4h0rIiLBSOpNcGY2EzgJaGJma4FbgJPMrD1et9Iq4KqkBSgiUkMlNTk457rG2D0p4YGIiMh2Uq5bSUREkk/JQUREoig5iIhIFCUHERGJouQgIiJRlBxERCSKkoOIiERRchARkShKDiIiKW5H1oKuKK0hLSKSwnZ0LeiKUstBRCSF7eha0BWl5CAiksJ2dC3oilJyEBFJYTu6FnRFKTmIiKSwHV0LuqKUHEREUtiOrgVdUZqtJCKS4rp1Cz4ZlKSWg4iIRFFyEBGRKEoOIiISRclBRESiKDmIiEgUJQcREYmi5CAiIlGUHEREJIqSg4iIRFFyEBFJomQs5OOHymeIiCRJshby8UMtBxGRJEnWQj5+KDmIiAQsXtdRshby8UPdSiIiASqt6ygjw9suKeiFfPxIasvBzCab2Q9mtiRiX2Mze9nMvgw/7p7MGEVEdkRpXUfJWsjHj2R3K00FziixbwiwwDnXGlgQ3hYRqZJK6zpK1kI+fiQ1OTjn3gQ2lth9HvBI+PkjwPkJDUpEpBKVtQZ0t26wahUUFXmPqZAYIPkth1j2dM59CxB+3CPWQWbWx8wWmtnC9evXJzRAEZF4Sg4+n3VWEruOnIP58+G778r91lRMDr4453Kccx2ccx2aNm2a7HBERLYNPufmet/LubnwyCPQo0eCu46KimDOHOjQwctOOTnlPkUqzlb63syaOee+NbNmwA/JDkhExI94g8/z5nldRoErKICZM2HECPjiC2jdGiZPrlAmSsWWw7NAj/DzHsAzSYxFRMS3pN238NtvMGECHHCA10zZaSd47DFYtgx69fK2yynZU1lnAu8CB5rZWjO7AhgJnGZmXwKnhbdFRFJeWYPPle6XX2DsWNhvP+jXD/bcE559FhYvhi5dIC2twqdOareSc65rnJdOSWggIiKVIDt7+xveIKDB502b4L77YNw4yMuDk0+G6dO9R7NK+YhU7FYSEamSAr9v4YcfYOhQryly881wzDHw7ruwYAGcckqlJQZIzQFpEZEqq1u3AGYirVkDY8bAQw/B779D584wbBi0a1fJH/QnJQcRkVS1YgWMHAnTpnlzYy+/HIYM8QaeA6ZuJRGRHRDIYj1LlsCll8KBB8KMGd5AxooV3rTUBCQGUMtBRKRcQiHvfobVq6FxY9iyBfLzvdd2eLGeDz6A4cPhmWegQQMYNAiuuw722qvS4vdLLQcREZ9K3gGdl/dnYihW7sV6nIPXX4fTToOjjoI334Rbb/U+ZNSopCQGUMtBRMS3WHdAx+LrprfiukfZ2fDOO949CqNHQ9++0LDhDse6o5QcRER88nunc6k3vRUWwpNPet1Hixd7B99/P/TuDTvvXClxVgZ1K4mI+OTnTue4N70VFHhV+A45BC6+2GuCTJniDTRffXVKJQbwkRzMrJWZ1Q0/P8nMBphZo+BDExFJLbFWbqtTB9LTS7np7bffYPx4b5ZRz55Qty7MmgWff+5t16mTwCvwz0/LYQ5QaGb7A5OAlsCjgUYlIpJCiqerXn659wd+ZDKYMgU2bIixWM/PP8Odd0LLlpCZ6Q0sP/ccfPyx13LYgbpHieBnzKHIOfeHmXUC7nHO3WdmHwcdmIhIspQ2XTUvz2s9TJ8eZ7rqjz/+Wfdo40avrMWjj8JJJ1VqeYug+Wk5FJhZV7zy2XPD+1KzHSQisoMqPF31+++9u5ebN4dbboHjjoP33oNXXoGOHatUYgB/LYdeQF8g2zn3tZm1BGYEG5aISHKUe7rq6tVe3aOHH/ayyMUXe8XxDj000DiDVmZycM59bmaDgYzw9tdojQURqab8Tlc9sdmXcEW47hFA9+5ey6F16+CCSyA/s5X+CSwGXghvtzezZ4MOTEQkGcqarnoInzE7rSuvfnuQN5bQty+sXAmTJlWbxAD+xhxuBY4ENgE45xbjzVgSEanyIgvnNWnizTwqqU4dOG3X93mGc/mMQzm/zlxq3XSjNz3pvvsCXOotefyMOfzhnPvJth9McQHFIyKSMMWDz8VjDHl5JY9wnLfr64zfN5tmSxd4U5cG/oc611wDu++e6HATyk9yWGJmlwJpZtYaGAC8E2xYIiLBiz/47Dib5xnGcI7d/C7k7eXds3DVVV611BrAT7fSNUBb4HdgJrAZuDbIoEREghYKedNVI9WikM7M5mMOYy7/5C+sI5MH4euv4YYbakxiAH+zlbYCWeEfEZEqr7g7qVhtCuhGiCGM5CD+xxccSA+m8iiXsnfzOlAvebEmS9zkYGbPUcrYgnPu3EAiEhEJUCgEPXp4xVHr8Su9mcxNjKY5q/mY9nRmNk9yAUWkxS+iVwOU1nK4M2FRiIgkQHGLYefCLfRlAjdwF3vxPW9zLP0Yz3zOJD3dcBuheYaXGCq0ols1EDc5OOfeSGQgIiJBGz1kI4O23sdAxtGYH3mZU7mEx3iDvwNG8+be7FQpvVtptnPuYjP7jBjdS865qn1vuIjUGHMe+I71w8by383jacjPPM15DGcYH3LktmNqchdSLKV1Kw0MP56TiEBERCpdbi6vnzOGs5ZMYifymUUXRjCUJfx1u8PS0mKsw1DDxZ3K6pz7Nvw00zmXG/kDZCYmPBGRCli+nJUn9qKgxf4cuySHEN04kP/RjUejEkP9+t4CbUoM2/Nzn8NpMfadWdmBiIhUVCjklb5oZ58wy7pQdOBBNHtrFg+SSStWciUPs5L9o94Xc+U2AUofc+iH10LYz8w+jXipIfB20IGJiJQlFIKBA2H/vPeYQjb/ZC6bacgoBnM317GePeK+V4PPpSttzOFRYD4wAhgSsX+Lc25joFEBZrYK2AIU4tV36hD0Z4pI1ZHZz7F8wqvMYjin8Cp5NObf3Mb99GcTpdc9MtPgc1lKm8r6E/AT0NXM0oA9w8c3MLMGzjmfVc93SEfnXIwaiSJSU4VmOOZlzuWaLdkczfusoxnXcxc59OEXyi5vYeZV2VZXUunKLJ9hZv3xynZ/DxSFdztAU1lFJHEKC5l05hN0eHk43fiUr2lBX8YzlZ787rO+RXq6t7SzEkPZ/FRlvRY40DkXVcw2YA54ycwcMNE5l5PgzxeRVJCfz7v9Z7DHpJFcUfQlyziI7jzCTLryh8/l7JUUys9PcliD172UaMc559aZ2R7Ay2b2hXPuzeIXzawP0AcgoxoutCFS4/36K0yaxMahoznm5zV8xGFcyBM8RSdcGRMta9WCoiJv0Lkml8DYEX6Sw1fA62b2PF7ZbgCcc2MDi8o7/7rw4w9m9hTeanRvRryeA+QAdOjQQYsPiVRxxTOP8vM205cJXM9Y9uJ7Puc4spnIC5wBWNz3F48lPPhg4mKuzvwkh9Xhn53CP4Ezs12AWs65LeHnpwO3JeKzRSSxQiFvDZ26v+QxgHsZwL3sziZe5HQuJou3OLHMc6jbqPL5Wc/hP4kIpIQ9gafCS5PWBh51zr2QhDhEJECZmfDU+G+5hbH0YzwN+IWnOJ/hDGMhR5T5frUWguNntlJT4Ca81eC2TQlwzp0cVFDOua+AdkGdX0SS7/8uz6XtjNGMZRJ1KOAxLmEEQ1nKIb7er9ZCsPx0K4WAWXgF+PoCPYD1QQYlItXXs2P+x+ZhI7jljxAOYyo9Gc1NMctbxNOvn1oLQfOTHNKdc5PMbGB4jYc3zExrPYhI+SxezKKLhnPOyif4jXo8wNXcySC+YZ9ynUaJITH8JIeC8OO3ZnY2sA7K+V9TRGqud97x5pPOm8f+7MpIhnAP15Za9ygWdSMllp/kcIeZ7QbcANwH7ApcF2hUIlKlhWY45mQu4Jot2XTkdTaQzt3cwQNczU80KvW9DRrAhAlKAsnmZ7bS3PDTn4COwYYjIlVaURGvD5pL63uyedJ9wDf8hesYSw592MouZb5dXUapw89spSnEXia0dyARiUjVU1jIfwfMZveJIzip8DO+oiV9mMgj9CCfur5OocSQWvx0K82NeF4P6IQ37iAiNV1+PjP+MZ2jXh/J8azgcw7mMqbzGJdQ6OvrRd1IqcpPt9KcyG0zmwm8ElhEIpLyHpvyK4v7P0zm1jFcxhoW8TcuYA5Pc36ZdY8iqbWQuvyl9u21BlTpTqQm2ryZh//2IP9ceTeX8ANvcTx9yOFF/kFpdY9iUWJIbX7GHLbgjTlY+PE7YHDAcYlICnl8Qh4rB4zjqoL7+BebeIF/kE0W/+WEcp9LU1KrBj/dSg0TEYiIpJZQCG65ch19fx1LXybQgF+YwwWMYCiLKN+qvWolVD2lJgcz2xnoBrQJ71oIPOGcyw86MBFJvOKy2Q3yVjGYUSxlMmkUMpOujGAoy7Z9FfhTq5ZXcVWJoeqJO3JkZn8FlgEnAKuAXOAfwNtm1sjM7khIhCISqMxMr7qpGdx+2RfcldeDFexPbyYzlZ4cwHK6M913YqhXD2bMAOegsFCJoaoqreVwL3Clc+7lyJ1mdiqwBFgaZGAiEpziNRR++cXbbs/HDGM4FzKH36jHfVzDnQxiHXv7PqfGEqqX0pJDs5KJAcA594qZFeDd7yAiVUgoBL17Q364Y/hY3iaLbM5iPj+xK8MZxjgGsoGmvs95yinwiia3VzulTUiuZWZRtzaaWT2gwDm3NbiwRKSyhUJw+eWQn+84lZd5jZN4m+M5gg8ZRjYZrObf3OE7MTRo4HUfKTFUT6Ulh2nAHDNrUbwj/Hw2MD3IoESkcoVC0P2yIs51T/M+R/Eyp7M/K7iWu2nBKkYwjM3sVuZ5ihOCc7Bli7qQqrO4ycE5dwfwAvCmmW0wsw3AG8DLzrnbExWgiFRcKAT1d/qD5y97lMW042k6kU4eV5JDK1YyjmvLLIhn5k1FVUKoWUqdyuqcux+438wahre3JCQqEdkhmZnw8Ph8ujONTxnJ/qxkKW3oxgxm0cVX3SMNMNdsvspnKCmIVB2HH7yV4794iJXcyb6s5UM6cD5P8Sznllr3SAXwJFJFaiuJSIoJheDGPj/RY+uDzOdu9mA9b3AiVzCJlzmN0uoeabaRxKLkIFLFXXDiBg57axyfcx+N+In5nEE2WbzN8WW+V2UtJB4/hffq4y0RmuGcu9LMWgMHRqwQJyJJMLTHOppOu5PpTGRnfuVJLmA4w/iYv/l6vxKDlMZPy2EKsAg4Jry9Fnic7RcBEpEEyMyE+eO/ZjCjuJUppFHIo1zKSIaUq7zFww9rbEFK52dVjlbOudFAAYBz7lfKW7hdRCosFPIGi9vY5xw9vjtf0ppeTGEKvTiA5fRgmq/EkJ7u3aPw669KDFI2Py2H/HB1VgdgZq2A3wONSkS2lbpom/8Rj5BNJ57iV3ZmHAMZy/W+6x61aQNLVQlNyslPcrgF72a4fc0sBBwH9AwyKJGaLhSCCZf9l6fJ5kxeYBO7kU0W4xhIHk18ncMM+vbVuIJUjJ/Ffl42s4+Ao/G6kwY65zYEHplIDZTZz7Fywktkkc1bvMUPNGUow3mQTF/lLUAtBakccZODmZWc8vBt+DHDzDKccx8FF5ZIzREKQd8+RZy69RmGMZwjWMga9mEA43iYf/Er9X2fSzOQpLKU1nK4q5TXHHByJcciUuP07/sHP06cxbuM4BCWsoJW/IuHmEZ3CtjJ1znUfSRBiJscnHMdExlISWZ2BjAOSAMeds6NTGY8IpUlMxMmjf+dHjzCYEbRiq9YQlsuJcRsLvZV9wjUSpBglTmV1czqmdn1Zvakmc0xs2vDazoExszSgAeAM/HWr+5qZuVbvFYkxYRCsGvaL+w0/h5W0oocriKPdM7jaQ7lU2ZyaZmJIXIJTiUGCZKfP1GmAVuA+8LbXfHWc+gcVFDAkcAK59xXAGb2GHAe8HmAnylSqSKX4tyVn7iaB1jJ3TRlA6/zd3oxhVc4Fb+3DakGkiSSn+RwoHOuXcT2a2b2SVABhe0NrInYXgscFXmAmfUB+gBkZGQEHI6If5FLcTZhPUO5h/7cz25s5nnOYjjDeIfjfJ9PYwqSDH7ukP7YzI4u3jCzo4C3gwvJ+5gY+9x2G87lOOc6OOc6NG3qf71bkcpWfAezmfdz2WXQJP8bxnIdq2jBUEbwEqdzGB9xDs/7TgzFC+wUFSkxSOL5aTkcBXQ3s9Xh7QxgmZl9Bjjn3KEBxLUW2Ddiex9gXQCfI1JhkS2EYvuxkpsYTU+mkkYhM7iMUQzmCw72fV4NNEsq8JMczgg8imgfAq3NrCXwDXAJcGkS4hCJ6dRTYcGCP7fbsJShjKArMymgDpO4gtHcRC4tfJ9Ti+1IKvFzh3Sume2O95d87Yj9gd0E55z7w8z6Ay/iTWWd7JzTPZ+SdKEQ9OgBhYXe9t9YRBbZXMBT/Mwu3M113MUNfEczX+fTILOkKj/rOdyOV0tpJX/2+wd+E5xzbh4wL8jPEPEjMxPGj99+3wm8yTCGcwYv8iONuI1/M46BbCTd1znVSpBU56db6WK8st35ZR4pUo2UbCWA4x+8SBbZnMB/+Z49GMxIxtOPLexa6rm0hoJUNX6SwxKgEfBDwLGIpIzIMQWjiPN5mmEMpwOLWM2+XMO9TOKKMuseqYUgVZWf5DACbzrrEiLWcXDOnRtYVCIJFnnDWrE0/qArMxnKCNqwjC/Zn95MYgaXxa17pBaCVBd+ksMjwCjgM6Ao2HBEEq/kzKOd+J2eTGUwo9iPr/mUv3IJM3mczhSRFvMctWvD1KlKClJ9+EkOG5xz9wYeiUiClRxTqM8v9CGHQdzJ3qzjfY7kWu5hLufgSrlfVDOOpDrykxwWmdkI4Fm271bSeg5S5WRmemMALuJ++93YRH/u51ruoQl5vEpHujONVzmZ0uoeaTxBqjM/yeGw8OPREfu0noNUKdEzj7y6R9dxN1fzALuxmbmcTTZZvMcxcc+jVoLUFH5ugkvqug4iOyLWPQp7s5ZB3EkfcqjHbzzBRQxnGJ/QPu55NKYgNY2vVUXM7GygLbBtHQfn3G1BBSWyo2IlhVasYDCj6MEjGI4ZXMZIhrCcA+OeRxVRpabyc4f0BKA+0BF4GLgI+CDguETKLdZ0VIC2LGEoI7iExyigDg9xJWO4MW7dI7USRPy1HI51zh1qZp865/5jZncBTwYdmIhfsVoJAB34kGEMpxNPs4UG3MUN3M11pdY90piCiMfPeg6/hh+3mtlfgAKgZXAhifgTCkHduiUTg+NE3uBFTudDjuTvvMGt3EJzchnM6LiJoUEDb/lNJQYRj5+Ww1wzawSMAT7Cm6n0UKBRiZSh5I1r4DiDF8gim+N5m+/Yk5sYxXj68TMNY55DrQSR+MpsOTjnbnfObXLOzQGaAwc5524OPjSRaJmZ3iBxZN2jC3mCRRzOfM4ig9X05z5a8jVjuClmYlArQaRscZODmR1hZntFbHcHZgO3m1njRAQnAtsvw1nchVSbAi5nGktpyxN0pgE/04vJ7M8KHqA/v7HzdueoXdtLCM7Bli0abBYpS2kth4lAPoCZnQiMBKYBPwE5wYcmNVWsNZmLZyDV5TeuYgLLOYBp9CCfnejCYxzMMqbSK6ognpm37GZBgRKCSHmUNuaQ5pzbGH7eBcgJdy3NMbPFwYcmNVG8mUe78DNXMZEbuIu/8C3vcRQDuJe5nEOsEhdah1lkx5SaHMystnPuD+AUoI/P94n4Fu/ehGKN+JH+3M9AxtGEPBZwMpcxg9foiJKCSHBK+5KfCbxhZhvwprO+BWBm++N1LYlUWFlJoSk/bKt7tCtbeI5zyCaL97cr8fUnraMgUrniJgfnXLaZLQCaAS85t62OZS3gmkQEJ9VTvK4jgH1Yw42M4Uoeoi6/M5uLGcFQPqVdzOOVFESCUWr3kHPuvRj7lgcXjlRXsUplR2rFCoYwku5Mw3BM53JGMoQvOSDm8SqXLRIsjR1I4EprKRzCZwxlBF2YRQF1yKEPY7iR1TSPOrZWLa8rSmMKIsFTcpBAhUKxE8MRfEAW2ZzHs2yhAXcyiLu5ju/Za7vj0tLgkUfUQhBJND+1lUR8y8z0/sKPvEfhT46/8zovcRofcBQn8Ba3cCvNyWUIo6ISQ3q6EoNIsqjlIJUiFILevSE/P9arjrOYxzCGcxzv8C17MYgxTOSq7cpbaBqqSOpQcpAKCYVg4EDIy4t/TC0KuYAnGcZwDmMxuWSQyQNMpje//7lulAaXRVKQupWk3DIzve6ieImhNgV05xGW0pbHuZj6bKUnU9ifFYwnc1tiSE/36h2p1pFI6lHLQXzx01Koy2/0ZjI3MZoW5LKYdnRmNk9yAUWkbTtup51g8mQlBJFUpuQgZQqFoFcvr3hdLLvwM32ZwA3cRTO+4x2O4WoeYB5nUbLEhbqQRKqGlOtWMrNbzewbM1sc/jkr2THVRKEQNGny54yjWIlhdzbyb24jl+bcyY0spS0deZXjeJt5nE1kYlAXkkjVkqoth7udc3cmO4iaxk/XEcAefM/1jCWTB2nIzzzDuQxnGB9w1LZj1EIQqdpSNTlIgpV2F3OxfVnNjYzhXzzMTuRvq3v0GYduOyY9HcaNU1IQqepSrlsprL+ZfWpmk81s92QHU50Vdx+Vlhhas5xJ9GYlrejLBB7lUg7iCy5lJuvSD922wppzsGGDEoNIdZCUloOZvQIlbof1ZAHjgdsBF368C+gd4xx9CK8xkZGREVis1VlZxfD+yqcMYzideZx8dmI8/biTQawhg1q1YMY0JQKR6spcvG+GFGBmLYC5zrlDSjuuQ4cObuHChQmJqboIheDyy2MnhiN5nyyyOZfn2ExDHiSTu7mOH9gT0FRUkerCzBY55zrEei3lupXMrFnEZidgSbJiqa5CIejRo2RicHTkVV7mVN7naI7jbf7NbTQnl6GM3JYY0tOVGERqglQckB5tZu3xupVWAVclN5zqJboryXE2z5NFNsfwHt+yFzdwJxO5CmvQQDOORGqolEsOzrnLkx1DdREKQVYW5OZ69ytEthRqUciFzGEYw2nPJ6yiOf14kCn0okF6PSZqxpFIjZZy3UpSfqEQtGjhJYDatb3HWrW8m9dyc71jihNDbQroyRQ+pw2z6UI9fqMHU2nNl0y0fvTuV08zjkQk9VoO4l+sm9YKC73HkgPN9fh1W92j5qzmY9pzEY/zFJ0oIo20NJiutRNEJEzJoYoprasolgZs2Vb3aC++522OpR/jmc+ZFJe3MNOiOiKyPSWHKiJWK6G0xLA7GxnAvQzgXhrzIy9xGl3I4k1OJLLmkRn07avEICLbU3KoAkIh6NMHtm4t+9g9+Y7rGUs/xtOQn3ma8xjOMD7kyKhjVepCROJRckgxxd1Gq1dD48bevrIK4QFkkLut7lEdCphFF0YwlCX8NepYJQURKYuSQwqIN47gJykcwP8YwkguYwYOYxo9GMlgVqXtT2EhNG8O2dlKBCJSPkoOSVayy8hvNZND+WRb3aPfqMeMhpk0umMQ/xqwL/8KLlwRqSGUHJIgsuuoVq0/p5/6cTTvkkU25/A8W6why84ZTNuHr6PXHnsEF7CI1DhKDglWsqXgLzE4TuZVssjmZF5jY610Prngdto91J+2jRoFGa6I1FC6QzrBsrL8zTryOM7hOd7lGBZwKgfbFyy69C4a/7SKdo//HygxiEhAlBwSbPXqso9Jo5AuPMZnae15jnPZgx8Y1ngCr0/6isND13trcIqIBEjJIWDFdY9q1fIei6enlpSWBjuRz6D0yWzc62AeoyuHHFAA06axX8FyhuddRdde9RIZuojUYBpzCFDJ8YXcXKhTx1ssJz//z+Ma7/wrL3SexBGvjYY1a+Cww+D+J6BTJy+riIgkmL55AhRrfKGgABo29O4/2JXNjGg0irV1WnDEtGsgIwPmzYNFi+DCC5UYRCRp1HIIUNzxhbw8Vl19L9x7L2zaBKef7mWSE09MaHwiIvHoT9NKVNb4wl58y2huJNeaw223QceO8MEH8OKLSgwiklLUcqgkpY0vNMtfxU2MpjeTqUMBq4/pSsucodC2bXKDFhGJQ8mhksQaX9iv4AtuqTuSiwjhMJ5o0JNd/jOY865vlZwgRUR8UnKoJJHjC+1YTBbZXMgcfvu9HnUGXg2DBnHpPvskL0ARkXLQmEMFlBxbCIW8iUbH8A5zOZvFHMbpvMQIhnL8Prlwzz2gxCAiVYhaDuUUPbbgmHnFAl5Kz+YAXmcD6WRxBw9wNQX1G5EzMrnxiohUhFoO5VQ8tmAUcS7P8B5HM/f309j1++Us6jaWE/bNZYRl0ah5I3JytI6CiFRNSg4lxOoyirQ2t5BLmMkntOMZzqcp6+nDRFoUfsXhM65j2UGUDRsAAAl0SURBVOpdKCqCVauUGESk6lK3UoRY01H79PGed+ucD9On82XtkbT8YwVLacNlTOcxLqGQ2jRvnry4RUQqm5JDhFjTUd3WrXx5zcMwZAysXcuuLQ/nkm+eZHb+ebhww6t+fW8pThGR6qLadiuV1T0US+R01IZsZjAjWUULbv1xILRsCS+8QPrKD/nn5E5kNK+FmVcjSWMLIlLdmPO7aHEK69Chg1u4cOG27ZLdQ+D9dV/Wl3iLFvBz7gYGcC/XcB+7s4n5nMHkPYfx+HcnBHcBIiJJYGaLnHMdYr1WLVsOsbqHtm719se1bh3z29zAKlpwM7fzKidzOAu5qP58zr9LiUFEapZqmRziVUONuf/rr6FfP2jZkoNfGsf64zpxarOldLY55DU/XF1GIlIjJSU5mFlnM1tqZkVm1qHEa0PNbIWZ/c/M/lGR82dk+Ni/bBn06AGtW8PkydCzJyxfTsv/TueVdW00HVVEarRktRyWABcAb0buNLM2wCVAW+AM4EEzSyvvybOzvTGGSNtmFH38MVx0kVcR9YknYMAA+OormDgR9tuvgpcjIlK9JGUqq3NuGYCZlXzpPOAx59zvwNdmtgI4Eni3POcv/ms/K8vrSsrIgJweb3N6KBvmz4fddoNhw2DgQGjadEcvR0Sk2km1MYe9gTUR22vD+6KYWR8zW2hmC9evXx/1erdusOprR9ELL7Gq+d85/bbj4cMPYfhw7+62O+5QYhARiSOwloOZvQLsFeOlLOfcM/HeFmNfzLm2zrkcIAe8qaxRB+TlwZlneglh7729yqhXXhnd3yQiIlECSw7OuVMr8La1wL4R2/sA6yoUQOPG3o1rV14J3btD3boVOo2ISE2UauUzngUeNbOxwF+A1sAHFTqTGcyaVYmhiYjUHMmaytrJzNYCxwDPm9mLAM65pcBs4HPgBeBq51xhMmIUEanJkjVb6SngqTivZQMqYycikkSpNltJRERSgJKDiIhEUXIQEZEoSg4iIhJFyUFERKIoOYiISJRqsRKcma0HcpMdxw5qAmxIdhAJUpOuFWrW9epaq5bmzrmYReaqRXKoDsxsYbzl+qqbmnStULOuV9dafahbSUREoig5iIhIFCWH1JGT7AASqCZdK9Ss69W1VhMacxARkShqOYiISBQlBxERiaLkkELMbIyZfWFmn5rZU2bWKNkxBcXMOpvZUjMrMrNqOR3QzM4ws/+Z2QozG5LseIJkZpPN7AczW5LsWIJmZvua2Wtmtiz8//DAZMcUBCWH1PIycIhz7lBgOTA0yfEEaQlwAfBmsgMJgpmlAQ8AZwJtgK5m1ia5UQVqKnBGsoNIkD+AG5xzBwNHA1dXx/+2Sg4pxDn3knPuj/Dme3hraFdLzrllzrn/JTuOAB0JrHDOfeWcywceA85LckyBcc69CWxMdhyJ4Jz71jn3Ufj5FmAZsHdyo6p8Sg6pqzcwP9lBSIXtDayJ2F5LNfwCqenMrAVwGPB+ciOpfElZJrQmM7NXgL1ivJTlnHsmfEwWXtM1lMjYKpufa63GLMY+zRuvRsysATAHuNY5tznZ8VQ2JYcEc86dWtrrZtYDOAc4xVXxm1DKutZqbi2wb8T2PsC6JMUilczM6uAlhpBz7slkxxMEdSulEDM7AxgMnOuc25rseGSHfAi0NrOWZrYTcAnwbJJjkkpgZgZMApY558YmO56gKDmklvuBhsDLZrbYzCYkO6CgmFknM1sLHAM8b2YvJjumyhSeWNAfeBFvwHK2c25pcqMKjpnNBN4FDjSztWZ2RbJjCtBxwOXAyeF/p4vN7KxkB1XZVD5DRESiqOUgIiJRlBxERCSKkoOIiERRchARkShKDiIiEkXJQVKGmaVHTA38zsy+CT/fZGafJziW9pHTE83s3IpWVjWzVWbWJMb+3cxsmpmtDP+EzGz3HYk7zufHvRYzu9XMBlX2Z0rVp+QgKcM5l+eca++caw9MAO4OP28PFFX255lZaRUC2gPbvlCdc88650ZWcgiTgK+cc62cc62AFXjVTStbIq5FqhklB6kq0szsoXD9/JfMbGcAM2tlZi+Y2SIze8vMDgrvb25mC8JrYywws4zw/qlmNtbMXgNGmdku4bUIPjSzj83svPAdzbcBXcItly5m1tPM7g+fY8/wehufhH+ODe9/OhzHUjPrU9rFmNn+wOHA7RG7bwPamdmBZnaSmc2NOP5+M+sZfn5zON4lZpYTvmMXM3vdzEaZ2QdmttzMTijrWkrEFO932Tn8WZ+YWbUssS7RlBykqmgNPOCcawtsAi4M788BrnHOHQ4MAh4M778fmBZeGyME3BtxrgOAU51zNwBZwKvOuSOAjsAYoA5wMzAr3JKZVSKWe4E3nHPtgL8BxXc+9w7H0QEYYGbppVxPG2Cxc66weEf4+cfAwWX8Lu53zh3hnDsE2BmvFlex2s65I4FrgVvC5cJLu5ZI8X6XNwP/CF/vuWXEJtWECu9JVfG1c25x+PkioEW4KuaxwOPhP54B6oYfj8FbTAhgOjA64lyPR3wpnw6cG9HvXg/IKCOWk4HusO0L/afw/gFm1in8fF+8hJYX5xxG7Cqtsaq5ltTRzG4C6gON8ZLTc+HXiovALQJa+DiX96Gl/y7fBqaa2eyI80s1p+QgVcXvEc8L8f5irgVsCo9LlCXyi/iXiOcGXFhy4SEzO6o8wZnZScCpwDHOua1m9jpeoolnKXCYmdVyzhWFz1ELOBT4CC9BRbbs64WPqYf3F30H59waM7u1xOcU/54KKd+/77i/S+dc3/Dv42xgsZm1d87FS3pSTahbSaqscA39r82sM3jVMs2sXfjld/AqoQJ0A/4b5zQvAtdE9NsfFt6/Ba8IYiwLgH7h49PMbFdgN+DHcGI4CG/5yNJiX4HXhfR/Ebv/D1jgnFsN5AJtzKyume0GnBI+pjgRbAj/tX9RaZ/j41qK44n7uzSzVs65951zNwMb2L4UuVRTSg5S1XUDrjCzT/D+Gi9einMA0MvMPsWroBlvEfjb8cYYPjWzJfw5QPwa3pfzYjPrUuI9A/G6dj7D675pC7wA1A5/3u14y7yWpTdeWe8VZrYeL6H0BXDOrQFmA5/ijZl8HN6/CXgI+Ax4Gq80eFlKu5ZI8X6XY8zss/Dv503gEx+fKVWcqrKKpAAzOxCYhzcgPC/Z8YgoOYiISBR1K4mISBQlBxERiaLkICIiUZQcREQkipKDiIhEUXIQEZEo/w+hDraJRQbL2QAAAABJRU5ErkJggg==\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:                Mon, 24 Feb 2020   Deviance:                     0.087389\n",
      "Time:                        22:49:06   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:                Mon, 24 Feb 2020   Deviance:                   1.0215e-05\n",
      "Time:                        22:49:06   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.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
