Skip to content

Instantly share code, notes, and snippets.

@pr0nstar

pr0nstar/covid.alt.soliz.ipynb

Created December 31, 2021 14:57
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save pr0nstar/a9e7dcd4dbb03e376042ea27b8713ab4 to your computer and use it in GitHub Desktop.
Save pr0nstar/a9e7dcd4dbb03e376042ea27b8713ab4 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
covid.alt.soliz.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "smaller-bookmark",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"import datetime\n",
"import unidecode\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"\n",
"import sklearn.linear_model\n",
"import statsmodels.formula.api as smf\n",
"\n",
"from matplotlib import pyplot\n",
"\n",
"warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "pregnant-sweet",
"metadata": {},
"outputs": [],
"source": [
"soliz_df = pd.read_excel('https://github.com/daliagachc/cov-alt/raw/master/data/soliz_2021.xlsx')\n",
"soliz_df.columns = [_.lower().replace('/', '_').replace(' ', _) for _ in soliz_df.columns]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "identified-carroll",
"metadata": {},
"outputs": [],
"source": [
"soliz_df = soliz_df[~soliz_df['altitude'].isna()]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "environmental-infection",
"metadata": {},
"outputs": [],
"source": [
"def get_binned(bin_width, stat):\n",
" soliz_binned = soliz_df.assign(count=1)\n",
" soliz_binned = soliz_binned.groupby(\n",
" np.round((soliz_df['altitude']) / bin_width + .4999) * bin_width\n",
" ).agg({\n",
" 'cases_pop_den': stat,\n",
" 'count': 'count'\n",
" })\n",
" soliz_binned = soliz_binned.reset_index()\n",
" \n",
" return soliz_binned"
]
},
{
"cell_type": "markdown",
"id": "proved-monthly",
"metadata": {},
"source": [
"### Original"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "functional-violin",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>OLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>np.log(cases_pop_den)</td> <th> R-squared: </th> <td> 0.603</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>OLS</td> <th> Adj. R-squared: </th> <td> 0.594</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 63.82</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Fri, 31 Dec 2021</td> <th> Prob (F-statistic):</th> <td>5.80e-10</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>10:52:34</td> <th> Log-Likelihood: </th> <td> -85.227</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 44</td> <th> AIC: </th> <td> 174.5</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 42</td> <th> BIC: </th> <td> 178.0</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 1</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 8.5537</td> <td> 0.505</td> <td> 16.934</td> <td> 0.000</td> <td> 7.534</td> <td> 9.573</td>\n",
"</tr>\n",
"<tr>\n",
" <th>altitude</th> <td> -0.0016</td> <td> 0.000</td> <td> -7.989</td> <td> 0.000</td> <td> -0.002</td> <td> -0.001</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>14.013</td> <th> Durbin-Watson: </th> <td> 1.693</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.001</td> <th> Jarque-Bera (JB): </th> <td> 16.319</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td>-1.096</td> <th> Prob(JB): </th> <td>0.000286</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 5.024</td> <th> Cond. No. </th> <td>4.91e+03</td>\n",
"</tr>\n",
"</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 4.91e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" OLS Regression Results \n",
"=================================================================================\n",
"Dep. Variable: np.log(cases_pop_den) R-squared: 0.603\n",
"Model: OLS Adj. R-squared: 0.594\n",
"Method: Least Squares F-statistic: 63.82\n",
"Date: Fri, 31 Dec 2021 Prob (F-statistic): 5.80e-10\n",
"Time: 10:52:34 Log-Likelihood: -85.227\n",
"No. Observations: 44 AIC: 174.5\n",
"Df Residuals: 42 BIC: 178.0\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 8.5537 0.505 16.934 0.000 7.534 9.573\n",
"altitude -0.0016 0.000 -7.989 0.000 -0.002 -0.001\n",
"==============================================================================\n",
"Omnibus: 14.013 Durbin-Watson: 1.693\n",
"Prob(Omnibus): 0.001 Jarque-Bera (JB): 16.319\n",
"Skew: -1.096 Prob(JB): 0.000286\n",
"Kurtosis: 5.024 Cond. No. 4.91e+03\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 4.91e+03. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"original_data = get_binned(bin_width=100, stat='sum')\n",
"original_model = smf.ols(\n",
" 'np.log(cases_pop_den) ~ altitude', \n",
" data=original_data, \n",
").fit()\n",
"original_model.summary()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "occupied-study",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PolyCollection at 0x7f45ef7c9280>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4PUlEQVR4nO3deXyVZ5nw8d+dQMhGSAghkBAgLKWlC7Sle21r94Uuo9bq1LFuU5dx3D5O1dHXGf2M7zjWcV59tY7VUes4b53RLpbQBVrZ1EKBQoFCKSEh7CQsgSwnOTnJ9f5xPYdzCElOzn6Sc30/n3zIeZKc53nS9Fznvu7rvm4nIhhjjDFDyUn3BRhjjMl8FiyMMcZEZMHCGGNMRBYsjDHGRGTBwhhjTERj0n0B0Zo0aZLMnDkz3ZdhjDEjysaNG4+KSEWsPz/igsXMmTPZsGFDui/DGGNGFOdcUzw/b2koY4wxEVmwMMYYE5EFC2OMMRFZsDDGGBORBQtjjDERWbCIoLXTn+5LMMaYtLNgMYSO7gB1Ww7R3h1I96UYY0xaWbAYwvaDp/jas9vYcfBUui/FGGPSyoLFIDq6Azy+ZjcAP13TYKMLY0xWs2AxiO0HT7F8ezMAy7YfsdGFMSarWbAYQPioIshGF8aYbDbiekOlQnt3Dx+6upaHrqo9fcw56OjuoXic/cqMMdnHXvkGUFlSQGVJQbovwxhjMoaloVLI1mwYY0YqCxYpYms2jDEjmQWLFLE1G8aYkcyCRQrYmg1jzEhnwSIFbM2GMWaks2CRZLZmwxgzGljpLFqlVFqYl5TntjUbxpjRIOtfrYJVSvddXJ2UF29bs2GMGQ2yPg1lVUrGGBNZSoKFc+7nzrlm59y2sGMTnXPLnXO7vH/LUnEt4axKyRhjhidVI4tfArf3O/Zl4BURmQu84j1OKatSMsaY4UlJsBCR1cDxfofvBZ7wPn8CuC8V1xJkVUrGGDN86ZzgrhSRQ97nh4HKwb7ROfcw8DDA9OnTE3Jyq1Iyxpjhy4hXRRER55wM8fXHgccBFi1aNOj3RcOqlIwxZvjSWQ11xDk3FcD7tzmN12KMMWYI6QwWzwEPeZ8/BPw+jddijDFmCKkqnX0SeBWY55zb75z7KPBt4Bbn3C7gZu+xMcaYDJSSOQsRef8gX7opFec3xhgTn+xawd3ZCX196b4KY4wZcTKiGiplNm6EQADOOQcqK2FMdt2+McbEKrteLQMBDRBbtui/s2dDdTWMG5fuKzPGmIyWXcECIC8PJk/WwFFfD2+/DdOnw4wZUFSU7qszxpiMlH3BImjMGCgv1zmMgwehqQmmTIHaWigtTffVGWNMRsneYBGUkwNlZSACra3w6qswYQLMnavBJGfgGoBkbphkjDGZJruqoYbiHJSUaIqqtxc2bIDVq+HAAU1ZhQlumGRNB40x2cKCxUAKCzVo5OXpZPiKFdDQAN3dgG2YZIzJPpaGGsq4cWdNhndOqeaJ17Tb+k/XNHBeVYl1qTXGjHo2shiO4GR4eTl7t+6ibflKzj9cz9rXd7PjwMl0X50xxiSdBYso+AJ9/LahnWOFEyju7uSSAzt45RfP0rHvoK0MN8aMapY/iUKnP8DdC6pZfFH16WO5XT66171GUfkEmDNHy2/Hjk3jVRpjTOJZsIhCeXE+5cX5/Y6W6j9+P7z5JuzYoWs1pk2DAttcyRgzOliwSJS8PKio0MnwhgadEJ82TVeGl5Sk++qMMSYu2RMs2tp04V2yjRkDkybpHEZzM+zbp49nz4aJE3U9hzHGjDDZEyw+8xlYtgwWL4a774apU5N7vpycUNuQ9nZYtw6Ki3Vl+OTJkJub3PMbY0wCZU811O2364v0449rsPjkJ2HpUvD5kn/u4mJtiZ6bC5s36yK/PXtOL/IzxphM5yQVqZkEWrRokWzYsCG2H16xQvs/vfQS1NVpK4/CQrj5Zg0gCxfGnSZq6+phfH6EaqhAQK8DoKZGu94WF8d1XmOMGYpzbqOILIr557MuWBQW6ryCiL7LX7IEXn5Zd9GbNk3TVHfdFVOayucPsGJnCzfMq6AwbxgZvr4+DRqBgI56Zs3S1JXNaxhjEsyCRTTCg0U4nw/+8Acdbaxfr8cuu0wDx403DrsE9s0DJ/nS01v5l3ddyPnVE4Z/XSI6r+HzaeXUnDk6KW7zGsaYBIk3WGTPnMVQCgp0NPHjH8Nzz8EnPqF7XPzDP8Btt8E3vwmvvz5kNZXPH+DpTfsBeGbTATr9UXSkdQ7Gj9fRhYiea9Uq3WPD74/37owxJm5pr4Zyzn0e+BggwFbgwyLSlbYLqqqCj30MPvpR2LRJRxsvv6xBpLo6lKaqqjrjxxpaOljXeAKAtY3HaWzpiG50EVRQoB89PfDWW7Bzp85p1NTYTn7GmLRJ68jCOVcNfAZYJCIXALnA+9J5Tac5B5dcAl//uk6If+MbGiB+8hO45x4dfdTVgc93xqgiKOrRRX9jx2oqqrRU12qsXq3Bq7U1NetFjDEmTFrnLLxgsRZYAJwCngV+ICLLBvuZpMxZROPQIS25rauD/fuhsJCu629g3zU3037uBad31nMOasoGag8SIxFdWOjz6U5+Nq9hjInCiJ/gds59FvgW4AOWiciDA3zPw8DDANOnT7+0qakptpMlIlgEicAbb2g11fLlWk1VXa0pqrvu0s+TxefTwDFunK4MnzpV240YY8wgRnSwcM6VAU8BDwCtwG+B34nIrwf7mbSPLAbi8+lzB6upRODSS3Xtxo036jmToacHTnr7aUyfrh82r2GMGcBIDxb3A7eLyEe9xx8ErhSRTw32MxkZLML1T1MVFOiiv7vu0jmQnCRME/X2wqlTGjxsvYYxZgDxBot0V0PtBa50zhWiaaibgBgjQYaYOjVUTfXGGxo0li/XdFVVlQaNxYsTm6bKzYWystC8xtq1Woo7e7b1oTLGJEQmzFl8A01DBYBNwMdEZNCmSRk/shhIV9eZi/5EdJRx991w000xpakithXp6tLRxtixOtKoqoL8BE22G2NGnBGdhorFiAwW4Q4fDqWp9u3TNNVNN+loY5hpqqjaigQCOq/R16drNWpqbH8NY7LQSE9DZZ8pUzRF9ZGPnJmmqqsLpanuukv7VA2ioaWDx1buZsbEwsgL/8aMgfJyDRaHD8PevbqvRnB/jWTMoRhjRp3sGlls3qwT0M7pO/qiosyYBO7qClVTvfZaKE21eLGOOsIqnHz+AN9dtpN1jSe4snYiX7j1nOE1LQzX0aEf+fm2b7gxWSJlaSjn3LuAfwEmA877EBFJaU4jrmABWuba2qq9n44e1XfceXnaIjyd6amgw4fh+ec1cOzdqy/oN92k8xuXXMKbh9r40tNbT3971E0Lw/n9mqJyDmbO1NGMld4aMyqlMljUA3eLyI5YT5YIcQeLcIGATgIfOaIjju5uTcsUFQ2702zSiMCWLRo0li2Djg76pkzhj/Ou4FcTL+RwySSA2EcX4fqX3tbWanVVJoy6jDEJkcpg8ScRuSbWEyVKQoNFOBFNzRw/rpsiBTcnys/X4JHE8tNhVTatXIn/988xdsN6nAjt8y/k2Dtv5eQ111NdNTExbUWCv4POTr3nOXOgosJSVMaMAqkMFt8HpqD9m06XtorI07GePBZJCxb9dXdriubwYR15BAKapiouTmhrjag3TDpyJFRNFUxT3XijpqkuvTRxE9bB0tvcXB1pVFcnbyW6MSbpUhksfjHAYRGRj8R68likLFiE6+vTF86jR3Wuo7MzNEleWBhXuiauDZO2bg31pmpv14nq4KK/mpqYr+kMwVRdcDc/S1EZMyLZOot06Ow8c5JcRJv6FRVFNUmekMom0FHAqlUaONat0+tZuFCDxs03J2Z/7/Dd/CxFZcyIk8qRxTnAj4FKEbnAOXcRcI+I/FOsJ49FRgSLcMFmfs3NOknu92vqpqgo4orp4KgiKK7KpqAjR+CFFzRwNDVpEAumqRYtSkyaKjxFNXOmpqhSWEXV2umntDB9XXbTfX5jYpHKYLEK+DvgJyJysXdsm7dpUcpkXLAI19en776PHdNJ8vZ2PT7Amo7wUUVQQiqbgkRg2zYNGsuWJSdN1T9FNXOmLvRLYoqqozvAM5sOcN/F1RSPS32pc7rPb0ysUhks1ovIZc65TWHBYrOILIz15LHI6GDR3xBrOo51Bdh3ouuMTe8SvmFSUDBNVVenaaq+vsSmqcKrqAoKNEU1eXJS9thY33ic+3/yKr/9+FVcVjsx4c+f6ec3JlapbPdx1Dk3G90rG+fce4BDsZ44KwT30546NdSjqbkZDh6k3O+nvHB46aq45efDbbfpR3NzaNHfP/0TPPpo/Gkq5zTgFBdrFdm2bV7k83pRjR+fkNvo6A7w+JrdAPx0TQPnVZWk9N19us9vTDpFM7KYBTwOXA2cABqBD4jInqRd3QBG1MhiMMFW4v3TVYWFcVdXRXUNb76paaqXXkp8mip8oV9ZWagXVRzrVYLv6oNS/e4+3ec3Jh4pr4ZyzhUBOSLSFutJ4zEqgkV/Ph+cOKET5C0tMVdXxay7G1av1sCxdq2mqRYs0NFGItJUwV5UeXlaejt1atQr5Du6A3zuvzexfHvz6WO3zq/kew8sTMm7+3Sf35h4JT1YOOe+MNTXReR7sZ48FqMyWIQLr646eFAf5+ZqKmcYcwCDrQaPuEo8qKUllKZqbAxVUy1erGmqeFayB1Nxvb0aMGbMGPaOfkdO+ahv7jhrjmfO5CIqS5LfmiXd5zcmXqkIFv/gfToPuAx4znt8N/CaiHwg1pPHYtQHi3B9fZquOnpUt2jt7Az1rsrPP+tFdrDV4FGvEodQmqquTtNUbW1QWRlKU02fHvt9BdNwXV16L7NmJW1CPB5WImtGk1RWQ60G7gqmn5xz44GlInJdrCePRTYFi7NerDo6NF114ID2sHIu1LsqJ2fQ1eAxrxIPGixNtXgx3HJLfGmqri4NHKABaNq0jNicyUpkzWiTymqoSsAf9tjvHTNJ0NEdoG7LoTNfrIqK9GPatFDvqoMHobkZX5efF/7YRG5fL89sOkBtRRGFeWPw+QM8vWk/wBnHozJunAaFW245M031rW/Bd78L73ynBo7LLos+TZWfrx99fTpn09SkwWLWLFqLJlBakp5+VNsPnuJrz25jXuV4m8Q2huiCxa+A15xzz3iP7wN+megLMirii9W4cZq6mTwZenvZsWUPLxzfT6Xfx65tDeyrLWTenKk0HO06vfBvbeNxGls64lslXlEBDz0EH/wgbN8eSlO9+KKmqe68UwPHjBnRPW9Ojs5fAPh8+NZvZENjK1fdeAlFM4dffpuI1JGVyBpztmEX1YvIt4APo2WzJ4APi8g/B7/unCtL/OVlp/4vVu3dgaG/PyD8eOtxdk+q4c8zFrBx2nx+2ZzLyVOdLF+1lYmdJxnXo42Cn9l0gE7/0M83LM7B+efDl76k7UW+/W1djPfEE/Dud+u2sc88EyoLjkZBAQ0U8oONLezbtAP++EdNfx0+rJPkgwiOxiL9viLZfvDU6aqnZduPsOPgqbiez5jRIKq3SyLyOvD6IF9+Bbgk2gtwzpUCPwMuQBf8fUREXh3yh0a5gV6shkqFtHf38KGra3noqtrTx5yD46X5XFE1l6tbW8k7fJAxp05BDnSdOEnh5AS25Rg3Tktsb75ZJ+ODvamCaaobbtAy3GGmqYKps76cHH7X2MkX5k6jsMcPmzZpKfHMmVpN1W+uJBGpo/BAHWSjC2OiDBYRxPrK833gRRF5j3MuD8jqTRNiebGqLCkYtHyzdlIxMANYoPMcJ07oPEdLi84TjBunL7oJ2typrXgC4//qr+ADH4AdO0KL/l56SVNmd92lHzNnhn6mX1lvQ0vHwKmzwkIdWezZA/X1uthv5kyYNImOXhKSOhos8HZ091iwMFktYS3KnXOvi0hUIwvn3ARgMzBLhnkho70aKmX1/MH1HIcP68RyIKDtxsePP70QcKi1GQN9bdASXb9fq6nq6uDVV3WdxUUXweLF+K5/JysOdp3+magaLHZ2aporN5ctuaW8v66JjnH6XsNWVxtzpozZzyLGYLEQbSGyHVgAbAQ+KyIdg/3MaA8W6dDa5qM00KXtzQ8cAL8fX5+w6mA3111QddaL9GBBYVgluuFpqoYG+vLyWFN1PjMeeoCZd76TY76eqBss+jq7eOy5TWzZc4y2cUXsm1DJokvn8uiDi2w0YIwnk4LF6W60UfzMImAtcI2IrPO2bj0lIv+r3/c9DDwMMH369EubmpoScs1mgPUEInDqFG9squd//2QZ37x1NvOqJmiqymt4OFBQiLSR01kjERG6tr7J1h//J/O2vEpJdyd9FRXkBBf9haWpIjnW3nU6wOR0+cjt7ICcHKacP5tJ82ZrKa7t7GeyXCrXWeCcuwS4Fp2I/pM34R10Uwzn3w/sF5F13uPfAV/u/00i8jg6AmHRokUja2u/DHfWpLBzdOQX8cO9faybfiE/7Cnh27OqKTpySNdzBPpYsu4gcOa6jUHnGdBAsmbX0TNHIs6xu7yGb8y/mzHz7uDy/W/ytx07GP+f/wm//CVceKEGjVtvjVg2W14cPuoo1X96e6HtBPz5z6FV4hUVOkdjjInasEtnnXNfB54AyoFJwC+cc18Lfl1Ejkd7chE5DOxzzs3zDt2EpqRMCgxWonu6Gss5ljS0sX1MGVx7LVx/PTvLZ7Bu70nKO1t5+81Gmppa8HX3nF74FxReotvQ0sFjK3fT2BLKLoYvFgzkjuHPMxbw/bs+Reezz8HnPqer1f/5n7Wt+t//fWiuY7hyc3XdxuTJ+vm2bfCHP8Dmzdrtt68v5t+bMdkompHFg8ACEekCcM59G52cjndb1b8F/surhGpA13KYFBioRHd+Vcmg1VguN48f1fvYVH0e4wJ+SnxtPP3WMR7OCfAX0/K4Z04tfQXaYt05DQgOBlxB3ukPcPeCahZfVH36PM6Bb3w+hR/4ADz4oFZTBRf9LVumI4M779Qy3CjSVKdXiYuEuvvm5elzTJmSsC1hrZeUGc2i6Q21AvgLEWn1HpcCT4vIjUm7ugHYBHdiDNZy++uL59N0vHPAaixg4EqtCWOoDPhCJbkiWuZaVMSbB0/Fv8+43w9r1mjg+POfdYRxwQWhNFUsvaTCO+CWlp4uwWXsMDrzDsB6SZlMl8pGgs+iXWeXo3MWtwCvofMOiMhnYr2IaFiwSIyklej6/bqV7P79+PYf4ier6nn1UCedY/MRlxP/PuNHj2prkSVLYPduHSFcf72ONq64Irb1IsES3JwcqKrS3lsTJkS1a6Btt2oyXSqDxUNDfV1Enoj1IqJhwWLkOHLsFE31Bxhz+BB5x45CnyAF+VRVTaQ83gaBIrBzpwaNF1/UUcKkSaHeVLNmRf+cwZbw3d2atpoxQ/tdRUhThY/SbEMkk6lSWjrrzSuc4z3cKSI9sZ44VhYsRqieHh1xHDig6zn6+nS3PK+9elz8fu0ftWRJKE11/vk62og1TRVctNjXp6OMYJpqgD03bLtVMxKkcmRxA1oNtQdt7VEDPCQiq2M9eSwsWIwCAwUOb44j7vUQx46FFv3t3q1zEDfcoKONK66IbZtany/UEHHKFE1TlZVBbq5tt2pGjFQGi43AX4rITu/xOcCTInJprCePhQWLUaanRyuU9u8P9atKROAYLE11xx064oglTSWiQcPn06BTU0NzUSm7unKRsNZott2qyUSpDBZbROSiSMeSzYLFKOb3nxk4QANHYWF8gaOnJ5Sm+tOfNE01f34oTTUhhv09AgGd3+jp0fmN2lot7U1QGa4xiZbKYPFzoA/4tXfoQSBXRD4S68ljYcEiS3R369ax+/drBZRz2nKkIM5368eOhaqp6us1TXX99ZqmuvLK2NJUfj+cOqWjopISnRifNOl0exRjMkEqg8U44G/Qdh8Aa4DHRKQ71pPHwoJFFurq0hf5vXt1riMn54xeVcMxUG8qdu7UtRsvvKBpqvLyUDXV7NmxXWtwfkNERxo1NTBxYszrN4xJlHRUQ52HjjB2iog/wo8knAWLLNfZqSONvXs1DZSbq+/mh3gxHrR1elBPj6anlizRdFUwTbV4sbYbiSVNJaItSzo7NbhNmQLV1acnxo1JtVSOLO4C/h3YjVZD1QIfF5EXYj15LCxYmNPa26G5GZqadPTRbz+OoGG1Tg86flzTVHV18Pbb+pzXXafzG7Gmqfr69Fq7ujRQVFfrTn+lpfGXDRszTKkMFm8Bi0Wk3ns8G1gqIufGevJYWLAwZ/HaqnPokM5xBCedi4vxBfqGbJ0+pPA0VWurpqnuuENHHHPmxHatvb06IvL7NRDV1OjCvwkT4i4btt5UZiipDBbrReSysMcOeC38WCpYsDBD6u0NreE4dIidB0/ytZcb8Y0ZB87F1psqmKaqq9MeVb29cN55oTRVaWls1xpeUTVunAaOyZNj2n/DelOZSFIZLH6Mbub8P2hvqPuBvcDLACLydKwXEQ0LFma4Otp9/K+fr6R+41uU+U4hOOafU81n7rog9t5UJ06EqqneflvTUsE01VVXxZamAg0cp07pv/n5MH26TpCPHz+swGG9qUwkqQwWvxjiy5KqEloLFma4wpsl5nT5GHPiOPn7m6gZ10dZSZG+g4/1xR0SkqYacJ/znh4dcQQDx4wZGjiKiwcMHNabygxHJm2r+hUR+eeEPNkQLFiYuIjoC/GRI9DURFu7j/ETivWFONbJ5mCaaulSWL162GmqiFVaoHMbbW36nAUFoRFHWOAYqb2pbI4lteINFoksxbg/gc9lTHI4pyOKuXPpuOY6XiiooaN4gq7jaGnRdRLRCvafevRRTVF98YtaAfXoo3D77fDIIzrXEQic8WMD7SB4lrw8HbFMnqzzGvX1Wt67ahU0NNDRcpzHV9ef8SPhux5mqo7uAHVbDmX8dZqQRI5V4+wAZ0xqbT/cziMrDlD78au4bMEFun6jsVHLcXNztUIp2jRVWRm87336sWuXpqmef163dA1LU/mmzxxwB8EhBQMH6Iijvh5/m4/P5AoP31xNoHwSvUXFuBxHR3dPRqeiztr73WS8RP41JSafZUwK9N9//LwHFlJcVaWbH7W3w+HDun7D79d5g2FONJ9h7lz4/Ofhb/9WW6cvWQK/+Q38+tcw+xwqyi9gfO2lrG2ExpaO6Kq0vMBRVg5lPT3Q1gqHj+q1TpsGMhb6xmXkOo6zfvdVJRkd2IyykYXJSgPtP376HW5xsU5Qz5ql1U/79mnwAA0a0fZ8ClZMXXcdtLbiX/o8rb95ik/sfpqPbvg9r007n23+m6j95LspLIyhn9TYsdpSBHT+ZM8eTVfl5WngmDw56p3/kmnI373JWIkMFr9N4HMZkzTh72yDBnyHm5OjaZ/ycm1s2NKiL8TNzfoCXVISfeuO0lLa7n0XR667k9bGBspXLOPy1S9zzZP/Rt+Lvwz1ppo7N7abCw8cgYAGuoYGPV5VpW1HJkxIW8uRYf/uTcaJpnT2O8A/AT7gReAi4PMi8ushfzDB4qmGsuoLA3HuPx6spjpwQF+Ie3u1LXk8rckDgVCaKjgRPm+ert24/fbYF/31P0d7u6bVcnM1aFRV6XPHUz4cpaTt/W4iSuU6i80istA59xfAYuALwGoRWRDrycOeOxfYABwQkcVDfW+swcJWuJqECwS0imrPHu0pFeukeLjWVnjpJQ0cb72lz/WOd+ho45prEvPC3turTQ67unT0VFERanI4wLaxZnSIN1hE85cX/N67gN+KyEkX7xaYIZ8FdgAxbJY8PFZ9YRJuzBjt61RZqS++hw9r4PD7Y9/tr7QUHnhAP+rrQ9VUK1boi3lwp78o0lRnLfwLduotKdES37Y22LRJv1ZeHgoc8e4dYkaVaGa86rxmgpcCrzjnKoCueC/AOTcNDUA/i/e5BtO/+sJqu03CFRXpHhg33ACXXaaT5C0tWo7b0xPbc86ZA5/7nAaLf/s3uPhi+J//gfe/Hx58UCurWluHfAqfP8CaXUfp9A/yNx/cG6SiQjds8vlgyxbaX1quCw337g3tzxGj1s6U72QwLJl6XZlq2MFCRL4MXA0sEpEeoBO4NwHX8H+AR9A9MgbknHvYObfBObehJbjdZhQGqr4wJilyc/VF99JLNXDMnasvwM3N+g4+lhfdYCrqO9/RRX+PPKLn+e53dU7ji1+ElSvPWvQHw1z4F+QcFBXhm1DG6tYcDTBvvaXzKKtWwe7duklU36D/q54lUxffZep1ZbJhBwvnXCHwKeDH3qEqIOb8l/eci4FmEdk41PeJyOMiskhEFlVUVER1jsGqL+yPxCRdQYGW315/PVx+uZbdtrToPEeE0UZb1yBfLy2F974XfvUrHVm8//2wdasGjDvugH/9V21wiI4qwhf+DTq66Od0gGnrPXP1+O7dOhG/YgXs2KHzNL29Qz5XMP2baW/QhrouG3EMLJo5i18AG9HRBcABtFy2Lo7zXwPc45y7E8gHSpxzvxaRD8TxnGdo7+7hQ1fX8tBVtaePOUfGr3A1o0h4Ca7Pp/tuNDYOOrcRTB0N2TMKNE312c/C3/wNrF2r8xu/+x08+SSccw6t193MjrYqyC9mbePxYS386x9gTq8sD189HgjAwYO6aDEnR4NJsLIqbII8UxffDXVdwRGHFcKcLZrfxmwRecA5934AEel0cc5wi8hXgK8AOOduAL6YyEABUFlSYCV5JnMERxszZ+o786ams9qLBN/Zz5hYOLxV3WPGwLXX6sfJk7BsGX3PLWHqzx7jVy6H9dPm88qcy/n9+hJqK+YPGYAaWjpY13gCYPAAM2ZMqJy3r0/nTQ4d0oBXXq6Bo6yM7Ue6MnLx3VCLAq0QZnDRBAu/c64Ar62Ht1Ned1KuyphR5qw1Pjk5OrcxaZLu0+2NNnwdPurWHQSi6BkVbsIEuP9+TtxxNy1vvEXZimUsWvkyV634OYGNT9Gz41Z41326jqOf8FFFUMRryMnR9Nr48Tof4/PBtm34unt4Zd0hprXm0FpQQnteQUaMLoZaFOggI0dCmSKadRa3AF8D5gPL0BTSh0RkZdKubgDWotyMNMNe49Pby6bX6/nH/7uUCb42Ajm5fPV9V3D+9DPf4Q64B8ZQAoFQmmrVKp0vOeccXbtx++2nV3wfa+9i34musxbM1ZTlU14cXRuSY+1d7D9yipyOdlxvHzJ2DP7KqdScO4PJNVPStoJ8qEWBe4/5RmSr9+FK6X4Wzrly4Eq0D9RaETka64ljZcHCjDTD3cUufBOjQr+PyrZj3FPq5xPvqKWgvAzy84e3B8ZQvDQVS5bA9u36on3ttRo4rr1W24IkQ/8V5MF5jgkTdPI8zcJ/90GjbSOpVK7gvgbYLCIdzrkPAJcA3xeRplhPHgsLFmYkiWYXu4He9eYEepib66Oi+SC0tfHWsS7+bvkevv3uBdHvJd7f7t2hRX/HjukLd3Cnv3nzol9QOFx9faEV5KDzH9XVOsIpLEzeeYeQDW1IUhkstgAL0J5QvwD+A3iviFwf68ljYcHCjCQJ28VOhI7mY3znZ8tofGMXl8yYyMN3X0xhUQJeyAZKU82dq0HjjjtCjQmTxefT4NHXpwUA1dW6SLCkJGM65Y4GqWz3ERARcc7dC/xIRP7DOffRWE9szGiX0A6rzrG9M4cn2iaQN30B9e3HuX3fMc6dmKfvxouLY7/QAaqpWLJEV43/4Afak+ruu5OXpiooCLUW8fu1Zcru3ZquqqyEqVN11GN9q9IqmpHFKrTb7IeB64Bm4A0RuTB5l3c2G1mYkSKRqY2Bcuq3nVvBv94yneIDezWNFGvb9ME0NITSVEeP6gv27bdr4EhmmioovOFheFluaWl8XX6zVCrTUFOAvwTWi8ga59x04AYR+VWsJ4+FBQuTjSIGnvZ2bZve1KTpnFg2aRpMIADr1mngWLlS01Rz5oTSVMHFeskULMvt6NDPCws1XTVpkqWrhiml1VCZwIKFMUPo6dFFfrt36wtrQYGmqBI1Cjh1StNUdXWwbZuOYq6+WgPHO96RulSR368BMhDQNNrUqZqysnTVoFI5srgS+L/AeUAekAu0i0icJRnRsWBhzDCI6JawTU3aOj0Re23019ioQWPp0lCa6rbbNE117rmpq2oKpqu6vTXCZWWhNutpqq5Khng3b0tlsNgAvA/tB7UI+CBwjteyI2UsWBgTpc5OTVE1NiY+RQVnpqlWrdJ3/bNnh3b6mzQpceeKJJiu6uzUe83P133Ig+mqNC0GjFciNm9LabAQkUXOuS0icpF3bJOIXBzryWNhwcKYGAUC2vW2vl5TOPn5GjgS+c47U9JUQT092ho+EMjIxYDDNdyFnUNJZelsp3MuD9js7cd9iOg2TzLGpFMwtz9lijb/27Mn8SmqkhJ4z3v0Y88eLcF9/nndEyMdaaqxY0PrRIJNDw8f1hFIaakGjokTEzuvk2CZ0r03mpHFDOAIOl/xeWAC8JiI1Cfv8s5mIwtjEig8RdXbqy/2Q6Soou5LBfq8r72mgWPlSk1TzZqlQeOOO1KbpgrX1RXaBTAvTwPp5MmJn9uJU6IWdqYyDVUE+ESkz3ucC4wTkc5YTx4LCxbGJEEgoFVU9fWDVlHF3ZcKNCUUTFNt3aqjmquu0jTVddelr5IpvHdVsCPw1Kk6+igsTM81kdieVakMFmuBm0Wk3XtcDCwTkauH/snEsmBhTBIFq6gaGzV4hC30e/PASb709Fb+5V0Xxt+XCjRNtXSpfjQ363luu00Dx/z56UsLiWjA9PnOXtMxfnxKJ8kTubAzlcFis4gsjHQs2SxYGJMiHR2wbx/s3UtXt5/vvXaEP+9r58raiXzh1nNiH13019sL69eH0lTd3ZqmWrwY7rwzfWmqoPA1HcFJ8mALkkRWlSVZKie4O5xzl4jI696JLwV8sZ7YGJPhiop0InrWLN5av5M36rcyMeDnjbe7aLy4OjGjC9AX4Cuv1I9gmmrpUu1L9cMfaprq7ru1miodFUx5eWdPkh86pI9LSkIdc4uLR/VK8mhGFpcBvwEOovtZTAEeEJGNybu8s9nIwpjUCubNX37zMKW+Nma0HuKdFWP461vOo7C8LHkvkME01fPPw5EjmgIKpqnOPz8zqpe6unQE1turKbspUzJ2JXmqNz8aCwT3Y9wpIj1hX7tFRJbHeiHDZcHCmNQaKG8+pr2NOV1HmXTymAaLZFYQDZSmqq0NpakqKob88ZgquGIRCIRWkjunK8i9/cgpKkp7cMuY3lDOuddF5JKEPNkQLFgYk0G6uuDgQe1FFQhELL2NW3s7LF+ugWPLFg1UV16J77Y7KLj5xrPSVAmp4IrFQCvJq6pC+3SkoTQ3k4JFSlZzW7AwJgMFV4fv2qXvruPdY2M4mpqgro6+pUvJaW5Gxo/H9UtTJbyCK1YZUJqbScHCRhbGZLtg6W1DgwaPsWM1RZXEid839x7nyR8/wxf9b1O69o+n01Q9d9zJD/Pn8soxEl/BFY/+pblFRTrqSHJpbiqroRLOOVcD/AqoBAR4XES+n85rMsbEwTmtDJo4Ud9J792r5begQSPBO+35/AGe3nKIzVXz+GHtVXzhq1+hcPVKWLKEsY/9iM84xzuqzuXlOZez5/wK5tcOPb+REs7pqCs48uru1nUt9fUaVCsrdaI8w0pzEzmyeFpE3hXlz0wFporI68658cBG4D4R2T7Yz9jIwpgRprtbS03r67WxXwLnNYJppqBgusnnD/DzJ1dRvvJlbmxYz+SOVnz5hYy543bG3ntP5lRT9dfbq/McXV066ggvzY2z6WMqF+UNFAhOAltFpHmAr0V/Mc79HvjhUFVVFiyMGaF6e0Ndb9va4t6YyecP8N1lO1nXeOL0sWC6yecPsO9El1Zw9fVRvG0z5SuWUbZ2Da67G2bODFVTTZ6cmPtLhmD/qr4+nbxftEgDSAxSGSyWAlcBK7xDN6AjgVrgmyLyn7FehPf8M4HVwAUicqrf1x4GHgaYPn36pU1NTfGcyhiTTiK6sK2hQdt8jBmjE71Rzmsca+8KBQSPc1BTlk958SAjl/Z2ePllraZ6443T1VQsXgzXX5/ZbctbWuCSS2IObqkMFi8BHxSRI97jSnS+4f3AahG5IOaL0D5Tq4BvicjTQ32vjSyMGUX6z2uUlg5ZVprQNRP79mlDw7q60KK/W2/V1eKZmKZKc7CIJpTXBAOFp9k7dhzoGeRnIvIW+j0F/FekQGGMGWWKi7Vp4A03wLx5mp46ckTTL/34/AHW7DpKpz+QmHPX1MAnP6mjjMceg2uv1cDxoQ/B/ffDL3+pIx8DRFcNtdI5V4duqwrwHu9YEdAay8mdcw74D2CHiHwvlucwxowC48bpPEJNjb5A796t/4at12ho6eCxlbuZMbEwsWsmcnLg8sv140tf0jTV0qXal+qxx+CKK0JpqgyqTkq1aNJQDngXcK136E/AUxJHOZVz7lpgDbAV6PMO/72IPD/Yz1gaypgs0G+9ho8cvrv2EOuaTqZuzcS+fRo06up0d73i4lCa6oILUp+mSnMaati/bRER59wfAT+6JuK1eAKF95x/RJsSGmNMSPh6jbY26v+8hfrtTZTl5PDa7l4aWxLY9XYwNTXwiU/Aww/Dxo2hLWKffhpmzNCgkenVVAk07DkL59x7gdfQ9NN7gXXOufck68KMMSNLa6c/Kc/bkVfAD47ksXbGRewpncr47k5eXr2VzraOpJzvLDk5cNll8M1vwksvwde/DuXlmqa66y749KfhxRcHnGcZTaJJQ70B3BJcU+GcqwBeFpEFSby+s1gaypjM09Ed4JlNB7jv4uqot/uM5Kyut4EAecdamHPqIBNdn7bLKCpK6DmHZf/+UJrq0CG9hmCa6sILE5+mGilpKCCn3+K7Y0RXTWWMGaW2HzzF157dxrzK8VxWOzGhz11ZUjDAFqJToO98OH5cF/k1N+skeUlJ6uYSpk2Dj38c/vqv4fXXNU31wgvwzDMwfXooTVVZmZrrSbJogsWL3lqLJ73HDwCDTkQbY7JDR3eAx9fsBuCnaxo4r6ok4aOLAQW7t06aBCdP6mZJBw/qOo0JE1K3V3ZOjq6sXrQIHnkEXnlFRxs/+pFWU11+uQaOG24Y0dVU0W5+9G7gGu/hGhF5JilXNQRLQxmTWdY3Huf+n7x6+vFvP35VwkcXw9bZqVVMTU1aUZWE5oXDFkxTLV2qQSyYplq8GC66KPoR0EhZwZ0pLFgYkzmCW64u3x7KUN86v5LvPbAwNaOLwfj9oeaFfn/yN2UaSl+fpqnq6nQNR1eXpqmCvammTBne82R6sHDOtaGlsmd9Ca2oja2rVYwsWBiTOQbactU5mDO5aIB5hjTo37wwFZsyDaWjA/7wB53feP11/WVddpmmqd75zqEDWqYHi0xjwcIYE7XgIr/du+HoUcjL0xRVOvs/DZSmuuUWHXEsWHD2tVmwiI4FC2NMXE6d0jmN/ft1Ery0NHWT4QPp64NNm0JpKp9PFwQuXqzrOIJpKgsW0bFgYYxJCJ9PA0Zjo75gT5igI4506uwMpak2bgylqRYv1rUbV19twWK4LFgYYxKqp0d7P9XX6+RzOifDwx04EEpTHTigm0U99RTccUdMTzei9+A2xpi0GztW0z5VVaHJ8H4db9Oiulr7Un3sY7B5sy72mz8/bZdjwcIYY0DnLaZM0RXXwcnw5ub0T4bn5OhcRU2Nji7SxIKFMcaEC+94m2mT4WlkwcIYYwZTUqITy7NnhybDIb0rw9PEgoUxxkRSWAjnnKO7+WXKyvAUs2BhjDHDlZenGx9Nm6bzGbt26b/papOeQhYsjDEmWrm5MHWqTogfO5a+NukpZMHCGGNi5dyZbdIbGzVNNXaszmvkjJ4tfyxYGGNMIkyYAAsXwty5sHevfoBWUI0Z+S+1aQ97zrnbnXM7nXP1zrkvp/t6jDEmLkVFcN55utnR3LlaftvcDN3d6b6yuKQ13DnncoEfAbcA+4H1zrnnRGR7Oq/LGGPiNm4czJqle1ccPqyT4SdP6qrwwsJ0X13U0j02uhyoF5EGAOfcb4B7AQsWxpjRYcwYrZ6qqtLJ8F274MiRUDuRETIZnu5gUQ3sC3u8H7giTddijDHJk5MDFRU6Gd7aCg0Nmp4aIZPh6Q4Ww+Kcexh4GGD69OlpvhpjjImDc1BWBpdeqrv37d2r+4Y7l9GT4ekOZQeAmrDH07xjZxCRx0VkkYgsqqioSNnFGWNMUo0fD+efr5Phs2bpnEZLi64OzzDpDmHrgbnOuVo0SLwP+Mv0XpIxxqRYfr5WTs2cGdpbo7U1o9qJpDVYiEjAOfdp4CUgF/i5iLyZzmsyxpi0GWpvjTRL98gCEXkeeD7d12GMMRljoL01Tp5M6yWlPVgYY4wZRPjeGm1taU1JWbAwxpiRYPz4tJ4+3dVQxhhjRgALFsYYYyKyYGGMMSYiCxbGGGMismBhjDEmIgsWxhhjIrJgYYwxJiILFsYYYyKyYGGMMSYiCxbGZLHWzsxrhW0ykwULY7JUR3eAui2HaO8OpPtSzAhgwcKYLLX94Cm+9uw2dhw8le5LMSOABQtjslBHd4DH1+wG4KdrGmx0YSKyYGFMFtp+8BTLtzcDsGz7ERtdmIgsWBiTZcJHFUE2ujCR2H4WxmSZ9u4ePnR1LQ9dVXv6mHPQ0d1D8Th7STADs78MY7JMZUkBlSUF6b4MM8JYGsoYY0xEFiyMMcZElLZg4Zx71Dn3lnNui3PuGedcabquxRhjzNDSObJYDlwgIhcBbwNfSeO1GGOMGULagoWILBORYK3eWmBauq7FGGPM0DJlzuIjwAuDfdE597BzboNzbkNLS0sKL8sYkw7W4DDzJDVYOOdeds5tG+Dj3rDv+SoQAP5rsOcRkcdFZJGILKqoqEjmJRtj0swaHGampK6zEJGbh/q6c+5DwGLgJhGRZF6LMWZkCDY4nFc5nstqJ6b7cownndVQtwOPAPeISGe6rsMYkzmswWHmSuecxQ+B8cBy59xm59y/p/FajDEZwBocZq60tfsQkTnpOrcxJvMM1uDwvKoS61mVAey/gDEmI1iDw8xm/wWMMRnBGhxmtkxZZ2GMMSaDWbAwxhgTkQULY4wxEVmwMMYYE5EFC2OMMRG5kdZlwznXAjTF+OOTgKMJvJyRJJvvHbL7/rP53iG77z/83meISMzN9UZcsIiHc26DiCxK93WkQzbfO2T3/WfzvUN2338i793SUMYYYyKyYGGMMSaibAsWj6f7AtIom+8dsvv+s/neIbvvP2H3nlVzFsYYY2KTbSMLY4wxMbBgYYwxJqKsCRbOududczudc/XOuS+n+3oSwTn3c+dcs3NuW9ixic655c65Xd6/Zd5x55z7gXf/W5xzl4T9zEPe9+9yzj2UjnuJlnOuxjm3wjm33Tn3pnPus97xbLn/fOfca865N7z7/4Z3vNY5t867z/92zuV5x8d5j+u9r88Me66veMd3OuduS9MtRc05l+uc2+Scq/MeZ8W9O+f2OOe2epvGbfCOJf/vXkRG/QeQC+wGZgF5wBvA/HRfVwLu6zrgEmBb2LHvAF/2Pv8y8C/e53cCLwAOuBJY5x2fCDR4/5Z5n5el+96Gce9TgUu8z8cDbwPzs+j+HVDsfT4WWOfd1/8A7/OO/zvwSe/zTwH/7n3+PuC/vc/ne/8/jANqvf9PctN9f8P8HXwB+H9Anfc4K+4d2ANM6ncs6X/32TKyuByoF5EGEfEDvwHuTfM1xU1EVgPH+x2+F3jC+/wJ4L6w478StRYodc5NBW4DlovIcRE5ASwHbk/6xcdJRA6JyOve523ADqCa7Ll/EZF27+FY70OAG4Hfecf733/w9/I74CbnnPOO/0ZEukWkEahH/3/JaM65acBdwM+8x44sufdBJP3vPluCRTWwL+zxfu/YaFQpIoe8zw8Dld7ng/0ORvzvxksrXIy+u86a+/fSMJuBZvR/9t1Aq4gEvG8Jv5fT9+l9/SRQzsi9//8DPAL0eY/LyZ57F2CZc26jc+5h71jS/+5tp7xRTETEOTeqa6Odc8XAU8DnROSUvmFUo/3+RaQXWOicKwWeAc5N7xWlhnNuMdAsIhudczek+XLS4VoROeCcmwwsd869Ff7FZP3dZ8vI4gBQE/Z4mndsNDriDTPx/m32jg/2Oxixvxvn3Fg0UPyXiDztHc6a+w8SkVZgBXAVmmYIvgkMv5fT9+l9fQJwjJF5/9cA9zjn9qAp5RuB75Md946IHPD+bUbfJFxOCv7usyVYrAfmetUSeegk13NpvqZkeQ4IVjY8BPw+7PgHveqIK4GT3rD1JeBW51yZV0Fxq3cso3k55/8AdojI98K+lC33X+GNKHDOFQC3oPM2K4D3eN/W//6Dv5f3AH8Qnel8DnifVzFUC8wFXkvJTcRIRL4iItNEZCb6//IfRORBsuDenXNFzrnxwc/Rv9dtpOLvPt0z+6n6QKsC3kbzul9N9/Uk6J6eBA4BPWjO8aNoLvYVYBfwMjDR+14H/Mi7/63AorDn+Qg6uVcPfDjd9zXMe78Wzd1uATZ7H3dm0f1fBGzy7n8b8HXv+Cz0Ba8e+C0wzjue7z2u974+K+y5vur9XnYCd6T73qL8PdxAqBpq1N+7d49veB9vBl/LUvF3b+0+jDHGRJQtaShjjDFxsGBhjDEmIgsWxhhjIrJgYYwxJiILFsYYYyKyYGHMILzunpOcc6XOuU+FHa9yzv3O+3yhc+7OGJ77H51zX0zk9RqTTBYsjImsFO1cCoCIHBSR4OKvhej6DmNGNQsWxgDOuWe9xmxvhjVnC/o2MNvbP+BR59xM59w2rxvAN4EHvK890H/E4H3fTO/zrzrn3nbO/RGYF/Y9s51zL3rnX+Ocy4oeT2ZksUaCxqiPiMhxr3XGeufcU2Ff+zJwgYgshNNdbhERv3Pu6+iq2E97X/vHgZ7cOXcp2ppiIfr/3evARu/LjwOfEJFdzrkrgMfQfkfGZAwLFsaozzjn/sL7vAbtE5RI7wCeEZFOAOfcc96/xcDVwG/DOuaOS/C5jYmbBQuT9bw21zcDV4lIp3NuJdpPKBYBzkzvRnqeHHQfhoUxns+YlLA5C2O0ZfUJL1Cci24/Ga4N3bp1IP2/tgfd6hZvv+Na7/hq4D7nXIHXNfRuABE5BTQ65+73fsY55xbEf0vGJJYFC2PgRWCMc24HOpm9NvyLInIM+JM3Wf1ov59dAcwPTnCj+2tMdM69CXwa7XSM6Baw/412C30BbZsf9CDwUedcsJPoiN/y14w+1nXWGGNMRDayMMYYE5EFC2OMMRFZsDDGGBORBQtjjDERWbAwxhgTkQULY4wxEVmwMMYYE9H/B90MrAx0yngYAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"original_data['log_cases_pop_den'] = np.log(original_data['cases_pop_den'])\n",
"\n",
"ax = sns.scatterplot(\n",
" data=original_data,\n",
" x='altitude',\n",
" y='log_cases_pop_den',\n",
" marker='^'\n",
")\n",
"\n",
"prediction = original_model.get_prediction(original_data['altitude']).summary_frame()\n",
"\n",
"ax.plot(\n",
" original_data['altitude'], \n",
" prediction['mean'],\n",
" color='red'\n",
")\n",
"ax.fill_between(\n",
" original_data['altitude'], \n",
" prediction['mean_ci_lower'],\n",
" prediction['mean_ci_upper'],\n",
" color='red',\n",
" alpha=.2\n",
")"
]
},
{
"cell_type": "markdown",
"id": "suited-falls",
"metadata": {},
"source": [
"### Ponderado"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "linear-insert",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>WLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>np.log(cases_pop_den)</td> <th> R-squared: </th> <td> 0.135</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>WLS</td> <th> Adj. R-squared: </th> <td> 0.115</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 6.568</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Fri, 31 Dec 2021</td> <th> Prob (F-statistic):</th> <td>0.0141</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>10:53:03</td> <th> Log-Likelihood: </th> <td> -106.87</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 44</td> <th> AIC: </th> <td> 217.7</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 42</td> <th> BIC: </th> <td> 221.3</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 1</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>t</th> <th>P>|t|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 2.7469</td> <td> 0.234</td> <td> 11.763</td> <td> 0.000</td> <td> 2.276</td> <td> 3.218</td>\n",
"</tr>\n",
"<tr>\n",
" <th>altitude</th> <td> -0.0008</td> <td> 0.000</td> <td> -2.563</td> <td> 0.014</td> <td> -0.001</td> <td> -0.000</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>13.859</td> <th> Durbin-Watson: </th> <td> 2.099</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.001</td> <th> Jarque-Bera (JB): </th> <td> 31.003</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td> 0.654</td> <th> Prob(JB): </th> <td>1.85e-07</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td> 6.899</td> <th> Cond. No. </th> <td>1.02e+03</td>\n",
"</tr>\n",
"</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 1.02e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" WLS Regression Results \n",
"=================================================================================\n",
"Dep. Variable: np.log(cases_pop_den) R-squared: 0.135\n",
"Model: WLS Adj. R-squared: 0.115\n",
"Method: Least Squares F-statistic: 6.568\n",
"Date: Fri, 31 Dec 2021 Prob (F-statistic): 0.0141\n",
"Time: 10:53:03 Log-Likelihood: -106.87\n",
"No. Observations: 44 AIC: 217.7\n",
"Df Residuals: 42 BIC: 221.3\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 2.7469 0.234 11.763 0.000 2.276 3.218\n",
"altitude -0.0008 0.000 -2.563 0.014 -0.001 -0.000\n",
"==============================================================================\n",
"Omnibus: 13.859 Durbin-Watson: 2.099\n",
"Prob(Omnibus): 0.001 Jarque-Bera (JB): 31.003\n",
"Skew: 0.654 Prob(JB): 1.85e-07\n",
"Kurtosis: 6.899 Cond. No. 1.02e+03\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 1.02e+03. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n",
"\"\"\""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"weighted_data = get_binned(bin_width=100, stat='mean')\n",
"weighted_model = smf.wls(\n",
" 'np.log(cases_pop_den) ~ altitude', \n",
" data=weighted_data, \n",
" weights=weighted_data['count']\n",
").fit()\n",
"weighted_model.summary()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "provincial-function",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PolyCollection at 0x7f45ef86a0a0>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1JklEQVR4nO29eZhc5XXn/zndrV61tPa90QJCCCEJ1EISEpJMsCEGY0LsYCeODSRRPLGTOJ4ZJ44zTn4z41/seJ4k9ng8tryC7WCMDTbBYCMTuiWMWSQkNgntG5JQNxJSq1u995k/zr2qUqu3qq6qW8v5PM99qurWct+3uvp+73nPJqqK4ziO4wxEUdQDcBzHcbIfFwvHcRxnUFwsHMdxnEFxsXAcx3EGxcXCcRzHGZSSqAeQKBMmTNBZs2ZFPQzHcZycYuvWrW+p6sRk359zYjFr1iy2bNkS9TAcx3FyChE5NJz3+zKU4ziOMyguFo7jOM6guFg4juM4g+Ji4TiO4wyKi4XjOI4zKC4WGeT0uY6oh+A4jpMULhYZoqW9i0dfPk5ze1fUQ3Ecx0kYF4tBSJU1sONYE3/301fZeawpJZ/nOI6TSVwsBiBV1kBLexcbNu8D4Bub97t14ThOzuFiMQCpsgZ2HGti444GAJ7YccKtC8dxco6sEAsRKRaRbSLyaNRjCUmVNRD/OSFuXTiOk2tkS22ovwR2AqOjHkhIX9bAstnjEv6c5vZO7rpuNh9ZOfv8PhFoae9kZFm2fP2O4zgDE/nZSkRmALcAnwM+GfFwgP6tgSumjU74BD95dAWTR1ekcniO4zgZJ3KxAP4V+BQwqr8XiMh6YD1ATU1N2gfk1oDjOM6FRHrmE5FbgQZV3Soi6/p7napuADYA1NbWarrH5daA4zjOhUTt4F4F3CYiB4EfAjeIyPejHZLjOJnAKxrkFpGKhap+WlVnqOos4APAf6jqh6Ick+M46ccrGuQeUVsWjpNz+BXx8PGKBrlH1oiFqtap6q1Rj8NxBsKviIePVzTITbJGLBwnF/Ar4uHjFQ1yExcLxxkifkU8fLyiQe7iSQOOM0RSldVfyHgOU+7ifx3HGQKpzOovZDyHKXfxX7njDAG/InYKHf+VO84Q8Ctip9BxB7fjOI4zKC4WjuM4zqC4WDiO4ziD4mLh5BxebsNxMo+LhZNTeLkNJ1X4RUdiuFg4OYWX23BSgV90JI6LhZMzeLkNJ1X4RUfiFJZY7N0LjY3Q0xP1SJwk8AJ0Tirwi47kKCyxOHIEnn8eNm+G48ehuzvqETlDxAvQOalisIsO92X0TdQ9uMuBTUBZMJYfq+rfp/WgEyaYSGzfDmVlcOmlMHUqjBiR1sM6w8PLbTipYLAaX6Ev4/arp/vvqhdRfxvtwA2q2iwiI4CnReRxVX02rUctK4NJk6CzE3bsgNdfh7lzYfp0KC9P66Gd5PByG04qGOyiI/RlXD55VNZVFD59roPqytLIjh+pWKiqAs3BwxHBphkbwIgRMHEidHXBvn2wZw/MmAE1NTB6dMaG4WQnUf9zOqlnoIuO3r6MbKoonA0WT+Q+CxEpFpHtQAOwUVWfy/ggSkpg/HjbTpyAX//afBsnT4JmTruc7MFDKwuPbA6gyIborcjFQlW7VXUJMAO4VkQW9n6NiKwXkS0isqWxsTF9gykqgupqW6JqbTXB2LQJ3ngDOtzpVUhkwz+nkzmyOYAiW6K3ssPGAlT1tIg8BdwMvNrruQ3ABoDa2trMXOqPHGlbezu89pptM2faNmpURobgREM2L0c46SGbAyiypUNj1NFQE4HOQCgqgHcCX4hyTBdRVmZbdzccOwaHDsHYsTBnji1bFRdHPUInxWTLP2e+kAu+n2wNoMimDo1RXy5NBe4VkWJsSexHqvpoxGPqm+JiEwmAc+dg61bzdcycCdOmmbUhEu0Y84woTjLZ9M+ZD2SDYzaXySaLJ+poqJeBqzNysDNn+s3cPtvWyajyBPIsKitt6+42f8aBAyYWs2ZZdFVZWWrGXMBEdZLJpn/OfCCbQ1FzgWyyeArn1//nfw6PPQZr18INN8DSpTBiBK0dXWze8xbrLp9IZWmCX0e8tdHWZn4NVROMmTPtOU/2S4qoTjLD+efMheWWTOK+n/yicP5yt99u/obHH4eHHoKqKli1ircWXct3jo7kknHLuHL6mOQ/v7zcNlVobrZlKhGLrJo506KsXDiGRC6eZHy55WLc95NfFM6v+o477Eq/uBi2bYO6OrS+nplPPMG/FRVz+PkFdNzx25Te8A6zDJJFJBZJpQpnz8KWLRaWO2mS+TfGjPFM8QEYzkkmqqt7X265EPf95B+F91crK4PVq2H1anZ8+M/47jceZcXhV1h5+BVKv/hP8MV/goULbblq3TrzQyTruI4Xjp4eaGqypD9V83FMnw7jxtn9oshTXrKC4Zxkorq6z0VLKN247yf/EM2xDOXa2lrdsmVLcm9+6ilzTJeU0NrRxf96YhfPHXjbnlPl1pFt3NN9iNLNm6xmFFjpj3XrbFu4MHUn9bY2aGkxJ3lJiVkzkyZZmZGqqoKNrDrR1MrehpYLEudF4NJJVYP6El44cIr3f/03PPinKzN6dR8eNyTTx3ecoSAiW1W1Ntn3F6zEn+vo4j2Lp3Prounn94nA2bHXM/5P/hgaGqC+3rYf/ADuu8/yKtasMatj2bLhRT2FPg4wwThzxqyOnp4LxWPUKBO4PM7niF86StbBHNXVfTYtt7iD3UknBWtZJERzs9WLqquDZ54xi6CyElauNItj1arUFh7s7rZyI62ttmRVVGR+jvHjzVFeVQUVFXlhfbS0d/HwtqPDXjqK6up+OJZQKknV9+jkL25ZZIKRI+Gmm2zr6DCHdWh1PPmkXfUvXWoWx9q1MGXK8I5XXBzzdYAJRnu7RXPt3Wv7SkpMQKqr7bay0gQkUSGMmFQ4hqO8us+WOHh3sDvpxi2L4dDTY76NujrbDh60/fPnx/wcc+emxwLo6TEBaW+3vhwhFRVm5YQWSLjcVZp9yxMt7V184oFtbNzRwLsWTOaf71yS1Mk9W67uoyJV36OT37hlESVFReb0XrgQPv5xE4u6OqtU+7Wv2TZ9uonG2rWweHHqfA9FRSYMFb1Ohp2dFnX11lsXto0tKTERGTXKtvLyWN2riIQkVXH42XJ1HxXpyGdw/4fTGxeLVDJrFtx1l21vvWWiUVcHP/qROcmrq+H66008li9PT67FiBF9J/91d9sS2vHjcPjwhX06iorMChk5MnZbWhrbRoxIuXWUTY7hXCYd36MnGDp94ctQmaClxRzj9fXw9NPmMC8vhxUrTDhWrzYhiYqeHrNIwi20SMLfhohZMJWVJibh8lYoTEkISqEvHaWKdHyPUYUgO+nFl6ESobgYTp2KneDKyjIjHFVV8M532tbZaaVAQgd5XZ2Na8mS2HLVtGnpH1M8RUWxJam+ULXWs62ttsTV2XlxB0ER+07DSK1wC4UkfhMp+KWjVJHq79ETDJ3+KCzLornZrvLPnrWT3pkz5iBWtZNdeNIsL8+MiKiagzwUjf37bf+8ebEM8nnzciNEVtUskq6uC7fwObB5qJqAlJdfKCplZSYkJSUX3ubC3PMITzDMX4ZrWRSWWPRFV5cJRmuricnbb8Pp07YvJDy5pbv0+JEjMeF46SU7sU6dGrM4lizJudDYPgmFJF5c4svHh6IiYoIRRnSFWxgiHC8s4eYkTXxUVYhHV+UPLhbporPTBOTcOROQkydNTABEONsjjBo3Jn0nqFOnzEFeXw/PPWfO6TFjzL+xbp35O3pHQuUj8aIS3obC0tdSWLicFgpLeL+3qISbWy7ncT9SfpPTYiEiM4H7gMmAAhtU9UsDvSdjYtEXXV3Q0sK5U2d48pmd3DheqCA4cZWWmgM4HWXIz52DZ5814di82ZbQysrg2mtNONasifXVKGTil8Lib3t6LhYWiFkuocDEC01pad/iUlzsAuPkJLnu4O4C/rOqvigio4CtIrJRVXdEPK6+CbKmXzvVzZ+/0sWD61ewbGql+UBOnrTaTm+/HbvCrapKTV5FZaU1bLrhBjsBbt8eW67avNl8LYsXx/wcM2YM/5i5iEjiy1GhoITLkKHI9GW1hITBEWGOSn8CU1zsFoyTN2TVMpSI/Az4iqpu7O81kVoWDCFbtq3NrvzffNPEo6vLThphuGkqUYXdu0006uvtPljWeJhBPn9+Rk9UCbeozTVC66X3FgpM6G+JJ96CibdkwmWyeFGJF5o8Lh7pZJ6cXoaKR0RmAZuAhara1Ou59cB6gJqamqWHDh3K/AADEooW6ekxq+Ptt+HoURMREROOysrUD+7o0Vgi4LZtdvzJk22Zat06q1+VRidwa0cXT+1qTK5Fbb7Tl7j09PS/RAZmMYaJkaGwhJZMaMX0FpriYu+N4vRJxsRCRO4AvgBMAiTYVFWHXW5VREYC9cDnVPWhgV4bpWUx7GiRtjZzXB85YrehcKSjguzp05YAWF9vCYHt7VbmY9UqE46VK+3YSdKXBfHa0TP89UOv8IU7rhpei1rHCMUk3vfS3zJZPCUlMZGJd/CHt/EC40tlBUMmxWIv8B5V3Znswfr53BHAo8AvVfWfB3t9lGKR0miR3sJRVGRlNtJRAqSt7UIH+enTdtJYtizmIJ8wYcgf15cFEd9MasXscXzyXfPcuoiKUFT62nrnvMTTW2RCC6a3yPQWG7dkcoJMisWvVXVVsgfq5zMFuBc4paqfGMp7ovZZpIW2NmhstBLkZ8/Giv6lY8mouxtefjlWKffoUTtx9G4lOwB9WRDhvhC3LnKQwUQmtDx6+2aKiy9cHhtouSz+1n0yGSWTYvElYArwU+B8xtpgy0aDfOZqYDPwCoQxqPytqj7W33vyUixCVE0sjh83i6OzM1aPKR1LBKqwb1/MQb4zMBpnzYJ162hZuZqqqxddcOXYlwUhcGGLWnDropDoLTJDXS4rKrrQ4d/bqulruSwUGV8yS5hMisV3+titqnpPsgdPhrwWi3i6u80xfuiQtXgtLrakvMDaSEvU0Ztvnq9ZpVu3It3d9EyYQFHoIK+t5bXG1ossiCljyjjydttFy3Mzx5YzfmR5+sbr5Da9I8vi/TOq/QtNGF0WXxk5XnD688sU+JJZ3kRDDZWCEYt4zp2DY8esX0ZnJ62l5Tx1pCWtUUc7d73Bz7/2Y/6k+zBjtr0Ara1oZSU7Zi/k5+MuZ8uMBZwrrRiSBeFRUk7K6WupLLRoBjqnhUtm8ctkQ10yy3FrJmNJeSIyD/i/wGRVXSgii4DbVPV/JntwZ4hUVsKll8Ls2XDyJAc2b+eBn29hTsli5l8+I+VXS60dXfx4zxmem1NL2+x38cnP/Q8qX9pG+8Ynmff0Zq587Xl6SkpoXriYM+WraD86isrZM/v9vP2NLXy1bh+XjKvMOj+GWzw5SrI+j/gIs+bmC5fQ4uuT9aY/aya87SsAILzNE2smkcu8bwD/Ffg6gKq+LCL/BrhYZIriYlrGjONfWsaydcaV3H+4g78e3UBl+YgLlqiGy/7GlvP+h2cPnOLA1dO5cvVqylevtn+qV16hqL6e0fX1jP76l+HrX4YFC2KJgLNnn78Ka+3o4qFtbwDw8LajzJ5YFYl10ZcotHZ0sXnPW27xFBJFRbYl+r8SLpn19FidttbWoVsz8b6ZeHGJz/rvb9ksi6yZRL6xSlV9Xi4cfFeKx5PzpLsd5fkWmmWV3NsE75l3NbWlbeao7u6OtUxNkviTe8gFJ/mw98aSJfAXfxFrJVtXB1/9qm01NRZZtXYt+8fVXCg8jS0Zty76E4VstnicLCMsJZMMoTUTlpVpaYklZXZ39x3GHB6zuPjCzP/LLrMQ+whIZPZvichcrOAfIvI+4HhaRpWjpLsdZV8tNDc8f4z5dy5hZE2NOcL37rXbsKNdgpzr6OI9i6dz66Lp5/eJ2An3oqtvEbMiZs+Gu++244aVcu+/H773PS6pGs2fT1vAb2qu4qWp8yKxLvoShWyxeDKBL7VFTLLWDFxovRw9CtOn54RYfAzYAMwXkaPAAeBDaRlVjrLjWBN/99NXuXzyqLQ0jGlu7+Su62bzkZWzz+8TgZb2TkaOrrAOe1OnWlHDUDQqKszaGCLjR8YimBJm0iR43/tsa27m7H/U0/Ufddz44nPctOdZusvLOXv1Mro71sGN77BckjTTnyhctNQWgcWTCXypLccJfTNhif0IGfLRVXU/cKOIVAFFqno2fcPKPTLRjnJILTRFLBt7/HjL1N6/30SjtNT8GplaAx05klG33QK33XK+lWxxXR3V9fXwm83whf8frrkmlgg4ZUpahtGXKMyZWDXwUlse4UttTqoYNHRWRD450PNDKdGRSrI1dDar21GePWu+hTfesCuUMWOii9Do6Ym1kq2vj7WSnT8/5iCfOzclohafQBiyYvY4/uT6ORxvGjgvJB/wEix5RmOjXWBNmpTU2zMROhuuYVwOLAMeCR6/B3g+2QPnE335ErKq2f2oUXDVVTBnjiX5HTpkJm11deZFo6jISossXAgf+5iNJezN8fWvw9e+ZuuyocWxaFHS5nd//peSYlgyszol08lm0rHU5v6PwiWRDO5NwC3h8lPQrOjnqromjeO7iGy0LHKhHeUFUVqtrXaSPnjQBjp2bHbU6Tl50god1tXB88/HWsmuWWPisWJFegot5iH9WVXDsS48uTJicsCyCJkMdMQ97gj2FTxD8iVEyEVRWhUVtuwzaxYcOcLZnbsYVV5qlkaUojF+PNx+u20tLfCb38Ssjn//dwsdXLHCLI7rr7fxOn2SUFTbEHH/R2GTyK/mPuB5EXk4eHw78N1UD8hJPf1GaZWX01Izm38/rrx3fA9VbxyyM0rUogEW9nvjjbZ1dcGLL8YKHtbX23LWkiXn8zkKtpVsPwwrqq0PCinU2OmbIS9Yq+rngLuBt4PtblX9x/B5ERmb+uE5w6V3lFZz+4V5lDuONfG3j+1hR+WkWPb16dO2JNTdnfkB90VJCVx7LXzqU/Doo/D978M991jnwX/5F7NEPvhB83m8/vrA2bRZxtm2zqiHMCT68n84hUVClwaq+iLwYj9PPwlcM+wROSnlfMY38MSOE+w81nTeurgo3PfOJYy87DLLwD58GA4csA/JBksjRMSW0ObPh49+1CK8QovjW9+Cb3zDwnBDB/nVV0cen94fuZIDMWhWv1MQpPIvnT1FTBxg8CitfoUkLCsQikYY3potjvB4ZsyAD33ItrffjjnIf/pTeOABiwRbvTrWSjYdvc+TJFd8AOnwfzi5Ryr/0knZ/iLybeBWoEFVF6ZwPAXPQBnfAoOH+4aiMXOmRU8dOGC+grFjs7OS5tixcNtttrW2XthK9vHHLTFx+XKzOtasgXHR5cHkkg8g1f4PJzfJhl/nd4GvYA50J4UMFKV1oqm1/9IhvXNDysvh8svN0ghDbktKok3uG4yKCnjHO2zr6oLt22ORVZs322QXLYolAs7sv8R6OhgsB8LzGZxsI/JlKFXdJCKzUjiOfkl3RdhcIqlw3zDk9pJLzMo4dMgywqurs6qU8kWUlEBtrW2f/CTs2ROrlPulL9k2Z46Jxtq1Vm49xfOJP/kP5gPIFV+GU1gk9EsUkWuA1diS068Dh3fIb6VyYKkm3RVhC4qKCjuhXnKJ+TPeeCPztaeSRQTmzbNt/XrrQLhpkwnHvffCt79tSU9hSO7SpSaIw6D3yX8wH0Cu+DKcwiKRTnmfBd4PPBTs+o6IPBh2ylPVU2kYX3js9cB6gJqamqQ+I90VYYdDzlo8VVVWRmTWLOunceyYCUkGqsmmjGnT4AMfsO3MGXj66VgS4IMPWjnoVatiDvIkykP3PvkP5APIJV+GU1gk8iv8A2CxqrYBiMjnge1koFOeqm7AyqNTW1ubsCM9ExVhkyXVFk8kwjNqlCXIzZkDu3fDiRMmJBHV3U+aMWPglltsa2uDF14w4di0CX75S7Mwli2LWR0TJgz6kYme/AuldLqTeyTinTwGxF8OlQFHUzuc9NBXiGi2EFo8qRhTKDy9E+8yxujR5hdYudIiqU6cgHPnohnLcCkvt5Ii/+2/wS9+Ad/8Jtx5p4US/+M/ws03W8On737XHP79kEgyW3++jHMd3pDSiZ5ELmXPAK+JyEbMZ/FOrPzHlwFU9S+SGYCI3A+sAyaIyBvA36vqt5L5rL7I5oqwqbZ4smapbexYC1E9dcoyqhsaTEhytQhgfCvZv/xL89OEiYBf+YptNTWxyKqFC6GoKOFkNs9ncLKZRKrOfmSg51X13pSMaBASrTqbzRVhU9kDo6W9i088sI2NOxp414LJ/POdSyIXQ8BKbzQ0wK5dVhxwzBizOvKFEydiDvItW6xEyvjxsGYNTcuv48DsK+kZEVsWzMe+GSEe7ptmcqXqrKreKyKlwLxg1y5VzfrCNtlaETbVFs9AZT0iRQQmT4aJE+3E+vrr5kgeO3bYUUZZweTJ8P7323b2LPz61yYcv/wlox9+mMWVlXDddWZxrFqVUIvbXMLDffOfRKKh1gH3AgexnIqZIvIRVd2UlpHlOQP2005QLLJ5qe08RUXWH3ziRIua2r3brsKrq7O2dlPCjBplvoybb7ZeHPEO8l/9ypazli6N5XNMzp8K/x7um/8ksgy1Ffh9Vd0VPJ4H3K+qS9M4vovIxuZHUZPNS2390tkJR47A3r22VDVuXPZmgw+Xnh549VXzcTz1lDnJwXJVwsiqFLWSjQJv35ohIl6GSkQsXlbVRYPtSzcuFnlGe7tlgx84kDuJfcPl4MFYBvmrr9q+GTNiDvKrrsq+go0D8NrRM/z1Q6+cf/yFO65y6yId5IrPAtgiIt8Evh88/gPAz9rO8CgrsxIiNTUWZXTkiO0bk8cnm1mz4K67bGtstGWq+nr44Q+tV8fYsRa2u26d9fHI4igyL19eOCRiWZQBH8PKfQBsBr6qqu1pGlufuGWR55w9a7Wb3nwzNxP7hkNzMzzzjAnH009b9Fh5ueWtrF1rAjJmTFZFHZ1sbuPI220XLYHma8RXpOTKMlRwsFLgCqAHi4bqGOQtKcfFokA4fdoip06dshyNiiz1vaSLzk7YujWWz9HYCMXFdC9ewuvzr2HO776biksyWynXiZhcEQsRuQX4GrAPi4aaDfypqj6e7MGTwcWigFC19q47d9pVd77laAyVnh77DurqaNv4JOVvBA7yyy+PdQS87LL89/UUOjkkFq8Dt6rq3uDxXODnqjo/2YMnQyGJRc4WGEw1PT2xHI22Noucypdw2wQIo46OvLSb32vexztO7KD4lVdMVKdNiwnH4sUF+f3kPTnk4D4bCkXAfuBssgd2BsZLqscRn6Nx9KjlaPT0ZGeb1zRyvs7U6In86+iJTP34n3BleZc1c6qvh5/8BO6/3yyw0EG+YkVWO8gTJZv8NYVGotFQjwE/wmpDvR94QUTuAFDVhwZ6s5MYWVPnKZsoKbEeGlOnWq7C3r0mFtXV+ZujEdBv1NG75lF5++1w++1WtPHZZ2OJgI8+ast2y5ebcFx/vQlsjuJZ4tGSyDdeDpwA1gaPG4EK4D2YeLhYpIhsLqmeFZSWwqWXwvTpsY59eZ6jMaQig5WVcMMNtnV1wbZtsUTATZtMUBcvji1XzZgRzWSSxLPEoyWhaKgBP0jk06r6jyn5sAEoBJ9FKgsMFgTNzRZue/x44YXbDgVVK+QY9iDfs8f2z50bSwScPz+rhXY4WeJ5s3QVsc8ilbb7+1P4WQVLf3WeIutRkQuMHAlXX20F+8rLzRne2hr1qLIHERODP/1T82n87GfWi3zMGPjOd+AP/xBuvRW+8AV47jmzSrKMwfqCnG3ru6ZpuHTlPUGGTyrXNrL3siSHSGWBwYKjutoynk+ehB07rDR6dbUtUTkxpk+H3/99206fjmWQP/KItZIdNQpWr7blqpUrzVqLkMGyxAfyZfjSVepI5dknNetZBU62llTPGUSs3enq1WZh7NwZK4nu4aQXU10Nt91mW1tbzEG+eTM8/riVkb/22piDfAitZFPNYP6a/gTB+5mnlsgtCxG5GfgSUAx8U1U/n8IxOYVKX+G2qnZyLKBw24QoL4/5MLq64KWXzOKor4fPfc7O0FddFSuxfsklGRnW+JH9lw4ZSBC8n3lqSaVYPJjoG0SkGPg/WIvWN7BQ3EdUdUcKx+UUMmG47ZQpVu31wAG7Ws7jyKmUUFJivTeWLoW/+isLUw4d5F/+sm2zZ8eEY8GCSMKX+xMEL3CYehLJ4P4n4H8CrcAvgEXAX6nq9wd848CfuRL4B1W9KXj8aYCBoqoKIRrKSSPnzsG+fVbdtqLC6k45ifHmm7Fcjq1brYnVxImxkNylSzPSBTE+QiokjJRq7ejKvwKHOVTuY7uqLhGR3wFuBT4JbFLVxUkfXOR9wM2q+sfB4z8Elqvqx3u9bj2wHqCmpmbpoUOHkj2k4xhNTbY01dBgDt3KyqhHlJucORNrJfub31gUWlWVtZBdt84i1NIUylxwFW9zqNxH+NpbgAdV9YxkyIxX1Q3ABjDLIiMHdfKb0aOhttaq2u7cac7w6urCLFQ4HMaMgXe/27a2NmslW19vVscTT9hy1rJlsY6AEyem7NAD+TKc1JOIWDwaFBNsBf6TiEwE2oZ5/KNAfJ3lGcE+x8kM48ZZeGhDQyxyqkALFQ6b8nKLmLr+eluaeuWVmJ/j85+37corY070WbPcb5RDJNrPYhxwRlW7RaQKGKWqbyZ9cJESYDfwW5hIvID1+X6tv/e4z8JJG93d8MYbHjmValQtsCDszfFa8O9dUxNzkF91Vd7X9xo2ubIMJSKVwJ8BNZj/YBpwOfBosgdX1S4R+TjwSyx09tsDCYXjpJXi4lihwkOHzBFeUmKi4VfAySMCc+bYds89ZsVt2mTi8YMfwH33wfjxsGaNCceyZb4cmIUk4uB+ANgKfFhVFwbi8YyqLknj+C7CLQsnY7S2xiKnyss9ciodNDfHHOTPPGOtZCsqzDG+dq0lV/r3buSKZQHMVdU7ReSDAKp6TjLl4XacKKiogIULbblk925zgnvkVGoZORJuusm2jg7YsiUWlvvkk2btLV0ac5BPmRL1iAuWRCyLZzDfwq9V9ZqgU979qnptOgfYG7csnMgII6fOnPHIqXTT02P1verqbDt40PbPnx9zkM+dW1jLgzmUZ/FO4O+ABcATwCrgLlWtS/bgyeBi4URKT08scqqAW7xmnIMHYxbHyy/bvunTL2wlm+/BCLkiFsHBxgMrsDpQz6rqW8keOFlcLJyswCOnouOtt2IO8hdegM5O+/7DVrLLl+dVK9nz5IpYiMgqYLuqtojIh4BrgC+pakbTqV0snKyio8Ouevfv95pTUdDSYpnjTz1ljvLmZhOKFStMOFavNiHJB3LIwf1/gcUishgr9fEt4D5ibVYdp/AoLYV586xF6b59Zm1UVJgj3Ek/VVVw4422dXbCiy/G8jnq6szaW7Ikls8xbVq0481hErEsXgwc258Fjqrqt8J96R3ihbhl4WQ1TU3WwrSx0SOnokTV/Eqhg3z/fts/b17MzzFvXm5ZgTm0DFWPVZu9G1gDNAAvqepVyR48GVwsnKxHNRY5dfasd+vLBo4ciVkbL71kf6OpU2PCsWRJ9gcq5JBYTAF+H3hBVTeLSA2wTlXvS/bgyeBi4eQMPT2xbn0dHd6tL1s4dSrWSva55+xvM2bMha1kK7KwW2WuiEW24GLh5BxdXbHIKTDR8DpI2cG5c9ZKtr7eWsk2NVn+TNhKds0a+3tlA7ni4BaRFcD/Bq4ASrFaTs2q6n0KHWcgSkqswurUqVZQb/9+OyF55FT0VFbCDTfY1tUF27fHHOSbN5uoL1oUc5DPnDnIB+YviSxDbQE+gLVPrQU+DMxT1U+nb3gX45aFk/O0tMCePXDsmEXzpKk5kDMMVM0SDIUjtArnzo35Oa64IrNinyvLUCKyRVVrReRlVV0U7Numqlcne/BkcLFw8obTp+H1120NfcyY/EwkyxeOHYs5yLdtM3/U5MmxSrmZaCWbK8tQwDkRKQW2B/24jwO+8Oo4yVJdbdnGjY2xbn3jxmWkf7WTINOmwQc/aNvp0/D00yYejzwCDz5o1uHq1WZxrFxpFmOekYhlcQlwAvNX/BUwBviqqu5N3/Auxi0LJy/p7objx83S6O42p6qXD8l+2tosoqquznwcp0+b2C9bFnOQT5iQmmPl0DJUFdCqqj3B42KgTFXPJXVgkfcD/4A5zK9V1SEpgIuFk9d0dsLhw+bT8MZLuUV3txU5DBMBjx61v91VV5lohK1kkyWHxOJZ4EZVbQ4ejwSeUNXrkjqwyBVAD/B14L+4WDhOHK2tFjV1+LA3XspFVK38S+gg37nT9s+aFSuxvmBBYiHUOeSzKA+FAkBVm4NueUmhqjsBvH+S4/RBRQVceaU1Xtq1yxsv5RoicOmltv3xH8Obb5po1NfD974H3/2uLU+FFkdtbdZn+SciFi0ico2qvgggIkuB1vQMy3EcwASitjZWPqShwSKnvPFSbjFlCtx5p21NTbFWso8/Dg89ZA7x+FayWRhOnYhYfAJ4UESOYf0spgB3DvQGEflV8LrefEZVfzbUA4vIemA9QE1NzVDf5jj5w7hxFmUTNl5qavLyIbnK6NHw279tW3u79eQIrY6NG+1vWlsbayWb5LJTqkm0+dEI4PLg4S5V7Yx77p2qujHhAYjU4T4Lxxk6YfmQPXtsbdzLh+QH3d3w6qsxP8fhw7Z/wQJbqlq0CG6/3fI7kiBrakMlW67cxcJxkqS93cuH5CuqsVaydXXw2mu2/2c/g9tuS+ojM+ngHnQsCb1Y5HewWlMTgZ+LyHZVvSmF43Gc/KasDObPt3pFXj4kvxCB2bNtu/tui4R67DFbnoqIVIpFQiaKqj4MPJzC4ztOYVJVZf0YZs2ypL6GBlsX9/Ih+cPEifDud0fqo/KFTsfJF8LyIUuXml+jocGS/BwnBaRSpg6m8LMcx0kGEYueGT/elqV27fLyIU5KSKSfxR197D4DvKKqDara1/OO40RBcbH5MiZPhkOHLJvYy4c4wyARy+KPgJXAU8HjdcBWYLaI/HdV/V6Kx+Y4znApLYXLLoMZM0wwDh+2LPBRo6IemZNjJCIWJcAVqnoCQEQmA/cBy4FNgIuF42QrFRWwcOGF5UNGj87OXtNOVpKIg3tmKBQBDcG+U4B70RwnFxg92sIvly+3WP6GBujoiHpUTg6QiGVRJyKPYm1VAd4X7KsCTqd6YI7jpAkRc4CvWmUF7l5/Hc6c8fIhzoAk8sv4GHAHsDp4fC/wE7UU8HekemCO46SZoiLrADdxIhw5Yol9RUXmBPfyIU4vhiwWqqoi8jTQgSXgPa+pqhXiOE50jBgBc+aYcOzfb9FTYfkQxwkY8uWDiPwe8Dy2/PR7wHMi8r50DcxxnAxTXm5F61avNqE4cQJaWqIelZMlJLIM9Rlgmao2AIjIROBXwI/TMTDHcSJi1CjLAvceGk4ciSxMFoVCEXAywfc7jpNLhD00rr7aKtw2NloZEacgScSy+IWI/BK4P3h8J/BY6ofkOE7WUFRkXd4mTICjR2H3bu+hUaAk4uD+ryLyu8CqYNeGoHKs4zj5TkkJXHKJCYf30ChIEgqqVtWfAD9J01gcx8l2vIdGwTKoWIjIWfruVSFYRO3olI/KcZzsJr6Hxs6dFjk1Zoz30MhjBhULVU1LxTER+SLwHixvYx9wt6qeTsexHMdJE9XVsGKFRUzt2AFNTeYY90zwvCNKD9VGYKGqLgJ2A5+OcCyO4ySLiJVCX7PGihU2NcHJk9ZHw8kbIhMLVX1CVcM4vGeBGVGNxXGcFBD20Fi71npHnzoFb79t0VNOzpMtsW/3AI/396SIrBeRLSKypbGxMYPDchwnYcIeGmvWWNe+hgZobo56VM4wSatYiMivROTVPrb3xr3mM0AX8IP+PkdVN6hqrarWTpw4MZ1DdhwnVVRWwqJFVt22vNyc4G1tUY/KSZK0eqFU9caBnheRu4Bbgd/yooSOk6eMGQPXXgtvvWVO8BMnzAk+YkTUI3MSILKQBRG5GfgUsFZVz0U1DsdxMoCIlUJfvRqOH7ceGt3dlgleXBz16JwhEGV821eAMmCjWAbos6r60QjH4zhOuikutn7gkyfDwYPWF7ykxEJwPRM8q4lMLFT10qiO7ThOxIwYYU7wGTNMMA4fNh/HqLSkdTkpIFuioRzHKUQqKiw3Y/Vqywo/cQJaW6MeldMHLhaO40TP6NFQWwvLl1teRkMDdHREPSonDhcLx3GyAxEYP95CbRctgnPnLILKe2hkBV7AxXGc7KKoCKZPt4S+I0esum1RkUVOuRM8MlwsHMfJTkaMgDlzYNq0mBO8vNyWrJyM48tQjuNkN+XlcOWV5gQfNcqc4Oc8NSvTuFg4jpMbjBoVc4KLuBM8w7hYOI6TW4wfD9ddZ82XWluhsdGd4BnAfRaO4+QeRUUwdaqVEDl82JzgIuYEL/Jr4HTgYuE4Tu5SUhJzgu/fD4cOWZ/wMWOiHlne4RLsOE7uU14OCxbA9ddbtJQ7wVOOi4XjOPnDyJHmBF+xwpajGhqgvT3qUeUFLhaO4+Qf48bBypXmBG9vdyd4CnCfheM4Wcfpcx1UV5YO70PineBvvAG7d9t+d4InhX9jjuNkFS3tXTz68nGa21NkCZSUwKxZsHYtzJxp9aZOn07NZxcQkYmFiPwPEXlZRLaLyBMiMi2qsTiOkz3sONbE3/30VXYea0rtB5eVwRVXmBN87Fhzgre0pPYYeUyUlsUXVXWRqi4BHgU+G+FYHMfJAlrau9iweR8A39i8P3XWRTwjR8I115hPo6TERMOd4IMSmVioavxlQxWgUY3FcZzsYMexJjbuaADgiR0nUm9dxDN2rAnGNde4E3wIROrgFpHPAR8GzgDviHIsjuNES7xVEfKNzfu5YtpoRpal6VQlAlOmwIQJcPSoOcFV3QneB6Kavgt6EfkVMKWPpz6jqj+Le92ngXJV/ft+Pmc9sB6gpqZm6aFDh9IxXMdxIuREUyt7G1qIPyWJwKWTqpg8uiIzg2hvhwMHLBs8zATPlh4ajY1mBU2alNTbRWSrqtYme/i0isWQByFSAzymqgsHe21tba1u2bIlA6NyHKdgaWkxK+PNN6Gy0vwcUROxWEQZDXVZ3MP3Aq9HNRbHcZwLqKqCq6+2TPARIywTvK0t6lFFSpQ+i8+LyOVAD3AI+GiEY3GcgiQlyW/5TOgEb2iAHTugqcmyw0sKL585shmr6u9GdWzHcWLJb7dfPT19DuR8QAQmTzYneJgJXoBO8MKZqeM4F5C25Ld8pbgYLrkE1qyBmhrzIZw+DVng980ELhaOU4BkJPktXykrg/nzrXzI+PG2RNXcHPWo0o6LheMUIBlNfstXqqqsqu1110FpqWWC57ET3MXCcQqM/pLf3LpIkupqi5pauhQ6OszSyMNMcPdqOU6B0dzeyV3XzeYjK2ef3ycCLe2d7uhOlgJwgvsvw3EKjMmjKzKXEV1ohE7wKVPg4EHLBC8tza5M8CRxsXAcx0k1ZWVw+eUwYwbs2QPHjpmPIxsywZMkP+wjx3GcbCSPnOAuFo7jOOkmD5zgvgzlOI6TCeKd4EePwq5dOeUEd7FwHMfJJMXFlgEeOsH37csJJ7iLheM4ThSUlsK8eTB9ek44wbPf9nEcx8ln4p3gZWVZWw7dxcJxHCcbqK6G5cvNCd7ZmXVOcF+GchzHyRZErBPe+PG2LPX66+YEr66OemQuFo7jOFlHcTHMnGnRU6ETPOKlqciXoUTkP4uIisiEqMfiOI6TVYRO8LVr4bLL7HFERGpZiMhM4F3A4SjH4TiOk9VUVsLChZEOIWrL4l+ATwGF0WrKcRwnR4lMLETkvcBRVX1pCK9dLyJbRGRLY2NjBkbnOI7jxJPWZSgR+RUwpY+nPgP8LbYENSiqugHYAFBbW+tWiOM4ToZJq1io6o197ReRq4DZwEti6e0zgBdF5FpVfTOdY3Icx3ESJxIHt6q+AkwKH4vIQaBWVd+KYjyO4zjOwETt4HYcx3FygKxIylPVWVGPwXEcx+kftywcx3GcQRHV3AouEpFG4FCSb58AFKpfpJDnDoU9/0KeOxT2/OPnfomqTkz2g3JOLIaDiGxR1dqoxxEFhTx3KOz5F/LcobDnn8q5+zKU4ziOMyguFo7jOM6gFJpYbIh6ABFSyHOHwp5/Ic8dCnv+KZt7QfksHMdxnOQoNMvCcRzHSQIXC8dxHGdQCkYsRORmEdklIntF5G+iHk8qEJFvi0iDiLwat2+ciGwUkT3B7dhgv4jIl4P5vywi18S95yPB6/eIyEeimEuiiMhMEXlKRHaIyGsi8pfB/kKZf7mIPC8iLwXz//+C/bNF5Llgng+ISGmwvyx4vDd4flbcZ3062L9LRG6KaEoJIyLFIrJNRB4NHhfE3EXkoIi8IiLbRWRLsC/9v3tVzfsNKAb2AXOAUuAlYEHU40rBvNYA1wCvxu37J+Bvgvt/A3whuP9u4HFAgBXAc8H+ccD+4HZscH9s1HMbwtynAtcE90cBu4EFBTR/AUYG90cAzwXz+hHwgWD/14D/FNz/M+Brwf0PAA8E9xcE/w9lWCXofUBx1PMb4nfwSeDfgEeDxwUxd+AgMKHXvrT/7gvFsrgW2Kuq+1W1A/gh8N6IxzRsVHUTcKrX7vcC9wb37wVuj9t/nxrPAtUiMhW4CdioqqdU9W1gI3Bz2gc/TFT1uKq+GNw/C+wEplM481dVbQ4ejgg2BW4Afhzs7z3/8Hv5MfBbYv0B3gv8UFXbVfUAsBf7f8lqRGQGcAvwzeCxUCBz74e0/+4LRSymA0fiHr8R7MtHJqvq8eD+m8Dk4H5/30HOfzfBssLV2NV1wcw/WIbZDjRg/+z7gNOq2hW8JH4u5+cZPH8GGE/uzv9fsZbMPcHj8RTO3BV4QkS2isj6YF/af/dZUXXWSQ+qqiKS17HRIjIS+AnwCVVtsgtGI9/nr6rdwBIRqQYeBuZHO6LMICK3Ag2qulVE1kU8nChYrapHRWQSsFFEXo9/Ml2/+0KxLI4CM+Mezwj25SMnAjOT4LYh2N/fd5Cz342IjMCE4geq+lCwu2DmH6Kqp4GngJXYMkN4ERg/l/PzDJ4fA5wkN+e/CrhNrGnaD7Hlpy9RGHNHVY8Gtw3YRcK1ZOB3Xyhi8QJwWRAtUYo5uR6JeEzp4hEgjGz4CPCzuP0fDqIjVgBnArP1l8C7RGRsEEHxrmBfVhOsOX8L2Kmq/xz3VKHMf2JgUSAiFcA7Mb/NU8D7gpf1nn/4vbwP+A81T+cjwAeCiKHZwGXA8xmZRJKo6qdVdYZaH5wPYHP5Awpg7iJSJSKjwvvY7/VVMvG7j9qzn6kNiwrYja3rfibq8aRoTvcDx4FObM3xj7C12CeBPcCvgHHBawX4P8H8X8Ha2Iafcw/m3NsL3B31vIY499XY2u3LwPZge3cBzX8RsC2Y/6vAZ4P9c7AT3l7gQaAs2F8ePN4bPD8n7rM+E3wvu4DfjnpuCX4P64hFQ+X93IM5vhRsr4Xnskz87r3ch+M4jjMohbIM5TiO4wwDFwvHcRxnUFwsHMdxnEFxsXAcx3EGxcXCcRzHGRQXC8fph6C65wQRqRaRP4vbP01EfhzcXyIi707is/9BRP5LKsfrOOnExcJxBqcaq1wKgKoeU9Uw+WsJlt/hOHmNi4XjACLy06Aw22txxdlCPg/MDfoHfFFEZonIq0E1gP8O3Bk8d2dviyF43azg/mdEZLeIPA1cHveauSLyi+D4m0WkIGo8ObmFFxJ0HOMeVT0VlM54QUR+Evfc3wALVXUJnK9yi6p2iMhnsazYjwfP/UNfHy4iS7HSFEuw/7sXga3B0xuAj6rqHhFZDnwVq3fkOFmDi4XjGH8hIr8T3J+J1QlKJdcDD6vqOQAReSS4HQlcBzwYVzG3LMXHdpxh42LhFDxBmesbgZWqek5E6rB6QsnQxYXLu4N9ThHWh2FJksdznIzgPgvHsZLVbwdCMR9rPxnPWax1a1/0fu4g1uqWoN/x7GD/JuB2EakIqoa+B0BVm4ADIvL+4D0iIouHPyXHSS0uFo4DvwBKRGQn5sx+Nv5JVT0J/DpwVn+x13ufAhaEDm6sv8Y4EXkN+DhW6Ri1FrAPYNVCH8fK5of8AfBHIhJWEs35lr9O/uFVZx3HcZxBccvCcRzHGRQXC8dxHGdQXCwcx3GcQXGxcBzHcQbFxcJxHMcZFBcLx3EcZ1BcLBzHcZxB+X/CLuCctuNw2wAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"weighted_data['log_cases_pop_den'] = np.log(weighted_data['cases_pop_den'])\n",
"\n",
"ax = sns.scatterplot(\n",
" data=weighted_data,\n",
" x='altitude',\n",
" y='log_cases_pop_den',\n",
" marker='^'\n",
")\n",
"\n",
"prediction = weighted_model.get_prediction(weighted_data['altitude']).summary_frame()\n",
"ax.plot(\n",
" weighted_data['altitude'], \n",
" prediction['mean'],\n",
" color='red'\n",
")\n",
"ax.fill_between(\n",
" weighted_data['altitude'], \n",
" prediction['mean_ci_lower'],\n",
" prediction['mean_ci_upper'],\n",
" color='red',\n",
" alpha=.2\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment