{ "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": { "iopub.execute_input": "2026-07-28T14:15:34.212335Z", "iopub.status.busy": "2026-07-28T14:15:34.212196Z", "iopub.status.idle": "2026-07-28T14:15:37.113170Z", "shell.execute_reply": "2026-07-28T14:15:37.112575Z" } }, "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": { "iopub.execute_input": "2026-07-28T14:15:37.114923Z", "iopub.status.busy": "2026-07-28T14:15:37.114663Z", "iopub.status.idle": "2026-07-28T14:15:37.123751Z", "shell.execute_reply": "2026-07-28T14:15:37.120627Z" } }, "outputs": [], "source": [ "raw = StringIO(\"\"\"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\"\"\")" ] }, { "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": { "iopub.execute_input": "2026-07-28T14:15:37.126323Z", "iopub.status.busy": "2026-07-28T14:15:37.125993Z", "iopub.status.idle": "2026-07-28T14:15:37.142533Z", "shell.execute_reply": "2026-07-28T14:15:37.141480Z" } }, "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": { "iopub.execute_input": "2026-07-28T14:15:37.144538Z", "iopub.status.busy": "2026-07-28T14:15:37.143943Z", "iopub.status.idle": "2026-07-28T14:15:37.189833Z", "shell.execute_reply": "2026-07-28T14:15:37.189422Z" } }, "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: Tue, 28 Jul 2026 Deviance: 6.1260\n", "Time: 14:15:37 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": { "iopub.execute_input": "2026-07-28T14:15:37.191674Z", "iopub.status.busy": "2026-07-28T14:15:37.191498Z", "iopub.status.idle": "2026-07-28T14:15:37.488822Z", "shell.execute_reply": "2026-07-28T14:15:37.488411Z" }, "lines_to_next_cell": 1 }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-packages/statsmodels/base/model.py:1315: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n", " predict_results = self.model.predict(self.params, exog, *args, **kwargs)\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAGxCAYAAABhi7IUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWQRJREFUeJzt3XlcVOX+B/DPgDCIsgiKoKAimoaKO+SSa6hoZlqGqHWtTNN2uy20GZV567bY1bKy5Xe9VqZlmTdFccml3BUIsTLFFUQFWRSBgTm/P7gzOTAznJk5M3POnM/79br35Zx5zpnnPAHnO8/zfZ5HIwiCACIiIiKV8nJ3BYiIiIjcicEQERERqRqDISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVWMwRERERKrWxN0VUAK9Xo/8/HwEBARAo9G4uzpEREQkgiAIKC8vR5s2beDlZbn/h8GQCPn5+YiKinJ3NYiIiMgOZ86cQWRkpMX3GQyJEBAQAKCuMQMDAyW/vk6nw6ZNmzBq1Cj4+PhIfn21YDtKg+3oOLahNNiO0lBzO5aVlSEqKsr4HLeEwZAIhqGxwMBApwVD/v7+CAwMVN0PqpTYjtJgOzqObSgNtqM02I5oNMVFkQnUFy5cwJ49e1BWVib6nPz8fGRmZuLKlStOrBkREREpjaKCocOHD2PKlCno1q0bBgwYgEOHDjV6TlVVFZKTkxETE4PJkyejdevWWLp0qQtqS0REREqguGBowoQJ2Lt3r+hzXn31VezatQt//vknjh07huXLl+Ohhx7C/v37nVhTIiIiUgpF5Qzdd999AICzZ8+KPufTTz/FrFmz0LZtWwDAHXfcgdjYWHz22Wfo37+/U+pJREREyqGoniFb5efn4/z58w2CnoSEBBw+fNhNtSIiIiI5UVTPkK2Ki4sBACEhISbHQ0NDUVRUZPG8qqoqVFVVGV8bErV1Oh10Op3k9TRc0xnXVhO2ozTYjo5jG0qD7SgNNbej2Hv26GDIMIXw+sAGAK5du2Z1euHChQuRlpbW4PimTZvg7+8vbSWvk5GR4bRrqwnbURpsR8exDaXBdpSGGtuxoqJCVDmPDoaioqLg5eWF/Px8k+P5+flo3769xfNSU1Mxb94842vDok2jRo1y2jpDGRkZSExMVO0aEFJgO0qD7eg4tqE02I7SUHM7il2Cx+OCoRMnTuDKlSuIi4uDv78/brrpJqxbtw7Tpk0DUBclbt68GS+++KLFa2i1Wmi12gbHfXx8nPqD5OzrqwXbURpsR8exDaXBdpSGGttR7P0qKhi6ePEijh8/josXLwIAcnNz4efnh8jISOOeI6+//jr27NmDnJwcAMBrr72GUaNG4cUXX8SAAQPwr3/9CyEhIZg1a5bb7oOIyKBWL2BfXjEulFciLMAP8dEh8PbihtBErqSoYGjfvn149dVXAdTNCFu+fDmWL1+OmTNnYubMmQCAmJgYXLt2zXjO8OHDsWXLFixZsgQ7duxAjx498NlnnzlluIuIyBbpOQVIW5eLgtJK47GIID/MHx+LMd0j3FgzInVRVDA0btw4jBs3zmqZ1NTUBseGDBmCIUOGOKtaREQ2S88pwJwVhyDUO36+tBJzVhzC0ul9GBARuYhHrzNERCRHtXoBaetyGwRCAIzH0tblolZvrgQRSY3BEBGRi+3LKzYZGqtPAFBQWol9ecWuqxSRijEYIiJysQvllgMhe8oRkWMYDBERuVhYgJ+k5YjIMQyGiIhcLD46BBFBfrA0gV6Dulll8dEhFkoQkZQYDBERuZi3lwbzx8cCQIOAyPB6/vhYrjdE5CIMhoiI3GBM9wgsnd4H4UGmQ2HhQX6cVk/kYopaZ4iIyJOM6R6BxNhwrkBN5GYMhojIIdxOwjHeXhoMiAl1dzWIVI3BEBHZjdtJEJEnYM4QEdnFsJ1E/cUDDdtJpOcUuKlmRES2YTBERDbjdhJE5EkYDBGRzbidBBF5EgZDRGQzbidBRJ6EwRAR2YzbSRCRJ2EwREQ243YSRORJGAwRkc24nQQReRIGQ0RkF24nQUSegosuEpHduJ0EEXkCBkNE5BBuJ0FESsdhMiIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVWMwRERERKrGYIiIiIhUjcEQERERqRqDISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVVNcMPTxxx+jW7duaNmyJYYPH46DBw9aLV9WVobHH38cN9xwA1q1aoW+ffti2bJlLqotERERyZ2igqEvvvgCjz76KF588UUcPHgQ3bp1w8iRI5Gfn2/xnNmzZ2PDhg344osvkJubi8cffxxz587FqlWrXFhzIiIikitFBUMLFy7EvffeiylTpqB9+/b417/+haZNm2Lp0qUWz9m9ezdSUlLQv39/tGrVCnfffTdiY2Oxe/duF9aciIiI5EoxwVBJSQmOHDmCkSNHGo95eXlhxIgR2LVrl8XzJk6ciHXr1uHcuXMQBAFbt27FiRMncNttt7mi2kRERCRzTdxdAbEMQ2FhYWEmx8PCwqzmDb399tuYNm0aIiMj0aRJE3h5eeHjjz/G8OHDLZ5TVVWFqqoq4+uysjIAgE6ng06nc+Q2zDJc0xnXVhO2ozTYjo5jG0qD7SgNNbej2HtWTDBk4OXl1eC1IAgWy8+aNQuZmZn45ZdfEBMTg02bNmHWrFkICQnB+PHjzZ6zcOFCpKWlNTi+adMm+Pv7O3YDVmRkZDjt2mrCdpQG29FxbENpsB2locZ2rKioEFVOI1iLJGTk0qVLaNWqFdasWYOJEycaj99zzz3Iy8vDzp07G5xTXFyM0NBQfP3117jrrruMx6dOnYpz585h+/btZj/LXM9QVFQULl26hMDAQAnvqo5Op0NGRgYSExPh4+Mj+fXVgu0oDbaj49iG0mA7SkPN7VhWVoaWLVuitLTU6vNbMT1DLVu2RKdOnbBjxw6TYGj79u0mgY45vr6+DV5biwG1Wi20Wm2D4z4+Pk79QXL29dWC7SgNtqPj2IbSYDtKQ43tKPZ+FZNADQCPP/44Pv30U2zfvh2VlZV49dVXceHCBTz44IPGMg8//DDi4+MBACEhIRg4cCBef/11nDx5Enq9Hps3b8bq1astDpERERGRuiimZwgAHnroIRQVFWHixIkoLS1F586d8cMPPyAmJsZYprKy0mSMcNWqVZg3bx569+6Nq1evIiwsDE8//TSefPJJd9wCERERyYyigiEAeOmll/DSSy9Bp9OZ7f56//33odfrja/btm2Lr7/+GgBQU1ODJk0Ud8tERETkRIqNDCyNA5rL9TFgIERERET1KSpniIiIiEhqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFSNwRARERGpGhfeISLyYLV6AfvyinGhvBJhAX6Ijw6Bt5fG3dUikhUGQ0REHio9pwBp63JRUFppPBYR5If542MxpnuEG2tGJC8cJiMi8kDpOQWYs+KQSSAEAOdLKzFnxSGk5xS4qWZE8sNgiIjIw9TqBaSty4Vg5j3DsbR1uajVmytBpD4MhoiIPMy+vOIGPULXEwAUlFZiX16x6ypFJGMMhoiIPMyFcsuBkD3liDwdgyEiIg8TFuAnaTkiT8dgiIjIw8RHhyAiyA+WJtBrUDerLD46xJXVIpItBkNERB7G20uD+eNjAaBBQGR4PX98LNcbIvofBkNERB5oTPcILJ3eB+FBpkNh4UF+WDq9D9cZIroOF10kIvJQY7pHIDE2nCtQEzWCwRARkQfz9tJgQEyou6tBJGscJiMiIiJVYzBEREREqsZgiIiIiFSNwRARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFSNwRARERGpGoMhIiIiUjXF7VqfnZ2Njz76CIWFhejRowcee+wxBAcHWz2npqYGy5cvx9atW+Hv74/7778fCQkJrqkwERERyZqieob27duHhIQE6PV6jB8/Hhs3bsTgwYNx7do1i+dUVlZixIgReOuttzB06FAMHToUTz/9NA4dOuTCmhMREZFcKapnKDU1FaNGjcLSpUsBABMmTEBkZCQ+/fRTPPzww2bPef3113H06FEcPXoULVu2BACkpKTg6tWrLqs3ERERyZdieoYqKyuxfft2TJo0yXgsODgYI0eORHp6usXzPv/8c9x9993GQAgAvLy8EBAQ4NT6EhERkTIopmfo9OnTqK2tRWRkpMnxqKgobNu2zew5xcXFOHv2LHr16oV//OMfOHjwINq0aYO7774b/fr1s/hZVVVVqKqqMr4uKysDAOh0Ouh0OgnuxpThms64tpqwHaXBdnQc21CcWr2AA6cu40J5FcICtOjXvgW8vTTG99mO0lBzO4q9Z8UEQ9XV1QAAf39/k+P+/v7G9+qrqKgAADz77LO48847MXnyZOzduxc33XQT1qxZg9tuu83seQsXLkRaWlqD45s2bWrw+VLKyMhw2rXVhO0oDbaj49iGlmUVabDmpBdKqv8KfoJ9BUzqoEfPUMGkLNtRGmpsR0Mc0BiNIAhC48Xc7+zZs4iKisJ///tfjBs3znj8/vvvR3Z2Nvbv39/gnCtXriAgIAATJ07EmjVrjMenTp2KvLw87N692+xnmesZioqKwqVLlxAYGCjhXdXR6XTIyMhAYmIifHx8JL++WrAdpcF2dBzb0LqNRwrxyMos1H/4GMKixVN6YnS31mxHiai5HcvKytCyZUuUlpZafX4rpmeobdu2CA0NRVZWlkkwlJmZid69e5s9p3nz5ujUqRNiYmJMjnfs2BF79+61+FlarRZarbbBcR8fH6f+IDn7+mrBdpQG29FxbMOGavUCFmz4vUEgBAAC6gKiBRt+R1JcWxhaju0oDTW2o9j7VUwCtUajwd13341PPvkERUVFAOq6/A4dOoR77rnHWG7x4sV46KGHjK/vu+8+/PDDD8a8nytXrmDt2rUYPHiwa2+AiIiwL68YBaWVFt8XABSUVmJfXrHrKkWqp5ieIQB49dVXkZWVhS5duuCGG27A4cOH8corr2DIkCHGMllZWdizZ4/x9d///ndkZmYiJiYGPXr0wJEjR9C1a1e888477rgFIiJVu1BuORBqWE76tAQicxQVDDVv3hxbt25FdnY2CgsL0a1bN7Rp08akzKOPPmrSU+Tj44Ovv/4ax44dw6lTp9CuXTvccMMNrq46EREBCAvwk7QckRQUFQwZxMXF2fxe586d0blzZ2dViYiIRIiPDkFEkB/Ol1aazRvSAAgP8kN8dAj0tTWurh6plGJyhoiISPm8vTSYPz4WwF+zxwwMr+ePjzVZb4jI2RgMERGRS43pHoGl0/sgPMh0KCw8yA9Lp/fBmO4RbqoZqZUih8mIiEjZxnSPQGJsOPblFeNCeSXCAuqGxtgjRO7AYIiIiNzC20uDATGh7q4GEYfJiIiISN0YDBEREZGqMRgiIiIiVWMwRERERKrGYIiIiIhUjcEQERERqRqDISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIiIiVWMwRERERKrGYIiIiIhUjcEQERERqRqDISIiIlK1Ju6uABERkSeo1QvYl1eMC+WVCAvwQ3x0CLy9NO6uFonAYIiIiGRP7oFGek4B0tbloqC00ngsIsgP88fHYkz3CDfWjMRgMERERLK28UghFmz4XbaBRnpOAeasOASh3vHzpZWYs+IQlk7vI4t6kmXMGSIiItnKKtLgkZVZJoEQ8FegkZ5T4Kaa1anVC0hbl9sgEAJgPJa2Lhe1enMlSC4YDBERkSzV6gWsOekl60BjX15xg0DtegKAgtJK7Msrdl2lyGYMhoiISJYOnLqMkmrLeUFyCDQulFsOhOwpR+7BYIiIiGTpQnmVyHLuCzTCAvwkLUfuwWCIiIhkKSxAK7Kc+wKN+OgQRAT5wVL/lQZ1yd7x0SGurBbZiMEQERHJUr/2LRDsK8g60PD20mD++Fhjfa5neD1/fKyslgGghhgMERGRLHl7aTCpgx6AvAONMd0jsHR6H4QHmfZQhQf5cVq9QnCdISIFk/tCdESO6hkqYPGUng3WGQqX0TpDQF1AlBgbzt9HhWIwRKRQXPGW1GJ0t9ZIimsr+0DD20uDATGh7q4G2YHBEJECccVbUhsGGuRMissZunbtGr777jt8+OGH+Pnnn20699ixY1i0aBG2b9/upNoROR9XvCUikpaigqHz58+jV69eeOmll/Dzzz9jwoQJmDFjhqhzKysrcccdd+D555/Hd99959yKEjkRV7wlIpKWoobJnn32WTRt2hR79uyBn58fsrOz0bt3b0ycOBETJkyweu7jjz+OoUOHwstLUfEfUQNc8ZaISFqKiQz0ej2+/fZbzJgxA35+ddMX4+LiMGjQIKxatcrqud9++y127tyJN9980xVVJXIqrnhLrlKrF7D7eBHWZp7D7uNFHHolj6WYnqHTp0/jypUr6Nq1q8nxrl27Yv/+/RbPO3XqFObOnYuNGzeiadOmoj6rqqoKVVV/LQNfVlYGANDpdNDpdHbU3jrDNZ1xbTVRSzv2jgxAeKAWhWVVZvOGNADCg7ToHRlgV1uopR2dyRPacOORQry2/jecL/vrb2F4oBYvjO2K0d1au6QOntCOcqDmdhR7z4oJhsrLywEAwcHBJsdbtGhhfK++mpoapKSk4KmnnkKvXr1Ef9bChQuRlpbW4PimTZvg7+8v+jq2ysjIcNq11UQN7Tg2XIPPygwdu9dPLxYgAEhqXYGN6Rsc+gw1tKOzKbUNs4o0+OyPhj9f58sq8fDKTNx3gx49Q13XS6TUdpQbNbZjRUWFqHKKCYYMQUj9wKesrMxigPLFF18gNzcXkydPxqJFiwAAFy9exOHDh7Fo0SI89thj0GgarlORmpqKefPmmXxGVFQURo0ahcDAQInu6C86nQ4ZGRlITEyEj4+P5NdXCzW141gAfcx8c48I8sPzSY59c1dTOzqLktuwVi9g4ds7AJjbJFUDDYANhf54etoQp6/zo+R2lBM1t6NhZKcxigmG2rVrB61Wi+PHjyMxMdF4/Pjx4+jcubPZc2JiYjBjxgycOnXKeKy6uhplZWU4efIkBEEwGwxptVpotQ03CPTx8XHqD5Kzr68WamnHW3tFOnUhOrW0ozMpsQ0PHC8yCbDrq5utWIXDZ8tdtu6PEttRjtTYjmLvVzHBkI+PD8aOHYsvv/wSs2bNgpeXF06fPo2ffvoJn376qbHc5s2bcf78eUyfPh2DBw/G4MGDTa7z008/YejQocaeIiIl40J0JDXOViQ1UkwwBABvvPEGBg4ciKSkJCQkJODLL7/EkCFDMHXqVGOZlStXYs+ePZg+fboba0pEpEycrUhqpKhgqHPnzsjJycEXX3yBwsJCvPzyy5gyZQq8vb2NZRITEy0OmwHAtGnT0KlTJ1dUl4jI5ZvpOvp58dEhiAjyw/nSSiuzFeuuS+QpFBUMAUDr1q1NkpvrS05Otnr+U089JXWViIjM2niksMFu687cTFeKzXu9vTSYPz4Wc1YcggYwCYgMIdX88bGy2ySVyBGKWXSRiEhJsoo0eGRlVoOtUwyb6abnFEj6eYbNe6X4vDHdI7B0eh+EB5kOhYUH+XETYPJIiusZIiKSu1q9gDUnvSxupqtB3Wa6ibHhkvSwNLZ5rz2fN6Z7BBJjw106xEfkLgyGiIgkduDUZZRUWw4art9MV4rZgLZs3mvL59WfrWjYnoPBEXkaBkNERBK7UG55nR7TctJMT3fFdHgp8pGI5Io5Q0REEgsLaLhoq/ly0kxPd/Z0eCnzkYjkiMEQEZHE+rVvgWBfAZYGkDSo61WRanq6YTq8Mz6vsXwkoC4fiTvak5IxGCIikpi3lwaTOugBoEGA4ozp6Ybp8M74PFvykYiUSnQwVFJSIvp/RERq1zNUwOIpPV02Pd1Z0+G5PQepgegE6hYtWoi+qCCwu5SIaHS31k7dTLc+Z0yH5/YcpAaig6H9+/c7sx5ERB7J1ZvpSv153J6D1EB0MNSvXz9n1oOIiGSI23OQGjCBmoiIrOL2HOTp7F50cfXq1Vi9ejVOnz6Nmpoak/cOHDjgcMWIiEg+uD0HeTK7eobee+89PPTQQ+jYsSP27t2LW2+9FWFhYTh48CCH04iIPJQhH2lCr7YYEBPKQIg8hl3B0AcffIDVq1fjH//4BwDg5Zdfxvr16/HPf/4T+fn5klaQiIiIyJnsCoZOnDiBAQMGAAC0Wi2uXLkCAJg5cyZ++uknySpHRERE5Gx2BUM1NTXw9fUFALRr1w4HDx4EABQUFMDb21u62hERERE5mcO71t9zzz2YMmUKRo8eja1bt2LChAlS1IuIiIjIJewKhs6cOWP89/PPP49WrVph9+7dmDt3Lh5//HGp6kZERETkdHYFQ5GRkcZ/azQazJ49G7Nnz5asUkRERESuYlcwdP78eavvh4eH21UZIpK/Wr3AtWaIyKPYFQxFRFhfbZQbtRJ5pvScAqSty0VB6V87lEcE+WH++FiuQkxEimVXMPTrr7+avNbr9Th27Biee+455gwReaj0nALMWXGowWad50srMWfFIW7LQESKZVcw1L179wbH4uLi0L59ezzyyCOYM2eOwxUjItcyDH8VlFzFiVINavUCfK57L21drtldywXUbdiZti4XibHhHDIjIsVxeGr99W688Ubk5ORIeUkicoGGw1/e+ObtHXj5tm4Y0z2iLki6bmisPgFAQWkl9uUVY0BMqEvqTEQkFcl2rdfpdFi0aBHatm0r1SWJyAUMw1/1g53CsirMWXEI6TkFuFBuORC6nthyRERyYlfPUHBwcINjV65cQdOmTfHll186WicichGxw19v3dlT1PXCAvykrB4RkUvYFQx98sknDY61aNECvXv3RkhIiMOVIiLXEDv8BU3drLHzpZVmAycNgPCgumn2RERKY1cwdOedd0pdDyJyA7HDWpeuVGH++FjMWXEIGsAkIDKkS88fH8vkaSJSJNHBUGMLLV6Piy4SKYPYYa2wAD8MiAnF0ul9GqwzFM51hohI4UQHQ40ttHg9LrpIpAzx0SE2DX+N6R6BxNhwrkBNRB5FdDB09OhR4783b96MBQsW4Pnnn0f//v0BAPv378eCBQvwwgsvSF9LInIKby+NzcNf3l4aTp8nIo8iOhjq2rWr8d/Tpk3DN998g0GDBhmPJSQkoFevXnjsscfw0EMPSVvLeoqKinDx4kV06NABfn7iuvnPnDmDJk2a2NTDRaQGY7pHWBj+0mL++G4c/iIij2fXOkNHjx5Ft27dGhzv1q0bfvvtN4crZYlOp8Pf/vY3tGnTBomJiQgLC8Nnn31m9Zz33nsP7dq1w6BBg9CzZ0906tQJW7ZscVodiZRoTPcI7HpmBL564Ca8M7kHHo6txbZ5QxgIEZEq2BUMtWvXDosXL25wfMmSJWjfvr3DlbJkwYIF2LRpE37//XecOXMGH374IR544AEcPnzYbPna2lrk5eVh9+7dOH36NM6fP49JkyZh4sSJuHDhgtPqSaREhuGv8XER6BwkMA+IyMPU6gXsPl6EtZnnsPt4EWr1zO81sGtq/aJFi3D77bdj9erV6N+/PwRBwIEDB3Ds2DGsXbtW6joaLVu2DDNnzkSHDh0AAFOnTsVrr72GTz75BO+//36D8t7e3li0aJHxtZeXF+bNm4d//vOfOHjwIJKSkpxWVyIiIrnIKtJg4ds7cL6syngsgjNBjezqGRozZgz+/PNPjB8/HsXFxbh8+TLGjx+PP//8E6NGjZK6jgCAgoIC5OfnIz4+3uT4TTfdhEOHDom+TnZ2NgA4tQeLiIhILjYeKcRnf3iZBEIAcL600rjljtrZvVFrZGQkFixYIGVdrCoqKgIAhIaazmIJDQ3FpUuXRF2jpKQEDz/8MG677TbExsZaLFdVVYWqqr9+aMrKygDU5SzpdDpbq94owzWdcW01YTtKg+3oOLahNNiOjqvVC3h1vflc3r+23DmCYZ1DPXJoXOzPjuhgqKSkBEDdvmSGf1tibu8yR/n4+ACASZBieG14z5qrV6/i1ltvRbNmzfDvf//batmFCxciLS2twfFNmzbB39/fhlrbJiMjw2nXVhO2ozTYjo5jG0qD7Wi/Y6UaFJZ546/FMkzVbblThSVfp6NzkOflEFVUVIgqJzoYatGiBYC6BRUN/7bEGYsutm3bFhqNBgUFpt15+fn5iIqKsnru1atXMXbsWJSXl2Pr1q2NBmupqamYN2+e8XVZWRmioqIwatQoBAYG2n0Pluh0OmRkZCAxMVFUYEfmsR2lwXZ0HNtQGmxHx63LLgByf220XMduvTA2zvNyhwwjO40RHQzt37/f7L9dpXnz5oiPj8f69esxdepUAEBlZSW2bNmC1NRUY7kzZ86goqICXbp0AVAXFY4bNw4lJSXYsmVLg2E2c7RaLbRabYPjPj4+Tv2FdPb11YLtKA22o+PYhtJgO9ovIriZ6HKe2MZi70l0MNSvXz+z/zaoqKhA06ZNodE4b8zxlVdewbhx43DjjTdiwIABWLRoEQICAjB79mxjmbS0NOzZswc5OTmoqanB+PHjcfToUXz55Zc4e/Yszp49C6BueYCQEO6wTUREnis+OgThgVqcL6uEuaGy+lvuqJVds8l+//13k96YJ598Es2bN0dUVJRxtpYzjBo1CuvXr8fevXuRmpqKsLAw7Nq1C0FBQcYy7dq1M66WfeXKFRQVFSEiIgJPPvkkZsyYYfzf9u3bnVZPIiIiOfD20uCFsXXPxPqhkOH1i+NuxL68YlWvP2TXbLInnngCjz/+OADg2LFj+PDDD/Hdd99h69atePrpp5Geni5lHU0kJiYiMTHR4vsvvfSS8d/BwcHIzMx0Wl2IiIjkbnS31rjvBj3Wn/c3mV4fHuSH23pG4NUfj5psxaPG9YfsCoZ+/vlnrF69GgCwceNGTJgwARMmTMCQIUPQsWNHSStIROpWqxewL68YF8orERZQ153viVOAiZypZ6iAp6cNweGz5cbfpctXq/HQl4dQvx/IsP7Q0ul9VBMQ2RUM+fr6orCwEB07dsSGDRswceJEAEBNTQ18fX0lrSARqVd6TkGDDWTV+K2VSAqGLXeAui8Zg9/Y2iAQAq5ffygXibHhqvjyYVfO0Lhx43DnnXdi1qxZ2LFjB8aPHw8A2LJlC0aOHClpBYlIndJzCjBnxSGTQAjgqrlEUtiXV9zgd+t6desPVWJfXrHrKuVGdgVDS5YswejRo3HlyhV8//33aN26NQBg27ZtePnll6WsHxGpUK1eQNq6XIvfWoG6b61qTPSUE278qVwXyi0HQvaUUzq7hsmaN2+OhQsXNjj+0UcfOVwhIhLPU/NpbPnWauj2J9fiEKayhQX4SVpO6ezem6ympgYHDx7EiRMnkJKSAgAoLS01meZORM7jyQ8jfmuVN8MQpqcn3nrqlw2gbv2hiCA/nC+tNNsDq7b1h+waJjt79iz69euHm2++2bgaNADcc889+PHHHyWrHBGZ5+n5NPzWKl9qGcJMzynA4De2ImXZHjy2MhMpy/Zg8BtbFf+7ZeDtpcH88XUblltaf2j++FiPCf4aY1cw9MQTT6Bv374N9vx4+umnzQ6fEZF0lPYwsievxPCt1dKfYQ3qesHU8q3VEVLn9agh8dbTv2wYjOkegaXT+yA8yPRLRXiQn8f07oll1zDZTz/9hNzcXPj5mTZgXFwcDhw4IEnFiMg8JeXT2DuUZ/jWOmfFIWgAk8BPjd9a7eWMoVRPH8Js7MuGp005H9M9Aomx4R47HCiWXT1D165dQ5MmdXHU9XuRnT9/Hv7+/tLUjIjMUsrDyNFv1/zW6hhn9W5IMYQp51loauj5qs+w/tCEXm0xICZUdYEQYGfP0JAhQ/Dxxx/jmWeeMQZD165dw9NPP40RI0ZIWkEiMqWEfBqpvl17wrdWdyThOrN3w9HEW7kn/ivlywZJy65g6K233sLQoUOxYcMGCIKAlJQU7NixAzqdDj///LPUdSSi6yhhFoiUQ3nXr5qrNBuPFGLBht9d/uB35lCqI0OYSpiFpoQvGyQ9u4bJYmNjkZ2djaFDhyIpKQklJSWYMWMGsrKy0LlzZ6nrSETXUcIsEH67BrKKNHhkZZZbknCd3f72DGEqJfGfyfvqZPc6QxEREUhLS2twPCcnB927d3eoUkRkneFhVH+4IVwmww1q/3Zdqxew5qSX25JwXdH+tg5hKiXxn8n76mRXMFRTUwO9Xm+yKeuZM2fw0ksvYfny5aitrZWsgkRknpzzacQM5bUO1EIvCFibeU5WdZfCgVOXUVJt+V6c/eB31VCqLUOYSuotlPuXDZKeTcHQxYsXcf/992Pjxo2oqanBxIkT8e9//xuffvopnn32WbRr1w6rV692Vl2JqB655tM09u1aAFBZo8e0T/Yaj8spiVYMa4nRF8qrRF3DWQ9+OfZuKK23UM5fNkh6NgVDzz77LI4cOYIFCxYAAD744AOMHj0av/76K9566y3MmjXLOOWeiNTN0rfrIH8flFToUFKhMykvpyTaxjQ2IyosQCvqOs588Mutd0MJif/1mfuy4clbdKiZTZHLxo0b8d///he9evUCAIwYMQJ9+/ZFeno6Ro8e7Yz6EZGC1f923bKZFk+uzgKga1BWKQvaiZkRNaxzKIJ9BZRWa9z64JdT74Yce6tsJfdlAch+Ns0mKygoQFxcnPF1z549AYBrCxGRRdcv6OblpcH5MuUuaCd2RhQATOqgB+D+GX9yWlBPyQtpqmWLDrWyqWdIr9fDy+uv+Mnb2xsA4OPjI22tiEhWpBoaUFISrTliZ0QdOHUZPUMFLJ7Ss8E6Q2pPwpVTb5VYatuiQ41sTvDp1KlTo8f+/PNP+2tERC53fbAT6t8E1y/1IuXQgNKSaOsTH8xVwRvA6G6tkRTXVlEPfleQa+K/JUpZFoDsZ1Mw9PzzzzurHkTkJuaCnWBfb/h0KESTJt6SrhislCRaSz1h4oM5LYr+92+lPfipIaX3aFLjbAqGXnvtNWfVg4jcwFIycEk18MjKLAT5+0g6NKCEJFprPWGJseGigrl+7Vtg41GXVZmcTGwQfPJShZNrQs5i13YcROR8lnb2lmrHb2t5EEDdLKj609+vZ2+ys5yTaBtLks3IPS/7rVBIevHRIQgPbHy5hJX7T7t9OxGyDxcFIpKh9dn5eGFtDoqv/hWMRAT54baeEfghq0CS/J3G8iDEsmdoQI5JtGKTZHc9M6LR9Xt0OstBJCmPt5cGKfHt8O7mY1bLMW9IuRgMEcnMwvW5+GhHXoPjBaWVZo/bm78jVX6DvcnOYnJpXLnAnS1JsnIM5si5OrRsJqoc84aUicEQkYyszy4wG/BYY2/+jqMztpyd7OzqBe5sTZJlYrS6KH0mJFnHnCEimajVC3hhbY5d59qTv2OY2WUpdNIAaOHvY/x3/fcA5+XHuGOBOz7syBoxvy8RMpgJSfZhMERuJVUysCfYl1eM4qvVDl3Dli56w8wuoGGwY5jntXBSD3zo4mRnsas8S/2zwocdWWPt94XJ88rHYTJyG+7zY0qKXANbey0sbeYZ7Au8Nqmn8b+DK/Nj3LXAnRKm/ZN7yW3zW5IOgyFyCzGbXbr7D4urd6d2dPjF3l6L+snAof5NcDF3D0Z3a20s48r8GHcucMeHHTWGyfOeicGQm1jb/sDTKWGfH3f0WhmGaeyd7v7iuBvtbq/rgx2dTof1blww0N25O3zYUWOYPO95GAy5gbXtD27tFenGmpkndQ+J3Pf5cVev1fXDNPbExi2aNb4onBLIYcsOT3vYubqXk0hpFJdAvXz5cvTt2xeRkZFISkpCdna2U85xFkuzZAzbHzhjlowj0nMKMPiNrUhZtgePrcxEyrI9GPzGVofqKed9ftyVvGtgGKaJCLK918NT1jdhomrjbJl44IzfYSJPo6hgaNWqVZg5cyYeeeQRbNmyBZGRkRg+fDgKCwslPcdZGtv+AHDug9ZWzpre7O5hEGts6bWSmuEBV1Wjx1t39sQXMxPw3pReeHHcjaLO96Qp33LessPdbAlu3LFEAZESKSoYWrBgAe69917MmDEDXbp0wYcffogmTZrggw8+kPQcZ3Hng9ZWzuwhkfMUZnf1WtV/wE37dC/+vjoL2iZemDEoWrbtZQtbl1EY0z0Cu54Zga8euAnvTemFrx64CbueGaH6QEhscOPuXk4iMeSyvIpicoZKS0uRnZ2NF154wXjM29sbI0aMwM6dOyU7x5nkPDxUnzPzeuQ8hdkdvVZicpTk2l5i2ZuQ7mm5O46wdeKB3HPziOS0vIpigqH8/HwAQOvWrU2Ot27dGpmZmZKdAwBVVVWoqqoyvi4rKwNQN8vGkQ0YQ/3FNXdhaQXWHDyNsAAt+rVv4ZaHXEHJVdHldLpAm68/sktLLJ7SE6+t/w3ny/5q6/AgLZ5P6oqRXVra3NaG8o78N+odGYDwQC0Ky6qsJO9q0TsyQJLNOGv1Al7+4UgjD7gj2DZviOTtZYkU7Xi9jUcK8cjKLIvB3uIpPU2m8XsCqdsQAPaKDG52/3kBCdEhTv8ddoXKqmocK9Xg+8NnERHs77a/h0rnjJ9HR7nq74LYe1ZMMCQIdU3m7e1tcrxJkybQ6/WSnQMACxcuRFpaWoPjmzZtgr+/v031vp5eqJs1VlINmFvzFxCgAfD6hj+MR4J9BUzqoEfPUMe6DvUCcLxMgzIdEOgDxAQKsPY35USpBoC35QKGckcysf7sYbvr9Uxs/XpdRe2pg1h/yu5LIiMjw/6TAYwN1+CzMsMI8vWNJEAA0CugAq+vSBfVjo05VqrB+TLL7Vz3gKvCkq/T0TlIcEp7WeJoOwJ1P3dph7z/9wfPtKGE//3/C2syoTtZ61A7ypUUbWhw8JK438lNO/ei6Kjgst9he4j5e5RVpMGak14oqfYGcnMBSPf3UK2k/Hl0hCv/LlRUVIgqp5hgqFWrVgCAixcvmhy/ePGi8T0pzgGA1NRUzJs3z/i6rKwMUVFRGDVqFAIDHfsG5dOhLhoGUC8irusHqP8rXlqtwed/eDsUJW88UoiF9XsUArV4YWxXi9es1Qv45u0djfaQPJw8RDbf1HQ6HTIyMpCYmAgfHx+7rzMWQJ8jhQ16YYL9fQBBgw1n//qm0Vg7NmZddgGQ+2uj5Tp264WxcbZ1G9fqBRw4dRkXyqts6mWUqh2But6Mkj0HrJTQoKQaaBV7ExJknvNkCynb0CA0rxjLj1lryzqjbk5AQnSIbH+Hxfw92nikEJ/vbthrIMXfQzVyxs+jI1z5d8EwstMYRQVD0dHR2LVrF26//Xbj8R07dmDSpEmSnQMAWq0WWm3DNVt8fHwc/kG6tVckmjTxbjBOWj8XxMAwVLJgw+9Iimtr8x+t9JwCs12RhWVVeGRllsWZOT4AXr6tWyN5Kt3gp/W1qT62sHdtFKn+OyXFtTV+/slLFVi0+Q+b27ExEcHNRJez5Z6kGIuXoh2LKmpEl5PDH2mpSdGGBgM6hYlaf2lApzB4e2lk8Ttcn5i/R4mx4Viw4Xen/D1UOyl/Hh3hyr8LYs9X1GyyRx99FJ988gn27t2LmpoavPnmm8jPz8fs2bONZZ544gncfPPNNp3javVnyTyXdMP/BsjMs3eWmaOzSdw5vVkOa6MYkndvjWuDlftPK2ZmnZymU8t5GQWlsWf9JTktUSD279GeE0WKmXVL9pHj3wXF9AwBwGOPPYYLFy7glltuQVVVFdq2bYs1a9bghhtuMJYpLy/H5cuXbTrHHa6fJbPm4GlR59g6y0yK2SSu2Jqgfg/Q5avVeOhL+exbpqSZdWJmHL38wxEE+Png0pUqp69GLIfVpD2JPXunyWV7EbG/R7uPF4m6nhxm3ZJ95Ph3QVHBkEajweuvv47XXnsNV69eRUBAQIMyixYtQk1NjU3nuFtYgLhtFGyNkqWayu/M6c3mhnO8NNaHDB3dt8zW4TdnL4kg5eagYh4458uqMO2TvcZjtgyf2dp2cl5GQansCW7ksESB+N8PcT2s7E1ULjn+XVBUMGTg5eVlMahp3ry5zee4W7/2LRDsK6C0umECNWB/lCzHrsjrWVpfx9pok6Nro5jdF66pD+4dFI2HR3Qy+8vninaU6tu7PQGZ2B43e/OQuBO89OQQ3NhK7O/HgI4t8e2hc7LqNSDpye3vgiKDIU/j7aXBpA56fP6Ht6RRshy7Ig2sb03SOHse+paCr5JrOry7+Q98/kse/jGpR4NfQle1oxQPOHsCsvo9buY4unmtXIZqyH3E/h7dFBMqu14Dcg45/V1QVAK1J+sZKmDxlJ6SJjrKecPLxoZzGmPrQ19M8FVSocODZhKM5dyO9TWWkG2JtaRUqbZ1MAR7E3q1xYCYUFm0F7mOLb9Hhl6D1oGmKQTcm87zyOXvAnuGZGR0t9Ym07mliJLl1hVpYG9+jb29MLYEX+ZykuTajvVZG4sXo+6/i+laWtzWgaRiy+/RmO4RGNY5FEu+TkfHbr0QEdyMvYnkNAyGZMYZuQBy6oo0sGc4x5FeGFuCL0sPdjm2ozmWHjhimPvvoqQ99ayxd90qkpYtv0feXhp0DhIwNi5CFuvjkOdiMKQScku4bCx/AKibVXb9yIsjvTBSzcSTWztakhgbjgCtD3afuARAg4ToEDz1TTYKyxrPe9LXmi6IJvdEfDHktCEkKef3iNSDwRC5hZiplUtSeqNFM60k3+QNwZfYnhI5P9gbY+7B/+0hP0zoFYGPd+Q1mpSqrzW9npwT8cVwNPnb07HHjIjBELmRK/Nwrg++rOXRyP3B3hhrD/6Pd+Rh1pBo/JBVYFN7y3FNELHELELp6LpVSsYeM6I6DIbIrVyZh2MIvp5d8ytKKnQN3pf7g70xYh78P2QVYPtTw3Hw1GWb2ttdCeSO9lqITf7+v5/zMGNQtCL/u9tLyT1m7M0iqTEYIrdzZf6AIfhasvUYPv/5JEquXbcDvcK/EYt98B88ddmu9nZ1ArkUvRZik7pf/fEo3v/pT7w2oTvGxrWxqZ7mHsxy5+4eM0eCGfZmkTMwGCLV8fbS4LFbbsDDIzp71LdLV8z6clXgKlWvhS25X8VXdZj75WHMPluC1LGxoutp7sH8fFIX0Z/rDu5cLsGRYEbJvVkkb1x0kVRLLot9ScUTZn0B0i3yCNi3COVHO/KwPrug0XKGB3P9oOJ8aSUeWZmFrCL5/jy5a7kEa202x8yCp9eT8ueCqD4GQ0QeorEHvwZ138DlPoxjS69FY6ytemzNi2tzrD5UxTyY15z0ku2D2R2Bs6PBjJQ/F0T1MRgi8hCObBtSqxew+3gR1maew968Yqub5Tqb1L0WhuTv+lvdWFN0tdrqQ1XMg7mkWoMDpy6L/kxHXf/fcPfxIquBmDsCZ0eDGU9Z/JPkiTlDRB7Enllf5nI4gn294dOhELf2inRJva/njF4LQ/L3//2ch1d/PCrqHGsPVfEP5ipR5Rxlax6OO5ZLcDSY8ZRhYJInBkNEHsaWWV+WElJLqoFHVmahSRNvyRJSxc4gctYij95eGswYFI33f/oTxVcbLq1Qn7WHqvgHs7bxQg6yN6nY1cslOBrMKH3xT5I3BkNEHkjMrC9rORyG/gGpplfb0nPhzF4Lby8NXpvQHXO/PGy1XGNDRGIezEG+Avq1b2FzHW3h6BR5Vy6X4Ggwo+TFP0n+mDNEpFKuSki1ZwaRpTyf8CA/sz0d1TV6fLrzBF5am4NPd55AdY3eYn3GxrXB7CHRFt/XoPGHqpj8rEkd9E5/MEvx39BVsyodyWkzsPXngkgs9gwRqZQrElId6bkQ22uxcH0ulu3MM0n6XrD+KB64OdriekGpY2PRM7IFXlibg+Kr1cbjtizeZ22Y6fmkLqg9dbDRazhKaUnFUgzNuXrxT1IHBkNEKuWKhFRHF/drbLhv4fpcfLQjr8FxvQDjcUsB0di4CIzu7thD1dKDWV9bg/WnRF/GbkpMKpYimOGu9yQ1BkNEKuWKhFRn9lxU1+ixbGfDQOh6y3bm4clRXeHbxHxGgBQPVXPX0Nc6dEnRlJpUzGCG5IY5Q0QqZX1BwrpHq6MJqc7sufjP7pONroekF+rKeSop8nCIiMEQkapZSkgN9gUWT+lpNodDLov7nSqukLScNbbcs6sxqZjIcRwmI1K5+jkcof5NcDF3D0Z3a92grJwW92sf4i9pOUuUsEs6k4qJHMOeISIymV6dEB0Cc89QezfZdFbPxd0DOpit5/W8NHXl7OXIxqKu5mkbDxO5EnuGiKhRclzcz7eJFx64OdrsbDKDB26Otpg83RhH75mIlIPBEBE1ytEp8oBzZhAZps3XX2fISwOr6wyJIcU9y5nY7VGI1IDBEBE1Ss6L+6WOjcWTo7riP7tP4lRxBdqH+OPuAR3s7hEykPM9O0oJeVBErsRgiIgaJffF/XybeOH+mztKek2537O97N3YlciTMYGaiBrlzCnycuWJ99xYHhRQlwclp6UDiFyBwRARNUqui/s5c/0fd96zs+7LVZvzEikNh8mISBQpNtmUkivyXtxxz868L0/OgyJyBIMhIhJNLov7uTLvxZX37Oz78tQ8KCJHKS4Y+u2337Bs2TIUFhaiR48eeOihh9C8eXOL5auqqrBy5Ur88ssvaNKkCQYPHozk5GR4eXGEkMgezt5ks7Ep3+5Y/8cVG4u64r6UurErkbMpKiI4dOgQ+vbti+LiYtx888345ptvcPPNN6Oqqspseb1ej9jYWPz000/o27cvOnfujKeeegqTJk2CIDBBkEhu0nMKMPiNrUhZtgePrcxEyrI9GPzGVpOVnj0178UV9yXX3C8id1NUz1BqaiqGDRuGzz//HABw5513IioqCp9//jkefPDBBuU1Gg127NiBtm3bGo8lJCRg4MCBOHjwIPr16+eyuhORdWKHiDw178VV9yW33C8iOVBMMFRVVYWtW7fi448/Nh4LDQ3FyJEjsX79eovB0PWBEABERkYCAC5fvuzcChORaLYMEXlq3osr70suuV9EcqGYYOj06dOoqalBu3btTI63a9cO27dvF32dxYsXIzg4GAkJCRbLVFVVmQy9lZWVAQB0Oh10Op2NNW+c4ZrOuLaasB2l4Y523CtyiGj3nxfQr30LhAdqUVhWZSXvRYvekQEO3UOtXsCBU5dxobwKYQFa9GvfQnSwYE8b9o4McMl9Xa9fu0AAgQAAfW0N9LWSXFYy/J2WhprbUew9uzUYev/997Ft2zarZZYsWYLw8HBjcOLv72/yfvPmzVFZKa7beNWqVXj77bfx1VdfITAw0GK5hQsXIi0trcHxTZs2Nfh8KWVkZDjt2mrCdpSGK9vx4CUNAO9Gy23auRdFRwWMDdfgszJDyuP1AYoAAUBS6wpsTN9gd32yijRYc9ILJdV/XTvYV8CkDnr0DBWfb2hrGzr7vpSKv9PSUGM7VlRUiCrn1mAoISEBrVu3tlomICAAABAUFASg4fBWUVERgoODG/2stWvX4u6778b777+Pu+66y2rZ1NRUzJs3z/i6rKwMUVFRGDVqlNUgyl46nQ4ZGRlITEyEj4+P5NdXC7ajNNzRjqF5xVh+7ECj5UbdnICE6BCMBdDnSCFeW/8bzpf91YsbEeSH55O6YnQ3639XrNl4pBCf785q0DtTWq3B5394Y/GUno1e3942dOZ9KRF/p6Wh5nY0jOw0xq3BUL9+/UQnMUdGRqJFixbIzs7G2LFjjcezs7MRFxdn9dx169YhOTkZ7777rtncovq0Wi20Wm2D4z4+Pk79QXL29dWC7SgNV7bjgE5hoqZ8D+gUZhyqurVXJJLi2kqa91KrF7Bgw+9Wc5cWbPgdSXFtRX2OPW3ojPtSOv5OS0ON7Sj2fhUztV6j0WDatGn49NNPUVJSAgDYvn079u/fj+nTpxvLffjhhya9Oj/++CMmT56Md955B3PnznV1tYlIBHunfBvW/5nQqy0GxIQ6HDDIZdq+1PdFRNYpJoEaABYsWIBDhw7hxhtvxI033oi9e/fiueeew4gRI4xlDhw4gD179gCo6x674447EBgYiK1bt2Lr1q3GcnPmzMHIkSNdfg9EZJ4cpnx76rR9IrJOUcFQYGAgdu3ahf3796OwsBDdu3dHdHS0SZk5c+YgOTkZAODn54cVK1aYvVZMTIzT60tEtnH3lG9PnbZPRNYpKhgC6obL4uPjLb7ft29f4799fX1x5513uqJaRCQRV2x9YQm3q3BMY1upEMmV4oIhIiJnMeQuzVlxCBrAJCCSYrsKTw4W0nMKGgxxRnBVa1IIBkNERNdxVu6SJwcLYrdSIZIrBkNERPVInbvkycGCLVupeEovGHkeBkNERGZIlbvUWLAAKDtYsGU5AnflghE1RjHrDBERKVFjwQJQFyws2XrMRTWSFpcjcK1avYDdx4uwNvMcdh8vQq1e/PYwZBl7hoiInEhsEPDu5mPoEh6guOEyLkfgOp6cd+Zu7BkiInIiW4KAtHW5ivumb1iOwNIAnwZ1D2wuR+AYQ95Z/V5GQ95Zek6Bm2rmGRgMERE5kSFYEMMVW31Izd6tVEg8sXlnSguk5YTBEBGRE10fLIihxNwaw3IE4fWCvvAgP0XPlJMLueyZ58mYM0RE5GRjukfgiVtuwLub/2i0rFJza9y9lYonY5K68zEYIiJygYdHdMJX+07hfFmV2fc9YasPd26l4smYpO58HCYjInIBby8NXr6tGzRgbg3ZhknqzsdgiIjIRZhbQ/ZgkrrzcZiMiMiFmFtD9nDWnnlUh8EQEZGLMbeG7MFA2nkYDBERESkEA2nnYM4QERERqRqDISIiIlI1BkNERESkaswZIiKyQ61eYCIrkYdgMEREZKP0nIIGU5wjOMWZSLE4TEZEZIP0nALMWXGowcaZBaWVeHDFIazPzpf8M2v1AnYfL8LazHPYfbyIu5MTSYw9Q0REItXqBaSty4W1UOThrw5j0V1xkn0me6GInI89Q0REIu3LK27QI1SfXgAe/TobWUWO5w9Z6oU6X1qJOSsOIT2nwOHPICIGQ0REol0otx4IXW/NSS+HhrOs9UIZjqWty+WQGZEEGAwREYkUFuDXeKH/KanW4MCpy3Z/VmO9UALq8pT25RXb/RlEVIfBEBGRSPHRIYgIEh8QXSivsvuzxPZC2dJbRUTmMRgiIhLJ20uD+eNjRZcPC9Da/Vlie6Fs6a0iIvMYDBER2WBM9wh8MLU3rK2vqAEQ7CugX/sWdn+OoRfK0sdoUDerLD46xO7PIKI6DIaIiGw0Nq4NlqT0MfueIXiZ1EHv0IrU1/dC1b+K4fX88bFc9ZpIAgyGiIjsMDYuAh9O79Mghyg8yA+Lp/REz1DHZ3mN6R6BpdP7INzMZyyd3ofrDBFJhIsuEhHZaUz3CCTGhjfYo0xfW4P1p5z7GewRIpKO4oKh6upqbNmyBYWFhejRowf69u0r+ty8vDxs3LgRcXFxGDhwoBNrSURq4e2lwYCYUJNj+lrnf4YUuNmsa7G95UtRwdDFixcxYsQIVFdXIy4uDk888QSSk5Px4YcfNnpuVVUV7rjjDvz222+YNWsWgyEiUjVu8+FabG95U1TOUGpqKgDg8OHDWL16NTZv3oyPP/4YP/74Y6PnPvXUU+jfvz9uuOEGZ1eTiEjWuM2Ha7G95U8xwZBer8eqVatw7733wt/fHwDQt29fDBw4ECtXrrR67g8//ICNGzfinXfecUVViYhki9t8uBbbWxkUM0x25swZlJeXIzbWdMGz2NhYHDhwwOJ5Z8+exezZs7Fu3To0a9ZM1GdVVVWhquqvlWPLysoAADqdDjqdzo7aW2e4pjOurSZsR2mwHR0n5zbcK3Kbj91/XkCCm9cwknM7iiWH9vaEdrSX2Ht2azC0bds2/P7771bLTJ06FYGBgcaAJDg42OT9Fi1aGN+rr7a2FtOmTcOjjz6Kfv36ia7XwoULkZaW1uD4pk2bjL1SzpCRkeG0a6sJ21EabEfHybEND17SAPButNymnXtRdFQevRVybEex5NTeSm5He1VUVIgq59Zg6PTp08jMzLRa5o477gAANG3aFABQXl5u8n5ZWZnFAOXLL79EVlaWSZJ1UVERfv31V3z44YeYPXs2NJqGmfypqamYN2+eyWdERUVh1KhRCAwMFH1/Yul0OmRkZCAxMRE+Pj6SX18t2I7SYDs6Ts5tGJpXjOXHLPemG4y6OUEWPUNybUex5NDentCO9rLUWVKfW4Ohv/3tb/jb3/4mqmz79u3h6+uLvLw8k+N5eXno1KmT2XOioqIwZcoUZGdnG49du3YNFy9eRGZmJgRBMBsMabVaaLUN9xTy8fFx6g+Ss6+vFmxHabAdHSfHNhzQKQwRQX44X1ppNo9Fg7pFHQd0CpPNtG9b2lFu09fl1N5y/Hl0NrH3q5icIR8fH4wZMwYrV67EAw88AI1Gg/z8fGzbtg0fffSRsdz27dtRWFiIu+66C8OGDcOwYcNMrrNnzx4MGzYMixYtcu0NEBHJgGGbjzkrDkEDmDyglb7Nhxynr3tye3sSxcwmA4A33ngDmZmZmDBhAl5//XWMGDECCQkJmD59urHMf/7zH7zyyiturCURkbx54jYfcp6+7ont7WkU0zMEAF27dkV2djaWL1+OwsJC/P3vf8c999yDJk3+uo1hw4ahffv2Fq8xefJkrjVERKrnSdt8NDZ9XYO66euJseFuuz9Pam9PpKhgCAAiIyPx3HPPWXz/+l4ic55//nmpq0REpEjO2ubD1faJnL6+L6/YrffrKe3tiRQ1TEZERFTfhXLLgZA95Uh9GAwREZGihQX4NV7IhnKkPgyGiIhI0eKjQxAR5AdL2Tca1M0qi3fzukkkXwyGiIhI0QzT1wE0CIg4fZ3EYDBERESKx+nr5AjFzSYjIiIyh9PXyV4MhoiIyGNw+jrZg8NkREREpGoMhoiIiEjVGAwRERGRqjEYIiIiIlVjAjURkQrV6gXOuiL6HwZDREQqk55TgLR1uSabm0YE+WH++FhZrcdTqxdwrFSDddkFiAhuxoCNnIbBEBGRiqTnFGDOikMQ6h0/X1qJOSsOyWaBwvScArz8wxGcL/MGcn8FIM+AjTwDc4aIiFSiVi8gbV1ug0AIgPFY2rpc1OrNlXAdQ8B2vqzK5LghYEvPKXBTzchTMRgiIlKJfXnFJkNj9QkACkorsS+v2HWVqkcpARt5FgZDREQqcaHcciBkTzlnUELARp6HwRARkUqEBfg1XsiGcs6ghICNPA+DISIilYiPDkFEkB8szcfSoC5JOT46xJXVMqGEgI08D4MhIiKV8PbSYP74WABoEBAZXs8fH+vW6etKCNjI8zAYIiJSkTHdI7B0eh+EB5n2rIQH+cliWr0SAjbyPFxniIhIZcZ0j0BibLhsV6A2BGx16wz9Nb0+nOsMkZMwGCIiUiFvLw0GxIS6uxoWjekegWGdQ7Hk63R07NaLK1CTUzEYIiIiIzntWebtpUHnIAFj4yLg4+PjljqQOjAYIiIiAMrZs4xIakygJiIi4xYY9Rc85BYYpAYMhoiIVI5bYJDaMRgiIlI5boFBasdgiIhI5bgFBqkdgyEiIpXjFhikdgyGiIhUjltgkNoxGCIiUjlugUFqp8hg6OrVqzh79ixqa2ttOu/SpUsoLy93Uq2IiJRL7nuWETmTohZdrK2txaOPPopPPvkEzZo1g7e3NxYvXowpU6ZYPW/79u145JFHcOrUKfj7+2PkyJH44IMPEBgY6KKaExHJn9z3LCNyFkX1DL355ptYtWoVMjMzUVxcjAULFmD69OnIycmxeM7hw4cxevRoTJkyBUVFRSgoKMD48eNx/PhxF9aciEgZDHuWTejVFgNiQhkIkSooKhhaunQpZs6ciRtvvBEAMGvWLHTs2BEff/yxxXNeeOEF9O/fH8899xyaNKnrCEtOTkbv3r1dUmciIiKSN8UEQ4WFhThz5gwGDBhgcnzQoEE4cOCA2XN0Oh22bNmCSZMmoba2FufOnYNOp3NFdYmIiEgh3JozdOHCBZSVlVkt0759e/j4+ODSpUsAgNDQUJP3W7ZsiZ9//tnsuZcuXUJVVRUKCgoQHR2NmpoaFBUVITk5GR988AGaN29u9ryqqipUVVUZXxvqqNPpnBJMGa7JQM0xbEdpsB0dxzaUBttRGmpuR7H37NZg6O2338a3335rtczWrVvRrl07eHnVdWLVv7Hq6mp4e3ubPVejqRvr/s9//oOdO3eiU6dOyMvLw5AhQ5CamorFixebPW/hwoVIS0trcHzTpk3w9/dv9L7slZGR4bRrqwnbURpsR8exDaXBdpSGGtuxoqJCVDmNIAiK2HmvvLwcgYGB+Oqrr0xmjyUnJ6O4uNjsf+Ta2lo0b94cDz74IN59913j8dTUVHz77bf4448/zH6WuZ6hqKgoXLp0ySkz0HQ6HTIyMpCYmAgfHx/Jr68WbEdpsB0dxzaUBttRGmpux7KyMrRs2RKlpaVWn9+KmVofEBCAvn37YuPGjcZgqLq6Gps3b8a8efOM5S5evIjKykpERUXB29sbQ4cOxeXLl02uVVxcbLVRtFottFptg+M+Pj5O/UFy9vXVgu0oDbaj49iG0mA7SkON7Sj2fhUTDAHA/PnzMWnSJPTq1QsDBgzAO++8A19fXzz44IPGMqmpqdizZ49xuv38+fNxyy23ICEhAYMGDcLevXvx73//G0uXLnXXbRAREZGMKGY2GQCMHz8eq1evxtq1a3HvvfcCAHbs2GGSVB0WFoZ27doZXw8YMAAbNmzAf//7X0ydOhXfffcdvvzyS+P5REREpG6K6hkCgNtvvx233367xfdff/31BseGDBmCIUOG2P2ZhrSqxma+2Uun06GiogJlZWWq68KUEttRGmxHx7ENpcF2lIaa29Hw3G4sPVpxwZA7GPYzi4qKcnNNiIiIyFbl5eUICgqy+L5iZpO5k16vR35+PgICAozT9aVkmK125swZ7pfmALajNNiOjmMbSoPtKA01t6MgCCgvL0ebNm2MS/SYw54hEby8vBAZGen0zwkMDFTdD6ozsB2lwXZ0HNtQGmxHaai1Ha31CBkoKoGaiIiISGoMhoiIiEjVGAzJgFarxfz5880u9EjisR2lwXZ0HNtQGmxHabAdG8cEaiIiIlI19gwRERGRqjEYIiIiIlVjMERERESqxmBIpgRBwIkTJ3Dw4EFcuHDB3dVRtJKSEhw8eNBp26mohSAI+OWXX3D48GF3V0WRdDodcnJykJeXB71e7+7qyF51dTUOHTqE33//3d1VUSzDc+TIkSOorKx0d3VkjcGQDGVnZyM2NhYDBw7ErFmzEB0djTvuuAPXrl1zd9UUpbq6Gg8++CDatGmDWbNmoXv37njzzTfdXS3Fevfdd3HzzTfj7rvvdndVFKWyshLPPPMMwsPDkZKSgsGDB+PGG2/EL7/84u6qydamTZvQtm1b3HnnnRg4cCD69OmDc+fOubtaivL5558jJiYGt9xyCyZPnoyIiAh8+OGH7q6WbDEYkqE5c+agffv2OHv2LA4ePIijR49iy5YteP/9991dNUV56KGHkJ6ejiNHjuDgwYM4ceIEgoOD3V0tRTp48CDee+89TJ8+3d1VUZzS0lKEhobi1KlT+PXXX3H69GkMGzYMEyZMQFVVlburJzvFxcW46667MHfuXJw4cQL5+fnw9/fHvffe6+6qKUpBQQG2bduGEydOIDc3F0uXLsXcuXOxZ88ed1dNlhgMydDFixfRv39/NGlSt1tKu3btEBkZiYsXL7q5Zspx8uRJfPbZZ1i4cCGio6MBAE2aNMGsWbPcXDPlKS8vR0pKCj766CO0atXK3dVRnNatW+Ppp59G8+bNAQDe3t6YNWsWLl26xCEgM7799ltUVVXhqaeeAlC3Rs5TTz2FjIwMnDlzxs21U47nnnsO7du3N76eMmUKgoKC2CNpAfcmk6FXXnkFf//739GhQwe0b98emzZtQkVFBebOnevuqinGtm3bIAgCxo0bh1OnTqGkpASdOnVCs2bN3F01xXnwwQcxZswYjBkzBps3b3Z3dTzC/v374e3tjQ4dOri7KrJz+PBhdOnSxRg8AkB8fDwAIDMzE1FRUe6qmqIdO3YMpaWl6NSpk7urIksMhlzg2LFjKCwstFqmX79+8PPzAwCMGDEC8fHxSE1NRWRkJE6cOIFnn30W7dq1c0V1ZWvXrl1W3w8KCkKPHj0AAPn5+QgODsZjjz2GLVu2ICgoCCdOnMD8+fPx9NNPu6K6snTx4sVGeyM6deqE8PBwAHV5B1lZWThw4IArqqcYOTk5KCkpsVpm0KBB0Gg0DY4fP34czz//PB577DFVbprZmOLiYoSGhpocM7wuLi52R5UUr7q6GjNmzEDv3r0xduxYd1dHlhgMucDq1auxfv16q2W+/vprtG3bFoIgICkpyZgz5OvrizNnziA+Ph41NTV44YUXXFRr+Xn22Wetvh8XF4cPPvgAAODj44PLly+jefPmOHXqFDQaDX744QdMmDABCQkJGDp0qCuqLDuHDx/GK6+8YrXMU089hQkTJuDixYt45JFH8OabbxqDoXPnzqGiogK7du1C9+7dVZuD9emnn2L//v1Wy2zduhW+vr4mx/Lz8zF69GgMHDgQ//jHP5xZRcXy8fFpMPPJMHmkfntS42pqajBlyhTk5+djx44dxvQLMsXtOGTm3LlziIyMxPr165GUlGQ8Pnv2bBw+fBj79u1zY+2UY9WqVUhOTkZ2draxtwgA2rZtizlz5qg6qBTr9OnTmDp1qsmxkydP4vLly+jduzfefvttJCQkuKl2ypOfn4/hw4ejU6dO+O677/hgt+CFF17AihUrcPLkSeOx3NxcdOvWDbt27cKgQYPcVzmFqampQUpKCvbv34+ffvqJw7JWMIFaZkJCQuDl5YWzZ8+aHD9z5gyTV20wfPhw+Pj4mEzHvXr1KkpKStiOIrVr1w67du0y+d+UKVMQHR2NXbt2MRCyQUFBAYYPH46YmBisWbOGgZAViYmJOHXqFLKzs43H1q5di+DgYPTv39+NNVOW2tpaTJ06lYGQSOwvk5mmTZvigQcewHPPPYfa2lp07NgRmzZtQnp6On788Ud3V08xWrVqhWeffRZz5szBa6+9huDgYPzrX/9Cq1atkJyc7O7qkYqUlpZi5MiRqK2txZNPPmkyvKbmoUZLhg4diqSkJCQnJ+OVV15BYWEhXnnlFbz55psMIm1w7733Yt26dVi2bBnOnj1r/IIdGRnJwMgMDpPJUG1tLZYvX46NGzeiqKgI7du3x8yZM3HTTTe5u2qKs3z5cnz77bfQ6XTo3bs3nnjiCbRs2dLd1VKs999/H3v37sXy5cvdXRXF+PPPPzFjxgyz73Go0bxr167hnXfewfbt2+Hv74+UlBR+ibHRrbfeajbJPyUlBQ899JDrKyRzDIaIiIhI1ZgzRERERKrGYIiIiIhUjcEQERERqRqDISIiIlI1BkNERESkagyGiIiISNUYDBEREZGqMRgiIod99913OH36tLurISslJSVYuXIldDodAODSpUtYuXIl9Hq9m2tGRPUxGCKiRv3444/49ddfLb7/wAMP4JdffnFhjeTv5MmTSElJwdWrVwEAv/32G1JSUlBdXS3q/AsXLmDlypXgurhEzsdgiIga9dhjj+Hrr7+2+P6kSZPQvn17F9ZIeQz74nl7e4sqn5ubi5SUFNTW1jq5ZkTEjVqJyGFJSUlo27at8fXq1asxaNAgaDQaZGVlITg4GPHx8fDyMv3+VVVVhd27d+PKlSvo3r17gw0k09PTUVJSAo1Gg7Zt26J3795o1qyZSRnDZ1VXVyMzMxMxMTHo0aNHgzoayun1emRnZyMkJAQJCQnQaDSirnXhwgXs27cP/v7+6NOnj9kNVg8fPoxz587hxhtvbPBeaGgobr/99gbB0KVLl7Bv3z40a9YMN910E7RaLUpLS7F161YAwKpVq+Dl5YWYmBjjru3Hjx/Hr7/+iuDgYAwcONBkA9Pz589jx44dmDx5Mvbv34/Tp09j5MiRaNGiRYM6EdH/CEREjYiJiRGef/55i++HhoYKX331lfG1VqsVkpKShJiYGOHWW28VwsLChNGjRwt6vd5YZs+ePUKbNm2Efv36CePGjRNCQkKEJ554wuS6f//734Xk5GRh8uTJQq9evYQ2bdoI+/btMymj1WqFMWPGCB06dBAmTpwofP3112brqNVqhdGjRwuRkZFCUlKSEBISIowbN07Q6XSNXuutt94SgoKChFGjRgnDhg0TQkJChB9++MHk+vfdd5/QvHlzYfTo0UJ0dLSQlJQkABAuX74sCIIg7Ny5UwAgXLt2zXjOP//5T6Fp06bCwIEDhWHDhglxcXHCyZMnhTNnzgjDhw8XAAh33XWXkJycLHz00UeCIAjC448/LjRr1kwYNWqU0LlzZ6FDhw7Cb7/9Zrzmhg0bBC8vLyEpKUno27evkJycLBw/ftzifzsiEgQGQ0TUKHuCoaFDhwqVlZWCIAjC6dOnBV9fX2HTpk2CIAjC1atXhfDwcGHZsmXGc86ePSuEhoYK69ats/g5qampwk033WRyTKvVCv379xeuXr1q9R60Wq3QpUsXoaSkxFinFi1aCB9//LHVa23ZskVo0aKF8McffxiPrVy5UggJCRFKS0sFQRCEH3/8UfD19RVycnIEQRCEa9euCTfddJPVYGjjxo2CRqMRNmzYYLzub7/9Jhw9elQQBEHYtm2bAMAkWNu4caPg7e0tHDhwQBAEQdDpdEJSUpIwfPhwY5kNGzYIAKz+9yIiUxwmIyKnuO+++6DVagEAUVFR6NixI37//XckJiYiPT0dRUVFCAwMxDfffAOh7osZoqOjsW3bNtx6663G65w8eRK//fYbSktL4e/vj4MHD6K2ttZkuGnmzJnw9/dvtE4PPPAAgoKCjHVKSUnBqlWr8MADD1i81v/93/+ha9euyMrKQmZmJgRBgF6vR0lJCbKysnDzzTdj1apVGDt2LLp16wYA8PPzw6OPPoo9e/ZYrMu///1vjBw5EmPGjDEe69Kli9X6f/XVVxg9ejT69u0LAGjSpAmeeeYZDBs2DBcuXEBYWJix7KOPPtpoexBRHQZDROQUISEhJq+1Wi0qKysB1AU4vr6+WLNmjUmZmJgYdOzYEQAgCAJmzpyJVatWIT4+HqGhoSgrK4NOp0NJSQlCQ0ON50VERIiqU/2cpOjoaGRkZJgcq3+tkydPoqioCN98843J8cmTJxtzdU6fPo2ePXs2uLY1p0+fRlxcnKh6G5w6dcoYcBnExMQY3zMEQ97e3iaBERFZx2CIiFwuMDAQOp0O//nPf+Dj42O2zM6dO7FixQocP34ckZGRAIDNmzdj48aNDaabX58Ebc3ly5cbvG7ZsqXVawUGBqJly5ZYuXKlxeuGhoaavbY1wcHBKCoqElNto5YtW6K4uNjkmOF1/fsgIvE4tZ6IXO6WW25BbW0tPv/8c5PjOp0OFy9eBFA3KyogIABt2rQxvl+/d8ZW33//vfHfer0e33//PQYNGmT1nDFjxiA9PR0nT540OX7+/HnjAoqDBw/Gxo0bjT1fABr0etU3atQobNiwwXi/AFBbW2sMopo3bw4AJtc0fI5h7SKgbgZcVFQU2rVrZ/XziMgy9gwRkSi5ubkNekfi4uIQGxtr87U6dOiAN998Ew8//DCysrLQp08fnD59Gt988w0WL16MESNGYOjQoaitrcW0adOQmJiInTt3mgQz9ti7dy+mTZuGoUOHYs2aNbh06RKeeuopq+fMmjULa9euxYABA/DII4+gVatWyMzMREZGBo4cOQIvLy/MmjULS5YswciRI3HPPfcgMzOz0brOmjULq1evRkJCAubOnQtfX1989dVXeO+99xAfH48bbrgBAQEBSE1NxcCBA9GpUyfMnj0by5Ytw7Bhw3D//ffj2LFjWLx4Mb766ivR6xcRUUPsGSKiRt16663w9fXF999/b/K/33//HUDDRRfvuusuk3WHgLoelq5duxpfz5s3D7t27UKzZs2wc+dO+Pj4YN26dRgxYgQAoHXr1ti7dy/atm2L7du3o0uXLti6dSuSk5ONidmWPsuSjz/+GIMGDcKBAwcQFxeH/fv3m+TWmLuWr68v0tPT8e677+L06dM4cOAAevTogczMTOMQX9OmTbF7926MHDkS+/btQ6dOnfDTTz8hOTnZmFdUf9FFrVaLLVu24MUXX8TRo0dx6tQpLFmyBPHx8QDqhue2bNkCb29vrFu3DocPH4ZWq8Uvv/yClJQU7N69GzU1NdixYwfuuOMOY30jIiKQnJwsqj2IqI5GqD/4TkTkgfz8/PDNN9+YzFQjIgLYM0REREQqx2CIiFTBluE0IlIXDpMRERGRqrFniIiIiFSNwRARERGpGoMhIiIiUjUGQ0RERKRqDIaIiIhI1RgMERERkaoxGCIiIiJVYzBEREREqsZgiIiIiFTt/wHVq4SfSs/EUwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "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": { "iopub.execute_input": "2026-07-28T14:15:37.491378Z", "iopub.status.busy": "2026-07-28T14:15:37.490724Z", "iopub.status.idle": "2026-07-28T14:15:37.498085Z", "shell.execute_reply": "2026-07-28T14:15:37.497640Z" }, "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": { "iopub.execute_input": "2026-07-28T14:15:37.500493Z", "iopub.status.busy": "2026-07-28T14:15:37.499871Z", "iopub.status.idle": "2026-07-28T14:15:37.549387Z", "shell.execute_reply": "2026-07-28T14:15:37.548536Z" } }, "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: Tue, 28 Jul 2026 Deviance: 7.2134\n", "Time: 14:15:37 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": { "iopub.execute_input": "2026-07-28T14:15:37.553018Z", "iopub.status.busy": "2026-07-28T14:15:37.552848Z", "iopub.status.idle": "2026-07-28T14:15:37.762970Z", "shell.execute_reply": "2026-07-28T14:15:37.762588Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/hostedtoolcache/Python/3.14.6/x64/lib/python3.14/site-packages/statsmodels/base/model.py:1315: FutureWarning: linear keyword is deprecated, use which=\"linear\"\n", " predict_results = self.model.predict(self.params, exog, *args, **kwargs)\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'Residual')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAGxCAYAAABx6/zIAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQ+NJREFUeJzt3Xl8lPW5///3JCQTtgQShIRF2QQNyL6IqCAUCCBSrBZBy696qkes2mpPXYr9Alartlbt0SMudWnlIAouyEGCiFqwIhD2sIkIApqwBJNAMGHI3L8/0okMyWQmk5m5l3k9Hw8fknvuTK75zHJf81muj8swDEMAAAAOkGB2AAAAAJFCYgMAAByDxAYAADgGiQ0AAHAMEhsAAOAYJDYAAMAxSGwAAIBjkNgAAADHaGR2ALHm9Xr17bffqnnz5nK5XGaHAwAAQmAYho4fP662bdsqISFwv0zcJTbffvutOnToYHYYAAAgDAcOHFD79u0D3h53iU3z5s0lVTVMamqqydH8wOPx6IMPPtDo0aOVlJRkdji2QbuFj7YLH20XPtoufPHedqWlperQoUP1dTyQuEtsfMNPqamplktsmjRpotTU1Lh8wYaLdgsfbRc+2i58tF34aLsqwaaRMHkYAAA4BokNAABwDBIbAADgGCQ2AADAMUhsAACAY5DYAAAAxyCxAQAAjkFiAwAAHIPEBgAAOIbtKg+fPn1aO3bs0IkTJ9StWzdlZGSYHRKABqr0Glq795gOHy9X6+YpGtQpXYkJbFILoP5sldi88cYbmjFjhtLS0tSoUSNt2bJF06dP1xNPPGF2aADClJtfoNmLt6ugpLz6WFZaimZOyFZOzywTIwNgR7ZKbCoqKvTZZ5+pdevWkqTPP/9cQ4YM0ejRo5WTk2NydADqKze/QNPnbpBx1vHCknJNn7tBc27oR3IDoF5sNcdm2rRp1UmNJPXt21eNGjXSkSNHTIwKQDgqvYZmL95eI6mRVH1s9uLtqvTWdgYA1M5WPTaSdPToUeXl5am0tFSvvvqqLrnkEv3kJz8JeH5FRYUqKiqqfy4tLZVUtUuqx+OJeryh8sVipZjsgHYLn9ltt2bvMb/hp7MZkgpKyrX6y8Ma3Ck9doGFwOy2szPaLnzx3nahPm6XYRi2+jq0fv16zZgxQ0ePHtW+ffs0a9Ys3X777QHPnzVrlmbPnl3j+Lx589SkSZNohgqgDuuPuvSP3YlBz5t2fqX6t7LVxxSAKDh58qSmTp2qkpISpaamBjzPdonNmdatW6fLL79czz//vKZNm1brObX12HTo0EFHjx6ts2FizePxaPny5Ro1apSSkpLMDsc2aLfwmd12a/Ye0w0v5wU9b+5NAyzZY8PrLjy0Xfjive1KS0vVqlWroImN7YaizjRw4EANHDhQK1asCJjYuN1uud3uGseTkpIs+cKwalxWR7uFz6y2G9K1tbLSUlRYUl7rPBuXpMy0FA3p2tqyS7953YWPtgtfvLZdqI/ZNpOHKysrdfz4cb9j5eXl+uqrr9SmTRuTogIQrsQEl2ZOyJZUlcScyffzzAnZlk1qAFiTbXpsPB6PBg8erKlTpyo7O1vFxcV66aWXZBiG7rzzTrPDAxCGnJ5ZmnNDvxp1bDKpYwMgTLZJbFJSUrRq1So9//zzev3119W0aVP99Kc/1U033aTmzZubHR6AMOX0zNKo7EwqDwOICNskNpKUkZGh3/3ud2aHASDCEhNcGtKF7VEANJxt5tgAAAAEQ2IDAAAcg8QGAAA4BokNAABwDBIbAADgGCQ2AADAMUhsAACAY5DYAAAAxyCxAQAAjkFiAwAAHIPEBgAAOAaJDQAAcAwSGwAA4BgkNgAAwDFIbAAAgGOQ2AAAAMcgsQEAAI5BYgMAAByDxAYAADgGiQ0AAHAMEhsAAOAYJDYAAMAxSGwAAIBjkNgAAADHILEBAACO0cjsAAAgHlR6Da3de0yHj5erdfMUDeqUrsQEl9lhAY5DYgMAUZabX6DZi7eroKS8+lhWWopmTshWTs8sEyMDnIehKACIotz8Ak2fu8EvqZGkwpJyTZ+7Qbn5BSZFBjgTiQ0AREml19Dsxdtl1HKb79jsxdtV6a3tDADhILEBgChZu/dYjZ6aMxmSCkrKtXbvsdgFBTgciQ0ARMnh44GTmnDOAxAciQ0AREnr5ikRPQ9AcCQ2ABAlgzqlKystRYEWdbtUtTpqUKf0WIYFOBqJDQBESWKCSzMnZEtSjeTG9/PMCdnUswEiiMQGAKIop2eW5tzQT5lp/sNNmWkpmnNDP+rYABFGgT4AiLKcnlkalZ1J5WEgBkhsACAGEhNcGtIlw+wwAMdjKAoAADgGiQ0AAHAMEhsAAOAYJDYAAMAxSGwAAIBjkNgAAADHILEBAACOQWIDAAAcgwJ9AOBglV6DiseIKyQ2AOBQufkFmr14uwpKyquPZaWlaOaEbPaogmMxFAUADpSbX6Dpczf4JTWSVFhSrulzNyg3v8CkyIDoIrEBAIep9BqavXi7jFpu8x2bvXi7Kr21nQHYG4kNADjM2r3HavTUnMmQVFBSrrV7j8UuKCBGSGwAwGEOHw+c1IRzHmAnJDYA4DCtm6dE9DzATkhsAMBhBnVKV1ZaigIt6napanXUoE7psQwLiAkSGwBwmMQEl2ZOyJakGsmN7+eZE7KpZwNHIrEBAAfK6ZmlOTf0U2aa/3BTZlqK5tzQjzo2cCwK9AGAQ+X0zNKo7EwqDyOukNgAgIMlJrg0pEuG2WEAMcNQFAAAcAwSGwAA4BgMRQFx4uxdnvu2b252SAAQcSQ2QByobZfnzFS3xmW6NM7EuAAg0hiKAhwu0C7Ph0or9PIXCVq27ZBJkQFA5JHYAA4Wyi7PDy/dyS7PABzDVkNRJ0+e1GuvvabPPvtMjRo10qWXXqpp06YpMTHR7NAASwq2y7PkUkFJhdbuPcaSYACOYJseG6/Xqx49emjTpk0aOXKkBg4cqNmzZ2vChAnyer1mhwdYErs8A4g3tumxcblcWrNmjVq3bl19rF+/fho8eLDy8vI0aNAgE6MDrIldngHEG9v02LhcLr+kRpIyMzMlSaWlpWaEBFhesF2eJUNZaW52eQbgGLbpsanNE088ofT0dA0ePDjgORUVFaqoqKj+2ZcEeTweeTyeqMcYKl8sVorJDmi34GaM7a475m+WS/KbROz7+b7R58tbeVreSnPisyNed+Gj7cIX720X6uN2GYZhy+UQr732mm688UYtXLhQP/7xjwOeN2vWLM2ePbvG8Xnz5qlJkyZRjBCwjs1FLr29L0HFp37ou2mRbOjqjl71zrDlRwCAOHPy5ElNnTpVJSUlSk1NDXieLRObBQsW6Prrr9fzzz+vG2+8sc5za+ux6dChg44ePVpnw8Sax+PR8uXLNWrUKCUlJZkdjm3QbqGr9BrK+/o7HT5eodbN3erdtpk+WvEhbRcGXnfho+3CF+9tV1paqlatWgVNbGw3FPXWW2/phhtu0HPPPRc0qZEkt9stt9td43hSUpIlXxhWjcvqaLfgkiRd2q1N9c++bl3aLny0Xfhou/DFa9uF+phtM3lYkt555x1NnTpVc+bM0U033WR2OAAAwGJs02NTWlqqyZMnKy0tTW+//bbefvvt6tvuvPNOjR492sToAACAFdgmsWncuLFfMnOmCy+8MMbRAAAAK7JNYpOUlKQrr7zS7DAAAICF2WqODQAAQF1IbAAAgGOQ2AAAAMcgsQEAAI5BYgMAAByDxAYAADgGiQ0AAHAMEhsAAOAYJDYAAMAxSGwAAIBjkNgAAADHILEBAACOQWIDAAAcg8QGAAA4BokNAABwDBIbAADgGCQ2AADAMRqZHQAAa6v0Glq795gOHy9X6+YpGtQpXYkJLrPDwr/x/AD+SGwABJSbX6DZi7eroKS8+lhWWopmTshWTs8sEyODxPMD1IahKAA1VHoN/fXD3bp17ga/i6YkFZaUa/rcDcrNLzApOkhVSc10nh+gBhIbAH5y8ws09NEVevLDL2q93fj3/2cv3q5Kr1HrOYiuSq+h2Yu3q7bW5/lBvCOxAVDN1wtQWFpR53mGpIKScq3deyw2gcHP2r3HavTUnInnB/GMxAaApLp7AQI5fDzwxRXRE2q78/wgHpHYAJAUvBegNq2bp0QpGtQl1Hbn+UE8IrEBIKl+3+5dqlp9M6hTevQCQkCDOqUrKy1FgRZ18/wgnpHYAJBU/2/3MydkUy/FJIkJLs2ckC1JNZIb3888P4hXJDYAJAXvBfDJSkvRnBv6USfFZDk9szTnhn7KTPNPSDN5fhDnKNAHQNIPvQDT526QS6p1EvFdPzpft484n54Ai8jpmaVR2ZlUHgbOQGIDoJqvF4BqtvaRmODSkC4ZZocBWAaJDQA/9AIAsDMSGwA10AsAwK6YPAwAAByDHhsAgO1Ueg2GS1ErEhsAgK3k5hcwwR0BMRQFALAN30atZ2//UVhSrulzNyg3v8CkyGAVJDYAAFuoa6NW37HZi7er0lufrVzhNCQ2AABbCLZRqyGpoKRca/cei11QsBwSGwCALYS6UWt9NnSF85DYAABsIdSNWuu7oSuchcQGAGALwTZqdalqddSgTumxDAsWQ2IDALAF30atkmokN76fZ07Ipp5NnCOxAYAQVXoNrd5TpMVbCrS7xMXqGxP4NmrNTPMfbspMS9GcG/pRxwYU6AOAUNQsCpeohX9ZqVlX9eBiGmNs1Iq6kNgAQBC+onBn988cKq3Q9Lkb6CkwARu1IhCGogCgDhSFA+yFxAYA6kBROMBeSGwAoA4UhQPshcQGAOpAUTjAXkhsAKAOFIUD7IXEBgDqQFE4wF5IbAAgiMBF4dws9QYshjo2ABCCM4vCFRSX6attm3T75MuV4k42OzQAZyCxAYAQ+YrCeTypev/gRoafAAtiKAoAADgGiQ0AAHAMhqIAxESl12DTQgBRF3JiU1xcHPKdtmjRIoxQADhVzZ2xq2q/zJyQzYoiABEVcmLTsmXLkO/UMNgMDkCVQDtjF5aUszM2gIgLObFZt25dNOMA4EDBdsZ2qWpn7FHZmQxLAYiIkBObAQMGRDMOABFmhTkt9dkZe0iXjNgFBsCxmDwMOJBV5rSwMzaAWAs7sVmwYIEWLFig/fv36/Tp03635eXlNTgwAOGx0pwWdsYGEGth1bH561//ql/+8pfq3Lmz1qxZoyuvvFKtW7fW+vXrYzJkVVFRoYMHD6qioiLqfwuwskqvodV7irRo0zdavadIp05765zTIlXNaan0xmaCPztjA4i1sHpsnn32WS1YsEDDhg3TY489plmzZkmSHn/8ca1cuTKS8fnZu3ev/ud//kevvfaaDh8+rI8//ljDhw+P2t8DrKy24ab0psk6VnYq4O+cOadlwLmpUY/RtzP29Lkb5JL8Eq5o7IxthXlFAMwVVmLz1VdfaciQIZIkt9utEydOqFmzZvrFL36hBx98MKIBnumdd95RmzZttGzZMvXt2zdqfwewukDDTXUlNWeqmtMS/cRG+mFn7LOTsMwIz/mxyrwiAOYKK7E5ffq0kpOrdrQ999xztX79eg0bNkwFBQVKTEyMaIBnuvvuuyVJBw8ejNrfAKyuriXUoYr1nJYzd8aORm+KleYVATBXg1dFTZs2Tdddd53GjBmjjz76SBMnToxEXAACCLaEui4uVfWUDOqULm/l6aDnR5JvZ+xIo1YOgDOFldgcOHCg+t8zZszQOeeco9WrV+u2227Tr3/960jFFhEVFRV+k4xLS0slSR6PRx6Px6ywavDFYqWY7CAe262guCys3/Nd0meM7S5v5WnHtN2aEGvlrP7ysAZHaJKyU9rODLRd+OK97UJ93C7DhvsfHDx4UB06dAhp8vCsWbM0e/bsGsfnzZunJk2aRClCIHp2l7j0zPbgQ75NGxkqO/1DD0WLZENXd/Sqd4bt3vJ1Wn/UpX/sDt4e086vVP9WznrsQDw5efKkpk6dqpKSEqWmBp4jGFaPTWFhYZ23Z2ZmhnO3UXH//fdXz82RqnpsOnTooNGjR9fZMLHm8Xi0fPlyjRo1SklJSWaHYxvx2G6VXkML/7JSh0orah1+qRpucuvDX1+mjQeKdfh4hVo3d2vAeS39hmKc0nYZe4/pH7uD184afdngiPbYOKHtzEDbhS/e28434hJMWIlNVlbdk/Cs1AnkdrvldrtrHE9KSrLkC8OqcVldPLVbkqRZV/UIsoS6h5o2duvSbm2C35/N225I19bKSktRYUl5HYleioZ0bR3xOTZ2bzsz0Xbhi9e2C/Uxh1Wgb+vWrX7/bd68WQsXLlS3bt307LPPhnOXITl58qQOHjxY3WN05MgRHTx4MOQsDnAK3xLqzDT/1U2ZaSlxtwLIVytHUo1CgNGolQPA2sLqsenZs2eNY7169dJ5552nO+64Q9OnT29wYLVZsmSJ7rrrLklSu3btqv999913+w03AfEg2kuo7SRWtXIAWF9EN8G88MILlZ+fH8m79HPttdfq2muvjdr9A3YTrSXUdkSiB0CKYGLj8Xj01FNPqV27dpG6SwCoFxI9AGElNi1atKhx7MSJE2rcuLHmzZvX0JgAAADCElZi87e//a3GsZYtW6pv375KT2eXXiAUbNgIAJEXVmJzzTXXRDoOIK6wYSMAREfIiU2wonxnslKBPsBqnLZhIz1PAKwk5MQmWFG+M1mpQB/Mx4XvB07bsJGeJwBWE3Jis2PHjup/f/jhh3r44Yc1Y8YMDRw4UJK0bt06Pfzww3rggQciHyVsiwufv2A7c/s2bFy795jlV/c4recJgDOEnNhccMEF1f++/vrrtXDhQg0dOrT62ODBg9WnTx/96le/0i9/+cvIRglb4sJX0+HjgZOacM4zi9N6ngA4R1hbKuzYsUM9evSocbxHjx7auXNng4OC/QW78ElVF75Kb3wNW7ZunhL8pHqcZ5b69DwB4aj0Glq9p0iLNn2j1XuK4u6zAuELa1XUueeeq6efflq///3v/Y4/88wzOu+88yISGOzNSUMukTSoU3pIGzYOitAu1NHilJ4nWBND2GiIsBKbp556Sj/+8Y+1YMECDRw4UIZhKC8vT7t379aiRYsiHSNsiAtf7XwbNta9M7f1N2x0Ss8Tqlhpgj9D2GiosBKbnJwcffnll5ozZ462b98uSZowYYJuu+02tlSAJC58dXHCho1O6XmCtXpHmLuFSAh7r6j27dvr4YcfjmQscBAufHWz+4aNTul5indW6x1hCBuREPLk4eLiYhUXF/v9O9B/gO/CJ/1wofPhwlfFt2HjxD7tNKRLhu3awtfzlJnm3+uWmZbCcIENWHGCP0PYiISQe2xatmwpqar4nu/fgVCgD5IzhlxQN7v3PMUzK/aOMISNSAg5sVm3bl2t/wbqwoXP+Xw9T7AXK/aOMISNSAg5sRkwYECt//Y5efKkGjduLJeLCxb8ceEDrMeKvSPM3UIkhFWgb9euXbr//vurf/7Nb36jZs2aqUOHDtqyZUvEggOAUFDMrf58vSOBUgSXqlZHxbp3hLlbaKiwVkXddddd+vWvfy1J2r17t5577jm98847+uijj3TPPfcoNzc3kjECQEBWWq5sJ1buHWEIGw0RVmLzr3/9SwsWLJAkLVu2TBMnTtTEiRN1+eWXq3PnzhENEAACsdpyZbux8gR/hrARrrASm+TkZB06dEidO3fW0qVLNWnSJEnS6dOnlZycHNEAAaA2FHOLDHpH4DRhJTbjx4/XNddcowEDBmjlypV6+eWXJUkrVqzQyJEjIxogANTGisuV7YreEThJWJOHn3nmGY0ZM0YnTpzQu+++qzZt2kiSPv74Y82aNSuS8QFAray4XBmA+cLqsWnWrJkeeeSRGseff/75BgcEAKGw4nJlAOYLq8dGqppPs2bNGr3++uvVx0pKSiISFAAEY9XlygDMFVZic/DgQQ0YMECXXXaZpk6dWn182rRpWrJkScSCA4BA2I8sfNT9cbZ4f37DrmPTv39/ff7552rcuHH18XvuuUf33nuvxo8fH7EAASAQKy9Xtqr61v2p9BqsmLIR6jqFmdh88skn2r59u1JS/Meue/Xqpby8vIgEBgChYLly6Opb94eLpL1Q16lKWENR33//vRo1qsqJztwbqrCwUE2aNIlMZAAQIt9y5Yl92mlIlwySmloEq/sjVdX98Q1b+C6SZy+p910kc/MLohsw6qW+z6+ThZXYXH755XrhhRck/ZDYfP/997rnnns0YsSIyEUHAIiI+tT94SJpP/V5fp0urKGoxx9/XMOGDdPSpUtlGIamTJmilStXyuPx6F//+lekYwQAhOHM+TG7D50I6XcOHy+n+KENUdfpB2ElNtnZ2dqyZYuee+45NW3aVMXFxfr5z3+u22+/XVlZzh+/AwCrq21+TChaN0/hImlDodZr2n3ouFbvKXL0PLSwEhtJysrK0uzZs2scz8/PV8+ePRsUFAAgfIEmkdbFparVZIM6pYc8XEHxQ+vw1XUqLCmv83l/5uM9eubjPY6eBB7WHJvTp0/r1KlTfscOHDigG2+8Ub17945IYACA+qtrfkwgZ9f9ofih/dRV16k2Tp4EXq/E5siRI7rqqqvUtGlTNW7cWNdcc43Kysr03//93+revbtWr16tBQsWRCtWwPHivbAWGi7Y/JjaZKal+C0FpvihPfnqOmWmBe9Jc/Ik8HoNRd13333atm2bHn74YUnSs88+qzFjxmjr1q16/PHHdcstt1QvA0d4Kr2G8vYUUY8jDsW6Zkil19DuEpcWbylQVoumvNYcItR5L7df0UXnt2ke8HOG4of2dGZdp399eVTPfPxlwHOdOgm8XlnIsmXL9H//93/q06ePJGnEiBHq37+/cnNzNWbMmGjEF1c2F7n0yF9WqrC0ovqYk8dB8YNYF9bKzS/QrPe2qbA0Udq+VRKvNacIdd7L0K7nBL2YUfzQnnx1neJ1Eni9hqIKCgrUq1ev6p9982moXdNwy7Yd0stfJPglNZKzx0GdIBJDR7GuGeJLonitOVOk58dQ/NC+Qk1ynTYJvF49Nl6vVwkJP+RCiYmJkqSkpKTIRhVnKr2GHnp/Z623Gar6IJq9eLtGZWfyoRJhDdkHJ1JDR7GsGRIsieK1Zn+++THT526QS/J7rpkfE1+CrZQ6cyWck9R7QkzXrl2DHvvyy8Bjeqhp7d5j//72XPsHjVPHQc3WkMQkkkNHsewupvBafGB+DKT4TXLrldjMmDEjWnHEtXgdBzVTQxKTSPd6xLK7mNda/GB+DKT4THLrldg89NBD0YojrsXrOKhZGpqYRLrXI5bdxbzW4otvfgziW7wluWEV6ENkDeqUrsxUtxSgpBbFsCKroZvFRbrXI5Y1Qyi8BsSneJoETmJjAYkJLj0w7gJJFMOKhYYmJtHo9QhUWOvswmkNReE1AE5HNT2LGNOjjW7q5tX7hU38luE6eRzULA1NTKI1dBSr7mJfElVVx4bXGgBnIbGxkN4Zhu65/nJtPHg8LsZBzdLQxCSaKw1iNScip2eWhp+foWfeyFXnHn2oPGxRDSlHAMQrEhuLYbJf9EUiMXHCSoPEBJfOTzM0rlcWtagsKNZbbABOQWKDuBSJxCTeVhogdmK9xQbgJCQ2iFuRSEzoYUOkUR06Nhjmcy4SG8Q1EhNYDdWho49hPmdjuTdgA5HYbBP2QHXo6PIN852dPBaUlOtWNoF1BHpsAIvj22V8oTp09NQ1zOdz39tbGeazOXpsAAsL9O3SN4mUb5fOQ3Xo6Ak2zCdJxSc9euYjNnK2MxIbwKKCTSKVqiaRMizlLFSHjp5Qh+9e+Wwv7ysbI7EBLKqhe1rBvmK1xUa8CXX4rvikh/eVjTHHBrAoJpHGN+okRd6gTulq0ThJxd97gp7L+8q+SGwAi2ISKShHEFmJCS7dOLSjnvxwd9BzeV/ZF0NRgEUxiRSIvNtHnK8WTQJvIcL7yv5IbBAUNVTMwSRSWIHT3v+JCS49evVFtX5h4H3lDAxFoU7UUDGXEzbbDAfl7q3Bqe//ut5Xvx+frbTGyVq06RteezZFYoOA2IjPGuJtEmk0LqYkSvXn9Pd/be+r78oq9Iclzkvk4g2JDWrFRnzWEi+TSKNxMXVqr0M0xcv7/8z3VW5+gX45b6NjE7l4whwb1IoaKoi1aBQkpHJzeOLt/U8xTGexXWJTVFSkOXPmaNasWXrrrbfk9XrNDsmRqKGCWIv0xZSLVfji7f0fb4mc09kqsdm7d68uuugivfHGGyorK9NvfvMbTZgwgeQmCqihgliL9MWUi1X44u39H2+JnNPZKrG599571bFjR61YsUJ//vOf9dFHH+mDDz7QG2+8YXZojkMNFcRaqBfJVk3dIZ3HxSp88fb+j7dEzulsk9icPn1aixcv1g033KDExERJUufOnTVs2DC9/fbbJkfnPNRQQawFu5j6/GbB5pDmxnCxCl+8vf/jLZFzOtusitq/f7/Ky8vVtWtXv+Ndu3bV6tWrA/5eRUWFKioqqn8uLS2VJHk8Hnk8wfcLiRVfLFaKaWT3Vnr6ut566P2dKiz9oQ0z09yaMfYCjezeyvR4rdhudmHFtpsxtrvumL9ZLqnWuTGSdKi0auLv09f11pgebQLeV9/2zZWZ6tah0opa78ulqtdy3/bN690GVmy7SIvW+9+qbRfotec643Zv5Wl5K00I7t+s2naxEurjdhmGYYuZc1u3blWvXr20evVqXXzxxdXH7733Xr311lv68ssva/29WbNmafbs2TWOz5s3T02aNIlavE7iNaQ9pS6VeqTUJKlLqiGHfFFDlDTkNbO5yKW39iaoxFPXLxhqkSzN7FdZ5/1uLnLp5S98HdNnnlj1sXdTN696Z9jiI9A08fT+31zk0tv7ElR86ocH2CLZ0NUdeZ1YwcmTJzV16lSVlJQoNTU14Hm26bFp1qyZJKmkpMTveHFxcfVttbn//vt19913V/9cWlqqDh06aPTo0XU2TKx5PB4tX75co0aNUlJS4H1M4I92C1+02m7ZtkN65Oxv+aluPTDugjp7WHzGSRq2p0j/36vr6zjLpeJT0jnZF2twHcMD4yT123aoRq9DVlqKZowNLZ7a8LoLn5Xbbpyke7yG8r7+ToePV6h1c7cGnNfSMkNuVm67WPCNuARjm8Tm3HPPVdOmTbVr1y6NGTOm+viuXbt04YUXBvw9t9stt7vmZMOkpCRLvjCsGpfV0W7hi2Tb5eYX6I75m2sM/RwqrdAd8zeHXOSsuDy0/v6ik6eDxn5ln/Ya26tdVCoP87oLn1XbLknSpd3CS3hjxaptF+0K36E+ZtskNomJiZo0aZL+8Y9/6NZbb1VycrK2b9+uTz/9VG+++abZ4QFxL5LVaiM98TdeKjfbDVtdOIeVKnzbJrGRpEcffVSXXXaZhgwZov79++u9997TNddco0mTJpkdGhD36lM3JliS4VulUlhSXsfEX1ap2NmybYf08NJdfq+Z9KbJemhiT43rxdYFdmK1fcVss9xbktq1a6ctW7boN7/5jS688ELNnTtX8+fPl8tFhg+YLZJ1Y+JtuXG82Vzk0h3zN9dIhI+VndJt8zbo4SXbTYoM9WXFCt+26rGRqiYRT5061ewwANvydf8XFJfpqxKXKr2GIjFaH+nho5yeWZpzQ78a3duZDtzAMp6GZCq9ht7elxBwOb8kvbhqryRDM8b3iFVYCFMke2ojxXaJDYDw1RwHT9TCv6zUrKt6NDhRiMbwUU7PLI3KznT0Rd9KcxNiIe/r7/yWUwfy4qp96tuhpcb1ahuDqBAuK1b4ttVQFIDwBdrp+lBpRUR2uo7W8JFv4u/EPu00pEuG45KaeNt9/PDxiuAn/dsDi/LZpNTirFjhm8QGiAOxGgf3DR9lpvl/iGWmpcR8AqHVWXFuQiy0bh7aXl+SdKzMwyalFmfF7SgYikJci5e5DbEcB4+H4aNIsOLchFgYcF5LNWtk6MTp0F4PbFJqbb6e2ulzNwTcjiLWE/1JbBC34mluQ6zHwakbE1wsnxMrJfCJCS5d29mrV75IDOl8Nim1PqtN9CexQVyyWt2FaLPiOHi8i9VzYsUEvk+Gof8Yep5e+tfXdZ7Hjtr2YaWeWubYIO7E49wGK46Dx7tYPCdWnpx8X0533XxZx4C3u0StIruxykR/EhvEnfrMbXCKeCh4V+k1tHpPkRZt+kar9xRZPjGN9nNihwR+xvgeenZqX6U39a+klMVkczQAQ1GIO1asuxAJweZR5PTM0v9M7asHFuXrWJmn+nhmmlszJzS8jo2ZrDjcEopozk2wy+Tkcb3aakzPLEsMYcAZSGwQd5w43ySUC3tufoH+sGSHX1LTrJGh+3O6W/riH4zd50tFa26CnRJ4JpsjkkhsEHectsFiKBd2SbWec+K09Ks3tig5qZGlL/6BRHJHcTNF48LuxAQeCAVzbBB3nDTfJNR5FLPe2xZgbx5X9TlWn5NSm3icLxWqhk5OttucJcCHHhvEJavVXQhXqBf2ulhlrkU47DTcEmsNKZxm1zlLgERigzhmpboL4YrkBduOF3+GW+oWTgIfqzlLlV5DeXuKbPveg3WR2CCu2X3SYiQv2Ha8+DttvlQ01CeBj9Wcpc1FLj3yl5UqLP1hQ0x6hBApzLEBbCzUeRSZqW5HFudz0nypaAq1cFos5iwt23ZIL3+R4JfUSNYoGghnILEBTBCpiZmhXthnXdWj1nN8My/sfPFnR/HIifacpUqvoYfe31nrbVYpGgj7YygKiLFIT8wMdR5Fbee0SJYeurq37S/+TpgvZQXRnrO0du+xf/fU1P682HkiO6yDxAaIoWhNzAzlwn72ORlNGunI9s81pkebBj4qa7D7fCkriPacJVaxIRZIbIAYifbEzFAu7Gee4/F49P6Oev8ZOFhDloiHglVsiAXm2AAxQjE52EE05ywN6pSuzFS3VEe5SLtOZId10GMDxAjd8PEr2AalVhOtOUuJCS49MO4C3T5/U1R6hACJxAZB2O0DOZbq2zZ0w8cnu1bxjdacpTE92uimbl69X9jEb8m33ap+w7pIbBCQXT+QYyGctqGYXPyx+87j0dI7w9A911+ujQeP86UJEcccG9TK94F89pwQimiF3zYUk4svoW5QGq81W0ItGgjUF4kNauADObCGtg3F5IJzyq7STBYHzMFQFGqozwdyvNUNiUTbUEwuMLOGP6Mxl4zJ4oA5SGxQAx/IgUWqbSgmV5NZ81FqS6bSmybroYk9Na5X+H/PqpPFWRAApyOxQQ1W/UC2AtomOmK1q/TZAiVTx8pO6bZ5G3TzgY6aMb5HWPdtxcniLAhAPGCODWoIdcfoeFy9Q9tEhxnzUepKpnxeXLVPDy/ZHtb9W22yOAsCEC9IbFCD1T6QrYS2iQ4zhj+DJVM+L67aq/e3hHfRt8pkcRYEIJ4wFIVahbpjdDyibSLPjCG++iRJv1+UrzE9wxsGs8JkcRYEIJ6Q2CAgK3wgWxVtE1lmzEepT5JUVHZKr/5rr34+tFPUNiiNJhYEIJ6Q2KBOZn8gWxltEznR3lW6NoM6pSu9abKOlZ0K6fw/LNmhv326VzMnZGtk91YRiyMW9h09GdJ5THqHEzDHBoAlxHo+SmKCSw9N7Fmv3/FNtF227VBEY4mm3PwCPfXhF0HPS2+apMLSclsXRQQkemwAhCFatVDqO8TX0DjG9crSzQc66sVV+0I633e5/3+Lt2tG/XIiU4Sy8svnWJlHd72xSRJLwGFvJDaIKoqBOU+0a6GEOsQXqcJ6VXVqXHpx1d6Qf+dYmUf/b32i3J0O6co+7UP+vVgLdeXX2eJ9k07YG4kNooZiYM5jld2qgxXW+8+DnXT/uOyQ72/G+Gz17dBSDyzKD3nOTdlp6Y75m9WoUaJlX8/hTgaOZlFEINqYY4OooBiY81ilFkoowyvPr9yr97d8W6/7HdcrS+tm/Ei/H39hiL9RdbG3cv2XhkwGZpPO6HHKRq9WRY8NIi6UC+Dv3tmq7z1eZaYyPGUXVqmFEurwygOL8jWmZ1a9XluJCS79fGgn/e3TvQGXnp/J6vVfgi2jDwVLwCOLnuzoo8cGERfKhcc3UXHKi5/r0sc+ogfHBqxSCyXU+z9W5gmrt+HM6tKRjilc4X7Dr6tSdqhYAh459GTHBokNIq6+H/K8qe2hPtWBo9nVXp8LbbgJh2/peXrT5IjHVF+5+QW69LGPNOXFz/Wr+fX/MhBoGX1WWopaNEli37MYscpQbjxgKAoRV98PeSYq2kOo1YG/K6vQpY99FLWu9voU1mtIwpHTM0sjLmijix/5UMfKPLWeE+0duiM1WTvQMvrl2wtjWhQxnlllKDce0GODiAu2A3ZtmKhofaFsAHpV7yz9ct7GqHa1h1pY7+zehnB6kZIbJeiPky6SS7UN5VT9frQu/pH+hu9bRj+xTzsN6ZKhxASXZTbpjAdWGcqNB/TYIOLqKo8fDG9qa6trA9Dfj8/WH5YEvhBHslduXK8s/efBTnp+Ze21Z1zyTzgaMmEz0GNukSw9dHXvqF38Y/UNn33PYsOMjV7jFYkNoiLQxSAY3tTWF+hCGOuu9vvHZat3+xb/rj3zw1DR2QlLJIZzzn7MGU0a6cj2zzWmR5sGP45AYvkNn33Pos+MjV7jFYkNoubMi0Fhyff6w5Id+q7sFG9qB6jtQmhGV/u4Xm01pmdWwN6GYMM59elFOvMxezwevb8jYg+jVnb9hk+18dqZsdFrvCKxQVSdeTFonJzIm9rBzLoQ19XbYOcJm3b8hk+NlrrVNZRLG0UOiQ1ihje1s1nxQmznCZt2+4Zvle02rI45TdFHYoOY4k3tXFa8ENt1OMfHLl8G6jPkB+Y0RRuJDWKON7VzWe1CbMVepPqyw5eB+gz5DTg3NXaBIS6R2ACIqGhdiMOZlGrFXqRwWP3LQP2G/EhsEF0kNgAiLtIX4mjUobHacI6d2X3ID85CYgPA0qJRh8aKwzl2Vp8hP2/l6ViHhzjDlgoALCuS2wrUtqUAIiOU7TbsMOQHZyCxAWBZ9ZmUCnOx7xSsgqEoAJZl5zo08YghP1gBiQ0Ay2JSqv1YfQUXnI+hKACW5ZuUGuj7vktVq6OsXIcGQGyR2ACwLCalAqgvEhsAlsakVAD1wRwbAJbHpFQAoSKxAWAL8TwpNZztJIB4ZavEpqysTK+//rqee+457dy5U0uXLtVll11mdlhR+9Bx0oeZkx6LVZndxmb/faeqbTuJFo2TdOPQTrp9RFfaGDiLrRKbP/zhDzpy5IgeeOABTZo0SZWVlWaH1KA9bMy4XzM46bFYldltbPbfd6pA20kUf+/Rkx9+oVc+26tHr76INgbOYKvJw48++qheeuklDRgwwOxQJP3woXN2ZVTfHja5+QWWul8zOOmxWJXZbWz233equraT8Ck+6dGttDHgx1aJjZVEcg+bWNyvGZz0WKzK7DaO9N+v9BpavadIizZ9o9V7iuL6tRFsO4kz8T4CfmCroahwVFRUqKKiovrn0tJSSZLH45HH4wn7fteEuIfN6i8Pa3AIxcN8sXy+50hE79dMkW6j2vjarSHPpZ01pI0j0XaRfI6XbTukh97fqcLSH96vmaluPTDuAo3p0SbsGKMhFq+7guKy0M+1yWeCFNm2q/Qayvv6Ox0+XqHWzd0acF5LR885ivfPu1Aft6mJza233qq5c+fWec6WLVvUuXPnsP/GI488otmzZ9c4/sEHH6hJkyZh3+/6oy5JiUHP+2DVGhXtCP2b1Eer10flfs0QrTaqzfLlyxv0+3YViTZuSNtF6jneXOTSy1/4OpB/uDAVlpbr9vmbdFM3r3pnWO/1Hs3X3VclobWtjx0+E87U0LbbXOTS2/sSVHzqh9dLi2RDV3e05mslkuL18+7kyZMhnecyDMO0V0BFRUXQDKxp06Zyufwz8IMHD6pDhw76+OOPNXz48KB/4+wemw4dOujo0aNKTU0NO/Y1e4/phpfzgp4396YBIffYLF++XGnnD9DP/7EpYvdrpki3UW187TZq1CglJSWFdR921pA2jkTbReI5rvQaGv6XlX49NWdyScpMc+vjuy+3zLfxWLzugrXL2ezwmSBFpu2WbTukO+ZvrjEE6nt1PH1db8v18kVCvH/elZaWqlWrViopKanz+m1qj43b7Zbb7TblbyQlJTXohTGka2tlpaWosKS81vkFVR/GKRrStXW9Powv7nJOVO7XDNFqo9o09Pm0q0i0cUPaLhJ/P29PUZ0X76rhrAptPHjccnVsovm6S5I066oeta6KOpOdPhPOFG7bVXoNPbx0V8B5XS5JDy/dpbG92tmqPeojXj/vQn3MTB4OU7T2sHHS3jhOeixWZXYbR+LvHz4e2gTZUM9zEt92Ei2a1P6BHo/vo2CTqn3zutbuPRa7oGAptkps/vd//1fNmjVT9+7dJUljx45Vs2bN9Mc//tGUeKK1h42T9sZx0mOxKrPbuKF/v3XzlDpvr+95TpPTM0vrHxilu350vlo09k9w4vF9RCKMYGy1Kmry5MmaOHFijePJyckmRFMlWnvYOGlvHCc9Fqsyu40b8vcHdUoPaThrkA3mj0RLYoJLv/pRN90+4vy4fx+RCCMYWyU2jRo1UrNmzcwOo4Zo7WHjpL1xnPRYrMrsNg737/uGs6bP3SCX5JfcxONQS13Mfo6tgEQYwdhqKAqAM5k9nAb7MHteGazPVj02AJzL7OG0aGBj0OjwJcJn70+Wyf5kEIkNAAtx0lALG4NGlxMTYUQGiQ0ARFigXbl9G4MyvBYZTkqEETnMsQGACDJ7Y1Ig3pHYAEAEUUAOMBeJDQBEEAXkAHMxxwYAIogCcpBYEWcmEhsAiCAKyEWGnRMDVsSZi8QGACKISsoNZ+fEgBVx5mOODQBEGJWUw+dLDM6egO1LDHLzC0yKLDhWxFkDPTYAEAUUkKu/YImBS1WJwajsTEu2Y31WxFF/J3pIbAAgSiggVz92TwxYEWcNDEUBACzB7okBK+KsgcQGAGAJdk8MfCviAg2SuVQ1CZoVcdFFYgMAsAS7Jwa+FXGSajwGVsTFDokNAMASnJAYsCLOfEweBgBYhi8xOLuOTaZN6thIrIgzG4kNAMBSnJAYsCLOPCQ2AADLITFAuJhjAwAAHIPEBgAAOAaJDQAAcAwSGwAA4BgkNgAAwDFYFQUAqJdKr2HrpdhwNhIbAEDIcvMLahTPy7JR8Tw4H0NRAICQ5OYXaPrcDX5JjSQVlpRr+twNys0vMCky4AckNgCAoCq9hmYv3i6jltt8x2Yv3q5Kb21nALFDYgMACGrt3mM1emrOZEgqKCnX2r3HYhcUUAsSGwBAUIePB05qwjkPiBYSGwBAUK2bp0T0PCBaSGwAAEEN6pSurLQUBVrU7VLV6qhBndJjGRZQA4kNACCoxASXZk7IlqQayY3v55kTsqlnA9OR2AAAQpLTM0tzbuinzDT/4abMtBTNuaEfdWxgCRToAwCELKdnlkZlZ1J5GJZFYgMAqJfEBJeGdMkwOwygVgxFAQAAxyCxAQAAjkFiAwAAHIPEBgAAOAaJDQAAcAwSGwAA4BgkNgAAwDFIbAAAgGPEXYE+wzAkSaWlpSZH4s/j8ejkyZMqLS1VUlKS2eHYBu0WPtoufLRd+Gi78MV72/mu277reCBxl9gcP35cktShQweTIwEAAPV1/PhxpaWlBbzdZQRLfRzG6/Xq22+/VfPmzeVyWWdvk9LSUnXo0EEHDhxQamqq2eHYBu0WPtoufLRd+Gi78MV72xmGoePHj6tt27ZKSAg8kybuemwSEhLUvn17s8MIKDU1NS5fsA1Fu4WPtgsfbRc+2i588dx2dfXU+DB5GAAAOAaJDQAAcAwSG4twu92aOXOm3G632aHYCu0WPtoufLRd+Gi78NF2oYm7ycMAAMC56LEBAACOQWIDAAAcg8QGAAA4BomNhR08eFB5eXk6ePCg2aHYjmEY2rlzp3bt2hW0/DZqKisr06effqrdu3ebHYptVFZWaufOnfriiy90+vRps8OxrK+//lp5eXnVVeARGo/Ho/z8fH311VeqrKw0OxxLI7GxoMLCQg0dOlQXXXSRbr31VvXs2VNDhw5VYWGh2aHZwieffKKuXbtq1KhRuu6663TxxRdr3759ZodlK7fccouGDRumhx9+2OxQbOGPf/yj2rdvr0mTJiknJ0fnnXeeFi1aZHZYlnLy5ElNmDBB2dnZ+tnPfqbMzEy98MILZodleadOndKMGTPUtm1bXXfddbr88svVvXt3/fOf/zQ7NMsisbGgGTNm6NixYzpw4IDy8vK0f/9+HT16VA888IDZoVnezp07NW7cOP385z/XgQMHtHHjRj3//PP69ttvzQ7NNl5++WXt27dPl1xyidmh2EZZWZm2bt2qHTt26KuvvtL06dN13XXX6ZtvvjE7NMuYMWOG8vPztXfvXu3YsUMvvfSSbr31Vm3atMns0CytrKxMTZs21VdffaX8/Hzt379f48aN06RJk3TixAmzw7Mklntb0FVXXSW3260FCxZUH7vmmmvk8Xj4FhjEtGnTtGHDBm3dutVSe4HZxc6dO3XFFVfos88+04033qiOHTvq1VdfNTss2ykqKlKrVq20aNEiXXXVVWaHYzqv16tWrVrp3nvv1b333lt9vFu3bho3bpyeeuop84KzoW3btqlnz55as2aNBg0aZHY4lhN3e0XZwX333aerr75aTzzxhHr37q2NGzfq008/1bvvvmt2aJa3YsUK/exnP1NFRYW2bdum1q1bq3379iQ5IaioqNDkyZP12GOPqVOnTmaHY2vr1q2TJHXp0sXkSKxh3759+u6779S/f3+/4wMHDtTGjRtNisq+1q1bp4SEBN6nAZDYxMDBgweDzvHo0aOHWrZsKUnq27evrrnmGj344IPq3LmzvvrqK/3sZz9T7969YxCttezfv1/79++v85yLLrpIaWlpMgxDhYWF+uabb9S9e3elp6fr66+/Vrdu3TR//nx17NgxNkFbxN69e4MOhfTp00fNmjWTJN1999264IILNG3atFiEZ2lffvll0Dlt/fv3V+PGjWscLyoq0i9/+UtdffXV6tGjR7RCtJVjx45JkjIyMvyOZ2RkaMuWLWaEZFv79+/XPffco1tvvVXnnHOO2eFYEolNDPzzn//UnDlz6jznT3/6U/WchptuuklffPGFvv76a6WlpamkpEQjRozQzTffrLlz58YiZMtYtmyZ/v73v9d5ztNPP62+ffvK5XIpMTFRS5YsUV5enjp37qyysjKNHj1aN998s5YvXx6jqK1h8eLFevPNN+s856WXXlL37t21cuVKvfrqq5o/f74+/fRTSVJJSYkOHz6sTz/9VBdffLEaNYqfj4u33npLixcvrvOcefPm6dxzz/U7VlpaqnHjxikjI0OvvPJKNEO0laSkJElSeXm53/Hvv/9eycnJZoRkS4cOHdLo0aPVt29fPfHEE2aHY1nMsbGg5s2b68EHH9Rdd91VfezJJ5/UzJkzVVpaamJk1tetWzf17dtXb7zxRvWxOXPm6O6779bJkycZkgrggw8+0IMPPuh3bOvWrUpJSdH555+vJUuWKC0tzaTo7OH48eMaM2aMTp06pQ8//FAtWrQwOyTLKC4uVsuWLbVgwQJdc8011cfHjRsnt9utd955x8To7OHw4cO64oor1LZtWy1evFgpKSlmh2RZ8fMVzEbOOeecGrVrDhw4QLdjCMaMGVNjzP7gwYNq1aoVSU0dRo8erdGjR/sdGz58OJOHQ3T8+HHl5OTo1KlTWr58OUnNWVq0aKEBAwbovffeq05sSkpK9Mknn+jxxx83OTrrO3LkiEaMGKGsrCy99957JDVBkNhY0J133qn77rtP55xzjvr376/169fr2Wef1aOPPmp2aJZ33333qW/fvrrzzjt15ZVXatu2bXrqqaf02GOPmR0aHOr06dMaP368vvjiC73yyivatm1b9W1du3ZVZmamidFZxx//+EeNGzdOHTt2VP/+/fXkk0/q3HPP1Y033mh2aJZ24sQJjRw5UidOnNC9996r9evXV9+WnZ2t9PR0E6OzJoaiLGrx4sVasGCBvv32W7Vt21bXXnutJkyYYHZYtrBv3z79+c9/1q5du5SVlaUpU6Zo3LhxZodlO3fccYeysrL0u9/9zuxQLO3EiRPKycmp9bbf/va3mjhxYowjsi7ffMMjR46oT58+1V/gENj+/fs1derUWm975JFHdNlll8U4IusjsQEAAI5B5WEAAOAYJDYAAMAxSGwAAIBjkNgAAADHILEBAACOQWIDAAAcg8QGAAA4BokNAD/vvPNO0B3V401xcbHmz58vj8cjSTp69Kjmz58vr9drcmQAzkZiA8SZJUuWaOvWrQFvv/nmm/XZZ5/FMCLr27dvn6ZMmaKysjJJ0s6dOzVlyhSdOnUqpN8/fPiw5s+fL+qhAtFHYgPEmV/96ld+u5+f7eqrr9Z5550Xw4js55xzztHkyZOVmJgY0vnbt2/XlClTVFlZGeXIALAJJgA/Y8eOVbt27ap/XrBggYYOHSqXy6XNmzerRYsWGjRokBIS/L8XVVRUaPXq1Tpx4oR69uypjh07+t2em5ur4uJiuVwutWvXTn379lXTpk39zvH9rVOnTmnTpk3q0qWLLrroohox+s7zer3asmWL0tPTNXjwYL8d3Ou6r8OHD2vt2rVq0qSJ+vXrV+tu3Bs3btQ333yjCy+8sMZtGRkZ+vGPf1wjsTl69KjWrl2rpk2b6uKLL5bb7VZJSYk++ugjSdKbb76phIQEdenSRQMHDpQk7dmzR1u3blWLFi10ySWXKDk5ufr+CgsLtXLlSl177bVat26d9u/fr5EjR6ply5Y1YgLwbwaAuNKlSxdjxowZAW/PyMgwXn/99eqf3W63MXbsWKNLly7GlVdeabRu3doYM2aM4fV6q8/5/PPPjbZt2xoDBgwwxo8fb6Snpxt33XWX3/3+13/9lzF58mTj2muvNfr06WO0bdvWWLt2rd85brfbyMnJMTp27GhMmjTJeOONN2qN0e12G2PGjDHat29vjB071khPTzfGjx9veDyeoPf1+OOPG2lpacbo0aON4cOHG+np6cZ7773nd/833XST0axZM2PMmDFGp06djLFjxxqSjO+++84wDMNYtWqVIcn4/vvvq3/nz3/+s9G4cWPjkksuMYYPH2706tXL2Ldvn3HgwAHjiiuuMCQZP/3pT43Jkycbzz//vGEYhvHrX//aaNq0qTF69Gjj/PPPNzp27Gjs3Lmz+j6XLl1qJCQkGGPHjjX69+9vTJ482dizZ0/A5w6AYZDYAHEmnMRm2LBhRnl5uWEYhrF//34jOTnZ+OCDDwzDMIyysjIjMzPTePHFF6t/5+DBg0ZGRoaxePHigH/n/vvvNy6++GK/Y2632xg4cKBRVlZW52Nwu91G9+7djeLi4uqYWrZsabzwwgt13teKFSuMli1bGl988UX1sfnz5xvp6elGSUmJYRiGsWTJEiM5OdnIz883DMMwvv/+e+Piiy+uM7FZtmyZ4XK5jKVLl1bf786dO40dO3YYhmEYH3/8sSHJL/FatmyZkZiYaOTl5RmGYRgej8cYO3asccUVV1Sfs3TpUkNSnc8XAH8MRQEI6qabbpLb7ZYkdejQQZ07d9auXbs0atQo5ebmqqioSKmpqVq4cKGMqi9M6tSpkz7++GNdeeWV1fezb98+7dy5UyUlJWrSpInWr1+vyspKvyGdX/ziF2rSpEnQmG6++WalpaVVxzRlyhS9+eabuvnmmwPe16uvvqoLLrhAmzdv1qZNm2QYhrxer4qLi7V582ZddtllevPNNzVu3Dj16NFDkpSSkqI777xTn3/+ecBY/v73v2vkyJHKycmpPta9e/c643/99dc1ZswY9e/fX5LUqFEj3XvvvRo+fLgOHz6s1q1bV5975513Bm0PAFVIbAAElZ6e7vez2+1WeXm5pKpkJTk5WW+//bbfOV26dFHnzp0lSYZh6Be/+IXefPNNDRo0SBkZGSotLZXH41FxcbEyMjKqfy8rKyukmM6ew9OpUyctX77c79jZ97Vv3z4VFRVp4cKFfsevvfba6rkt+/fvV+/evWvcd13279+vXr16hRS3z9dff12dPPl06dKl+jZfYpOYmOiX5ACoG4kNgAZJTU2Vx+PRa6+9pqSkpFrPWbVqlebOnas9e/aoffv2kqQPP/xQy5Ytq7EE+swJwHX57rvvavzcqlWrOu8rNTVVrVq10vz58wPeb0ZGRq33XZcWLVqoqKgolLCrtWrVSseOHfM75vv57McBIHQs9wbQID/60Y9UWVmpV155xe+4x+PRkSNHJFWt7mnevLnatm1bffvZvSb19e6771b/2+v16t1339XQoUPr/J2cnBzl5uZq3759fscLCwuri+1deumlWrZsWXWPlKQavVFnGz16tJYuXVr9eCWpsrKyOiFq1qyZJPndp+/v+GrjSFUruTp06KBzzz23zr8HIDB6bIA4tH379hq9Fr169VJ2dna976tjx47605/+pNtvv12bN29Wv379tH//fi1cuFBPP/20RowYoWHDhqmyslLXX3+9Ro0apVWrVvklJuFYs2aNrr/+eg0bNkxvv/22jh49qt/+9rd1/s4tt9yiRYsWaciQIbrjjjt0zjnnaNOmTVq+fLm2bdumhIQE3XLLLXrmmWc0cuRITZs2TZs2bQoa6y233KIFCxZo8ODBuu2225ScnKzXX39df/3rXzVo0CB169ZNzZs31/33369LLrlEXbt21X/+53/qxRdf1PDhw/Uf//Ef2r17t55++mm9/vrrIdfHAVATPTZAnLnyyiuVnJysd9991++/Xbt2SapZoO+nP/2pX10bqarn44ILLqj++e6779ann36qpk2batWqVUpKStLixYs1YsQISVKbNm20Zs0atWvXTv/85z/VvXt3ffTRR5o8eXL1pORAfyuQF154QUOHDlVeXp569eqldevW+c1Fqe2+kpOTlZubqyeffFL79+9XXl6eLrroIm3atKl6GK1x48ZavXq1Ro4cqbVr16pr16765JNPNHny5Op5OGcX6HO73VqxYoV+//vfa8eOHfr666/1zDPPaNCgQZKqhsBWrFihxMRELV68WBs3bpTb7dZnn32mKVOmaPXq1Tp9+rRWrlypn/zkJ9XxZmVlafLkySG1B4AqLuPsAW4AsLiUlBQtXLjQb8UVAEj02AAAAAchsQFgO/UZsgIQXxiKAgAAjkGPDQAAcAwSGwAA4BgkNgAAwDFIbAAAgGOQ2AAAAMcgsQEAAI5BYgMAAByDxAYAADgGiQ0AAHCM/x//4BIY1iL0WAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }