3_SLR_Adequecy_Residuals_fkingfisher.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa13cd0f-de9a-4a9a-9177-7c53d9aea186",
   "metadata": {
    "tags": []
   },
   "source": [
    "The purity of oxygen produced by a fractional distillation process is thought to be related to the percentage of hydrocarbons in the main condenser of the processing unit. Twenty samples are shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cd90d453-fe1f-44f2-a453-51c7779907ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Purity (%)</th>\n",
       "      <th>Hydrocarbon (%)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>86.91</td>\n",
       "      <td>1.02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>89.85</td>\n",
       "      <td>1.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>90.28</td>\n",
       "      <td>1.43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>86.34</td>\n",
       "      <td>1.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>92.58</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>87.33</td>\n",
       "      <td>0.95</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>86.29</td>\n",
       "      <td>1.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>91.86</td>\n",
       "      <td>0.87</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Purity (%)  Hydrocarbon (%)\n",
       "0       86.91             1.02\n",
       "1       89.85             1.11\n",
       "2       90.28             1.43\n",
       "3       86.34             1.11\n",
       "4       92.58             1.01\n",
       "5       87.33             0.95\n",
       "6       86.29             1.11\n",
       "7       91.86             0.87"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Intercept      77.863284\n",
       "Hydrocarbon    11.801028\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "dat46 = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/46table.csv\")\n",
    "display(dat46.head(8))\n",
    "\n",
    "# Creating X and y\n",
    "X1 = dataset['Hydrocarbon (%)']\n",
    "y1 = dataset['Purity (%)']\n",
    "dat46 = pd.DataFrame(X1).join(pd.DataFrame(y1))\n",
    "# Change the column names\n",
    "dat46.columns =['Hydrocarbon', 'Purity']\n",
    "\n",
    "# calculate least square estimate coefficients using linear model\n",
    "model = sm.formula.ols('Purity ~ Hydrocarbon', dat46)\n",
    "lm46 = model.fit()\n",
    "display(lm46.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b0ddae9-7fe7-4246-8a4e-0e18da05c52b",
   "metadata": {
    "tags": []
   },
   "source": [
    "**a. Construct a normal probability plot of the residuals. Does there seem to be any problem with the normality assumption?**\n",
    "\n",
    "Looking at the normal probability plot, the residuals indicate a moderately skewed/non-normal distribution, therefore there is evidence the data set may violate the assumption of normally distributed errors required for model adequacy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a82ebcb5-6eb6-48d1-abe2-aebb34206c1e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept      77.863284\n",
       "Hydrocarbon    11.801028\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHpCAYAAABeNIDUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABNX0lEQVR4nO3deXicVfn/8fddlhZkCZZFkKX9atlLWxqQfZVVBCwiKCiIgsBXjAEVcCOKqCA/46igX1REZFWURQVEFFkExRSKLQIiUqSAUpYUCm1pm/P745lkJmuTNrPm/bquXOmcZ+aZO9M0/eTM/ZwTKSUkSZIklc6oShcgSZIk1TtDtyRJklRihm5JkiSpxAzdkiRJUokZuiVJkqQSM3RLkiRJJbZyJZ88Ii4FDgGeTyltmx97M3AtMA6YDbwvpfRy/tjZwEeApcAnUkq/XdZzrLvuumncuHGlKF+SJEnqMn369BdSSuv1dSwquU53ROwBzAcuLwrdFwAvpZS+HhFnAeuklM6MiK2Bq4EdgY2A24HNU0pLB3qOxsbG1NbWVtKvQ5IkSYqI6Smlxr6OVbS9JKV0F/BSj+HDgJ/k//wT4PCi8WtSSotSSk8C/yQL4JIkSVJVq8ae7g1SSs8B5D+vnx9/K/B00f3m5MckSZKkqlaNobs/0cdYn70xEXFSRLRFRNvcuXNLXJYkSZI0sIpeSNmP/0bEhiml5yJiQ+D5/PgcYJOi+20MPNvXCVJKlwCXQNbT3fP44sWLmTNnDgsXLhzeylUyY8aMYeONN2aVVVapdCmSJElDVo2h+ybgOODr+c83Fo1fFRHfJLuQcgJw//I8wZw5c1hzzTUZN24cEX1NoKuapJR48cUXmTNnDuPHj690OZIkSUNW0faSiLgauA/YIiLmRMRHyML2fhHxOLBf/jYppYeBnwF/B24F/ndZK5f0Z+HChYwdO9bAXSMigrFjx/rOhCRJqlkVnelOKb2/n0P79nP/84DzhuO5Ddy1xb8vSZJUy2rpQkpJkiSpJhm6B6HnBkLDtaHQnDlzOOyww5gwYQJve9vbaGpq4o033uCyyy7j4x//+LA8x4q44YYb+Pvf/951+4tf/CK33357BSuSJEmqTYbuZWhpaaG5ubkraKeUaG5upqWlZYXOm1Ji2rRpHH744Tz++OP84x//YP78+Xzuc58bhqp7W7JkyZAf0zN0f/nLX+ad73zncJYlSZI0Ihi6B5BSor29nVwu1xW8m5ubyeVytLe3r9CM9x/+8AfGjBnDhz/8YQBWWmklWltbufTSS3n99dd5+umnOfDAA9liiy340pe+BMBrr73Gu971LiZNmsS2227LtddeC8D06dPZc889mTp1KgcccADPPfccAHvttRef/exn2XPPPTnvvPMYN24cHR0dALz++utssskmLF68mB/84AfssMMOTJo0iSOOOILXX3+de++9l5tuuolPf/rTTJ48mSeeeILjjz+e6667DoDf//73TJkyhYkTJ3LCCSewaNEiAMaNG8c555zD9ttvz8SJE3n00UcBuPPOO5k8eTKTJ09mypQpvPrqq8v92kmSJNUaQ/cAIoLW1laamprI5XKMGjWKXC5HU1MTra2tK3Rx38MPP8zUqVO7ja211lpsuummLFmyhPvvv58rr7ySGTNm8POf/5y2tjZuvfVWNtpoIx566CFmzZrFgQceyOLFiznttNO47rrrmD59OieccEK32fL29nbuvPNOzjnnHCZNmsSdd94JwK9+9SsOOOAAVlllFaZNm8Zf//pXHnroIbbaait+9KMfscsuu3DooYfyjW98gxkzZvC2t72t65wLFy7k+OOP59prr2XmzJksWbKE733ve13H1113XR544AFOOeUULrzwQgAuvPBCLrroImbMmMHdd9/NaqutttyvnSRJUq0xdC9DZ/AutqKBG7JZ9L7O0Tm+3377MXbsWFZbbTWmTZvGPffcw8SJE7n99ts588wzufvuu1l77bV57LHHmDVrFvvttx+TJ0/mK1/5CnPmzOk631FHHdXtz52z49dcc03XsVmzZrH77rszceJErrzySh5++OEBa3/ssccYP348m2++OQDHHXccd911V9fxadOmATB16lRmz54NwK677srpp5/Ot7/9bdrb21l55WpcIl6SJKk0DN3L0NlSUqy4x3t5bbPNNrS1tXUbe+WVV3j66adZaaWVegXyiGDzzTdn+vTpTJw4kbPPPpsvf/nLpJTYZpttmDFjBjNmzGDmzJncdtttXY9705ve1PXnQw89lFtuuYWXXnqJ6dOns88++wBw/PHH893vfpeZM2dyzjnnLHM97GV97aNHjwaylpnOXvKzzjqLH/7whyxYsICddtqpq+1EkiRpJDB0D6C4h7upqYmOjo6uVpMVDd777rsvr7/+OpdffjkAS5cu5YwzzuD4449n9dVX53e/+x0vvfQSCxYs4IYbbmDXXXfl2WefZfXVV+fYY4/lU5/6FA888ABbbLEFc+fO5b777gOyLe77m6leY4012HHHHWlqauKQQw5hpZVWAuDVV19lww03ZPHixVx55ZVd919zzTX77L3ecsstmT17Nv/85z8B+OlPf8qee+454Nf7xBNPMHHiRM4880waGxsN3ZIkaViVarW54WLoHkBE0NDQ0K2Hu7PHu6GhYYVaTCKC66+/np///OdMmDCBzTffnDFjxvDVr34VgN12240PfvCDTJ48mSOOOILGxkZmzpzJjjvuyOTJkznvvPP4/Oc/z6qrrsp1113HmWeeyaRJk5g8eTL33ntvv8971FFHccUVV3RrOzn33HN5xzvewX777ceWW27ZNX700UfzjW98gylTpvDEE090jY8ZM4Yf//jHHHnkkUycOJFRo0Zx8sknD/j1futb32Lbbbdl0qRJrLbaahx00EHL+9JJkiR1U6rV5oZTVNtvAcOtsbEx9WzjeOSRR9hqq60GfY6e/df99WOrtIb69yZJkupfz86E1tbWXrfLldsiYnpKqbGvY17NNgh99VdLkiSp8ooXvcjlcuRyOYCyB+5lsb1EkiRJNa1Uq80NJ0O3JEmSalqpVpsbToZuSZIk1axSrjY3nOzpliRJUs3qb7U5YIVXmxtOhm5JkiTVtJaWlm6ry3UG72oJ3GB7ifKeffZZ3vve9w7Lufbaa6+u3TYPPvhg2tvbh+W8kiRJ/an21eac6a6wlBIpJUaNWv7ff5YsWcLKK6/YX+VGG23Eddddt0Ln6MvNN9887OeUJEmqNc50V8Ds2bPZaqutOPXUU9l+++15+umn+cY3vsEOO+zAdtttxznnnNN133PPPZctt9yS/fbbj/e///1ceOGFQDab/NnPfpY999yTXC7H9OnT2XPPPZk6dSoHHHAAzz33HADf/va32Xrrrdluu+04+uijAbjzzjuZPHkykydPZsqUKbz66qvMnj2bbbfdFoCFCxfy4Q9/mIkTJzJlyhTuuOMOAC677DKmTZvGgQceyIQJE/jMZz6zzK913LhxvPDCC11f84knnsg222zD/vvvz4IFC4Bsi/gDDzyQqVOnsvvuu7tFvCRJqjvOdAOlePdhWRfKPvbYY/z4xz/m4osv5rbbbuPxxx/n/vvvJ6XEoYceyl133cXqq6/OL37xCx588EGWLFnC9ttvz9SpU7vO0d7ezp133snixYvZc889ufHGG1lvvfW49tpr+dznPsell17K17/+dZ588klGjx7d1eZx4YUXctFFF7Hrrrsyf/58xowZ0622iy66CICZM2fy6KOPsv/++/OPf/wDgBkzZvDggw8yevRotthiC0477TQ22WSTQb0mjz/+OFdffTU/+MEPeN/73scvfvELjj32WE466SS+//3vM2HCBP7yl79w6qmn8oc//GGQr7QkSVL1M3Sz7IBcCpttthk77bQTALfddhu33XYbU6ZMAWD+/Pk8/vjjvPrqqxx22GGsttpqALz73e/udo6jjjoKyAL8rFmz2G+//QBYunQpG264IQDbbbcdxxxzDIcffjiHH344ALvuuiunn346xxxzDNOmTWPjjTfudt577rmH0047DYAtt9ySzTbbrCt077vvvqy99toAbL311jz11FODDt3jx49n8uTJAEydOpXZs2czf/587r33Xo488siu+y1atGhQ55MkSaoVhu4KedOb3tT155QSZ599Nh/72Me63afnzkr9nSOlxDbbbMN9993X6z6/+c1vuOuuu7jppps499xzefjhhznrrLN417vexc0338xOO+3E7bff3m22e6D1LEePHt3155VWWoklS5YM/IUO8NgFCxbQ0dFBQ0MDM2bMGPR5JEmSao093VXggAMO4NJLL2X+/PkAPPPMMzz//PPstttu/OpXv2LhwoXMnz+f3/zmN30+fosttmDu3LldoXvx4sU8/PDDdHR08PTTT7P33ntzwQUX0N7ezvz583niiSeYOHEiZ555Jo2Njb16qPfYYw+uvPJKAP7xj3/w73//my222KIkX/taa63F+PHj+fnPfw5kgf+hhx4qyXNJkiRVijPdVWD//ffnkUceYeeddwZgjTXW4IorrmCHHXbg0EMPZdKkSWy22WY0NjZ2tXYUW3XVVbnuuuv4xCc+wbx581iyZAmf/OQn2XzzzTn22GOZN29e125NDQ0NfOELX+COO+5gpZVWYuutt+aggw7quvAS4NRTT+Xkk09m4sSJrLzyylx22WXdZqmH25VXXskpp5zCV77yFRYvXszRRx/NpEmTSvZ8kiRJ5RbVsjVmqTQ2NqbONaM7PfLII2y11VYVqmho5s+fzxprrMHrr7/OHnvswSWXXML2229f6bIqopb+3iRJ0sgTEdNTSo19HXOmu8qddNJJ/P3vf2fhwoUcd9xxIzZwS5Ik1TJDd5W76qqrKl2CJEmSVtCIvZCy3ttq6o1/X5IkqZaNyNA9ZswYXnzxRYNcjUgp8eKLL/baxEeSJKlTz1xXbTlvRLaXbLzxxsyZM4e5c+dWuhQN0pgxY3pt4iNJkgTQ0tJCe3s7ra2tRES3VdtaWloqXR4wQkP3Kquswvjx4ytdhiRJklZQSon29nZyuRyQbS7Y3NxMLpejqamJlBIRUeEqR2joliRJUn2IiK5dvHO5XFf4bmpq6pr5rgYjcp1uSZIk1ZeUEqNGFS5X7OjoKHvgHmid7hF5IaUkSZLqR2cPd7Hm5uaqupjS0C1JkqSa1Rm4O3u4Ozo6aGpqIpfLVVXwtqdbkiRJNSsiaGho6NbD3dnj3dDQYE93udjTLUmSVP96rlJSiVVL7OmWJElSXesZsKtlhruToVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6JUmSpBIzdEuSJEklZuiWJEmSSszQLUmSJJWYoVuSJEkqsaoN3RHRFBGzIuLhiPhkfmxyRPw5ImZERFtE7FjhMiVJkqRlqsrQHRHbAicCOwKTgEMiYgJwAfCllNJk4Iv525IkSVJVW7nSBfRjK+DPKaXXASLiTuA9QALWyt9nbeDZypQnSZIkDV61hu5ZwHkRMRZYABwMtAGfBH4bEReSzdLv0teDI+Ik4CSATTfdtBz1SpIkSf2qyvaSlNIjwPnA74BbgYeAJcApQHNKaROgGfhRP4+/JKXUmFJqXG+99cpUtSRJktS3qgzdACmlH6WUtk8p7QG8BDwOHAf8Mn+Xn5P1fEuSJElVrWpDd0Ssn/+8KTANuJqsh3vP/F32IQvikiRJUlWr1p5ugF/ke7oXA/+bUno5Ik4EchGxMrCQfN+2JEmSVM2qNnSnlHbvY+weYGoFypEkSZKWW9W2l0iSJEn1wtAtSZIklZihW5IkSSoxQ7ckSZJUYoZuSZIkqcQM3ZIkSVKJGbolSZKkEjN0S5IkSSVm6JYkSZJKzNAtSZIklZihW5IkSSoxQ7ckSZJUYoZuSZIkqcQM3ZIkSVKJGbolSZKkEjN0S5IkSSVm6JYkSZJKzNAtSZIklZihW5IkSSoxQ7ckSZJUYoZuSZIkqcQM3ZIkSVKJGbo1oqWUBrwtSZI0HAzdGrFaWlpobm7uCtopJZqbm2lpaalsYZIkqe4YujUipZRob28nl8t1Be/m5mZyuRzt7e3OeEuSpGG1cqULkCohImhtbQUgl8uRy+UAaGpqorW1lYioZHmSpBEopdTt/5+et1Xbot5n9BobG1NbW1uly1CVSikxalThDZ+Ojg5/wEmSyq6lpYX29vauiZ/Od2AbGhpse6whETE9pdTY1zHbSzRidf5AK1bc4y1JUjnY8jgyGLo1IhX/QGtqaqKjo4OmpqZuP/AkSSqHzpbHzv+HRo0a1fX/ky2P9cP2Eo1YvpUnSaomtjzWvoHaS7yQUiNWS0tLt4tUOmca/AEnSSq3/loe/X+pftheohGt5w8yf7BJksrNlseRwZluSZKkCooIGhoauvVwdy5r29DQ4IRQnbCnW5IkqQq4Tnftc8lASZKkKmfLY30zdEuSJEklZuiWJEmSSszQLUmSJJWYoVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6JUmSpBIzdEuSJEklZuiWJElSzUsp5T93v10tDN2SJEmqaWeeeT4TJvyNCBg1Kgvczc3NtLS0VLq0LlUbuiOiKSJmRcTDEfHJovHTIuKx/PgFFSxRkiRJFbJ0KXzlKxABF1xwJk88MYnNNvs7//1vFrhzuRzt7e1VM+O9cqUL6EtEbAucCOwIvAHcGhG/ATYGDgO2Syktioj1K1imJEmSyuzmm+Fd7yrcXm89+PWvE1ddlQXtDTbIxpuammhtbSUiKlNoD1Et6b9YRBwJHJBS+mj+9heARUAjcElK6fbBnquxsTG1tbWVplBJkiSV3D//CdOmwcyZhbEf/xiOOy6b6YaspWTUqEITR0dHR9kDd0RMTyk19nWsWttLZgF7RMTYiFgdOBjYBNgc2D0i/hIRd0bEDhWtUpIkSSXx2mtwwglZqJ4wIQvc//u/sGBBdrHk8cd3D9zNzc3dHt/c3Fw1rSVQpaE7pfQIcD7wO+BW4CFgCVk7zDrATsCngZ9FH7/CRMRJEdEWEW1z584tX+GSJElabinBd7+bhek11shms3feGZ56qnBszJiejyn0cDc1NdHR0UFTUxO5XK6qgndV9nQDpJR+BPwIICK+CswBtgJ+mbJX7/6I6ADWBeb2eOwlwCWQtZeUs25JkiQNzT33wIEHZrPbAKusArfeCvvss+zHRgQNDQ3derhbW1sBaGhosKd7WSJi/ZTS8xGxKXAbsDNwFLBRSumLEbE58Htg0zTAF2FPtyRJUvV55hl4//vh7rsLY62t0NRUaBsZipRSt4Dd83Y5DNTTXbUz3cAvImIssBj435TSyxFxKXBpRMwiW9XkuIECtyRJkqrHokVw9tlZuO50zDFw8cWw1lordu6eAbtaZrg7VW3oTint3sfYG8CxFShHkiRJy+mKK+CDHyzc3moruP562GKLytVUblUbuiVJklS7ZsyAQw7J2kg6XX89HH54pSqqrKpcvUSSJEm158UX4dBDs57sKVOywH3OObB4cbb6yEgN3OBMtyRJUr+q4eK8ard0KZx3XhauO73rXXDZZbDuuhUrq+oYuiVJkvrQ0tJCe3t71zJ0netBNzQ00NLSUunyKu7Xv4Z3v7twe8MNs7Htt69cTdXM9hJJkqQeUkq0t7d322ClcwOW9vb2qtlwpdwefxy23TZrH+kM3D/5SdY68uyzBu6BONMtSZLUQ/EGK7lcjlwuB9BtA5aR4tVXs+3Xf/rTwthpp8EFF/TeHVL9q9rNcYaLm+NIkqTllVJi1KhCY0BHR8eICNwpwXe+k21U02nXXeHqq2GTTSpXV7UbaHMc20skSZL60NlSUqyz1aRe3XknrLYajBqVBe7Ro+GOO7IQfs89Bu4VYeiWJEnqobiHu6mpiY6ODpqamrr1eNeLOXNgt92yPu299oKFCyGXg46O7M977VXpCuuDPd2SJEk9RAQNDQ3derg7e7wbGhpqvsVk0SL4zGfg298ujH3wg3DRRbDmmpWrq57Z0y1JktSPelun+/LL4bjjCre32QZ++UvYfPPK1VRPBurpdqZbkiSpHz0Ddi0G7gceyDar+c9/CmM33pjtHKnysadbkiSpzrzwAhxySNanPXVqFri/9CVYsiS7KNLAXX7OdEuSJNWBJUvg3HPhy18ujL373fDjH8PYsZWrSxlDtyRJUg278UY4/PDC7be+FX71K5gypWIlqQ+2l0iSJNWYxx6DrbbK2kc6A/fll2etI3PmGLirkaFbkiSpBrz6Khx7bBa0t9wSHn0UmpuztbRTypb8U/WyvUSSJKlKpQTf+hacfnphbPfds+3Y3/rWipWl5WDoliRJqjJ33AEHHACLF2e33/QmuOWWLHCrNtleIg2DnptM1fumU5Kk4ffvf8POO2ftI/vskwXu73wn2459/nwDd60zdEsrqKWlhebm5q6gnVKiubmZlpaWyhYmSap6CxfCxz+eBe3NNoM//znbMfLVV7PWks5jqn2GbmkFpJRob28nl8t1Be/m5mZyuRzt7e3OeEuSekkpWzs7AlZbDS66CCZOhH/8Izt22WWwxhqVrlLDzZ5uaQVEBK2trQDkcjlyuRwATU1NtLa21uR2wZKk0mhrg4MPhrlzC2O//nW2RbvqX9T7TFxjY2Nqa2urdBmqcyklRo0qvHHU0dFh4JYkMXcufOhDcOuthbFzz4Wzz4aVVqpcXSqNiJieUmrs65jtJdIK6mwpKVbc4y1JGlmWLIHPfz5rH1l//SxwH3YYvPBC1j7y+c8buEciQ7e0Aop7uJuamujo6KCpqalbj7ckaWS4/vosaK+yCpx3HmyyCcyYkQXtG26AsWMrXaEqyZ5uaQVEBA0NDd16uDt7vBsaGmwxkaQ698gj2Sz2448Xxq68Ej7wgcrVpOpkT7c0DFJK3QJ2z9uSpPrxyitw8snZrpCdzjgDvvpVWHXVytWlyhuop9uZbmkY9AzYBm5Jqi8dHfDNb8KnP10Y23NPuOoq2GijytWl2mHoliRJ6sfOO2cb1nRac81sO/Zdd61cTapNXkgpSZJU5MorswsiIwqBe4ststnuV14xcGv5ONMtSZJGvOee67tN5MknYdy4spejOuRMtyRJGpFSKsxoFwfuiy/OjqVk4NbwMXRLkqQR5fjjs6BdtJEwEyYUgvYpp1SsNNUx20skSVLd+8MfYN99e48vXAijR5e/Ho08hm5JklSXXn0V1lqr93hbG0ydWv56NLLZXiJJkurKmDFZ+0hx4D7zzEL7iIFbleBMtyRJqnnnngtf/GLv8TrfeFs1xNAtSZJq0qxZMHFi7/GXX4aGhrKXIw3I9hJJklQzFi8uLPNXHLhvvrnQPmLgVjUydEuSpKo3dWoWtFddtTB2xBGFoH3QQZWrTRoM20skSVJV+slPsjW1e+royAK4VEsM3ZIkqWo88wxsvHHv8aeegk03LX890nCxvUSSJFVU8XbsxYH7kksK7SMGbtU6Q7ckSaqID3yg93bs22xTCNonnli52qThZnuJJEkqm9/9Dvbfv/f4okXdL5KU6o2hW5IkldS8eX0v4zdjBkyaVO5qpMqo2vaSiGiKiFkR8XBEfLLHsU9FRIqIdStUniRJWoZRo7L2keLA/bnPFdpHDNwaSaoydEfEtsCJwI7AJOCQiJiQP7YJsB/w78pVKEmS+vKFLxQuiizegr0zaH/lK5WrTaqkqgzdwFbAn1NKr6eUlgB3Au/JH2sFPgOk/h4sSZLK56GHCkG7OFS3txfCtjTSVWvongXsERFjI2J14GBgk4g4FHgmpfTQQA+OiJMioi0i2ubOnVuOeiVJGlHeeKMQtCdPLozfdlshaK+9dsXKk6rOkC+kjIgtgU2BdYEFwPPAzJTSK8NVVErpkYg4H/gdMB94CFgCfA7o45rnXo+/BLgEoLGx0d+vJUkaJhMnwqxZ3ceOPhquvroy9Ui1YlChOyL2AT4CvJMsbPfUEREPAtcBl6aUXljRwlJKPwJ+lH/+rwL/BY4BHops79eNgQciYseU0n9W9PkkSVLffvjDvtfMdjt2afAGDN0RMQ04D9gcCOAZ4EbgP8BLwGrAWGBLYDLQCHwpIi4HvphS+u/yFhYR66eUno+ITYFpwM4ppVzR8dlA43AEfEmS1N3TT/e9C+TTT/e9TbukgfUbuiPiLmA34BHgbOCalFK/K4ZExKrA3sBxwLHA0RHxwZTSTctZ2y8iYiywGPjflNLLy3keSZI0CCl13x2y02WXwXHHlb0cqa4MNNO9JnD4YENzSukN4LfAbyNifeCzwBbLW1hKafdlHB+3vOeWJEkFRx4J113XfWzKFHjggcrUI9WjfkN3SmnK8p40pfQ88MnlfbwkSSqtW26Bgw/uPf7GG7DKKuWvZyhSSkRRM3nP21I1qtYlAyVJ0jBrby8s81ccuP/2t8Iyf9UeuFtaWmhubiblF/9OKdHc3ExLS0tlC5OWYYVCd0SMiYht8h9jhqsoSZI0fDqD9jrrFMa+9KVC0J44sXK1DUVKifb2dnK5XFfwbm5uJpfL0d7e3hXEpWo05HW6ASJiZbJVTU4DRueHF0VEDvhCfhdJSZJUIWefDV//evexVVeFRYsqU89wiAhaW1sByOVy5HLZomZNTU20trbaYqKqFsvzW2FEfAs4GbgCmA6MAd4N7Am0ppQ+NYw1rpDGxsbU1tZW6TIkSSq56dOhsbH3+CuvwJprlr+eUkkpMapomZWOjg4Dt6pCRExPKfXxr3D520s+CJyRUvpoSul7KaXWlNI+ZDtIfnB5C5UkSUOzaFGhfaQ4cP/+94X2kXoL3M3Nzd3Ginu8pWo1YOiOiHsjoq9l/9YAHu9j/J/5Y5IkqYQmTMiC9piiK6o+9KFC0N5nn8rVVirFPdxNTU10dHTQ1NTUrcdbqlbL6uleAsyIiC8BF6SUOvLj9wEXRsTHgAfJ+rrfTbYxzr2lKlaSpJHs+9+HU07pPT5StmOPCBoaGrr1cHf2eDc0NNhioqq2zJ7uiDiN7KLJx4ATUkozI2I74BbgLcV3Jdsm/oCU0t9LVO+Q2dMtSapls2fD+PG9x599FjbcsOzlVAXX6Va1WqGe7pTSd4DtgHagLSJagL8Dbye7mPJbQA44EZhQTYFbkqRa1DlzHdE9cF9xRaF9ZKQGbqBXwDZwqxYMasnAlNJsYL+IOBG4AHgP2az3D0pYmyRJI8phh8FNN3Uf22knuO++ytQjafgMafWSfMjeFngauC8ivh4Ro5fxMEmS1I9f/aowq10cuBcvzma0DdxSfRh06I6IsQAppWdSSocAHwE+CjwUEbuUqD5JkurOSy8VgvahhxbGH3640D6y8nJtXyepWi1rycBVIuL8iHgVeD4iXo2ICyJilZTST4FtgFnAXRHxrYhYvRxFS5JUizqD9tixhbGvfrUQtLfeunK1SSqtZc10fx74NHA/8I385zOALwCklP6bUnovcHT+Y2ZE7F26ciVJqi1nnFEI253WXLMQtM8+u3K1SSqfZb15dSxwW0rpwM6BiPgtcAzwxc6xlNJ1EfEH4DvA7cBKJahVkqSacP/98I539B5/9VVYwy3kpBFpWTPd6wIze4z9LT/eTUrppZTSMWSb5EiSNKIsXFiY0S4O3H/8Y2FW28AtjVzLCt0zgPdGxKYAEbExcER+vE8ppZuHqzhJkqrduHFZ0F5ttcLYRz9aCNp77lmx0iRVkWW1l5wB/B54IiLmAusBC4CjSl2YJEnV6rvfhdNO6z0+UrZjlzR0A4bulFJbRGwFfBDYFPg3cEVK6ZlyFCdJUrV44gl4+9t7j//nP7DBBuWvR1JtWeYqoCmlZ4Hzy1CLJElVpaMDVupjaYCrr4ajjy5/PZJq15B2pJQkaSQ48MCsTaQ4cO+xR6FP28Ataaj6Dd0RccSKnDgiNoyInVfkHJIklcv11xdWH/ntbwvjS5ZkQfvOOytXm6TaN9BM988jYnpEHBURowd7wojYIiJagX8C71zhCiVJKpEXXigE7WnTCuOPPlqY1e6rvUSShmqgnu59gVbgamBeRNwI/AloA54DXgbGAGOBLYGdgAOARuAN4NvAt0pVuCRJy6uvFUYuuAA+/eny1yJpZOg3dKeU7oiIKcD7gf8FPkS2ikl/AmgHckAupfTUMNYpSdIKOe20bKm/YmPHZrPdklRqy1oyMAFXAVdFxBZk7SK7kS0fOJZsze7nyXap/CPwh5TSglIWLEnSYN13H+yyS+/x116D1Vcvfz2SRq5lLhnYKaX0GPAYcFHpypEkacUsWNB3oL7nHth11/LXI0ngkoGSpDqx4YZZr3Zx4D7llMIFkQZuSZU06JluSZKqzTe/CWec0Xs8pfLXIkkDMXRLwyClRBQth9DztqTh8/jjsPnmvceffx7WW6/89UjSYNheIq2glpYWmpubSfmptZQSzc3NtLS0VLYwqY4sXVpYT7s4cF93XaF9xMAtqZoZuqUVkFKivb2dXC7XFbybm5vJ5XK0t7d3BXFJy2fffbOgvXLR+7LvfGchaB+xQnsnS1L52F4irYCIoLW1FYBcLkculwOgqamJ1tZWW0yk5XDddXDkkb3Hlyxxd0hJtSvqfSausbExtbW1VboM1bmUEqNGFd446ujoMHBLQ/D887DBBr3H//EPmDCh/PVI0vKIiOkppca+jg26vSS/O6WkHjpbSooV93hL6ltKhT7t4sDd2lpoHzFwS6oXQ+npnh4Rf4mIEyLCfbwk6NbD3dTUREdHB01NTd16vCV1d/LJWdAuenOIjTYqBO1PfrJipUlSyQylp/tm4ADgB8A3I+KnwCUppZklqUyqARFBQ0NDtx7uzh7vhoYGW0ykvHvugd137z3++uuw2mrlr0eSym1IPd0R8VbgROAEYGMgAX8G/g/4WUppYSmKXBH2dKscXKdb6u2112CNNXqP33cf7LRT+euRpFIblp5ugJTSMymlFmAccBhwC7Aj8GPgmYhojYitVqxcqfb0DNgGbo1kb35z1j5SHLg/8YlC+4iBW9JItFzrdKeUOlJKv0opHUIWwL8MvAF8ApgVEX+MiPcOX5mSpGp2wQWFiyJffrkw3hm086tpStKINRzrdG8DbAeMBQJ4Adgd2D0iZgBHpJRmD8PzSJKqyKOPwlZ9vLf5wgswdmz565GkarZcM90RsX5EnBURT5C1mBwO/BGYBrwFeDtZn/dk4OLhKFSSVHlLlhRmtIsD9/XXF2a1DdyS1NuQZrojYl/gY2T93KsALwPfAr6XUvpn0V2fBE6NiNHA+4anVElSpey+e7YCSbGDD4bf/KYy9UhSrRl06I6Ix4H/IWshaSObwb5mGSuWPA68aYUqlCRVxNVXwwc+0Ht86dLua2xLkpZtKDPdbwUuAy5OKU0f5GOuBO4balGSpMr4z39gww17jz/xBPzP/5S/HkmqF0MJ3RullNqHcvKU0tPA00OqSJJUVin1PXP9ne/Axz9e/nokqR4N5Q3Cb0bEoQPdISIOiYhLV7CmznM1RcSsiHg4Ij6ZH/tGRDwaEX+LiOsjomE4nkuSRqKPfKT3duzjxxcuiDRwS9LwGUroPp5sNZKBTAKOW95iOkXEtmQ7X+6YP+chETEB+B2wbUppO+AfwNkr+lySNJL88Y+F1UcuLZoiWbAgC9r/+lfFSpOkujbcl8KMBpYOw3m2Av6cUno9pbQEuBN4T0rptvxtyLaf33gYnkuS6tr8+YWgvffehfH77y/Mao8ZU7n6JGkkGGroTv0dyC8PuAfwnxWqKDML2CMixkbE6sDBwCY97nMC2RrhfdVyUkS0RUTb3Llzh6GcoUkpDXhbksrhTW/KgvaaaxbGzjijELR32KFytUnSSDPghZQR0fONxuaI+HAfd10JWI9spvv7K1pUSumRiDifrJ1kPvAQ0DnDTUR8Ln/7yn4efwlwCUBjY2NZE29LSwvt7e20trYSEaSUaG5upqGhgZaWlnKWImkEOu88+Pzne4/7u78kVdayZrpHka3LHWSz3NHPx2JgJnA+8OnhKCyl9KOU0vYppT2Al8jW/CYijgMOAY5JVTaFnFKivb2dXC5Hc3NzV+DO5XK0t7c74y2pJB5+uNA+Uhy4X3yxMKstSaqsAWe6U0rjOv8cER1Aa0rpy6UuKv9866eUno+ITcm2l985Ig4EzgT2TCm9Xo46hiIiaG1tBSCXy5HL5QBoamrqmvmWpOGwZAmsskrv8V/9Cg45pPz1SJIGFoOdfY2IPYHZKaWnSltS1/PdDYwlm0U/PaX0+4j4J1kLy4v5u/05pXTyQOdpbGxMbW1tpS22h5QSo4rW4Oro6DBwSxoWO+wAPX+kHX44XH99RcqRJBWJiOkppca+jg16c5yU0p3DV9Kgnm/3PsbeXs4alkdnS0mx5uZmZ7olLbef/hQ+9KHe427HLkm1o9/QHRGdP+KvTym9WnR7mVJKl69wZTWouIe7s6Wk8zZg8JY0aM88Axv3sSjqk0/CuHFlL0eStIIGmum+jOziyT8DrxbdHkjnBZcjMnRHBA0NDd16uDt7vBsaGgzckgbU33bs3/8+fOxj5a9HkjR8+u3pjojjyQL0L/Mz3Z23lyml9JPhKnBFVaqnuzhg97wtScWOPRau7LEA6hZbwKOPVqYeSdLyWa6e7pTSZQPdVv96BmwDt6Sebr8d9tuv9/jChTB6dPnrkSSV1qAvwYmISyOiedn3lCT15ZVXCutpFwfu6dML62kbuCWpPg3luvcPAOuXqhBJqlcrr5wF7bXXLoydfXYhaG+/feVqkySVx6CXDARmY+iWpEFpaYEvfan3uLtDStLINJTQfRVwckSsk1J6uVQFSVKt+tvfYNKk3uMvvwwNDWUvR5JURYbSXvI1oA24IyIOiYgNSlSTJNWMxYsLfdrFgfuWWwrtIwZuSdJQZroX5j8HcCP0uypHSikN5bySVHMmTcpmtosdeST87GeVqUeSVN2GEo7vZpDrdEtSPbr0UvjIR3qPd3RkM92SJPVn0KE7pbRXCeuQpKr09NOw6aa9x//9b9hkk/LXI0mqTUPp6ZakESGlQp92ceD+4Q8LfdoGbknSUBi6JSnvqKOyoD2q6CfjxImFoN1Xa4kkSYMx6PaSiPjiIO+aUkrnLmc9klRWv/0tHHhg7/FFi2DVVctfjySpPg3lQsqWAY51XmAZ+T8buiVVrfZ2WGed3uMzZvS9zrYkSStqKKF7737GG4AdgE8AvwG+v4I1SVJJ9LXCyBe+AF/+cvlrkSSNLENZveTOAQ7fGBHXAvcD16xwVZI0TD77Wfja17qPjRoFS5dWph5J0sg0bJvYpJRmRsSNwGfJb54jSZXw4IOw/fa9x+fNg7XWKn89kiQN9+ol/wa2HeZzStIyvfFGYZm/4sD9u98VVh8xcEuSKmW4Q/c7gAXDfE5J6tdWW2VBe/TowtgHPlAI2u98Z+VqkySp01CWDOxjT7auc2wCnAjsBvxsGOqSpH793//BySf3Hnc7dklStRpKT/dsCksD9iWAx4FPrUhBktSXRx6BrbfuPT5nDrz1reWvR5KkoRhK6L6cvkN3B/Ay2colN6aUFg1HYZLU0QErrdR7/Cc/gQ99qPz1SJK0vIayZODxJaxDkrpssAE8/3zv8TTQe22SJFWxQV1IGRGbRsQRETEtIjYpdVGSRp4f/KCw+khx4F60qHBRpCRJtWqZM90RcSHwSbKebYAUEa0ppU+XsjBJ9e+//4W3vKX3+J/+BLvsUv56JEkqlQFnuiPiA8DpZIH7UeCx/J9Pj4j3l748SfWoc0a7OHAfc0xhRtvALUmqN8tqL/kIsAR4Z0ppm5TS1sABZBdPfqTUxUmqHwccUAjbxTqD9hVXVKYuSZLKYVmhezvghpTSHZ0DKaXbybZ5n1zCuiTVgd/9rhC0b7utMP7yy/ZpS5JGlmWF7nXIWkp6ehRoGPZqpGGUeiS6nrdVGq+/Xgja++9fGL/22kLQbmioWHmSJFXEsi6kHAUs7mN8MYULK6Wq09LSQnt7O62trUQEKSWam5tpaGigpaWl0uXVpb52gtx2W5g5s/y1SJJUbQazZKDTg6opKSXa29vJ5XI0Nzd3Be5cLkd7e7sz3sOoubnvPu2OjmxG28AtSVImBgogEdHB0EN3SikNZafLkmpsbExtbW2VLkNlVhy0OzU1NXXNfGv5zZoFEyf2Hv/Xv2D8+PLXI0lStYiI6Smlxj6PDSJ0D1lKaVCb7pSDoXvkSikxalThW7Gjo8PAvZz62479a1+Ds84qfz2SJFWjgUL3gOE4pTRqeT5K82VIg9c5012ss9VEg9fQkLWO9AzcnRdEGrglSRocA7LqTnFrSVNTEx0dHTQ1NXXr8Vb/Lr640Kc9b15h/I03XOZPkqTlVTW919JwiQgaGhq69XC3trYC0NDQYItJH557DjbaqPf4X/4CO+5Y/nokSao3A/Z01wN7ukeulFK3gN3ztvpe5u/44+HHPy57KZJUlfy/REOx3D3dUi3r+UPRH5KZvfceeDt2A7ckZVpaWrq1JXa2L7rfg5aHoVsaAW6+uRC0//jHwnh7u33aktQX93zQcLOnW6pTr70Ga6zRe/wXv4Bp08pfjyTVkuLrgXK5XNe+D+75oOVlT7dUZ/r6f2DqVPCfgSQNnXs+aCjs6Zbq3Mc/PvB27AZuSRo693zQcDJ0SzVqxoxC0L7oosL4U08V+rSdjJGk5eOeDxpu9nRLNWTpUli5j3+13/gGfOpT5a9HkuqVez5ouFVtT3dENAEnAgH8IKX0rYh4M3AtMA6YDbwvpfTyQOexp1v1YPTobEfInqr0n68k1Q3X6dZQ1FxPd0RsSxa4dwQmAYdExATgLOD3KaUJwO/zt6W69O1vF9pHigP34sUu8ydJ5eKeDxou1dpeshXw55TS6wARcSfwHuAwYK/8fX4C/BE4swL1SSUxZw5ssknv8ba2bAUSSZJUm6pyphuYBewREWMjYnXgYGATYIOU0nMA+c/rV7BGaUh6tnIV3+6c0S4O3CeeWJjRNnBLklTbqnKmO6X0SEScD/wOmA88BCwZ7OMj4iTgJIBNN920JDVKQ9HS0kJ7e3vXxTgpJTba6Cn+859xve5r20j/7K2UJNWqap3pJqX0o5TS9imlPYCXgMeB/0bEhgD5z8/389hLUkqNKaXG9dZbr3xFS30o3kp42rTvEQGjRkW3wP3KK/ZpL0tLS0u3Zbo6l/NqaWmpbGGSJA1C1YbuiFg//3lTYBpwNXATcFz+LscBN1amOmnwFi0KcrlvAYkbbji1a/zd7/4BHR2JlGDNNStWXk0o/sWlM3h3rp/b3t7uermSpKpXzUsG3g2MBRYDp6eUfh8RY4GfAZsC/waOTCm9NNB5XDJQlbLbbvCnP/Uc/TxwHuBWwkNVHLQ7Fa+fK0lSpQ20ZGDVhu7hYuhWOf34x3DCCd3H1l8f/vMfA+NwSCkxalThDTp/cZEkVZOaW6dbqiWzZxdWHykO3C+9lPVoFwdutxJefp0z3cV8/SRJtcLQLS2HpUsLQXv8+ML4bbcVLohcZ51srL+thJuamtxKeJCKW0v8xUWSVIuqcslAqVoddxxcfnn3sY9+FH7wg4Ef19LS0m15u87gbeAenP5+cQH8xUWSVBPs6ZaW4dZb4aCDeo8vXQqjfK+orFynW5JUzQbq6XamW+rDiy/Cuuv2Hn/qKXC/pcrpGbAN3JKkWuE8nZSXEjQ0ZH3axYH7pz8t9GkbuCVJ0vIwdGvEa2khv0skzJuXje27byFoH3tsRcuTJEl1wPYSjUhtbbDDDr3HFy6E0aPLX48kSapvhm6NGAsWwOqr9x5/6CHYbrvy1yNJkkYO20tU97bdNmsfKQ7cX/96oX3EwC1JkkrNmW7VpfPPh7PO6j628cbw9NOVqUeSJI1shm7VjUcega237j3+wgswdmz565EkSepke4lq2pIlhe3YiwP3jTcW2kcM3JIkqdIM3apJu+ySBe1VVimMHXJIIWgfemjlapMkSerJ9hLVjKuugmOO6T3uduySJKnaGbpV1Z57DjbaqPf4v/4F48eXvx5JkqTl4fygqk5KhT7t4sB90UWF9hEDtyRJqiWGblWN448vbMfe6e1vLwTtU0+tWGmSJEkrxPYSVdQdd8A++/Qedzt2SZJUTwzdKrvFi2HVVXuPt7XB1Knlr0eSJKnUbC9R2fziF1n7SHHg/vSnC+0jBm5JklSvnOlWSf3979ma2U88URi76ip4//srV5MkSVK5OdOtYTdvHhx9dDarvc02WeA+4wxYtCib0TZwS5KkkcaZbg2Ljg648EI488zC2N57w5VXwoYbVq4uSZKkamDo1gq5/XbYf/9sBhtgzTXhlltg110rW5ckSVI1sb1EQzZ7NuywQ9Y+st9+WeD+3vey2e5XXjFwS5Ik9WTo1qAsWAAnn5wF7fHjs+X9PvIRmD8/C92dxyRJktSb7SXqV0rwwx/CSScVxqZMgeuug//5n8rVJUmSVGsM3erlL3+Bgw6Cl18ujN1yCxx4YOVqkiRJqmW2lwiA//4368+OgJ12ygL3174GS5dmM94GbkmSpOVn6B7BFi+Gs87KgvZb3pKtRHLEEfDSS1nQPussGOV3iCRJ0gqzvWQE+vnP4X3vK9wePx5uvBEmTqxcTZIkSfXM0D1CzJqVbcf+5JOFsWuugaOOqlxNkiRJI4XNA3WsvR2OPDJrH5k4MQvcn/lMYTt2A7ckSVJ5ONNdZzo64IIL4OyzC2P77gtXXJH1bUuSJKn8DN114rbb4IADCrfXWQduvjlbiUSSJEmVZXtJDfvXv2D77bP2kc7A/X//l812v/SSgVuSJKlaGLprzOuvZztERsDb3gYPPggnngivvZb1aXcekyRJUvWwvaQGpJTNYJ9ySmFs6tRs6b/x4ytXlyRJkgbH0F3F7rsv24593rzC2K23du/dliRJUvWzvaTKPPcc7LNP1iKyyy5Z4D7//MJ27AZuSZKk2mPorgJvvJGtnx0BG20Ed9yRra/98stZ0P7MZ9yOXZIkqZbZXlJB11wD739/4fbb3pZtx77NNpWrSZIkScPP0F1mM2fCu98NTz1VGPvZz7KZbUmSJNUnmxbK4OWX4YgjsvaR7bbLAvfZZ2dtJSkZuCVJkuqdM90lsnRpdgHk5z5XGNtvP/jpT2GDDSpXlyRJksrP0F0iU6ZkrSRjx8JvfgPveEelK5IkSVKlVG3ojohm4KNAAmYCHwa2BL4PjAGWAKemlO6vWJEDePDBbMURd4eUJElSVfZ0R8RbgU8AjSmlbYGVgKOBC4AvpZQmA1/M365KK61k4JYkSVKmKkN33srAahGxMrA68CzZrPda+eNr58ckSZKkqlaV7SUppWci4kLg38AC4LaU0m0R8TTw2/yxUcAufT0+Ik4CTgLYdNNNy1S1JEmS1LeqnOmOiHWAw4DxwEbAmyLiWOAUoDmltAnQDPyor8enlC5JKTWmlBrXW2+9cpUtSZIk9akqQzfwTuDJlNLclNJi4Jdks9rH5f8M8HNgxwrVJ0mSJA1atYbufwM7RcTqERHAvsAjZD3ce+bvsw/weIXqkyRJkgatWnu6/xIR1wEPkC0N+CBwSf5zLn9x5ULyfduSJElSNavK0A2QUjoHOKfH8D3A1AqUI0mSJC23am0vkSRJkuqGoVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6JUmSpBIzdEuSJEklZuiWJEmSSszQLUmSJJWYoVtllVIa8LZUK/xeliQNhaFbZdPS0kJzc3NXOEkp0dzcTEtLS2ULk4bI72VJ0lAZulUWKSXa29vJ5XJdYaW5uZlcLkd7e7uzhKoZfi9LkpZH1Pt/EI2Njamtra3SZQi6hZNOTU1NtLa2EhEVrEwaGr+XJUl9iYjpKaXGPo8ZulVOKSVGjSq8wdLR0WFIUU3ye1mS1NNAodv2EpVN5+xgseK+WKlW+L0sSRoqQ7fKovjt+KamJjo6OmhqaurWFyvVAr+XJUnLY+VKF6CRISJoaGjo1vfa2toKQENDg2/Lq2b4vSxJWh72dKusUkrdQknP21Kt8HtZktSTPd2qGj1DiSFFtcrvZUnSUBi6JUmSpBIzdEtl4rbhkiSNXIZuqQzcNlySpJHN0C2VmNuGS5IklwyUSqx4SblcLte1dbjbhkuSNHK4ZKBUJm4bLklSfXPJQKnC3DZckqSRzdAtlZjbhkuSJHu6pRJz23BJkmRPt1QmbhsuSVJ9s6dbqgJuGy5J0shl6JYkSZJKzNAtSZIklZihW5IkSSoxQ7ckSZJUYoZuSZIkqcQM3ZIkSVKJGbol1Yye+wrU+z4DkqT6YeiWVBNaWlpobm7uCtopJZqbm2lpaalsYZIkDYKhW1LVSynR3t5OLpfrCt7Nzc3kcjna29ud8ZYkVb2VK12AJC1LRLD22mszefJkcrkcuVwOgMmTJ7P22mu7u6ckqeo50y2p6qWUmDdvHjNmzOg2PmPGDObNm+dMtySp6hm6JVW9iOCb3/wmkydP7jY+efJkvvnNbzrTLUmqeoZuSVUvpcTpp5/e50z36aef7ky3JKnqGbolVb3inu5i9nRLkmqFoVtS1Svu6W5qaqKjo4OmpiZ7uiVJNcPVSyRVvYigoaGBpqYmWltbiQhaW1sBaGhocKZbklT1olpniCKiGfgokICZwIdTSgsj4jTg48AS4Dcppc8MdJ7GxsbU1tZW8nollV5KqVvA7nlbkqRKiojpKaXGvo5V5Ux3RLwV+ASwdUppQUT8DDg6Ip4CDgO2Syktioj1K1qopLLqGbAN3JKkWlHNPd0rA6tFxMrA6sCzwCnA11NKiwBSSs9XsD5JkiRpUKoydKeUngEuBP4NPAfMSyndBmwO7B4Rf4mIOyNih0rWKUmSJA1GVYbuiFiHrI1kPLAR8KaIOJZs9nsdYCfg08DPoo/3lyPipIhoi4i2uXPnlrFySZIkqbeqDN3AO4EnU0pzU0qLgV8CuwBzgF+mzP1AB7BuzwenlC5JKTWmlBrXW2+9shYuSZIk9VStofvfwE4RsXp+Jntf4BHgBmAfgIjYHFgVeKFSRUqSJEmDUZWrl6SU/hIR1wEPkC0N+CBwCdnygZdGxCzgDeC4VK1rHkqSJEl5VRm6AVJK5wDn9HHo2HLXIkmSJK2Iam0vkSRJkuqGoVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6JUmSpBIzdEuSJEklZuiWJEmSSszQLUmSJJWYoVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6VVYppQFvS5Ik1SNDt8qmpaWF5ubmrqCdUqK5uZmWlpbKFiZJklRihm6VRUqJ9vZ2crlcV/Bubm4ml8vR3t7ujLckSaprK1e6AI0MEUFraysAuVyOXC4HQFNTE62trUREJcuTJEkqqaj3GcbGxsbU1tZW6TKUl1Ji1KjCGywdHR0GbkmSVBciYnpKqbGvY7aXqGw6W0qKFfd4S5Ik1StDt8qiuIe7qamJjo4OmpqauvV4S5Ik1St7ulUWEUFDQ0O3Hu7OHu+GhgZbTCRJUl2zp1tllVLqFrB73pYkSapV9nSravQM2AZuSZI0Ehi6JUmSpBIzdEuSJEklZuiWJEmSSszQLUmSJJWYoVuSJEkqMUO3JEmSVGKGbkmSJKnEDN2SJElSiRm6JUmSpBIzdEuSJEklZuiWJEmSSszQLUmSJJWYoVuSJEkqMUO3JEmSVGKRUqp0DSUVEXOBpypdRxmsC7xQ6SLqnK9xefg6l56vcXn4Opeer3F5+DoP3mYppfX6OlD3oXukiIi2lFJjpeuoZ77G5eHrXHq+xuXh61x6vsbl4es8PGwvkSRJkkrM0C1JkiSVmKG7flxS6QJGAF/j8vB1Lj1f4/LwdS49X+Py8HUeBvZ0S5IkSSXmTLckSZJUYobuGhIRB0bEYxHxz4g4q4/j60TE9RHxt4i4PyK2rUSdtSwiLo2I5yNiVj/HIyK+nf87+FtEbF/uGuvBIF7nLSPivohYFBGfKnd99WAQr/Ex+e/hv0XEvRExqdw11oNBvM6H5V/jGRHRFhG7lbvGWres17jofjtExNKIeG+5aqsng/he3isi5uW/l2dExBfLXWOtM3TXiIhYCbgIOAjYGnh/RGzd426fBWaklLYDPgTkyltlXbgMOHCA4wcBE/IfJwHfK0NN9egyBn6dXwI+AVxYlmrq02UM/Bo/CeyZ/3lxLvZsLq/LGPh1/j0wKaU0GTgB+GEZaqo3lzHwa9z5f+T5wG/LUVCduoxlvM7A3SmlyfmPL5ehprpi6K4dOwL/TCn9K6X0BnANcFiP+2xN9gOelNKjwLiI2KC8Zda2lNJdZIGvP4cBl6fMn4GGiNiwPNXVj2W9ziml51NKfwUWl6+q+jKI1/jelNLL+Zt/BjYuS2F1ZhCv8/xUuHjqTYAXUg3RIH4uA5wG/AJ4vvQV1adBvs5aAYbu2vFW4Omi23PyY8UeAqYBRMSOwGb4H+lwG8zfg1RrPgLcUuki6lVEvCciHgV+QzbbrWEUEW8F3gN8v9K1jAA7R8RDEXFLRGxT6WJqjaG7dkQfYz1nTL4OrBMRM8h+638QWFLiukaawfw9SDUjIvYmC91nVrqWepVSuj6ltCVwOFkrj4bXt4AzU0pLK11InXuAbIvzScB3gBsqW07tWbnSBWjQ5gCbFN3eGHi2+A4ppVeAD0N2wR9Zz+aT5SpwhFjm34NUKyJiO7Ie44NSSi9Wup56l1K6KyLeFhHrppReqHQ9daQRuCb7b491gYMjYklK6YaKVlVn8hmj8883R8TFfi8PjTPdteOvwISIGB8RqwJHAzcV3yEiGvLHAD4K3FX8j0TD4ibgQ/lVTHYC5qWUnqt0UdJQRcSmwC+BD6aU/lHpeupVRLw9PwlCfrWjVQF/wRlGKaXxKaVxKaVxwHXAqQbu4RcRbyn6Xt6RLEP6vTwEznTXiJTSkoj4ONmV2SsBl6aUHo6Ik/PHvw9sBVweEUuBv5O9ZawhiIirgb2AdSNiDnAOsAp0vcY3AwcD/wReJ//OgoZmWa9zRLwFaAPWAjoi4pPA1v4SOXiD+F7+IjAWuDj//+iSlFJjZaqtXYN4nY8g+0V9MbAAOKrowkoNwiBeYw2DQbzO7wVOiYglZN/LR/u9PDTuSClJkiSVmO0lkiRJUokZuiVJkqQSM3RLkiRJJWboliRJkkrM0C1JkiSVmKFbkoZZRMyOiNmVrqOUIqIlIlJE7FXpWoYiIlaNiMcj4jclfI4zImJxRGxZqueQVHsM3ZLqRj4EDrgOaj4Qp4gYV6ayVF0+AbydbJ3yLhGxVn6HvTkR8WJE/Coi3tbXCSLixHyontzPc1wMPA9cOJyFS6pthm5J0ogQEW8CPgf8LqU0vcfhy4CTgNuBa4F9gd9HxOo9zvFW4BvA11NKM/p6npTSAiAHvCsidhnOr0FS7TJ0S5JGig8ADWQBu0tEbAC8B/hySun4lNKpwCnAZsAhPc7xf8Ac4NxlPNcVwFLg1BWuWlJdMHRLGtEiYp2IeD0inoj8fuh93OfX+ZaUqUVjEREfj4iHI2JhRDwTEd+NiLX7Ocfx+XMcHxEHRsQfI2JecTtMRKwdEV+LiMfy53w5In4bEe8coP79860Qz0fEooh4OiJuLH5Mvo/54xFxc0Q8lb/fSxFxe0Qc1M95Z+c/1oqIb+b/vDgiWvq473ER8WBELMjXcWlEvKWf806IiMvzr9cbEfFs/vaEPu7b1TceEe+NiPvzf1cvRcQ1+VnnofgI8AZwQ4/xzfKf7y8au7/HMSLig8CBwAkppTcGeqKU0rPA3cB7I2KtIdYpqQ4ZuiWNaCmll4FrgP8BeoXbiNiYLGhN79GS8C3gO8A6wCX5cxxI1p6w6gBP+V7g18CrwPeBn+WfpwG4FzgLmJc//y+AnYHbIuJjfdT2JeC3wF75z/8P+D2wFXBs0V3fTNbusCbwO+CbwE3AFODmiPhoP7WuCvwBOBy4LX+OJ3vcpzn/dTyUr/kx4MPAvRGxXo96dwDa8rX9lazn+c/AMUBbRDT2U8epZDPHs4GLgFnAUcDtETG6n8d0k/9lqBF4IKX0eo/D/85/nlo01lnLU/nHb5D/+lpTSsXhfCB/AkYDewzy/pLq2MqVLkCShltfs7FFGvoYu5gsKH6MLJQW+yiwEllbQef5dyG7IO8JYMeU0kv58c8BdwAbkg9rfTgYODildGuP8fOBrckC/MkppZQ/5/lkQfXbEfHblNLs/Pj+ZBcDPgnsnlJ6pvhk+V8WOr0MbJZSmtPjPmuTBcMLIuLKfC9ysQ2BvwN7ppRe6+frOQh4R0rpwaLztgKfBL5ONrtM/l2Ey4G1gGNTSlcW3f8osl9aroiIrVNKHT2e40Bgh5TSzKLHXAW8HziM/C8uy7Az2d9jW88DKaX/RMSNwDn5iycXAseRhfHOVU4uAl6kxwWYy/DX/Oc9yH7RkjSSpZT88MMPP+riA0hD+BjX47F/BRYDbykaWwl4GngFWKNo/Af5c3y4jxr2yh+b3WP8+Pz49X08ZhXgNbLZ7zf3cfzc/GO/WDT2q/zYe1bwNTs9f549eozPzo9P6udxLfnjP+rj2NpAO7AAGJ0f2zV//3v7Od/dPesoeo6v9HH/vfPHLhzk13lS/v6f7ef42mS/8DwLvEQWtifkj70X6AB2J5u5/k7+Pm8AfwS27uec78g/5zWV/rfhhx9+VP7D9hJJdSelFP190P8M9MVk7/6dUDR2MLAxcEVKaX7R+Pb5z3f2cZ67gSUDlNdXa8KWwOrAQyk/a97DH/KfpxSN7UQW6HrOmPcpIraJiMsi4l/53uvO5RX/X/4uffVHLwT+toxT93oNUkrzgBnAGLJWFyi8Zn/oef8e41P6ONZrdprslyHI2nsGY2z+88t9HUwpzUspnZRS2iil9OaU0rtSSo9HxFiyWe6LUkp3k70jcRLZLwSH5s97a0SM6eO0nX+X6w6yRkl1zPYSScpcQxZAT4yIr6esxaGzj/r/ety382LJ//Y8SUppaUS8OMDz/KePsc7zPdfPYzrHG4rGGoCXU++WkF4iYieyULsyWc/3TWSz9x3AZLIWjb56o59PKQ247jl9vAZ5nV/n2j0+D+Vr7NTex1jnLzYrDVBbsc7Xqa9wPJBvA68DZ+eXHDwF+GlK6dsAEfEacBfZyiiX9njsaj2eW9IIZuiWJLK1lSPiMrILA/ePiFlkvcR/SSk91OPu8/KfNwD+VXwgIlYim/18hr71FWI7z9fnih9kvdXF94MsiI6NiNUGEbw/TxYA904p/bFHvWeThe7B1trTBv2Md34t83p8HsrXOJyez38eO+C9ikTEIWRher+U0vyI2I7s4tIHiu7WeXHtNn2covO5nu/jmKQRxvYSSSr4HlnQ/Bh9XEBZpDN07dnHsd0Z+oTGY2SzqZMjoq92ib17PC9kq34E2S8Gy/J24KWegTuvr69hKHo9Pn+B5mSy9pRH8sOdF1ru1c95Oscf6Of4iupskxnU1uz5r+H7wA9TSrd3Duc/F78rMNDMeedzzRhkjZLqmKFbkvJSSo+TtV8cApxMNpt8bR93vSz/+XMR8ebOwXxf79eW43nfAK4E1gC+XHwsv5rGJ8gu8vxp0aHv5D//v77Wq+4xNht4c36mtvg+HwEOGGq9PXwwInr2YbeQtZNcnVJalB/7E9kvF7tFxHt71PFeshU+/gHcs4L19OdhYC5ZL/xgfDP/+VNFY/8ku3iyeMOcdxedv6fO57pjkM8pqY7ZXiJJ3V1Mtl73BsB3Uu81nUkp/SkivgOcBsyKiOvIQvFhZBfq9de3PJCzyGbJP55fz/oOsgvw3ke2vvbHU0pda2SnlG6LiHOBLwCPRMQNZBcXbgDsRjYTfnz+7t8iC9f3RMTPyFo4GvP3u45sdY7ldQvwp/x5n8ufczeyoH9WUb0pIo4jW5Lx2vwSfY8CW5CtA/4q8KHUe7nAYZF//uuBkyJim5RSXyEZgIjYj+yC2kPyF4V2nuO1iLgIaI6IW8lC+IfJXverepxjFNlW8o+llGYN/1ckqdY40y1J3d0EvJD/c1+tJZ2ayEL3PLJ2lPeTbVDzTrLZ0CHJr1qyM3ABWS/w6cCRZKudHJhSuriPx3wReBfZpjqHkM3KHkDW0nF50f1uJZuR/TvZpjIfARaRta38hhXTSrZ5zWSytbm3JHsnYJeUUrde5pTSX4AdyALqzsCngV2Aq8nW4f7LCtayLJ2v4Yf6u0NErEG2JOQVKaW+XpuzyTYJmkrWgvQXsr+fhT3u906yFWG+v6JFS6oPsewL0yVp5IiI/yGbwfxTSmn3Stej4RURvwUmAeMHs/LLCjzPL8j63d9WPFsuaeRypluSuvsU2QVz3610ISqJT5G17ZxaqieIiMnAe4AWA7ekTvZ0SxrxImJTsqXhJpD16D4E/LyiRakkUkozI+IEsj75UtmQrNfe1hJJXWwvkTTiRcReZBcuvk62esYpKaV/DfQYSZKGwtAtSZIklZg93ZIkSVKJGbolSZKkEjN0S5IkSSVm6JYkSZJKzNAtSZIklZihW5IkSSqx/w/DhI9JDrJJqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAHwCAYAAABzBnP9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABYVUlEQVR4nO3dd5gUVfq38fshIyKIAgZQzAgqqCNiThjAHEDRXbOIOaHiuuqublB/pl0DJkyvac05IeqaAxgQBBURAUVAUQRB4nn/qNYFnIEBZqYm3J/rmqu7T1dXfacs4eH0qXMipYQkSZKkslMr7wCSJElSdWORLUmSJJUxi2xJkiSpjFlkS5IkSWXMIluSJEkqYxbZkiRJUhmzyJYkSZLKmEW2JKlEEfFKRBybdw5JqmossiUpRxExOiImRESj+dqOjYhXcoxVahHRLiKeiIgpETE1Il6KiM6l/OyREfFxREyPiG8j4oaIaFLemSWpIlhkS1L+6gCnLetOIlNhf65HxDrAG8DHwFrAasBjwICI6LSYz54FXAacDTQBOgNtgBciom75pZakimGRLUn5+z+gT0Q0Le7NiNg6It4r9Ba/FxFbz/feKxHx94h4A5gOrB0RKSJOjIjPC73Ll0TEOhHxVkT8FBEPRES9wudXjIinImJSRPxQeN6qlLn/AryVUjo/pTQ5pTQ1pfRv4G6yArpYEbEC8FfglJTScyml2Sml0UAPsmL90FIeX5IqLYtsScrfIOAVoM/Cb0REM+Bp4N/ASsBVwNMRsdJ8m/0R6AU0Br4qtO0BbE7WQ3wOcDNwGNAa2AjoWdiuFnA7sCawBjADuK6UuXcFHiym/QFgu4hoUMLntgYaAI/M35hSmgY8C+xWyuNLUqVlkS1JlcOFwCkR0Xyh9j2Bz1NK/y+lNCeldB8wAth7vm3uSCkNK7w/u9B2WUrpp5TSMGAo8EJKaVRKaQpZIbspQErp+5TSwyml6SmlqcDfgR1KmXllYHwx7eOB2kCzRXzuu5TSnBI+u/A5kKQqxyJbkiqBlNJQ4Cmg70Jvrcb/eqd/9RWw+nyvxxazywnzPZ9RzOvlASJiuYi4KSK+ioifgFeBphFRuxSxvwNWLaZ9VSAB30XEdhExrfAzbL7PrRwRdUr47KRSHFuSKjWLbEmqPC4CjmPBAvobsqEc81sD+Hq+12kZjnkWsAGwZUppBWD7QnuU4rMvAt2Lae8BvJ1SmpVSei2ltHzhp33h/beAmcAB83+oMMNKV+C/S/F7SFKlYpEtSZVESmkk8B/g1PmanwHWj4hDI6JORBwMtCPr9S4Ljcl6tn8sjP++aAk++1dg68KNl80ionFEnAIctaj9FIas/BW4NiL2iIi6EdGGbHz3d8A9S/m7SFKlYZEtSZXLxcBvc2anlL4H9iLrcf6e7CbGvVJK35XR8a4BGpIVt28Dz5X2gymlz4FtgQ7AaOBH4BJg/5TSgMV89nLgT8AVwFTgS2A5oEtK6ecl/B0kqdKJlJblW0ZJkjKFqf/eBi5KKfVfws8eTda7vU1KaUx55JOkilTcTSeSJC2xlNK4iOgK7BsRyxem5CvtZ2+LiNlk0/tZZEuq8uzJliRJksqYY7IlSZKkMmaRLUmSJJWxajkme+WVV05t2rTJO4YkSZKqscGDB3+XUip2ldpqWWS3adOGQYMG5R1DkiRJ1VhELLwi728cLiJJkiSVMYtsSZIkqYxZZEuSJEllzCJbkiRJKmMW2ZIkSVIZs8iWJEmSyphFtiRJklTGLLIlSZKkMmaRLUmSJJUxi2xJkiSpjFlkS5Ik6TdHH300LVq0YKONNvqt7eCDD6Zjx4507NiRNm3a0LFjx/wCVhF18g4gSZKkyuPII4/k5JNP5vDDD/+t7T//+c9vz8866yyaNGmSR7QqxSJbkiRJv9l+++0ZPXp0se+llHjggQd46aWXKjZUFeRwEUmSJJXKa6+9RsuWLVlvvfXyjlLpWWRLkiTVYGPGjOHUU0+lqKiIvffem2effbbEbe+77z569uxZgemqLotsSZKkGmrMmDFstdVWNGzYkOuvv54ePXpw4okncu+99/5u2zlz5vDII49w8MEH55C06nFMtiRJUg11xRVXcNhhh3HZZZcBsOWWW7L55puz3Xbbscoqqyyw7Ysvvkjbtm1p1apVHlGrHHuyJUmSaqi33nqLAw44YIG2Sy65hClTpvDZZ5/RqlUr+vfvD8D9999fuYaKzJsHN94Ip52Wd5Ji2ZMtSZJUQ6266qp89tlndO7c+be2W265heeee45PP/2UFi1a/NZ+xx135JBwEa65Bs46C7p0gZkzoX79vBMtwJ5sSZKkGurEE0/kwgsv5JNPPgFg2rRpnHrqqeyxxx4LFNiVxpw5MG5c9vzYY+HOO+GFFypdgQ0W2ZIkSTXWHnvswfnnn89OO+3EhhtuyBprrMGMGTO4+eab8472ex99BJ07w+67w+zZsMIKcPjhEJF3smI5XESSJKkGO+644zj88MP5/PPPadGiReXrwZ45E/72N7j0UmjWDK67DupU/hK28ieUJElSuapfvz4bbbRR3jF+b8wY2GMPGD4867W+6ipYaaW8U5WKw0UkSZJUOa26Kqy/PjzzTDb+uooU2GCRLUmSpMrkxRdhm23ghx+gbl147DHo2jXvVEvMIluSJEn5+/FHOOYY2HVXmDQJvvkm70TLJLciOyJaR8TLETE8IoZFxO9mEo/MvyNiZEQMiYjN8sgqSZKkcvTYY9CuXTYkpG/fbCaR9u3zTrVM8rzxcQ5wVkrp/YhoDAyOiAEppU/m26YrsF7hZ0ugX+FRkiRJ1UFKcOut0LIlPPUUbFY9+lRzK7JTSuOB8YXnUyNiOLA6MH+RvS9wV0opAW9HRNOIWLXwWUmSJFVFKcHdd2djr9deG+66Cxo3zsZgVxOVYkx2RLQBNgXeWeit1YGx870eV2iTJElSVfTVV9CtWzYl3/XXZ23NmlWrAhsqQZEdEcsDDwOnp5R+WvjtYj6SSthPr4gYFBGDJk2aVNYxJUmStCzmzcuK6o02gtdeg2uvhf/7v7xTlZtci+yIqEtWYN+TUnqkmE3GAa3ne90KKPZW05TSzSmlopRSUfPmzcs+rCRJkpbe1VfDySfDVlvB0KHZ81q59/eWm9zGZEdEAP2B4Smlq0rY7Ang5Ii4n+yGxymOx5YkSaoi5syB8eOhdWvo1QtWWQUOPRSiuMEK1Uues4tsA/wR+DgiPiy0/QlYAyCldCPwDNANGAlMB46q+JiSJElaYh9+mM17/csv2fPGjeGww/JOVWHynF3kdYofcz3/Ngk4qWISSZIkaZn98gtccglcdhmsvHI2Drua3dRYGnn2ZEuSJKk6GTMGdt8dRoyAI4+EK6/MZg6pgSyyJUmStGxSysZZr7YabLghXHNNVmzXYNX3lk5JkiSVv+efz2YMmTwZ6tSBRx6p8QU2WGRLkiRpaUyenA0J2WMPmDIFvv0270SVikW2JEmSlszDD0O7dtnS6OefDx98kL3WbxyTLUmSpNJLCe68Mxt//dxz0LFj3okqJYtsSZIkLVpKcNddsO22sM46WZHduHE2BlvFcriIJEmSSjZ6dDbu+sgjoV+/rG3FFS2wF8MiW5IkSb83bx5cey1stBG8+Wa2qMzll+edqsqwyJYkSdLvXXUVnHoqbLcdDB0KJ54ItSwdS8t+fkmSJGVmz4ZvvoE114Tjj4fVV4dDDskWmtES8Z8jkiRJgvffhy22gK5ds2K7cWPo2dMCeylZZEuSJNVkM2ZA377QqRNMnAj/+AfUrZt3qirP4SKSJEk11Vdfwa67wuefwzHHwBVXQNOmeaeqFiyyJUmSapqUsmEgq68OHTpkU/PtskveqaoVh4tIkiTVJM8+C1tuCZMnZ3NdP/igBXY5sMiWJEmqCb7/Hg4/HLp1g59/hgkT8k5UrVlkS5IkVWcpwQMPwIYbwn33wQUXZDOJbLhh3smqNcdkS5IkVXd33w1rrAEDBmRjsFXuLLIlSZKqm5Tgtttgxx1hnXXgrrtg+eWzMdiqEA4XkSRJqk5GjYIuXeDYY+HGG7O2pk0tsCuYRbYkSVJ1MHcuXHMNbLwxvPdeVmBfdlneqWos/0kjSZJUHVx9NZx9Nuy5Z1Zgt2qVd6IazSJbkiSpqpo1C775Btq0geOPz25u7N49W2hGuXK4iCRJUlX03ntQVJTNez17NjRuDD16WGBXEhbZkiRJVcn06XDOOdC5c7bAzGWXQd26eafSQhwuIkmSVFV89VU2c8jIkXDccfB//wdNmuSdSsWwyJYkSarsUsqGgay+Omy6Kdx0E+y8c96ptAgOF5EkSarMnnoqG3v9/ffZXNcPPGCBXQVYZEuSJFVGkybBoYfC3ntns4hMmpR3Ii0Bi2xJkqTKJCW47z5o1w4eegj+8hcYPBjats07mZaAY7IlSZIqm//8B9ZeG/r3h402yjuNloJFtiRJUt5SgltvhZ12gnXXhTvvhOWXh9q1806mpeRwEUmSpDx98QXssgv06gU335y1NWligV3FWWRLkiTlYe5cuPJK2HjjbMz1LbdkC8uoWrDIliRJysNVV0GfPtniMp98Asce65Lo1YhFtiRJqtGOPvpoWrRowUYL3WB47bXXssEGG9C+fXvOOeecsjnYzJnw5ZfZ8xNOyGYPefzxbJEZVSsW2ZIkqUY78sgjee655xZoe/nll3n88ccZMmQIw4YNo0+fPst+oHfegc03h27dYPbs7MbGAw+097qassiWJEk12vbbb0+zZs0WaOvXrx99+/alfv36ALRo0WLpD/Dzz3DmmbDVVjBlClxxBdStuyyRVQVYZEuSJC3ks88+47XXXmPLLbdkhx124L333lu6HY0eDZtsAldfDb17w7BhsOeeZZpVlZPzZEuSpBpn6tSpPPDAA4wZM4aioiLatWu3wPtz5szhhx9+4O233+a9996jR48ejBo1iijt0I6UsmEgrVpBp05w++2w/fbl8Juossq1JzsibouIiRExtIT3d4yIKRHxYeHnworOKEmSqpdPPvmEDTfckKeffpqUEhdffDGHHHII8+bN+22bVq1accABBxARdOrUiVq1avHdd9+V7gBPPJGNvf7+e6hTJ1si3QK7xsl7uMgdwB6L2ea1lFLHws/FFZBJkiRVY8cffzwXXHABjzzyCBdffDHvvPMOq6yyygJF9H777cdLL70EZENHZs2axcorr7zoHU+cCIccAvvuC3PmQGmLclVLuRbZKaVXgcl5ZpAkSTXHxIkTGTp0KMccc8xvbYcddhhvvfUWkyZNolWrVvTv35+jjz6aUaNGsdFGG3HIIYdw5513ljxUJCW45x5o1w4efRQuuQQGDYINNqig30qVUVUYk71VRHwEfAP0SSkNyzuQJEmqmiKClBIppd/a7rvvPt5++22OO+44Pv7449/a77777tLv+MEHYb31oH//rNhWjZf3cJHFeR9YM6XUAbgWeKykDSOiV0QMiohBkyZNqqh8kiSpCmnevDmbbrop/fr1+61t7ty5XH755fTo0aP0O5o3D/r1g88/z25wvOsueP11C2z9Jub/l1wuASLaAE+llDYqxbajgaKU0iIHORUVFaVBgwaVTUBJklStfP755+y66660bduWjTfemGeffZZVVlmFJ598koYNG5ZmB9kS6K++CueeC5deWv6hVSlFxOCUUlFx71XqnuyIWCUKA6AiohNZ3u/zTSVJkqqy9dZbj+HDh3PEEUew8sor869//YsBAwYsvsCeMwcuvzyb93rIELjtNvjnPysmtKqcXMdkR8R9wI7AyhExDrgIqAuQUroROAg4ISLmADOAQ1LeXe+SJKnKa9iwIT179lyyD119ddZzvf/+cP31sOqq5RNO1ULuw0XKg8NFJElSmZg5E77+GtZeG6ZNgxdfzKboK+2iNKrWquxwEUmSpNy89RZsuil06wazZ8Pyy8N++1lgq1QssiVJkuY3bRqcfjpss032/OqroW7dvFOpiqkK82RLkiRVjNGjYaedsseTTspubGzcOO9UqoIssiVJklLKhoG0bg1bb53Ne73ddnmnUhXmcBFJklSzPfoodOwI330HtWtnS6RbYGsZWWRLkqSa6dtvoXt3OOAAqFULJk/OO5GqEYtsSZJUs6SUDQdp1w6efBL+8Q94911Yf/28k6kacUy2JEmqeR57LCuyb70V2rbNO42qIYtsSZJU/c2bB/36wa67Zj3Wd94JjRplw0SkcuCVJUmSqrdPP4UddoCTT4bbb8/aGje2wFa58uqSJEnV0+zZcOml0KEDDBsGd9yRjb+WKoBFtiRJqp6uvhrOOw/23hs++QSOOMIl0VVhHJMtSZKqj19+gXHjYN114cQTs5sa99kn71SqgezJliRJ1cPrr2dDQ/baC+bMgeWXt8BWbiyyJUlS1TZ1anZT43bbwaxZcO21UMcv65Uvr0BJklR1ffkl7LgjjB0Lp54Kf/971oMt5cwiW5IkVT0pZTcxrrFGNj1f796w9dZ5p5J+43ARSZJUtTz0EGyyCXz3HdSunS2RboGtSsYiW5IkVQ3jx8OBB0L37lCvHvzwQ96JpBJZZEuSpMotpWwhmXbt4OmnswVm3nkH1lsv72RSiRyTLUmSKrcIePJJ2HhjuPVWWH/9vBNJi2WRLUmSKp+5c+H662H33WGDDbKe7EaNoJZfwqtq8EqVJEmVy/Dh2ZzXp52W3dQI0LixBbaqFK9WSZJUOcyenc1z3bEjfPYZ3H03/O1veaeSlopFtiRJqhyuvhr+/GfYf3/45BM47LBsPLZUBTkmW5Ik5WfGDBg3Lpsp5KSToH172HPPvFNJy8yebEmSlI9XX4UOHWCvvWDOnOzGRgtsVRMW2ZIkqWL99FPWa73DDllxfcMNUMcv11W9eEVLkqSKM2oU7LhjNkTk9NOzGxsbNco7lVTmLLIlSVKZmjJlCm+++SYrrLACW221FbVq1YJ587Ip+NZcE3beGXr3hs6d844qlRuLbEmSVGZuvPFGzjvvPDbbbDMmTZrEzF9+YWDv3rTq3x9eeQWaN88WlpGqOYtsSZJUJt5++20uueQSBg0axDrrrEP6+mvG7LUXrc46i1RUREyZkhXZUg3gjY+SJKlM3HHHHZx++umss/bacOutRPv2rDFiBFetuiqvX345rLtu3hGlCmORLUmSysSUKVNo2bJltoDMc89Bx47Exx/z/MYb8+O0aXnHkyqURbYkSVp2c+dyekq81K8f8+bNy8Zdv/QSX9Wty7vvvsu2226bd0KpQjkmW5IkLZthw+CYY9jynXf4bI012G233TjyyCOZOHEi11xzDRdffDErrrhi3imlCmWRLUmSls6sWXDppdlc102awL330n2//Zh177088cQTNGnShHvvvddebNVIkVLKO0OZKyoqSoMGDco7hiRJ1dtll0HfvtCzJ/zrX84cohonIganlIqKe8+ebEmSVHrTp2erNa6/Ppx8MmyyCXTtmncqqdLxxkdJklQ6r7ySFdV77w1z5mTLoVtgS8WyyJYkSYs2ZQocfzzstFP2+sYboY5fhkuLkuv/IRFxG7AXMDGltFEx7wfwL6AbMB04MqX0fsWmlCSpavnhhx94+umnmTNnDt26daNFixZLv7NRo2D77WH8eDjrLLj4YlhuubILK1VTefdk3wHssYj3uwLrFX56Af0qIJMkSVXWI488wjrrrMMjjzzCs88+ywYbbED//v2XfEfz5mWPa64Ju+0Gb70FV1xhgS2VUq492SmlVyOizSI22Re4K2VToLwdEU0jYtWU0viKSShJUtUxadIkjj32WAYOHMimm24KwMiRI9lqq63YYYcdWLc0y5qnBPffn/VY//e/0KIF3HZbOSeXqp+8e7IXZ3Vg7HyvxxXaJEnSQh599FH22GOP3wpsgHXXXZdDDz2U//znP4vfwbhxsM8+cOih2bzXU6eWY1qpeqvsRXYU01bsxN4R0SsiBkXEoEmTJpVzLEmSKp+ZM2fSqFGj37U3atSImTNnlvzBlODmm6F9exg4EK66Ct54A9ZZpxzTStVbZS+yxwGt53vdCvimuA1TSjenlIpSSkXNnQxfklQD7bnnnjz66KOMH/+/UZU//vgjd999N/vss0/JH4yAF1+EoiIYOhTOOANq166AxFL1Vdnn33kCODki7ge2BKY4HluSpOKtvfba9OnTh6KiIo455hjq1q3LbbfdRo8ePSgqWmhRujlzslUa99oLNtgAbr89u6kxivsSWdKSynsKv/uAHYGVI2IccBFQFyCldCPwDNn0fSPJpvA7Kp+kkiRVDX379mW33XbjgQceYPr06dx9991svfXWC2708cdwzDHw3nvZHNgXX5wtLCOpzEQ2cUf1UlRUlAYNGpR3DEmSKpeZM+Ef/8h+VlwRrrsOune391paShExOKVUVNx7lX1MtiRJKivXXJP1Wh9yCAwfDj16WGBL5aSyj8mWJEnL4uefs6n5NtgATjkFNt00W1xGUrmyJ1uSpOrqpZdgk01g772zGx2XW84CW6ogFtmSJFU3P/4Ixx0Hu+ySTcV3661Qxy+vpYrk/3GSJFUnX3wB228P334L55wDf/kLNGyYdyqpxrHIliSpOpg7N+u1btMGunaF3r2zxWUk5cLhIpIkVWUpwT33QLt2MHHi/4aHWGBLubLIliSpqho7Nlux8Q9/gGbNYNq0vBNJKrDIliSpqkkJ+vXLeq9feSVbHv3112HttfNOJqnAMdmSJFU1EVlxvdVWcNNNsNZaeSeStBCLbEmSKkBKiUGDBvH111+z+eab07p16yXbwZw5cNVV2ZzXG24It9+ezRriio1SpeRwEUmSytm3337LVlttxWGHHcZtt91Ghw4dOOOMM0gplW4HH30EW24J554L99+ftS23nAW2VIlZZEuSVM6OOeYYdt55Z0aMGMETTzzBqFGjeOONN7jjjjsW/cFffoE//zmbKeTrr+Ghh+Cvf62QzJKWjUW2JEnlaMKECbz11ltceOGF1KqV/bXbtGlTLrroIm6//fZFf/iaa+Dvf4fDDoNPPoEDDyz/wJLKhGOyJUkqSCnx448/0qhRI+rVq1cm+5w2bRqNGjWifv36C7Q3b96cn376qbgPwLhx0LYtnHpq1ovdpUuZZJFUcezJliQJePzxx2nXrh1rrrkmLVu25KyzzmLWrFnLvN+11lqLBg0a8OKLLy7Qfvvtt7P77rsvuPGAAbDxxrDPPtmNjsstZ4EtVVEW2ZKkGu+NN96gd+/eXH/99UyZMoWhQ4fy2Wefcfrppy/zvmvVqsV1113HYYcdxp///GfuvvtuevTowcsvv0yfPn2yjX74AY4+GnbbDerXh9tugzp+2SxVZVHqO5urkKKiojRo0KC8Y0iSqoju3bvTpUsXjj/++N/aJk+ezDrrrMOoUaNYccUVl/kYn376Kbfccgtff/01W265JUcffTQrrLACfPEFbLstTJoE55wDF14IDRos8/Eklb+IGJxSKiruPf+ZLEmq8UaNGsVmm222QFuzZs1YddVV+eabb8qkyN5ggw244oor/tcwd272uNZa2fCQ3r1h002X+TiSKgeHi0iSaryOHTsyYMCABdpGjx7NhAkTaNOmTdkeLCW4887sxsYJE6BWrWzVRgtsqVqxJ1uSVOOdffbZbL/99jRq1IgDDjiATz/9lLPOOos+ffrQqFGjsjvQV1/B8cfD88/DNtvA9Ollt29JlYo92ZKkGq9t27a8+OKLvPbaa2y55Zace+65nH766fTt27dsDpASXHcdtG8Pr78O114Lr76aDRWRVC3Zky1JErDJJpvw0EMPlc/OI7Liettts6Eha65ZPseRVGksUZEdESsCrVNKQ8opjyRJ1cPs2XDllbDvvrDhhnD77dmsIRF5J5NUARY7XCQiXomIFSKiGfARcHtEXFX+0SRJqqI++AA6dYLzzoP//Cdra9jQAluqQUozJrtJSukn4ADg9pTS5oDLT0mStLBffskK6y22gPHj4eGH4S9/yTuVpByUpsiuExGrAj2Ap8o5jyRJVdc118Cll8Lhh8Pw4XDAAXknkpST0ozJvhh4Hng9pfReRKwNfF6+sSRJqiKmToWxY6FdOzj11GyYyM47551KUs5cVl2SpKX1/PPQqxfUq5f1XNdx0i6pJlmqZdUj4lqgxAo8pXRqGWSTJKnqmTwZzjgD7rorW7mxf38LbEkLWNSfCHYFS5K0sJEjs9UaJ0+GP/85+6lfP+9UkiqZEovslNKdFRlEkqRKbc6crLd67bWzGxp794YOHfJOJamSKs082c0j4oqIeCYiXvr1pyLCSZKUu5SyhWTatoUJE6BWLejXzwJb0iKVZgq/e4DhwFrAX4HRwHvlmEmSpMrhyy9h993h6KNhtdVgxoy8E0mqIkpTZK+UUuoPzE4p/TeldDTQuZxzSZKUn5Tg3/+GjTaCt96CG26AV16BNm3yTiapiijNrdCzC4/jI2JP4BugVflFkiQpZxHw5puwww5w442wxhp5J5JUxZSmyP5bRDQBzgKuBVYAzijXVJIkVbTZs+Gyy2D//aF9+2wcdoMGWcEtSUtosUV2SunXpdSnADuVbxxJknIweHA27nrIEJg3LyuyGzbMO5WkKmyxRXZE3E4xi9IUxmZLklR1zZgBf/kLXHEFtGwJjz0G++6bdypJ1UBphos8Nd/zBsD+ZOOyJUmq2v71L7j8cjj2WPi//4OmTfNOJKmaKM1wkYfnfx0R9wEvllsiSZLK008/wdix2ZCQ006DrbbKbnCUpDJUmin8FrYe4G3WkqSq55lnsmn59t8/W8GxYUMLbEnlojQrPk6NiJ9+fQSeBM4ti4NHxB4R8WlEjIyIvsW8v2NETImIDws/F5bFcSVJNcx338Ef/wh77gmNG8Ndd2VLpEtSOSnNcJHG5XHgiKgNXA/sCowD3ouIJ1JKnyy06Wsppb3KI4MkqeqZPn06V155JQ8//DARwUEHHcSZZ55Jw5JmAxk5ErbeGn74AS68EP70J6hfv2JDS6pxSiyyI2KzRX0wpfT+Mh67EzAypTSqcLz7gX2BhYtsSZIAmDdvHnvttRcrrrgi/fr1I6XE5Zdfzj777MPzzz9PrVrzfUE7Z07WW7322tC9Oxx/PGyySX7hJdUoi+rJvrLw2AAoAj4CAtgEeAfYdhmPvTowdr7X44Ati9luq4j4iGxGkz4ppWHF7SwiegG9ANZwZS5JqpZeeOEFfvjhB1588cXfCuqHH36YTTfdlJdeeokuXbpkS6Lfeiv885/Zqo2rrALXX59zckk1TYljslNKO6WUdgK+AjZLKRWllDYHNgVGlsGxi1tCa+H5uN8H1kwpdSBbbfKxReS9uZCxqHnz5mUQT5JU2bzzzjt069ZtgR7r2rVr061bN95991344gvYZRfo1QvWXBNmzcoxraSarDSzi7RNKX3864uU0lCgYxkcexzQer7XrVho/u2U0k8ppWmF588AdSNi5TI4tiSpClpjjTUYNuz3X2gO/fhjdhs6FDbeOFu98aabYOBA8JtNSTkpTZE9PCJuLcz0sUNE3AIML4NjvwesFxFrRUQ94BDgifk3iIhVIiIKzzsV8n5fBseWJFVB3bt3Z9CgQdx8883MmTOHOXPm0K9fPz4aMoSOc+ZkvdjDhmU92bWWZpZaSSobpZm/6CjgBOC0wutXgX7LeuCU0pyIOBl4HqgN3JZSGhYRvQvv3wgcBJwQEXOAGcAhKaXfLfEuSaoZll9+eZ5//nmOO+44Luzbl9NmzOCTDTfkhRdeoM5aa2WzhkRxoxElqWJFdaxZi4qK0qBBg/KOIUkqL++9x+wjjqDu8OFw8cVwwQV5J5JUA0XE4JRSUXHvLWoKvwdSSj0i4mN+f0MiKSXnQZIkVazp0+Gii+Cqq6i76qrwxBOw9955p5Kk31nUcJFfh4e4EIwkqXL497/hiiuyOa8vuwyaNMk7kSQVq8QiO6U0vvD0O2BGSmleRKwPtAWerYhwkiQxZQqMGZPNHHLaabDNNrDddnmnkqRFKs2t168CDSJidWAg2Y2Qd5RnKEmSAHjySWjXDg44IFvBsWFDC2xJVUJpiuxIKU0HDgCuTSntD7Qr31iSpBpt0iQ49FDYZx9YaSW4995siXRJqiJK8ydWRMRWwGHAMUvwOUmSltznn8NWW8FPP2Uzh5x7LtSrl3cqSVoipSmWTwfOAx4tzGO9NvByuaaSJNU8s2dD3bqwzjpZL/bxx0P79nmnkqSlstjhIiml/6aU9gGuK7welVI6tdyTSZKqjalTp/Lhhx/y3Xff/f7NefOyZdDXWw++/TZbqfHf/7bAllSlLbbIjoitIuITCkupR0SHiLih3JNJkqq8lBKXXHIJa665Jocffjjrrbcexx13HDNnzsw2GDkyWwq9d++sB3vWrHwDS1IZKc2Nj9cAuwPfA6SUPgK2L8dMkqRqon///jzyyCMMGTKEIUOG8OWXXzJx4kT6nntuNt/1xhvDBx/ArbfCiy/CGmvkHVmSykSpbmBMKY2NiPmb5pZPHElSdXLDDTdw5ZVX0qpVKwCaNm3KDTfcQPv27blijz2ovfvucMMNsNpqOSeVpLJVmiJ7bERsDaSIqAecSmHoiCRJi/Ltt9+y3nrrZS9mzoR//IPVDjqIefPm8dM117Biy5awYCeOJFULpSmyewP/AlYHxgEvACeWZyhJUvWwzTbb8PDDD3N6585wzDHwySeMHDOG1q1b09QCW1I1ttgiO6X0Hdkc2QBExIpkRfbfyzGXJKka+MvZZ/PajjuSZsxgVsuWvHDCCRz3yCPceuuthAW2pGqsxCI7IloDFwCrAY8C9wEXA4cXnkuStEjtBw6k/YwZvNKuHRc3bEjLH37gscceo3PnznlHk6Rytaie7LuA/wIPA3sAbwPDgE1SSt9WQDZJUg6mTZvGFVdcwUMPPURE0L17d/r06cNyyy1Xuh38+COMGQObbAJnnAE77MCOW2/NjuUZWpIqmUUV2c1SSn8pPH8+IiYAW6SUZpZ/LElSHubOncuee+7JKquswu233868efO4/PLL2XvvvRkwYAC1ai1m5tfHH4cTToBGjWDECGjQALbeumLCS1Ilssg/LSNixYhoFhHNgG+B5eZ7LUmqZp577jmmTZvGfffdxxZbbMGWW27JAw88wKRJkxg4cGDJH5w4EQ4+GPbbD1q0gPvvh9q1Kyy3JFU2i+rJbgIMBua/M+X9wmMC1i6vUJKkfAwaNIiuXbsu0GNdu3ZtunbtyqBBg9h1111//6HPP4fOnWHaNPjb3+Ccc6Bu3QpMLUmVT4lFdkqpTQXmkCRVAmussQaPP/7479o//vhjDjvssAUbZ82CevWy5dD/+Efo1QvataugpJJUuZVmWXVJUg3RvXt3Bg0aRL9+/Zg9ezazZ8/m2muvZdiwYRxwwAHZRvPmZas0rrsujB8PtWrBNddYYEvSfCyyJUm/WX755RkwYAD33XcfLVq0oHnz5jzyyCMMGDCAhg0bwmefwY47wkknQdu2MGdO3pElqVIqzYqPkqQaZMMNN+TVV19l4sSJRATNmzeHlOCyy+Cii6BhQ7j9djjiCFdslKQSLGoxmkXOIJJSmlz2cSRJlUWLFi3+9yICPv4Y9twTrrsOVl01v2CSVAUsqid7MNksIgGsAfxQeN4UGAOsVd7hJEk5mjkzmy2kRw/YeGO47bbsRkdJ0mItanaRtQAi4kbgiZTSM4XXXYEuFRNPkpSLN9+EY47JFpRZbrmsyLbAlqRSK82Nj1v8WmADpJSeBXYov0iSpNxMmwannQbbbgvTp8Nzz8F55+WdSpKqnNIU2d9FxJ8jok1ErBkR5wPfl3cwSVIOrr02+znpJBg6FHbfPe9EklQllabI7gk0Bx4t/DQvtEmSqoMffoCPPsqen3FGNlTk2muhceN8c0lSFbbYKfwKs4icFhHLp5SmVUAmSVJFeeSRrNd6+eWz8dcNGmRLpEuSlslie7IjYuuI+AT4pPC6Q0TcUO7JJEnl59tv4aCD4MADs+n4HnwQatfOO5UkVRulWYzmamB34AmAlNJHEbF9uaaSJJWfzz7LequnT4d//hPOOgvq1s07lSRVK6Va8TGlNDYWXNVrbvnEkSSVm1mzsmn41l0XjjwSjj8eNtgg71SSVC2V5sbHsRGxNZAiol5E9AGGl3MuSVJZmTcvW6VxnXXgm2+gVi246ioLbEkqR6UpsnsDJwGrA+OAjoXXkqTKbsQI2H57OOUUaN8eUso7kSTVCKUZLtIwpXTY/A0RsUo55ZEklYWU4NJL4a9/zVZsvPNO+OMfYcGhf5KkclKanuwvI+K+iGg4X9szJW4tScpfBHzyCey9NwwfDocfboEtSRWoNEX2x8BrwOsRsU6hzT+pJamymTED/vQnGDIke92/fzY1X8uW+eaSpBqoNMNFUkrphoj4CHgyIs4FHNQnSZXJ66/DMcdk0/M1aQKbbJLNJCJJykVperIDIKX0BrALcDbQtjxDSZJKaepUOPlk2G67bIq+AQPg3HPzTiVJNV5piuxuvz5JKY0Hdgb2KLdEkqTSu/56uOEGOO00+Phj6NIl70SSJBYxXCQi/pBSuhvoGcXfLPPqsh48IvYA/gXUBm5NKV260PtReL8bMB04MqX0/rIeV5KqtO+/h7FjoWNHOP102Hln6NQJgHnz5hERlPDntiSpgiyqJ7tR4bFxCT/LJCJqA9cDXYF2ZMV8u4U26wqsV/jpBfRb1uNKUpWVEjz0ELRrB927w9y50KABdOrEsGHD6NatG/Xq1aNp06aceuqpTJs2Le/EklRjlVhkp5RuKjz+tbifMjh2J2BkSmlUSmkWcD+w70Lb7AvclTJvA00jYtUyOLYk5S6lxIABAzjzzDO58MIL+fTTT0veePx4OPDArLhu3TortmvXBuDbb7+lS5cudOvWjZ9++olhw4YxefJkevToUUG/iSRpYYsaLvLvRX0wpXTqMh57dWDsfK/HAVuWYpvVgfEL7ywiepH1drPGGmssYzRJKl/z5s3j8MMP54MPPuAPf/gDkydPZrvttuOKK67g8MMPX3Djzz6DLbeEX36Byy+HM86AOv/74/vWW29l33335eSTTwZgueWW44477mDttddmyJAhbLLJJhX5q0mSWPQUfoMLj9uQDef4T+F19/neWxbFDRhceGrA0myTNaZ0M3AzQFFRkVMMSqrUnnzySYYNG8bgwYNp0KABAEcffTRbb701++67L02aNIGZM6F+fVh33Wx6vl69YP31f7evESNGsNtuuy3QVqdOHTp16sSIESMssiUpB4saLnJnSulOsvHQO6WUrk0pXUs2jV/HMjj2OKD1fK9bAd8sxTaSVOU88cQTHHvssb8V2AAbbrghnTt35qUBA+Caa2DtteHrr6FWLbjiimILbID27dvz6qsL3os+a9Ys3nrrLdq3b1+ev4YkqQSlmcJvNRa80XH5Qtuyeg9YLyLWioh6wCHAEwtt8wRweGQ6A1MK0whKUpVWv359pk+f/rv25pMmseMFF2RDQjp2LNVS6MceeyzPPvssl19+OZMnT+bzzz/n0EMPpVOnThbZkpST0hTZlwIfRMQdEXEH8D7wj2U9cEppDnAy8DwwHHggpTQsInpHRO/CZs8Ao4CRwC3Aict6XEmqDHr27Mn111/PhAkTsoaU+OKoo7hl8GCaTpoEd98NTz0Fqy2+T6N58+a88sorvPfee7Ru3ZrtttuOddZZh3vvvbecfwtJUkkipZKHL0dELaAzWaH7602J76SUvq2AbEutqKgoDRo0KO8YkrRIf/vb37j66qvZc889+f777zl84EB22mYbWtx3H7RokXc8SdJiRMTglFJRse8tqsgufPitlNJW5ZKsnFhkS6oSZsxgyllnMbBFC2a1bcteu+/O8iuumHcqSVIpLarILs1wkRci4sBw+TBJKjv//S906ECTfv04YLnlOOSQQyywJakaWdQUfr86k2z1xzkR8QvZtHoppbRCuSaTpOrop5/g3HPhxhuz2UMGDsyWRZckVSuLLbJTSsu8hLokqeCGG+Dmm+HMM+Hii6FRo7wTSZLKQWl6somIFcnmy/5tQteU0qslf0KS9JvvvoMxY2CzzbKp+bp0gaJih/BJkqqJxRbZEXEscBrZQjAfks028hbg95uStCgpwX/+A6ecAk2bwogR2QqOFtiSVO2V5sbH04AtgK9SSjsBmwKTyjWVJFV1X38N++0HPXvCWmvBI49A7dp5p5IkVZDSDBf5JaX0S0QQEfVTSiMiYoNyTyZJVdWnn0KnTjB7Nlx5JZx2mgW2JNUwpSmyx0VEU+AxYEBE/AB8U56hJKlK+uUXaNAA1l8feveGXr1gnXXyTiVJysFih4uklPZPKf2YUvoLcAHQH9ivnHNJUtUxdy5cdVU2Jd/XX0MEXHaZBbYk1WAl9mRHRLNimj8uPC4PTC6XRJJUlQwdCsccA+++C3vtBbVKc6uLJKm6W9RwkcFAIlt8Zg3gh8LzpsAYYK3yDidJlVZK2TzXf/87NGkC994LhxyS9WJLkmq8EovslNJaABFxI/BESumZwuuuQJeKiSdJlcvUqVO5/fbbee211+jz8cesv+uurHjHHdC8ed7RJEmVSGm+19zi1wIbIKX0LLBD+UWSpMrpx2++4eG11mL0449zwAEH8GLPnrQdNIiH/vvfvKNJkiqZ0swu8l1E/Bm4m2z4yB+A78s1lSRVNi+/zLwDD+TIH36Arl2z+a+BXbt2Zf/992ffffelbt26OYeUJFUWpenJ7gk0Bx4lm8avRaFNkqq/KVOyqfh23pmfZ8zg/auugj59fnu7U6dONGnShGHDhuUYUpJU2Sy2JzulNJls1UdJqnn69YP+/eHsszltyBC6r7IKm8339uzZs5k8eTJNmjTJLaIkqfJZbJEdEesDfYA282+fUtq5/GJJUo4mTYIxY2DzzeGMM2C33WCzzTjs4Ye58MIL2WWXXWjRogUpJS699FLatm3LWms54ZIk6X9KMyb7QeBG4FZgbvnGkaQcpQT33QenngorrggjRkD9+rBZ1nd9wAEHMGTIEDbYYAM6d+7MF198QZMmTXj00UdzDi5JqmwipbToDSIGp5Q2r6A8ZaKoqCgNGjQo7xiSqpKxY+GEE+Dpp2HLLbMhIu3bF7vpxIkTeeedd1hllVUoKioinBtbkmqkQp1cVNx7penJfjIiTiS78XHmr42FsdqSVPV9+ilssUW2PPrVV8Mpp0Dt2iVu3qJFC/bee+8KDChJqmpKU2QfUXg8e762BKxd9nEkqQLNmAENG8L668NJJ8Fxx8Ha/tEmSVp2pZldxLt5JFVpc+bM4fXXX2f69Olst912NG7YMOuxvvJKGDQIWrWCf/4z75iSpGqkND3ZRMRGQDugwa9tKaW7yiuUJJWVwYMHc+CBB7LyyivTpEkT/nHIITy28sqs/OWXsM8+ixwWIknS0lrsYjQRcRFwbeFnJ+ByYJ9yziVJy2zmzJnsu+++XHnllQx67z0GbrMNr82YAaNHM/aKK+Cxx2DVVfOOKUmqhkrTk30Q0AH4IKV0VES0JJvOT5IqtRdeeIF1112XAw88MGsYO5bo2ZMbV1qJqRMncpmzgkiSyklpllWfkVKaB8yJiBWAiXjTo6QqYNqECfT55hv48MOs4ZZb4K67WGGttZgyZUqu2SRJ1VtpiuxBEdEUuAUYDLwPvFueoSRpmQ0cSPdLLmGvzz/np8cey9rq1GHOnDncfffd7L777rnGkyRVb6WZXeTEwtMbI+I5YIWU0pDyjSVJS+nHH6FPH+jfnzrrrcf/O+44LrrrLk5bcUWaNGnCrbfeSsuWLZ3nWpJUrhZbZEfEwJTSLgAppdELt0lSpXLjjXDHHdC3L1x4IX9s2JBWL7/M3XffzfTp0+nduzcHH3wwdeqUanIlSZKWSonLqkdEA2A54GVgR+DXO4RWAJ5NKW1YEQGXhsuqSzXMhAkwZky2auPMmTB8OHTsmHcqSVI1t7TLqh8PnA6sRjYW+9ci+yfg+rIMKElLJSX4f/8PTj8dVloJRoyA+vUtsCVJuSvxxseU0r8Kqz32SSmtnVJaq/DTIaV0XQVmlKTfGzMGunWDI46Atm3hiSdcWEaSVGmU2JMdEVsAY1NK1xZeHw4cCHwF/CWlNLliIkrSQkaMyIaGpAT//jeceKIFtiSpUlnUFH43AbMAImJ74FLgLmAKcHP5R5OkhUyfnj1usAGcdhoMHQqnnGKBLUmqdBZVZNeer7f6YODmlNLDKaULgHXLP5okFcyZA5deCmutBePGQQT87W/Qpk3eySRJKtYii+yI+HU4yS7AS/O959xXkirGhx/CllvCeefBtttC3bp5J5IkabEWVWTfB/w3Ih4HZgCvAUTEumRDRiSp/KQEf/5zNvb666/hoYfg4YehZcu8k0mStFgl9kinlP4eEQOBVYEX0v8m1K4FnFIR4STVYBHwzTdw2GFw1VXQrFneiSRJKrVFDvtIKb1dTNtn5RdHUo02bVrWe33kkdlc17fc4k2NkqQqybHVkiqHF16AXr2y+a/XXDMrsi2wJUlV1KLGZJebiGgWEQMi4vPC44olbDc6Ij6OiA8jwnXSpepo8mQ46ijYfXdo0ABeew3OOCPvVJIkLZNcimygLzAwpbQeMLDwuiQ7pZQ6lrQuvKQq7uabs6XRzz8fPvyQwQ0asOeee9KkSRPWX399rrnmGubNm5d3SkmSlkheRfa+wJ2F53cC++WUQ1Ievv0W3n03e37mmfDBB/C3vzH8yy/p2rUr++67L6NGjeLuu+/mP//5D+eff36+eSVJWkLxv0lDKvCgET+mlJrO9/qHlNLvhoxExJfAD0ACbkoplbjSZET0AnoBrLHGGpt/9dVXZZ5b0jJKCe68MxsOsvLK2fLo8427Pu6442jTps0CRfWECRNo27YtX375JU2bNs0htCRJxYuIwSWNtii3nuyIeDEihhbzs+8S7GablNJmQFfgpMLy7sVKKd2cUipKKRU1b958mfNLKmOjR2fjro86CjbaCJ566nc3Ng4dOpQdd9xxgbaWLVuy5pprMmrUqIrLKknSMiq32UVSSl1Kei8iJkTEqiml8RGxKjCxhH18U3icGBGPAp2AV8slsKTyM2IEFBVlc19ffz307g21fv9v/PXWW493332XbbbZ5re2H374ga+++oo111yzIhNLkrRM8hqT/QRwROH5EcDjC28QEY0iovGvz4HdgKEVllDSsvv55+xxgw2ysddDh8KJJxZbYAOcfvrp/POf/+TRRx9l7ty5jBo1ip49e9KzZ09WWmmlCgwuSdKyyavIvhTYNSI+B3YtvCYiVouIZwrbtARej4iPgHeBp1NKz+WSVtKSmT0b/vEPWGstGDs268G++OJs/utF2Gyzzbj//vv55z//SYMGDdhiiy3YfPPN+de//lVBwSVJKhu5LEaTUvoe2KWY9m+AboXno4AOFRxN0rJ6/304+mj46CPo3h3q11+ij++88868++67zJo1i7p16xIR5RRUkqTyk1dPtqTqJiU47zzo1AkmTIBHHoEHHoAWLZZqd/Xq1bPAliRVWRbZkspGBEycCEceCZ98Avvvn3ciSZJyk8twEUnVxNSp8Kc/ZdPybbZZtnrjQtPySZJUE9mTLWnpPPsstG+fTcn32mtZmwW2JEmARbakJfX993D44dCtGyy/PLzxBpx2Wt6pJEmqVCyyJS2ZW26B++6DCy6ADz6ArbbKO5EkSZWOY7IlLd748TBmDGy5ZbaozF57ZUujS5KkYtmTLalkKcFtt8GGG8If/whz50K9ehbYkiQthkW2pOKNGgW77grHHAMdOsDTT3tjoyRJpeRwEUm/N3w4FBVlRXW/ftCrF9Ty3+SSJJWWRbak/5k2LZsxpG1bOOecbHn01q3zTiVJUpVj15QkmDULLrkE2rTJbnCMgIsussCWJGkp2ZMt1XSDBmXjrocMgUMOgYYN804kSVKVZ0+2VFOlBOeem03L99138Pjj2fzXzZvnnUySpCrPIluqqSJg8uSsF/uTT2CfffJOJElSteFwEakm+ekn6Ns3K6w33xxuuslZQyRJKgf+7SrVFE8/De3bZ4X1G29kbRbYkiSVC/+Glaq7SZPgsMOypdCbNIE334RTT807lSRJ1ZpFtlTd3XYbPPhgNiXf++9nNzpKkqRy5ZhsqTr6+utsvuuttoIzzoC994Z27fJOJUlSjWFPtlSdpAS33JIV1IcfDnPnQr16FtiSJFUwi2ypuvjiC9hlF+jVK5s55LnnoHbtvFNJklQjOVxEqg6GD88K67p14eab4dhjs3mwJUlSLiyypaps6lRo3Bjatv3f/Nerr553KkmSajyHi0hV0axZ8Ne/Qps28NVXWa/1hRdaYEuSVEnYky1VNe++m/VYDx0Khx4KjRrlnUiSJC3EnmypqkgJ+vTJpuX74Qd48km45x5YeeW8k0mSpIXYky1VFRHw00/Z7CGXXQYrrJB3It58800eeughUkoccMABbLvttoQ3XEqSZE+2VKn9+CMcfzwMGpS9vvFG6NevUhTYF110ET179qRZs2Y0b96cI444gr59++YdS5KkSsGebKmyeuIJOOEE+PZb2HhjKCqCWpXj38UjRozgpptuYujQoaxcGK5ywgknsPHGG3PooYfSoUOHnBNKkpSvyvE3tqT/mTgRDjkE9t0XVloJ3nkHTj4571QLeOqpp+jevftvBTbAiiuuSM+ePXnyySdzTCZJUuVgkS1VNrffDo8+Cpdckg0TKSrKO9HvNGzYkKlTp/6uferUqTRs2DCHRJIkVS6RUso7Q5krKipKg34dwypVBWPHZj9bbw2zZ2dLpLdtm3eqEk2YMIENN9yQl156iY4dOwIwbNgwtt9+ez788ENat26db0BJkipARAxOKRXbG+aYbClP8+bBLbfA2WfDKqtky6PXrVupC2yAli1bcsstt7DTTjux/fbbU6tWLV555RWuv/56C2xJkrDIlvLz+edw3HHw3//CLrvAzTdD7dp5pyq1Aw88kF122YVnn32WlBL9+/enWbNmeceSJKlSsMiW8jB8OGy2GdSvD/37w1FHZfNgVzFNmzalZ8+eeceQJKnSsciWKtJPP2VzXLdtC+efD0cfDautlncqSZJUxpxdRKoIM2fCBRdAmzbw1VdZr/Wf/2yBLUlSNWVPtlTe3noLjjkmGyJy+OGw/PJ5J5IkSeXMnmypvKQEZ5wB22wD06bBM8/AnXdmC8xIkqRqzSJbKi8RMGMGnHgiDBsGXbvmnUiSJFWQXIrsiOgeEcMiYl5ElLicXUTsERGfRsTIiOhbkRmlpfLjj9m0fL8uhtSvH1x3HTRunGssSZJUsfLqyR4KHAC8WtIGEVEbuB7oCrQDekZEu4qJJy2Fxx6Ddu2yZdHfeSdrq4LT8kmSpGWXS5GdUhqeUvp0MZt1AkamlEallGYB9wP7ln86aQlNmAA9esD++0OLFlmBfdJJeaeSJEk5qsxjslcHxs73elyhrVgR0SsiBkXEoEmTJpV7OOk3d9wBjz8Of/87vPcebL553okkSVLOym0Kv4h4EVilmLfOTyk9XppdFNOWSto4pXQzcDNAUVFRidtJZWLMmOxn223hzDOzXuz11887lSRJqiTKrchOKXVZxl2MA1rP97oV8M0y7lNaNvPmZTcz9u0Lq66azX1dt64FtiRJWkBlHi7yHrBeRKwVEfWAQ4Ancs6kmuzTT2GHHeDkk2HrreGFF6B27bxTSZKkSiivKfz2j4hxwFbA0xHxfKF9tYh4BiClNAc4GXgeGA48kFIalkdeiU8+gQ4dsvmu77gDnnsuWyJdkiSpGJFS9Ru+XFRUlAb9Ok+xtCymTIEmTbLVGy+9FI46ClYp7lYDSZJU00TE4JRSsWu+VObhIlJ+fvkFzj8/660ePTqb7/q88yywJUlSqZTbjY9SlfXGG3DMMdkY7COPhBVWyDuRJEmqYuzJln6VEpx2Gmy3XdaT/fzz2eqNzZrlnUySJFUxFtnSryJg1qxs9pChQ2G33fJOJEmSqiiHi6hmmzwZ+vSBE06ALbaAG27Iiu1qYu7cuQwcOJCxY8eyxRZbsMkmm+QdSZKkGsGebNVcDz8M7drBXXfBr7PRVKMCe+zYsXTo0IE//elPvP7663Tr1o1DDz2U2bNn5x1NkqRqzyJbNc/48XDggXDQQbDaalmBfcIJeacqc8cddxwHH3wwgwYN4vbbb+eLL75g4sSJXHfddXlHkySp2rPIVs3z//4fPP10Nu/1u+9Cx455Jypz3333HW+//TZnn332b23169fnz3/+M/fcc0+OySRJqhkck62aYfRoGDcOtt0WzjgDDjgA1l03tzjff/89Q4YMYfXVV2f99dcv8/3PnDmTunXrUq9evQXaGzduzIwZM8r8eJIkaUH2ZKt6mzcPrr0WNtoom/t63jyoWze3AjulxEUXXcS6667LX/7yF3bccUd23313Jk+eXKbHWW211WjVqhUPPvjgAse+4YYb2Hvvvcv0WJIk6ffsyVb1NXw4HHssvPkm7LEH3HQT1Mr335X33nsvjz32GCNGjKBly5bMnj2bM888k169evHQQw+V2XEightvvJG9996bgQMHsvHGG/PUU08xYcIEXn755TI7jiRJKl6klPLOUOaKiorSoF9ni1DNNGwYbLYZLL88XHMN/OEPlWLmkJ122okzzjiDffbZ57e2n3/+mdVXX52RI0ey8sorl+nxxo8fz5133snYsWPp1KkTPXr0oGHDhmV6DEmSaqqIGJxSKiruPXuyVb388AOsuGI2Nd/FF2fLordsmXeq33z//fesvvrqC7Q1atSIFVZYgSlTppR5kb3qqqvSt2/fMt2nJElaPMdkq3qYMQP69oU2beDLL7Ne63PPrVQFNsDOO+/8u9k9XnvtNWrXrk2bNm3yCSVJksqcPdmq+l59NRt7/fnn2c2NTZvmnahE55xzDltvvTVTp05ln332Yfjw4Vx11VXceOON1K5dO+94kiSpjNiTraorJTj5ZNhhB5gzB158EW69NRsuUkmtttpqvPfee7Ru3Zp+/foxcuRInnnmGfbbb7+8o0mSpDLkjY+q2k46CerVg7/9DRo1yjuNJEmqQRZ146M92apavv8+u5nxnXey19ddB1dfbYEtSZIqFYtsVQ0pwYMPZrOG3HMPfPBB1l4JpuWTJElamEW2Kr9vvsmWQe/RA1q3hsGDoXfvvFNJkiSVyCK7HPz4448cdNBBtG3blg033JC33nor70hV2913w3PPweWXw9tvwyab5J1IkiRpkZzCrxycdtpp7LHHHjz00EPMmjWL6dOn5x2p6hk1CsaNg+23hzPOgAMPhHXWyTuVJElSqdiTXcZ++uknXn31VY455hgA6tWrR9NKPG9zpTN3brYM+sYbw3HHwbx5ULdumRbYv/zyC9dddx3dunXjoIMO4vHHH6c6zrIjSZLyY5FdxkaNGkXz5s056qij2HTTTTn22GP5+eef845VYd58803+8Ic/sPPOO/OnP/2Jb7/9tvQfHjYMttkm67neaScYOBBqle0lOmvWLPbYYw+eeeYZjjvuOPbaay/+9Kc/cf7555fpcSRJUs1mkV0GfvnlF5577jmeffZZpk2bxvvvv88JJ5zABx98QKNGjbj00kvzjlgh7r//frp3707nzp0577zzmDp1Kp06deLrr79e/IeHDYNNN4UvvoB774Unn4RWrco840MPPcS8efN46qmn2H///TnyyCN59dVXuemmmxgzZkyZH0+SJNVMjsleRs8//zx//OMfadu2LbVq1eLjjz9mpZVWYssttwTgoIMOqhFF9pw5c+jTpw+PP/44RUXZnOy77ror9evX5/LLL+df//pX8R+cPBmaNcum5vvHP+CII6B583LL+corr3DIIYdQa74e8pVWWokuXbrw2muvcdhhh5XbsSVJUs1hT/Yy+O677zj00EN55JFHePXVV3nllVd4+umnmTx5Mm+88QYAAwcOpF27djknLX9ffPEF9erV+63A/tXBBx/Mf//7399/YPp0OOccWGst+PLLbL7rPn3KtcCGrKAursd6zJgxrLTSSuV6bEmSVHPYk70MHnzwQbp27cq22277W1vnzp3p2rUrPXv2pGnTpqy99trcfvvtOaasGM2aNeOHH37g559/ptF8qy+OHj2aFi1aLLjxK69kNzWOHJk9NmtWYTmPPPJItt12Ww466CCKiopIKXHbbbcxYcIEunTpUmE5JElS9WaRvQx++uknmhfT89q2bVu22GILLrjgghxS5aN58+Z06dKFPn36cM0111C/fn3GjBnDBRdcwN///vdso5TgpJOgXz9Ye+3sxsadd67QnBtssAE33XQTe+65J61bt2bKlCk0aNCAp556ijp1/N9BkiSVjaiOU5cVFRWlQYMGlftxPvjgA/bZZx+GDRvGCiusAMC0adPYaKONePDBB9liiy3KPUNl8uOPP3L44Yfz9ttv06ZNG0aOHEnfvn0555xz/rfRySdD/fpwySWw3HK5ZZ05cyaDBg2iUaNGdOjQgXB5dkmStIQiYnBKqajY9yyyl80pp5zCSy+9xEknnUREcMMNN7DNNttw4403VsjxK6PRo0czfvx4NtpoIxr/8ks2Jd/JJ0PnzllvtgWtJEmqBiyyy1FKiWeffZYHHniAlBIHHXQQe+21lz2jKcH998Opp8KUKXD99dn4a0mSpGrCIlsV6+uvoXdveOop6NQJ+veHjTbKO5UkSVKZWlSR7RR+Knv33pvd1HjVVfDmmxbYkiSpxnE6BZWNkSOzHuwddsjGYB90UDYHtiRJUg1kT7aWzdy5cOWVsMkmcPzxMG8e1KljgS1Jkmo0i2wtvaFDYautspUau3TJhojU8pKSJElyuIiWzrBhsNlm0LRpNotIjx5OzSdJklRgt6OWzHffZY/t2sGll8Inn8DBB1tgS5IkzcciW6Xz889w5pnZWOtRo7Ki+swzYeWV804mSZJU6ThcRIv30kvZQjKjRsEJJ1hYS5IkLUYuPdkR0T0ihkXEvIgodgLvwnajI+LjiPgwIlxdpqKllM0YsssuULs2vPIK3HADrLBC3skkSZIqtbyGiwwFDgBeLcW2O6WUOpa0mo7KUQTUrw/nnAMffZTNgS1JkqTFymW4SEppOEB4s1zlM3EinH46nHJKNj3fv/7lTY2SJElLqLLf+JiAFyJicET0yjtMtZYS3HNPNmvIww9nU/SBBbYkSdJSKLee7Ih4EVilmLfOTyk9XsrdbJNS+iYiWgADImJESqnYISaFIrwXwBprrLFUmWussWOhd2945hno3Bn698+KbUmSJC2VciuyU0pdymAf3xQeJ0bEo0AnShjHnVK6GbgZoKioKC3rsWuU++/Pbmq85ho4+eTsJkdJkiQttUo7XCQiGkVE41+fA7uR3TCpsvDZZ/Dyy9nzM87IFpU57TQLbEmSpDKQ1xR++0fEOGAr4OmIeL7QvlpEPFPYrCXwekR8BLwLPJ1Sei6PvNXKnDlw+eXQoQOceCLMmwd16sCaa+adTJIkqdrIa3aRR4FHi2n/BuhWeD4K6FDB0aq3jz6Co4+G99+H/feH66+HWpX2ywxJkqQqyxUfa4qhQ6GoCJo1gwcfhAMPdOYQSZKkcmI3ZnU3aVL22L49XHEFDB8OBx1kgS1JklSOLLKrq2nTskVl1l4bvvgiK6pPOy3ryZYkSVK5crhIdTRgAPTqBaNHw0knQYsWeSeSJEmqUezJrk7mzYNjj4XddoN69eDVV+G666Bx47yTSZIk1SgW2dVJrVqw/PLQt282k8h22+WdSJIkqUayyK7qvv0WevSAt97KXl99Nfzzn9CgQb65JEmSajCL7KoqJbjrLmjXDp54Ips1BJw1RJIkqRKwyK6KvvoKunaFI47IiuwPP8wWmZEkSVKlYJFdFT30ELz+Olx7bXZzY9u2eSeSJEnSfJzCr6r49FP4+mvYeedsvusePaB167xTSZIkqRj2ZFd2s2dnNzJ26JDNeT1vHtSpY4EtSZJUiVlkV2YffABbbgl/+hPsvTe8/HI2TZ8kSZIqNYeLVFZDh8IWW8DKK8PDD8MBB+SdSJIkSaVkt2hlM2FC9ti+fTbn9fDhFtiSJElVjEV2ZTF1KpxyCqyzDowcmc13fcopsOKKeSeTJEnSEnK4SGXw/PPQqxeMHZsV1qusknciSZIkLQOL7DzNmwfHHgu3357Ndf3667D11nmnkiRJ0jJyuEieatWCpk3h/POzmUQssCVJkqoFi+yKNn48HHQQvPlm9vqqq+Bvf4MGDfLNJUmSpDJjkV1RUsqGhbRrB089la3gKEmSpGrJIrsijB4Nu+8ORx8NG28MQ4bAUUflnUqSJEnlxCK7Ijz0ELz1FtxwA7zyCqy/ft6JJEmSVI6cXaS8DB8OX38NXbrA6afDIYdAq1Z5p5IkSVIFsCe7rM2end3I2LEjnHpqNk1fnToW2JIkSTWIRXZZGjwYiorgggtgv/3g5ZezafokSZJUozhcpKx8/DF06gQtW8Kjj2ZFtiRJkmoki+yystFG8O9/w2GHZQvMSJIkqcayyC4rEXDSSXmnkCRJUiXggGFJkiSpjFlkS5IkSWXMIluSJEkqYxbZkiRJUhmzyJYkSZLKmEW2JEmSVMYssiVJkqQyZpEtSZIklTGLbEmSJKmMWWRLkiRJZcwiW5IkSSpjFtmSJElSGculyI6I/4uIERExJCIejYimJWy3R0R8GhEjI6JvBceUJEmSlkpePdkDgI1SSpsAnwHnLbxBRNQGrge6Au2AnhHRrkJTSpIkSUshlyI7pfRCSmlO4eXbQKtiNusEjEwpjUopzQLuB/atqIySJEnS0qoMY7KPBp4tpn11YOx8r8cV2iRJkqRKrU557TgiXgRWKeat81NKjxe2OR+YA9xT3C6KaUuLOF4voBfAGmusscR5JUmSpLJSbkV2SqnLot6PiCOAvYBdUkrFFc/jgNbzvW4FfLOI490M3AxQVFRUYjEuSZIklbcovr4t54NG7AFcBeyQUppUwjZ1yG6K3AX4GngPODSlNKwU+58EfFV4uTLwXVnkrkE8Z0vOc7bkPGdLznO25DxnS8fztuQ8Z0uuOpyzNVNKzYt7I68ieyRQH/i+0PR2Sql3RKwG3JpS6lbYrhtwDVAbuC2l9PelONaglFJR2SSvGTxnS85ztuQ8Z0vOc7bkPGdLx/O25DxnS666n7NyGy6yKCmldUto/wboNt/rZ4BnKiqXJEmSVBYqw+wikiRJUrVSE4rsm/MOUAV5zpac52zJec6WnOdsyXnOlo7nbcl5zpZctT5nuYzJliRJkqqzmtCTLUmSJFWoaldkR8T/RcSIiBgSEY9GRNMSttsjIj6NiJER0beCY1YqEdE9IoZFxLyIKPEu34gYHREfR8SHETGoIjNWNktwzrzOCiKiWUQMiIjPC48rlrBdjb/OFnfdRObfhfeHRMRmeeSsTEpxznaMiCmF6+rDiLgwj5yVSUTcFhETI2JoCe97nS2kFOfM62whEdE6Il6OiOGFvzdPK2abanmtVbsiGxgAbJRS2oRsnu3zFt4gImoD1wNdgXZAz4hoV6EpK5ehwAHAq6XYdqeUUsfqPOVOKS32nHmd/U5fYGBKaT1gYOF1SWrsdVbK66YrsF7hpxfQr0JDVjJL8P/aa4XrqmNK6eIKDVk53QHssYj3vc5+7w4Wfc7A62xhc4CzUkobAp2Bk2rKn2nVrshOKb2QUppTePk22UqRC+sEjEwpjUopzQLuB/atqIyVTUppeErp07xzVCWlPGdeZwvaF7iz8PxOYL/8olRqpblu9gXuSpm3gaYRsWpFB61E/H9tKaSUXgUmL2ITr7OFlOKcaSEppfEppfcLz6cCw4HVF9qsWl5r1a7IXsjRwLPFtK8OjJ3v9Th+/x9cv5eAFyJicET0yjtMFeB1tqCWKaXxkP2hC7QoYbuafp2V5rrx2lpQac/HVhHxUUQ8GxHtKyZaleZ1tnS8zkoQEW2ATYF3FnqrWl5ruSxGs6wi4kVglWLeOj+l9Hhhm/PJvqK4p7hdFNNWradZKc05K4VtUkrfREQLYEBEjCj8q75aKoNz5nX2P+cvwW5q1HVWjNJcNzXu2lqM0pyP98mWP54W2WrCj5F9Na2SeZ0tOa+zEkTE8sDDwOkppZ8WfruYj1T5a61KFtkppS6Lej8ijgD2AnZJxc9ROA5oPd/rVsA3ZZew8lncOSvlPr4pPE6MiEfJvqKttsVPGZwzr7P5RMSEiFg1pTS+8DXgxBL2UaOus2KU5rqpcdfWYiz2fMz/l3pK6ZmIuCEiVk4pfVdBGasir7Ml5HVWvIioS1Zg35NSeqSYTarltVbthotExB7AucA+KaXpJWz2HrBeRKwVEfWAQ4AnKipjVRQRjSKi8a/Pgd3Ibv5TybzOFvQEcETh+RHA774N8DoDSnfdPAEcXrgjvzMw5dehODXUYs9ZRKwSEVF43ons77/vKzxp1eJ1toS8zn6vcD76A8NTSleVsFm1vNaqZE/2YlwH1Cf7mhng7ZRS74hYDbg1pdQtpTQnIk4GngdqA7ellIblFzlfEbE/cC3QHHg6Ij5MKe0+/zkDWgKPFs5pHeDelNJzuYXOWWnOmdfZ71wKPBARxwBjgO4AXmcLKum6iYjehfdvBJ4BugEjgenAUXnlrQxKec4OAk6IiDnADOCQEr7prDEi4j5gR2DliBgHXATUBa+zkpTinHmd/d42wB+BjyPiw0Lbn4A1oHpfa674KEmSJJWxajdcRJIkScqbRbYkSZJUxiyyJUmSpDJmkS1JkiSVMYtsSZIkqYxZZEvSEoqIlSLiw8LPtxHxdeH5jxHxSQVn2S8i2s33+uKIWOKFlCKiTUQUOyd5RLSPiJci4rOI+CIi/hoRZf73x6J+l4h4JSKKyvqYklReLLIlaQmllL5PKXVMKXUEbgSuLjzvCMwr6+NFxKLWNNgP+K0wTSldmFJ6sQyP3ZBsoYhLU0rrAxuTrcJ5WlkdYz77UY6/iyRVJItsSSpbtSPilogYFhEvFIpUImKdiHguIgZHxGsR0bbQvmZEDIyIIYXHNQrtd0TEVRHxMnBZcZ+PiK2BfYD/K/Skr1P43EGFfWwREW9GxEcR8W5ENC70WL8WEe8XfrZezO9zKPBGSukFgMJKuicDZxeO8ZeI6PPrxhExNCLaFJ4/Vsg7LCJ6zbfNtIj4eyHX2xHRcnG/y/wiYreIeKuQ/8GIWL7QfmlEfFI4l1cs+X86SSo7FtmSVLbWA65PKbUHfgQOLLTfDJySUtoc6APcUGi/DrgrpbQJcA/w7/n2tT7QJaV0VnGfTym9SdbLfHahZ/2LXz8Y2fLi/wFOSyl1ALqQrUA3Edg1pbQZcPBCxytOe2Dw/A2F4zSMiKaL+ezRhbxFwKkRsVKhvRHZarwdgFeB4xb1u8wvIlYG/lw4L5sBg4AzI6IZsD/QvnAu/7aYbJJUrqrjsuqSlKcvU0ofFp4PBtoUelq3Bh4sLBkPUL/wuBVwQOH5/wMun29fD6aU5i7m8yXZABifUnoPIKX0E0BENAKui4iOwFyyQn5RAihuaeAopm1hp0bE/oXnrcn+AfI9MAt4qtA+GNi1FPv6VWeyISVvFM5FPeAt4CfgF+DWiHh6vv1LUi4ssiWpbM2c7/lcoCHZt4Y/FsZtL878Be3Phccl+fyvSiqOzwAmAB0K+/1lMfsZBmy/wI4j1ga+Syn9GBFzWPBb0QaFbXYk6z3fKqU0PSJe+fU9YHZK6ddsc1myv4sCGJBS6vm7NyI6AbsAh5ANadl5CfYrSWXK4SKSVM4KvchfRkR3gMh0KLz9JllRCHAY8PoSfn4q0LiYw44AVouILQqfaVy4gbIJWQ/3POCPQO3FxL8H2Ha+WT4akg0xuajw/mhgs8J7mwFrFdqbAD8UCuy2ZD3Qi1PS7zK/t4FtImLdwjGXi4j1C739TVJKzwCnk92EKkm5sciWpIpxGHBMRHxE1ju8b6H9VOCoiBhCVvSWNGtHSZ+/Hzg7Ij6IiHV+3TilNItszPW1hc8MIOtJvgE4IiLeJhsq8jOLkFKaQXZD4vkR8RnwHdmNkPcUNnkYaBYRHwInAJ8V2p8D6hR+r0vIiuPFKfZ3WSjPJOBI4L7Cvt8G2pIV508V2v5L1mMvSbmJ/31jJ0nSokXEfsBVwE4ppa9yjiNJlZZFtiRJklTGHC4iSZIklTGLbEmSJKmMWWRLkiRJZcwiW5IkSSpjFtmSJElSGbPIliRJksqYRbYkSZJUxv4/zgcKNBddZ1gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/46table.csv\")\n",
    "\n",
    "# Creating X and y\n",
    "X1 = dataset['Hydrocarbon (%)']\n",
    "y1 = dataset['Purity (%)']\n",
    "dat46 = pd.DataFrame(X1).join(pd.DataFrame(y1))\n",
    "# Change the column names\n",
    "dat46.columns =['Hydrocarbon', 'Purity']\n",
    "\n",
    "# calculate least square estimate coefficients using linear model\n",
    "model = sm.formula.ols('Purity ~ Hydrocarbon', dat46)\n",
    "lm46 = model.fit()\n",
    "display(lm46.params)\n",
    "\n",
    "from statsmodels.stats.outliers_influence import OLSInfluence as olsi\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# basic scatterplot with least squares estimate linear regression line\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.scatter(X1,y1,label='Observations', marker='x', color='black')\n",
    "\n",
    "#Create and plot the regression line\n",
    "yhat = 11.8010285*X1 + 77.863284\n",
    "plt.plot(X1, yhat, lw=1, c='blue', label = 'regression line')\n",
    "plt.xlabel('Hydrocarbon (%)', fontsize = 20)\n",
    "plt.ylabel('Purity (%)', fontsize = 20)\n",
    "plt.legend()\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as scipystats\n",
    "\n",
    "studentized_residuals = olsi(lm46).resid_studentized\n",
    "leverage_pts = olsi(lm46).hat_matrix_diag\n",
    "cook_dist = olsi(lm46).cooks_distance\n",
    "\n",
    "sorted_student_residuals = pd.Series(studentized_residuals)\n",
    "sorted_student_residuals = sorted_student_residuals.sort_values(ascending = True)\n",
    "df = pd.DataFrame(sorted_student_residuals)\n",
    "df.columns = ['sorted_student_residuals']\n",
    "df['theoretical_quantiles'] = scipystats.probplot(df['sorted_student_residuals'], dist = 'norm', fit = False)[0]\n",
    "rankings = abs(df['sorted_student_residuals']).sort_values(ascending = False)\n",
    "top3 = rankings[:3]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_figheight(8)\n",
    "fig.set_figwidth(12)\n",
    "x = df['theoretical_quantiles']\n",
    "y = df['sorted_student_residuals']\n",
    "ax.scatter(x,y, edgecolor = 'k',facecolor = 'none')\n",
    "ax.set_title('Normal Q-Q')\n",
    "ax.set_ylabel('Standardized Residuals')\n",
    "ax.set_xlabel('Theoretical Quantiles')\n",
    "ax.plot([np.min([x,y]),np.max([x,y])],[np.min([x,y]),np.max([x,y])], color = 'r', ls = '--')\n",
    "for val in top3.index:\n",
    "    ax.annotate(val,xy=(df['theoretical_quantiles'].loc[val],df['sorted_student_residuals'].loc[val]))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "028faa54-4f4a-45ac-abcc-5ecdfba2208c",
   "metadata": {
    "tags": []
   },
   "source": [
    "**b. Construct and interpret a plot of the residuals versus the predicted response.**\n",
    "\n",
    "To calculate the residuals and standardized residuals by hand:\n",
    "\n",
    " * Residuals/Error values are the difference between response values from the data (y) and predicted response from linear model (y_hat, given a value x)\n",
    " * Standardized residuals are the residuals/error values divided by the square root of the Mean Squared Error (MSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b2caab6c-d127-404a-85bd-e96c34287949",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Hydrocarbon (%)</th>\n",
       "      <th>Residual</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.87</td>\n",
       "      <td>3.729821</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.95</td>\n",
       "      <td>-3.874261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.95</td>\n",
       "      <td>-1.744261</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.98</td>\n",
       "      <td>1.131708</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.99</td>\n",
       "      <td>7.303698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1.01</td>\n",
       "      <td>5.217677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.01</td>\n",
       "      <td>2.797677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.02</td>\n",
       "      <td>-2.990333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.02</td>\n",
       "      <td>-0.040333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-4.622426</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-1.112426</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-4.672426</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1.15</td>\n",
       "      <td>-4.124467</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.40</td>\n",
       "      <td>-0.734724</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.43</td>\n",
       "      <td>0.871245</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.43</td>\n",
       "      <td>-4.458755</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1.46</td>\n",
       "      <td>1.637214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1.55</td>\n",
       "      <td>3.265122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1.55</td>\n",
       "      <td>2.505122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1.55</td>\n",
       "      <td>-0.084878</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Hydrocarbon (%)  Residual\n",
       "7              0.87  3.729821\n",
       "18             0.95 -3.874261\n",
       "5              0.95 -1.744261\n",
       "19             0.98  1.131708\n",
       "17             0.99  7.303698\n",
       "16             1.01  5.217677\n",
       "4              1.01  2.797677\n",
       "0              1.02 -2.990333\n",
       "9              1.02 -0.040333\n",
       "3              1.11 -4.622426\n",
       "1              1.11 -1.112426\n",
       "6              1.11 -4.672426\n",
       "15             1.15 -4.124467\n",
       "14             1.40 -0.734724\n",
       "8              1.43  0.871245\n",
       "2              1.43 -4.458755\n",
       "10             1.46  1.637214\n",
       "11             1.55  3.265122\n",
       "12             1.55  2.505122\n",
       "13             1.55 -0.084878"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Hydrocarbon (%)</th>\n",
       "      <th>Standardized Residual</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.02</td>\n",
       "      <td>-0.831443</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-0.309303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.43</td>\n",
       "      <td>-1.239728</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-1.285235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.01</td>\n",
       "      <td>0.777876</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.95</td>\n",
       "      <td>-0.484980</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.11</td>\n",
       "      <td>-1.299137</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.87</td>\n",
       "      <td>1.037052</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.43</td>\n",
       "      <td>0.242244</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.02</td>\n",
       "      <td>-0.011214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1.46</td>\n",
       "      <td>0.455217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1.55</td>\n",
       "      <td>0.907846</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1.55</td>\n",
       "      <td>0.696533</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1.55</td>\n",
       "      <td>-0.023600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.40</td>\n",
       "      <td>-0.204285</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>1.15</td>\n",
       "      <td>-1.146781</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1.01</td>\n",
       "      <td>1.450741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.99</td>\n",
       "      <td>2.030746</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.95</td>\n",
       "      <td>-1.077213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.98</td>\n",
       "      <td>0.314664</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Hydrocarbon (%)  Standardized Residual\n",
       "0              1.02              -0.831443\n",
       "1              1.11              -0.309303\n",
       "2              1.43              -1.239728\n",
       "3              1.11              -1.285235\n",
       "4              1.01               0.777876\n",
       "5              0.95              -0.484980\n",
       "6              1.11              -1.299137\n",
       "7              0.87               1.037052\n",
       "8              1.43               0.242244\n",
       "9              1.02              -0.011214\n",
       "10             1.46               0.455217\n",
       "11             1.55               0.907846\n",
       "12             1.55               0.696533\n",
       "13             1.55              -0.023600\n",
       "14             1.40              -0.204285\n",
       "15             1.15              -1.146781\n",
       "16             1.01               1.450741\n",
       "17             0.99               2.030746\n",
       "18             0.95              -1.077213\n",
       "19             0.98               0.314664"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/46table.csv\")\n",
    "\n",
    "# Creating X and y\n",
    "X1 = dataset['Hydrocarbon (%)']\n",
    "y1 = dataset['Purity (%)']\n",
    "dat46 = pd.DataFrame(X1).join(pd.DataFrame(y1))\n",
    "\n",
    "# Calculate residuals/error values manually\n",
    "yhat = 11.8010285*X1 + 77.863284\n",
    "Residuals = y1 - yhat\n",
    "dfErr = pd.DataFrame(X1).join(pd.DataFrame(Residuals)).sort_values(by='Hydrocarbon (%)')\n",
    "dfErr.columns = ['Hydrocarbon (%)', 'Residual']\n",
    "display(dfErr)\n",
    "\n",
    "# Calculate standardized residuals manually\n",
    "SS_error = np.sum(np.square(Residuals))\n",
    "MSE = SS_error/(20-2) #approximate variance or sigma_squared of the sample\n",
    "StandardizedResiduals = Residuals/(np.sqrt(MSE))\n",
    "dfStdRes = pd.DataFrame(X1).join(pd.DataFrame(StandardizedResiduals))\n",
    "dfStdRes.columns = ['Hydrocarbon (%)', 'Standardized Residual']\n",
    "display(dfStdRes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97408316-6354-48cc-8640-80ed191957c0",
   "metadata": {
    "tags": []
   },
   "source": [
    "Looking at the residuals plotted against predicted values, they aren't distributed randomly. The large gap in the middle seems to indicate some pattern. Therefore the data violates the assumption of constant variance (Homoscedasticity of errors) required for an adequate regression model. Standardized residuals show several outliers - observations 7:|-1.299|, 17:|1.45|, 18:|2.03| - which could be analyzed for leverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8e5db4e8-fa47-459e-9d6b-242380f70907",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHgCAYAAABaYIDJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkk0lEQVR4nO3df3Rf6V0f+Pcn8tfHEnVi0RmRXex0oiUHbQzhR92cZckRgnS3JG2d/mC34VAosHskNk0C2SVAysYl3mZ355R0Q0lKNDsks8nSZBsgrLMMS0JmVQHTpMwkmRAThTNHJLULqWZAmXGQjGXl2T+kMfZkbF9PLF1Jfr3O0bm6z/fX+1x/j/XWo/t9brXWAgAAXN8z+g4AAAC7hfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHe3rO8CNuO2229odd9zRdwwAAPa4Bx988NHW2u1PHt9V5fmOO+7IAw880HcMAAD2uKr67FONO20DAAA6Up4BAKAj5RkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADra13cAuFnmFpYyO7+YM8srOTI6kpnJ8UxNjPUdCwDYQ8w8syfMLSzlxKnTWTp3PoeGB1k6dz4nTp3O3MJS39EAgD1EeWZPmJ1fzGCoMrJ/X6o2toOhyuz8Yt/RAIA9RHlmTzizvJLhwdAVY8ODoZxdXukpEQCwFynP7AlHRkeyurZ+xdjq2noOj470lAgA2IuUZ/aEmcnxrK23rFy4mNY2tmvrLTOT431HAwD2EOWZPWFqYiwnjx/N2MEDeWx1LWMHD+Tk8aNW2wAAbipL1bFnTE2MKcsAwJYy8wwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQUa/luaoOVdUvVtVCVX2qqr6lzzwAAHAt+3p+/Z9J8v+21r6rqvYnGek5DwAAXFVv5bmqnplkMsn3J0lr7UKSC33lAQCA6+nztI3xJI8keUdVfayq7q6qr+gxDwAAXFOf5Xlfkm9O8nOttW9K8qdJfuLJd6qq6ap6oKoeeOSRR7Y7IwAAXNJneT6b5Gxr7SOb+7+YjTJ9hdbaXa21Y621Y7fffvu2BgQAgMv1Vp5ba59LcqaqvnZz6MVJfq+vPAAAcD19r7bxqiS/sLnSxmKSH+g5DwAAXFWv5bm19vEkx/rMAAAAXbnCIAAAdKQ8AwBAR8ozAAB0pDwDAEBHyjMAAHSkPAMAQEfKMwAAdKQ8AwBAR8ozAAB0pDwDAEBHyjMAAHSkPAMAQEfKMwAAdKQ8AwBAR8ozAAB0pDwDAEBHyjMAAHSkPAMAQEf7+g4AO83cwlJm5xdzZnklR0ZHMjM5nqmJsb5jAQA7gJlnuMzcwlJOnDqdpXPnc2h4kKVz53Pi1OnMLSz1HQ0A2AGUZ7jM7PxiBkOVkf37UrWxHQxVZucX+44GAOwAyjNc5szySoYHQ1eMDQ+GcnZ5padEAMBOojzDZY6MjmR1bf2KsdW19RweHekpEQCwkyjPcJmZyfGsrbesXLiY1ja2a+stM5PjfUcDAHYA5RkuMzUxlpPHj2bs4IE8trqWsYMHcvL4UattAABJLFUHX2JqYkxZBgCekplnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoqPfyXFVDVfWxqvp/+s4CAADX0nt5TvLDST7VdwgAALieXstzVR1O8teT3N1nDgAA6KLvmec3J/mxJF/sOQcAAFxXb+W5qv5GkqXW2oPXud90VT1QVQ888sgj25QOAAC+VJ8zz9+a5HhVfSbJe5J8R1X9n0++U2vtrtbasdbasdtvv327MwIAwCW9lefW2utaa4dba3ckeXmS+1prf7+vPAAAcD19n/MMAAC7xr6+AyRJa20uyVzPMQAA4JrMPAMAQEc7YuZ5p5pbWMrs/GLOLK/kyOhIZibHMzUx1ncsAAB6Yub5KuYWlnLi1OksnTufQ8ODLJ07nxOnTmduYanvaAAA9ER5vorZ+cUMhioj+/elamM7GKrMzi/2HQ0AgJ4oz1dxZnklw4OhK8aGB0M5u7zSUyIAAPqmPF/FkdGRrK6tXzG2uraew6MjPSUCAKBvyvNVzEyOZ229ZeXCxbS2sV1bb5mZHO87GgAAPVGer2JqYiwnjx/N2MEDeWx1LWMHD+Tk8aNW2wAAuIVZqu4apibGlGUAAC4x8wwAAB0pzwAA0JHyDAAAHSnPAADQkfIMAAAdKc8AANCR8gwAAB1Z5xkAgB1lbmEps/OLObO8kiOjI5mZHN8x194w8wwAwI4xt7CUE6dOZ+nc+RwaHmTp3PmcOHU6cwtLfUdLojwDALCDzM4vZjBUGdm/L1Ub28FQZXZ+se9oSZRnAAB2kDPLKxkeDF0xNjwYytnllZ4SXUl5BgBgxzgyOpLVtfUrxlbX1nN4dKSnRFdSngEA2DFmJseztt6ycuFiWtvYrq23zEyO9x0tifIMAMAOMjUxlpPHj2bs4IE8trqWsYMHcvL40R2z2oal6gAA2FGmJsZ2TFl+MjPPAADQkfIMAAAdKc8AANCR8gwAAB0pzwAA0JHyDAAAHSnPAADQkXWe4UnmFpYyO7+YM8srOTI6kpnJ8R271iQAsL3MPMNl5haWcuLU6SydO59Dw4MsnTufE6dOZ25hqe9oAMAOoDzDZWbnFzMYqozs35eqje1gqDI7v9h3NABgB1Ce4TJnllcyPBi6Ymx4MJSzyys9JQIAdhLlGS5zZHQkq2vrV4ytrq3n8OhIT4kAgJ1EeYbLzEyOZ229ZeXCxbS2sV1bb5mZHO87GgCwAyjPcJmpibGcPH40YwcP5LHVtYwdPJCTx49abQMASGKpOvgSUxNjyjIA8JTMPAMAQEfKMwAAdOS0DQAAdpSdfLVfM88AAOwYO/1qv8ozAAA7xk6/2q/yDADAjrHTr/arPAMAsGPs9Kv9Ks8AAOwYO/1qv8ozAAA7xk6/2q+l6gAA2FF28tV+zTwDAEBHyjMAAHSkPAMAQEfKMwAAdKQ8AwBAR8ozAAB0pDwDAEBHyjMAAHTUW3muqiNV9f9V1aeq6nRV/XBfWQAAoIs+rzB4Mcn/0Fr7aFUdTPJgVX2wtfZ7PWYCAICr6m3mubX2R621j25+fy7Jp5J8dV95AADgenbEOc9VdUeSb0rykZ6jAADAVfVenqvqLyT5pSQ/0lp7/Clun66qB6rqgUceeWT7AwIAwKZey3NVDbJRnH+htfbLT3Wf1tpdrbVjrbVjt99++/YGBACAy/S52kYl+fkkn2qt/bO+cgAAQFd9zjx/a5LvTfIdVfXxza+X9pgHAACuqbel6lprv5Wk+np9AAC4Ub1/YBAAAHYL5RkAADpSngEAoCPlGQAAOurtA4MAAHvV3MJSZucXc2Z5JUdGRzIzOZ6pibG+Y3ETdJp5rqqfrqqjWx0GAGC3m1tYyolTp7N07nwODQ+ydO58Tpw6nbmFpb6jcRN0PW1jIcldVfWRqvqhqnrWVoYCANitZucXMxiqjOzfl6qN7WCoMju/2Hc0boJO5bm1dndr7VuTfF+SO5J8oqr+ZVV9+1aGAwDYbc4sr2R4MHTF2PBgKGeXV3pKxM3U+QODVTWUZGLz69EkDyX576vqPVuUDQBg1zkyOpLVtfUrxlbX1nN4dKSnRNxMXc95/mfZOHXjpUn+59baX26t3dla+5tJvmkrAwIA7CYzk+NZW29ZuXAxrW1s19ZbZibH+47GTdB1tY1PJvkfW2tP9feGF97EPAAAu9rUxFhOZuPc57PLKzlstY09pWt5/p7W2tsvH6iqD7XWXtxae2wLcgEA7FpTE2PK8h51zfJcVQeSjCS5rapGk9TmTc9M8h9vcTYAANhRrjfzPJPkR7JRlD962fjjSd66RZkAAGBHumZ5bq39TJKfqapXtdZ+dpsyAQDAjnS90za+o7V2X5J/X1V/58m3t9Z+ecuSAQDADnO90za+Lcl9Sf7mU9zWkijPAADcMq532sY/rqpnJPm11tq/2qZMAACwI133IimttS8meeU2ZAEAgB2t6+W5P1hVP1pVR6rqK5/42tJkAACww3S9SMoPbm7/4WVjLYnrTLIt5haWMju/mDPLKzniSk0AQE86lefW2nO3OghczdzCUk6cOp3BUOXQ8CBL587nxKnTOZko0ADAtupUnqvq+55qvLX2zpsbB77U7PxiBkOVkf0bb9eR/fuycuFiZucXlWcAYFt1PW3jr1z2/YEkL87GFQeVZ7bcmeWVHBoeXDE2PBjK2eWVnhIBALeqrqdtvOry/ap6VpJ3bUkieJIjoyNZOnf+0sxzkqyurefw6EiPqQCAW1HX1TaebCXJ825mELiamcnxrK23rFy4mNY2tmvrLTOTPq8KAGyvruc8vz8bq2skG4X7+UlcNIVtMTUxlpPZOPf57PJKDlttAwDoSddznn/6su8vJvlsa+3sFuSBpzQ1MaYsAwC9u2Z5rqoDSX4oydck+d0kP99au7gdwQAAYKe53szz/5FkLclvJnlJNk7X+OGtDgUAcLO54BY3w/XK8/Nba1+fJFX180n+7dZHAgC4uVxwi5vleqttrD3xjdM1AIDd6vILblVtbAdDldn5xb6jsctcb+b5G6rq8c3vK8nw5n4laa21Z25pOgCAm8AFt7hZrlmeW2tD2xUEAGCruOAWN8vTvUgKAMCu4YJb3CzKMwCw501NjOXk8aMZO3ggj62uZezggZw8ftSHBblhXS+SAgCwq7ngFjeDmWcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI56Lc9V9Z1V9emqeriqfqLPLAAAcD29leeqGkry1iQvSfL8JN9dVc/vKw8AAFzPvh5f+4VJHm6tLSZJVb0nycuS/F6PmaCTuYWlzM4v5szySo6MjmRmcjxTE2N9xwIAtlifp218dZIzl+2f3Ry7qs9+9rN5//vfnyS5ePFipqenc++99yZJzp8/n+np6XzgAx9IknzhC1/I9PR07rvvviTJ5z//+UxPT2d+fj5J8uijj2Z6ejr3339/kuRzn/tcpqen85GPfGQjzNmzmZ6ezoMPPpgk+cxnPpPp6ek89NBDSZKHH34409PTOX36dJLk05/+dKanp/PpT386SXL69OlMT0/n4YcfTpI89NBDmZ6ezmc+85kkyYMPPpjp6emcPXs2SfKRj3wk09PT+dznPpckuf/++zM9PZ1HH300STI/P5/p6el8/vOfT5Lcd999mZ6ezhe+8IUkyQc+8IFMT0/n/PnzSZJ7770309PTuXjxYpLk/e9/f6anpy8dy/e97315xStecWn/ve99b1796ldf2n/3u9+d17zmNZf23/Wud+W1r33tpf177rknr3vd6y7t33333Xn9619/af9tb3tb3vCGN1zaf8tb3pI3vvGNl/bf/OY3584777y0/6Y3vSlvetObLu3feeedefOb33xp/41vfGPe8pa3XNp/wxvekLe97W2X9l//+tfn7rvvvrT/ute9Lvfcc8+l/de+9rV517vedWn/Na95Td797ndf2n/1q1+d9773vZf2X/GKV+R973vfpf3p6elL770PffIP82OveWXOfuK3c2h4kP+w/Hh+7DWvzJvv2Xi895733la99/y/573nvee9l3jvbed776n0OfNcTzHWvuROVdNJppPkWc961lZnguv633/rD1KV7N83lKrKyP59+dNKfuNTS/mRvsMBwB6w8EeP51/e9eH8uz8+l2f80eP52j98PC/tO9Smau1L+ur2vHDVtyT5qdbaX9vcf12StNb+l6s95tixY+2BBx7YpoTw1F505305NDxI1Z///tday2Ora/nNH/+OHpMBwO43t7CUE6dOZzBUGR4MZXVtPWvrLSePH93WUySr6sHW2rEnj/d52sbvJHleVT23qvYneXmSUz3mgU6OjI5kdW39irHVtfUcHh3pKREA7B2z84sZDG38ZfeJv/AOhiqz84t9R0vSY3lurV1M8sokv57kU0n+VWvtdF95oKuZyfGsrbesXLiY1ja2a+stM5PjfUcDgF3vzPJKhgdDV4wND4Zydnmlp0RX6vOc57TW7k1yb58Z4EZNTYzlZDZ+Mz67vJLDVtsAgJvmyOhIls6dz8j+P6+pO+kvvL2WZ9itpibGlGUA2AIzk+M5cep0Vi5cvOKc553yF16X5wYAYMeYmhjLyeNHM3bwQB5bXcvYwQPb/mHBazHzDADAjrKT/8Jr5hkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADra13cAAODWNbewlNn5xZxZXsmR0ZHMTI5namKs71hwVWaeAYBezC0s5cSp01k6dz6HhgdZOnc+J06dztzCUt/R4KqUZwCgF7PzixkMVUb270vVxnYwVJmdX+w7GlyV8gwA9OLM8kqGB0NXjA0PhnJ2eaWnRHB9yjMA0IsjoyNZXVu/Ymx1bT2HR0d6SgTXpzwDAL2YmRzP2nrLyoWLaW1ju7beMjM53nc0uCrlGQDoxdTEWE4eP5qxgwfy2Opaxg4eyMnjR622wY5mqToAoDdTE2PKMruKmWcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoyGobwLabW1jK7Pxiziyv5MjoSGYmx33aHoBdwcwzsK3mFpZy4tTpLJ07n0PDgyydO58Tp05nbmGp72gAcF3KM7CtZucXMxiqjOzfl6qN7WCoMju/2Hc0ALgu5RnYVmeWVzI8GLpibHgwlLPLKz0lAoDulGdgWx0ZHcnq2voVY6tr6zk8OtJTIgDoTnkGttXM5HjW1ltWLlxMaxvbtfWWmcnxvqMBwHUpz8C2mpoYy8njRzN28EAeW13L2MEDOXn8qNU2ANgVLFUHbLupiTFlGYBdycwzAAB0pDwDAEBHyjMAAHSkPAMAQEfKMwAAdKQ8AwBAR8ozAAB0pDwDAEBHyjMAAHSkPAMAQEfKMwAAdKQ8AwBAR/v6DsDuNbewlNn5xZxZXsmR0ZHMTI5namKs71gAAFvGzDNPy9zCUk6cOp2lc+dzaHiQpXPnc+LU6cwtLPUdDQBgyyjPPC2z84sZDFVG9u9L1cZ2MFSZnV/sOxoAwJZRnnlaziyvZHgwdMXY8GAoZ5dXekoEALD1lGeeliOjI1ldW79ibHVtPYdHR3pKBACw9ZRnnpaZyfGsrbesXLiY1ja2a+stM5PjfUcDANgyyjNPy9TEWE4eP5qxgwfy2Opaxg4eyMnjR622AQDsaZaq42mbmhhTlgGAW0ovM89V9U+raqGqPlFV76uqQ33kAACAG9HXaRsfTPJ1rbUXJPn9JK/rKQcAAHTWS3lurX2gtXZxc/fDSQ73kQMAAG7ETvjA4A8m+bW+QwAAwPVs2QcGq+o3kjz7KW76ydba/715n59McjHJL1zjeaaTTCfJc57znC1ICgAA3WxZeW6t/dVr3V5V/yDJ30jy4tZau8bz3JXkriQ5duzYVe8HAABbrZel6qrqO5P8eJJva625njMAALtCX+c8vyXJwSQfrKqPV9XbesoBAACd9TLz3Fr7mj5eFwAAvhw7YbUNAADYFZRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjvb1HQB2o7mFpczOL+bM8kqOjI5kZnI8UxNjfccCALaYmWe4QXMLSzlx6nSWzp3PoeFBls6dz4lTpzO3sNR3NABgiynPcINm5xczGKqM7N+Xqo3tYKgyO7/YdzQAYIspz3CDziyvZHgwdMXY8GAoZ5dXekoEAGwX5Rlu0JHRkayurV8xtrq2nsOjIz0lAgC2i/IMN2hmcjxr6y0rFy6mtY3t2nrLzOR439EAgC2mPMMNmpoYy8njRzN28EAeW13L2MEDOXn8qNU2AOAWYKk6eBqmJsaUZQC4BZl5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICOlGcAAOhIeQYAgI6UZwAA6Eh5BgCAjpRnAADoSHkGAICO9vUdgN1rbmEps/OLObO8kiOjI5mZHM/UxFjfsQAAtoyZZ56WuYWlnDh1OkvnzufQ8CBL587nxKnTmVtY6jsaAMCWUZ55WmbnFzMYqozs35eqje1gqDI7v9h3NACALdNrea6qH62qVlW39ZmDG3dmeSXDg6ErxoYHQzm7vNJTIgCArddbea6qI0n+iyT/rq8MPH1HRkeyurZ+xdjq2noOj470lAgAYOv1OfP8vyX5sSStxww8TTOT41lbb1m5cDGtbWzX1ltmJsf7jgYAsGV6Kc9VdTzJv2+tPdTH6/Plm5oYy8njRzN28EAeW13L2MEDOXn8qNU2AIA9bcuWqquq30jy7Ke46SeT/KMk/2XH55lOMp0kz3nOc25aPr58UxNjyjIAcEup1rb3rImq+vokH0ryxCfLDif5wyQvbK197lqPPXbsWHvggQe2OCFwq7FmOQBPVlUPttaOPXl820/baK39bmttrLV2R2vtjiRnk3zz9YozwFawZjkAN8I6z8AtzZrlANyI3i/PvTn7DNCLM8srOTQ8uGLMmuUAXI2ZZ+CWZs1yAG6E8gzc0qxZDsCNUJ6BW5o1ywG4Eb2f8wzQN2uWA9CVmWcAAOhIeQYAgI6UZwAA6Mg5zwDALWFuYSmz84s5s7ySI6MjmZkc93kHbpiZZwBgz5tbWMqJU6ezdO58Dg0PsnTufE6cOp25haW+o7HLKM8AwJ43O7+YwVBlZP++VG1sB0OV2fnFvqOxyyjPAMCed2Z5JcODoSvGhgdDObu80lMidivlGQDY846MjmR1bf2KsdW19RweHekpEbuV8gxsu7mFpXz3XR/Oi+68L99914edcwhsuZnJ8aytt6xcuJjWNrZr6y0zk+N9R2OXUZ6BbeVDO0AfpibGcvL40YwdPJDHVtcydvBATh4/arUNbpil6oBtNTu/mAsX1/PHX7iYC+tfzP6hZ+TggX2ZnV/0QwzYUlMTY/6f4cumPAPb6vf/w+N5/PzFPCOVoapcXG/54z+9kIvrj/cdDQCuS3kGttXaekuSPOMZlSSpSr74xZYLm+MAsJM55xnYVvv3PSNpyRdbS0vLF1tL2uY4AOxwfloB2+p5Ywdz28H92feMyvoXW/Y9o3Lbwf153tjBvqMBwHUpz8C2mpkcz2BoKM9+1oF87VcdzLOfdSCDoSHLRQGwKyjPwLayXBQAu5kPDALbznJRAOxWZp4BAKAj5RkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADpSngEAoCPlGQAAOlKeAQCgI+UZAAA6Up4BAKAj5RkAADpSngEAoKNqrfWdobOqeiTJZ/vOscVuS/Jo3yFuAY7z1nOMt4fjvPUc4+3hOG89x/jG/KXW2u1PHtxV5flWUFUPtNaO9Z1jr3Oct55jvD0c563nGG8Px3nrOcY3h9M2AACgI+UZAAA6Up53nrv6DnCLcJy3nmO8PRznrecYbw/Hees5xjeBc54BAKAjM88AANCR8tyTqvrOqvp0VT1cVT/xFLePVtX7quoTVfVvq+rr+si5m1XV26tqqao+eZXbq6r++ea/wSeq6pu3O+Nu1+EYT1TVv6mqP6uqH93ufHtFh+P8PZvv4U9U1f1V9Q3bnXG363CMX7Z5fD9eVQ9U1Yu2O+NecL3jfNn9/kpVrVfVd21Xtr2iw3t5qqoe23wvf7yqTmx3xt1Oee5BVQ0leWuSlyR5fpLvrqrnP+lu/yjJx1trL0jyfUl+ZntT7gn3JPnOa9z+kiTP2/yaTvJz25Bpr7kn1z7Gf5Lk1Ul+elvS7F335NrH+Q+SfNvm/xf/U5zX+HTck2sf4w8l+YbW2jcm+cEkd29Dpr3onlz7OD/xM/LOJL++HYH2oHtynWOc5Ddba9+4+XVyGzLtKcpzP16Y5OHW2mJr7UKS9yR52ZPu8/xs/Ged1tpCkjuq6qu2N+bu1lqbz0Z5u5qXJXln2/DhJIeq6j/annR7w/WOcWttqbX2O0nWti/V3tPhON/fWlve3P1wksPbEmwP6XCMv9D+/ENCX5HEB4aehg7/LyfJq5L8UpKlrU+093Q8xnwZlOd+fHWSM5ftn90cu9xDSf5OklTVC5P8pfiBeLN1+XeA3ea/SfJrfYfYi6rqb1fVQpJfzcbsMzdZVX11kr+d5G19Z9njvqWqHqqqX6uqo32H2W2U537UU4w9eRbjf00yWlUfz8Zv4R9LcnGLc91quvw7wK5RVd+ejfL8431n2Ytaa+9rrU0k+VvZOD2Gm+/NSX68tbbed5A97KPZuOz0NyT52SS/0m+c3Wdf3wFuUWeTHLls/3CSP7z8Dq21x5P8QLLxwbZsnNP4B9sV8BZx3X8H2C2q6gXZOA/3Ja21P+47z17WWpuvqv+kqm5rrT3ad5495liS92z82MttSV5aVRdba7/Sa6o9ZLNfPPH9vVX1L7yXb4yZ5378TpLnVdVzq2p/kpcnOXX5Harq0OZtSfLfJpm//A3PTXEqyfdtrrrxnyV5rLX2R32HghtVVc9J8stJvre19vt959mLquprNicysrkyz/4kfkm5yVprz22t3dFauyPJLyZ5heJ8c1XVsy97L78wG13Qe/kGmHnuQWvtYlW9MhufJB5K8vbW2umq+qHN29+W5D9N8s6qWk/ye9n4Uyw3oKrenWQqyW1VdTbJP04ySC4d43uTvDTJw0lWsjnTT3fXO8ZV9ewkDyR5ZpIvVtWPJHm+XwRvTIf38okkfzHJv9j8mXixtXasn7S7U4dj/Hez8cv2WpLVJH/vsg8Q0lGH48yXqcMx/q4k/11VXczGe/nl3ss3xhUGAQCgI6dtAABAR8ozAAB0pDwDAEBHyjMAAHSkPAMAQEfKM8A2qKovPGn/+6vqLTfw+Duq6pM3P9k1X3Ouqix5B3AZ5RlgF6uqLVmvv6qGtuJ5AXY75RmgR1V1sKr+oKoGm/vPrKrPVNWgqv5yVT1UVf8myT+87DHfX1Xvrar3J/lAVX1lVf1KVX2iqj68eanuVNVfqKp3VNXvbt72dzfHf66qHqiq01X1hsue9zNVdaKqfivJf7U5/Per6v6q+uTm1chyjdf7qap6++aM9WJVvXo7jiHAdnKFQYDtMVxVH79s/yuTnGqtnauquSR/PcmvJHl5kl9qra1V1TuSvKq19q+r6p8+6fm+JckLWmt/UlU/m+RjrbW/VVXfkeSdSb4xyeuzcdn5r0+SqhrdfOxPbj5uKMmHquoFrbVPbN52vrX2os37/1CSr2it/edVNZnk7Um+LskbrvJ6STKR5NuTHEzy6ar6udba2pd15AB2EDPPANtjtbX2jU98ZeOS2k+4O39+efgfSPKOqnpWkkOttX+9Of6uJz3fB1trf7L5/YueuL21dl+Sv7j5+L+a5K1PPKC1trz57X9dVR9N8rEkR5M8/7Ln/b+e9Drv3nzsfJJnVtWha7xekvxqa+3PWmuPJllK8lXXPiwAu4vyDNCz1tpvJ7mjqr4tyVBr7ZNJKkm7xsP+9LLv66me9qmeo6qem+RHk7y4tfaCJL+a5MBVnveJ53mq532q10uSP7tsbD3+wgnsMcozwM7wzmzM8r4jSVprn0/yWFW9aPP277nGY+efuL2qppI82lp7PMkHkrzyiTttnrbxzGwU5Meq6quSvOQ6uf7e5mNflI1TQB67xusB7HlmBAB2hl9I8k+yeZrEph9I8vaqWkny69d47E9l41SPTyRZSfIPNsf/SZK3bi5xt57kDa21X66qjyU5nWQxyW9fJ9dyVd2fjdL9g9d5PYA9r1q71l8FAdgOVfVdSV7WWvvevrMAcHVmngF6trlaxkuSvLTvLABcm5lnAADoyAcGAQCgI+UZAAA6Up4BAKAj5RkAADpSngEAoCPlGQAAOvr/AQstRoWAC0iyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([-0.86439713, -0.31816538, -1.31228081, -1.32206231,  0.81008566,\n",
       "       -0.5114297 , -1.33636284,  1.1193882 ,  0.25642111, -0.01165879,\n",
       "        0.48590043,  1.00060757,  0.76770304, -0.0260111 , -0.21467098,\n",
       "       -1.17718719,  1.51081239,  2.12273817, -1.13596082,  0.32958687])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAHgCAYAAABw0HFmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApKklEQVR4nO3dfZRdaV0n+u9vQhhKBAJ2C51A2+1M33iRAdIGBGEp+DIRBlcHRIWrgjjSF+XNe8dcyThX0PGu4a44M8qLvFwHsK8uHB1CaKW1UFBaVEbSBAjQZuyFOHQFpXkJyKXUdPzdP+pUm2SqKqc6+9Q5Vfl81jqrzn72Pmf/zk5117eeep5nV3cHAAAYzj+adgEAALDVCNkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAA7vHtAuYhMsuu6yvuuqqaZcBAMAWdsstt3y6uy9fad+WDNlXXXVVjh49Ou0yAADYwqrqL1bbZ7gIAAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAA7vHtAuAjXTk2EIOzZ/IyVOL2bljLgf27c7+PbumXRYAsMUI2VwyjhxbyMHDx7N4+kySZOHUYg4ePp4kgjYAMCjDRbhkHJo/cVfAXrZ4+kwOzZ+YUkUAwFYlZHPJOHlqcV3tAAB3l5DNJWPnjrl1tQMA3F1CNpeMA/t2Z277tnPa5rZvy4F9u6dUEQCwVZn4yCVjeXKj1UUAgEkTsrmk7N+zS6gGACbOcBEAABiYkA0AAAMTsgEAYGBCNgAADEzIBgCAgU0tZFfVQ6rq96rq1qr6SFW9eIVjqqpeUVW3VdWHquraadQKAADrMc0l/O5M8q+6+/1VdZ8kt1TV73T3R8865klJrhk9vj7Ja0ZfAQBgZk2tJ7u7P9nd7x89/+sktyY5fwHj65Lc0Evem2RHVV2xwaUCAMC6zMSY7Kq6KsmeJP/1vF27knzirO3b8z8GcQAAmClTD9lV9eVJ3pLkR7v7C+fvXuElvcr7XF9VR6vq6B133DF0mQAAMLaphuyq2p6lgP0r3X14hUNuT/KQs7YfnOTkSu/V3a/v7r3dvffyyy8fvlgAABjTNFcXqST/Kcmt3f0fVjnsxiTPGq0y8pgkn+/uT25YkQAAcDdMc3WRxyX5/iTHq+oDo7Z/neTKJOnu1ya5KcmTk9yW5EtJnrPxZQIAwPpMLWR393uy8pjrs4/pJM/fmIoAAGAYU5/4CAAAW42QDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMDuMc2TV9Ubkjwlyae6+2Er7H9Ckrcl+fNR0+Hu/ukNKxBWceTYQg7Nn8jJU4vZuWMuB/btzv49u6ZdFgAwI6YaspO8KcmrktywxjF/0N1P2Zhy4MKOHFvIwcPHs3j6TJJk4dRiDh4+niSCNgCQZMrDRbr75iSfnWYNsF6H5k/cFbCXLZ4+k0PzJ6ZUEQAwazbDmOzHVtUHq+q3quprp10MnDy1uK52AODSM+sh+/1Jvqq7H5HklUmOrHZgVV1fVUer6ugdd9yxUfVxCdq5Y25d7QDApWemQ3Z3f6G7vzh6flOS7VV12SrHvr6793b33ssvv3xD6+TScmDf7sxt33ZO29z2bTmwb/eUKgIAZs20Jz6uqaoelOSvurur6tFZ+qXgM1Mui0vc8uRGq4sAAKuZ9hJ+b07yhCSXVdXtSV6aZHuSdPdrkzw9yQ9X1Z1JFpM8o7t7SuXCXfbv2SVUAwCrmmrI7u5nXmD/q7K0xB8AAGwaMz0mGwAANiMhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwO4x7QK2giPHFnJo/kROnlrMzh1zObBvd/bv2TXtsgAAmBIh+yIdObaQg4ePZ/H0mSTJwqnFHDx8PEkEbQCAS5ThIhfp0PyJuwL2ssXTZ3Jo/sSUKgIAYNqE7It08tTiutoBANj6hOyLtHPH3LraAQDY+oTsi3Rg3+7Mbd92Ttvc9m05sG/3lCoCAGDaTHy8SMuTG60uAgDAMiF7APv37BKqAQC4i+EiAAAwMCEbAAAGJmQDAMDAhGwAABiYkA0AAAObasiuqjdU1aeq6sOr7K+qekVV3VZVH6qqaze6RgAAWK9pL+H3piSvSnLDKvuflOSa0ePrk7xm9BUAgEvckWMLM3uvkqn2ZHf3zUk+u8Yh1yW5oZe8N8mOqrpiY6oDAGBWHTm2kIOHj2fh1GI6ycKpxRw8fDxHji1Mu7Qksz8me1eST5y1ffuoDQCAS9ih+RNZPH3mnLbF02dyaP7ElCo616yH7FqhrVc8sOr6qjpaVUfvuOOOCZcFAMA0nTy1uK72jTbrIfv2JA85a/vBSU6udGB3v76793b33ssvv3xDigMAYDp27phbV/tGm/WQfWOSZ41WGXlMks939yenXRQAANN1YN/uzG3fdk7b3PZtObBv95QqOtdUVxepqjcneUKSy6rq9iQvTbI9Sbr7tUluSvLkJLcl+VKS50ynUgAAZsnyKiKzurpIda84xHlT27t3bx89enTaZQAAsIVV1S3dvXelfbM+XAQAADYdIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMbF0hu6ruX1UPn1QxAACwFVwwZFfV71fVfavqAUk+mOSNVfUfJl8aAABsTuP0ZN+vu7+Q5GlJ3tjdX5fkWydbFgAAbF7jhOx7VNUVSb47yW9OuB4AANj0xgnZP51kPslt3f2+qvrqJH822bIAAGDzuseFDujuX0/y62dtfyzJd06yKAAA2MxWDdlV9cokvdr+7n7RRCoCAIBNbq2e7KMbVgVsMkeOLeTQ/ImcPLWYnTvmcmDf7uzfs2vaZQEAM2LVkN3dv7SRhcBmceTYQg4ePp7F02eSJAunFnPw8PEkEbQBgCRjjMmuqsuT/HiShya513J7d3/zBOuCmXVo/sRdAXvZ4ukzOTR/QsgGAJKMt7rIryS5NcnVSX4qyceTvG+CNcFMO3lqcV3tAMClZ5yQ/RXd/Z+SnO7ud3f3DyZ5zITrgpm1c8fcutoBgEvPOCH79OjrJ6vqX1TVniQPnmBNMNMO7Nudue3bzmmb274tB/btnlJFAHBpOnJsIY97+bty9Uvense9/F05cmxh2iXd5YJjspP8TFXdL8m/SvLKJPdN8r9NtCqYYcvjrq0uAgDTM+sLEVT3qkthb1p79+7to0etQAgAsFU97uXvysIK86F27ZjLH75kY9bnqKpbunvvSvvGWV3kjVnhpjSjsdkAALDhZn0hgnGGi/zmWc/vleSpSU5OphwAALiwnTvmVuzJnpWFCC448bG733LW41eSfHeSh02+NAAAWNmsL0QwTk/2+a5JcuXQhQAAwLhmfSGCccZk/3WWxmTX6OtfZukOkAAAMDX79+yamVB9vguG7O6+z0YUAgAAW8WqIbuqrl3rhd39/uHLAQCAzW+tnux/P/p6ryR7k3wwS0NGHp7kvyZ5/GRLAwCAzWnV1UW6+4nd/cQkf5Hk2u7e291fl2RPkts2qkAAANhsLriEX5Kv6e7jyxvd/eEkj5xYRQAAsMmNs4TfrVX1i0l+OUuri3xfklsnWhUAAGxi44Ts5yT54SQvHm3fnOQ1E6sIAAA2uXGW8PubJP9x9AAAAC5grSX8fq27v7uqjmdpmMg5uvvhE60MAAA2qbV6speHhzxlIwoBAICtYq0l/D45evrpJJ/o7r9I8o+TPCLJyQ2oDQAANqVxlvC7Ocm9qmpXkndmaSLkmyZZFAAAbGbjhOzq7i8leVqSV3b3U5M8dLJlAQDA5jVWyK6qxyb53iRvH7WNs/QfAABcksYJ2T+a5GCSt3b3R6rqq5P83kSrAgCATWycdbLfneTdVXXv0fbHkrxo0oUBAMBmdcGe7Kp6bFV9NKNbqVfVI6rqFyZeGQAAbFLjDBf5uST7knwmSbr7g0m+cYI1AQDApjZOyE53f+K8pjMTqAUAALaEcVYJ+URVfUOSrqp7Zmk89q2TLQsAADavcXqyn5fk+Ul2Jbk9ySOT/MgEawIAgE3tgiG7uz/d3d/b3Q/s7q9M8sIkPzzEyavq26vqRFXdVlUvWWH/E6rq81X1gdHjJ4c4LwAATNKqIbuqHlJVr6+q36yqf1lVX1ZVP5vkRJKvvNgTV9W2JK9O8qQs3UHymVW10p0k/6C7Hzl6/PTFnhcAACZtrZ7sG5KcTPLKJA9L8t4sDRl5eHe/eIBzPzrJbd39se7+uyS/muS6Ad4XAACmaq2Jjw/o7peNns9X1V8leVR3/+1A596V5OxVS25P8vUrHPfYqvpglgL/j3X3RwY6PwDA1Bw5tpBD8ydy8tRidu6Yy4F9u7N/z65pl8VA1lxdpKrun6RGm3+Z5MvOuvPjZy/y3LVCW5+3/f4kX9XdX6yqJyc5kuSaVWq9Psn1SXLllVdeZGkAAJNz5NhCDh4+nsXTS6siL5xazMHDx5NE0N4i1houcr8kt5z1uG+WQu8tSY4OcO7bkzzkrO0HZ6m3+i7d/YXu/uLo+U1JtlfVZSu9WXe/vrv3dvfeyy+/fIDyAAAm49D8ibsC9rLF02dyaP7ElCpiaKv2ZHf3VRM+9/uSXFNVVydZSPKMJP/L2QdU1YOS/FV3d1U9Oku/FHxmwnUBAEzUyVOL62pn8xnnZjQT0d13VtULkswn2ZbkDd39kap63mj/a5M8PckPV9WdSRaTPKO7zx9SAgCwqezcMZeFFQL1zh1zU6iGSZhayE7uGgJy03ltrz3r+auSvGqj6wIAmKQD+3afMyY7Sea2b8uBfbunWBVDmmrIBgC4FC1PbrS6yNa1asiuqges9cIBVhcBALhk7d+zS6jewtbqyb4lS0vqVZIrk3xu9HxHkv+e5OpJFwcAAJvRqkv4dffV3f3VWZqY+B3dfVl3f0WSpyQ5vFEFAgDAZrPWOtnLHjWaoJgk6e7fSvJNkysJAAA2t3EmPn66qv5Nkl/O0vCR74u1qgEAYFXj9GQ/M8nlSd46elw+agMAAFZwwZ7s0SoiL66qL1++xTkAALC6C/ZkV9U3VNVHk3x0tP2IqvqFiVcGAACb1DjDRf5jkn0ZjcPu7g8m+cZJFgUAAJvZOCE73f2J85rOrHggAAAw1uoin6iqb0jSVXXPJC9KcutkywIAgM1rnJD9vCQ/n2RXktuTvCPJ8ydZFNwdR44t5ND8iZw8tZidO+ZyYN9ut6sFAKZinJA9193fe3ZDVT1oQvXA3XLk2EIOHj6exdNLI5kWTi3m4OHjSSJoAwAbbpwx2X9eVW+uqrmz2m5a9WiYgkPzJ+4K2MsWT5/JofkTU6oIALiUjROyjyf5gyTvqap/MmqryZUE63fy1OK62gEAJmmckN3d/QtZmvD4G1X1HVm6vTrMjJ075tbVDgAwSeOE7EqS7v7DJN+S5ECSr5lkUbBeB/btztz2bee0zW3flgP7dk+pIgDgUjbOxMcnLz/p7k9W1Tcn+YbJlQTrtzy50eoiAMAsWDVkV9X3dfcvJ3lm1YpDsG+eWFVwN+zfs0uoBgBmwlo92fcefb3PRhQCADBt7rnAUFYN2d39utHXn9q4cgAApsM9FxjSWsNFXrHWC7v7RcOXAwAwHWvdc0HIZr3WWl3kltHjXkmuTfJno8cjk5xZ/WUAAJuPey4wpLWGi/xSklTVDyR5YnefHm2/Nsk7NqQ6AIANsnPHXBZWCNTuucDdMc462Ttz7uTHLx+1AQBsGe65wJDGWSf75UmOVdXvjba/KcnLJlYRAMAUuOcCQ1ozZFfVP0pyIsnXjx5J8pLu/stJFwYAsNHcc4GhrBmyu/vvq+rfd/djk7xtg2oCAIBNbZwx2e+oqu+sVW77CAAAnGucMdn/e5bu/nhnVf1NkkrS3X3fiVYGAACb1AVDdne7rToAAKzDOD3Zqar7J7kmSzemSZJ0982TKgoAADazC4bsqvqhJC9O8uAkH0jymCR/nOSbJ1oZAABsUuNMfHxxkkcl+YvufmKSPUnumGhVAACwiY0Tsv+mu/8mSarqH3f3nyZx6yMAAFjFOGOyb6+qHUmOJPmdqvpckpOTLAoAADazcVYXeero6ctGt1a/X5LfnmhVAACwia0asqvqASs0Hx99/fIkn51IRQAAsMmt1ZN9S5LO0s1nrkzyudHzHUn+e5KrJ10cAABsRqtOfOzuq7v7q5PMJ/mO7r6su78iyVOSHN6oAgEAYLMZZ3WRR3X3Tcsb3f1bSb5pciUBAMDmNs7qIp+uqn+T5JezNHzk+5J8ZqJVAQDAJjZOT/Yzk1ye5K1ZWsbvK0dtAADACsZZwu+zWbrrIwAAMIYLhuyq+p+S/FiSq84+vru/eXJlAQDA5jXOmOxfT/LaJL+Y5MxkywEAgM1vnJB9Z3e/ZuKVAADAFjHOxMffqKofqaorquoBy4+JVwYAAJvUOCH72UkOJPmjLN0F8pYkR4c4eVV9e1WdqKrbquolK+yvqnrFaP+HquraIc4LAACTNM7qIhO5fXpVbUvy6iTfluT2JO+rqhu7+6NnHfakJNeMHl+f5DWjrwAAMLPGGZOdqnpYkocmuddyW3ffcJHnfnSS27r7Y6Nz/GqS65KcHbKvS3JDd3eS91bVjqq6ors/eZHnhok7cmwhh+ZP5OSpxezcMZcD+3Zn/55d0y4LANgI3b3mI8lLk/xekr9K8sYkf5nkv1zodWO879OT/OJZ29+f5FXnHfObSR5/1vY7k+xd5f2uz9IwlqP3u9/9Okt3p+wkffTo0T569Og5bS996Uu7u/uKK664q+3aa6/t7u7nPve55xy7sLDQN9544zltr3vd63oU/u96POUpT+nu7qc85SnntHd3v+51rzun7cYbb+yFhYVz2p773Od2d/e11157V9sVV1zR3d0vfelLfaZN+Jke9Oyf6wc9++e21Gfaiv9OPpPP5DP5TD7T5v9MD3/Ba/pnf/mmjf5MR3uVrFtLr1tdVR1P8ogkx7r7EVX1wCyF4+9Y84UXUFXflWRfd//QaPv7kzy6u1941jFvT/Lvuvs9o+13Jvk/uvuWtd577969ffToIMPG4W553MvflYVTi/9D+64dc/nDl1hiHgAu1pFjCzl4+HgWT//DCtNz27fl3z3tn23YX46r6pbu3rvSvnEmPi52998nubOq7pvkU0m+eoC6bk/ykLO2H5zk5N04BmbOyRUC9lrtAMD6HJo/cU7ATpLF02dyaP7ElCo61zgh+2hV7Ujy/2RpZZH3J/mTAc79viTXVNXVVXXPJM9IcuN5x9yY5FmjVUYek+TzbTw2m8DOHXPragcA1mfWO7QuGLK7+0e6+1R3vzZLK4E8u7ufc7En7u47k7wgyXySW5P8Wnd/pKqeV1XPGx12U5KPJbktSyH/Ry72vLARDuzbnbnt285pm9u+LQf27Z5SRQCwtcx6h9YFVxepqnd297ckSXd//Py2i9HdN2UpSJ/d9tqznneS51/seWCjLY8Fs7oIAEzGgX27VxyTPSsdWquG7Kq6V5IvS3JZVd0/SY123TfJzg2oDTa1/Xt2CdUAMCGz3qG1Vk/2/5rkR7MUqG/JP4TsL2TpJjIAADA1s9yhtWrI7u6fT/LzVfXC7n7lBtYEAACb2qoTH6vqUVX1oOWAXVXPqqq3VdUrquoBG1ciAABsLmutLvK6JH+XJFX1jUlenuSGJJ9P8vrJlwYAAJvTWmOyt3X3Z0fPvyfJ67v7LUneUlUfmHhlAACwSa3Vk72tqpZD+LckeddZ+y649B8AAFyq1grLb07y7qr6dJLFJH+QJFX1T7M0ZAQAAFjBWquL/F9V9c4kVyR5x+jGMMlS7/cLN6I4AADYjNYc9tHd712h7b9NrhwAANj81hqTDQAA3A1CNgAADEzIBgCAgQnZAAAwMCEbAAAGJmQDAMDA3LkRAJhpR44t5ND8iZw8tZidO+ZyYN/u7N+za9plwZqEbABgZh05tpCDh49n8fSZJMnCqcUcPHw8SQRtZprhIgDAzDo0f+KugL1s8fSZHJo/MaWKYDxCNgAws06eWlxXO8wKIRsAmFk7d8ytqx1mhZANAMysA/t2Z277tnPa5rZvy4F9u6dUEYzHxEcAYGYtT260ugibjZANAMy0/Xt2CdVsOoaLAADAwIRsAAAYmJANAAADE7IBAGBgJj4CM+vIsQUrCgCwKQnZwEw6cmwhBw8fv+t2ygunFnPw8PEkEbQBmHmGiwAz6dD8ibsC9rLF02dyaP7ElCoCgPEJ2cBMOnlqcV3tADBLhGxgJu3cMbeudgCYJUI2MJMO7Nudue3bzmmb274tB/btnlJFADA+Ex+BmbQ8udHqIgBsRkI2MLP279klVAOwKRkuAgAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICB3WMaJ62qByT5z0muSvLxJN/d3Z9b4biPJ/nrJGeS3NndezeuSgAAuHum1ZP9kiTv7O5rkrxztL2aJ3b3IwVsAAA2i6n0ZCe5LskTRs9/KcnvJ/nxKdXChB05tpBD8ydy8tRidu6Yy4F9u7N/z65plwUAMDHT6sl+YHd/MklGX79yleM6yTuq6paqun6tN6yq66vqaFUdveOOOwYul7vryLGFHDx8PAunFtNJFk4t5uDh4zlybGHapQEATMzEQnZV/W5VfXiFx3XreJvHdfe1SZ6U5PlV9Y2rHdjdr+/uvd299/LLL7/o+hnGofkTWTx95py2xdNncmj+xJQqAgCYvIkNF+nub11tX1X9VVVd0d2frKorknxqlfc4Ofr6qap6a5JHJ7l5IgUzESdPLa6rHQBgK5jWcJEbkzx79PzZSd52/gFVde+qus/y8yT/PMmHN6xCBrFzx9y62gEAtoJpheyXJ/m2qvqzJN822k5V7ayqm0bHPDDJe6rqg0n+JMnbu/u3p1Itd9uBfbszt33bOW1z27flwL7dU6oIAGDyprK6SHd/Jsm3rNB+MsmTR88/luQRG1waA1teRcTqIgDApWRaS/hxCdm/Z5dQDQBcUtxWHQAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwIRsAAAYmJANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMLB7TLsA2KqOHFvIofkTOXlqMTt3zOXAvt3Zv2fXtMsCADaAkA0TcOTYQg4ePp7F02eSJAunFnPw8PEkEbQB4BJguAhMwKH5E3cF7GWLp8/k0PyJKVUEAGwkIRsm4OSpxXW1AwBbi5ANE7Bzx9y62gGArUXIhgk4sG935rZvO6dtbvu2HNi3e0oVAQAbycRHmIDlyY1WFwGAS5OQDROyf88uoRoALlGGiwAAwMCEbAAAGJiQDQAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADCwqYTsqvquqvpIVf19Ve1d47hvr6oTVXVbVb1kI2sEAIC7a1o92R9O8rQkN692QFVtS/LqJE9K8tAkz6yqh25MeQAAcPdN5Y6P3X1rklTVWoc9Oslt3f2x0bG/muS6JB+deIEAAHARZnlM9q4knzhr+/ZRGwAAzLSJ9WRX1e8medAKu36iu982zlus0NZrnO/6JNcnyZVXXjlWjQAAMAkTC9nd/a0X+Ra3J3nIWdsPTnJyjfO9Psnrk2Tv3r2rhnEAAJi0WR4u8r4k11TV1VV1zyTPSHLjlGsCAIALmtYSfk+tqtuTPDbJ26tqftS+s6puSpLuvjPJC5LMJ7k1ya9190emUS8AAKzHtFYXeWuSt67QfjLJk8/avinJTRtYGgAAXLRZHi4CAACbkpANAAADE7IBAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABjaVdbK5tBw5tpBD8ydy8tRidu6Yy4F9u7N/z65plwUAMDFCNhN15NhCDh4+nsXTZ5IkC6cWc/Dw8SQRtAGALctwESbq0PyJuwL2ssXTZ3Jo/sSUKgIAmDwhm4k6eWpxXe0AAFuBkM1E7dwxt652AICtQMhmog7s25257dvOaZvbvi0H9u2eUkVw9x05tpDHvfxdufolb8/jXv6uHDm2MO2SAJhRJj4yUcuTG60uwmZnEi8A6yFkM3H79+wSQtj01prE6/sbgPMZLgIwBpN4AVgPIRtgDCbxArAeQjbAGEziBWA9jMkGGINJvACsh5ANMCaTeAEYl+EiAAAwMCEbAAAGJmQDAMDAhGwAABiYkA0AAAMTsgEAYGBCNgAADMw62QAAI0eOLbjpFIMQsoGZ5YcdsJGOHFvIwcPHs3j6TJJk4dRiDh4+niT+38O6GS4CzKTlH3YLpxbT+YcfdkeOLUy7NGCLOjR/4q6AvWzx9Jkcmj8xpYrYzIRsYCb5YQdstJOnFtfVDmsRsoGZtLDKD7XV2gEu1s4dc+tqh7UI2cBM2la1rnaAi3Vg3+7Mbd92Ttvc9m05sG/3lCpiMzPxEZhJZ7rX1Q5wsZYnN5pwzRCEbGAm7doxt+LQkF3+bAtM0P49u4RqBmG4CDCT/NkWgM1MTzYwk/zZFoDNTMgGZpY/2wKwWRkuAgAAAxOyAQBgYEI2AAAMTMgGAICBCdkAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAQBgYFMJ2VX1XVX1kar6+6rau8ZxH6+q41X1gao6upE1AgDA3XWPKZ33w0meluR1Yxz7xO7+9ITrAQCAwUwlZHf3rUlSVdM4PQAATNSsj8nuJO+oqluq6vppFwMAAOOYWE92Vf1ukgetsOsnuvttY77N47r7ZFV9ZZLfqao/7e6bVznf9UmuT5Irr7zybtUMAABDqO6e3smrfj/Jj3X3BSc1VtXLknyxu392jGPvSPIXF13gbLssibHqk+UabwzXefJc443hOk+ea7wxXOfxfVV3X77SjmlNfLygqrp3kn/U3X89ev7Pk/z0OK9d7cNuJVV1tLtXXZmFi+cabwzXefJc443hOk+ea7wxXOdhTGsJv6dW1e1JHpvk7VU1P2rfWVU3jQ57YJL3VNUHk/xJkrd3929Po14AAFiPaa0u8tYkb12h/WSSJ4+efyzJIza4NAAAuGizvroIq3v9tAu4BLjGG8N1njzXeGO4zpPnGm8M13kAU534CAAAW5GebAAAGJiQPcOq6tur6kRV3VZVL1lh//2r6q1V9aGq+pOqetg06tzsquoNVfWpqvrwKvurql4x+nf4UFVdu9E1bnZjXOOvqao/rqq/raof2+j6tooxrvP3jr6HP1RVf1RV5r2s0xjX+LrR9f1AVR2tqsdvdI1bwYWu81nHPaqqzlTV0zeqtq1ijO/lJ1TV50ffyx+oqp/c6Bo3OyF7RlXVtiSvTvKkJA9N8syqeuh5h/3rJB/o7ocneVaSn9/YKreMNyX59jX2PynJNaPH9UleswE1bTVvytrX+LNJXpTkguvgs6Y3Ze3r/OdJvmn0/4x/G+Mu7443Ze1r/M4kj+juRyb5wSS/uAE1bUVvytrXefnn5P+dZH4jCtqC3pQLXOMkf9Ddjxw9xlpGmX8gZM+uRye5rbs/1t1/l+RXk1x33jEPzdL/0NPdf5rkqqp64MaWufmN7iL62TUOuS7JDb3kvUl2VNUVG1Pd1nCha9zdn+ru9yU5vXFVbT1jXOc/6u7PjTbfm+TBG1LYFjLGNf5i/8Nkp3snMfHpbhjj/8tJ8sIkb0nyqclXtPWMeY25CEL27NqV5BNnbd8+ajvbB5M8LUmq6tFJvip+aE7COP8WsNn8yyS/Ne0itqLRvSD+NMnbs9SbzcCqaleSpyZ57bRr2eIeW1UfrKrfqqqvnXYxm42QPbtqhbbze0RenuT+VfWBLP1GfyzJnROu61I0zr8FbBpV9cQshewfn3YtW1F3v7W7vybJ/iwNy2F4P5fkx7v7zLQL2cLen6Vbhj8iySuTHJluOZvPzN5Wndye5CFnbT84ycmzD+juLyR5TrI0OS9L4y3/fKMKvIRc8N8CNouqeniWxgk/qbs/M+16trLuvrmq/klVXdbdn552PVvM3iS/uvSjL5cleXJV3dndR6Za1RYyyhjLz2+qql/wvbw+erJn1/uSXFNVV1fVPZM8I8mNZx9QVTtG+5Lkh5LcfPZ/FAzmxiTPGq0y8pgkn+/uT067KFivqroyyeEk39/d/23a9WxFVfVPR50eGa1EdM8kfpkZWHdf3d1XdfdVSf5Lkh8RsIdVVQ8663v50VnKjL6X10FP9ozq7jur6gVZmjW9LckbuvsjVfW80f7XJvmfk9xQVWeSfDRLf/5lnarqzUmekOSyqro9yUuTbE/uus43JXlyktuSfCmjvx4wvgtd46p6UJKjSe6b5O+r6keTPNQvjeszxvfyTyb5iiS/MPrZeWd3751OtZvTGNf4O7P0S/npJItJvuesiZCMaYzrzEUa4xo/PckPV9WdWfpefobv5fVxx0cAABiY4SIAADAwIRsAAAYmZAMAwMCEbAAAGJiQDQAAAxOyAWZEVX3xvO0fqKpXreP1V1XVh4evbM1z/n5VWQYQ4DxCNsAWV1UTuSdCVW2bxPsCbAVCNsCMq6r7VNWfV9X20fZ9q+rjVbW9qr6uqj5YVX+c5PlnveYHqurXq+o3kryjqh5QVUeq6kNV9d7R7dVTVV9eVW+squOjfd85an9NVR2tqo9U1U+d9b4fr6qfrKr3JPmuUfP3VdUfVdWHR3eGyxrne1lVvWHUA/6xqnrRRlxDgI3mjo8As2Ouqj5w1vYDktzY3X9dVb+f5F8kOZLkGUne0t2nq+qNSV7Y3e+uqkPnvd9jkzy8uz9bVa9Mcqy791fVNye5Ickjk/yfST7f3f8sSarq/qPX/sTodduSvLOqHt7dHxrt+5vufvzo+OcluXd3f0NVfWOSNyR5WJKfWuV8SfI1SZ6Y5D5JTlTVa7r79EVdOYAZoycbYHYsdvcjlx9Zug36sl9M8pzR8+ckeWNV3S/Jju5+96j9/z3v/X6nuz87ev745f3d/a4kXzF6/bcmefXyC7r7c6On311V709yLMnXJnnoWe/7n887z5tHr705yX2rasca50uSt3f333b3p5N8KskD174sAJuPkA2wCXT3Hya5qqq+Kcm27v5wkkrSa7zs/zvrea30tiu9R1VdneTHknxLdz88yduT3GuV911+n5Xed6XzJcnfntV2Jv6qCmxBQjbA5nFDlnqN35gk3X0qyeer6vGj/d+7xmtvXt5fVU9I8unu/kKSdyR5wfJBo+Ei981SkP58VT0wyZMuUNf3jF77+CwNPfn8GucDuCToPQDYPH4lyc9kNDxj5DlJ3lBVX0oyv8ZrX5alISYfSvKlJM8etf9MklePlv47k+SnuvtwVR1L8pEkH0vyhxeo63NV9UdZCuc/eIHzAVwSqnutvzQCMCuq6ulJruvu7592LQCsTU82wCYwWh3kSUmePO1aALgwPdkAADAwEx8BAGBgQjYAAAxMyAYAgIEJ2QAAMDAhGwAABiZkAwDAwP5/7ou3ynA9L7wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.formula.api import ols\n",
    "import seaborn as sns\n",
    "\n",
    "dataset = pd.read_csv(\"https://raw.githubusercontent.com/forestwhite/Regression/main/data/46table.csv\")\n",
    "\n",
    "# Creating X and y\n",
    "X1 = dataset['Hydrocarbon (%)']\n",
    "y1 = dataset['Purity (%)']\n",
    "dat46 = pd.DataFrame(X1).join(pd.DataFrame(y1))\n",
    "dat46.columns = ['Hydrocarbon', 'Purity']\n",
    "\n",
    "#define figure size\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "#view residual plot by itself\n",
    "sns.residplot(x='Hydrocarbon', y='Purity', data=dat46)\n",
    "\n",
    "# alternately, produce four diagnostic plots, including residuals\n",
    "# import statsmodels.api as sm\n",
    "# fig = sm.graphics.plot_regress_exog(lm46, 'Hydrocarbon', fig=fig)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "#fit simple linear regression model\n",
    "lm46 = ols('Purity ~ Hydrocarbon', data=dat46).fit()\n",
    "#create influence object containing studentized residuals\n",
    "influence = lm46.get_influence()\n",
    "#obtain studentized residuals, which are nearly equal to standardized residuals\n",
    "standardized_residuals = influence.resid_studentized_internal\n",
    "#display standardized residuals, not quite equal to our manual calculations above\n",
    "display(standardized_residuals)\n",
    "\n",
    "#define figure size\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "\n",
    "plt.scatter(X1, standardized_residuals)\n",
    "plt.xlabel('Hydrocarbon')\n",
    "plt.ylabel('Standardized Residuals')\n",
    "plt.axhline(y=0, color='black', linestyle='--', linewidth=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cad2d31-8bcb-47e3-9797-53d00057cf86",
   "metadata": {
    "tags": []
   },
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}