{
 "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:                Sun, 16 Aug 2020   Deviance:                       4078.8\n",
      "Time:                        18:00:44   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": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABCHklEQVR4nO2de3hU5bXwfyvDQAICQUCUAIIoIFQRRdFSq6IVrVWpWqmX2mNr+/kd7alWqfj1glhbaKmXHq21HrU3bb2XiqhUC9WKpYhctEHCQZFLUAQhXCcwSdb3x8yEyWTvPXsuezKZrN/z5GFm73e/e82E7PW+6yqqimEYhmGkUtbWAhiGYRjFiSkIwzAMwxFTEIZhGIYjpiAMwzAMR0xBGIZhGI6YgjAMwzAcMQVhGICIDBYRFZFOPsb+h4i8nuH8p4pITfYStpjLt6yGkQumIIx2h4h8ICL7RaRPyvHl8Qfn4DYSLfnhvTvpZ4Wq/kNVhyeN+0BEzvKY53QRaYpfv0tEakTk6izkuU1EHs328xgdG1MQRntlLXBZ4o2IHANUtJ04rahU1YPiP6OznGOTqh4E9ABuAf5HREbmT0TD8MYUhNFe+QNwVdL7rwK/Tx4gIj1F5PciskVE1onI90WkLH4uJCI/F5GtIvI+cJ7DtQ+LyIciUisid4hIKFth4zuCjfHXfwAGAXPiO4Tvel2rMWYD24FWCkJE+ovIcyKyTUTWiMg34sfPAf4fMDmxk8lWfqNjYjZMo72yCPiKiBwNrAYmA58B7kgacy/QEzgC6A38FfgQeBj4BvAFYAywB3gmZf7fAZuBI4FuwPPABuDXuQquql8RkVOBa1T1lXTj40rtQqASeMdhyJ+AaqA/MAJ4WUTeV9WXROQnwJGqemWuchsdD9tBGO2ZxC7ic8AqoDZxIr7anwzcqqq7VPUD4E7gK/EhlwL3qOoGVd0GzEi6th9wLnCDqu5R1Y+Bu4EvZyDbVhGpi//cnOXn6y8idcBWYBrwFVVt4egWkYHEFOMtqlqvqsuBhzjwOQ0ja2wHYbRn/gC8BgwhxbwE9AE6A+uSjq0DquKv+xPbESSfS3A4EAY+FJHEsbKU8enoo6oNiTcicnoG1ybYpKoD0ozpD2xT1V1Jx9YBY7O4n2G0wBSE0W5R1XUishb4PPD1lNNbgSixh/3K+LFBHNhlfAgMTBo/KOn1BmAfKQ/5PJOvMsqbgINFpHuSkkj+nFau2cgaMzEZ7Z2vAxNUdU/yQVVtBJ4Efiwi3UXkcOA7QCLk80ngv0RkgIj0AqYmXfshMX/FnSLSQ0TKRGSoiJyWR7k3E/ON5ISqbgDeAGaISLmIHEvsO3ks6T6DE855w8gE+09jtGtU9T1VXeJy+lvEHNDvA68DfwQeiZ/7H2AesAJYCjybcu1VxExUK4lFDz0NHJZH0WcA38/RR5HgMmAwsd3En4Fpqvpy/NxT8X8/EZGlOd7H6GCINQwyDMMwnLAdhGEYhuGIKQjDMAzDEVMQhmEYhiOmIAzDMAxH2l0eRJ8+fXTw4MFtLYZhGEa74q233tqqqn0zuabdKYjBgwezZIlbVKNhGIbhhIisSz+qJWZiMgzDMBwxBWEYhmE4YgrCMAzDcMQUhGEYhuGIKQjDMAzDkXYXxWQYhtERmL2sllnzathUF6F/ZQVTJg5n0piq9BfmEVMQhmEYRcbsZbXc+uw7RKKNANTWRbj12Vi32UIqCTMxGYZhFBmz5tU0K4cEkWgjs+bVuFwRDKYgDMMwioxNdZGMjgeFKQjDMIwio39lRUbHg8IUhGEYRpExZeJwKsKhFscqwiGmTBxeUDnMSW0YhlFkJBzRFsVkGIZhtGLSmKqCK4RUzMRkGIZhOGI7CMMwjCJm9rJaps+pZvveKACVFWFuu2BUQXYXpiAMwzCKlNnLapny9Aqijdp8rC4SZcpTK4Dgk+ZMQRiGUdIEXbLCaX7Ij4N51ryaFsohQbRJmTWvxtecn+zex50vr8743mAKwjCMEibokhVO8095egVo7CGe6z29EuPSJc1FG5v4/T/Xcc8rq9m7v9FzrBvmpDYMo2QJumSF0/zRRm1WDrne0ysxzuvcq6u3cM49r/Gj51dy3MBK5t1wasb3BttBGIZRwgRdsiKTebK555SJw1v5IADCZeKYNLd26x5+PHclr7z7MYN7d+Whq8Zy5tGHICIZ3xtMQRiGUcL0r6yg1uHBnK+SFW7zu43NlIRJKl0U0676KPctWMMjr6+lc6iMqeeO4Orxg+nSKeQ4r19MQRiGUbJMmTi8hY8A8luywmn+cEha+CByvadXwlxTk/L00o387KUatu7ex5dOGMCUc4ZzSPfyrO6ViikIwzBKlqBLVrjNH+Q9Eyxdv53pz1WzYuMOxgyq5OGvjmX0wMq83kNUW4dQFTNjx47VJUuWtLUYhmEYbcJHO+r56Uur+POyWvr16MLUc0dw4egqysq8/Qwi8paqjs3kXraDMAzDaAfURxt5+PW1/HLBGhqalOvOGMp/nn4k3boE9xg3BWEYhlHEqCrzqjfz4xdWsmFbhImj+vG9z49kUO+ugd/bFIRhGEaRUvPRLqbPqeaN9z5hWL+DeOyacYw/sk/B7m8KwjAMo8io27ufu15ezaOL1tG9PMztF47i8pMG0SlU2NxmUxCGYRge5LuWk9d8DY1N/HHxeu56eTU7I1GuPPlwbjxrGL26dc7Xx8kIUxCGYRgu5LuWk9d8h3TvwvQ5K6nZvItTjujNtAtGMuLQHnn6JNlhCsIwDMMFr1pO2VZndZpv6rNvUx9tYkCvCh648ngmjjo06/IY+cQUhGEYhgv5ruXkdl19tImbzx7GNaceQXk4t/IY+cSquRqGYbjgVj8p21pObtcd2qOc6yccVVTKAUxBGIZRIsxeVsv4mfMZMnUu42fOZ/ayWl/nvJgycTgVKQ/tXOoqTT5xIKkJzxXhEFPPHZHVfEFjJibDMNo9Xs5fIGtHc75qOX28q55ZL9Xw1Fsb6V7eiZAIdZEoVQHVacoXgdZiEpFzgF8AIeAhVZ2Zcr4n8CgwiJiy+rmq/sZrTqvFZBhGKuNnzncsu10VN+m4nVs4dUKgcu1vaOI3C9dy7/w17Gto5Gvjh3D9hCPpXh4O9L5OFFUtJhEJAb8EPgdsBN4UkedUdWXSsOuAlap6voj0BWpE5DFV3R+UXIZhlB7ZOJPz1TTICVVl/qqPuWPuu6zduoczRxzC9847miP6HhTYPYMgSBPTScAaVX0fQEQeBy4EkhWEAt0lFs91ELANaAhQJsMwipRcEtLSNQYKsmlQKms+3s3tz6/ktdVbGNq3G7+9+kROH35IIPcKmiAVRBWwIen9RmBcypj7gOeATUB3YLKqNqVOJCLfBL4JMGjQoECENQyj7cg1IS1dY6AgmwYl2BGJ8t9/+19+98YHVHQO8YMvjOSqUw4nXODyGPkkSAXhlOWR6vCYCCwHJgBDgZdF5B+qurPFRaoPAg9CzAeRf1ENw2hLck1I8+NMnjWvhtq6CCGR5rmTr82WxiblySUb+Pm8Grbt3c+XTxzITWcPp89BXXKatxgIUkFsBAYmvR9AbKeQzNXATI15yteIyFpgBLA4QLkMwygyck1IS2eeSrzOZ9kMgMVrt3Hbc9Ws/HAnJw0+mN+dP5JPVfXMaq5iJEgF8SZwlIgMAWqBLwOXp4xZD5wJ/ENE+gHDgfcDlMkwjCIknQ/BC7/mqXyWzaitizDjhXd5/u0P6d+znHsvG8MXjj2sKMpj5JPAFISqNojI9cA8YmGuj6hqtYhcGz//APAj4Lci8g4xk9Qtqro1KJkMw8g/+ah2ms6H4HVfJ8Xi9ODPR9mMyP5Gfv3aezzw6nuowrfPPIprTxtKRefiyoDOF4EmyqnqC8ALKcceSHq9CTg7SBkMwwiOfFU7zTQhLfW+TqQ++HPZpagqc9/5kBkvrKK2LsJ5xx7GreeOYECv4Lu6tSWWSW0YRtZkarZx2m0k5tlUF6FnRZjKrmE21UU8nchO900l9cGfzS4FoHrTDqbPWcnitds4+rAe3HXpaMYd0dvzmlLBFIRhGFmTidnGabcx5ekVoBBtigUn1kWizeO9diPpzEJOD/5Mdymf7N7Hz/+6msffXE+vrp35yRePYfKJAwmlFlMqYUxBGEY7Jt/dzjIlE7ON06o/2ugdte62G3G7L+BZ32jSmKq030+0sYnf/3Md97yymsj+Rq7+9BC+feZR9Oxa+PIYbY0pCMNop+S721k2eJltUpWX2wM9HU67Bbf7zrjomBaf3UuBOp3r1a0zt8+p5r0tezj1qD5MO38kRx7SPSu5SwFTEIbRTsl3t7NscDPbQOucA6F1pqwfnHYjfsxFmVZ4/c6Ty2lSGNy7Kw9dNZYzjz6k5MJWM8UUhGG0U/Ld7SyZTExXTmab8TPnt1JeCq2URKhMaGxyVxteTuR05iIvBZp4nUyTQo/yTsy78bN06VSaYauZ0n6LhBhGByff3c4SJFbetXURlAMrb79NdsBdSSkxH4EAvbqGHR9A3TqHkPi4VJNRJngpUDdz1876BlMOSZiCMIx2Sr67nSVIt/JOx+xltZS5mGYSPRjWzjyPrp07NUcvJVPZtTNrZ57HwqkTcjKVuSlKrxpJoQ5uUkrFFIRhtFMmjalixkXHNK/Ic11xJ8jFdJXYfTQ6NCJLVV5BmsjAWYGGRNiye5/rNU5yd2TMB2EY7Rg/YZuZkkvGsVsCW0iklfLK5T5+SNzrZy+tYtOOegQQgetOH8qzS2v5cEd9q2uq8nTvUsF2EIZhtCAX05Xb6r9JtZUiC8pElkBVKQ+HCIViZqOzR/Vj/k2nM2XiCG45Z0Sg9y4VbAdhGEYLMs04TiaTXUEu90lHzUe7uP35ahau+YRh/Q7isWvGMf7IPgW5dykh2s5sbmPHjtUlS5a0tRiGYTjgVETPKYEtKOr27ueul1fz6KJ1dC8Pc9PZw7j8pEF0CpW1edZ5WyMib6nq2EyusR2EYRh5o61W5g2NTfxx8Xruenk1OyNRrjz5cG48axi9unUGiiPrvD1iCsIwjLyuroNwnHvxxpqtTJ+zkprNuzjliN5Mu2AkIw7t0WJMMWSdt0dMQRhGB6e9rq43bNvLj+e+y0vVHzGgVwUPXHk8E0cd6lgeI+iQ2lLFFIRhdHDa2+p6z74GfvX393jwH+8TEuHms4dxzalHUB52z4AOOqS2VDEFYRgdnGJYXfsxcakqf1m+iRkvvsvmnfuYdFx/pp57NIf2LE87f7bNgjo6piAMo4PT1qtrx0ZCT61g+pxq6vZG6V9ZweQTB/L3mo9Zur6OYwf05P4rjueEww/2fQ8La80OC3M1jHZMPpzLQYWm+pVt/Mz5vnpFdC/vxA++MJJLjh9AWQfq6pYvLMzVMDoQ+XIuB7G6zkQ2v42EDurSiUvHDsxaJiNzTEEYRjslXdXVTB746UJTE7uB2roIIREaVR1beyaPS8XN8Z2YLx0f7qhn9rJaMwsVEFMQhtFOcXMiJ1br+QpbTd0NJB7mqfM6mar8yJxJBdX2EH5bSmRUrE9EykSkR/qRhmFkw+xltYyfOZ8hU+cyfuZ8zyY9Xk7kTPo5fH/2Owy99QUGT53L0Ftf4Puz32lx3q1Ca+q8XuO8ZD6sR/ooJKf7GcGTVkGIyB9FpIeIdANWAjUiMiV40QyjY5FpJzenaqheOK3evz/7HR5dtL55Fd+oyqOL1jPyBy823zdduGvifLpxqWGljU3KnxavZ9e+hlZjwyF3J7QltxUOPzuIkaq6E5gEvAAMAr4SpFCG0RHJtJNbcsMgPzit3v/0rw2OY/dGm5qVU7pw18R5r3GpzYwWr93G+fe+zq3PvsPIw3pw89nDWjQ+mnXJaNfPZclthcOPDyIsImFiCuI+VY2KSPuKjTWMdkA2CWsJ5/KQqXPx+qN0Swrzsv8nlJNTkpnTvG7JaMmKobYuwowX3uX5tz+kf89y7r1sDF849jBEhOsnHNVqfktua1v8KIhfAx8AK4DXRORwYGeQQhlGRySXhDW3awEEuPgE5yildBFEm+oiLcJgvaKYvMJlI/sb+fVr7/HAq++hCt8+8yiuPW0oFZ3dTWSW3Nb2ZJUoJyKdVLW14bAAWKKcETRt1Tcgl4S1dBFEVZUVLJw6odXxhA/CDbfr/KKqzH3nQ2a8sIraugjnHXsYt547ggG9ujp+BlMGwRFIopyI9AN+AvRX1XNFZCRwCvBwdmIaRvHSlpVNc1kxJ8bc8MRyx/NuZqo7Jh0DwGP/Wk/qWjFXc071ph1Mn7OSxWu3cfRhPbjr0tGMO6K349j2WlG21Em7gxCRF4HfAN9T1dEi0glYpqrHFELAVGwHYQSJW9mHXFfShSIX+XNZwSdf269HOUP6dGPR2k/o1bUzN589nMknDiTkUR6jvX/v7YFsdhB+opj6qOqTQBNA3LTkHexsGO2UYqhsmgtOoa9+dwKTxlSxcOoE7p58HAA3PrE8bS4GtA7P/WhnPf98/xM+e1RfFtx0OpePG+SpHKD9f++lih8FsUdEekMsSEJETgZ2BCqVYbQRbg7h9hJamRz6mggZzaToXqa5GOCeILfm49307Br2dd/2/r2XKn6imL4DPAcMFZGFQF/gkkClMow2olj7BmRi/sml5WemzYPWbt3jGj1VWxdh/Mz5vmQu1u+9o5NWQajqUhE5DRhOLGKuRlWjgUtmGG1AMYZWFtKB69fUs6s+yn0L1vDI62sRcMzBEA5Uak0nczF+74a/KKarUg4dLyKo6u8Dkskw2pRcVuBBUMiWoOlyMZqalKeXbuRnL9Wwdfc+LjlhAMdW9WTGi6tayOikNNLJXGzfu+HPxHRi0uty4ExgKWAKwjAKQCEduGeM6OuYF3HGiL4sXb+d6c9Vs2LjDsYMquThr45l9MBKAHpUhFus/t3MTuZ0bl/4MTF9K/m9iPQE/uBnchE5B/gFEAIeUtWZDmNOB+4BwsBWVT3Nz9yG0VHIJMM612SzBau2OB5/5q1aHl20nn49unD35NFcOLqqRVe31NW/W9iql9PZEuVaUgzfR0blvuPsBVoXTUlBRELAL4FzgZHAZfEku+QxlcD9wAWqOgr4UhbyGEZJ4zd0NZsIpFTcVviRaCPXnTGU+TedzhfHpG/5mWm4bT5kLyWK5fvwU+57jog8F/95HqgB/uJj7pOANar6vqruBx4HLkwZcznwrKquB1DVjzMT3zBKH7+hq5lWg3XCbYXfr3sXpkwcQbcu/nqMZRpu6yb7TU+u8NUbo9TIx+8yH/j5bf886XUDsE5VN/q4rgpIriW8ERiXMmYYsWqxfwe6A79wcn6LyDeBbwIMGjTIx60No3jIh6nAjwM3H76KK8cNYtZfa2hK8jBXhEPc+vmjm9/7/TyZOJ3dZHTrXlfqFEvioB8fxKtZzu20B00NbOgEnEDM8V0B/FNEFqnq6hQZHgQehFipjSzlMYyCU8gQ1Vyqwdbt3c9dL6/m0UXrKA+HCIfK2BGJtqrYOntZLVOeWkG06cCDe8pTK4DcPo+XYztBUJFbxUguv8t84qogRGQX7uHNqqrpWo9uBAYmvR8AbHIYs1VV9xDL2H4NGA2sxjACIgjnX2LO1HLYe/Y1OJoKbnhiOTc8sZzKijC3XTAqLw+9bJLNGhqb+OPi9dz18mp2RqJcefLh3HjWMHp16+w4/rbnqpuVQ4Jok3Lbc9XNfamz+W69ek4k01GioIolcdBVQahq9xznfhM4SkSGALXAl4n5HJL5C3BfvABgZ2ImqLtzvK9huBLEij51zmSzSDrqItG8rMCTr/f7gH5jzVamz1lJzeZdnHJEb6ZdMJIRh3qv++oizjmydZFoTt9tquxlLn0qOkrpjWJJHPTncQJE5BBieRAAJBzLbqhqg4hcD8wjFub6iKpWi8i18fMPqOq7IvIS8DaxYoAPqeq/s/gchuGLbJLO0q2K3WoR+SXapEyfU+14/0xW5OnGzl5Wy/Q51Wzfe+Ahf3C3zjxw5fFMHHUoIt6RSenINaEv2Wfh1hujI5XeKIbEQT+Z1BcAdwL9gY+Bw4F3gVHprlXVF4j1sU4+9kDK+1nALP8iG0b2ZOr887MqzofZY/ve2Ao89YE+5ekVRBuT7P1PO+820sk5e1ktNz+1nIamlvfduXc/9dEm38qhV9dwCwWTfDyfjtViWUF3dPzkQfwIOBlYrapDiDmUFwYqlWEERKZVQ/2EG6Yze/TqGqbKh2kkNYRx+pzqZuWQINoY221kIqdq7JpU5QDQoK3v68W080cRDrVUJuGQMO38UXmvyJooP7525nksnDrBlEMb4EdBRFX1E6BMRMpUdQFwXLBiGUYwZJrA5WdV7DRn8tzTzh/FwqkTuGfyca0erl73clqpJ46n5ga4yVlbF+HiX73hOldijBOzl9Uyfub8FveaNKaKWZeMbpHfMOuS0UwaU5VTLwqjOPHjg6gTkYOA14DHRORjYvkQhpFXClFaIFPThZ9ww+Q5U6OYkudO/Hvjk8tbtfdMnTMdydm1XnIC/O/Huz3nEnA0b3mZrKwia8fAT8vRbkCE2G7jCqAn8Fh8V1FwrOVoaeLmlMyk2U0xy/X92e/wp39toFEVkdhDOTlaVIBPDz2YDz6JND9ct+3ZRyTqYBdKoaqygsG9K1j43rZW50Ye1p1VH+4i3SyprT3dailVVoRZPu3stDIZxUdQLUe/CfRX1QZV/Z2q/ndbKQejdCmW0gKp5NqhDWLK4dFF65vDNlVbKgeI7QgWvretRe2dhib19QdaWxdxVA4Q2z2kVzGtTVRuJqtEOKvRMfBjYuoBzBORbcTqKT2tqpuDFcvoaBRLaQEncg03/NO/NqQf5EC0UenVNUzXzp08cwPSzeGHVPOWl8kq22zmYqhOamRG2gWKqk6PV1q9jlio66si8krgkhkdilLuSZzpQz2Zur3R5kieOy8d7enkzhYnR7KXYzkbpV0s1UmNzMik3PfHwEfAJ8AhwYhjdFRKOQImlEMCWisFmedKZG4ms0ljqujVNexPJh8UqwnR8MZPotz/BSYDfYGngW+o6sqgBTM6FoWMgCm0qePkI3o5+gjKwNM/kFCQyXWe8olAC8d0KtPOH5W3bOZiNiEa7vjxQRwO3KCqywOWxejgFKK0QCGrqybut/iD7Y7nLj95UHNkUyohEWZcdAyAryJ22ZBuJ5BPpV0s1UmNzPBT7ntqIQQxjEKQa72gZNwquA7uXcEb721Law1asGoLd1462jOM9pQZf8tZOYRDAkqLKqxCrM90OvKltIulOqmRGb6L9RlGKZAvU4dXBVe/pqBNdRHHJLtItJGfvbSKZeu38+GO+ozkSkagedW/ZN02Hlu0vllpKbE+02MPP7ggkUSWRNc+MQVhdCjyZerItYJr8j0TD8lkhbNpRz2/++c6yqR1zkSCkEfYa2ri26x5Na12NIVuwFMM1UmNzMgkiskw2j1uZpVte/ZlFHKZq3M1XCYtzCtuCqdHeSfCZa2joMIh4c5LR3PP5ON8RX+Zk9jIBlcFISK7RGSn208hhTSMfLFg1RbH45FoU0Zx+bk6V2d9aXTzajra2ORqltoRaeCg8tYb/WijNq/+/WR6l3KeiREcaTvKicjtxPIf/kDMrHkFkGu3OcNoE7xWzMlx+els5X5bZLqRuE+vbp253aF8d4L+lRVpV/+pjXZmzavhxieWt5DdnMRGNvgxMU1U1ftVdZeq7lTVXwEXBy2YYQRBuhVzIuw1XcZv8sodDiTDVVVWMH7owaRLjauti/CdJ5fz1UcW09ikXPOZIZR3avnnmHiA+139e2Ur56OmlNHx8OOkbhSRK4jVYVLgMiD/QdmGkSN+EuDSrfwTUUTJRKKN3PDEcm54YjkAnUNCty6dqNsbpX9lBfdMPq5VqexEVdbKrmH2RRvZ61CVtUljPoZ5N36WLp1CfKqqp6v8TjLv3d/Q/PCfvayWm55c0cppneyINiexkSl+FMTlwC/iP0qsm9zlQQplGJniNwEu8Tq1NzPEVux+TEb7G5X98WudWnsmy+HVqAdgV30DXTqFmq/36rNw23PV1EUOzLd9b5Rbn32HJeu28cxbta4RTZk4oq2gnpGMn2J9H6jqharaR1X7quokVf2gALIZhm8yqfUzaUwVy354NvdMPq6VycVPa9BUUn0Xmfgl/DqJJ42poluX1uu5SLSRP/1rg+c9y0R8Od/zWVDPqRud0f7wU4tpGPAroJ+qfkpEjgUuUNU7ApfOMHySTRin24o9G+dz4j6ZrNYzdRK7zZ2uWmyjqq9yIvnKMi90ORMjOPw4qf8HuBWIAqjq28CXgxTK6Nhks/rMVxhnqjPXL4n7+L1fotZSJg9Mt7n9VIv1Uzk1X7kSVrm1dPCjILqq6uKUY9aT2giEbM0cmZYL91JCk8ZUNfdguPLkQb7kPmNEX1SVoYd0Szu2IhzizktHZ7yadvuMl40b2Oq4E14P+tnLailzUTSZKllLyisd/Dipt4rIUOKV6EXkEuDDQKUySg6/zs9MzBypc158QhULVm1Je490JpDUeccPPTht8b2/Vm9m8dptrN68u8Vxp17T2Tp+veoZjT384Objbp3n3B70ie/D6ZpsciWscmvp4EdBXAc8CIwQkVpgLbFkOcPwRSY2ab+rT6c5H120nq7hMnpWhNlUF2k2aaSGoHqFg6YWtauti7Btz37ujoeyDpk611FRfLxrH1t272t1XIFF72+nSTUvUUFOfpNUhXbGiL4881at76Q4N8d6NmYwsMqtpYQfE9M6VT2LWMOgEar6GVVdF7BcRgmRiU3ary/B7aG2N9pEXSTqaJ7yWikTH5+sHJxkdZOvW+cQbr7iRtXA2mw6meSeeauWi0+o8p0U56aUm1Sz3ulYUl5p4GcHsVZEXgKeAOYHLI9RgmRik/a7+vRrz042T02fU+0ZneRVHTVxvykTh/Pdp99mf+OBxLdOZRAOleEnf9QrKiibHASnzxSJNrJg1ZYW1VwTPhenuYMwCVlSXmngZwcxHHiFmKlprYjcJyKfCVYso5TIJMIo1+JzTmyqizB7WW3apDWvcNH+lRVs2LaXl/79EfsbmwjFK6xWVnRCkBYJbH7kSSUb57zXZ0q+R7q5S7kfuJEbfjrKRYAngSdFpBexjOpXgfRhE0aHwWv1m6lN2s/qM5Nief0rK3IOsRzZvwdn3vUqIRFuPnsY15x6BOXhEONnzqcu4hzUJ+Dor3BSbtnkIHh9puR7pDPxJc4nd8WzDGoDfPaDEJHTROR+YClQDlwaqFRGuyLdCjUIm3Rizl5dw57jEooo1xDLl1duBtVY5vLiDdwxdyXjZ8737B6n0LzTSCDEvp/U0NpsQkO9ziUrX7dxyYUJIbaDSnxfphwM8JdJvRZYTmwXMUVV9wQtlNG+8LP6DcImnZgzefeSWhyvS7xCqpud3Qk3X8T+xgNtRR9dtN7XXE1NsRV5bV2kxY4iNZIrGz+A2zWVFeEW37XbOLfChIXsMmcUN547CBEJAb9R1S+q6p9MORhOtHViVHJi27TzR6FJOdB1kVhBuzNG9PWVTAbpS1dkQmKmyoqwa3TU7GW17N3f2kyVzg/g5ju47YJRvsblo7ifUdp4KghVbQTOKJAsRjulmLqVue1mFqza0sLMVVkRxqGTZyDU1kVcndiJnUSqs7myIpzWDOfXdOc2zq0woSW0GQn8hLm+ISL3EQtzbd5BqOrSwKQy2hX5SIzKtsx06nVuZqRNdZFWmciVXcPsikRpyN+GIWOczDwA3bp08vX5/Zru/BYmtOglIxk/CuLT8X9vTzqmwASHsUYHxKsEhB+yrf7pdJ1X5JBTr4bkTUTvbp35ZM9+XzLnA6/+E4Uw8+T6ezNKH9E82lsLwdixY3XJkiVtLYaRR8bc/lfHeP6qyooWyV6puEURpSqJinCIGRcdw6x5NY7je5aHWfKDswiHyjznFIl1gcuWXl3DdO3cqcXD2E0mp7H24DZyQUTeUtWxmVyTNsxVRPqJyMMi8mL8/UgR+Xq2QhpGMn6TvZxwMycpONrl3cbvrI/GM6HdHbp3Tz7OtZRGKuEyIRxq6eCoCIeYdv6o5h7TiVpRZ4zoSzjFGVImsLu+IS+NewwjF/yYmH4L/Ab4Xvz9amL+iIcDkskoEgrRftJvspcTbuGoIZFWO4/Fa7cRDgnRRu8qp15ml0xW+05zAK1MYk8s3kBqt+omjdVBSsbCT422wI+C6KOqT4rIrQCq2iAivtpticg5xDKvQ8BDqjrTZdyJwCJgsqo+7U90I0gy9Qtkq0y8chO8nKWzl7n3YG5Uba47dEiPLvTvWcGyDXVUdg2zZ19DCyWRcMr6kf+MEX1bFfNL7Ay8ekknGD9zfiufQzQDm5WFnxqFxk8m9R4R6c2BfhAnAzvSXRTPofglcC4wErhMREa6jPspMC8DuY2AyaQCa7ZNfmYvq3Xt2paa7OV0PzcS2coKbN65j2Ub6pg4qh//nHomsy4Z3cr8BKSVf/ayWp55q7aFchDg4hP8JwDm+oC38FOj0PjZQXwHeA4YKiILiZX9vsTHdScBa1T1fQAReRy4EFiZMu5bwDPAiX6FNoLHqzzD+JnzW6yws+1lPGtejWPEkUCrZK/U67xqMDnN+e/anVR0DjmGezqt7FPld7qnAgtWbXGUIbEjqa2LeFaJ9YuFnxptgZ9ifUtF5DRiVV0FqFFVP6Urq4ANSe83AuOSB4hIFfBFYiGzrgpCRL4JfBNg0CB/LSCN3PDKKUg1N2WbSe12XvEOb81mJe51jdfnTHe9V2XWhEJJpxxCZUKjh6kp28Y9hpErfqKYvgRUqGo1MAl4QkSO9zG3k/Ug9a/gHuCWeMa2K6r6oKqOVdWxffv29XFrI1econmS8dNEJ51JxO28W4av33kzvSbk0os5+XgmnzHdDieV7l06eX7mbBv3GEau+PFB/EBVd8V7QEwEfgf8ysd1G4GBSe8HAJtSxowFHheRD4iZre4XkUk+5jYCZtKYKi4+ocrVRwAtm+hk008gl+syId2cXs7u5Hv6lTXTHc6OSJSFUydY6Quj6PCjIBJLofOAX6nqX4DOPq57EzhKRIaISGfgy8R8Gc2o6hBVHayqg4Gngf9U1dl+hTeCZcGqLY72/ASJB1e25bxzua5bZ3+F9/yYZ9wezCGRrEqWZ/pAr4yXLLfGPUax4cdJXSsivwbOAn4qIl3woVji4bDXE4tOCgGPqGq1iFwbP/9ADnJ3GPJVoyibHAavlXDqgyuTct65yLZ26x6u/+NS9uz3Z8LxY55xaz7UqNrC1+L3M2bSzAhoTsCz0hdGseFHQVwKnAP8XFXrROQwYIqfyVX1BeCFlGOOikFV/8PPnB2JfNYo8nNdKl59BLJ1mmYr2676KPctWMND/1jr6NAtcymD4Wc1n7jvTU+uaGVuyiZBLflB7yeKaUdSpVfr5WwUE36imPbGfQTnxhPfFqrqXwOXzMgpfDST69xW9G5VWr2UQ7rdQaayNTUpzyzdyM/m1bBl1z46hwSndXmP8jD7Gpqyrkw6aUwVNz6x3PFcNlFTbuG0mTYFMoy2xE8U0w+JOaZ7A32A34jI94MWzMi+EU82IZlOSWKZ+gj8JMxlItvS9dv54v0LmfL02wzoVcF3zhrW3NUtlR2RaM5tTYPua2E+BqO94cfEdBkwRlXrAURkJrHe1HcEKVgxUojaRMlk04Yy0+vSregzMXm4zXXDE8uZNa+muVBdOtk+2lHPT19axZ+X1dKvRxeuHDeI+as+5q5XVrveW+P3z+V3kq6vRa6/f/MxGO0NPwriA6AcqI+/7wK8F5RAxUq+7PqZkG0jnkyuy2e70HTJaLc++w4Xn1DFM2/VOspWH23k4dfX8ssFa2hoVK47YyiDenXltjkrfTl8c/2deD3A8/X7Nx+D0Z5wVRAici+xhdk+oFpEXo6//xzwemHEKx6y9QdkQ/JKtWdFmPJwGXV7o75XnJmsVLPdpTjhlX0NLVt/Jst289nDKA+H+Nzdr7JhW4SJo/rxvc+PZFDvro5lMLzI9Xfi9gAv5O/fMIoFrx1EoivPW8Cfk47/PTBpiph8rrS9SF2p1kWizf0IglipZrNLycSpnUqi9WdCtpqPdnH789UsXPMJw/odxGPXjGP8kX1ajM8ULyWVLYX6/RtGMeGqIFT1dwAiUg4cSWz38F7CF9HRyOdK24tCr1QztYv7MbW49U2AA99X3d793PXyah5dtI7u5WFuv3AUl580iE6hslbjM33gu5XOyIVC/f4No5hwbTkqIp2AnwBfA9YRi3gaQLx5kM+CfXmnrVqOpj4YIX3IZzYMmTrXtcLp2pnn5e0+2eIWqllZEaZbl07NZrFoY1OrZLaKcIg7Jn2KvfsbuPPl1eyMRPn00D6s+Xg3m3fWOyqn2ctqufGJ5Z4Z3U5cefIgFqza4qr0MnU45+v3X+hAB8NIkE3LUS8T0yygOzBEVXfFb9AD+Hn859vZCtoeKVQESrGvVN1MKnWRKHXxhK+6SOu1Q2VFmMvHDeLB196nZvMuTjmiN6ce1Yd756/x3I1MGlPFDS75CV48umh98+vUebNxOOfj998WgQ6GkQteO4j/BYZpyoB4g59VqnpUAeRrRVvtIApFPlaqQZbncNtBpKM8XEZ9tIkBvSr4/nlHM3HUoXzmpwsc56qqrGjRMtTtnr26hqmPNvl2YifmdZsv9b75pq3uaxiQ3Q7CK1FOU5VD/GAjzj1ZjDyQbQG7BLl0d/NzXboy4G7UR5voUd6Jjdsj/Oj5d/nL8k2+Hb9uCWbTzh/V3BHOD4l528rhbI5uo73hpSBWishVqQdF5EpgVXAiGZPGVLFw6gTWzjyPhVMnZGR+yKRVaDbXOSmwXvFqpOnYWd8AHFA+FWHn/36p5jQvpTlpTFXa/hGp87qZ65TYKj+dMs2WoDO1DSPfePkgrgOeFZGvEQt1VWJd3yqIdYEzihC/q9RUc5Kb2chpvtQQWiezWDq8xm7bs6+51IfbPZPxE16bHLrrNT5Iv0C2iY+G0VZ4hbnWAuNEZAIwilggzYuq+rdCCWdkjh8nt5OzVHC2G6Ze5+SjmDSmip2RKDNfWsVen2W4vYhEm5jy1ArA30N60pgqlqzbxmOL1rf4DInPVJXiT0kXjhtUWLGV2jDaG36quc4H5hdAFsMnXs7kKROHM+WpFUSTal+Hy6TFKtXJnKTQSkmk1iFyisBpaGzikz37uXf+GqKNTRzUpRO79zXk/BmjTZrRQ9qpuVFCOTg5gBOKzS2suLYuwuxltUyfU832vbGorMqKMLddMCqnB7qV2jDaE346yhlFhC9ncmqeWMp7NzNU4oHq5Bx381Hc8uw7zHhxFeOGHMy8Gz7LHp/KwU8qWybO22wdwF72/5ueWtGsHCAWvjvlqRWB+SgMo9jwU6zPCIBsQ1HTZVrPmldDNKUkdrRRm53Ns+bVuIageYVbuvkoGpuU3159InV7o3zl4cW+wtsqXRLpUsnEeZtt/siUicNdE/GcGhNlurMxjPaM7SDagGxDUSH9StntfOIebg96L2fp92e/4ypP/57l1O2Nes6dyr6G9Moh1SyWjmx7LUwaU5VxzLaFpRodBdtBBIzTTiGXekvpVspebULdonxSnbip8idnJScjwHfPGcFtz1X7jmDykiNBNrZ+pzafyWG6XnNVZVjvycJSjY6CKYgAcXPsuj0g/axM04VKup13u6eAZxbvHXNXup5LrLydSmu43cutN3M+ak0llECm5SycvrNwmdBEazNTpjsbw2jPmIkpQNx2Cm7VRv2sTNNlWrudd0smc7tnbV2E6/+4lK2797vKUlVZkTYBL4FbGG06OTIlm0RBp+9s1pdGc+eXRrdIAqysCDPrS6PN/2B0GGwHESBuO4JGVcIhaeFMziVhasm6ba3MWE67Aj/tNFNNLeEyaREym0A44OD1g5dyyGeyWLbRTG7hp0ErA6vuahQztoMIEM9VscaKzWVab8nJwf3oovVpHd5eO4/Zy2qZ+szbjnb4aJMSKmu54xHgipMHMWlMVc4r/0xrTaWjPZWzyCVYwTAKge0gAsSrpEO0SenauRPLfnh2RnM6mVBSiUQbmT6nunl8YnV6xoi+juN/PPdd6huaXOdralKqKitcE/Oy6dcAMYd1vlfM7amchbUxNYod20EESGLV7kY24ZJ+r9m+N8qUp1d47jSmPvM2lz24iC2793nOpcQevP3jSmLWvJrmVW42YaIJGlXzvmLOtRpuIbHqrkaxYzuIgEkkr+WrCVAmLThTE+ZSqW9o4p/vf0K3LiH27HPflQjekUFuYaIhEdeopQRBrJjbSzmLYm8OZRi2gygA2SZxzV5Wy/iZ8xkydW5zGeogTCU/nnRM2h4PXpFBbp/vsnEDCZelL6rRUVfM2f6/MIxCYQoiB5we4E5kY/Zwc2ACvvsv+KGqsqKFfE647QESD3a3z3fHpGM4qDz9JrVMJO13WIq0J3OY0TFxbTlarBRLy9F8NbF3w6s9pZ/+B35IyAstndl79jX4Sn7r1TWc1snuVi01nUz5bK9qoaSGkf+Wo4YH2XZu84ubn6G2LtJq5VlZEW4VippK184hLj6+qtVqFWi1U/GbGb27viHtit/Lnu6UMJjJd+gnTNRCSQ0je8xJnSVBR6C4OXgTD9VkR+z4mfNdH+oicOOZw/ivs45yPD9+5vysdyJ+Kpu6hZ3OuOgY1yQ7v9+hnzBRCyU1jOwxBZElQUeguEX/NKr6bhcKcOeXRnPR8QNcz+eq0PxkKINzF7Vco7v8KGkLJTWM7DEFkSVBJ2S5hY726hpuFXLqRkiEi44f4GmDd1MwvbqG6dq5U/M1e/c3tGiek8Bv/Sin1Xqu36EfJW2hpIaRPeaDyJKgI1DcQiBVW4ecupHYbXjZ4N3uM+38USycOoG1M89j4dQJTDt/VN5DMnP9Dv2EiVooqWFkj0UxFTFOK/9MylokwlbdoqESBf38RvnkIxoo3xFFFsVkGP7IJorJTExthJ+HlpNp5qcvruLDnfW+7uFVbTXZBu838zjXDGW3/hiJuYOivWRWFxpTnEY6zMTUBmQTetnQ2MQf/vkBO+r9haBWVoQ9q622hQ0+36HBFsKaPfbdGX4IVEGIyDkiUiMia0RkqsP5K0Tk7fjPGyIyOkh5ssFvtnQmZPqgfGPNVs7779f5wV+qGT2gku9OHN5st+/VNdyqnEVFOMRtF4wCissGn++IoqBzUUoZ++4MPwRmYhKREPBL4HPARuBNEXlOVZN7WK4FTlPV7SJyLvAgMC4omTIlKJOI3wflhm17+fHcd3mp+iMG9KrggSuPZ+KoQxER/vOMI1vI6WYq8AozLTT5jiiyENbsse/O8EOQPoiTgDWq+j6AiDwOXAg0KwhVfSNp/CLAPWC/DQgqycrtQanEEtf+a8KRbNge4cF/vE9IhJvPHsY1px5BuUtBvXQ29mKxwec7NNhCWLPHvjvDD0GamKqADUnvN8aPufF14EWnEyLyTRFZIiJLtmzZkkcRvQlqleVk9klQWxfhlmff4b4Fa/j8pw5lwc2nc/2EoygPhwIxdxWSfIcGF9J81t6/+1SKyfRoFC9B7iCcigM5RmiKyBnEFMRnnM6r6oPEzE+MHTu2YHG5Qa2yks0+boluZcBflm/izQ+2N//Rppq7bnxiOUvWbeOOSe5NiYqNfO5mCmU+a6voqyApJtOjUbwElgchIqcAt6nqxPj7WwFUdUbKuGOBPwPnqurqdPMWMg8i24qtyT6Byq5hVGFHJOr4Rzh46ty0clSEQ5SHyxwzmRP9oRes2uLrD93JXwH2oPDCq7JuIpfEMIqdYsuDeBM4SkSGALXAl4HLkweIyCDgWeArfpRDoclmlZWqVJIf6skrz88fcxi/WbgWwb3fQoJItNE1e1qBxxatb57Da3XrtBKe8vQK0FjhvXTXd1TMoWt0VAJTEKraICLXA/OAEPCIqlaLyLXx8w8APwR6A/dLrEppQ6Yazg+5JARlahJxcmwnE4k28r0/v8NNT62gsUkJlwmqSkMOG7nUS90c6U6yObUltWqnLTGHrtFRCTSTWlVfAF5IOfZA0utrgGuClKHQ9mM/q8o9+w88pKNNSjgkVHbuxI5IlDIffZyzlSOTFa+tjg8QdGFGwyhWSj6TutAJQdmsKqONSl3cR3HZuIEZXevWJshJjkxkcxpbapE8frHWoEZHpeRrMRXafpxLO9DaugjPvFVLt86hFrsMN6oqKzhjRF+eeavW1+rWSbZwSFr4INyuL8VInkwollwSwygkJa8gCm0/TjxE7pi7kq279wMQKhMam5Qqj74KCSLRRiorwlSEvct6J0fQjD38YF8+Fjenu9MxP/4L81UYRmlT8gqi0Pbj2roIr7y7ma2799O/Zzm3fv5ovnDsYcSd8I6hs6nsiES5e/JxzXkSqZFOqfJnsrp1G5vueovkMYyOR8kriEIlBEX2N/Lr197jgVffQxW+feZRXHvaUCo6t8xW9ZMk17+yosWDvBjKMlskj2F0PKxhUI6oKnPf+ZAZL6yiti7CcQMr+XBHhI937vOVtJZNIl5b0J5kNQyjNcWWKFfyVG/awfQ5K1m8dhtHH9aDi8ZU8dDra307cttTuYP2JKthGPnBdhBZ8Mnufdz58moeX7yenhVhbp44nC+fOIjP/myBlWQwDKMosR1EwEQbm/j9P9dxzyuriexv5D8+PYRvn3kUPbuGAXPkGoZRWpiC8Mmrq7dw+5xq3tuyh1OP6sO080dy5CHdW4zp6I7cYnCmG4aRP0o+kzpX1m7dwzW/e5OvPrKYhibloavG8vuvndRKOYBzjX0h5oso5szjfGRIW49jwyg9bAfhwq76KPctWMMjr6+lc6iMqeeO4Orxg+nSybnRD7QOYU3OXyjWzON8ZUhbIp1hlB6mIFJoalKeWbqRn82rYcuufVxywgC+e85wDule7nqNk2nFKc8hEm3kpidXcOMTy4vGBJOvB7v5Xwyj9DAFkcTS9duZ/lw1KzbuYMygSh66aiyjB1Z6XuO2AnfLlE5Uai2WHUW+Huwd3f9iGKWI+SCAzTvrufGJ5Vx0/xt8tLOeuyeP5plrP51WOQDc9ly14wo8JG51VluOC6qqrF/cHuCZPtitx7FhlB4degdRH23k4dfX8ssFa2hoVK47Yyj/efqRdOvi72uZvayWuohz4b1GVSrCobRVXdvaBJOvWlWWSGcYpUeHVBCqyrzqzfz4hZVs2BZh4qh+nDj4YH6z8APuX/Ce74eb1+o/JMKMi45pfmC6NQJqaxNMPh/sVhLbMEqLDqcgaj7axe3PV7NwzScM63cQj10zji279mUcyTN7Wa1rsT2glTLoUdGJ3fUNafsutAX2YDcMw4kOoyDq9u7nrpdX8+iidXQvDzP9glFcMW4QnUJljJ85P6NInoRj2ovKinALpbN9bzTWWrQizI549zgzwRiGUcyUvIJoaGziT4vXc+fLq9kZiXLlyYdz41nD6NWtc/OYTCN5nEJDkwmXCSKtG/5EG5VuXTqxfNrZWXyS/GJZz4ZhpKOkFcQba7Yyfc5Kajbv4pQjejPtgpGMOLRHq3GZhmh6OZYrK8LcdsEobnxiueP52roIQ6bObdOHckdvH2oYhj9KMsx1w7a9XPuHt7j8oX+xZ38DD1x5PH/8xjhH5QCZh2i6KY6qygqWTzubSWOqPJ3PbV2Kwis5zjAMI0FJKYg9+xr4+bwazrzrVV5dvYWbzx7GK985jXM+daDlpxOTxlQx46JjqKqsQIg96L0a4fhRKE5jUmmrh7JlPRuG4YeSMDGpKn9ZvokZL77L5p37mHRcf6aeezSH9nQvj5FKpn2dwTs0NHWMW9eNtngoW9azYRh+aPcK4u2Nddz2XDVL19dx7ICe3H/F8Zxw+MGB39ePQkkeM37m/KJ5KOcrOc4wjNKm3SqIj3fVM+ulGp5eupHe3brws0uO5ZLjB1BWlr7ERVtQTA9ly3o2DMMP7U5BqMKvX32Pe+evYV9DI9889Qiun3Ak3cvDbS2aJ8X2ULbkOMMw0tHuelJ3Hzhce19xF2eOOITvnXc0R/Q9qK1FMgzDKHo6RE9qAX579YmcPvyQthbFMAyjpGl3Ya5HHdLdlINhGEYBaHcKwkebBcMwDCMPtDsFYRiGYRQGUxCGYRiGI+0uiklEtgDrcpymD7A1D+IEQTHLBsUtXzHLBsUtXzHLBiZfLiRkO1xV+2ZyYbtTEPlARJZkGu5VKIpZNihu+YpZNihu+YpZNjD5ciEX2czEZBiGYThiCsIwDMNwpKMqiAfbWgAPilk2KG75ilk2KG75ilk2MPlyIWvZOqQPwjAMw0hPR91BGIZhGGkwBWEYhmE4UrIKQkTOEZEaEVkjIlMdzo8QkX+KyD4RubkI5btCRN6O/7whIqOLSLYL43ItF5ElIvKZQsnmR76kcSeKSKOIXFIssonI6SKyI/7dLReRHxZKNj/yJcm4XESqReTVYpJPRKYkfXf/jv9+g+8Q5k+2niIyR0RWxL+7qwshVwby9RKRP8f/dheLyKfSTqqqJfcDhID3gCOAzsAKYGTKmEOAE4EfAzcXoXyfBnrFX58L/KuIZDuIA/6rY4FVxfTdJY2bD7wAXFIssgGnA88X8v9bhvJVAiuBQfH3hxSTfCnjzwfmF4tswP8Dfhp/3RfYBnQuIvlmAdPir0cAf0s3b6nuIE4C1qjq+6q6H3gcuDB5gKp+rKpvAtEile8NVd0ef7sIGFBEsu3W+P8yoBu4ttxuE/nifAt4Bvi4CGVrK/zIdznwrKquh9jfSZHJl8xlwJ8KIpk/2RToLiJCbBG1DWgoIvlGAn8DUNVVwGAR6ec1aakqiCpgQ9L7jfFjxUKm8n0deDFQiQ7gSzYR+aKIrALmAl8rkGzgQz4RqQK+CDxQQLnA/+/1lLgZ4kURGVUY0QB/8g0DeonI30XkLRG5qmDSZfB3ISJdgXOILQIKgR/Z7gOOBjYB7wDfVtWmwojnS74VwEUAInIScDhpFp6lqiCcioIXUzyvb/lE5AxiCuKWQCVKuqXDsVayqeqfVXUEMAn4UdBCJeFHvnuAW1S10WFskPiRbSmxmjijgXuB2UELlYQf+ToBJwDnAROBH4jIsKAFi5PJ3+35wEJV3RagPMn4kW0isBzoDxwH3CciPYIVqxk/8s0kpvyXE9thLyPNDqfddZTzyUZgYNL7AcS0erHgSz4RORZ4CDhXVT8pJtkSqOprIjJURPqoaiGKlfmRbyzweGynTx/g8yLSoKqz21o2Vd2Z9PoFEbm/yL67jcBWVd0D7BGR14DRwOoikS/BlymceQn8yXY1MDNufl0jImuJ2foXF4N88f97VwPEzWBr4z/uFMKBUugfYorvfWAIBxw2o1zG3kbhndRp5QMGAWuATxehbEdywEl9PFCbeF8M8qWM/y2Fc1L7+e4OTfruTgLWF9N3R8xE8rf42K7Av4FPFYt88XE9idn3uxVCrgy+u18Bt8Vf94v/XfQpIvkqiTvNgW8Av083b0nuIFS1QUSuB+YR8+4/oqrVInJt/PwDInIosAToATSJyA3EvP473eYtpHzAD4HewP3xlXCDFqBapE/ZLgauEpEoEAEma/x/XZHI1yb4lO0S4P+KSAOx7+7LxfTdqeq7IvIS8DbQBDykqv8uFvniQ78I/FVju5yC4FO2HwG/FZF3iJl8btHC7Az9ync08HsRaSQWqfb1dPNaqQ3DMAzDkVJ1UhuGYRg5YgrCMAzDcMQUhGEYhuGIKQjDMAzDEVMQhmEYhiOmIIySI17hc3nSz2AReSN+brCIXJ409jgR+XwW9/i7iOQcdpyveQwjCExBGKVIRFWPS/r5QFU/HT83mFhBugTHARkrCMPoCJiCMDoEIrI7/nImcGp8Z3ELcDswOf5+soh0E5FHRORNEVkmIhfGr68QkcfjtfSfACoc7nGuiDyZ9P50EZkTf/0rifXOqBaR6WlkREQuEZHfxl/3FZFn4jK9KSLj48dPS9olLROR7vn4rgwjQUlmUhsdnop4QTKAtar6xaRzU4mVVvkCgIhsBsaq6vXx9z8h1mPgayJSCSwWkVeA/wPsVdVj4zWyljrc92Xg1yLSLZ7lOxl4In7ue6q6TURCwN9E5FhVfdvn5/kFcLeqvi4ig4hlyx4N3Axcp6oLReQgoN7nfIbhC1MQRikSUdXjsrz2bOACOdBlsJxYXazPAv8NoKpvi0irh3u83MFLwPki8jSxiqjfjZ++VES+Sexv7jBitfn9KoizgJHxkisAPeK7hYXAXSLyGLEeDhsz+6iG4Y0pCMNoiQAXq2pNi4Oxh7OfujRPANcRKyb3pqruEpEhxFb7J6rq9rjpqNzh2uT5k8+XAaeoaiRl/EwRmUvMh7JIRM7SWCMYw8gL5oMwOhq7gO4e7+cB34qXQ0ZExsSPvwZcET/2KWKtVp34O7EKt9/ggHmpB7AH2CGxDl7nuly7WUSOFpEyYgXpEvwVuD7xRkSOi/87VFXfUdWfEis8OcJlXsPIClMQRkfjbaBBYh3dbgQWEDPfLBeRycQqcoaBt0Xk3xxohvQr4KC4aem7uNT411iToueJKYHn48dWEGvOUg08Qsw05MTU+DXzgQ+Tjv8XMDbuIF8JXBs/foOI/FtEVhCrDFuoroNGB8GquRqGYRiO2A7CMAzDcMQUhGEYhuGIKQjDMAzDEVMQhmEYhiOmIAzDMAxHTEEYhmEYjpiCMAzDMBz5/4d1PVYu0Ow9AAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7O0lEQVR4nO2df5wdZXXwv2c3N8kmCBtKtLAlJlINikgiq6CxFqySggUjoFSxWmuLvv56oZoarG8JiiVKKX2rb1WsVlspBoWuYKxBTdSKRU2aDSGa1B/8XKgkkkVIlmR/nPePmdnM3p0fz8ydO3fu3fP9fPLJ3rnz48xzZ855nvOccx5RVQzDMAwjC12tFsAwDMNoP8x4GIZhGJkx42EYhmFkxoyHYRiGkRkzHoZhGEZmzHgYhmEYmTHjYbQUEblYRG5P+P7bIvKnBVznDBF5sNHzVAERWSsiX2i1HC50UrsbUzHjYTgjIveKyIiIPCEi/yMinxORIxo5p6reoKpnFSVjXkRERWS/f2+/EpFvichFrZarHRCRPxaRcb/tfi0igyLyBznO8zkRuaoZMhrFY8bDyMq5qnoEsAxYDlzeWnEK5RT/3pYCnwM+LiJXtFaktuE//bbrBT4D3CQiR7dWJKOZmPEwcqGq/wNsxDMiAIjI6SLyfREZFpHtInJG6Ls/FpFfiMjjInKPiFwc2v690H6vEJFdIvKYiHwckNB3U9w1IrLYHzHM8j+/WUR+4l/jFyLy1pz3tldV/wX4X8DlIvIb/vmPEpHPiMjDIjIkIleJSHfoPu4QkY/5su8Skd8LyZp27PdE5G9EZJ/fPmeHjl0iIt/x7+sbwDFheVPa/dsi8iFftsdF5HYROSb0/UtCxz4gIn/sb5/jy3O/iPxSRD4pIj0ObTcBfBboAZ5R/72IPNuXaVhEdorIef72S4CLgb/wRzC3pV3LaC1mPIxciMhvAWcDP/M/9wEbgKuAo4H3AjeLyEIRmQ/8PXC2qj4FeDEwGHHOY4CbgQ/gKcifAysyiPUI8AfAkcCbgetE5Pl57s/nK8As4IX+588DY8Bv4426zgLC8zGnAb/wZb8CuCXU+3Y5drd/7EeBz4hIYDj/Fdjqf/ch4E3BQUntHjr36/Ha46nAbH8fRGQR8O/Ax4CFeB2BQf+YjwDP8rf9NtAH/FVia3nnnOXf1xPAT+u+qwG3Abf7srwLuEFElqrq9cANwEdV9QhVPTftWkZrMeNhZGVARB4HHsBT1oFb5w3A11T1a6o6oarfALYA5/jfTwDPFZEeVX1YVXdGnPsc4Meq+mVVHQX+DvgfV8FUdYOq/lw9voOnpH4nz0365xsF9gJHi8jT8Izlpaq6X1UfAa4D/jB0yCPA36nqqKquxzMGr3Q89j5V/bSqjuMZmmOBp/kK/gXA/1HVg6r6XTwFHJDW7gD/pKr/raojwE0cHi1eDHxTVW/0Zf6Vqg76RuvPgMtU9VFVfRz46zp56zldRIbxfq/XAa9W1cfq9wGOANap6iFV3QR81d/faDNmtVoAo+1YparfFJHfxesRHwMMA08HXiMi4R5jDdisqvv9yef34vWo7wDeo6q76s59HJ5RAkBVVUQewBHf1XMFXo+5C5gH7Mh6g6Hz1fB65I/i3V8NePjwgICusLzAkE6tNHof3j25HDtpJFX1gL/fEXjtu09V99ed93j/79h2jzo3cMA/L/45fh5x6wvx2m5rSF4BuiP2DbhTVV+S8D34v6/v2gq4D29UY7QZZjyMXKjqd0Tkc8DfAKvwFOG/qOqfxey/Edjo+82vAj7N9FHBwxxWivg94OND3+/HU2oBvxnadw6ey+uNwFdUdVREBgjNmeTgVXiuph/iuXsOAseo6ljM/n0iIiEDsgi4Fa9t0o6N42FggYjMDxmQRUBwjcR2T+EBDrvkwuwFRoCTVHUox3njeAg4XkS6QgZkEfDf/t9W4ruNMLeV0Qh/B7xCRJYBXwDOFZGVItItInPFi/H/LRF5moic5899HMTzh49HnG8DcJKInO/7zt9NyEDg+eNfKiKLROQopkZ6zQbmAHuAMX8UkisEWESOFm9C//8BH/HdOQ/jucGuFZEjRaRLRE7wR2ABTwXeLSI1EXkN8Gw8l5LLsZGo6n14bqgrRWS2iLwECI8yYtvd4VZvAF4uIq8VkVki8hsissxX7J/GmzN6qt8mfSKy0uGcSfwArwPwF34bneHfyxf9739JxCS7UU3MeBi5UdU9wD/j+eMfwOupvx9PgT8ArMZ7xrqA9+D1PB8Ffhd4e8T59gKvAdYBvwKeCdwR+v4bwHrgLrwJ5K+Gvnscz9jcBOzDmyS+NeMtbReRJ/CCAP4Uz+cfniR+I56R+rF/jS/jzU0E/MCXeS/wYeBCVf2V47FJvB5vQv1RPLfcPwdfpLR7Iqp6P97cyHv8cw8Cp/hfvw+vHe4UkV8D38QLYc6Nqh4CzsOb/9kL/APwxpD78jPAc/xIrIFGrmU0H7HFoAyjcfwQ1z918PsbRkdgIw/DMAwjM2Y8DMMwjMyY28owDMPIjI08DMMwjMx0RJ7HMccco4sXL261GIZhGG3F1q1b96rqwvQ9p9MRxmPx4sVs2bKl1WIYhmG0FSJyX95jzW1lGIZhZMaMh2EYhpEZMx6GYRhGZsx4GIZhGJkx42EYhmFkpiOirYz2YGDbENds3M1DwyMc19vD6pVLWbXclnIwjHbEjIdRCgPbhrj8lh2MjHqV2IeGR7j8Fm+dJjMghtF+mNvKKIVrNu6eNBwBI6PjXLNxd4skMgyjEcx4GKXw0PBIpu2GYVQbMx5GKRzX25Npu2EY1caMh1EKq1cupafWPWVbT62b1SsbWpzOMIwW0VLjISKfFZFHROTu0La1IjIkIoP+v3NaKaNRDKuW93H1+SfT19uDAH29PVx9/sk2WW4YbUqro60+B3yc0JrMPtep6t+UL47RTFYt7zNjYRgdQktHHqr6XeDRVspgGIZhZKeqcx7vFJG7fLfWgqgdROQSEdkiIlv27NlTtnyGYRgzmioaj08AJwDLgIeBa6N2UtXrVbVfVfsXLsy1lolhGIaRk8oZD1X9paqOq+oE8Gngha2WyTAMw5hK5YyHiBwb+vhq4O64fQ3DMIzW0NJoKxG5ETgDOEZEHgSuAM4QkWWAAvcCb22VfIZhGEY0LTUeqvq6iM2fKV0QwzAMIxOVc1sZhmEY1ceMh2EYhpEZMx6GYRhGZsx4GIZhGJkx42EYhmFkxoyHYRiGkRkzHoZhGEZmzHgYhmEYmTHjYRiGYWTGjIdhGIaRGTMehmEYRmbMeBiGYRiZMeNhGIZhZMaMh2EYhpEZMx6GYRhGZsx4GIZhGJlp6WJQhmEYSQxsG+Kajbt5aHiE43p7WL1yKauW97VaLAMzHoZhVJSBbUNcfssORkbHARgaHuHyW3YAmAGpAOa2Mgyjklyzcfek4QgYGR3nmo27WySREcaMh2EYleSh4ZFM241yMbdVG2P+YKOTOa63h6EIQ3Fcb08LpDHqsZFHmxL4g4eGR1AO+4MHtg21WjTDKITVK5fSU+uesq2n1s3qlUtbJJERxoxHm2L+YKPTWbW8j6vPP5m+3h4E6Ovt4erzT7bRdUUwt1WbYv5gYyawanmfGYuKYiOPNiXO72v+YMMwysCMR5ti/mDDMFqJua3alGAob9FWhmG0AjMebYz5gw3DaBXmtjIMwzAyYyMPw4jAEjANIxkzHoZRhxXkM4x0Wuq2EpHPisgjInJ3aNvRIvINEfmp//+CVspozDwsAdMw0mn1nMfngN+v27YG+JaqPhP4lv/ZMErDEjANI52Wuq1U9bsisrhu86uAM/y/Pw98G3hfeVKVh/nVq4kV5DOMdFo98ojiaar6MID//1OjdhKRS0Rki4hs2bNnT6kCFoEVNqwuloBpGOlU0Xg4oarXq2q/qvYvXLiw1eJkxvzq1cUK8hlGOlWMtvqliByrqg+LyLHAI60WqBmYX73aWAJmY5hLtvOp4sjjVuBN/t9vAr7SQlmahhU2NDoVc8nODFodqnsj8J/AUhF5UETeAqwDXiEiPwVe4X/uOMyvbnQq5pKdGbQ62up1MV/9XqmCtAArbGh0KuaSnRlUcc5jxmB+daMTsVDnmUEV5zwMw2hjzCU7M7CRh2EYhWIu2ZmBGQ+jbbDwz/bBXLKdT6rbSkT+t4gcKR6fEZH/EpGzyhDOMAIs/NMwqoXLnMefqOqvgbOAhcCb6dDwWaO6WPinYVQLF+Mh/v/nAP+kqttD2wyjFCz80zCqhYvx2Coit+MZj40i8hRgorliGcZULCPfMKqFi/F4C96aGi9Q1QPAbDzXlWGUhoV/Gka1iI22EpHn1216hoh5q4zWYOGfhlEtkkJ1r034ToGXFSyLYSRi4Z+GUR1ijYeqnlmmIIZhGEb74JQkKCLPBZ4DzA22qeo/N0sowzDag7jETUvo7HxSjYeIXIG3pvhzgK8BZwPfA8x4GMYMJkjcDPJvgsTNLfc9ys1bh6ZtB8yAdBAu0VYX4pVI/x9VfTNwCjCnqVIZhlF54hI3b/zBA5bQOQNwcVuNqOqEiIyJyJF4y8I+o8lyGR2IuTI6i7gEzXHVTPsb7YmL8dgiIr3Ap4GtwBPAD5splNF5xLk4wFwZ7Urcuh3dIpEGxBI6O4tUt5Wqvl1Vh1X1k3jLwr7Jd18ZhjNWm6rziEvcfN1px1tC5wzAZcL8pVHbVPW7zRHJ6ESsNlXnkZS42f/0o81F2eG4uK1Wh/6eC7wQz31lSYI5mKl+f1uatL1wfU7jEjctobPzSTUeqnpu+LOIHA98tGkSdTAz2e+/euXSKfcOnefK6JSOwUx+Tg138qxh/iDw3KIFmQnMZL//quV9XH3+yfT19iBAX28PV59/cscoo05arGomP6eGOy5zHh/Dq2UFnrFZBmxvokwdy0z3+3eyKyNJ4bbbPc/059RwwylUN/T3GHCjqt7RJHk6GvP7dy6dpHDtOTVccAnV/Xzo3w1mOPJja1J0Lp20WJU9p4YLSet57OCwu2oaqvq8pkjUwVR5TYpOmextFVUOCMj621b5OTWqg2hMKQERebr/5zv8///F//9i4ICqfrDJsjnT39+vW7ZsSd/RiKQ+ugY8xddJE9plUEUDbL+tkYSIbFXV/lzHxhmP0MnvUNUVadtaiRmPxlixblOkj7uvt4c71nRuOk/Zyr4VxmWm/raGG40YD5cJ8/ki8hJV/Z5/sRcD8/NczKgmnTTZ60rZuQytyp2Yib+tUQ4ueR5vAf6fiNwrIvcC/wD8SVOlMkqlWZO9A9uGWLFuE0vWbGDFuk2VynkoO5ehVbkTnTSRb1QLlwzzrcApfjl2UdXHmi8W+IbqcWAcGMs7tOoUmunyaMZkb9WzlMvukRdxvTzPwJknLuQLd94fud0wGiEp2uoNqvoFEfnzuu0AqOrfNlk2gDNVdW8J16k0zVbEzYiuqXrSXNm5DI1eL+8zsHnXnkzbZwJVDGxoR5JGHsG8xlPKEMSIJ04RX3nbzsIe+qKzv6vuay87tLbR6115285cxrjqv0PZVH1E3E7EGg9V/ZT//5XliTNVBOB2EVHgU6p6ffhLEbkEuARg0aJFLRCvPOJe9H0HRhnYNpT40Leql1WFLOWkew/+X3vrToZHRgGYW8tT6s2NRkZ3A9uG2HdgNPK7NCNQhd+hSlR9RNxOpL4tIvJRETlSRGoi8i0R2SsibyhBthWq+nzgbOAd9euKqOr1qtqvqv0LF3a2/zbpRU+acG1lsb5WZym73vvBsYnJv/cdGG1q+6xa3scda17GPeteyR1rXuasrJJ+4zQj0OrfoWrYSKw4XLpaZ6nqr4E/wKuo+yymrvHRFFT1If//R4B/w1tHZEaS9KInPfStrI7a6iq6LvfeLtVjk37jNCPg+jtUOTKuSCz6rDhc8jxq/v/n4BVFfDSYNG8WIjIf6FLVx/2/zwIqk9FeNquW901xr4RJeuhb3ctqZRVdl3tvdvuE3WZH9dQQgeEDo5ndh3Gup96eWiHtW/V5gCJdr1UuI9NuuIw8bhORXUA/8C0RWQg82VyxeBrwPRHZDvwQ2KCqX2/yNSvN2vNOyux+mMm9LJd7b2b71LvNhkdG2XdgNJf7MM71tPa8kzLLEXXtKo/Aina9tnpE3Em4VNVdA7wI6FfVUeAA8KpmCqWqv1DVU/x/J6nqh5t5vXYgz0M/k/3dLvfezPaJUshhsijnRhSei2HIMgIr273VDMOWd+7JmIrLYlDz8IojLsKLbjoOWAp8tbmiGfVkdQO1e3XUwF0xNDxCtwjjqvQ53oPLvTezfVxcX1ncY3ldgC6GIS0iK/w7CIdLbZfh3mq169WIx2XO45+ArcCL/c8PAl/CjEelqfcTX3fRsrYxGjDdDz/uF/DMorBcFG6z5mXiFHL9Ps3GJVQ3aR6g/neoL6Pa7DBXCzWuLi5zHieo6keBUQBVHQGaO2NuNEQnrKed5Papij8e4t04US6xMGW5D11cc0lusTT3GzR3FDCTXa9Vx2XkcUhEevA7HSJyAnCwqVIZDdHsRKgyEg/TFFKz3RYu9+gSpVREtFUjuLrm4kZgLu3czFFAu7teOxkX43EF8HXgeBG5AVgB/HEzhTIao5l+4rLCOtPcPs1UWK73mGakmx2q7GrE6+UIRksuyjjtdyhjFNDKkG8jHpdoq28A5+MZjBvxQnbvaa5YRiM0MwS1rLDOJLdPsxVW3D1eun5wimuqWUbaJaIpr2sy63FRv0Pgs57pYa4zJbEyjsSRh4i8COgDvquqG0TkecDfA78DHF+CfEYOmpkI1UyFWd+Lvvr8k3NHWzVC0r0MDY+w+kvbgeZM5hY16okj63Fluo2KdIc227Va9cTKMkgqyX4NXkmSQeB9IvJV4O3AX2OLQRVOkQ97M174QL64RYvrwzqzXDfuRbz6/JNbslRqmqtmdEJZe+tO1p53UuFG2lW55zXieY4rw21UpDIuQ7FbgcXkkccrgeWq+qSILAAeAp6nqj8tR7T2oAil34yHPeqFzytrvXz1xIV1ut5H1V7EqJFbPcMjo4UY6frfJM5o1Sv3vKOeqoa+FvkMZD1XnvfC8k+SjceIqj4JoKr7RGS3GY6pFKX0i1aeUS8DkFvWpHDNsBtpxbpNHbHmRNgopOVquPbKXX+TcBJemHrlntc1WdXaTkU+A1kz5vO8F1U1wmWSZDxOEJFbQ58Xhz+r6nnNE6s9KErpF/niRL0Ml60fZG6ti5HRiSn7joyO8+c3DXLp+kEAempdzK11TwsnjZNDYIpbKe99xBb+m1dzjgrKQ9p6H6uW97H8g7dHrqWxYF5t2rak60QpqDmzuqY9PwrTDEiUcs876qlq6GuRyjjLufK+w1U1wmWSZDzq61dd20xB2pGilH6RL07Uy6AwzXAETIS01MjoxOR+4R6Yq3x57yPqRax1C088OTapuIv2W7v2OK849yRWf3k7o+OHG6rWLVxxbnpRwoA4BRU3mlO8EV3WEFxXwscFBvSy9YMtNSRFKuMs58r7DlfVCJdJ0kqC3ylTkHakKKVf5ItTpKsn6IG5ypf3PqJexP0Hx6aVoB8ZHec9N22fckxe1t7qtqxrEUoi62/S19uTGihQ1bm2vHIVqYyznKuRd3im55+4JAkaMRSl9It8cVxqKmXhoeERJ/kCpTEyOp4rrLb+RVyyZkPkfuOqhSi4qLVRgMi2S1MSaQoz7jdZMK/Gk6MTmZ+fPEo/SsZmzLU1YoyKVMau5zL3U37MeDRA0b2lrNFRUd+vXrmUy9YPRk66zvPnPeLCbaMIemBJL2NUEcPgBcyrDJKMYKORWEkJjQKp68KHSVOYA9uGOHBobNpxPbXuSddX1ucnTzRRlIxxbrO8o9dGjFEZJW+iMPdTfsx4NEizhq4uSikuN+Li0xdxw533T5t0nVPr4kDM3EcUtS5x6oE1I9Q2LVy2Efdc0rEKmeROy7iPuofenhprzztpyuR8gEvpENdw3jQZgxFiPXkjhvLOH7Q64W6mu5/yklqeRESeJSKfFpHbRWRT8K8M4WYyaUop6furVp3MdRctm1YldTgiaiiJI+bOcnqpmhFqu2q5V+m1O2bJ40ZCItOOzSJ30r3HhTjPnxPdri6lQwa2DcWWtI67rzgZgxFimEZcNnnL4lR5JUMjHpeRx5eATwKfBpJrMxuFkaaQ076P6k255C2EcTU2WSYds7gngu1F+6TTRjVZ5E6696xG1WUEF5flL/59xd1PlIx9obmPVq4PXrU8n6y0yuXWalyMx5iqfqLpkhhTSFPIeaJEzjxxYaQ7a26tKzKXwbV3H3feeqWRxz3RqE866cW+8rad0+47q9xRClP8fbK6hVyUaNw+SnwbJin1oiepIdtvNbBtiK6C3Wdl0mqXWytxMR63icjbgX8jtI6Hqj7aNKmM1F5c1l7ewLYhbt46NEXBC3DBqX30P/3o3L37pPNGjXzyzI3kVXBpL3Ywd7T21p2T0Vdza9M9uUlyByG1Ucu0RinEWnf8PJJLhyApciuOMieFs/xWwe8T1U7tEvFUtdI6ZeJiPN7k/786tE2BZxQvjhGQ9sJnVQhxyYObd+3hqlUnZzqX63nrKcs9EV5zu56oF/vg2OEggn0HRqf1HNMmqAOFuWLdpnS3YEKom0uHYPXKpdOSFgGeeHIsMUqsipPCcXNC3SJtU+q93V1ujZBqPFR1SRmCGNNJe+GzKIQ8cyRJJCnouOvF9ZoVWLFuUyG94bQijvWypfUcgwnqKJ1fL7eLwhid0IbKn69a3jdlpORy3qoS114Tqm1zHzO5xlWq8RCRGvC/gJf6m74NfEpVs4XuGC2lyIfcRUFHnTdporooX3FSEcco2dKMX1IZ+uD4tDIuUcdA/HxM2v0/FpPg2G693U5QvDM5yTA1VBf4BHAq8A/+v1P9bUYbEbUiXN6HPE1Bx503CL/ti1EORYRnpinQsGwuYa8uCjmQ+8wTF8aeL0y3SO6VAMOyuW6vKlmfySqu3LdqeR8XnNo3GVLeLRI539eJuMx5vEBVTwl93iQi25slkNEcipw0TVKoaSVJgp71kjUbInv0jfSekyJ3omRLGlUcODTGkjUbEs8XZmh4hH+tiziLY1y1oYnWsnu7rqGoWUNWszyTVY1qCgJGgmdkXJWbtw7R//SjO96AuBiPcRE5QVV/DiAiz8DyPUqniFjyoiZNk/IGXFf+K9pl8YGBHdPChQN6at2RE7BJhioI4XUxHAGuuftdEu8uGxoeYcmaDanVdKG8pWFdlHZe5Z5W9ia4xygjXoWoJou2SmY1sFlEfoEXhfl04M1NlcqYQit7XXH1sxrt+RbZex7YNhRrOLqE2MidLEUku0WYUOWonhr7D41Ni3bKwkTKoYEba/WXtnPlbTunra8C5UVPuSrHZhdZjDPirZ7nmcnRVqlzHqr6LeCZwLv9f0tVdXOzBTMO06ryDXF+eWBy7iJc/iSLkgjPf4TPAWT2aye5n5IUdZTPPY4JVe5Z90oGrziLay48Jf2AAhidUPYdGM08J1IkrsqxaCXqEvgArZ/n6ZT5pzy4RFu9Bvi6qt4lIh8ArhCRq1T1v5ovngGt692kJcdl7VFGjWLCbq68I6y0dnj/LXdNue6ZJy5k8649PDQ8wlE9NebWuiZ791HriMBUZbBqeZ9zqZcuAdViSuUX4Q7J6v5s9kJgcbg821WIarJoq2T+j6o+LiIvAVYCn8eirUql2b2bcBTLsitvZ/kHb2fJmg2Zq7emXSMtuijvCCutHQ6MTky57hfuvH/y8/DIKE+OTnDdRcu4Y83LWHveSdNGI4JXgiWM66hlwjccq1cujY0yA2ILQNbTaEBB1givuPvcf3BsynFFRvNB/G/aLZJ7tNsM4kbQrZarDFyMR/A2vxL4hKp+BZjdPJGMeop+McPUK5ThkdFJV0kcRS2PW28Y8o6wsrifohgZHWftrTtZsW4Tl60fpD4lUIGbtw5NUZZRSuMNpy+aNBBhUxAo6TNPXEite7qRqHUJrzvteKd7OK63J3fIah7jHNxnffmT4ZHRKYanPgy7W2Ty3HlcbXHP/LWvPYV71r0y18i3Waxa3scda15WObmajcuE+ZCIfAp4OfAREZmDm9FpCBH5feD/At3AP6rqumZfs6oUGV1T77bYt/8gceubR+G60l29rHEGYGh4ZDJLO6/rI6nQoSvDI6OT7qqo9qgvh5+0QNd7btoeGRm0edcerrnwlClyhtf26H/60ZPnjpqY76l1c+aJC3MHTzSyXvc1G3dPa9t6N1pUFeS8wR1lRpQZ+RBNCUUUkXnA7wM7VPWnInIscLKq3t40oUS6gf8GXgE8CPwIeJ2q/jhq//7+ft2yZUuzxOkYXDLD4xBweoGjQmaTKveG97ng1D5u3jo0zX+cxQ2w/IO35zYgLvTUumPlS2tfAe5Z90rna31gYAc3/uABxlXpFm90snnXntxh0nG1t1yOjcvLqb+nRq5hlI+IbFXV/lzHJhkPEekC7lLV5+YVLg8i8iJgraqu9D9fDqCqV0ftf/TTn62veP9nS5SwPdl2/zCHxt1HGWFmd3dx/NE9HHPEnNh99j5xkJ/v2R/5XXeXoKqJ0U/BNR54dIRD4xNO16znnr37eeTxg+k7Fsjs7i6WL+p1at8TFs6fcj97nzgYeb97nzjIPXv3p4b1hjltydGJ38e1zVOfMoclx8xPPDbu3oJ7D/jBPfHFttPkqxpxv00ncdPbXpzbeCS6rVR1QkS2i8giVb0/n3i56AMeCH1+EDgtvIOIXAJcAnDEsSeUJ1kbk9dwBMfes9czDHEv0L2/OhB7/PiEcsLC+ZMvY9w1jjliTkMvaNoCVrO7u5hb6+LXT05fVzyJLokP+w3ux6V9w21YbyDCbfzAoyOZDMfs7nRPclzbuCz6dfzRPdOMWZd42+vliGuHbfcPF6KAy1DqSb9NpxmQvLjMeRwL7BSRHwKT3UpVPa9pUhFZImjKq6Sq1wPXg+e2Wv/WFzVRnOoTnmfonVdD1SugF3Y1OZUMJ15RTig8OTpBXFsvXrMh9px9vT186z1nAMmujUZ+x4FtQ7E937B7ZcW6TanGo9YlHDF31mQI75knLpx0IcXJ7dK+4TZcsW7TtHYOvh/NYOi7BD564fNSXXtLYn6f0fH43zSMS5jvwLYhVn9pO6MRD9Ch8QkeGn6Sd73smbnnLgLXYNhgN3rOKJJ+m07SNTe9Lf+xLsbjyvynz82DwPGhz78FPNQCORqirOUp633tYZ9/2qp3Ubz+tEWxGdsPDY/kuq/6NSmKjo0P2iCO8KR70gRx1NyO66JFru3rspRwlryQpFUEw2QNSEjLy4klIeq40VyVssqBzOTMcVdc1vP4ThmC1PEj4JkisgQYAv4QeH0L5MhNmSVF0rJxo1a9i6sXBN5CTnGK5qie2rT7unT9IGtv3cn82d3sPzRdjvmzu6etSRGWowjDmtYGBw4dXiwpa22upEWLLjjVi0S6bP0gx/X2cMGpfbGT2gEuSwmvXrmUy9YPOhVajJq2zFNWpn70+sSTY5MjCNfn95qNu1NLtzSigJNqghVJJ5SLbzapjlIROV1EfiQiT4jIIREZF5FfN1MoVR0D3glsBH4C3KSqO5t5zaIps6SIy8sYXvQpiEmfSKgXFBdnL0KkIh0eGeXQ2ATdXVO7nbVu4cOvPnna/kXHxqe1QbBC4MC2ocx5M3HnDiqohpPubt7qnb+3J3pZWIEpI5U4OVYt7+Pi0xc5l3gPk6esTP0x+w6MTnM9uTy/Ls9iIwo4LpnSNcnSlWbmVnUKLm6rj+P1/L8E9ANvxKt11VRU9WvA15p9nWZR5rDXxcUR9cL2zqtFhrUe19sTOzrwkuiiGZ1QempdHFImw0svesHxToYhqqccdf24c7m0wcjoOO+5aTvXvvYUrj7/5IbPHSTC1V9j7a072X8oek7l4tMXOS8lfNWqk6fkfsyLGdm97rTDHt6kPJOksjKutaTSnt+036FRBRxXIDFL9WMXLM8kHRfjgar+TES6VXUc+CcR+X6T5Wp7yhz2pvnao8prDGwb4omISeNat0y+3FGVW9NqOoUT7JLWNggbi/qEuKHhEVZ/eTsozm4TVzfPuCqX37KDq88/2TnvIM7dE9feUbWxABbMq02uFx+QZanhgW1D/Pn6wWml379w5/1s3rWHM09cOGVtiXqSFL9rpybt+Y1qq2AZ37S1XlzoS3A5Fk1ZlYvbFZdM8QMiMhsYFJGPishlQHJQuFHqsLe+VMa82tSfNaq8xjUbd0dGxMyfPSvxhclaCiTK1RFVEqXeTz46rpncJquW9znND6SdJ+7cUe6erArLJSQ2iWs27o5dM2RoeIQb7rw/89LALt8FuDy/UW113UXLuLcg96S5k6qDy8jjj/CMzDuBy/CioC5oplCdQNkL9oSvc+DQGAfqSmyMjI5z5W07J/eLU7TB+thp62u//5a7pl0jjvperauLxOVcYeJ6pVnPA+6RRlEjkrhsehcFnRTJliZzkvFMU7CrVy5l9Ze3TzHi3V3CU+bMmhbynUYze+zmTqoOLtFW94lID3CsqrYibLdtaeZLFCiZoeGRSbcAJEed7Dswmlq6Iyi8lxQpFtQ6OpDT1dHIvE+SAnYNlU07j2uk3KrlfWy579EpJUQuONWrUZUnFDntunnLuneLuJV4qbM+XTBZd6tKmDupGrhEW50LDAJf9z8vE5FbmyyXkUDY7QPJPc4sBAqukQq4cecM08i8T/3cDRwuKX/Z+kHmzOpiwbzalEq3Wd0crpFycetXQ77FstKum6d6sADXvvYUp2vXuwlHJ7TpC44Z7YuL22ot8ELg2wCqOigii5sn0szFNfmuEbdPFPWJcXERVWGD4dILjpogHdg2xIGISKRwRneSMdy8a8+Uz/W99eGRUXpq3Vx30bLJ64YjllzcHEkVgJddefukG2f/wbFIZX/p+sFck8NpEXphl039iDMO1wTCZkcHlpUwa5SHi/EYU9XHpOA4amMqUS6Ly9YPsuW+R6dF6Li80L09NebPmcXQ8EhiXaaoxDiXSDGXCK/688ZVnQ2XJYf48iVw+N7Dbrt6okqFZ1FUSYYxiKRKM5x5kkKzROgJXsKmiDcRH5fwGZdv0si1s1JmwqxRHi7RVneLyOuBbhF5poh8DLBQ3YKJGk0ocMOd909bTCfthe6pdbP2vJNYvXIptW6JNRxx7huXiJYgqiYuOStKxrgR0/w5UyO8Vq9cGpscF56TSVLgYQObdfGkRheXCsga1RV13XCYdVSUWrAK4rWvPYVa1/RW239oLPV+465dVBRTmQmzRnm4GI93AScBB4F/BR4DLm2iTDOSuNGEwrSXLE7JwFT/elKpiKRJ1LSlNcNzDEf2zJqmtKKUzsC2IedlbeOyq5PmZOoJjFeepVeD+y+C8EgpzYCtWt7HBaf2TbnvcJh1Wl2nI+ZOdySMjrvNWzRzOVWrE9WZxLqtRGQu8Dbgt4EdwIv8siFGwQxsG4p1O8DUlyy8QFBAnH896eWcUHVOTquXtb4IY61b6O2pxYZ0ZilaGFCfXe0yJxMQNl55C+kFxrfRmkku0WthNu/aM20eI5A3aS5mYNtQbB6Jq5JuVhST1YnqTJLmPD4PjAL/AZwNPBsbcRROUsXWgOAl+8DADr5w5/RlVc48cWHkS5/ku8/74kYp49FxZf6cWQxecZbzMQFJrpE4ZZZ0X73+PMBl6wcTlX+UQq2f1A0ytvMGJ9S6JDV6zdXgp1XaTcquL0pJ5530bkYVZaP1JLmtnqOqb1DVTwEXAi8tSaYZRZoLJvyS3fiDByL3idsezHnUEyi1PORxQSR9l8c1Eueff8Ppizg4NsE+P2IriEiKol6hRrm3bt46xAWn9k26chbMq9FTc/H0ehwx15vLydJmcYo+UNZxczFxhqO+am6WuZ8wedx/Ac10iRmtI2nkMTkGVtUxi7ZqDmkZ0+HeXdaicMFxV962czI5sD6yKSt5XBBJJdDzyFEfshoUKIxarElhWkhrVK83bnSwedeeyaixNPdbPfsOjDKwbSi1zerrfNW6ZcpcVbjSLsClKW67gPDz02jEU6PraFhiX+eRZDxOCZVeF6DH/yyAquqRTZduBhCnWIIQy8D9snrlUrpj5kWSylEX/dLmcUHEFcuLSviLIsldEj5vnBENivIluVtcRgd58msuv2UHF5zaN839VesSDhwaY/GaDVOM2/DIKLUuYcG82uQqhmF5Xedi6kOlG1X+Nult1BNrPFS18VhFI5UoxVrrEvYfGpuSU3D5LTs4/RkLuOPn05dZDZfkbjZ5agsFZTzCqxMGUUThirtxZdnjesyuyjxukacwLiOqJEVZP1oICEYv4RLwQRXhYDRYf9TohDJv9iy2/dXUOaS4kjRx9xOmUeVvk95GPU4l2Y3mEaWMD4QUS8DI6Dj3/mqEN5y+aEotpdeddvy0JEJXskyANpohnBRFFOdWuWz9IHNrXVPKvIePc1F8rhOzZ564MDIY4dH9B51WIFy9cmmsO+mh4ZEpI8AV6zbFlm0PHxOmvn3C7jgXt1yjyr/Vk96WoV49zHhUgHrX0pI1GyL3GxoeYfOuPU61itLI4gMvIkM4Kcx0yZoNkaHKCtMMR/h8R/XUIpVwtwgTqpmUTH3Zk4CR0YnENeDD8xFx7qQ8hSHrj4lLIl0wr8YV556UqljjXIdDwyOsWLfJafQYyFG2ArcM9WpixqOCJIVkFvXiZPGBN+ovh+R7UrKvBBe4fuqpdQnXvCa7cU1S6HFrwNcrUNfeeZ7V9uLkC0aoaW65pLpYrs9Uqya9i3j+jOJxjzs0SiOtPEYRpR2y+MCLmCwtquQHHF5LPWqOYfasLq7ZuDtzOGqa+yZqDfj6xY2CDPEggCEo0R41CogKoQZvJBEVxpok36XrByfvNSkcN5C9r7cn1oVYRWyyvpqY8agg4bj4OBp9cZLyCRrZN456xZqXoKxKXDb1/kPjuXIR0oyb60JOUSXa66+/ankf82dHD/rnxazkmDa3MDQ8wuovbWf1l7en3n9WZdxIfkgRFPH8GcVjxqOihHuJUTT64mQphFdE0bx6xZqXoKyK6/279qgD4xZRW5CeWjdnnrgwVYFmKQD4WMyEeSOdgtEJnTYai7p+FmXcSHJgUdjSs9XEjEfFadaLkyXrt4gMYZew2iBLPGnE1TvPy3/J4gZzUciBcauvQNzbU5vM0yiyR++iwMM9/vfctD31HuKov36WZ6oKFXEtQ72a2IR5xWlmlEuWCdBGJ0uTFHj9YlQAy668PTKSKhi4RLXL/oNjkce4jFKSysVv3rXHacI2Szhs2uR6fYRRIyO2+utneaaqMt9gGerVw4xHG9AJL05SjkRUpFCcWye8vb5dBrYNsfpL26csp+pax6uIml1xBiFweUUp6jgFniebvdYlUBdIEDeicH2mLDnQiMOMh1EKq1cuzaTYcyut+jkLx/n5tOu51qbqnVdjzqyuyfL0i3+jZ0pmfX1YbJwCd+nZ17qF+bNnTSmFD8WOUludHGhUFzMebUDVsmtzy5NBsedRWlGLXwWLIaXJl3a9uO+i1jcJ1lCH6FLpLjkKccbMJQGyyGcjr9u0vthjsFxuFZ5foxjMeFSIrLWdWvEC5s32zarY8yitRvzzcdcLto2Mjk8WpgxXq12xblPihHLcTEWaTHHGrBUTxXEjpLhORP0zEp6HavXzaxSHGY+KEKeU58zqqlR2bd5s3zyKPetcT6P++ag5lPpJ6/ry6HkNVhA1Fr5WvSIOF1OsWo89qRORNl9j2eGdgYXqVoQ4pRxXQK9V2bV5lWUZiV5FhzW7hKkm3VfSvYWDp+JyKYDYbPZWk9Q2Ls+mZYe3P2Y8KkLWl6lV0S5H9dQybQ8oI9Gr6HwAF0OZdF9J9zY8Mjolh6PVuRRZSVsuN41OiNZqdeZ9qzG3VUWIc7ksmFfjydGJykS7xFUXSas6EqzpES4nH1X3qVGKDGt2cYOlzc2svXVn5OgxqGgL8TkcVe6dJ7VN1HxNmE6I1rJKvxUceYjIWhEZEpFB/985rZapDOJ6sFece1KlsmvjakrFbQ9wrfvUbLL0Fl1HS1HFEoPrxLkdXVL+qtw7T2qb+hFgb0+NBfNqlXh+i6IKmfetpqojj+tU9W9aLUSZpPVgq/Ky5Z2UrkJZ7ay9xUbCVJN63i7k6Z2XubiXy/NalWe2GVQl876VVNV4zEja4YXLmzSWtBhUXPZ10eQxYHl+kzzZ4ZBvEauAshf3Cvat+vPaLCzzvoJuK593ishdIvJZEVkQtYOIXCIiW0Rky5490avAGcWTd1I67qUKfP9lVGxttLfo6vJyOV/9FFFPrZtrX3tK7siqLG4Uc7k0jlX6bdHIQ0S+CfxmxFd/CXwC+BCeW/hDwLXAn9TvqKrXA9cD9Pf3N1bn28hEnh5n3DKoebKv89JIbzFLb91lpcALTu1j8649hY24yl7ca6bTzIKl7UJLjIeqvtxlPxH5NPDVJotjFEySPz28PU7BNkuJNVKnKYvLK8lQ9jVJyWQxjOZyKYaZ7LaDCs55iMixqvqw//HVwN2tlMfIRloPPfyyrVi3qVQl1khvMUtvvRW90iwVfa3YoVEEog2u7FY0IvIvwDK8jtq9wFtDxiSS/v5+3bJlS/OFM1KJMwhRpdejopKa3UPPSjCKihslxZWUjzpHsw3JBwZ2TMmjOf0ZC/iv+x+LrI8FM9vlYniIyFZV7c9zbOVGHqr6R62WwchP3h760PDIlDmQIpKuGlXaaSG3Lr31spLJovJovv/zR2PnlKpW7sRoP6oabWW0KVlrWAUJdn29PbGKLg8fGNjBZesHG4rkSgq5dY0yKyuyKeo6cT6FoeGRGVdKwygeMx5GoeQNYSxy8nxg29CUBZgCsirtuGsLOPfcy4psynq+ZoZEGzMDMx5GoeTJAxnYNhS7LlSeyfNrNu7OvY6Gy7WzyFRGNeGk88W1q+V1GI1ixsMonKhaT0lcedvOSGUvkCsCKMlAZFHaRSSClZVMFnedi09fFHuM5XUYjVC5CXOj/claY2lfTFFFJd+kclweQ1ZjVETIbVlhu0nX2bxrT8Mh0VVbCtloPZUL1c2DhepWh6gIpaTlU+NCe8EtDNZVBgEuPn0RV606OfP5iqZsRZz1Nyn6eKO6NBKqa24ro1CyRhcluU7yunai5l2uu2hZZQxH1KqBzZy8bnSRLKuFZURhbiujULJGF8W5mHp7ag31asPZ7EFP/7L1gy13ubSqNH0jpTSsFpYRhY08jELJGl0UN9G79ryTCpGnFT39JNpREZcVMWa0F2Y8jELJGl1U9Lrj9VTN5dKOitjKjxtRmNvKKJQ80UXNrE5atZ5+OxYltPLjRhRmPIzCqVKp6izlx8uIgmpXRVyl39SoBmY8jI7GtadfVgHD4HymiI12x+Y8jI7GdU6lanMjhlF1bORhdDwuPf2qzY0YRtWxkYdh0J5RUIbRSsx4GAadGY46sG2IFes2sWTNBlas22Ql2I1CMbeVYdC+UVBxlBkAYMxMzHgYhk8nRUG1qgyKMXMwt5VhdCAWAGA0GzMehtGBWACA0WzMeBhGB9KJAQBGtbA5D2NG06kr5LUyAKBT29SYihkPY8bS6RFJrQgA6PQ2NQ5jbisjkU7OFbCSJMVjbTpzsJGHEUun9yItIql4rE1nDjbyMGLp9F6kRSQVj7XpzMGMhxFLp/ciLSKpeKxNZw7mtjJiybKQUjvSaSVJqoC16cxBVLXVMjRMf3+/btmypdVidBz1cx7g9SKLXGPcMIzWISJbVbU/z7E28jBisV6kYRhxtMR4iMhrgLXAs4EXquqW0HeXA28BxoF3q+rGVshoeHRSsUDDMIqjVSOPu4HzgU+FN4rIc4A/BE4CjgO+KSLPUtXx6acwDMMwWkVLoq1U9SeqGhXv+Srgi6p6UFXvAX4GvLBc6QzDMIw0qhaq2wc8EPr8oL9tGiJyiYhsEZEte/bsKUU4wzAMw6NpbisR+SbwmxFf/aWqfiXusIhtkeFgqno9cD140Va5hDQMwzBy0TTjoaovz3HYg8Dxoc+/BTxUjESGYRhGUVQtVPdW4F9F5G/xJsyfCfww7aCtW7c+ISKdUTOjcY4B9rZaiIpgbXEYa4vDWFscJnfqf6tCdV8NfAxYCGwQkUFVXamqO0XkJuDHwBjwDsdIq915E106DRHZYm3hYW1xGGuLw1hbHEZEcmdXt8R4qOq/Af8W892HgQ+XK5FhGIaRhapFWxmGYRhtQKcYj+tbLUCFsLY4jLXFYawtDmNtcZjcbdERhRENwzCMcumUkYdhGIZRImY8DMMwjMy0lfEQkd8Xkd0i8jMRWRPxvYjI3/vf3yUiz2+FnGXg0BYX+21wl4h8X0ROaYWcZZDWFqH9XiAi4yJyYZnylYlLW4jIGSIyKCI7ReQ7ZctYFg7vyFEicpuIbPfb4s2tkLPZiMhnReQREbk75vt8elNV2+If0A38HHgGMBvYDjynbp9zgH/HK3NyOvCDVsvdwrZ4MbDA//vsmdwWof02AV8DLmy13C18Lnrx8qgW+Z+f2mq5W9gW7wc+4v+9EHgUmN1q2ZvQFi8Fng/cHfN9Lr3ZTiOPFwI/U9VfqOoh4It4VXjDvAr4Z/W4E+gVkWPLFrQEUttCVb+vqvv8j3filXrpRFyeC4B3ATcDj5QpXMm4tMXrgVtU9X4AVe3U9nBpCwWeIiICHIFnPMbKFbP5qOp38e4tjlx6s52Mh0vFXeeqvG1O1vt8C17PohNJbQsR6QNeDXyyRLlagctz8SxggYh8W0S2isgbS5OuXFza4uN4C9I9BOwA/reqTpQjXqXIpTerVtsqCZeKu85Vedsc5/sUkTPxjMdLmipR63Bpi78D3qeq414ns2NxaYtZwKnA7wE9wH+KyJ2q+t/NFq5kXNpiJTAIvAw4AfiGiPyHqv66ybJVjVx6s52Mh0vF3ZlSldfpPkXkecA/Amer6q9Kkq1sXNqiH/iibziOAc4RkTFVHShFwvJwfUf2qup+YL+IfBc4Beg04+HSFm8G1qnn+P+ZiNwDnIhDMdYOI5febCe31Y+AZ4rIEhGZjbdc7a11+9wKvNGPHjgdeExVHy5b0BJIbQsRWQTcAvxRB/Yqw6S2haouUdXFqroY+DLw9g40HOD2jnwF+B0RmSUi84DTgJ+ULGcZuLTF/XgjMETkaXgVZn9RqpTVIJfebJuRh6qOicg7gY14kRSfVa8K79v87z+JF0lzDt7ytQfwehYdh2Nb/BXwG8A/+D3uMe3ASqKObTEjcGkLVf2JiHwduAuYAP5RVSNDONsZx+fiQ8DnRGQHnuvmfaracaXaReRG4AzgGBF5ELgCqEFjetPKkxiGYRiZaSe3lWEYhlERzHgYhmEYmTHjYRiGYWTGjIdhGIaRGTMehmEYRmbMeBgzAr+a7mDo32IR+b7/3WIReX1o32Uick6Oa3xbRBoOhy7qPIbRTMx4GDOFEVVdFvp3r6q+2P9uMV7BwIBleHHvhmHEYMbDmLGIyBP+n+vwsq4HReR9wAeBi/zPF4nIfH9NhB+JyDYReZV/fI+IfNFfA2E9Xq2o+mucLSI3hT6fISK3+X9/QkS2+GtJXJkiIyJyoYh8zv97oYjc7Mv0IxFZ4W//3dDoapuIPKWItjKMetomw9wwGqRHRAb9v+9R1VeHvlsDvFdV/wBARH4J9KvqO/3Pfw1sUtU/EZFe4Ici8k3grcABVX2eX0fsvyKu+w3gUyIy368ndRGw3v/uL1X1URHpBr4lIs9T1bsc7+f/Atep6vf8UjQb8SrEvhd4h6reISJHAE86ns8wMmHGw5gpjKjqspzHngWcJyLv9T/PBRbhLbLz9wCqepeITFP8fpmMrwPnisiXgVcCf+F//VoRuQTvPTwWeA5e2RAXXg48J1Ql+Eh/lHEH8LcicgPeuh0PZrtVw3DDjIdhpCPABaq6e8pGT3G71PdZD7wDb0GeH6nq4yKyBG+U8AJV3ee7o+ZGHBs+f/j7LuBFqjpSt/86EdmAN2dzp4i8XFV3OchoGJmwOQ/DgMeBpyR83gi8S3xrISLL/e3fBS72tz0XeF7M+b+Ntwzon3HYZXUksB94zK/oenbMsb8UkWeLSBfeglYBtwPvDD6IyDL//xNUdYeqfgTYgldi3DAKx4yHYXiuojER2S4ilwGb8VxCgyJyEV711Rpwl4jc7X8G+ARwhO+u+gti1oFQ1XHgq3gG4qv+tu3ANmAn8Fk8d1MUa/xjNgHhMtnvBvr9yfofA2/zt18qIneLyHZghM5dQdJoMVZV1zAMw8iMjTwMwzCMzJjxMAzDMDJjxsMwDMPIjBkPwzAMIzNmPAzDMIzMmPEwDMMwMmPGwzAMw8jM/wewUhgyyXKe9wAAAABJRU5ErkJggg==\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+TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYdUlEQVR4nO3de5hcdX3H8ffHEOUuYBaMBFxFoCKt4XGNWNAilxoBDVipQKVJRUOrKFiqBm+FeouPCrRVsUF4CFdFBUEQJQbTiCK60YCJQVGJ3EKyXGISqWjIt3+c38BhMrNzdndmZ3/Zz+t55tlzP99z9sxnfnPOmRlFBGZmlp9ndLsAMzMbHge4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOCJpOWSDul2Hd0k6VhJ90raIOmAbtdTI+ksSZe1cXkXS/pY6n6VpF+2a9mldYSkF1Wcti3bJ+mLkj480uWMRYM9PyUdIum+Nq1nkaS3tWNZo2FcBLiklZIOrxs2S9Ittf6IeElELGqxnN70xNyqQ6V222eAUyNi+4j4WdWZyoGYm4j4fkTs2+062iEi/jkiPtrtOjqhyvNzPBoXAZ6LMfDC8HxgeZdraCtJE7pdgxXGwPG9xXGAJ+VWuqRpkvolrZO0WtI5abLF6e/adJrhlZKeIelDkn4naY2kSyQ9u7Tcf0zjHpb04br1nCXpa5Iuk7QOmJXWfauktZJWSfqcpGeWlheS3iHpLknrJX1U0l5pnnWSripPX7eNDWuV9CxJG4AJwO2SftNgXkk6N833e0l3SNpf0mzgH4D3pX3yzTT9HEm/STX+QtKxpWXNknSLpM9IelTS3ZJeVxr/Akn/m+ZdAEyqq+Wrkh5MdSyW9JLSuIslnS/pW5L+ALxG0gGSfpqW9xVg69L0T779lvTmtA21x+OSFqVxz0r13pOOiS9K2qa0nPem/9cDkt7aaP8PYfsOlPTDdAzcrnTqQNLxkvrrpn2PpOtK2147NbSzpOslDaR9fL2kKaX5FqVj5wepjpskTSqNP7hUw72SZlXZD3W1zUrLP1fSI8BZg80vaVKqc62kRyR9X9Iz0rjy82abtK2PSvoF8PK69T7t9NVQ9kvdcl6U/k+/l/RQOnbGlojY4h/ASuDwumGzgFsaTQPcCpyUurcHDkzdvUAAW5Xmeyvwa+CFadqrgUvTuP2ADcDBwDMpTlH8ubSes1L/MRQvptsALwMOBLZK61sBnF5aXwDXATsCLwEeBxam9T8b+AUws8l+aFpradkvajLva4ElwE6AgBcDk9O4i4GP1U1/HPC8tF1vBv5Qmn5W2u63U7xo/AvwAKDS/j8HeBbwamA9cFndduyQxp8HLC2Nuxj4PXBQWveOwO+A9wATgTeldX8sTX8IcF+D7d0x7ftTUv95ab/vktb9TeCTadx0YDWwP7AdcEWLfdl0+4DdgYeBI1P9R6T+HmDbNO3epWX9BDi+/v8APAf4uzTPDsBXgW+U5lsE/AbYh+K4WwTMTeP2TOs5Ie2z5wBTW+2HBts5C9gIvIvieN6mxX78JPDFtM6JwKtKx8RKnnrezAW+n5axB7Cs/D+s3/fD2C9vS91XAh9M/4etgYO7nWWb7eNuFzAqG1n88zcAa0uPx2ge4IuBs4FJdcvpZfMAXwi8o9S/L0VAbAV8BLiyNG5b4E88PcAXt6j9dOCauoPzoFL/EuD9pf7PAuc1WVbTWkvLbhY6hwK/onhxeUbduCefIINsx1JgRuqeBfy6br8E8FyK8NgIbFcafwWlAK9b7k5p3meXarmkNP7VlF4c0rAfMkiApyfs9cD5qV8UL0B7laZ5JXB36r6IFH6pf59m+7LV9gHvp/SimoZ9h/SiDFwGfCR1700RtNu2+j8AU4FHS/2LgA+V+t8BfDt1n1k+5krTDLofGkw/C7in6vzAfwDXNtlvK3nqefNbYHpp3GwqBnjF/VIL8EuAecCUwY7tbj7G0ymUYyJip9qD4oBt5mSKJ+Gdkn4i6ehBpn0eRQuv5ncU4b1bGndvbUREPEbRmiq7t9wjaZ/0tu5BFadVPkHdW2yK1l7N/zXo334YtQ4qIm4GPgd8HlgtaZ6kHZtNr+LU0dL0dngtReu0vB0Plpb9WOrcPtX4aET8oa7O2nInSJqr4vTMOoonNnXLLu/T5wH3R3pG1i+viY9TtM7enfprrd8lpe35dhpeW0d5nYMtf9Dto7gOcVxtPWldBwOT0/grKFrGACdStB4fo46kbSX9j4rTZesoGiU76enXBB4sdT/GU8fNHhSt83qt9kMj5f3Sav5PU7xDvEnSbyXNabLMoezvp6m4X2reR/Gi82MVd8EMemqsG8ZTgFcWEXdFxAnArsCngK9J2o7ilb3eAxRPuppaC2s1sAoon3fchuIt3NNWV9d/PnAnxdvkHYEPUBxE7TBYrS1FxH9FxMsoTt3sA7y3Nqo8naTnAxcApwLPSS+Yy6i2HauAndP+LtdZcyIwAzic4pRRb2215VLrlre7pPL48vKeRtLxFAH5poj4cxr8EMUL40tKjYBnR0Qt8FZRhF7L5VfYvnspWuA7lR7bRcTcNP4mYJKkqanOK5qs5wyKd1ivSMfRq2ubOEht5Rr2ajC81X5opPy/GHT+iFgfEWdExAuB1wP/KumwBststb8fo3ihqHluqbvyfomIByPi7RHxPOAU4AuqeGvoaHGANyDpLZJ6ImITxekWgCeAAWATxTnkmiuB96i4MLU9RYv5KxGxEfga8HpJf63iwuLZtH4C7QCsAzZI+guK88PtMlitg5L0ckmvkDSR4m3wHyn2CRQvAOV9UnuxG0jz/hNFC7yliPgd0A+cLemZkg6meDLX7EBx3v9hiifpJ1os8laKF6l3S9pK0huBaU228QDgvynerQ2UatpE8YJ0rqRd07S7S3ptmuQqigvQ+0naFvj3EWzfZRTHzGvTu42tVVxonZLmrx1Xn6Y4B7ygyap2oAjLtZJ2GaymBi4HDpf092mfPUfS1Ar7YVCt5pd0dLpwKIrnwBM8dYyVXQWcmS5ITqE4x162FDgx7b/pwN+UxlXeL5KOK13gfJTimG5UT9c4wBubDixXcWfGf1JcJPpjeqv6ceAH6S3ggRTnPy+leCt2N0WwvQsgIpan7i9TtBrWA2soAqiZf6NoZa6nONjbeeW7aa0V7JjqeZTiLevDFBdlAS4E9kv75BsR8QuKc/G3UoT7XwI/GEKdJwKvAB6heIJdUhp3SVr//RQXbH802IIi4k/AGynOxz5KcUH16iaTzwB2Bm7RU3ei3JjGvZ/i7f2P0lvv71K05IiIGykuzt2cprl5uNsXEfemOj5A8QJ4L8U7nfJz9QqKdyBfHeTF9zyKi4YPUeyjb7eo6UkRcQ/FRdQzUo1LgZem0U33Q0WDzb936t9Acex8IRrf+302xTFwN8U7kkvrxp9G8aK4luIOqW+Uxp1H9f3ycuC2lAPXAadFxN2tN3H01K7w2ihIrd61FKdHxtSBYGb5cQu8wyS9Pl042Y6ixfpznrrwZmY2bA7wzptBcfHwAYq3iMeH3/aYWRv4FIqZWabcAjczy9SofrnMpEmTore3dzRXaWaWvSVLljwUEZt9YGpUA7y3t5f+/v7WE5qZ2ZMkNfy0qU+hmJllygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpYpB7iZWaYc4GZmmXKAm5llalQ/iWnjT++cG4Y0/cq5R3WoErMtj1vgZmaZqhzg6fflfibp+tS/i6QFku5Kf3fuXJlmZlZvKC3w04AVpf45wMKI2BtYmPrNzGyUVArw9MvMRwFfKg2eAcxP3fOBY9pamZmZDapqC/w84H3AptKw3SJiFUD6u2ujGSXNltQvqX9gYGAktZqZWUnLAJd0NLAmIpYMZwURMS8i+iKir6dns+8jNzOzYapyG+FBwBskHQlsDewo6TJgtaTJEbFK0mRgTScLNTOzp2vZAo+IMyNiSkT0AscDN0fEW4DrgJlpspnAtR2r0szMNjOS+8DnAkdIugs4IvWbmdkoGdInMSNiEbAodT8MHNb+kszMrAp/EtPMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8tUlR813lrSjyXdLmm5pLPT8LMk3S9paXoc2flyzcyspsov8jwOHBoRGyRNBG6RdGMad25EfKZz5ZmZWTMtAzwiAtiQeiemR3SyKDMza63SOXBJEyQtBdYACyLitjTqVEl3SLpI0s5N5p0tqV9S/8DAQHuqNjOzagEeEU9ExFRgCjBN0v7A+cBewFRgFfDZJvPOi4i+iOjr6elpS9FmZjbEu1AiYi3Fr9JPj4jVKdg3ARcA09pfnpmZNVPlLpQeSTul7m2Aw4E7JU0uTXYssKwjFZqZWUNV7kKZDMyXNIEi8K+KiOslXSppKsUFzZXAKR2r0szMNlPlLpQ7gAMaDD+pIxWZmVkl/iSmmVmmHOBmZplygJuZZcoBbmaWqSp3oZg9qXfODd0uwcwSt8DNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlT2KOY/5UpVne3AI3M8tUlZ9U21rSjyXdLmm5pLPT8F0kLZB0V/rb8FfpzcysM6q0wB8HDo2Il1L8Av10SQcCc4CFEbE3sDD1m5nZKGkZ4FHYkHonpkcAM4D5afh84JhOFGhmZo1VOgcuaYKkpcAaYEFE3AbsFhGrANLfXZvMO1tSv6T+gYGBNpVtZmaVAjwinoiIqcAUYJqk/auuICLmRURfRPT19PQMs0wzM6s3pLtQImItsAiYDqyWNBkg/V3T7uLMzKy5Kneh9EjaKXVvAxwO3AlcB8xMk80Eru1QjWZm1kCVD/JMBuZLmkAR+FdFxPWSbgWuknQycA9wXAfrNDOzOi0DPCLuAA5oMPxh4LBOFGVmZq35k5hmZplygJuZZcoBbmaWKQe4mVmm/HWyNqYM9StuV849qkOVmI19boGbmWXKAW5mlikHuJlZphzgZmaZ8kXMLYh/49JsfHEL3MwsUw5wM7NMOcDNzDLlADczy5QvYlrW/MlNG8/cAjczy1SVn1TbQ9L3JK2QtFzSaWn4WZLul7Q0PY7sfLlmZlZT5RTKRuCMiPippB2AJZIWpHHnRsRnOleemZk1U+Un1VYBq1L3ekkrgN07XZiZmQ1uSOfAJfVS/D7mbWnQqZLukHSRpJ3bXZyZmTVXOcAlbQ98HTg9ItYB5wN7AVMpWuifbTLfbEn9kvoHBgZGXrGZmQEVA1zSRIrwvjwirgaIiNUR8UREbAIuAKY1mjci5kVEX0T09fT0tKtuM7Nxr8pdKAIuBFZExDml4ZNLkx0LLGt/eWZm1kyVu1AOAk4Cfi5paRr2AeAESVOBAFYCp3SgPjMza6LKXSi3AGow6lvtL8fMzKryJzHNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLl7wO3ccXfH25bErfAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFMOcDOzTDnAzcwy5QA3M8uUA9zMLFNVfhNzD0nfk7RC0nJJp6Xhu0haIOmu9HfnzpdrZmY1VVrgG4EzIuLFwIHAOyXtB8wBFkbE3sDC1G9mZqOkZYBHxKqI+GnqXg+sAHYHZgDz02TzgWM6VKOZmTUwpK+TldQLHADcBuwWEaugCHlJuzaZZzYwG2DPPfccUbHjzVC/+tTMxpfKFzElbQ98HTg9ItZVnS8i5kVEX0T09fT0DKdGMzNroFKAS5pIEd6XR8TVafBqSZPT+MnAms6UaGZmjVS5C0XAhcCKiDinNOo6YGbqnglc2/7yzMysmSrnwA8CTgJ+LmlpGvYBYC5wlaSTgXuA4zpSoZmZNdQywCPiFkBNRh/W3nLMzKwqfxLTzCxTDnAzs0w5wM3MMuUANzPLlAPczCxTDnAzs0w5wM3MMuUANzPLlAPczCxTQ/o6WbPxZjhf6bty7lEdqMRsc26Bm5llygFuZpYpB7iZWaYc4GZmmXKAm5llygFuZpapKj+pdpGkNZKWlYadJel+SUvT48jOlmlmZvWqtMAvBqY3GH5uRExNj2+1tywzM2ulZYBHxGLgkVGoxczMhmAk58BPlXRHOsWyc7OJJM2W1C+pf2BgYASrMzOzsuEG+PnAXsBUYBXw2WYTRsS8iOiLiL6enp5hrs7MzOoNK8AjYnVEPBERm4ALgGntLcvMzFoZVoBLmlzqPRZY1mxaMzPrjJbfRijpSuAQYJKk+4B/Bw6RNBUIYCVwSudKNDOzRloGeESc0GDwhR2oxczMhsCfxDQzy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUw5wM7NMOcDNzDLlADczy5QD3MwsUy2/Ttaa651zQ7dLsDFoqMfFyrlHdagS29K5BW5mlqmWAZ5+dX6NpGWlYbtIWiDprvS36a/Sm5lZZ1RpgV8MTK8bNgdYGBF7AwtTv5mZjaKWAR4Ri4FH6gbPAOan7vnAMe0ty8zMWhnuOfDdImIVQPq7a7MJJc2W1C+pf2BgYJirMzOzeh2/iBkR8yKiLyL6enp6Or06M7NxY7gBvlrSZID0d037SjIzsyqGG+DXATNT90zg2vaUY2ZmVVW5jfBK4FZgX0n3SToZmAscIeku4IjUb2Zmo6jlJzEj4oQmow5rcy1mZjYE/iSmmVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOBmZplygJuZZcoBbmaWqZY/6DAYSSuB9cATwMaI6GtHUWZm1tqIAjx5TUQ81IblmJnZEPgUiplZpkYa4AHcJGmJpNntKMjMzKoZ6SmUgyLiAUm7Agsk3RkRi8sTpGCfDbDnnnuOcHVmZlYzohZ4RDyQ/q4BrgGmNZhmXkT0RURfT0/PSFZnZmYlww5wSdtJ2qHWDfwtsKxdhZmZ2eBGcgplN+AaSbXlXBER325LVWZm1tKwAzwifgu8tI21mJnZELTjPvAtRu+cG7pdgo1DQz3uVs49qkOVWG58H7iZWaYc4GZmmXKAm5llygFuZpapbC5i+kKP2fD4ubPlcgvczCxTDnAzs0w5wM3MMuUANzPLVDYXMYfKn6q0LdVYO7aHU48vlLaHW+BmZplygJuZZcoBbmaWKQe4mVmmttiLmGY2PKNxkbTTnw4da8sfzjqqcAvczCxTIwpwSdMl/VLSryXNaVdRZmbW2kh+1HgC8HngdcB+wAmS9mtXYWZmNriRtMCnAb+OiN9GxJ+ALwMz2lOWmZm1MpKLmLsD95b67wNeUT+RpNnA7NS7QdIvR7DOTpoEPNTtIobINY+O3GrOrV5oUbM+1dmVD3P5Q9rPI9yG5zcaOJIAV4NhsdmAiHnAvBGsZ1RI6o+Ivm7XMRSueXTkVnNu9YJrHq6RnEK5D9ij1D8FeGBk5ZiZWVUjCfCfAHtLeoGkZwLHA9e1pywzM2tl2KdQImKjpFOB7wATgIsiYnnbKht9Y/40TwOueXTkVnNu9YJrHhZFbHba2szMMuBPYpqZZcoBbmaWKQd4iaRPS7pT0h2SrpG0U7drakXScZKWS9okaczehpXb1y5IukjSGknLul1LVZL2kPQ9SSvSMXFat2sajKStJf1Y0u2p3rO7XVNVkiZI+pmk67tZhwP86RYA+0fEXwG/As7scj1VLAPeCCzudiHNZPq1CxcD07tdxBBtBM6IiBcDBwLvHOP7+XHg0Ih4KTAVmC7pwO6WVNlpwIpuF+EAL4mImyJiY+r9EcW97WNaRKyIiLH66daa7L52ISIWA490u46hiIhVEfHT1L2eImB2725VzUVhQ+qdmB5j/q4KSVOAo4AvdbsWB3hzbwVu7HYRW4hGX7swZoNlSyCpFzgAuK3LpQwqnYpYCqwBFkTEmK43OQ94H7Cpy3WMvx90kPRd4LkNRn0wIq5N03yQ4u3o5aNZWzNVah7jKn3tgrWHpO2BrwOnR8S6btczmIh4ApiarjddI2n/iBiz1x0kHQ2siYglkg7pcjnjL8Aj4vDBxkuaCRwNHBZj5Cb5VjVnwF+7MEokTaQI78sj4upu11NVRKyVtIjiusOYDXDgIOANko4EtgZ2lHRZRLylG8X4FEqJpOnA+4E3RMRj3a5nC+KvXRgFkgRcCKyIiHO6XU8rknpqd3pJ2gY4HLizq0W1EBFnRsSUiOilOI5v7lZ4gwO83ueAHYAFkpZK+mK3C2pF0rGS7gNeCdwg6TvdrqleujBc+9qFFcBVY/1rFyRdCdwK7CvpPkknd7umCg4CTgIOTcfv0tRSHKsmA9+TdAfFi/yCiOjqbXm58Ufpzcwy5Ra4mVmmHOBmZplygJuZZcoBbmaWKQe4mVmmHOBmZplygJuZZer/AUeVqAoD5LCwAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApbUlEQVR4nO3dd3gU5drH8e9NaCIoElA5agARGx6xYC/HfmxHRUVFFASPCBHBglLynqNHDSAgig2IFAVWsGBFUbH3AooKoghKELFAEEEBE5Pn/WMmuGQ3ySRkdjfJ73NduXZndmb2Xsve+5S5H3POISIiEq1OsgMQEZHUo+QgIiIxlBxERCSGkoOIiMRQchARkRh1kx1AVWjevLlr3bp1ssMQEalW5s2bt9o51yLeazUiObRu3Zq5c+cmOwwRkWrFzHJLe03dSiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIxlBxERGqoSARat4Y6dbzHSCT4uTViKquIiGwpEoFevWDDBm87N9fbBujatfzz1XIQEamBsrL+SgzFNmzw9geR1ORgZpPM7GczWxC172Yz+97M5vt/pyczRhGR6mj58ortLynZLYcHgVPj7L/TOXeA//d8gmMSEan2MjIqtr+kpCYH59ybwJpkxiAiUhNlZ0OjRlvua9TI2x9EslsOpelrZp/53U47xDvAzHqZ2Vwzm7tq1apExyciktK6doWcHGjVCsy8x5ycYIPRAJbsNaTNrDUwyzm3n7+9E7AacMCtQEvnXM+yrtGxY0enwnsiIhVjZvOccx3jvZZyLQfn3E/OuULnXBHwAHBosmMSEaltUi45mFnLqM1OwILSjhURkXAk9SY4M5sOHAc0N7MVwE3AcWZ2AF630jLgymTFJyJSWyU1OTjnusTZPTHhgYiIyBZSrltJRESST8lBRERiKDmIiEgMJQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOIiIpbmvWgq4srSEtIpLCtnYt6MpSy0FEJIVt7VrQlaXkICKSwrZ2LejKUnIQEUlhW7sWdGUpOYiIpLCtXQu6spQcRERS2NauBV1Zmq0kIpLiunYNPxmUpJaDiIjEUHIQEZEYSg4iIhJDyUFERGIoOYiISAwlBxERiaHkICIiMZQcREQkhpKDiIjEUHIQEUmiZCzkE4TKZ4iIJEmyFvIJQi0HEZEkSdZCPkEoOYiIhKy0rqNkLeQThLqVRERCVFbXUUaGt11S2Av5BJHUloOZTTKzn81sQdS+ZmY2x8y+9h93SGaMIiJbo6yuo2Qt5BNEsruVHgROLbFvEPCKc64d8Iq/LSJSLZXVdZSshXyCSGpycM69Cawpsfts4CH/+UPAOYmMSUSkKpW3BnTXrrBsGRQVeY+pkBgg+S2HeHZyzv0A4D/uGO8gM+tlZnPNbO6qVasSGqCISGlKDj6ffnoSu46cg9mz4ccfK3xqKiaHQJxzOc65js65ji1atEh2OCIimwefc3O97+XcXHjoIejePcFdR0VFMHMmdOzoZaecnApfIhVnK/1kZi2dcz+YWUvg52QHJCISRGmDz88/73UZha6gAKZPh2HD4MsvoV07mDSpUpkoFVsOzwDd/efdgaeTGIuISGBJu29h0yYYNw723NNrptSvDzNmwKJF0KOHt11ByZ7KOh14D9jLzFaY2eXAcOBkM/saONnfFhFJeeUNPle533+H0aNh992hTx/YaSd45hmYPx8uvBDS0ip96aR2KznnupTy0okJDUREpApkZ295wxuENPi8di3ccw+MGQN5eXDCCTB1qvdoViVvkYrdSiIi1VLo9y38/DMMHuw1Rf77XzjiCHjvPXjlFTjxxCpLDJCaA9IiItVW164hzET67jsYORIeeAD++AM6d4YhQ6BDhyp+o78oOYiIpKolS2D4cJgyxZsbe+mlMGiQN/AcMnUriYhshVAW61mwAC6+GPbaC6ZN8wYylizxpqUmIDGAWg4iIhUSiXj3MyxfDs2awfr1kJ/vvbbVi/V8+CEMHQpPPw2NG8OAAXDttbDzzlUWf1BqOYiIBFTyDui8vL8SQ7EKL9bjHLz+Opx8Mhx2GLz5Jtx8s/cmt9+elMQAajmIiAQW7w7oeALd9FZc9yg7G95917tHYcQI6N0bmjTZ6li3lpKDiEhAQe90LvOmt8JCeOIJr/to/nzv4HvvhZ49YZttqiLMKqFuJRGRgILc6VzqTW8FBV4Vvv32gwsu8Jogkyd7A81XXZVSiQECJAcza2tmDfznx5lZPzNrGnpkIiIpJt7KbfXqQXp6GTe9bdoEY8d6s4wuuwwaNIBHHoEvvvC269VL4CcILkjLYSZQaGZ7ABOBNsDDoUYlIpJCiqerXnqp9wM/OhlMngyrV8dZrOe332DUKGjTBjIzvYHlZ5+FTz7xWg5bUfcoEYKMORQ55/40s07AXc65e8zsk7ADExFJlrKmq+blea2HqVNLma76yy9/1T1as8Yra/Hww3DccVVa3iJsQVoOBWbWBa989ix/X2q2g0REtlKlp6v+9JN393KrVnDTTXDUUfD++/Dyy3D88dUqMUCwlkMPoDeQ7Zz71szaANPCDUtEJDkqPF11+XKv7tGECV4WueACrzje/vuHGmfYyk0OzrkvzGwgkOFvf4vWWBCRGirodNVjW34Nl/t1jwC6dfNaDu3ahRdcAgWZrfQvYD7wgr99gJk9E3JcIiJJUd501f34nEfTuvDqD3t7Ywm9e8PSpTBxYo1JDBBszOFm4FBgLYBzbj7ejCURkWovunBe8+bezKOS6tWDk7f7gKc5i8/Zn3PqzaLOjTd405PuuSfEpd6SJ8iYw5/OuV9ty8EUF1I8IiIJUzz4XDzGkJdX8gjH2du9ztjdsmm58BVv6lL//1Hv6qthhx0SHW5CBUkOC8zsYiDNzNoB/YB3ww1LRCR8pQ8+O87gOYYwlCPXvQd5O3v3LFx5pVcttRYI0q10NdAe+AOYDqwDrgkxJhGR0EUi3nTVaHUopDOP8gkHMot/8TdWksn98O23cP31tSYxQLDZShuALP9PRKTaK+5OKlaXAroSYRDD2Zuv+JK96M6DPMzF7NKqHjRMXqzJUmpyMLNnKWNswTl3VigRiYiEKBKB7t294qgN2UhPJnEjI2jFcj7hADrzKE9wLkWklV5ErxYoq+UwKmFRiIgkQHGLYZvC9fRmHNdzBzvzE+9wJH0Yy2xOIz3dcGugVYaXGCq1olsNUGpycM69kchARETCNmLQGgZsuIf+jKEZvzCHk7iIGbzBPwCjVStvdqqU3a30qHPuAjP7nDjdS8656n1vuIjUGjPv+5FVQ0bz9rqxNOE3nuJshjKEjzh08zG1uQspnrK6lfr7j2cmIhARkSqXm8vrZ47k9AUTqU8+j3AhwxjMAv6+xWFpaXHWYajlSp3K6pz7wX+a6ZzLjf4DMhMTnohIJSxezNJje1DQeg+OXJBDhK7sxVd05eGYxNCokbdAmxLDloLc53BynH2nVXUgIiKVFYl4pS862Kc8YhdStNfetHzrEe4nk7Ys5QomsJQ9Ys6Lu3KbAGWPOfTBayHsbmafRb3UBHgn7MBERMoTiUD//rBH3vtMJpt/MYt1NOF2BnIn17KKHUs9V4PPZStrzOFhYDYwDBgUtX+9c25NqFEBZrYMWA8U4tV36hj2e4pI9ZHZx7F43Ks8wlBO5FXyaMZ/uIV76ctayq57ZKbB5/KUNZX1V+BXoIuZpQE7+cc3NrPGzrmAVc+3yvHOuTg1EkWktopMczyfOYur12dzOB+wkpZcxx3k0IvfKb+8hZlXZVtdSWUrt3yGmfXFK9v9E1Dk73aAprKKSOIUFjLxtMfpOGcoXfmMb2lNb8byIJfxR8D6Funp3tLOSgzlC1KV9RpgL+dcTDHbkDngJTNzwHjnXE6C319EUkF+Pu/1ncaOE4dzedHXLGJvuvEQ0+nCnwGXs1dSqLggyeE7vO6lRDvKObfSzHYE5pjZl865N4tfNLNeQC+AjBq40IZIrbdxI0ycyJrBIzjit+/4mAM5j8d5kk64ciZa1qkDRUXeoHNtLoGxNYIkh2+A183sObyy3QA450aHFpV3/ZX+489m9iTeanRvRr2eA+QAdOzYUYsPiVRzxTOP8vPW0ZtxXMdoduYnvuAoshnPC5wKWKnnF48l3H9/4mKuyYIkh+X+X33/L3Rmti1Qxzm33n9+CnBLIt5bRBIrEvHW0Gnwex79uJt+3M0OrOVFTuECsniLY8u9hrqNql6Q9Rz+l4hAStgJeNJfmrQu8LBz7oUkxCEiIcrMhCfH/sBNjKYPY2nM7zzJOQxlCHM5pNzz1VoIT5DZSi2AG/FWg9s8JcA5d0JYQTnnvgE6hHV9EUm+/7s0l/bTRjCaidSjgBlcxDAGs5D9Ap2v1kK4gnQrRYBH8Arw9Qa6A6vCDEpEaq5nRn7FuiHDuOnPCA7jQS5jBDfGLW9Rmj591FoIW5DkkO6cm2hm/f01Ht4wM631ICIVM38+884fyplLH2cTDbmPqxjFAL5n1wpdRokhMYIkhwL/8QczOwNYCRX8tykitde773rzSZ9/nj3YjuEM4i6uKbPuUTzqRkqsIMnhNjPbHrgeuAfYDrg21KhEpFqLTHPMzHyFq9dnczyvs5p07uQ27uMqfqVpmec2bgzjxikJJFuQ2Uqz/Ke/AseHG46IVGtFRbw+YBbt7srmCfch3/M3rmU0OfRiA9uWe7q6jFJHkNlKk4m/TGjPUCISkeqnsJC3+z3KDuOHcVzh53xDG3oxnofoTj4NAl1CiSG1BOlWmhX1vCHQCW/cQURqu/x8pv1zKoe9PpyjWcIX7MMlTGUGF1EY6OtF3UipKki30szobTObDrwcWkQikvJmTN7I/L4TyNwwkkv4jnkcxLnM5CnOKbfuUTS1FlJXsNS+pXaAKt2J1Ebr1jHhoPv519I7uYifeYuj6UUOL/JPyqp7FI8SQ2oLMuawHm/MwfzHH4GBIcclIinksXF5LO03hisL7uHfrOUF/kk2WbzNMRW+lqakVg9BupWaJCIQEUktkQjcdMVKem8cTW/G0Zjfmcm5DGMw86jYqr1qJVQ/ZSYHM9sG6Ars6++aCzzunMsPOzARSbzistmN85YxkNtZyCTSKGQ6XRjGYBZt/ioIpk4dr+KqEkP1U+rIkZn9HVgEHAMsA3KBfwLvmFlTM7stIRGKSKgyM73qpmZw6yVfckded5awBz2ZxINcxp4sphtTAyeGhg1h2jRwDgoLlRiqq7JaDncDVzjn5kTvNLOTgAXAwjADE5HwFK+h8Pvv3vYBfMIQhnIeM9lEQ+7hakYxgJXsEviaGkuoWcpKDi1LJgYA59zLZlaAd7+DiFQjkQj07An5fsfwkbxDFtmczmx+ZTuGMoQx9Gc1LQJf88QT4WVNbq9xypqQXMfMYm5tNLOGQIFzbkN4YYlIVYtE4NJLIT/fcRJzeI3jeIejOYSPGEI2GSznP9wWODE0bux1Hykx1ExlJYcpwEwza128w3/+KDA13LBEpCpFItDtkiLOck/xAYcxh1PYgyVcw520ZhnDGMI6ti/3OsUJwTlYv15dSDVZqcnBOXcb8ALwppmtNrPVwBvAHOfcrYkKUEQqLxKBRvX/5LlLHmY+HXiKTqSTxxXk0JaljOGacgvimXlTUZUQapcyp7I65+4F7jWzJv72+oREJSJbJTMTJozNpxtT+Izh7MFSFrIvXZnGI1wYqO6RBphrt0DlM5QURKqPg/fZwNFfPsBSRrEbK/iIjpzDkzzDWWXWPVIBPIlWmdpKIpJiIhG4odevdN9wP7O5kx1ZxRscy+VMZA4nU1bdI802kniUHESquXOPXc2Bb43hC+6hKb8ym1PJJot3OLrcc1XWQkoTpPBeI7wlQjOcc1eYWTtgr6gV4kQkCQZ3X0mLKaOYyni2YSNPcC5DGcInHBTofCUGKUuQlsNkYB5whL+9AniMLRcBEpEEyMyE2WO/ZSC3czOTSaOQh7mY4QyqUHmLCRM0tiBlC7IqR1vn3AigAMA5t5GKFm4XkUqLRLzB4n3tCw4f242vaUcPJjOZHuzJYrozJVBiSE/37lHYuFGJQcoXpOWQ71dndQBm1hb4I9SoRGRzqYv2+R/zENl04kk2sg1j6M9orgtc92jffWGhKqFJBQVJDjfh3Qy3m5lFgKOAy8IMSqS2i0Rg3CVv8xTZnMYLrGV7ssliDP3Jo3mga5hB794aV5DKCbLYzxwz+xg4HK87qb9zbnXokYnUQpl9HEvHvUQW2bzFW/xMCwYzlPvJDFTeAtRSkKpRanIws5JTHn7wHzPMLMM593F4YYnUHpEI9O5VxEkbnmYIQzmEuXzHrvRjDBP4NxtpFPhamoEkVaWslsMdZbzmgBOqOBaRWqdv7z/5ZfwjvMcw9mMhS2jLv3mAKXSjgPqBrqHuIwlDqcnBOXd8IgMpycxOBcYAacAE59zwZMYjUlUyM2Hi2D/ozkMM5Hba8g0LaM/FRHiUCwLVPQK1EiRc5U5lNbOGZnadmT1hZjPN7Bp/TYfQmFkacB9wGt761V3MrGKL14qkmEgEtkv7nfpj72IpbcnhSvJI52yeYn8+YzoXl5sYopfgVGKQMAX5iTIFWA/c4293wVvPoXNYQQGHAkucc98AmNkM4GzgixDfU6RKRS/FuR2/chX3sZQ7acFqXucf9GAyL3MSQW8bUg0kSaQgyWEv51yHqO3XzOzTsALy7QJ8F7W9Ajgs+gAz6wX0AsjIyAg5HJHgopfibM4qBnMXfbmX7VnHc5zOUIbwLkcFvp7GFCQZgtwh/YmZHV68YWaHAe+EF5L3NnH2uS02nMtxznV0znVs0SL4erciVa34DmYz7++SS6B5/veM5lqW0ZrBDOMlTuFAPuZMngucGIoX2CkqUmKQxAvScjgM6GZmy/3tDGCRmX0OOOfc/iHEtQLYLWp7V2BlCO8jUmnRLYRiu7OUGxnBZTxIGoVM4xJuZyBfsk/g62qgWVJBkORwauhRxPoIaGdmbYDvgYuAi5MQh0hcJ50Er7zy1/a+LGQww+jCdAqox0QuZwQ3kkvrwNfUYjuSSoLcIZ1rZjvg/ZKvG7U/tJvgnHN/mllf4EW8qayTnHO651OSLhKB7t2hsNDbPoh5ZJHNuTzJb2zLnVzLHVzPj7QMdD0NMkuqCrKew614tZSW8le/f+g3wTnnngeeD/M9RILIzISxY7fcdwxvMoShnMqL/EJTbuE/jKE/a0gPdE21EiTVBelWugCvbHd+uUeK1CAlWwng+CcvkkU2x/A2P7EjAxnOWPqwnu3KvJbWUJDqJkhyWAA0BX4ONxSR1BE9pmAUcQ5PMYShdGQey9mNq7mbiVxebt0jtRCkugqSHIbhTWddQNQ6Ds65s0KLSiTBom9YK5bGn3RhOoMZxr4s4mv2oCcTmcYlpdY9UgtBaoogyeEh4Hbgc6Ao3HBEEq/kzKP6/MFlPMhAbmd3vuUz/s5FTOcxOlNEWtxr1K0LDz6opCA1R5DksNo5d3fokYgkWMkxhUb8Ti9yGMAodmElH3Ao13AXszgTV8b9oppxJDVRkOQwz8yGAc+wZbeS1nOQaicz0xsDcFH322/PWvpyL9dwF83J41WOpxtTeJUTKKvukcYTpCYLkhwO9B8Pj9qn9RykWomdeeTVPbqWO7mK+9iedcziDLLJ4n2OKPU6aiVIbRHkJrikrusgsjXi3aOwCysYwCh6kUNDNvE45zOUIXzKAaVeR2MKUtsEWlXEzM4A2gOb13Fwzt0SVlAiWyteUmjLEgZyO915CMMxjUsYziAWs1ep11FFVKmtgtwhPQ5oBBwPTADOBz4MOS6RCos3HRWgPQsYzDAuYgYF1OMBrmAkN5Ra90itBJFgLYcjnXP7m9lnzrn/mdkdwBNhByYSVLxWAkBHPmIIQ+nEU6ynMXdwPXdybZl1jzSmIOIJsp7DRv9xg5n9DSgA2oQXkkgwkQg0aFAyMTiO5Q1e5BQ+4lD+wRvczE20IpeBjCg1MTRu7C2/qcQg4gnScphlZk2BkcDHeDOVHggzKJHylLxxDRyn8gJZZHM07/AjO3EjtzOWPvxGk7jXUCtBpHTlthycc7c659Y652YCrYC9nXP/DT80kViZmd4gcXTdo/N4nHkczGxOJ4Pl9OUe2vAtI7kxbmJQK0GkfKUmBzM7xMx2jtruBjwK3GpmzRIRnAhsuQxncRdSXQq4lCkspD2P05nG/EYPJrEHS7iPvmximy2uUbeulxCcg/XrNdgsUp6yWg7jgXwAMzsWGA5MAX4FcsIPTWqreGsyF89AasAmrmQci9mTKXQnn/pcyAz2YREP0iOmIJ6Zt+xmQYESgkhFlDXmkOacW+M/vxDI8buWZprZ/NAjk1qptJlH2/IbVzKe67mDv/ED73MY/bibWZxJvBIXWodZZOuUmRzMrK5z7k/gRKBXwPNEAivt3oRiTfmFvtxLf8bQnDxe4QQuYRqvcTxKCiLhKetLfjrwhpmtxpvO+haAme2B17UkUmnlJYUW/Ly57tF2rOdZziSbLD7YosTXX7SOgkjVKjU5OOeyzewVoCXwknOb61jWAa5ORHBSM5XWdQSwK99xAyO5ggdowB88ygUMYzCf0SHu8UoKIuEos3vIOfd+nH2LwwtHaqp4pbKjtWUJgxhON6ZgOKZyKcMZxNfsGfd4lcsWCZfGDiR0ZbUU9uNzBjOMC3mEAuqRQy9GcgPLaRVzbJ06XleUxhREwqfkIKGKROInhkP4kCyyOZtnWE9jRjGAO7mWn9h5i+PS0uChh9RCEEm0ILWVRALLzPR+4Uffo/AXxz94nZc4mQ85jGN4i5u4mVbkMojbYxJDeroSg0iyqOUgVSISgZ49IT8/3quO03meIQzlKN7lB3ZmACMZz5VblLfQNFSR1KHkIJUSiUD//pCXV/oxdSjkXJ5gCEM5kPnkkkEm9zGJnvzx17pRGlwWSUHqVpIKy8z0uotKSwx1KaAbD7GQ9jzGBTRiA5cxmT1YwlgyNyeG9HSv3pFqHYmkHrUcJJAgLYUGbKInk7iREbQml/l0oDOP8gTnUkTa5uPq14dJk5QQRFKZkoOUKxKBHj284nXxbMtv9GYc13MHLfmRdzmCq7iP5zmdkiUu1IUkUj2kXLeSmd1sZt+b2Xz/7/Rkx1QbRSLQvPlfM47iJYYdWMN/uIVcWjGKG1hIe47nVY7iHZ7nDKITg7qQRKqXVG053OmcG5XsIGqbIF1HADvyE9cxmkzupwm/8TRnMZQhfMhhm49RC0GkekvV5CAJVtZdzMV2Yzk3MJJ/M4H65G+ue/Q5+28+Jj0dxoxRUhCp7lKuW8nX18w+M7NJZrZDsoOpyYq7j8pKDO1YzER6spS29GYcD3Mxe/MlFzOdlen7b15hzTlYvVqJQaQmSErLwcxehhK3w3qygLHArYDzH+8Aesa5Ri/8NSYyMjJCi7UmK68Y3t/5jCEMpTOPkU99xtKHUQzgOzKoUwemTVEiEKmpzJX2zZACzKw1MMs5t19Zx3Xs2NHNnTs3MUHVEJEIXHpp/MRwKB+QRTZn8SzraML9ZHIn1/IzOwGaiipSU5jZPOdcx3ivpVy3kpm1jNrsBCxIViw1VSQC3buXTAyO43mVOZzEBxzOUbzDf7iFVuQymOGbE0N6uhKDSG2QigPSI8zsALxupWXAlUmNpoaJ7UpynMFzZJHNEbzPD+zM9YxiPFdijRtrxpFILZVyycE5d2myY6gpIhHIyoLcXO9+heiWQh0KOY+ZDGEoB/Apy2hFH+5nMj1onN6Q8ZpxJFKrpVy3klRcJAKtW3sJoG5d77FOHe/mtdxc75jixFCXAi5jMl+wL49yIQ3ZRHcepB1fM9760LNPQ804EpHUazlIcPFuWiss9B5LDjQ3ZOPmuketWM4nHMD5PMaTdKKINNLSYKrWThARn5JDNVNWV1E8jVm/ue7RzvzEOxxJH8Yym9MoLm9hpkV1RGRLSg7VRLxWQlmJYQfW0I+76cfdNOMXXuJkLiSLNzmW6JpHZtC7txKDiGxJyaEaiESgVy/YsKH8Y3fiR65jNH0YSxN+4ynOZihD+IhDY45VqQsRKY2SQ4op7jZavhyaNfP2lVcIDyCD3M11j+pRwCNcyDAGs4C/xxyrpCAi5VFySAGljSMESQp78hWDGM4lTMNhTKE7wxnIsrQ9KCyEVq0gO1uJQEQqRskhyUp2GQWtZrI/n26ue7SJhkxrkknT2wbw73678e/wwhWRWkLJIQmiu47q1Plr+mkQh/MeWWRzJs+x3pqw6MyBtJ9wLT123DG8gEWk1lFySLCSLYVgicFxAq+SRTYn8Bpr6qTz6bm30uGBvrRv2jTEaEWkttId0gmWlRVs1pHHcSbP8h5H8AonsY99ybyL76DZr8vo8Nj/gRKDiIREySHBli8v/5g0CrmQGXyedgDPchY78jNDmo3j9YnfcHDkOm8NThGRECk5hKy47lGdOt5j8fTUktLSoD75DEifxJqd92EGXdhvzwKYMoXdCxYzNO9KuvRomMjQRaQW05hDiEqOL+TmQr163mI5+fl/Hddsm4280Hkih7w2Ar77Dg48EO59HDp18rKKiEiC6ZsnRPHGFwoKoEkT7/6D7VjHsKa3s6Jeaw6ZcjVkZMDzz8O8eXDeeUoMIpI0ajmEqNTxhbw8ll11N9x9N6xdC6ec4mWSY49NZHgiIqXST9MqVN74ws78wAhuINdawS23wPHHw4cfwosvKjGISEpRy6GKlDW+0DJ/GTcygp5Moh4FLD+iC21yBkP79skNWkSkFEoOVSTe+MLuBV9yU4PhnE8Eh/F448vY9n8DOfu6tskJUkQkICWHKhI9vtCB+WSRzXnMZNMfDanX/yoYMICLd901eQGKiFSAxhwqoeTYQiTiTTQ6gneZxRnM50BO4SWGMZijd82Fu+4CJQYRqUbUcqig2LEFx/TLX+Gl9Gz25HVWk04Wt3EfV1HQqCk5w5Mbr4hIZajlUEHFYwtGEWfxNO9zOLP+OJntflrMvK6jOWa3XIZZFk1bNSUnR+soiEj1pORQQrwuo2grcgu5iOl8Sgee5hxasIpejKd14TccPO1aFi3flqIiWLZMiUFEqi91K0WJNx21Vy/vedfO+TB1Kl/XHU6bP5ewkH25hKnM4CIKqUurVsmLW0Skqik5RIk3HdVt2MDXV0+AQSNhxQq2a3MwF33/BI/mn43zG16NGnlLcYqI1BQ1tlupvO6heKKnozZhHQMZzjJac/Mv/aFNG3jhBdKXfsS/JnUio1UdzLwaSRpbEJGaxlzQRYtTWMeOHd3cuXM3b5fsHgLv1315X+KtW8Nvuavpx91czT3swFpmcyqTdhrCYz8eE94HEBFJAjOb55zrGO+1GtlyiNc9tGGDt79UK1cye9/rWUZr/sutvMoJHMxczm80m3PuUGIQkdqlRiaH0qqhxt3/7bfQpw+0acM+L41h1VGdOKnlQjrbTPJaHawuIxGplZKSHMyss5ktNLMiM+tY4rXBZrbEzL4ys39W5voZGQH2L1oE3btDu3YwaRJcdhksXkybt6fy8sp9NR1VRGq1ZLUcFgDnAm9G7zSzfYGLgPbAqcD9ZpZW0YtnZ3tjDNE2zyj65BM4/3yvIurjj0O/fvDNNzB+POy+eyU/johIzZKUqazOuUUAZlbypbOBGc65P4BvzWwJcCjwXkWuX/xrPyvL60rKyICc7u9wSiQbZs+G7beHIUOgf39o0WJrP46ISI2TamMOuwDfRW2v8PfFMLNeZjbXzOauWrUq5vWuXWHZt46iF15iWat/cMotR8NHH8HQod7dbbfdpsQgIlKK0FoOZvYysHOcl7Kcc0+XdlqcfXHn2jrncoAc8KayxhyQlwenneYlhF128SqjXnFFbH+TiIjECC05OOdOqsRpK4DdorZ3BVZWKoBmzbwb1664Arp1gwYNKnUZEZHaKNXKZzwDPGxmo4G/Ae2ADyt1JTN45JEqDE1EpPZI1lTWTma2AjgCeM7MXgRwzi0EHgW+AF4ArnLOFSYjRhGR2ixZs5WeBJ4s5bVsQGXsRESSKNVmK4mISApQchARkRhKDiIiEkPJQUREYig5iIhIDCUHERGJUSNWgjOzVUBusuPYSs2B1ckOIkFq02eF2vV59Vmrl1bOubhF5mpEcqgJzGxuacv11TS16bNC7fq8+qw1h7qVREQkhpKDiIjEUHJIHTnJDiCBatNnhdr1efVZawiNOYiISAy1HEREJIaSg4iIxFBySCFmNtLMvjSzz8zsSTNrmuyYwmJmnc1soZkVmVmNnA5oZqea2VdmtsTMBiU7njCZ2SQz+9nMFiQ7lrCZ2W5m9pqZLfL/G+6f7JjCoOSQWuYA+znn9gcWA4OTHE+YFgDnAm8mO5AwmFkacB9wGrAv0MXM9k1uVKF6EDg12UEkyJ/A9c65fYDDgatq4r9bJYcU4px7yTn3p7/5Pt4a2jWSc26Rc+6rZMcRokOBJc65b5xz+cAM4OwkxxQa59ybwJpkx5EIzrkfnHMf+8/XA4uAXZIbVdVTckhdPYHZyQ5CKm0X4Luo7RXUwC+Q2s7MWgMHAh8kOZQql5RlQmszM3sZ2DnOS1nOuaf9Y7Lwmq6RRMZW1YJ81hrM4uzTvPEaxMwaAzOBa5xz65IdT1VTckgw59xJZb1uZt2BM4ETXTW/CaW8z1rDrQB2i9reFViZpFikiplZPbzEEHHOPZHseMKgbqUUYmanAgOBs5xzG5Idj2yVj4B2ZtbGzOoDFwHPJDkmqQJmZsBEYJFzbnSy4wmLkkNquRdoAswxs/lmNi7ZAYXFzDqZ2QrgCOA5M3sx2TFVJX9iQV/gRbwBy0edcwuTG1V4zGw68B6wl5mtMLPLkx1TiI4CLgVO8P8/nW9mpyc7qKqm8hkiIhJDLQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOkjLMLD1qauCPZva9/3ytmX2R4FjOiS6mZma3mFmFb+ozs9alVSo1s/Zm9qqZLTazpWb2PzOr8v8ny/osZvZ6Ta2KK1tHyUFShnMuzzl3gHPuAGAccKf//ACgqKrfz8zKqhBwDl411eLY/uuce7kK33sbvJvihjvn9gT+jlesL4zyz+cQ4meRmknJQaqLNDN7wK+f/5L/5YqZtTWzF8xsnpm9ZWZ7+/tbmdkr/toYr5hZhr//QTMbbWavAbfHO9/MjgTOAkb6LZe2/nnn+9c4xMzeNbNPzexDM2vitxDeMrOP/b8jy/k8FwPvOOdeAvDviO8L3OC/x81mNqD4YDNb4Bd5w8ye8uNdaGa9oo75zcyy/bjeN7Odyvss0czsFDN7z4//Mb92EGY23My+8P9Zjqr4vzqpjpQcpLpoB9znnGsPrAXO8/fnAFc75w4GBgD3+/vvBab4a2NEgLujrrUncJJz7vp45zvn3sX7VX+D35JZWnyiXwrjEaC/c64DcBKwEfgZONk5dxBwYYn3i6c9MC96h/8+21j5izz19OPtCPQzs3R//7bA+35cbwJXlPVZoplZc+D//H8uBwFzgevMrBnQCWjv/7O8rZzYpIZQ4T2pLr51zs33n88DWvu/bI8EHvPK3QDQwH88Am8xIYCpwIioaz3mnCss5/zS7AX84Jz7CKC4GqeZbQvca2YHAIV4CagsRvwqrfGquZbUz8w6+c93w0uceUA+MMvfPw84OcC1ih2O1/X0jv/Poj5eOYx1wCZggpk9F3V9qeGUHKS6+CPqeSGwDV7Ld60/LlGe6C/i3/3HipxfrLQv9WuBn4AO/nU3lXOdhcCxW1zYbHdgtXNurZn9yZYt+4b+McfhtVaOcM5tMLPXi18DCqIq+RZSsf+/DZjjnOsS84LZocCJeMUD+wInVOC6Uk2pW0mqLf9X+7dm1hm8aplm1sF/+V28LzOArsDbFTx/PV4RxJK+BP5mZof45zTxB7a3x2tRFOEVZUsrJ/wIcHTUrKFt8LqibvJfXwYc5L92ENDG37898IufGPbG+8VfntI+S7T3gaPMbA//PRuZ2Z5+62p759zzwDV4kwOkFlBykOquK3C5mX2K92u8eCnOfkAPM/sM78u6tFlApZ0/A7jBzD4xs7bFB/tLfl4I3OOfMwfvl/v9QHczex+vS+l3yuCc24g3UJxlZouB1XgD1MULPM0EmpnZfKAP3priAC8Adf3PdSvel3p54n6WEvGsAi4DpvvXfh/YGy+pzPL3vYHXQpJaQFVZRVKAmZ0DjAaOd87lJjkcESUHERGJpW4lERGJoeQgIiIxlBxERCSGkoOIiMRQchARkRhKDiIiEuP/Ac3LuQUv5hqGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAApbUlEQVR4nO3dd3gU5drH8e9NaCIoElA5agARGx6xYC/HfmxHRUVFFASPCBHBglLynqNHDSAgig2IFAVWsGBFUbH3AooKoghKELFAEEEBE5Pn/WMmuGQ3ySRkdjfJ73NduXZndmb2Xsve+5S5H3POISIiEq1OsgMQEZHUo+QgIiIxlBxERCSGkoOIiMRQchARkRh1kx1AVWjevLlr3bp1ssMQEalW5s2bt9o51yLeazUiObRu3Zq5c+cmOwwRkWrFzHJLe03dSiIiEkPJQUREYig5iIhIDCUHERGJoeQgIiIxlBxERGqoSARat4Y6dbzHSCT4uTViKquIiGwpEoFevWDDBm87N9fbBujatfzz1XIQEamBsrL+SgzFNmzw9geR1ORgZpPM7GczWxC172Yz+97M5vt/pyczRhGR6mj58ortLynZLYcHgVPj7L/TOXeA//d8gmMSEan2MjIqtr+kpCYH59ybwJpkxiAiUhNlZ0OjRlvua9TI2x9EslsOpelrZp/53U47xDvAzHqZ2Vwzm7tq1apExyciktK6doWcHGjVCsy8x5ycYIPRAJbsNaTNrDUwyzm3n7+9E7AacMCtQEvnXM+yrtGxY0enwnsiIhVjZvOccx3jvZZyLQfn3E/OuULnXBHwAHBosmMSEaltUi45mFnLqM1OwILSjhURkXAk9SY4M5sOHAc0N7MVwE3AcWZ2AF630jLgymTFJyJSWyU1OTjnusTZPTHhgYiIyBZSrltJRESST8lBRERiKDmIiEgMJQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOIiIpbmvWgq4srSEtIpLCtnYt6MpSy0FEJIVt7VrQlaXkICKSwrZ2LejKUnIQEUlhW7sWdGUpOYiIpLCtXQu6spQcRERS2NauBV1Zmq0kIpLiunYNPxmUpJaDiIjEUHIQEZEYSg4iIhJDyUFERGIoOYiISAwlBxERiaHkICIiMZQcREQkhpKDiIjEUHIQEUmiZCzkE4TKZ4iIJEmyFvIJQi0HEZEkSdZCPkEoOYiIhKy0rqNkLeQThLqVRERCVFbXUUaGt11S2Av5BJHUloOZTTKzn81sQdS+ZmY2x8y+9h93SGaMIiJbo6yuo2Qt5BNEsruVHgROLbFvEPCKc64d8Iq/LSJSLZXVdZSshXyCSGpycM69Cawpsfts4CH/+UPAOYmMSUSkKpW3BnTXrrBsGRQVeY+pkBgg+S2HeHZyzv0A4D/uGO8gM+tlZnPNbO6qVasSGqCISGlKDj6ffnoSu46cg9mz4ccfK3xqKiaHQJxzOc65js65ji1atEh2OCIimwefc3O97+XcXHjoIejePcFdR0VFMHMmdOzoZaecnApfIhVnK/1kZi2dcz+YWUvg52QHJCISRGmDz88/73UZha6gAKZPh2HD4MsvoV07mDSpUpkoFVsOzwDd/efdgaeTGIuISGBJu29h0yYYNw723NNrptSvDzNmwKJF0KOHt11ByZ7KOh14D9jLzFaY2eXAcOBkM/saONnfFhFJeeUNPle533+H0aNh992hTx/YaSd45hmYPx8uvBDS0ip96aR2KznnupTy0okJDUREpApkZ295wxuENPi8di3ccw+MGQN5eXDCCTB1qvdoViVvkYrdSiIi1VLo9y38/DMMHuw1Rf77XzjiCHjvPXjlFTjxxCpLDJCaA9IiItVW164hzET67jsYORIeeAD++AM6d4YhQ6BDhyp+o78oOYiIpKolS2D4cJgyxZsbe+mlMGiQN/AcMnUriYhshVAW61mwAC6+GPbaC6ZN8wYylizxpqUmIDGAWg4iIhUSiXj3MyxfDs2awfr1kJ/vvbbVi/V8+CEMHQpPPw2NG8OAAXDttbDzzlUWf1BqOYiIBFTyDui8vL8SQ7EKL9bjHLz+Opx8Mhx2GLz5Jtx8s/cmt9+elMQAajmIiAQW7w7oeALd9FZc9yg7G95917tHYcQI6N0bmjTZ6li3lpKDiEhAQe90LvOmt8JCeOIJr/to/nzv4HvvhZ49YZttqiLMKqFuJRGRgILc6VzqTW8FBV4Vvv32gwsu8Jogkyd7A81XXZVSiQECJAcza2tmDfznx5lZPzNrGnpkIiIpJt7KbfXqQXp6GTe9bdoEY8d6s4wuuwwaNIBHHoEvvvC269VL4CcILkjLYSZQaGZ7ABOBNsDDoUYlIpJCiqerXnqp9wM/OhlMngyrV8dZrOe332DUKGjTBjIzvYHlZ5+FTz7xWg5bUfcoEYKMORQ55/40s07AXc65e8zsk7ADExFJlrKmq+blea2HqVNLma76yy9/1T1as8Yra/Hww3DccVVa3iJsQVoOBWbWBa989ix/X2q2g0REtlKlp6v+9JN393KrVnDTTXDUUfD++/Dyy3D88dUqMUCwlkMPoDeQ7Zz71szaANPCDUtEJDkqPF11+XKv7tGECV4WueACrzje/vuHGmfYyk0OzrkvzGwgkOFvf4vWWBCRGirodNVjW34Nl/t1jwC6dfNaDu3ahRdcAgWZrfQvYD7wgr99gJk9E3JcIiJJUd501f34nEfTuvDqD3t7Ywm9e8PSpTBxYo1JDBBszOFm4FBgLYBzbj7ejCURkWovunBe8+bezKOS6tWDk7f7gKc5i8/Zn3PqzaLOjTd405PuuSfEpd6SJ8iYw5/OuV9ty8EUF1I8IiIJUzz4XDzGkJdX8gjH2du9ztjdsmm58BVv6lL//1Hv6qthhx0SHW5CBUkOC8zsYiDNzNoB/YB3ww1LRCR8pQ8+O87gOYYwlCPXvQd5O3v3LFx5pVcttRYI0q10NdAe+AOYDqwDrgkxJhGR0EUi3nTVaHUopDOP8gkHMot/8TdWksn98O23cP31tSYxQLDZShuALP9PRKTaK+5OKlaXAroSYRDD2Zuv+JK96M6DPMzF7NKqHjRMXqzJUmpyMLNnKWNswTl3VigRiYiEKBKB7t294qgN2UhPJnEjI2jFcj7hADrzKE9wLkWklV5ErxYoq+UwKmFRiIgkQHGLYZvC9fRmHNdzBzvzE+9wJH0Yy2xOIz3dcGugVYaXGCq1olsNUGpycM69kchARETCNmLQGgZsuIf+jKEZvzCHk7iIGbzBPwCjVStvdqqU3a30qHPuAjP7nDjdS8656n1vuIjUGjPv+5FVQ0bz9rqxNOE3nuJshjKEjzh08zG1uQspnrK6lfr7j2cmIhARkSqXm8vrZ47k9AUTqU8+j3AhwxjMAv6+xWFpaXHWYajlSp3K6pz7wX+a6ZzLjf4DMhMTnohIJSxezNJje1DQeg+OXJBDhK7sxVd05eGYxNCokbdAmxLDloLc53BynH2nVXUgIiKVFYl4pS862Kc8YhdStNfetHzrEe4nk7Ys5QomsJQ9Ys6Lu3KbAGWPOfTBayHsbmafRb3UBHgn7MBERMoTiUD//rBH3vtMJpt/MYt1NOF2BnIn17KKHUs9V4PPZStrzOFhYDYwDBgUtX+9c25NqFEBZrYMWA8U4tV36hj2e4pI9ZHZx7F43Ks8wlBO5FXyaMZ/uIV76ctayq57ZKbB5/KUNZX1V+BXoIuZpQE7+cc3NrPGzrmAVc+3yvHOuTg1EkWktopMczyfOYur12dzOB+wkpZcxx3k0IvfKb+8hZlXZVtdSWUrt3yGmfXFK9v9E1Dk73aAprKKSOIUFjLxtMfpOGcoXfmMb2lNb8byIJfxR8D6Funp3tLOSgzlC1KV9RpgL+dcTDHbkDngJTNzwHjnXE6C319EUkF+Pu/1ncaOE4dzedHXLGJvuvEQ0+nCnwGXs1dSqLggyeE7vO6lRDvKObfSzHYE5pjZl865N4tfNLNeQC+AjBq40IZIrbdxI0ycyJrBIzjit+/4mAM5j8d5kk64ciZa1qkDRUXeoHNtLoGxNYIkh2+A183sObyy3QA450aHFpV3/ZX+489m9iTeanRvRr2eA+QAdOzYUYsPiVRzxTOP8vPW0ZtxXMdoduYnvuAoshnPC5wKWKnnF48l3H9/4mKuyYIkh+X+X33/L3Rmti1Qxzm33n9+CnBLIt5bRBIrEvHW0Gnwex79uJt+3M0OrOVFTuECsniLY8u9hrqNql6Q9Rz+l4hAStgJeNJfmrQu8LBz7oUkxCEiIcrMhCfH/sBNjKYPY2nM7zzJOQxlCHM5pNzz1VoIT5DZSi2AG/FWg9s8JcA5d0JYQTnnvgE6hHV9EUm+/7s0l/bTRjCaidSjgBlcxDAGs5D9Ap2v1kK4gnQrRYBH8Arw9Qa6A6vCDEpEaq5nRn7FuiHDuOnPCA7jQS5jBDfGLW9Rmj591FoIW5DkkO6cm2hm/f01Ht4wM631ICIVM38+884fyplLH2cTDbmPqxjFAL5n1wpdRokhMYIkhwL/8QczOwNYCRX8tykitde773rzSZ9/nj3YjuEM4i6uKbPuUTzqRkqsIMnhNjPbHrgeuAfYDrg21KhEpFqLTHPMzHyFq9dnczyvs5p07uQ27uMqfqVpmec2bgzjxikJJFuQ2Uqz/Ke/AseHG46IVGtFRbw+YBbt7srmCfch3/M3rmU0OfRiA9uWe7q6jFJHkNlKk4m/TGjPUCISkeqnsJC3+z3KDuOHcVzh53xDG3oxnofoTj4NAl1CiSG1BOlWmhX1vCHQCW/cQURqu/x8pv1zKoe9PpyjWcIX7MMlTGUGF1EY6OtF3UipKki30szobTObDrwcWkQikvJmTN7I/L4TyNwwkkv4jnkcxLnM5CnOKbfuUTS1FlJXsNS+pXaAKt2J1Ebr1jHhoPv519I7uYifeYuj6UUOL/JPyqp7FI8SQ2oLMuawHm/MwfzHH4GBIcclIinksXF5LO03hisL7uHfrOUF/kk2WbzNMRW+lqakVg9BupWaJCIQEUktkQjcdMVKem8cTW/G0Zjfmcm5DGMw86jYqr1qJVQ/ZSYHM9sG6Ars6++aCzzunMsPOzARSbzistmN85YxkNtZyCTSKGQ6XRjGYBZt/ioIpk4dr+KqEkP1U+rIkZn9HVgEHAMsA3KBfwLvmFlTM7stIRGKSKgyM73qpmZw6yVfckded5awBz2ZxINcxp4sphtTAyeGhg1h2jRwDgoLlRiqq7JaDncDVzjn5kTvNLOTgAXAwjADE5HwFK+h8Pvv3vYBfMIQhnIeM9lEQ+7hakYxgJXsEviaGkuoWcpKDi1LJgYA59zLZlaAd7+DiFQjkQj07An5fsfwkbxDFtmczmx+ZTuGMoQx9Gc1LQJf88QT4WVNbq9xypqQXMfMYm5tNLOGQIFzbkN4YYlIVYtE4NJLIT/fcRJzeI3jeIejOYSPGEI2GSznP9wWODE0bux1Hykx1ExlJYcpwEwza128w3/+KDA13LBEpCpFItDtkiLOck/xAYcxh1PYgyVcw520ZhnDGMI6ti/3OsUJwTlYv15dSDVZqcnBOXcb8ALwppmtNrPVwBvAHOfcrYkKUEQqLxKBRvX/5LlLHmY+HXiKTqSTxxXk0JaljOGacgvimXlTUZUQapcyp7I65+4F7jWzJv72+oREJSJbJTMTJozNpxtT+Izh7MFSFrIvXZnGI1wYqO6RBphrt0DlM5QURKqPg/fZwNFfPsBSRrEbK/iIjpzDkzzDWWXWPVIBPIlWmdpKIpJiIhG4odevdN9wP7O5kx1ZxRscy+VMZA4nU1bdI802kniUHESquXOPXc2Bb43hC+6hKb8ym1PJJot3OLrcc1XWQkoTpPBeI7wlQjOcc1eYWTtgr6gV4kQkCQZ3X0mLKaOYyni2YSNPcC5DGcInHBTofCUGKUuQlsNkYB5whL+9AniMLRcBEpEEyMyE2WO/ZSC3czOTSaOQh7mY4QyqUHmLCRM0tiBlC7IqR1vn3AigAMA5t5GKFm4XkUqLRLzB4n3tCw4f242vaUcPJjOZHuzJYrozJVBiSE/37lHYuFGJQcoXpOWQ71dndQBm1hb4I9SoRGRzqYv2+R/zENl04kk2sg1j6M9orgtc92jffWGhKqFJBQVJDjfh3Qy3m5lFgKOAy8IMSqS2i0Rg3CVv8xTZnMYLrGV7ssliDP3Jo3mga5hB794aV5DKCbLYzxwz+xg4HK87qb9zbnXokYnUQpl9HEvHvUQW2bzFW/xMCwYzlPvJDFTeAtRSkKpRanIws5JTHn7wHzPMLMM593F4YYnUHpEI9O5VxEkbnmYIQzmEuXzHrvRjDBP4NxtpFPhamoEkVaWslsMdZbzmgBOqOBaRWqdv7z/5ZfwjvMcw9mMhS2jLv3mAKXSjgPqBrqHuIwlDqcnBOXd8IgMpycxOBcYAacAE59zwZMYjUlUyM2Hi2D/ozkMM5Hba8g0LaM/FRHiUCwLVPQK1EiRc5U5lNbOGZnadmT1hZjPN7Bp/TYfQmFkacB9wGt761V3MrGKL14qkmEgEtkv7nfpj72IpbcnhSvJI52yeYn8+YzoXl5sYopfgVGKQMAX5iTIFWA/c4293wVvPoXNYQQGHAkucc98AmNkM4GzgixDfU6RKRS/FuR2/chX3sZQ7acFqXucf9GAyL3MSQW8bUg0kSaQgyWEv51yHqO3XzOzTsALy7QJ8F7W9Ajgs+gAz6wX0AsjIyAg5HJHgopfibM4qBnMXfbmX7VnHc5zOUIbwLkcFvp7GFCQZgtwh/YmZHV68YWaHAe+EF5L3NnH2uS02nMtxznV0znVs0SL4erciVa34DmYz7++SS6B5/veM5lqW0ZrBDOMlTuFAPuZMngucGIoX2CkqUmKQxAvScjgM6GZmy/3tDGCRmX0OOOfc/iHEtQLYLWp7V2BlCO8jUmnRLYRiu7OUGxnBZTxIGoVM4xJuZyBfsk/g62qgWVJBkORwauhRxPoIaGdmbYDvgYuAi5MQh0hcJ50Er7zy1/a+LGQww+jCdAqox0QuZwQ3kkvrwNfUYjuSSoLcIZ1rZjvg/ZKvG7U/tJvgnHN/mllf4EW8qayTnHO651OSLhKB7t2hsNDbPoh5ZJHNuTzJb2zLnVzLHVzPj7QMdD0NMkuqCrKew614tZSW8le/f+g3wTnnngeeD/M9RILIzISxY7fcdwxvMoShnMqL/EJTbuE/jKE/a0gPdE21EiTVBelWugCvbHd+uUeK1CAlWwng+CcvkkU2x/A2P7EjAxnOWPqwnu3KvJbWUJDqJkhyWAA0BX4ONxSR1BE9pmAUcQ5PMYShdGQey9mNq7mbiVxebt0jtRCkugqSHIbhTWddQNQ6Ds65s0KLSiTBom9YK5bGn3RhOoMZxr4s4mv2oCcTmcYlpdY9UgtBaoogyeEh4Hbgc6Ao3HBEEq/kzKP6/MFlPMhAbmd3vuUz/s5FTOcxOlNEWtxr1K0LDz6opCA1R5DksNo5d3fokYgkWMkxhUb8Ti9yGMAodmElH3Ao13AXszgTV8b9oppxJDVRkOQwz8yGAc+wZbeS1nOQaicz0xsDcFH322/PWvpyL9dwF83J41WOpxtTeJUTKKvukcYTpCYLkhwO9B8Pj9qn9RykWomdeeTVPbqWO7mK+9iedcziDLLJ4n2OKPU6aiVIbRHkJrikrusgsjXi3aOwCysYwCh6kUNDNvE45zOUIXzKAaVeR2MKUtsEWlXEzM4A2gOb13Fwzt0SVlAiWyteUmjLEgZyO915CMMxjUsYziAWs1ep11FFVKmtgtwhPQ5oBBwPTADOBz4MOS6RCos3HRWgPQsYzDAuYgYF1OMBrmAkN5Ra90itBJFgLYcjnXP7m9lnzrn/mdkdwBNhByYSVLxWAkBHPmIIQ+nEU6ynMXdwPXdybZl1jzSmIOIJsp7DRv9xg5n9DSgA2oQXkkgwkQg0aFAyMTiO5Q1e5BQ+4lD+wRvczE20IpeBjCg1MTRu7C2/qcQg4gnScphlZk2BkcDHeDOVHggzKJHylLxxDRyn8gJZZHM07/AjO3EjtzOWPvxGk7jXUCtBpHTlthycc7c659Y652YCrYC9nXP/DT80kViZmd4gcXTdo/N4nHkczGxOJ4Pl9OUe2vAtI7kxbmJQK0GkfKUmBzM7xMx2jtruBjwK3GpmzRIRnAhsuQxncRdSXQq4lCkspD2P05nG/EYPJrEHS7iPvmximy2uUbeulxCcg/XrNdgsUp6yWg7jgXwAMzsWGA5MAX4FcsIPTWqreGsyF89AasAmrmQci9mTKXQnn/pcyAz2YREP0iOmIJ6Zt+xmQYESgkhFlDXmkOacW+M/vxDI8buWZprZ/NAjk1qptJlH2/IbVzKe67mDv/ED73MY/bibWZxJvBIXWodZZOuUmRzMrK5z7k/gRKBXwPNEAivt3oRiTfmFvtxLf8bQnDxe4QQuYRqvcTxKCiLhKetLfjrwhpmtxpvO+haAme2B17UkUmnlJYUW/Ly57tF2rOdZziSbLD7YosTXX7SOgkjVKjU5OOeyzewVoCXwknOb61jWAa5ORHBSM5XWdQSwK99xAyO5ggdowB88ygUMYzCf0SHu8UoKIuEos3vIOfd+nH2LwwtHaqp4pbKjtWUJgxhON6ZgOKZyKcMZxNfsGfd4lcsWCZfGDiR0ZbUU9uNzBjOMC3mEAuqRQy9GcgPLaRVzbJ06XleUxhREwqfkIKGKROInhkP4kCyyOZtnWE9jRjGAO7mWn9h5i+PS0uChh9RCEEm0ILWVRALLzPR+4Uffo/AXxz94nZc4mQ85jGN4i5u4mVbkMojbYxJDeroSg0iyqOUgVSISgZ49IT8/3quO03meIQzlKN7lB3ZmACMZz5VblLfQNFSR1KHkIJUSiUD//pCXV/oxdSjkXJ5gCEM5kPnkkkEm9zGJnvzx17pRGlwWSUHqVpIKy8z0uotKSwx1KaAbD7GQ9jzGBTRiA5cxmT1YwlgyNyeG9HSv3pFqHYmkHrUcJJAgLYUGbKInk7iREbQml/l0oDOP8gTnUkTa5uPq14dJk5QQRFKZkoOUKxKBHj284nXxbMtv9GYc13MHLfmRdzmCq7iP5zmdkiUu1IUkUj2kXLeSmd1sZt+b2Xz/7/Rkx1QbRSLQvPlfM47iJYYdWMN/uIVcWjGKG1hIe47nVY7iHZ7nDKITg7qQRKqXVG053OmcG5XsIGqbIF1HADvyE9cxmkzupwm/8TRnMZQhfMhhm49RC0GkekvV5CAJVtZdzMV2Yzk3MJJ/M4H65G+ue/Q5+28+Jj0dxoxRUhCp7lKuW8nX18w+M7NJZrZDsoOpyYq7j8pKDO1YzER6spS29GYcD3Mxe/MlFzOdlen7b15hzTlYvVqJQaQmSErLwcxehhK3w3qygLHArYDzH+8Aesa5Ri/8NSYyMjJCi7UmK68Y3t/5jCEMpTOPkU99xtKHUQzgOzKoUwemTVEiEKmpzJX2zZACzKw1MMs5t19Zx3Xs2NHNnTs3MUHVEJEIXHpp/MRwKB+QRTZn8SzraML9ZHIn1/IzOwGaiipSU5jZPOdcx3ivpVy3kpm1jNrsBCxIViw1VSQC3buXTAyO43mVOZzEBxzOUbzDf7iFVuQymOGbE0N6uhKDSG2QigPSI8zsALxupWXAlUmNpoaJ7UpynMFzZJHNEbzPD+zM9YxiPFdijRtrxpFILZVyycE5d2myY6gpIhHIyoLcXO9+heiWQh0KOY+ZDGEoB/Apy2hFH+5nMj1onN6Q8ZpxJFKrpVy3klRcJAKtW3sJoG5d77FOHe/mtdxc75jixFCXAi5jMl+wL49yIQ3ZRHcepB1fM9760LNPQ804EpHUazlIcPFuWiss9B5LDjQ3ZOPmuketWM4nHMD5PMaTdKKINNLSYKrWThARn5JDNVNWV1E8jVm/ue7RzvzEOxxJH8Yym9MoLm9hpkV1RGRLSg7VRLxWQlmJYQfW0I+76cfdNOMXXuJkLiSLNzmW6JpHZtC7txKDiGxJyaEaiESgVy/YsKH8Y3fiR65jNH0YSxN+4ynOZihD+IhDY45VqQsRKY2SQ4op7jZavhyaNfP2lVcIDyCD3M11j+pRwCNcyDAGs4C/xxyrpCAi5VFySAGljSMESQp78hWDGM4lTMNhTKE7wxnIsrQ9KCyEVq0gO1uJQEQqRskhyUp2GQWtZrI/n26ue7SJhkxrkknT2wbw73678e/wwhWRWkLJIQmiu47q1Plr+mkQh/MeWWRzJs+x3pqw6MyBtJ9wLT123DG8gEWk1lFySLCSLYVgicFxAq+SRTYn8Bpr6qTz6bm30uGBvrRv2jTEaEWkttId0gmWlRVs1pHHcSbP8h5H8AonsY99ybyL76DZr8vo8Nj/gRKDiIREySHBli8v/5g0CrmQGXyedgDPchY78jNDmo3j9YnfcHDkOm8NThGRECk5hKy47lGdOt5j8fTUktLSoD75DEifxJqd92EGXdhvzwKYMoXdCxYzNO9KuvRomMjQRaQW05hDiEqOL+TmQr163mI5+fl/Hddsm4280Hkih7w2Ar77Dg48EO59HDp18rKKiEiC6ZsnRPHGFwoKoEkT7/6D7VjHsKa3s6Jeaw6ZcjVkZMDzz8O8eXDeeUoMIpI0ajmEqNTxhbw8ll11N9x9N6xdC6ec4mWSY49NZHgiIqXST9MqVN74ws78wAhuINdawS23wPHHw4cfwosvKjGISEpRy6GKlDW+0DJ/GTcygp5Moh4FLD+iC21yBkP79skNWkSkFEoOVSTe+MLuBV9yU4PhnE8Eh/F448vY9n8DOfu6tskJUkQkICWHKhI9vtCB+WSRzXnMZNMfDanX/yoYMICLd901eQGKiFSAxhwqoeTYQiTiTTQ6gneZxRnM50BO4SWGMZijd82Fu+4CJQYRqUbUcqig2LEFx/TLX+Gl9Gz25HVWk04Wt3EfV1HQqCk5w5Mbr4hIZajlUEHFYwtGEWfxNO9zOLP+OJntflrMvK6jOWa3XIZZFk1bNSUnR+soiEj1pORQQrwuo2grcgu5iOl8Sgee5hxasIpejKd14TccPO1aFi3flqIiWLZMiUFEqi91K0WJNx21Vy/vedfO+TB1Kl/XHU6bP5ewkH25hKnM4CIKqUurVsmLW0Skqik5RIk3HdVt2MDXV0+AQSNhxQq2a3MwF33/BI/mn43zG16NGnlLcYqI1BQ1tlupvO6heKKnozZhHQMZzjJac/Mv/aFNG3jhBdKXfsS/JnUio1UdzLwaSRpbEJGaxlzQRYtTWMeOHd3cuXM3b5fsHgLv1315X+KtW8Nvuavpx91czT3swFpmcyqTdhrCYz8eE94HEBFJAjOb55zrGO+1GtlyiNc9tGGDt79UK1cye9/rWUZr/sutvMoJHMxczm80m3PuUGIQkdqlRiaH0qqhxt3/7bfQpw+0acM+L41h1VGdOKnlQjrbTPJaHawuIxGplZKSHMyss5ktNLMiM+tY4rXBZrbEzL4ys39W5voZGQH2L1oE3btDu3YwaRJcdhksXkybt6fy8sp9NR1VRGq1ZLUcFgDnAm9G7zSzfYGLgPbAqcD9ZpZW0YtnZ3tjDNE2zyj65BM4/3yvIurjj0O/fvDNNzB+POy+eyU/johIzZKUqazOuUUAZlbypbOBGc65P4BvzWwJcCjwXkWuX/xrPyvL60rKyICc7u9wSiQbZs+G7beHIUOgf39o0WJrP46ISI2TamMOuwDfRW2v8PfFMLNeZjbXzOauWrUq5vWuXWHZt46iF15iWat/cMotR8NHH8HQod7dbbfdpsQgIlKK0FoOZvYysHOcl7Kcc0+XdlqcfXHn2jrncoAc8KayxhyQlwenneYlhF128SqjXnFFbH+TiIjECC05OOdOqsRpK4DdorZ3BVZWKoBmzbwb1664Arp1gwYNKnUZEZHaKNXKZzwDPGxmo4G/Ae2ADyt1JTN45JEqDE1EpPZI1lTWTma2AjgCeM7MXgRwzi0EHgW+AF4ArnLOFSYjRhGR2ixZs5WeBJ4s5bVsQGXsRESSKNVmK4mISApQchARkRhKDiIiEkPJQUREYig5iIhIDCUHERGJUSNWgjOzVUBusuPYSs2B1ckOIkFq02eF2vV59Vmrl1bOubhF5mpEcqgJzGxuacv11TS16bNC7fq8+qw1h7qVREQkhpKDiIjEUHJIHTnJDiCBatNnhdr1efVZawiNOYiISAy1HEREJIaSg4iIxFBySCFmNtLMvjSzz8zsSTNrmuyYwmJmnc1soZkVmVmNnA5oZqea2VdmtsTMBiU7njCZ2SQz+9nMFiQ7lrCZ2W5m9pqZLfL/G+6f7JjCoOSQWuYA+znn9gcWA4OTHE+YFgDnAm8mO5AwmFkacB9wGrAv0MXM9k1uVKF6EDg12UEkyJ/A9c65fYDDgatq4r9bJYcU4px7yTn3p7/5Pt4a2jWSc26Rc+6rZMcRokOBJc65b5xz+cAM4OwkxxQa59ybwJpkx5EIzrkfnHMf+8/XA4uAXZIbVdVTckhdPYHZyQ5CKm0X4Luo7RXUwC+Q2s7MWgMHAh8kOZQql5RlQmszM3sZ2DnOS1nOuaf9Y7Lwmq6RRMZW1YJ81hrM4uzTvPEaxMwaAzOBa5xz65IdT1VTckgw59xJZb1uZt2BM4ETXTW/CaW8z1rDrQB2i9reFViZpFikiplZPbzEEHHOPZHseMKgbqUUYmanAgOBs5xzG5Idj2yVj4B2ZtbGzOoDFwHPJDkmqQJmZsBEYJFzbnSy4wmLkkNquRdoAswxs/lmNi7ZAYXFzDqZ2QrgCOA5M3sx2TFVJX9iQV/gRbwBy0edcwuTG1V4zGw68B6wl5mtMLPLkx1TiI4CLgVO8P8/nW9mpyc7qKqm8hkiIhJDLQcREYmh5CAiIjGUHEREJIaSg4iIxFByEBGRGEoOkjLMLD1qauCPZva9/3ytmX2R4FjOiS6mZma3mFmFb+ozs9alVSo1s/Zm9qqZLTazpWb2PzOr8v8ny/osZvZ6Ta2KK1tHyUFShnMuzzl3gHPuAGAccKf//ACgqKrfz8zKqhBwDl411eLY/uuce7kK33sbvJvihjvn9gT+jlesL4zyz+cQ4meRmknJQaqLNDN7wK+f/5L/5YqZtTWzF8xsnpm9ZWZ7+/tbmdkr/toYr5hZhr//QTMbbWavAbfHO9/MjgTOAkb6LZe2/nnn+9c4xMzeNbNPzexDM2vitxDeMrOP/b8jy/k8FwPvOOdeAvDviO8L3OC/x81mNqD4YDNb4Bd5w8ye8uNdaGa9oo75zcyy/bjeN7Odyvss0czsFDN7z4//Mb92EGY23My+8P9Zjqr4vzqpjpQcpLpoB9znnGsPrAXO8/fnAFc75w4GBgD3+/vvBab4a2NEgLujrrUncJJz7vp45zvn3sX7VX+D35JZWnyiXwrjEaC/c64DcBKwEfgZONk5dxBwYYn3i6c9MC96h/8+21j5izz19OPtCPQzs3R//7bA+35cbwJXlPVZoplZc+D//H8uBwFzgevMrBnQCWjv/7O8rZzYpIZQ4T2pLr51zs33n88DWvu/bI8EHvPK3QDQwH88Am8xIYCpwIioaz3mnCss5/zS7AX84Jz7CKC4GqeZbQvca2YHAIV4CagsRvwqrfGquZbUz8w6+c93w0uceUA+MMvfPw84OcC1ih2O1/X0jv/Poj5eOYx1wCZggpk9F3V9qeGUHKS6+CPqeSGwDV7Ld60/LlGe6C/i3/3HipxfrLQv9WuBn4AO/nU3lXOdhcCxW1zYbHdgtXNurZn9yZYt+4b+McfhtVaOcM5tMLPXi18DCqIq+RZSsf+/DZjjnOsS84LZocCJeMUD+wInVOC6Uk2pW0mqLf9X+7dm1hm8aplm1sF/+V28LzOArsDbFTx/PV4RxJK+BP5mZof45zTxB7a3x2tRFOEVZUsrJ/wIcHTUrKFt8LqibvJfXwYc5L92ENDG37898IufGPbG+8VfntI+S7T3gaPMbA//PRuZ2Z5+62p759zzwDV4kwOkFlBykOquK3C5mX2K92u8eCnOfkAPM/sM78u6tFlApZ0/A7jBzD4xs7bFB/tLfl4I3OOfMwfvl/v9QHczex+vS+l3yuCc24g3UJxlZouB1XgD1MULPM0EmpnZfKAP3priAC8Adf3PdSvel3p54n6WEvGsAi4DpvvXfh/YGy+pzPL3vYHXQpJaQFVZRVKAmZ0DjAaOd87lJjkcESUHERGJpW4lERGJoeQgIiIxlBxERCSGkoOIiMRQchARkRhKDiIiEuP/Ac3LuQUv5hqGAAAAAElFTkSuQmCC\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:                Sun, 16 Aug 2020   Deviance:                     0.087389\n",
      "Time:                        18:00:44   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:                Sun, 16 Aug 2020   Deviance:                   1.0215e-05\n",
      "Time:                        18:00:44   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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
