{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quasi-binomial regression\n",
    "\n",
    "This notebook demonstrates using custom variance functions and non-binary data\n",
    "with the quasi-binomial GLM family to perform a regression analysis using\n",
    "a dependent variable that is a proportion.\n",
    "\n",
    "The notebook uses the barley leaf blotch data that has been discussed in\n",
    "several textbooks. See below for one reference:\n",
    "\n",
    "https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_sect016.htm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The raw data, expressed as percentages.  We will divide by 100\n",
    "to obtain proportions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "raw = StringIO(\n",
    "    \"\"\"0.05,0.00,1.25,2.50,5.50,1.00,5.00,5.00,17.50\n",
    "0.00,0.05,1.25,0.50,1.00,5.00,0.10,10.00,25.00\n",
    "0.00,0.05,2.50,0.01,6.00,5.00,5.00,5.00,42.50\n",
    "0.10,0.30,16.60,3.00,1.10,5.00,5.00,5.00,50.00\n",
    "0.25,0.75,2.50,2.50,2.50,5.00,50.00,25.00,37.50\n",
    "0.05,0.30,2.50,0.01,8.00,5.00,10.00,75.00,95.00\n",
    "0.50,3.00,0.00,25.00,16.50,10.00,50.00,50.00,62.50\n",
    "1.30,7.50,20.00,55.00,29.50,5.00,25.00,75.00,95.00\n",
    "1.50,1.00,37.50,5.00,20.00,50.00,50.00,75.00,95.00\n",
    "1.50,12.70,26.25,40.00,43.50,75.00,75.00,75.00,95.00\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The regression model is a two-way additive model with\n",
    "site and variety effects.  The data are a full unreplicated\n",
    "design with 10 rows (sites) and 9 columns (varieties)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(raw, header=None)\n",
    "df = df.melt()\n",
    "df[\"site\"] = 1 + np.floor(df.index / 10).astype(int)\n",
    "df[\"variety\"] = 1 + (df.index % 10)\n",
    "df = df.rename(columns={\"value\": \"blotch\"})\n",
    "df = df.drop(\"variable\", axis=1)\n",
    "df[\"blotch\"] /= 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the quasi-binomial regression with the standard variance\n",
    "function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                 blotch   No. Observations:                   90\n",
      "Model:                            GLM   Df Residuals:                       72\n",
      "Model Family:                Binomial   Df Model:                           17\n",
      "Link Function:                  Logit   Scale:                        0.088778\n",
      "Method:                          IRLS   Log-Likelihood:                -20.791\n",
      "Date:                Fri, 20 Mar 2026   Deviance:                       6.1260\n",
      "Time:                        11:18:46   Pearson chi2:                     6.39\n",
      "No. Iterations:                    10   Pseudo R-squ. (CS):             0.3198\n",
      "Covariance Type:            nonrobust                                         \n",
      "==================================================================================\n",
      "                     coef    std err          z      P>|z|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------\n",
      "C(variety)[1]     -8.0546      1.422     -5.664      0.000     -10.842      -5.268\n",
      "C(variety)[2]     -7.9046      1.412     -5.599      0.000     -10.672      -5.138\n",
      "C(variety)[3]     -7.3652      1.384     -5.321      0.000     -10.078      -4.652\n",
      "C(variety)[4]     -7.0065      1.372     -5.109      0.000      -9.695      -4.318\n",
      "C(variety)[5]     -6.4399      1.357     -4.746      0.000      -9.100      -3.780\n",
      "C(variety)[6]     -5.6835      1.344     -4.230      0.000      -8.317      -3.050\n",
      "C(variety)[7]     -5.4841      1.341     -4.090      0.000      -8.112      -2.856\n",
      "C(variety)[8]     -4.7126      1.331     -3.539      0.000      -7.322      -2.103\n",
      "C(variety)[9]     -4.5546      1.330     -3.425      0.001      -7.161      -1.948\n",
      "C(variety)[10]    -3.8016      1.320     -2.881      0.004      -6.388      -1.215\n",
      "C(site)[T.2]       1.6391      1.443      1.136      0.256      -1.190       4.468\n",
      "C(site)[T.3]       3.3265      1.349      2.466      0.014       0.682       5.971\n",
      "C(site)[T.4]       3.5822      1.344      2.664      0.008       0.947       6.217\n",
      "C(site)[T.5]       3.5831      1.344      2.665      0.008       0.948       6.218\n",
      "C(site)[T.6]       3.8933      1.340      2.905      0.004       1.266       6.520\n",
      "C(site)[T.7]       4.7300      1.335      3.544      0.000       2.114       7.346\n",
      "C(site)[T.8]       5.5227      1.335      4.138      0.000       2.907       8.139\n",
      "C(site)[T.9]       6.7946      1.341      5.068      0.000       4.167       9.422\n",
      "==================================================================================\n"
     ]
    }
   ],
   "source": [
    "model1 = sm.GLM.from_formula(\n",
    "    \"blotch ~ 0 + C(variety) + C(site)\", family=sm.families.Binomial(), data=df\n",
    ")\n",
    "result1 = model1.fit(scale=\"X2\")\n",
    "print(result1.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below shows that the default variance function is\n",
    "not capturing the variance structure very well. Also note\n",
    "that the scale parameter estimate is quite small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
 
 
 
 
    },
    "lines_to_next_cell": 1
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/genmod/generalized_linear_model.py:985: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAGwCAYAAACq12GxAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcrZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvhF0PpwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVbZJREFUeJzt3XtcVGX+B/DPcBvEYBSIi4pAdhFkvbIqqdtlA8VSy6103chKLVNzjX61mluKWdZ20U3Xy5rpFl18VZZahE7bRVPMVLAUUzOM0kESFUgCRub8/mBnZJwLczln5pw5n/fr1a5z5jlnnnmYy3ee5/s8j0YQBAFEREREKhXk7woQERER+RODISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVWMwRERERKoW4u8KyJ3JZMLJkycRGRkJjUbj7+oQERGRCwRBQH19Pbp06YKgIOd9PwyG2nHy5EkkJSX5uxpERETkgZ9++gndunVzWobBUDsiIyMBtDZmVFSUJI9hNBqxdetW5OTkIDQ0VJLHCHRsQ3GwHcXBdvQe21Acam7Huro6JCUlWb7HnWEw1A7z0FhUVJSkwVBERASioqJU92IVC9tQHGxHcbAdvcc2FAfbES6luDCBmoiIiFSNwRARERGpmqKCoW3btmHUqFHo0qULNBoNPvjgg3bP+eKLLzBgwACEh4fjiiuuwMqVK6WvKBERESmGooKh8+fPo0+fPli2bJlL5SsqKjBy5EgMGzYMpaWlePzxxzFz5ky89957EteUiIiIlEJRCdS5ubnIzc11ufzKlSvRvXt3LFmyBACQlpaGPXv24IUXXsCf/vQniWpJRERESqKoYMhdJSUlyMnJsTo2fPhwrFmzBkaj0W5mfVNTE5qamiy36+rqALRm5BuNRknqab6uVNdXA7ahONiO4mA7eo9tKA41t6M7zzmgg6GqqirEx8dbHYuPj8eFCxdw+vRpJCYm2pyzaNEiFBQU2BzfunUrIiIiJKsrAOj1ekmvrwZsQ3GwHcXBdvQe21AcamzHhoYGl8sGdDAE2K4vIAiC3eNmc+bMQX5+vuW2edGmnJwcSdcZ0uv1yM7OVu06EN5iG4qD7SgOtqP32IbiUHM7mkd2XBHQwVBCQgKqqqqsjlVXVyMkJAQxMTF2z9FqtdBqtTbHQ0NDJX8h+eIxAh3bUBxsR3GwHb3HNhSHGtvRnecb0MFQVlYWNm/ebHVs69atyMzMVN2LgojkqcUkYHfFGVTXNyIuMhwDU6MRHMRNoYl8SVHB0K+//orvv//ecruiogJlZWWIjo5G9+7dMWfOHJw4cQKvvfYaAGDq1KlYtmwZ8vPzMWXKFJSUlGDNmjV46623/PUUiIgsig8YULC5HIbaRsuxRF045o1Kx4gM25xGIpKGotYZ2rNnD/r164d+/foBAPLz89GvXz88+eSTAACDwYDKykpL+dTUVBQVFeHzzz9H37598dRTT+Hll1/mtHoi8rviAwY8WLjPKhACgKraRjxYuA/FBwx+qhmR+iiqZ+j666+3JEDbs27dOptj1113Hfbt2ydhrYiI3NNiElCwuRz2Ps0EABoABZvLkZ2ewCEzIh9QVM8QEVEg2F1xxqZHqC0BgKG2EbsrzviuUkQqxmCIiMjHqusdB0KelCMi7zAYIiLysbjIcFHLEZF3GAwREfnYwNRoJOrC4SgbSIPWWWUDU6N9WS0i1WIwRETkY8FBGswblQ4ANgGR+fa8UelMnibyEQZDRER+MCIjESvu6o8EnfVQWIIuHCvu6s91hoh8SFFT64mIAsmIjERkpydwBWoiP2MwREQe41YS3gsO0iCrh/29EonINxgMEZFHuJUEEQUK5gwRkdu4lQQRBRIGQ0Tklva2kgBat5JoMTneOoeISE4YDBGRW7iVBBEFGgZDROQWbiVBRIGGwRARuYVbSRBRoGEwRERu4VYSRBRoGAwRkVu4lQQRBRoGQ0TkNm4lQUSBhIsuEpFHuJUEEQUKBkNE5DFuJUFEgYDDZERERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFSNwRARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFRNccHQ8uXLkZqaivDwcAwYMADbt293Wv6NN95Anz59EBERgcTERNx7772oqanxUW2JiIhI7hQVDK1fvx6zZs3C3LlzUVpaimHDhiE3NxeVlZV2y3/55Ze4++67MWnSJBw8eBDvvPMOvv76a0yePNnHNSciIiK5UlQw9NJLL2HSpEmYPHky0tLSsGTJEiQlJWHFihV2y+/atQspKSmYOXMmUlNTMXToUDzwwAPYs2ePj2tOREREchXi7wq4qrm5GXv37sXs2bOtjufk5GDnzp12z7n22msxd+5cFBUVITc3F9XV1Xj33Xdx8803O3ycpqYmNDU1WW7X1dUBAIxGI4xGowjPxJb5ulJdXw3YhuJgO4qD7eg9tqE41NyO7jxnjSAIgoR1Ec3JkyfRtWtX7NixA9dee63l+DPPPIP//Oc/OHz4sN3z3n33Xdx7771obGzEhQsXMHr0aLz77rsIDQ21W37+/PkoKCiwOf7mm28iIiJCnCdDREREkmpoaMCECRNQW1uLqKgop2UV0zNkptForG4LgmBzzKy8vBwzZ87Ek08+ieHDh8NgMODRRx/F1KlTsWbNGrvnzJkzB/n5+ZbbdXV1SEpKQk5OTruN6Smj0Qi9Xo/s7GyHQRo5xzYUB9tRHGxH77ENxaHmdjSP7LhCMcFQbGwsgoODUVVVZXW8uroa8fHxds9ZtGgRhgwZgkcffRQA0Lt3b3Ts2BHDhg3DwoULkZiYaHOOVquFVqu1OR4aGir5C8kXjxHo2IbiYDuKg+3oPbahONTYju48X8UkUIeFhWHAgAHQ6/VWx/V6vdWwWVsNDQ0ICrJ+isHBwQBae5SIiIiIFBMMAUB+fj5eeeUVvPrqqzh06BAefvhhVFZWYurUqQBah7juvvtuS/lRo0Zhw4YNWLFiBX744Qfs2LEDM2fOxMCBA9GlSxd/PQ0iIiKSEcUMkwHAuHHjUFNTgwULFsBgMCAjIwNFRUVITk4GABgMBqs1h+655x7U19dj2bJleOSRR9CpUyfceOONeO655/z1FIiIiEhmFBUMAcC0adMwbdo0u/etW7fO5thDDz2Ehx56SOJaERERkVIpapiMiIiISGwMhoiIiEjVGAwRERGRqjEYIiIiIlVjMERERESqxmCIiIiIVI3BEBEREama4tYZIiIi17SYBOyuOIPq+kbERYZjYGo0goPsb2xNpGYMhoiIAlDxAQMKNpfDUNtoOZaoC8e8UekYkWG7STWRmnGYjIgowBQfMODBwn1WgRAAVNU24sHCfSg+YPBTzYjkicEQEVEAaTEJKNhcDsHOfeZjBZvL0WKyV4JInRgMEREFkN0VZ2x6hNoSABhqG7G74ozvKkUkcwyGiIgCSHW940DIk3JEasBgiIgogMRFhotajkgNGAwREQWQganRSNSFw9EEeg1aZ5UNTI32ZbWIZI3BEBFRAAkO0mDeqHQAsAmIzLfnjUrnekNEbTAYIiIKMCMyErHirv5I0FkPhSXowrHirv5cZ4joElx0kYgoAI3ISER2egJXoCZyAYMhIqIAFRykQVaPGH9Xg0j2OExGREREqsZgiIiIiFSNwRARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFSNwRARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaopLhhavnw5UlNTER4ejgEDBmD79u1Oyzc1NWHu3LlITk6GVqtFjx498Oqrr/qotkRERCR3If6ugDvWr1+PWbNmYfny5RgyZAhWrVqF3NxclJeXo3v37nbPufPOO3Hq1CmsWbMGV155Jaqrq3HhwgUf15yIiIjkSlHB0EsvvYRJkyZh8uTJAIAlS5Zgy5YtWLFiBRYtWmRTvri4GF988QV++OEHREdHAwBSUlKcPkZTUxOampost+vq6gAARqMRRqNRpGdizXxdqa6vBmxDcbAdxcF2dE2LScCeH8+iur4JcZFaZCZ3RnCQBgDbUCxqbkd3nrNGEARBwrqIprm5GREREXjnnXdw2223WY7/9a9/RVlZGb744gubc6ZNm4YjR44gMzMTr7/+Ojp27IjRo0fjqaeeQocOHew+zvz581FQUGBz/M0330RERIR4T4iISMX212iw4XgQzjVrLMc6hQkYm2JCnxhFfC2RzDU0NGDChAmora1FVFSU07KK6Rk6ffo0WlpaEB8fb3U8Pj4eVVVVds/54Ycf8OWXXyI8PBzvv/8+Tp8+jWnTpuHMmTMO84bmzJmD/Px8y+26ujokJSUhJyen3cb0lNFohF6vR3Z2NkJDQyV5jEDHNhQH21EcbEfnthw8hbUl+3FpyFPbrMHaI8FYOr4Pbrw6mm0oAjW/Fs0jO65QTDBkptForG4LgmBzzMxkMkGj0eCNN96ATqcD0DrUdvvtt+Nf//qX3d4hrVYLrVZrczw0NFTyF5IvHiPQsQ3FwXYUB9vRVotJwNMfH7YJhABAAKAB8PTHh3FT2jAAbEOxqLEd3Xm+iplNFhsbi+DgYJteoOrqapveIrPExER07drVEggBQFpaGgRBwM8//yxpfYmIyNbuijMw1DY6vF8AYKhtxJ4fz/quUqR6igmGwsLCMGDAAOj1eqvjer0e1157rd1zhgwZgpMnT+LXX3+1HDty5AiCgoLQrVs3SetLRES2qusdB0LW5ZraL0QkEsUEQwCQn5+PV155Ba+++ioOHTqEhx9+GJWVlZg6dSqA1nyfu+++21J+woQJiImJwb333ovy8nJs27YNjz76KO677z6HCdRERCSduMhwF8vZpisQSUVROUPjxo1DTU0NFixYAIPBgIyMDBQVFSE5ORkAYDAYUFlZaSl/2WWXQa/X46GHHkJmZiZiYmJw5513YuHChf56CkREqjYwNRqJunBU1TbazRvSAEjQhSMzuTO2HPJ17UitFBUMAa3T5adNm2b3vnXr1tkc69mzp83QGhER+UdwkAbzRqXjwcJ90ABWAZF5Ksy8UemW9YaIfEFRw2RERKR8IzISseKu/kjQWQ+ZJejCseKu/hiRkeinmpFaKa5niIiIlG9ERiKy0xOwu+IMqusbERcZjoGp0ewRIr9gMERERH4RHKRBVo8Yf1eDiMNkREREpG4MhoiIiEjVGAwRERGRqjEYIiIiIlVjMERERESqxmCIiIiIVI3BEBEREakagyEiIiJSNQZDREREpGoMhoiIiEjVGAwRERGRqjEYIiIiIlVjMERERESqxmCIiIiIVI3BEBEREakagyEiIiJStRB/V4CIiMiZFpOA3RVnUF3fiLjIcAxMjUZwkMbf1bKhlHqSLQZDREQkW1sOnsLTHx+GobbRcixRF455o9IxIiPRjzWzVnzAgILN5bKvJ9nHYTIiIpKl/TUaPPT2fqsAAwCqahvxYOE+FB8w+Klm1ooPGPBg4T7Z15McYzBERESy02ISsOF4EAQ795mPFWwuR4vJXgnfaTEJKNhcLvt6knMMhoiISHb2/HgW55od59sIAAy1jdhdccZ3lbJjd8UZmx6htuRST3KOwRAREclOdX2Ti+UcByK+4Orj+7ue5ByDISIikp24SK2L5cIlrok4j+/vepJzDIaIiEh2MpM7o1OYAEcDZRq0ztYamBrty2rZGJgajURduOzrSc4xGCIiItkJDtJgbIoJAGwCDfPteaPS/b6OT3CQBvNGpQOQdz3JOQZDREQkS31iBCwd3wcJOushpgRdOFbc1V826/eMyEjEirv6y76e5BgXXSRSMK54S4FueK945PbuKvvX+YiMRGSnJ8i+nmQfgyEiheKKt6QWwUEaZPWI8Xc12qWUepItDpMRKRBXvCUiEg+DISKF4Yq3RETiUlwwtHz5cqSmpiI8PBwDBgzA9u3bXTpvx44dCAkJQd++faWtIJHEuOItEZG4FBUMrV+/HrNmzcLcuXNRWlqKYcOGITc3F5WVlU7Pq62txd13340//vGPPqopkXS44i0RkbgUFQy99NJLmDRpEiZPnoy0tDQsWbIESUlJWLFihdPzHnjgAUyYMAFZWVk+qimRdLjiLflKi0lAybEabCw7gZJjNRx6pYClmNlkzc3N2Lt3L2bPnm11PCcnBzt37nR43tq1a3Hs2DEUFhZi4cKF7T5OU1MTmpou7olTV1cHADAajTAajR7W3jnzdaW6vhqoqQ37dYtEQpQWp+qa7OYNaQAk6LTo1y3S7fZQUztKKRDaccvBU1hY9B2q6i5+HiZEafH3kT0xvFe85I8fCG0oB2puR3ees2KCodOnT6OlpQXx8dZvwvj4eFRVVdk95+jRo5g9eza2b9+OkBDXnuqiRYtQUFBgc3zr1q2IiIhwv+Ju0Ov1kl5fDdTShiMTNHi1ztyx23YdEwECgNz4Bmwp/tjj66ulHaWm1HbcX6PBq0dsX19VdY2Y8XYZ7rvahD4xvuklUmobyo0a27GhocHlsooJhsw0GusFrARBsDkGAC0tLZgwYQIKCgpw9dVXu3z9OXPmID8/33K7rq4OSUlJyMnJQVRUlOcVd8JoNEKv1yM7OxuhoaGSPEagU1sbjgTQ384v90RdOObmev7LXW3tKBUlt2OLScCiF7cBsLdrvAYaAB+fisBjf/mDpAsKKrkN5UTN7Wge2XGFYoKh2NhYBAcH2/QCVVdX2/QWAUB9fT327NmD0tJSzJgxAwBgMpkgCAJCQkKwdetW3HjjjTbnabVaaLW2uyWHhoZK/kLyxWMEOjW14S19u0m2Mq+a2lFKSmzHPcdqrALsS7XOVmxC6c/1PllgUIltKEdqbEd3nq9igqGwsDAMGDAAer0et912m+W4Xq/HmDFjbMpHRUXh22+/tTq2fPlyfPrpp3j33XeRmpoqeZ2JpMYVb0lsnK1IaqSYYAgA8vPzkZeXh8zMTGRlZeHf//43KisrMXXqVACtQ1wnTpzAa6+9hqCgIGRkZFidHxcXh/DwcJvjRERSUdr+cZytSGqkqGBo3LhxqKmpwYIFC2AwGJCRkYGioiIkJycDAAwGQ7trDhER+cqWg6fw9MeHfbZ/nBiB18DUaCTqwlFV2+hktmLrtYkChaKCIQCYNm0apk2bZve+devWOT13/vz5mD9/vviVIiK6xP4aDdaW7LcJKMz7x624q7+oAZFYG/cGB2kwb1Q6HizcBw1gVX9zWDVvVLqse7eI3KWoRReJiJSgxSRgw/Egn+0fJ/bGvSMyErHirv5I0FkPhSXowkUP4ojkQHE9Q0REcrfnx7M41+y456Tt/nHeJsC3t3GvBq2BV3Z6glu9OSMyEpGdnqCofCciTzEYIiISWXW946np1uW8n5Hlzsa97gZebWcrKi0RnMgdDIaIiEQWF2m7Vpn9ct7PyPLFVHix8pGI5Io5Q0REIstM7oxOYQIc9Zto0BpMiDEjS+qp8GLnIxHJEYMhIiKRBQdpMDbFBAA2AZHYM7LMU+GlCLzay0cCxE0EJ/IXBkNERBLoEyNg6fg+ks/IMk+FB8QPvNzJRyJSMuYMERFJZHiveMn2j2vLPBX+0ryeBC/zerg1B6kFgyEiIgn5av84KabCc2sOUgsGQ0REAULswItbc5BaMGeIiIjskjIfiUhOXO4Zevnll12+6MyZMz2qDBERyYtU+UhEcuJyMLR48WKXymk0GgZDREQBhFtzUKBzORiqqKiQsh5ERCRjvkoEJ/IH5gwRERGRqnk8m+znn3/Gpk2bUFlZiebmZqv7XnrpJa8rRkREROQLHgVD//3vfzF69Gikpqbi8OHDyMjIwPHjxyEIAvr37y92HYmIiIgk49Ew2Zw5c/DII4/gwIEDCA8Px3vvvYeffvoJ1113He644w6x60hEREQkGY+CoUOHDmHixIkAgJCQEPz222+47LLLsGDBAjz33HOiVpCIiIhISh4FQx07dkRTUxMAoEuXLjh27JjlvtOnT4tTMyKSnRaTgJJjNdhYdgIlx2q4WzkRBQSPcoYGDx6MHTt2ID09HTfffDMeeeQRfPvtt9iwYQMGDx4sdh2JSAaKDxhsFt5L5MJ7RBQAPAqGXnrpJfz6668AgPnz5+PXX3/F+vXrceWVV7q8OCMRKUfxAQMeLNxnsz9VVW0jHizchxV39WdARESK5VEwdMUVV1j+HRERgeXLl4tWISLyjxaTgK8qzmDvaQ1iKs4g68o4BAdp0GISULC53O5GnQJa96gq2FyO7PQErkhMRIrEXeuJ6JIhsGC8dnSPZQhM1yHMamjsUgIAQ20jdlec4QrFRKRIHgVDQUFB0Ggc/wJsaWnxuEJE5FvtDYHdOyTFpetU1zsOmIiI5MyjYOj999+3um00GlFaWor//Oc/KCgoEKViRCQ9V4bANpaddOlacZHhYlaNiMhnPAqGxowZY3Ps9ttvR69evbB+/XpMmjTJ64oRkfR2V5xpdwis5nwzojuG4ux5o92gSQMgQde6izkRkRKJulHroEGD8Mknn4h5SSKSkKtDW7f17QqgNfBpy3x73qh0Jk8TkWKJFgz99ttvWLp0Kbp16ybWJYlIYq4Obd2UnoAVd/VHgs66fIIunNPqiUjxPBom69y5s1UCtSAIqK+vR0REBAoLC0WrHBFJa2BqNBJ14aiqbWx3CCw4SIPs9ATsrjiD6vpGxEVePE5EpGQeBUOLFy+2CoaCgoJw+eWXY9CgQejcubNolSMiaQUHaTBvVDoeLNwHDWAVENkbAgsO0nD6PBEFHI+CoXvuuUfkahCRv4zISMSKu/rbbLWRwK02iEglXA6GvvnmG5cv2rt3b48qQ0T+MSIjEdnpCSj5vhpbt3+FnGGDLCtQExEFOpeDob59+0Kj0UAQWjvS/bXo4vLly/H888/DYDCgV69eWLJkCYYNG2a37IYNG7BixQqUlZWhqakJvXr1wvz58zF8+HDJ6kekVMFBGgxKjUbNIQGDmAtEFFBaTAL2HKthvp8DLs8mq6iowA8//ICKigps2LABqampWL58OUpLS1FaWorly5ejR48eeO+99ySr7Pr16zFr1izMnTsXpaWlGDZsGHJzc1FZWWm3/LZt25CdnY2ioiLs3bsXN9xwA0aNGoXS0lLJ6khERCQn+2s0uP7Fbfjz6l3469tl+PPqXRj63KcoPmDwd9Vkw+WeoeTkZMu/77jjDrz88ssYOXKk5Vjv3r2RlJSEJ554ArfeequolTR76aWXMGnSJEyePBkAsGTJEmzZsgUrVqzAokWLbMovWbLE6vYzzzyDjRs3YvPmzejXr58kdSQiIpKLLQdP4dUjQQCarI6bt9vh0hitPEqg/vbbb5GammpzPDU1FeXl5V5Xyp7m5mbs3bsXs2fPtjqek5ODnTt3unQNk8mE+vp6REc7Xim3qakJTU0XXzR1dXUAWrccMRqNHtS8febrSnV9NWAbioPtKA62o/fYht5rMQl4qug7u/eZt9sp2HwQ118VE5BDZu68djwKhtLS0rBw4UKsWbMG4eGti7A1NTVh4cKFSEtL8+SS7Tp9+jRaWloQHx9vdTw+Ph5VVVUuXePFF1/E+fPnceeddzoss2jRIrv7q23duhURERHuVdpNer1e0uurAdtQHGxHcbAdvcc29NzRWg1O1QXDdu34VgIAQ20Tlq0vxlU6eyuNKVtDQ4PLZT0KhlauXIlRo0YhKSkJffr0AQDs378fGo0GH374oSeXdNmliduCIDhN5jZ76623MH/+fGzcuBFxcXEOy82ZMwf5+fmW23V1dUhKSkJOTg6ioqI8r7gTRqMRer0e2dnZCA0NleQxAh3bUBxsR3GwHb3HNvTe5m8MQPm37Za7oldfjOwdeENl5pEdV3gUDA0cOBAVFRUoLCzEd999B0EQMG7cOEyYMAEdO3b05JLtio2NRXBwsE0vUHV1tU1v0aXMm8e+8847uOmmm5yW1Wq10Gq1NsdDQ0Mlf0P64jECHdtQHGxHcbAdvcc29FxiJ9e+jxM7dQzINnbnOXkUDAFAREQE7r//fk9Pd1tYWBgGDBgAvV6P2267zXJcr9djzJgxDs976623cN999+Gtt97CzTff7IuqEhER+d3A1GgkRGlRVdcIe0NlbbfbUTuXg6FNmzYhNzcXoaGh2LRpk9Oyo0eP9rpi9uTn5yMvLw+ZmZnIysrCv//9b1RWVmLq1KkAWoe4Tpw4gddeew1AayB0991345///CcGDx5s6VXq0KEDdDqdJHUkIiKSg+AgDf4+sidmvF3mdLsdAChR+RpELgdDt956K6qqqhAXF+d06rxGo5Fs0cVx48ahpqYGCxYsgMFgQEZGBoqKiizT/g0Gg9WaQ6tWrcKFCxcwffp0TJ8+3XJ84sSJWLdunSR1JCIikovhveJx39UmFFVFoKru4kxp83Y7ADD0uU+ttuJJVOFWPC4HQyaTye6/fW3atGmYNm2a3fsuDXA+//xz6StERJJqMQnYXXFG1b9aibzRJ0bAY3/5A0p/rrd6H+nLq/Bg4T5cOo9MjWsQeZwzdKlz586hU6dOYl2OiAjFBww2G8iq8VcrkbeCgzTI6hFjud1iElCwudwmEALarkFUjuz0BFX8+HB5O462nnvuOaxfv95y+4477kB0dDS6du2K/fv3i1Y5IlKv4gMGPFi4zyoQAi7+auVWAkSe211xxua91VbrGkSN2F1xxneV8iOPgqFVq1YhKSkJQOtsrk8++QTFxcXIzc3Fo48+KmoFiUh92vvVCrT+am0xBd5CcUS+UF3vOBDypJzSeTRMZjAYLMHQhx9+iDvvvBM5OTlISUnBoEGDRK0gETkWqPk07vxqbdv1T74TqK89tYiLDBe1nNJ5FAx17twZP/30E5KSklBcXIyFCxcCaF0NWqqZZERkLZDzafirVd4C+bVnFujB3sDUaCTqwlFV22i3B1ZtaxB5FAyNHTsWEyZMwFVXXYWamhrk5uYCAMrKynDllVeKWkEismXOpwnUWSD81Spfgf7aA9QR7AUHaTBvVDoeLNzndA2iQAoAnfEoZ2jx4sWYMWMG0tPTodfrcdlllwFoHT5zNO2diMShtHyaFpOAkmM12Fh2AiXHalyql/lXq6OPYQ1av5zU8qvVU560fXvXU9JrzxNqStwfkZGIFXf1R4LO+kdFgi48IIJad3jUMxQaGor/+7//szk+a9Ysb+tDRO1QUj6Np7+w+avVe1L0bijptecJNU43H5GRiOz0hIAeEnSFRz1DAPD6669j6NCh6NKlC3788UcAwJIlS7Bx40bRKkdEtpSST+PtL2z+avWcVL0bSnnteUqt083NaxCN6dsVWT1iVBcIAR72DK1YsQJPPvkkZs2ahaefftqSNN2pUycsWbLE6capROQdJeTTiPULOxB+tfo6EVfK3g2xXntyTU4O9GCPHPMoGFq6dClWr16NW2+9Fc8++6zleGZmpt3hMyISjxJmgYg5nHLpyrlKsuXgKTz98WGfJuJKOZQlxmtPzsnJSvihQdLwaJisoqIC/fr1szmu1Wpx/vx5rytFRI6Z82kA2CQYyyWfhr+wgf01Gjz09n6fJ+JK2fbevvbknpzMxH318igYSk1NRVlZmc3xjz/+GGlpad7WiYjaIfd8GrX/wm4xCdhwPMgvs66kbntPX3tKmImmhB8aJA2PhskeffRRTJ8+HY2NjRAEAbt378Zbb72FZ555BmvWrBG7jkRkh5zzaZQwlCelPT+exblmx38HKWdd+aLtPXntKWUmmjnYu3QoL0EmQ3kkDY+CoXvvvRcXLlzAY489hoaGBkyYMAFdu3bF0qVLMWzYMLHrSEQOyDWfxpWp8U/cnCbLQE4M1fVNLpYTf5jQV8sSuPvaU9LQqZx/aJA0PAqGAGDKlCmYMmUKTp8+DZPJhJaWFjzzzDOYPn06fvvtNzHrSEQK5OwX9ug+iXjqo0OyTKJ1RXuzoeIitS5dR6phQjn2biht6FSuPzRIGm4FQ+fOncP06dOxdetWhIaGYvbs2ZgxYwYKCgrwwgsvID09Ha+++qpUdSUihbH3C/vs+WZMf1O52zm4MhsqM7kzOoUJqG3W+G2YUG69G4EwdCrXJQHIe24FQ48//ji2bduGiRMnori4GA8//DCKi4vR2NiIoqIiXHfddVLVk4gUqu0v7BaTgKHPfarYFX5d3ZcrOEiDsSkmrD0S7NcVtOXUu6H0VcXlvCQAec+t2WQfffQR1q5dixdeeAGbNm2CIAi4+uqr8emnnzIQIgpQYu5vpeQVft2dDdUnRsDS8X1kO+PPH+Q+C9IRuS8JQN5zq2fo5MmTSE9vnXZ4xRVXIDw8HJMnT5akYkTkG227/mMiQtA21hH717CSkmgv5U4gl9k9CgAwvFc8cnt35dBKG3IbvmuPGvcrUyO3giGTyYTQ0FDL7eDgYHTs2FH0ShGRb9gLdjqFBSM05RRCQoJdGhJyhxKSaB3lhbgXyEVZbstpqEoulNQmSlkSgLzjVjAkCALuueceaLWtMyUaGxsxdepUm4Bow4YN4tWQiCThKP/lXDPw0Nv7oYsIFf3XsNyTaJ31hCkhkCPxuRoE68urGAwpmFs5QxMnTkRcXBx0Oh10Oh3uuusudOnSxXLb/B8RecdRno5Y+TvOuv6B1hlQ5xqMDs/3NLdHziv8tpcXcvZ8M7dqUCFXg9tXdxxn7pCCudUztHbtWqnqQUT/U/TNSfx94wGcOX8xGEn839o8m/YbRMnfaa/r31We5PbIcQ0cV/JCnvqoHE/cnIbpb5a2OxvK1CJ1jclXzL2Z7b1fmDukbB4vukhE4ltUVI5V2ypsjhtqG+0e9zR/R6wEZW/2t2ovidaXa7q4mhfSuaNWdoEcScvcmzm1cJ/TcswdUjYGQ0QyUfSNwW7A44yn+Tve5rWIkdvjLInW12u6uJMcPaZvV0XNhiLvjchIxKQhKViz43i7ZeU4E5La59Gu9UQkrhaTgL9vPODRuZ7k75i7/p3lv3SOCLX8+9L7AOlye/yxpou7ydHmQG5M367I6hHDQEgFbkpPcKkcE+iVicEQ+Y2Yi/kp3e6KMzhzvtmra7jzi9RZIrM5G2bR2N9hpY8XyHN3YUOxuBIcMjla3fgaCWwcJiO/4NL21sToWnf3F6mjROZOYcDCsX0sfwdfDgn5a00XpW8VQdLjaySwMRgin3N1fyd/8vWGjN52rXv6i/TSROaYiBD8Ur4Lw3vFW8r4coE8f65QLcdZbiQvfI0ELgZDfuJsC4RApoSl7f3Ra+Xq9F1Hnrg5zeP2ahvsGI1GFB3y6DKi8PfChkrbKoJ8j6+RwMRgyA+cbYFwS99ufqyZLbF7SOS+tL2/eq3adsF7Ehd37qgVvU7+IIcVqpW0VYQrfN3LqQaB9hohBkM+194WCCEhwbLpapWih0TOG3X6u9fKURe8KwJlOi/zMlzjaoDD3Dwi1yhuNtny5cuRmpqK8PBwDBgwANu3b3da/osvvsCAAQMQHh6OK664AitXrvRRTW21twUCIM1MGU9INb3Z38MgzrjTayWVERmJ+PJvN+KtKYPxz/F98cTNaS6dF0jTec1BoS9nsSlJ8QEDhj73Kf68ehf++nYZ/rx6F4Y+96nNe9IfSxQQKZWieobWr1+PWbNmYfny5RgyZAhWrVqF3NxclJeXo3v37jblKyoqMHLkSEyZMgWFhYXYsWMHpk2bhssvvxx/+tOffF5/uQ8RmUnZQyKHYRBH/Nlr5eiXfotJwCtfVsiyvdzlznAN8zLsc3UY19+9nESukNMQrqKCoZdeegmTJk3C5MmTAQBLlizBli1bsGLFCixatMim/MqVK9G9e3csWbIEAJCWloY9e/bghRde8EswJOchorakDNrkPAzir16r9oYy5Npe7vBkuIZ5GdbcCXCU8sOL1EtuQ7iKCYaam5uxd+9ezJ492+p4Tk4Odu7cafeckpIS5OTkWB0bPnw41qxZA6PRiNDQUJtzmpqa0NTUZLldV1cHoHWWjdHoeBdvV8REuNbcp2obsGFvJeIitchM7uzzLzrDufMulzMao9y+/h+vicXS8X2wsOg7VNVdbOsEnRZzc3vij9fEut3W5vLe/I36dYtEQpQWp+qanPTCaNGvW6TXrwWzLQdP4aG39zv8pb90fB8M7xUvens5IkY7XsrV5xhIpGjHr1wMcEq+r0Z1fZPDcm15+h72BaPRCJMA7DhajTO/tfjt81DppHgtestXnwnuPGfFBEOnT59GS0sL4uOtGyg+Ph5VVVV2z6mqqrJb/sKFCzh9+jQSE22jz0WLFqGgoMDm+NatWxEREeHFMwBMQuussXPNgL11fwEBGgDPfHzEcqRTmICxKSb0ifFdHtEPtRoAwe2XO1iGop9LPX6cv6UDx+o0qDMCUaFAj6jzaPlxL4p+9PiS0Ov1np8MYGSCBq/WmVPp2v6NBAgA+kY24JnC4v/VV4A3n8smASjYF/y/DwTrCwn/+9+/byiD8XgLgjTStJcj3rajmbvPMdCI1Y4AsPe0a+/Lrdu/QlQoXCrr7XvYEybh0tex/ffR/hoNNhwPxrldZZZj/vg8DBRivha94cvPhIaGBpfLKiYYMtNoLmk8QbA51l55e8fN5syZg/z8fMvturo6JCUlIScnB1FR3v+CCk1pjYgBXBIVt3Z0X/oWr23WYO2RYK8j5RaTgD0/nkV1fVO7v7BaTALefXFbuz0kM8b9QTa/0oxGI/R6PbKzs+32+LlqJID+B0/Z9MJ0iggFBA0+/vniL42EKC3+PrKnx3+XryrO4NyuPU5KaHCuGbg8fTAGuZkT5M7fuy2x2tFMyucoZ2K3IwDEVJzBa0edtWWrnGGDkJncWZbv4S0HT2HRpT2cdt5HWw6ewtqS/f/7CXKRWJ+HaiLFa9EbvvxMMI/suEIxwVBsbCyCg4NteoGqq6tten/MEhIS7JYPCQlBTIz9cXKtVgut1nbNltDQUFFeSLf07YaQkGCbsdJL80HMzLkAT398GLm9u3r0weXu2GwogPmje7WTp9IL4dowt+viDk+S68T4O93Stxtye3e1PPbx0w1Y8skRm7/PqbomPPT2fo9nONU0XHC5nDvPSYyxeLFe71I9R6UQqx0BIOvKOJcmH2RdGYfgII0s3sNtFR8w2B0aufR91GIS8PTHhx32HHj7eahWYr4WveHLzwR3zlfM1PqwsDAMGDDApqtPr9fj2muvtXtOVlaWTfmtW7ciMzPTry+KS6dPP5579f8GyOzzZkq3p9Nr/T292dXpw1IxJ+/e0rsL3v66UpKNQ6VI2JbbdGo5L6WgNM4217WXTO/v93Bb7mzAK4clLkg6cv1MUEzPEADk5+cjLy8PmZmZyMrKwr///W9UVlZi6tSpAFqHuE6cOIHXXnsNADB16lQsW7YM+fn5mDJlCkpKSrBmzRq89dZb/nwaAKxnymzYW+nSOe7OMvN2eq0vpjfb6/3Rl1fJZu8yKWfliL3MgBynU8t5KQUlcndvLLksUeDO+0gps27JM3L9TFBUMDRu3DjU1NRgwYIFMBgMyMjIQFFREZKTkwEABoMBlZUXA4vU1FQUFRXh4Ycfxr/+9S906dIFL7/8sl+m1TsTF+naVgruRspifJFLOb3Z3nBOQpQWjRdMsvlCl/KDWexlBlz9e6/bUYHYSK3bX4yeDFvKeSkFpXI3wJHDEgXuvI/k2nNA4pDrZ4KigiEAmDZtGqZNm2b3vnXr1tkcu+6667Bv3z6Ja+WdzOTO6BQmoLbZNoEa8DxSlvMvLIeLx9U5nxIsxvooLSYBu36oQcmxGgACsq6IxeAeMXbffFJ/MIu5C7arf8enPrq4E6uruUTe5CFxp2/xySHAcYc77yO59hyQeOT4maC4YCgQBQdpMDbFhLVHgkWNlOX6C8v5tiSu8TSAKz5gwOwN3+Jcw8VZYcs+O4ZOEaF4duzvbN6EvvhgFmsow5O/oytDj2JsXiuX4RryD3feR217Di72B18sB7A3MRDI7TNBMQnUga5PjICl4/uImuxo/gBy9NLSoPXXva9/YbU3nOMKT774iw8YMLVwn1UgZHauwYipdhKM3U1a9ZT5l/6Yvl2R5aCXqj3t/b3taS8J3J3E1/aI8RxJmTxJ/l46vg86XTLZjfvTBRY5fSawZ0hGhveKt5rS7W2kLNexWW+G5TztiWkxCZi/qbzdcvM3HbTJR5Jjl649zv7ezrQdeszsbr2WFrd1ILG4+z4a3isexuMtuDx9MGoaLvi954ACG4MhmRE7F0COX+SeDst5E8DtrjiDqrr2g7Cquia7X+xy69J1xNHf2xWtQWqUnWOunitvctoUUq3cfR8FaYBBqdGyWB+HAhuDIRWQ2xe5K/kDuohQhIcEWwUw3gRw7nxZOyqrlKTVERmJuLFnPF4vOY4fzzRAEAS8vqv95RvsBalyzTtzl9w2hVQzpbyPSF0YDKmEnD6AXBm+e3bs70QN4Nz5spb7F3t77H3xB2la9wSyp+3Qo6nFenXYQJjZI0YCeCBjjxkRgyHyE1eH78QK4AamRiMhKrzdobKEKK2sv9jb4+iL31kgBFwcejS1WN8v17wzV8lxIUo5YY8ZUSsGQ+Q3vhy+a92rKR1TC52vOTV/dC/Ffim6smTBpT1Ergw9+ivvTIweC1cTwHf9UIMhV8Z6WWNlUWqPGXuySAoMhsivfDl8NyIjESvv6m+zzhAAh+sMKYkrSxaYBOCJm9PcXoHa13lnYvVYuJordv9re/D87b0xsncXt+uqRP7sMfMmmGFPFkmFwRCpivlL3dUVqJXE1S/+2EgtxvTt6vb1fRW4itlj4Wr+1/nmFkx7sxQP/HwOc0amu1Vfe1/ucuevJRO8CWaU2pNFysBgiFQnOEiDIVfGBtywSCDM/BK7x6K9BPBLrdpWgT7dOmNkb9e+VB19uc/Nvcal8/3FH0smeBPMMPeLpMYVqIkChFxXHHeHOz0Wrmi78rGrnth4wKUVtc1f7pfWt6q2EQ+9vR/7a+T7pezrwNnblczFfl0QXYrBEFGA8HTrkBaTgJJjNdhYdgJfVZxxOPPMF6TosTAngHfq4NrCfTXnm9v9UnXly33D8SCXgip/8HXg7G0wE0iLf5I8cZiMKIC4O/PL3jBPp7BghKacwi19u/ms3mZS9ViMyEhEpDYUf1nzlUvl2/tSdeXL/VyzBnt+PIuhV8e7U1WPuZOY7OslE7wNZgJhCJjkjcEQUYBxdeaXoxyOc83AQ2/vR0hIsKgJqa58WUu5yOPgHjGI7hiKM+dtN+q9VHtfqq5/uTe5VM5bniQm+3LJBG+DmUBY/JPkjcEQUQBqb+aX8zWJWgMUMRNSXf2ylrLHIjhIg4VjMjDtzVKn5VwZHnL9y13rcv085U1isq+WTPA2mFH64p8kf8wZIlIhXyakOks0frBwH4oPGKyOm3ssEnTWAUeCLtzr6dMje3fBA39IdXi/Bq59qbqSc9MpTEBmcmeP6+oKbxOTgYuB85i+XZEl0RITnuaztSXl64KIPUNEKuSrhFRPp0S722PRfMFk2Zg2OToCeVkpCAux/1tvzsh09OnWGX/feABnzjdbjruzeJ8rPRVjU0yS91T4a70gT4gxLCe3TacpcDAYIlIhXyWkevNl7eoij4uKyrF6e4XVLLiniw5hyrBUhwsojuydiOEZ3n2pOvtyn5t7DVp+3OvytTyltFlWYgQzctp0mgIHgyEiFfJVQqrUX9aLisqxaluFzXGTAMtxRwGRGF+qjr7cTS0XUPSjV5d2iRJnWTGYITlizhCRCjnL4TAP+oiRkCrll3XzBRNWb7cNhNpavb0CzRdMbl/bHb7IuXEkEBbaJJIDBkNEKuUoIbVTGLB0fB+7ORxtF2gsOVbT7qKCUn5Zv15yvN0FIk1CazlvuPucfUmMxGQi4jAZkapdOswTExGCX8p3YXgv24UCPVnLRsop0T+eaRC1nD1K2CXdl+sFEQUqBkNEKtc2h8NoNKLokG0Zb9eykeLLOjk6QtRyl1LSLumcZUXkHQZDROSUGDuGS/FlnZeVgqeLDjkdKgvStJZzlxJ3SWdiMpHnmDNERE6JtUCj2InGYSFBmDLM8QKKADBlWKrD9YacUcMu6XLOhSLyNfYMEZFTcl7Lxjxt/tJ1hoI0cLrOUHvk/JzFoIRcKCJfYjBERE7JfS2bOSPT8UhOT5dXoHaF3J+zN5SUC0XkKwyGiMgpJewYHhYShEnDrhDtekp4zp5QYi4UkS8wZ4iInFLjWjb+fM5S5vKoIReKyBPsGSKidslxLZsWkyDpVHJ/PGepc3kCPReKyFMMhojIJXJay8ZXCcC+fM6+yOUJ5FwoIm8wGCIil/liLZv2enx8nQDsq+fsi1yeQM2FIvIWgyEiko32enwCNQHYnVwebwIzKbdHIVIyxSRQnz17Fnl5edDpdNDpdMjLy8O5c+ccljcajfjb3/6G3/3ud+jYsSO6dOmCu+++GydPnvRdpYnIZeYen0uDAnOPT/EBQ8AmAPsyl8fRBr0JunBOqyfVUkzP0IQJE/Dzzz+juLgYAHD//fcjLy8Pmzdvtlu+oaEB+/btwxNPPIE+ffrg7NmzmDVrFkaPHo09e/b4supE1A5Xe3weG36NS9dTWgKwr3N55JT/RSQHigiGDh06hOLiYuzatQuDBg0CAKxevRpZWVk4fPgwrrnG9gNSp9NBr9dbHVu6dCkGDhyIyspKdO/e3e5jNTU1oampyXK7rq4OQGtPk9FoFOspWTFfV6rrqwHbUBz+asevXOzx+aX+N5euFxMR4tVzaDEJ2PPjWVTXNyEuUovM5M5uBQrutmO/bpFIiNLiVF2Tk1weLfp1ixT1b5PZPQpAFADA1HIBphbRLu01vqfFoeZ2dOc5KyIYKikpgU6nswRCADB48GDodDrs3LnTbjBkT21tLTQaDTp16uSwzKJFi1BQUGBzfOvWrYiI8Gz3a1ddGryR+9iG4vB1O+49rQEQ3G65n44eQqewIJxrBmxXAAIAAZ3CgF/Kd6HokGd12V+jwYbjQTjXfPH6ncIEjE0xoU+Me2v+uNOOIxM0eLXOnLnQ9rkJEADkxjdgS/HHbj1+IOB7WhxqbMeGhgaXyyoiGKqqqkJcXJzN8bi4OFRVVbl0jcbGRsyePRsTJkxAVFSUw3Jz5sxBfn6+5XZdXR2SkpKQk5Pj9DxvGI1G6PV6ZGdnIzQ0VJLHCHRsQ3H4qx1jKs7gtaPtD18P/8MgDGow4qG39wOwlwCswcKxfTC8V7xH9dhy8BTWluy36Z2pbdZg7ZFgLB3v2rU9aceRAPofPIWFRd+hqu5i73SiLhxzc3t6/JyUiu9pcai5Hc0jO67wazA0f/58u70wbX399dcAAI3G9legIAh2j1/KaDRi/PjxMJlMWL58udOyWq0WWq3W5nhoaKjkLyRfPEagYxuKw9ftmHVlnEtTvrOujENwkAYhIcGiL4bYYhLw9MeHneYtPf3xYeT27urykJm77XhL327I7d2VuTxt8D0tDjW2ozvP16/B0IwZMzB+/HinZVJSUvDNN9/g1KlTNvf98ssviI93/mvJaDTizjvvREVFBT799FPJeneIyHPuTvmWIgHYV9Pb2+OLdY2IyJpfg6HY2FjExsa2Wy4rKwu1tbXYvXs3Bg4cCAD46quvUFtbi2uvvdbheeZA6OjRo/jss88QE8MPGCK5cnf7C7GDBm5VQaReisgZSktLw4gRIzBlyhSsWrUKQOvU+ltuucUqebpnz55YtGgRbrvtNly4cAG333479u3bhw8//BAtLS2W/KLo6GiEhYX55bkQkWP+nPLNrSqI1EsRwRAAvPHGG5g5cyZycnIAAKNHj8ayZcusyhw+fBi1tbUAgJ9//hmbNm0CAPTt29eq3GeffYbrr79e8joTkfv8NUzErSq8I/XGuURSUkwwFB0djcLCQqdlBOHiR1hKSorVbSIiZ6TeqiKQgwVfbZxLJBXFBENERFJzN2/JVYEcLPh641wiKTAYIiJqQ+y8pUAOFgJ141xSHwZDRESXECtvqb1gAVB2sCCX5QiIvKWYXeuJiJRmz49nnQYLQGuwsOzToz6qkbi4HAEFCvYMERFJpLq+qf1CABZ/chTXJEQqbriMyxH4XiAn4vsTgyEiIonERdpu7eOIEofLuByBbwVyIr6/cZiMiEgimcmdkahzrVfEnFujJOblCICLyw+YibEcAV1kTsS/dNjVnIhffMDgp5oFBgZDREQSaRssuEKJuTXm5QgSLgn6EnThip4pJyeuJuK3mLi2nqc4TEZEJKERGYl4+KarsfiTI+2WVWpujT+3UVEDztqTHoMhIiKJzbjxSry1+0dU1dlPqA6E3Bp/baOiBpy1Jz0OkxERSSw4SIP5o3tBA+bWkPs4a096DIaIiHyAuTXkKfOsPUehsgats8qU3LPobxwmIyLyEebWkCek3kSYGAwREfkUc2vIE1JtIkytGAwREREpAHsWpcNgiIiISCHYsygNJlATERGRqjEYIiIiIlXjMBkRkZu4czhRYGEwRETkBu4cThR4OExGROQiRzuHG2obMbVwH/75yRFJNstsMQkoOVaDjWUnUHKshhtyEomMPUNERC5wtnO42eJPjuKt3T/h7yOvEe1x2RNFJD32DBERuaC9ncPNquoa8dDb+7G/xvscIkc9UVW1jXiwcB+KDxi8fgwiYjBEROQSd3cE33A8yKvhLGc9UeZjBZvLOWRGJAIGQ0RELnBnR3ABwLlmDfb8eNbjx2uvJ0pAa67S7oozHj8GEbViMERE5IL2dg63p7q+yePHc7Unyt0eKyKyxWCIiMgF5p3D3REXqfX48VztiXKnx4qI7GMwRETkIvPO4QlRzoMcDYBOYQIykzt7/Fjt9URp0DqrbGBqtMePQUStGAwREblhREYidsz+Ix6+6Wq795uDl7EpJq9WpW7bE3XpVcy3541K58rXRCJgMERE5KbgIA3+etNVWHlXfyTqrIepEnThWDq+D/rEeD/Ly9ITZecxVtzVn+sMEYmEiy4SEXloREYistMTbPYpM7VcQNGP0j4Ge4SIxMNgiIjIC8FBGmT1iLE6ZmqR/jG8xc1mfY9tLl8MhoiIVIZbfPge21zeFJMzdPbsWeTl5UGn00Gn0yEvLw/nzp1z+fwHHngAGo0GS5YskayORERyxy0+fI9tLn+KCYYmTJiAsrIyFBcXo7i4GGVlZcjLy3Pp3A8++ABfffUVunTpInEtiYjki1t8+B7bXBkUEQwdOnQIxcXFeOWVV5CVlYWsrCysXr0aH374IQ4fPuz03BMnTmDGjBl44403EBoa6qMaExHJD7f48D22uTIoImeopKQEOp0OgwYNshwbPHgwdDoddu7ciWuuucbueSaTCXl5eXj00UfRq1cvlx6rqakJTU0Xl9Cvq6sDABiNRhiNRi+ehWPm60p1fTVgG4qD7SgOubaj4dx5l8sZjVES18Y5ubahu/zd5oHSjp5w5zkrIhiqqqpCXFyczfG4uDhUVVU5PO+5555DSEgIZs6c6fJjLVq0CAUFBTbHt27dioiICJev4wm9Xi/p9dWAbSgOtqM45NaOP9RqAAS3X+5gGYp+LpW+Qi6QWxu6Sy5trvR29ERDQ4PLZf0aDM2fP99u4NHW119/DQDQaGynHwqCYPc4AOzduxf//Oc/sW/fPodl7JkzZw7y8/Mtt+vq6pCUlIScnBxERUnzS8loNEKv1yM7O5tDeR5iG4qD7SgOubZji0nAuy9uw6m6Jrs5LBoACTotZoz7g9+nfMu1Dd3l7zYPlHb0hHlkxxV+DYZmzJiB8ePHOy2TkpKCb775BqdOnbK575dffkF8fLzd87Zv347q6mp0797dcqylpQWPPPIIlixZguPHj9s9T6vVQqu13XcoNDRU8heSLx4j0LENxcF2FIfc2jEUwPzRvfBg4T5oAKsv54tbfPRCuDbM95VzwN02lNtaPnJpc7m9Fn3Bnefr12AoNjYWsbGx7ZbLyspCbW0tdu/ejYEDBwIAvvrqK9TW1uLaa6+1e05eXh5uuukmq2PDhw9HXl4e7r33Xu8rT0SkQOYtPi5d8yYhANa8ketaPoHc5oFCETlDaWlpGDFiBKZMmYJVq1YBAO6//37ccsstVsnTPXv2xKJFi3DbbbchJiYGMTHWK7aGhoYiISHBYcI1EZEaBOIWH+a1fC4dijKv5ePvvdwCsc0DiSKCIQB44403MHPmTOTk5AAARo8ejWXLllmVOXz4MGpra/1RPSIiRZFiiw9/aW8tHw1a1/LJTk/wa/ARSG0eaBQTDEVHR6OwsNBpGUFwvmiVozwhIiJSLnfW8mEwQvYoYtFFIiIiR6rrHQdCnpQj9WEwREREihYXGS5qOVIfBkNERKRoA1OjkagLh6NsIA1aZ5UNTI32ZbVIQRgMERGRogUHaTBvVDoA2AREF9fySefMLXKIwRARESmeeS2fBJ31UFiCLtzv0+pJ/hQzm4yIiMgZruVDnmIwREREAYNr+ZAnOExGREREqsZgiIiIiFSNwRARERGpGnOGiIhUqMUkMNGY6H8YDBERqUzxAQMKNpdb7eeVqAvHvFHpnIJOqsRhMiIiFSk+YMCDhftsNjatqm3Eg4X7UHzA4KeaWWsxCThaq8HmbwwoOVaDFpPzjbiJvMGeISIilWgxCSjYXA57YYWA1tWaCzaXIzs9wa9DZsUHDJi/6SCq6oKB8m8BsOeKpMWeISIildhdccamR6gtAYChthG7K874rlKXMPdcVdU1WR2XW88VBRYGQ0REKlFd7zgQ8qSc2NrruQJae644ZEZiYzBERKQScZHh7Rdyo5zYlNBzRYGJwRARkUoMTI1Goi7cZmd3Mw1ac3MGpkb7sloWcu+5osDFYIiISCWCgzSYNyodAGwCIvPteaPS/ZY8LfeeKwpcDIaIiFRkREYiVtzVHwk664AiQReOFXf19+tsLbn3XFHg4tR6IiKVGZGRiOz0BNmtQG3uuXqwcB80gFUitRx6rihwMRgiIlKh4CANsnrE+LsaNsw9V63rDF2cXp/AdYZIQgyGiIjIQg57lo3ISMT1V8Vg2fpiXNGrLxI7dZRFzxUFLgZDREQEQF57lgUHaXCVTsDI3okIDQ316WOT+jCBmoiIFLNnGZEUGAwREakcV34mtWMwRESkclz5mdSOwRARkcpx5WdSOwZDREQqx5WfSe0YDBERqRxXfia1YzBERKRyct+zjEhqDIaIiEjWe5YRSY2LLhIREQD57llGJDXF9AydPXsWeXl50Ol00Ol0yMvLw7lz59o979ChQxg9ejR0Oh0iIyMxePBgVFZWSl9hIiIFMu9ZNqZvV2T1iGEgRKqgmGBowoQJKCsrQ3FxMYqLi1FWVoa8vDyn5xw7dgxDhw5Fz5498fnnn2P//v144oknEB7OGRFERETUShHDZIcOHUJxcTF27dqFQYMGAQBWr16NrKwsHD58GNdcc43d8+bOnYuRI0fiH//4h+XYFVdc4ZM6ExERkTIoIhgqKSmBTqezBEIAMHjwYOh0OuzcudNuMGQymfDRRx/hsccew/Dhw1FaWorU1FTMmTMHt956q8PHampqQlNTk+V2XV0dAMBoNMJoNIr3pNowX1eq66sB21AcbEdxsB29xzYUh5rb0Z3nrIhgqKqqCnFxcTbH4+LiUFVVZfec6upq/Prrr3j22WexcOFCPPfccyguLsbYsWPx2Wef4brrrrN73qJFi1BQUGBzfOvWrYiIiPDuibRDr9dLen01YBuKg+0oDraj99iG4lBjOzY0NLhc1q/B0Pz58+0GHm19/fXXAACNxjaJTxAEu8eB1p4hABgzZgwefvhhAEDfvn2xc+dOrFy50mEwNGfOHOTn51tu19XVISkpCTk5OYiKimr/SXnAaDRCr9cjOzsboaGhkjxGoGMbioPtKA62o/fYhuJQczuaR3Zc4ddgaMaMGRg/frzTMikpKfjmm29w6tQpm/t++eUXxMfH2z0vNjYWISEhSE9PtzqelpaGL7/80uHjabVaaLVam+OhoaGSv5B88RiBjm0oDrajONiO3mMbikON7ejO8/VrMBQbG4vY2Nh2y2VlZaG2tha7d+/GwIEDAQBfffUVamtrce2119o9JywsDL///e9x+PBhq+NHjhxBcnKy95UnIiKigKCIqfVpaWkYMWIEpkyZgl27dmHXrl2YMmUKbrnlFqvk6Z49e+L999+33H700Uexfv16rF69Gt9//z2WLVuGzZs3Y9q0af54GkRERCRDigiGAOCNN97A7373O+Tk5CAnJwe9e/fG66+/blXm8OHDqK2ttdy+7bbbsHLlSvzjH//A7373O7zyyit47733MHToUF9Xn4iIiGRKEbPJACA6OhqFhYVOywiCYHPsvvvuw3333efx45qv6U4ilruMRiMaGhpQV1enujFdsbANxcF2FAfb0XtsQ3GouR3N39v2YoNLKSYY8pf6+noAQFJSkp9rQkRERO6qr6+HTqdzWkYjuBIyqZjJZMLJkycRGRnpcBq/t8zT93/66SfJpu8HOrahONiO4mA7eo9tKA41t6MgCKivr0eXLl0QFOQ8K4g9Q+0ICgpCt27dfPJYUVFRqnuxio1tKA62ozjYjt5jG4pDre3YXo+QmWISqImIiIikwGCIiIiIVI3BkAxotVrMmzfP7srX5Bq2oTjYjuJgO3qPbSgOtqNrmEBNREREqsaeISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRiSmSNHjmDMmDGIjY1FVFQUhgwZgs8++8zf1VKkjz76CIMGDUKHDh0QGxuLsWPH+rtKitTU1IS+fftCo9GgrKzM39VRlOPHj2PSpElITU1Fhw4d0KNHD8ybNw/Nzc3+rprsLV++HKmpqQgPD8eAAQOwfft2f1dJURYtWoTf//73iIyMRFxcHG699VYcPnzY39WSLQZDMnPzzTfjwoUL+PTTT7F371707dsXt9xyC6qqqvxdNUV57733kJeXh3vvvRf79+/Hjh07MGHCBH9XS5Eee+wxdOnSxd/VUKTvvvsOJpMJq1atwsGDB7F48WKsXLkSjz/+uL+rJmvr16/HrFmzMHfuXJSWlmLYsGHIzc1FZWWlv6umGF988QWmT5+OXbt2Qa/X48KFC8jJycH58+f9XTV5Ekg2fvnlFwGAsG3bNsuxuro6AYDwySef+LFmymI0GoWuXbsKr7zyir+ronhFRUVCz549hYMHDwoAhNLSUn9XSfH+8Y9/CKmpqf6uhqwNHDhQmDp1qtWxnj17CrNnz/ZTjZSvurpaACB88cUX/q6KLLFnSEZiYmKQlpaG1157DefPn8eFCxewatUqxMfHY8CAAf6unmLs27cPJ06cQFBQEPr164fExETk5ubi4MGD/q6aopw6dQpTpkzB66+/joiICH9XJ2DU1tYiOjra39WQrebmZuzduxc5OTlWx3NycrBz504/1Ur5amtrAYCvPQcYDMmIRqOBXq9HaWkpIiMjER4ejsWLF6O4uBidOnXyd/UU44cffgAAzJ8/H3//+9/x4YcfonPnzrjuuutw5swZP9dOGQRBwD333IOpU6ciMzPT39UJGMeOHcPSpUsxdepUf1dFtk6fPo2WlhbEx8dbHY+Pj2e6gIcEQUB+fj6GDh2KjIwMf1dHlhgM+cD8+fOh0Wic/rdnzx4IgoBp06YhLi4O27dvx+7duzFmzBjccsstMBgM/n4afudqO5pMJgDA3Llz8ac//QkDBgzA2rVrodFo8M477/j5WfiXq224dOlS1NXVYc6cOf6usiy52o5tnTx5EiNGjMAdd9yByZMn+6nmyqHRaKxuC4Jgc4xcM2PGDHzzzTd46623/F0V2eJ2HD5w+vRpnD592mmZlJQU7NixAzk5OTh79iyioqIs91111VWYNGkSZs+eLXVVZc3VdiwpKcGNN96I7du3Y+jQoZb7Bg0ahJtuuglPP/201FWVLVfbcPz48di8ebPVl09LSwuCg4Pxl7/8Bf/5z3+krqqsudqO4eHhAFoDoRtuuAGDBg3CunXrEBTE36GONDc3IyIiAu+88w5uu+02y/G//vWvKCsrwxdffOHH2inPQw89hA8++ADbtm1Damqqv6sjWyH+roAaxMbGIjY2tt1yDQ0NAGDzQRkUFGTp7VAzV9txwIAB0Gq1OHz4sCUYMhqNOH78OJKTk6Wupqy52oYvv/wyFi5caLl98uRJDB8+HOvXr8egQYOkrKIiuNqOAHDixAnccMMNlh5KBkLOhYWFYcCAAdDr9VbBkF6vx5gxY/xYM2URBAEPPfQQ3n//fXz++ecMhNrBYEhGsrKy0LlzZ0ycOBFPPvkkOnTogNWrV6OiogI333yzv6unGFFRUZg6dSrmzZuHpKQkJCcn4/nnnwcA3HHHHX6unTJ0797d6vZll10GAOjRowe6devmjyop0smTJ3H99deje/fueOGFF/DLL79Y7ktISPBjzeQtPz8feXl5yMzMRFZWFv7973+jsrKSuVZumD59Ot58801s3LgRkZGRlnwrnU6HDh06+Ll28sNgSEZiY2NRXFyMuXPn4sYbb4TRaESvXr2wceNG9OnTx9/VU5Tnn38eISEhyMvLw2+//YZBgwbh008/RefOnf1dNVKRrVu34vvvv8f3339vE0QyQ8GxcePGoaamBgsWLIDBYEBGRgaKiopU37PrjhUrVgAArr/+eqvja9euxT333OP7Cskcc4aIiIhI1Th4TURERKrGYIiIiIhUjcEQERERqRqDISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVWMwRESi0Gg0+OCDD/xdDVn4/PPPodFocO7cOQDAunXr0KlTJ7/WiYgcYzBERC655557cOuttzq832AwIDc313cVUpBx48bhyJEjLpdPSUnBkiVLpKsQEVnh3mREJAo5bDwqCAJaWloQEuL9R5uY1+rQoYNfNsdsbm5GWFiYzx+XSGnYM0REomg7THb8+HFoNBps2LABN9xwAyIiItCnTx+UlJRYnbNz50784Q9/QIcOHZCUlISZM2fi/PnzlvsLCwuRmZmJyMhIJCQkYMKECaiurrbcbx6O2rJlCzIzM6HVarF9+3abupnr8/bbb+Paa69FeHg4evXqhc8//7zdawmCgH/84x+44oor0KFDB/Tp0wfvvvuu1fWLiopw9dVXo0OHDrjhhhtw/Phxq/vtDZNt2rQJmZmZCA8PR2xsLMaOHQugdWPNH3/8EQ8//DA0Gg00Go3lnPfeew+9evWCVqtFSkoKXnzxRatrpqSkYOHChbjnnnug0+kwZcoU+38sIrImEBG5YOLEicKYMWMc3g9AeP/99wVBEISKigoBgNCzZ0/hww8/FA4fPizcfvvtQnJysmA0GgVBEIRvvvlGuOyyy4TFixcLR44cEXbs2CH069dPuOeeeyzXXLNmjVBUVCQcO3ZMKCkpEQYPHizk5uZa7v/ss88EAELv3r2FrVu3Ct9//71w+vRpm7qZ69OtWzfh3XffFcrLy4XJkycLkZGRlvKOrvX4448LPXv2FIqLi4Vjx44Ja9euFbRarfD5558LgiAIlZWVglarFf76178K3333nVBYWCjEx8cLAISzZ88KgiAIa9euFXQ6naU+H374oRAcHCw8+eSTQnl5uVBWViY8/fTTgiAIQk1NjdCtWzdhwYIFgsFgEAwGgyAIgrBnzx4hKChIWLBggXD48GFh7dq1QocOHYS1a9darpucnCxERUUJzz//vHD06FHh6NGjrv1xiVSOwRARucSTYOiVV16x3H/w4EEBgHDo0CFBEAQhLy9PuP/++62usX37diEoKEj47bff7D7G7t27BQBCfX29IAgXA5gPPvjAad3N9Xn22Wctx4xGo9CtWzfhueeec3itX3/9VQgPDxd27txpdb1JkyYJf/7znwVBEIQ5c+YIaWlpgslkstz/t7/9zWkwlJWVJfzlL39xWN/k5GRh8eLFVscmTJggZGdnWx179NFHhfT0dKvzbr31VictQUT2cJiMiCTTu3dvy78TExMBwDLMtXfvXqxbtw6XXXaZ5b/hw4fDZDKhoqICAFBaWooxY8YgOTkZkZGRuP766wEAlZWVVo+TmZnpUn2ysrIs/w4JCUFmZiYOHTrk8Frl5eVobGxEdna2VT1fe+01HDt2DABw6NAhDB482Go4q+3j2FNWVoY//vGPLtXZ7NChQxgyZIjVsSFDhuDo0aNoaWmxW38icg0TqIlIMqGhoZZ/m4MFk8lk+f8HHngAM2fOtDmve/fuOH/+PHJycpCTk4PCwkJcfvnlqKysxPDhw9Hc3GxVvmPHjh7XsW0Qc+m1zHX96KOP0LVrV6tyWq0WQGuitbs8SaYWBMGmrvYe25u2IFIrBkNE5Bf9+/fHwYMHceWVV9q9/9tvv8Xp06fx7LPPIikpCQCwZ88erx5z165d+MMf/gAAuHDhAvbu3YsZM2Y4LJ+eng6tVovKykpcd911Dstcur7Srl27nNajd+/e+O9//4t7773X7v1hYWFWvT3mx/nyyy+tju3cuRNXX301goODnT4eETnHYIiIXFZbW4uysjKrY9HR0ejevbvb1/rb3/6GwYMHY/r06ZgyZQo6duyIQ4cOQa/XY+nSpejevTvCwsKwdOlSTJ06FQcOHMBTTz3lVf3/9a9/4aqrrkJaWhoWL16Ms2fP4r777nNYPjIyEv/3f/+Hhx9+GCaTCUOHDkVdXR127tyJyy67DBMnTsTUqVPx4osvIj8/Hw888IBl+M+ZefPm4Y9//CN69OiB8ePH48KFC/j444/x2GOPAWidFbZt2zaMHz8eWq0WsbGxeOSRR/D73/8eTz31FMaNG4eSkhIsW7YMy5cv96pNiAicTUZErpk4caIAwOa/iRMnCoJgP4G6tLTUcv7Zs2cFAMJnn31mObZ7924hOztbuOyyy4SOHTsKvXv3tsyqEgRBePPNN4WUlBRBq9UKWVlZwqZNm6yua056NicqO2Kuz5tvvikMGjRICAsLE9LS0oT//ve/ljKOrmUymYR//vOfwjXXXCOEhoYKl19+uTB8+HDhiy++sJTZvHmzcOWVVwparVYYNmyY8OqrrzpNoBYEQXjvvfeEvn37CmFhYUJsbKwwduxYy30lJSVC7969Ba1WK7T9mH733XeF9PR0ITQ0VOjevbvw/PPPW13TXuI1EbVPIwgeDHgTESnI8ePHkZqaitLSUvTt29ff1SEimeFsMiIiIlI1BkNERESkahwmIyIiIlVjzxARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFTt/wFV7lOXBH9k1gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.grid(True)\n",
    "plt.plot(result1.predict(linear=True), result1.resid_pearson, \"o\")\n",
    "plt.xlabel(\"Linear predictor\")\n",
    "plt.ylabel(\"Residual\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative variance function is mu^2 * (1 - mu)^2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
 
 
 
 
    },
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "class vf(sm.families.varfuncs.VarianceFunction):\n",
    "    def __call__(self, mu):\n",
    "        return mu ** 2 * (1 - mu) ** 2\n",
    "\n",
    "    def deriv(self, mu):\n",
    "        return 2 * mu - 6 * mu ** 2 + 4 * mu ** 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the quasi-binomial regression with the alternative variance\n",
    "function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                 blotch   No. Observations:                   90\n",
      "Model:                            GLM   Df Residuals:                       72\n",
      "Model Family:                Binomial   Df Model:                           17\n",
      "Link Function:                  Logit   Scale:                         0.98855\n",
      "Method:                          IRLS   Log-Likelihood:                -21.335\n",
      "Date:                Fri, 20 Mar 2026   Deviance:                       7.2134\n",
      "Time:                        11:18:46   Pearson chi2:                     71.2\n",
      "No. Iterations:                    25   Pseudo R-squ. (CS):             0.3115\n",
      "Covariance Type:            nonrobust                                         \n",
      "==================================================================================\n",
      "                     coef    std err          z      P>|z|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------\n",
      "C(variety)[1]     -7.9224      0.445    -17.817      0.000      -8.794      -7.051\n",
      "C(variety)[2]     -8.3897      0.445    -18.868      0.000      -9.261      -7.518\n",
      "C(variety)[3]     -7.8436      0.445    -17.640      0.000      -8.715      -6.972\n",
      "C(variety)[4]     -6.9683      0.445    -15.672      0.000      -7.840      -6.097\n",
      "C(variety)[5]     -6.5697      0.445    -14.775      0.000      -7.441      -5.698\n",
      "C(variety)[6]     -6.5938      0.445    -14.829      0.000      -7.465      -5.722\n",
      "C(variety)[7]     -5.5823      0.445    -12.555      0.000      -6.454      -4.711\n",
      "C(variety)[8]     -4.6598      0.445    -10.480      0.000      -5.531      -3.788\n",
      "C(variety)[9]     -4.7869      0.445    -10.766      0.000      -5.658      -3.915\n",
      "C(variety)[10]    -4.0351      0.445     -9.075      0.000      -4.907      -3.164\n",
      "C(site)[T.2]       1.3831      0.445      3.111      0.002       0.512       2.255\n",
      "C(site)[T.3]       3.8601      0.445      8.681      0.000       2.989       4.732\n",
      "C(site)[T.4]       3.5570      0.445      8.000      0.000       2.686       4.428\n",
      "C(site)[T.5]       4.1079      0.445      9.239      0.000       3.236       4.979\n",
      "C(site)[T.6]       4.3054      0.445      9.683      0.000       3.434       5.177\n",
      "C(site)[T.7]       4.9181      0.445     11.061      0.000       4.047       5.790\n",
      "C(site)[T.8]       5.6949      0.445     12.808      0.000       4.823       6.566\n",
      "C(site)[T.9]       7.0676      0.445     15.895      0.000       6.196       7.939\n",
      "==================================================================================\n"
     ]
    }
   ],
   "source": [
    "bin = sm.families.Binomial()\n",
    "bin.variance = vf()\n",
    "model2 = sm.GLM.from_formula(\"blotch ~ 0 + C(variety) + C(site)\", family=bin, data=df)\n",
    "result2 = model2.fit(scale=\"X2\")\n",
    "print(result2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the alternative variance function, the mean/variance relationship\n",
    "seems to capture the data well, and the estimated scale parameter is\n",
    "close to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
 
 
 
 
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/genmod/generalized_linear_model.py:985: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAQHRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcrZGZzZzEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvhF0PpwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQQ9JREFUeJzt3Xt0VNX9///XJCYTAkkgREiQAClIIaSAQEHqpWjlEhShtn5V1B9YtQWBWmlrq/gpYL3WWmzpR+pHK2pZVJZXRDEab4AKpVyCchGFBrESoNwSBBOGzP79QScwJJNMJjNzLvN8rOWSzJxM3rPn5Jx39n7vvT3GGCMAAACHS7I6AAAAgGggqQEAAK5AUgMAAFyBpAYAALgCSQ0AAHAFkhoAAOAKJDUAAMAVzrA6gHjy+/3atWuXMjIy5PF4rA4HAACEwRijw4cPq1OnTkpKCt0fk1BJza5du5Sfn291GAAAIAJffPGFOnfuHPL5hEpqMjIyJJ1olMzMTIujOcnn8+nNN9/UiBEjlJKSYnU4jkLbRY62ixxtFznaLjKJ3m5VVVXKz8+vu4+HklBJTWDIKTMz03ZJTXp6ujIzMxPyZG0J2i5ytF3kaLvI0XaRod1OaKp0hEJhAADgCiQ1AADAFUhqAACAK5DUAAAAVyCpAQAArkBSAwAAXIGkBgAAuAJJDQAAcAWSGgAA4AoJtaIwAPup9RutLj+gvYer1SEjTYMLspWcxIazAJqPpAaAZUo2Vmj2ks2qqKyueywvK00zxxRqVFGehZEBcCLHDD/NmzdPffv2rdu3aejQoXr99detDgtAhEo2VmjygnVBCY0k7a6s1uQF61SyscKiyAA4lWOSms6dO+uBBx7QmjVrtGbNGl188cUaO3asNm3aZHVoAJqp1m80e8lmmQaeCzw2e8lm1fobOgIAGuaYpGbMmDEaPXq0evbsqZ49e+ree+9VmzZttGrVKqtDA9BMq8sP1OuhOZWRVFFZrdXlB+IXFADHc2RNTW1trZ577jkdOXJEQ4cODXlcTU2Nampq6r6uqqqSdGILd5/PF/M4wxWIxU4xOQVtFzkr267i0JGwj/P5MmMcTfNx3kWOtotMordbuO/bY4xxTP/uxx9/rKFDh6q6ulpt2rTRwoULNXr06JDHz5o1S7Nnz673+MKFC5Wenh7LUAE04rNKj/68ObnJ46YW1ursLMdcogDEyNGjRzV+/HhVVlYqMzP0HzqOSmqOHTumnTt36tChQ3rhhRf0xBNPaNmyZSosLGzw+IZ6avLz87Vv375GGyXefD6fSktLNXz4cKWkpFgdjqPQdpGzsu1q/UbDHl6uPVU1DdbVeCTlZnn17vQLbTm9m/MucrRdZBK93aqqqpSTk9NkUuOo4afU1FT16NFDkjRo0CD985//1B//+Ec99thjDR7v9Xrl9XrrPZ6SkmLLk8KucTkBbRc5K9ouRdKsy/to8oJ18khBiU0ghZk5po/SvKlxjau5OO8iR9tFJlHbLdz37JhC4YYYY4J6YgA4x6iiPM27boBys9KCHs/NStO86wawTg2AZnNMT82dd96p4uJi5efn6/Dhw3r22Wf13nvvqaSkxOrQAERoVFGehhfmsqIwgKhwTFKzZ88eXX/99aqoqFBWVpb69u2rkpISDR8+3OrQALRAcpJHQ7u3tzoMAC7gmKTmr3/9q9UhAAAAG3N0TQ0AAEAASQ0AAHAFkhoAAOAKJDUAAMAVSGoAAIArkNQAAABXIKkBAACuQFIDAABcgaQGAAC4AkkNAABwBZIaAADgCiQ1AADAFUhqAACAK5DUAAAAVyCpAQAArkBSAwAAXIGkBgAAuAJJDQAAcAWSGgAA4AokNQAAwBVIagAAgCuQ1AAAAFcgqQEAAK5AUgMAAFzhDKsDAAC3q/UbrS4/oL2Hq9UhI02DC7KVnOSxOizAdUhqACCGSjZWaPaSzaqorK57LC8rTTPHFGpUUZ6FkQHuw/ATAMRIycYKTV6wLiihkaTdldWavGCdSjZWWBQZ4E4kNQAQA7V+o9lLNss08FzgsdlLNqvW39ARACJBUgMAMbC6/EC9HppTGUkVldVaXX4gfkEBLkdSAwAxsPdw6IQmkuMANI2kBgBioENGWlSPA9A0khoAiIHBBdnKy0pTqInbHp2YBTW4IDueYQGuRlIDADGQnOTRzDGFklQvsQl8PXNMIevVAFFEUgMAMTKqKE/zrhug3KzgIabcrDTNu24A69QAUcbiewAQQ6OK8jS8MJcVhYE4IKkBgBhLTvJoaPf2VocBuB7DTwAAwBVIagAAgCuQ1AAAAFcgqQEAAK5AUgMAAFyBpAYAALgCSQ0AAHAFkhoAAOAKLL4HAC5V6zesZIyEQlIDAC5UsrFCs5dsVkVldd1jeVlpmjmmkD2n4FoMPwGAy5RsrNDkBeuCEhpJ2l1ZrckL1qlkY4VFkQGxRVIDAC5S6zeavWSzTAPPBR6bvWSzav0NHQE4G0kNALjI6vID9XpoTmUkVVRWa3X5gfgFBcQJSQ0AuMjew6ETmkiOA5yEpAYAXKRDRlpUjwOchKQGAFxkcEG28rLSFGritkcnZkENLsiOZ1hAXJDUAICLJCd5NHNMoSTVS2wCX88cU8h6NXAlkhoAcJlRRXmad90A5WYFDzHlZqVp3nUDWKcGrsXiewDgQqOK8jS8MJcVhZFQSGoAwKWSkzwa2r291WEAccPwEwAAcAWSGgAA4AoMPwEJoKHdmgHAbUhqAJcLtVvzjOJvWhgVAEQfw0+AizW2W/O0Zzdow35mwgBwD5IawKXC2a35xR1J7NYMwDVIagCXCme35kPHPFrz+cH4BQUAMURSA7hU+Ls118Q4EgCID8ckNffff7++/e1vKyMjQx06dNC4ceO0detWq8MCbCv83Zq9MY4EAOLDMUnNsmXLNGXKFK1atUqlpaU6fvy4RowYoSNHjlgdGmBL4ezW3DbVaFDXdvEMCwBixjFTuktKSoK+nj9/vjp06KC1a9fqwgsvbPB7ampqVFNzsmu9qqpKkuTz+eTz+WIXbDMFYrFTTE5B2zVuRvE3Ne3ZDfJIQQXDgUTnim5++WuPy+djFlRzcN5FjraLTKK3W7jv22OMceTUh23btunss8/Wxx9/rKKiogaPmTVrlmbPnl3v8YULFyo9PT3WIQK2sGG/Ry/uSNKhYycTl7apRld086tfe0f++gNIMEePHtX48eNVWVmpzMzMkMc5Mqkxxmjs2LE6ePCgVqxYEfK4hnpq8vPztW/fvkYbJd58Pp9KS0s1fPhwpaSkWB2Oo9B24an1G635/KD2Hq5RhwyvBnVtJ3/tcdouQpx3kaPtIpPo7VZVVaWcnJwmkxrHDD+daurUqfroo4/0/vvvN3qc1+uV11u/CDIlJcWWJ4Vd43IC2q5xKZLO79kx6LHAkBNtFznaLnK0XWQStd3Cfc+OS2qmTZumV155RcuXL1fnzp2tDgcAANiEY5IaY4ymTZuml156Se+9954KCgqsDgkAANiIY5KaKVOmaOHChVq8eLEyMjK0e/duSVJWVpZatWplcXQAAMBqjlmnZt68eaqsrNSwYcOUl5dX99+iRYusDg0AANiAY3pqHDhJCwAAxJFjemoAAAAaQ1IDAABcgaQGAAC4AkkNAABwBZIaAADgCiQ1AADAFUhqAACAK5DUAAAAVyCpAQAArkBSAwAAXIGkBgAAuAJJDQAAcAWSGgAA4AokNQAAwBVIagAAgCuQ1AAAAFc4w+oAANhXrd9odfkB7T1crQ4ZaRpckK3kJI/VYQFAg0hqADSoZGOFZi/ZrIrK6rrH8rLSNHNMoUYV5VkYGQJIOoFgJDUAgtT6jf78zjbNeevTes/trqzW5AXrNO+6ASQ2FiPpBOqjpgZAnZKNFTrvgbcbTGgkyfz3/7OXbFat3zR4DGKvZGOFJi9YF5TQSCeTzpKNFRZFBliLpAaApJM3yt1VNY0eZyRVVFZrdfmB+ASGILV+o9lLNquhlJKkE4mOpAZAozfKUPYerm76IETd6vID9XpoTkXSiURGUgOgyRtlQzpkpMUoGjQm3GSSpBOJiEJhAM26AXok5WadmGmD+As3mSTpRCKipwZAs2+AM8cUMnXYIoMLspWXlaZQre/RiVlQJJ1IRCQ1AJq8UQbkZaUxndtiyUkezRxTKEn1Pq/A1ySdSFQkNQAavVEG3HbJ2Xr/VxeT0NjAqKI8zbtugHKzgnvYckk6keCoqQEg6eSNkgXdnGFUUZ6GF+ayojBwCpIaAHW4UTpLcpJHQ7u3tzoMwDZIagAE4UYJwKmoqQEAAK5ATw0AwFHYnRyhkNQAAByD3cnRGIafAACOwO7kaApJDQDA9tidHOEgqQEA2B67kyMcJDUAANtjd3KEg6QGAGB77E6OcJDUAABsj93JEQ6SGgCA7bE7OcJBUgMAYar1G63cvl9LPqrQZ5UeZtrEGbuToyksvgcAYai/6Fuynn94uWZd3oebaRyx6SoaQ1IDAE0ILPp2er/MnqoaTV6wjl6COGPTVYTC8BMANIJF3wDnIKkBgEaw6BvgHCQ1ANAIFn0DnIOkBgAawaJvgHOQ1ABAI1j0DXAOkhoAaASLvgHOQVIDAE0Iveibl+ncgI2wTg0AhOHURd8qDh3RvzaVaepVFyrNm2p1aAD+i6QGAMIUWPTN58vU0n+vZ8gJsBmGnwAAgCuQ1AAAAFdg+AlAzNX6DRsQAog5khoAMVV/d+sT67rMHFPIrCEAUcXwE4CYCexuffreSbsrqzV5wTqVbKywKDIAbkRSAyAm2N0aQLyR1ACICXa3BhBv1NQALmSHwlx2twYQb2EnNX/605/CftGf/vSnEQUDoOXsUpjL7tYA4i3spGbOnDlhHefxeEhqgDhoqDemdPNuTV6wrl4dS6AwN577FAV2t95dWd1gXY1HUi67WwOIorCTmvLy8ljGAaAZGuqNyc1MU/Xx2pCFuR6dKMwdXpgbl6GowO7Wkxesk0cKiisWu1vbYcgNgLUcVSi8fPlyjRkzRp06dZLH49HLL79sdUhA3IWcJl1VrUNHfSG/z4rC3NC7W6dFtdeoZGOFzn/wHV3z+Crd+myZrnl8lc5/8B2mjAMJJuJC4X//+9965ZVXtHPnTh07dizouT/84Q8tDqwhR44cUb9+/XTDDTfoBz/4QUx+BmBnjU2TDle8C3NP3d06Fr0ogSTPDkNuAKwVUVLz9ttv6/LLL1dBQYG2bt2qoqIi7dixQ8YYDRgwINox1ikuLlZxcXHMXh+wu6amSYfDisLcwO7W0dbUWjjxHnIDYK2Ikpo77rhDP//5z3X33XcrIyNDL7zwgjp06KBrr71Wo0aNinaMEaupqVFNTU3d11VVVZIkn88nny90N328BWKxU0xOkWhtV3HoSMTfe6Iw16tzOmcE/Q44ue3+EeZaOCu37dWQKBYku6HtrELbRSbR2y3c9+0xxjS7JzsjI0NlZWXq3r272rVrp/fff199+vTRhg0bNHbsWO3YsaO5L9lsHo9HL730ksaNGxfymFmzZmn27Nn1Hl+4cKHS09NjGB0QG59VevTnzclhHBnopzj1a+lHPf3q1949K/iu3efRM5813R7/39m1GpjjnvcNJJqjR49q/PjxqqysVGZmZsjjIuqpad26dV0PSKdOnbR9+3b16dNHkrRv375IXjIm7rjjDk2fPr3u66qqKuXn52vEiBGNNkq8+Xw+lZaWavjw4UpJSbE6HEdJtLar9Rs9//By7amqCTlNOis9Rd5kj/YcPlnrlpeVphnFvTSyT8e6x9zQdu3LD+iZz9Y0edyIC4ZEvafG6W1nFdouMoneboGRlqZElNSce+65+uCDD1RYWKhLL71UP//5z/Xxxx/rxRdf1LnnnhvJS8aE1+uV1+ut93hKSootTwq7xuUEidJ2KZJmXd6n0WnSD1zxrWYV5jq57Yb26BDWWjhDe3SISU2Nk9vOarRdZBK13cJ9zxElNX/4wx/01VdfSToxxPPVV19p0aJF6tGjR9iL9AGITGCadL11ak5bNTgWhbl2E++1cADYW0RJzTe+8Y26f6enp+vRRx+NWkCN+eqrr7Rt27a6r8vLy1VWVqbs7Gx16dIlLjEAdhDradJOEm6SB8D9HLWh5Zo1a3TRRRfVfR2ol5kwYYKeeuopi6ICrBGradJORJIHQIowqUlKSpLHE/piUVtbG3FAjRk2bJgimKwFIAGQ5AGIKKl56aWXgr72+Xxav369nn766QanUAMAAMRaREnN2LFj6z32wx/+UH369NGiRYt04403tjgwwO3YgBEAoiuqNTVDhgzRzTffHM2XBFypoV228yhsBYAWidou3V9//bXmzp2rzp07R+slAVcKucv2fzdgdNLO0rV+o5Xb92tx2ZdauX2/av3UvAGwTkQ9Ne3atQsqFDbG6PDhw0pPT9eCBQuiFhycjyGWYG7agJHeJgB2E1FSM2fOnKCkJikpSWeeeaaGDBmidu3aRS04OBs3vfqa2mU7sAHj6vIDtp7JE+htOj05C/Q2zbtuQMJ+xgCsE1FSM3HixCiHAbfhptewvYdDJzSRHGcFN/U2AXCXsJOajz76KOwX7du3b0TBwB246YXWISMtqsdZwS29TbAvhq0RqbCTmv79+8vj8dQtfmfF4ntwBm56oQ0uyA5rA8bBUdxROtrc0NsE+2LYGi0R9uyn8vJy/etf/1J5eblefPFFFRQU6NFHH9X69eu1fv16Pfroo+revbteeOGFWMYLB+CmF1pgA0bp5IaLAU7ZgNENvU04yU4z2Nw0MxDWCLunpmvXrnX/vvLKK/WnP/1Jo0ePrnusb9++ys/P1//8z/9o3LhxUQ0SzsJNr3FO34DRDb1NOMFOvSIMWyMaIioU/vjjj1VQUFDv8YKCAm3evLnFQcHZuOk1zckbMAZ6myYvWCePFPQZO6W3CfYr5mfYGtEQ0eJ7vXv31j333KPq6pMnYE1Nje655x717t07asHBmdwwxBIPgQ0Yx/Y/S0O7t3dUewR6m3KzgnvbcrPSEnZmm5M01SsinegViedQFMPWiIaIemr+8pe/aMyYMcrPz1e/fv0kSRs2bJDH49Grr74a1QDhTE4fYkHTnNzblOjs2CvCsDWiIaKkZvDgwSovL9eCBQv0ySefyBijq666SuPHj1fr1q2jHSMcipue+wV6m+AsduwVYdga0RDxhpbp6en68Y9/HM1Y4ELc9AD7sWOvCLVaiIawk5pXXnlFxcXFSklJ0SuvvNLosZdffnmLAwOAcLBQW/PZtVeEYWu0VNhJzbhx47R792516NCh0SnbHo+HxfcAxIWdpiQ7iZ17RRi2RkuEPfvJ7/erQ4cOdf8O9R8JDYB4YKG2lrHzDDYnzwyEtSKuqTndoUOH1LZt22i9HACExEJt0UGvCNwmonVqHnzwQS1atKju6yuvvFLZ2dk666yztGHDhqgFBwANac6UZDSOXhG4SURJzWOPPab8/HxJUmlpqd566y2VlJSouLhYv/zlL6MaIACczo5TkgFYL6Lhp4qKirqk5tVXX9X/+3//TyNGjFC3bt00ZMiQqAYIAKez45RkANaLqKemXbt2+uKLLyRJJSUluuSSSyRJxhgKhQHEXGBKcqiBEo9OzIJioTYgsUSU1FxxxRUaP368hg8frv3796u4uFiSVFZWph49ekQ1QAA4HfuLRa7Wb7Ry+34tLvtSK7fvj+v+TkCsRTT8NGfOHHXr1k1ffPGFfve736lNmzaSTgxL3XLLLVENEAAawkJtzRfJuj4sbug8ifyZRZTUpKSk6Be/+EW9x3/2s5+1NB4ACBtTksMXWNfn9H6ZwLo+Da1Nw+KGzpPon1lEw0+S9Le//U3nn3++OnXqpM8//1yS9Mgjj2jx4sVRCw4AmsKU5KY1ta6PdGJdn1OHoljc0Hn4zCJMaubNm6fp06eruLhYhw4dqisObtu2rR555JFoxgcAaKHmrusTSRIEa/GZnRBRUjN37lw9/vjjmjFjhpKTk+seHzRokD7++OOoBQcAiMypBcEfbNsX1vcE1vVhcUPn4TM7IaKamvLycp1zzjn1Hvd6vTpy5EiLgwIARK6huopwBNb1YXFD5wn3s3j9v0NQbq09i6inpqCgQGVlZfUef/3119W7d++WxgQAiFCouorGnL6uD4sbOk+4n8UzKz/XNY+v0vkPvuPKGpuIemp++ctfasqUKaqurpYxRqtXr9bf//533XffffrrX/8a7RgBAGForK4ilIbW9Qksbri7srrB1/LoxNR5Fje0j6Y+s9M1NuvNySJKam644QYdP35ct99+u44eParx48frrLPO0ty5c3XBBRdEO0YgYSTy+hJouabqKhrS0Lo+gcUNJy9YJ48UdJNkcUN7auwza4hbd7OPKKmRpJtvvlk333yz9u3bJ7/fr9raWt13332aMmWKvv7662jGmDBq/UZrtu/nhpag4rm+RK3f6B/lB7R2n0ftyw9oaI8OnGsuEG5dxdSLuuvsjhmNXmdY3NB5Qn1moZxaPDy0e/vYBxgHzUpqDh06pClTpujNN99USkqKfv3rX2vq1KmaPXu2fv/736uwsFBPPvlkrGJ1tQ37Pbr/4eXaXVVT91giLZiU6CJZGK0lP+vkRS9Zz3y2hnPNJcKtqzivx5lh3cRY3NB5Tv3MXt9YoWdWft7k97ip4LtZhcJ33nmnli9frgkTJig7O1u33XabLrvsMq1YsUJLly7VP//5T11zzTWxitW13ti0R09+mhSU0EiJtWCSE0VrD514ri/B4lzuFouNPlnc0HkCn1lxmH+kuKngu1k9Na+99prmz5+vSy65RLfccot69Oihnj17suBeC9T6je5Z+kmDz7l1zNMuWlK/Es2houasL9GSLuKmkifONeejFganSsSC72b11OzatUuFhSd2xv3GN76htLQ03XTTTTEJLFGsLj/w3x6ahi8yibJgUryVbKzQ+Q++o2seX6Vbny1r1hTHaPd2xGtNEBbnSgyBuorcrOC/vnOz0lw30wWNS8Td7JvVU+P3+5WSklL3dXJyslq3bh31oBIJi1zFX0vqV2LR2xGvNUE41xIHtTAISLSC72YlNcYYTZw4UV6vV5JUXV2tSZMm1UtsXnzxxehF6HIschVfLU1KYjFUFK8uYs61xBKoqwASKcltVlIzYcKEoK+vu+66qAaTiAYXZCs306vdVdVqaAjKjWOeVmppUhKL3o541UEk4vg6gBMSJcltVlIzf/78WMWRsJKTPLprdC9NfbaMwr44aGlSEqvejnh0EVNECsDtIl58D9Ezsk9H/ainX0t3pwdN63brmKeVWpqUxLK3Ix5dxIk2vg4gsZDU2ES/9ka3X3uh1v/7sOvHPK3U0qQk1r0d8egiDiRPK7ft1Zsr/qERFwxhRWEArhDRLt2IDRa5ir1oTHF0w5TZ5CSPhhRka2CO0RCSZ1uK1uKOQCKhpwYJJxpDMIk0mwDxF899wAA3IalBQopGUpIoswkQX/HcByxRtWQ1cdgbSQ0SFkkJ7IatLGKPXjB3o6YGsDlqKxIHW1nEVqgtTioqqzWJDV1dgZ4awMb4qzKxsJVF7DTWCxbw6xc/phfM4eipAWwq2htnwv7YyiJ2muoFk6RDR3368zvb4hQRYoGkBrChpmorpBO1FQxFuUtgHaVQ/QQeneipYyuL5gu3d2v+h+X8XjkYSQ1gQ9RWJKZorKOEhoXbu3XoqI/fKwcjqQFsiNqKxOWGxR3taHBBttq2SgnrWH6vnItCYcCGqK1IbCzuGH3JSR7dcF43zXnrsyaP5ffKuUhqABuK5caZcAbWUYq+qRefrfkf7tCho74Gn+f3yvkYfgJsiNoK2IWb1klKTvLogSu+1WAhNr9X7kBPDRrFcuLWicYeVU7EOWcfblwnqbHfq/+5tFBZrVK1uOxLzj2HIqlBSG68oDlNotVWxOKcI0mKjJv3oGro9+rgkRr99jWud05HUoMGufmC5jSJUlsRi3OOxDwyibAH1am/VyUbKzRl4Xqudy5ATQ3qYeE3xFsszjlWZI5cIq2TxPXOXUhqUE8iXdBgD9E+57hRtUwirZPE9c5dSGpQTyJd0GAP0T7nuFG1TCKtk8T1zl0cl9Q8+uijKigoUFpamgYOHKgVK1ZYHZLrJNIFDfYQ7rmU09ob1nHcqFomkfag4nrnLo5KahYtWqSf/exnmjFjhtavX68LLrhAxcXF2rlzp9WhuUoiXdBgD02dcwE/f25DWLUw3KhaJpHWSeJ65y6Omv30hz/8QTfeeKNuuukmSdIjjzyiN954Q/PmzdP9999f7/iamhrV1NTUfV1VVSVJ8vl88vkaXlHSCoFY7BTTjOJvatqzG+SRguoSPKc87689Ln+tBcGdwo5t5xR2a7tQ59yp9lSdKPKde3U/jezTMeRrndM5Q7mZXu2pqmlkRWavzumcEdH7t1vbxcL3vpmjuVf30z1LP9HuqpPX0dwsr2YU99L3vpnjmrZzwvXOju0WT+G+b48xxhGVcseOHVN6erqee+45ff/73697/NZbb1VZWZmWLVtW73tmzZql2bNn13t84cKFSk9Pj2m8brBhv0cv7kjSoWMn/4Zpm2p0RTe/+rV3xGkDC/iNtL3KoyqflJkidc80CvcP+g37PXqhPEmVvsa+wahtqjRzQG2jr7thv0dPfhrojD71wBPn7o96ch6HoyWfp5NwvbO3o0ePavz48aqsrFRmZmbI4xyT1OzatUtnnXWWPvjgA33nO9+pe/y+++7T008/ra1bt9b7noZ6avLz87Vv375GGyXefD6fSktLNXz4cKWkhLeLbLzU+o3WfH5Qew/XqEOGV4O6trNVl7Od287uYtF2b2zaU/8v+0yv7hrdq9GelVN9uH2/Jjy1tsnjFvxokIY0MSTQUDx5/+1pCDeehnDeRc7ObWfn652d2y0eqqqqlJOT02RS46jhJ0nyeIJPMGNMvccCvF6vvN76hYUpKSm2PCnsGFeKpPN7Rn7xjxc7tp1TRKvtSjZWaNqzG+oN9+ypqtG0ZzeEvYDZoerw+vj3Hz3eZNyX9e+s4r5nxWxFYc67yNmx7ZxwvbNju0mxX7k73PfsmKQmJydHycnJ2r17d9Dje/fuVceO9j4JAbeL5gq00S7yTZQVmZ2m1m+0Zvt+tq9wATut3O2Y2U+pqakaOHCgSktLgx4vLS0NGo4CEH/RXBeG2Sjut2G/R8MeXq5rHl+lW58t0zWPr9K3731LSz9ilWensdvK3Y5JaiRp+vTpeuKJJ/Tkk09qy5Ytuu2227Rz505NmjTJ6tCAhBbNdWESaTpxInpj0x49+WlSUJ2TJB04cky3LFyne1/bbFFkaC47rtztmOEnSbrqqqu0f/9+3X333aqoqFBRUZGWLl2qrl27Wh0a4Ainjnu3Tz9D0brWRHvIaFRRnuZdN6Bel3auCzejTKRdxGv9Rvcs/aTRYx5fUS7JaMalfeITFCLWnB7aeA0BOyqpkaRbbrlFt9xyi9VhAI7T0Lh329RkpXTbo8v6d27RaweGjHZXVjeyLkzzhoxGFeVpeGGuq2/4dqpFiIfV5Qf+20PT+Gf4+IodOie/nUb37RSfwBARO67c7ajhJwCRCTXufeiYNO3Z8FbpbUyshowCRb5j+5+lod3buy6hsVMtQjw05+Z21+KNbDhqc3ZcuZukBnC5xsa9AylHNMa9A0NGuVnBF7DcrLSwp3MnCjvWIsRDc25uB4742HDU5uxY1O+44ScgWhKlliGe496JMGQUDXasRYiHwQXZyk5P0YGj4S15z4aj9hbooZ28YF3ILSbiXdRPUoOElEi1DPEe92ZdmKbF8zOxU/KenOTRrDG99dNFG9RUXY3EhqNOYLeifpIaJJxALcPpHfuBWga3DZXYcdw70cXrM7Fj8l5clKuLlq/XuxXJjR7HWkTOYaceWmpqkFASsZbBjuPeiS4en4mdC5HHdTP60Xe6hHzeI9Yichq7FPWT1CChRHPlW6dobGZSIJVz+g2k1m+0cvt+LS77Uiu377d9UhrrBQadkLzfUdxLj44/R9mtg/f0yaOwHC3A8BMSih3XVYiHUUV5+t/x5+iuxRt14MjJIs22qdI9V/Rz9A3EjkMs4YhlLYJTCpFH9+2kkUV5thi2gDuQ1CChuLW+pKli0JKNFfrta1uCEprs9BSN7VytkX2cuyGs0+ujYlWL4KTkncJyRBNJDRJKLFa+tVpTPRWhbvwHj/o0/9MkDdzU8hWFrRDNncGtFIubuluTd6Ap1NQgobhts8SmikGXflTRZG3Fva9/YvsalIYkYn1UuFpaiOy0GiUggJ4aJBy7rasQqXB6Kv5n8UbtP3KskVfxqKKyxvLaikg4aYgl3lqyKJpTa5QAiaQGCcpO6ypEKpyeisYTmpOceONniKVxkSTv8apR8hvpH+UHtP/ocUf+7sG+SGqQsJxeoBjNRMSJN3431kdFW3OS93jVKL2xaY9mr0vWoVVr6h6jJwjRQk0N4FDhJiLZrVMaWZDeKC/L68gbv9vqo2Il3EXR4lGjVLKxQtOe3aBDp3Ug2mFBQLgDSQ0QZ9Eqwgy3GPSesUV1X5/+vCTNKO7l2Bs/O4NHT6xrlIJ7goLPN7ssCAjnY/gJiKNoFmGGWww6qihP85I8DdRWeFXc8aij16mR3FEfZQexrlFyyoKAcDaSGiBOYlGEGW4xaEM3/nM6Z+iNktdb+rZswen1UXYQ6xolZqshHkhqgDiIZRFmuD0Vp9/4fT7f6S+FBNaSaeDhYLYa4oGaGiAOYl2EaZcdcuFssaxRYrd4xAM9NUAc0PWeuJral8tuYlWjdGpP0Mn+yROYrYZoIalBSE67GMdTc9uGrvfE5NTVeWNVozSqKE9zr+6nu14sC5rW7bTVvGFfJDVokFMvxvEQSduwUFzicfoO4rEysk9H+XbU6szCc1lRGFFHTQ3qaWqTxEReICvStmGhuMTSVGG4lNhrsiR5pCEF2dSAIepIahCEi3FoLW0bFooLjxt2iGYHccAaDD8hCAtkhRaNtmGhuMZZMewZi9oxCsMBa5DUIAgX49Ci1TYsFNcwK2pQGkqislun6p6xRRrdN/KfZdfCcIr/4XYkNQhi14uxHdA2sROvHaJPFSqJOnDkmG5ZuE4/+XeB7hhdGNFr27EwnOJ/JAJqahCEBbJCo21iJ941KI0lUQGPLS/Xq2VfRvT6disMp/gfiYKkBkHsdjG2E9omduI97NlUEhUwbVGZln4U2Q3fLoXhFP8jkTD8hHrC3SQxEdE2sRHvob1wkyNjpFsWrtNfkiJLQuxQGE7xPxIJSQ0aZIeLsV3RNtEX7xqU5iZHd770sS7u1VGpZzS/c9vqwnCK/5FISGoQktUXYzujbaIr1jtEn25wQbayW6fqwJFjTR8s6cARn869/23d9/0ix/XG7dh3NKzjKHCHG1BTA8AW4lmDkpzk0T1ji5r1PQeOHHNcUW3Jxgo98tanTR6X3TpFu6uqHbvYIRBATw0A22ju0F5L1l0Z3TdPP/l3gR5bXt6sGGcv2axhZ1/QrO+xQjgzvAIOHPHptkVlkpjmDWcjqUHMsNCXO8X6cw13aC8aC+fdMbpQ3+qUpWmLymTCuPsHimqfWfW5zrR5h0a4M7xOl+gbbsLZSGoQEyz05U52+VyjuXDeZf3PUlJSkm5ZuC7sn3/f65+qbWqyUrrt0WX9Ozcj8viJtPA3VosdAvFATQ2ijoW+3Mkun2u4C+ct/WhX2K85um+e/nLdAGW3Tgn7ew4dk6Y9u8G253NLCn/ZcBNORVKDqGpqoS+jE9NjX1rv3B2YE5GdFnALd1jlrsUbmxXPqKI8rbrjEmW3Tg3zO070YNh14bqmVsAOB9O8o88Nu9DbGcNPiKpwbjgUJTqPnRZwC/dGe+CIr9nxpJ6RpPu+X6TJC04MRTV1u4nH+460hqmxafLhYpp3dNll+NbN6KlBVDX3LzuGpJzBTgu4NedGG0k8oaaWR/vnhKNkY4XOf/AdXfP4Kt36bJmueXyVzn/wnbB/X0K9l7ysNLVNT2Efsziyy/Ct29FTg6hq7l92FCU6Q3O2MYj17KjmLJwXaU9DYGr5Ux+U67evbYnZz2lMqGLo5s5OCjVNvnTz7rgtdpjorNiFPlHRU4OoimQcn6JE+wt3h/KDR2pa1LMQjnAXzmtpT0NykkcTzyuwZGf2aNcwBabJj+1/loZ2b6/kJI9tNtxMBPHehT6RkdQgqhrbybopFCXaVzg7lF/eL09TFq6PS/f66L55+smFBSGf9yi4pyHS4szGz+cTrxGLHo143QRHFeXp/V9drL/ffK7+eHV//f3mc/X+ry4moYkyOw3fuh1JDaIukpoEiaJEu2vsL/v/HT9Ar2yoiOvsqDtGF+rR8efUm4add1pPQ6zqUtqmSnOv7heTBCCeN8GGenEQXfHehT6RUVODmDh1HH935df67WtbdPDIsbjswIzYCVWfYdXsqNF9O2lkUV7IGp5Y1aW0Tz9D/9m8SiP7dIzaezkVN0F3ifcu9ImMpAYxc+py961SkylKdImGtjGwsns91LYK0S7OPPXn+Hw+LW26fjhiTr4Jsj1KffHehT6RMfyEuKAo0d3s2LPg5OLMcGqY7HgTbOlQn5txDYwPemoQN83dgRnOYceeBacXZwZugqcv1pZr08XaojXU52ZcA2OPpAZxFe4OzHAWO3av27H3qLmcchNszlBfouMaGFsMPwGIilh3rzd3Wna4a+vYsS7lVE6YneTkoT64Cz01AKImVj0LkeyZY8feI7dq3lBfZmyDQUKjpwZAVEW7Z6Ele+ZQnBkfbhjqgzvQUwPAtqIxLdspdSlO1pxCcX/t8XiHhwRCTw0A24pWrYYT6lKczKlT0OE+JDUAbMvp07ITCUN9sAOGnwDYFrUazsJQH6xGUgPAtuy4qB8axzossBLDTwBsi1oNAM1BUgPA1qjVABAuhp8A2B61GgDCQVIDwBEStVaj1m9I5oAwkdS0UKwuOG66kLnpvdiV1W1s9c93q4a2h2jbKkU3nFegqRf3oI2B0zgmqbn33nv12muvqaysTKmpqTp06JDVIUW0H42Vr2sFN70Xu7K6ja3++W4V2B7i9Flfh772ac5bn2r+h+V64Ipv0cbAKRxTKHzs2DFdeeWVmjx5stWhSGrZfjRWvK4V3PRe7MrqNrb657tVY9tDBBw66tMk2hgI4pikZvbs2brtttv0rW99y+pQmtyPRjqxH02tv7FLUvxe1wpuei92ZXUbx+Ln1/qNVm7fr8VlX2rl9v0Je340tT3Eqfg9Ak5yzPBTJGpqalRTU1P3dVVVlSTJ5/PJ5/NF/Lr/CHM/mpXb9mpIGIuCBWJZtf0/UX1dK0W7jUIJtF1LPk+namkbt7Ttov0Zv7Fpj+5Z+ol2V538nc3N9Oqu0b00sk/HiGKMlVifdxWHjoR/rEOuCQHRartav9Gazw9q7+EadcjwalDXdq6uMUrka50U/vt2dVJz//33a/bs2fUef/PNN5Wenh7x667d55GU3ORxb674h/ZvCf8vqHdWro3J61ohVm0USmlpaYtfw2mi1caRtl00P+MN+z168tNAx/HJG9PuqmpNfbZMP+rpV7/29jvnY3Xe/asyvLYNcMI14XQtabsN+z16cUeSDh07ea60TTW6ops9z5NoSsRrnSQdPXo0rOMsTWpmzZrVYNJxqn/+858aNGhQRK9/xx13aPr06XVfV1VVKT8/XyNGjFBmZmZErylJ7csP6JnP1jR53IgLhoTdU1NaWqqLhw7UM5+VRe11rRTtNgol0HbDhw9XSkpKxK/jRC1t45a2XbQ+41q/0f0PL5dU08CzHnkkvb4nXbdfe6Ft/hKP9XlX6zd6/uHlQb1WjXHCNSGgpW33xqY9mr9yQ71hz8pjHs3/NFlzr+5nu569aEjka510cqSlKZYmNVOnTtXVV1/d6DHdunWL+PW9Xq+8Xm+9x1NSUlp0Ugzt0SGs/WiG9ujQrIvwud3PjMnrWiFWbRRKSz9TJ4pWG0fadtH6+Wu272/05n1iGKtG6/992Hbr1MTqvEuRNOvyPg3OfjqVk64Jp4uk7Wr9Rve+vjVkHZdH0r2vb1Vx37Mc1x7hSsRrnaSw37OlhcI5OTnq1atXo/+lpdlv991Y7Ufjpn1u3PRe7MrqNo7Wz997OLyC2HCPc4vA9hBt0xu+mCfi71FTBdSBOq7V5QfiFxRsxTGzn3bu3KmysjLt3LlTtbW1KisrU1lZmb766itL4onVfjRu2ufGTe/Frqxu42j8/A4Z4f3hEu5xbjKqKE9r7xqu2y45W21bBSc3ifh7RAKMpjimUPg3v/mNnn766bqvzznnHEnSu+++q2HDhlkSU6z2o3HTPjduei92ZXUbt/TnDy7IDmsYa7BDakaiLTnJo1sv6ampF5+d8L9HJMBoimOSmqeeekpPPfWU1WHUE6v9aNy0z42b3otdWd3GLfn5gWGsyQvWySMFJTaJOMQSitWfsR2QAKMpjhl+AuBeVg+jwRmsriOD/TmmpwaAu1k9jBYLbPQZfYEE+PT9xnLZbwwiqQFgI24aYmGjz9hxYwKM6CCpAYAoC7XDdmCjT4bUWs5NCTCih5oaAIgiqzcaBRIZSQ0ARBELxAHWIakBgChigTjAOiQ1ABBFLBAHWIdCYQCIIhaIA1P5rUNSAwBRxArJLefkpICp/NYiqQGAKGOBuMg5OSlgKr/1SGoAIAZYIK75nJwUNDWV36MTU/mHF+ZyDsQQSQ0AxAgLxIXP6UlBc6byc07EDrOfAACWc/r6PkzltweSGgCA5ZyeFDCV3x5IagAAlnN6UhCYyh9qYMyjEwXPTOWPLZIaAIDlnJ4UBKbyS6r3HpjKHz8kNQAAy7khKQhM5c/NCu5Nys1Ks/XMLTdh9hMAwBbcsL4PU/mtRVIDALANNyQFTOW3DkkNAMBWSAoQKWpqAACAK5DUAAAAVyCpAQAArkBSAwAAXIFCYQBAs9T6jaNnJ8G9SGoAAGEr2VhRbx2ZPAetIwN3Y/gJABCWko0VmrxgXb3dtHdXVmvygnUq2VhhUWTACSQ1AIAm1fqNZi/ZLNPAc4HHZi/ZrFp/Q0cA8UFSAwBo0uryA/V6aE5lJFVUVmt1+YH4BQWchqQGANCkvYdDJzSRHAfEAkkNAKBJHTLSmj6oGccBsUBSAwBo0uCCbOVlpSnUxG2PTsyCGlyQHc+wgCAkNQCAJiUneTRzTKEk1UtsAl/PHFPIejWwFEkNACAso4ryNO+6AcrNCh5iys1K07zrBrBODSzH4nsAgLCNKsrT8MJcVhSGLZHUAACaJTnJo6Hd21sdBlAPw08AAMAVSGoAAIArkNQAAABXIKkBAACuQFIDAABcgaQGAAC4AkkNAABwBZIaAADgCiQ1AADAFRJqRWFjjCSpqqrK4kiC+Xw+HT16VFVVVUpJSbE6HEeh7SJH20WOtoscbReZRG+3wH07cB8PJaGSmsOHD0uS8vPzLY4EAAA01+HDh5WVlRXyeY9pKu1xEb/fr127dikjI0Mej302X6uqqlJ+fr6++OILZWZmWh2Oo9B2kaPtIkfbRY62i0yit5sxRocPH1anTp2UlBS6ciahemqSkpLUuXNnq8MIKTMzMyFP1mig7SJH20WOtoscbReZRG63xnpoAigUBgAArkBSAwAAXIGkxga8Xq9mzpwpr9drdSiOQ9tFjraLHG0XOdouMrRbeBKqUBgAALgXPTUAAMAVSGoAAIArkNQAAABXIKkBAACuQFJjM59++qnGjh2rnJwcZWZm6rzzztO7775rdViO8dprr2nIkCFq1aqVcnJydMUVV1gdkqPU1NSof//+8ng8Kisrszoc29uxY4duvPFGFRQUqFWrVurevbtmzpypY8eOWR2aLT366KMqKChQWlqaBg4cqBUrVlgdku3df//9+va3v62MjAx16NBB48aN09atW60Oy7ZIamzm0ksv1fHjx/XOO+9o7dq16t+/vy677DLt3r3b6tBs74UXXtD111+vG264QRs2bNAHH3yg8ePHWx2Wo9x+++3q1KmT1WE4xieffCK/36/HHntMmzZt0pw5c/SXv/xFd955p9Wh2c6iRYv0s5/9TDNmzND69et1wQUXqLi4WDt37rQ6NFtbtmyZpkyZolWrVqm0tFTHjx/XiBEjdOTIEatDsycD2/jPf/5jJJnly5fXPVZVVWUkmbfeesvCyOzP5/OZs846yzzxxBNWh+JYS5cuNb169TKbNm0yksz69eutDsmRfve735mCggKrw7CdwYMHm0mTJgU91qtXL/PrX//aooicae/evUaSWbZsmdWh2BI9NTbSvn179e7dW88884yOHDmi48eP67HHHlPHjh01cOBAq8OztXXr1unLL79UUlKSzjnnHOXl5am4uFibNm2yOjRH2LNnj26++Wb97W9/U3p6utXhOFplZaWys7OtDsNWjh07prVr12rEiBFBj48YMUIffvihRVE5U2VlpSRxjoVAUmMjHo9HpaWlWr9+vTIyMpSWlqY5c+aopKREbdu2tTo8W/vXv/4lSZo1a5buuusuvfrqq2rXrp2++93v6sCBAxZHZ2/GGE2cOFGTJk3SoEGDrA7H0bZv3665c+dq0qRJVodiK/v27VNtba06duwY9HjHjh0ZWm8GY4ymT5+u888/X0VFRVaHY0skNXEwa9YseTyeRv9bs2aNjDG65ZZb1KFDB61YsUKrV6/W2LFjddlll6miosLqt2GJcNvO7/dLkmbMmKEf/OAHGjhwoObPny+Px6PnnnvO4ndhjXDbbu7cuaqqqtIdd9xhdci2EW7bnWrXrl0aNWqUrrzySt10000WRW5vHo8n6GtjTL3HENrUqVP10Ucf6e9//7vVodgW2yTEwb59+7Rv375Gj+nWrZs++OADjRgxQgcPHgzaWv7ss8/WjTfeqF//+texDtV2wm27lStX6uKLL9aKFSt0/vnn1z03ZMgQXXLJJbr33ntjHarthNt2V199tZYsWRJ0c6mtrVVycrKuvfZaPf3007EO1XbCbbu0tDRJJxKaiy66SEOGDNFTTz2lpCT+XjzVsWPHlJ6erueee07f//736x6/9dZbVVZWpmXLllkYnTNMmzZNL7/8spYvX66CggKrw7GtM6wOIBHk5OQoJyenyeOOHj0qSfUuiElJSXU9EYkm3LYbOHCgvF6vtm7dWpfU+Hw+7dixQ127do11mLYUbtv96U9/0j333FP39a5duzRy5EgtWrRIQ4YMiWWIthVu20nSl19+qYsuuqiud5CEpr7U1FQNHDhQpaWlQUlNaWmpxo4da2Fk9meM0bRp0/TSSy/pvffeI6FpAkmNjQwdOlTt2rXThAkT9Jvf/EatWrXS448/rvLycl166aVWh2drmZmZmjRpkmbOnKn8/Hx17dpVDz30kCTpyiuvtDg6e+vSpUvQ123atJEkde/eXZ07d7YiJMfYtWuXhg0bpi5duuj3v/+9/vOf/9Q9l5uba2Fk9jN9+nRdf/31GjRokIYOHar/+7//086dO6k/asKUKVO0cOFCLV68WBkZGXU1SFlZWWrVqpXF0dkPSY2N5OTkqKSkRDNmzNDFF18sn8+nPn36aPHixerXr5/V4dneQw89pDPOOEPXX3+9vv76aw0ZMkTvvPOO2rVrZ3VocKk333xT27Zt07Zt2+olgIzsB7vqqqu0f/9+3X333aqoqFBRUZGWLl2asD2p4Zo3b54kadiwYUGPz58/XxMnTox/QDZHTQ0AAHAFBn8BAIArkNQAAABXIKkBAACuQFIDAABcgaQGAAC4AkkNAABwBZIaAADgCiQ1AADAFUhqAATxeDx6+eWXrQ7DFt577z15PB4dOnRIkvTUU0+pbdu2lsYEIDSSGiDBTJw4UePGjQv5fEVFhYqLi+MXkINcddVV+vTTT8M+vlu3bnrkkUdiFxCAIOz9BCCIHTZiNMaotrZWZ5zR8ktUNF+rVatWlmwieOzYMaWmpsb95wJOQ08NgCCnDj/t2LFDHo9HL774oi666CKlp6erX79+WrlyZdD3fPjhh7rwwgvVqlUr5efn66c//amOHDlS9/yCBQs0aNAgZWRkKDc3V+PHj9fevXvrng8M87zxxhsaNGiQvF6vVqxYUS+2QDzPPvusvvOd7ygtLU19+vTRe++91+RrGWP0u9/9Tt/4xjfUqlUr9evXT88//3zQ6y9dulQ9e/ZUq1atdNFFF2nHjh1Bzzc0/PTKK69o0KBBSktLU05Ojq644gpJJzYg/Pzzz3XbbbfJ4/HI4/HUfc8LL7ygPn36yOv1qlu3bnr44YeDXrNbt2665557NHHiRGVlZenmm29u+MMCEMwASCgTJkwwY8eODfm8JPPSSy8ZY4wpLy83kkyvXr3Mq6++arZu3Wp++MMfmq5duxqfz2eMMeajjz4ybdq0MXPmzDGffvqp+eCDD8w555xjJk6cWPeaf/3rX83SpUvN9u3bzcqVK825555riouL655/9913jSTTt29f8+abb5pt27aZffv21YstEE/nzp3N888/bzZv3mxuuukmk5GRUXd8qNe68847Ta9evUxJSYnZvn27mT9/vvF6vea9994zxhizc+dO4/V6za233mo++eQTs2DBAtOxY0cjyRw8eNAYY8z8+fNNVlZWXTyvvvqqSU5ONr/5zW/M5s2bTVlZmbn33nuNMcbs37/fdO7c2dx9992moqLCVFRUGGOMWbNmjUlKSjJ333232bp1q5k/f75p1aqVmT9/ft3rdu3a1WRmZpqHHnrIfPbZZ+azzz4L78MFEhxJDZBgIklqnnjiibrnN23aZCSZLVu2GGOMuf76682Pf/zjoNdYsWKFSUpKMl9//XWDP2P16tVGkjl8+LAx5mQi8vLLLzcaeyCeBx54oO4xn89nOnfubB588MGQr/XVV1+ZtLQ08+GHHwa93o033miuueYaY4wxd9xxh+ndu7fx+/11z//qV79qNKkZOnSoufbaa0PG27VrVzNnzpygx8aPH2+GDx8e9Ngvf/lLU1hYGPR948aNa6QlADSE4ScATerbt2/dv/Py8iSpbvho7dq1euqpp9SmTZu6/0aOHCm/36/y8nJJ0vr16zV27Fh17dpVGRkZGjZsmCRp586dQT9n0KBBYcUzdOjQun+fccYZGjRokLZs2RLytTZv3qzq6moNHz48KM5nnnlG27dvlyRt2bJF5557btAw0ak/pyFlZWX63ve+F1bMAVu2bNF5550X9Nh5552nzz77TLW1tQ3GDyA8FAoDaFJKSkrdvwM3fb/fX/f/n/zkJ/rpT39a7/u6dOmiI0eOaMSIERoxYoQWLFigM888Uzt37tTIkSN17NixoONbt24dcYynJiOnv1Yg1tdee01nnXVW0HFer1fSiYLi5oqkaNgYUy/Whn52S9oCSFQkNQBaZMCAAdq0aZN69OjR4PMff/yx9u3bpwceeED5+fmSpDVr1rToZ65atUoXXnihJOn48eNau3atpk6dGvL4wsJCeb1e7dy5U9/97ndDHnP6+jyrVq1qNI6+ffvq7bff1g033NDg86mpqUG9L4Gf8/777wc99uGHH6pnz55KTk5u9OcBaBxJDZCAKisrVVZWFvRYdna2unTp0uzX+tWvfqVzzz1XU6ZM0c0336zWrVtry5YtKi0t1dy5c9WlSxelpqZq7ty5mjRpkjZu3Kjf/va3LYr/f//3f3X22Werd+/emjNnjg4ePKgf/ehHIY/PyMjQL37xC912223y+/06//zzVVVVpQ8//FBt2rTRhAkTNGnSJD388MOaPn26fvKTn9QNqzVm5syZ+t73vqfu3bvr6quv1vHjx/X666/r9ttvl3RiFtPy5ct19dVXy+v1KicnRz//+c/17W9/W7/97W911VVXaeXKlfrzn/+sRx99tEVtAkDMfgISzYQJE4ykev9NmDDBGNNwofD69evrvv/gwYNGknn33XfrHlu9erUZPny4adOmjWndurXp27dv3SwgY4xZuHCh6datm/F6vWbo0KHmlVdeCXrdQHFvoCA3lEA8CxcuNEOGDDGpqammd+/e5u233647JtRr+f1+88c//tF885vfNCkpKebMM880I0eONMuWLas7ZsmSJaZHjx7G6/WaCy64wDz55JONFgobY8wLL7xg+vfvb1JTU01OTo654oor6p5buXKl6du3r/F6vebUy+3zzz9vCgsLTUpKiunSpYt56KGHgl6zoQJjAE3zGBPBQDIAWGDHjh0qKCjQ+vXr1b9/f6vDAWAzzH4CAACuQFIDAABcgeEnAADgCvTUAAAAVyCpAQAArkBSAwAAXIGkBgAAuAJJDQAAcAWSGgAA4AokNQAAwBVIagAAgCv8/3r2oxpMLfMpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.grid(True)\n",
    "plt.plot(result2.predict(linear=True), result2.resid_pearson, \"o\")\n",
    "plt.xlabel(\"Linear predictor\")\n",
    "plt.ylabel(\"Residual\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
