{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Logistic Regression\n", "\n", "Want to know how different inputs affects a yes/no answer? Logistic regression is your new best friend! While it's best if you've read up on linear regression first, you'll probably survive even if you didn't." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "<p class=\"reading-options\">\n <a class=\"btn\" href=\"/regression/logistic-regression\">\n <i class=\"fa fa-sm fa-book\"></i>\n Read online\n </a>\n <a class=\"btn\" href=\"/regression/notebooks/Logistic Regression.ipynb\">\n <i class=\"fa fa-sm fa-download\"></i>\n Download notebook\n </a>\n <a class=\"btn\" href=\"https://colab.research.google.com/github/littlecolumns/ds4j-notebooks/blob/master/regression/notebooks/Logistic Regression.ipynb\" target=\"_new\">\n <i class=\"fa fa-sm fa-laptop\"></i>\n Interactive version\n </a>\n</p>" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "**We make scarves, and we're pretty lazy.** Since we're \ud83e\uddf6\ud83d\udc51Crafty Stylish People\ud83d\udc51\ud83e\uddf6, we might have some big questions about the scarves we knit, such as:\n", "\n", "* If we're knitting a **55-inch scarf**, how likely are we to finish it?\n", "* If we're knitting a **70-inch scarf**, how likely are we to finish it?\n", "* If we're knitting a **80-inch scarf**, how likely are we to finish it?\n", "\n", "We've been keeping records for the past few weeks, and we hope they look something like this:\n", "\n", "|\ud83d\udccf|\ud83e\udde3|\n", "|---|---|\n", "|55 inches|**Finished!**|\n", "|60 inches|**Finished!**|\n", "|65 inches|**Finished!**|\n", "|70 inches|Nope|\n", "|75 inches|Nope|\n", "|80 inches|Nope|\n", "\n", "That way we can says something like, \"if the scarf is over 65 inches, we're probably too lazy to finish it!\" Unfortunately our records do _not_ look like that. Instead, they look something like this:\n", "\n", "|\ud83d\udccf|\ud83e\udde3|\n", "|---|---|\n", "|55 inches|**Finished!**|\n", "|55 inches|**Finished!**|\n", "|55 inches|**Finished!**|\n", "|60 inches|**Finished!**|\n", "|60 inches|Nope|\n", "|70 inches|**Finished!**|\n", "|70 inches|Nope|\n", "|82 inches|**Finished!**|\n", "|82 inches|Nope|\n", "|82 inches|Nope|\n", "|82 inches|Nope|\n", "\n", "Right now our questions are about our chances of finishing a 55-inch, 70-inch and 80-inch scarf. We can _try_ to answer some of our questions with this data: for 55-inch scarves we've **finished all of them**, and we've finished **half of the 70-inch scarves**. There's a problem investigating our chances of finishing an 80-inch scarf, though: **we don't have any 80-inch scarves in our dataset!** \n", "\n", "To reason about our chances of finishing an 80-inch scarf, let's think about what we do know:\n", "\n", "* We rarely finish 82-inch scarves (a little longer than 80 inches)\n", "* We finish about half of the 70-inch and 60-inch scarves (a good amount shorter)\n", "* We're very good at finishing 55-inch scarves (much much shorter)\n", "\n", "It sounds like we're **generally bad at finishing longer scarves,** which is totally understandable. But _how_ bad?\n", "\n", "**This is where logistic regression rescues us.** Logistic regression can tell us two things:\n", "\n", "* What the effect of length is on our ability to complete a scarf\n", "* Given a length, what are our chances of completing it?\n", "\n", "Completing our scarf is a yes/no question, right? **Logistic regression is all about predicting categories** - in this case, the category yes vs the category no. If you're predicting a number instead, you use linear regression. That's it!\n", "\n", "Let's see what this looks like with a little code.\n", "\n", "> For the questioning souls: Yes, the \"chances of finishing\" are a number, but the inputs are yes/no or completed/incomplete. That makes it a logistic regression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Performing a logistic regression\n", "\n", "We'll start off with our data. We're going to use [pandas](https://pandas.pydata.org/), a super-popular Python library for doing data-y things." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>length_in</th>\n", " <th>completed</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>0</th>\n", " <td>55</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", " <td>55</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>2</th>\n", " <td>55</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>3</th>\n", " <td>60</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>4</th>\n", " <td>60</td>\n", " <td>0</td>\n", " </tr>\n", " <tr>\n", " <th>5</th>\n", " <td>70</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>6</th>\n", " <td>70</td>\n", " <td>0</td>\n", " </tr>\n", " <tr>\n", " <th>7</th>\n", " <td>82</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>8</th>\n", " <td>82</td>\n", " <td>0</td>\n", " </tr>\n", " <tr>\n", " <th>9</th>\n", " <td>82</td>\n", " <td>0</td>\n", " </tr>\n", " <tr>\n", " <th>10</th>\n", " <td>82</td>\n", " <td>0</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " length_in completed\n", "0 55 1\n", "1 55 1\n", "2 55 1\n", "3 60 1\n", "4 60 0\n", "5 70 1\n", "6 70 0\n", "7 82 1\n", "8 82 0\n", "9 82 0\n", "10 82 0" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "import numpy as np\n", "\n", "df = pd.DataFrame([\n", " { 'length_in': 55, 'completed': 1 },\n", " { 'length_in': 55, 'completed': 1 },\n", " { 'length_in': 55, 'completed': 1 },\n", " { 'length_in': 60, 'completed': 1 },\n", " { 'length_in': 60, 'completed': 0 },\n", " { 'length_in': 70, 'completed': 1 },\n", " { 'length_in': 70, 'completed': 0 },\n", " { 'length_in': 82, 'completed': 1 },\n", " { 'length_in': 82, 'completed': 0 },\n", " { 'length_in': 82, 'completed': 0 },\n", " { 'length_in': 82, 'completed': 0 },\n", "])\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our very tiny dataset has two columns:\n", "\n", "* Number of inches in the atttempted scarf\n", "* Whether we completed it or not. 1 means yes, 0 means no\n", "\n", "We want to ask a simple question using this data: **how does length affect our ability to finish a scarf?** Since finishing a scarf is a yes/no question, we get to use **logistic regression**.\n", "\n", "### Writing a regression\n", "\n", "To perform our logistic regression, we're going to use a library called [statsmodels](https://www.statsmodels.org), the same one we used for linear regression.\n", "\n", "We'll use `smf.logit` to create a new logistic regression, then use `model.fit()` to teach it all about scarf-making." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 0.531806\n", " Iterations 5\n" ] } ], "source": [ "import statsmodels.formula.api as smf\n", "\n", "# What effect does the length of the scarf have one whether it was completed?\n", "model = smf.logit(formula='completed ~ length_in', data=df)\n", "results = model.fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The technique of writing out the relationship like `\"completed ~ length_in\"` was stolen from R. It uses a library called [Patsy](https://patsy.readthedocs.io/), and you can read all about it in [the statsmodels documentation](https://www.statsmodels.org/stable/example_formulas.html).\n", "\n", "> When using statsmodels, linear regression is `smf.ols` while logistic is `smf.logit`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making predictions\n", "\n", "Now that we've taught our regression all about how bad we are at making scarves, let's **make some predictions** to see what it learned. Along with the 80-inch scarf, let's try out a few more sizes as well." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>length_in</th>\n", " <th>predicted</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>0</th>\n", " <td>55</td>\n", " <td>0.850526</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", " <td>65</td>\n", " <td>0.651815</td>\n", " </tr>\n", " <tr>\n", " <th>2</th>\n", " <td>75</td>\n", " <td>0.381148</td>\n", " </tr>\n", " <tr>\n", " <th>3</th>\n", " <td>80</td>\n", " <td>0.261047</td>\n", " </tr>\n", " <tr>\n", " <th>4</th>\n", " <td>90</td>\n", " <td>0.104122</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " length_in predicted\n", "0 55 0.850526\n", "1 65 0.651815\n", "2 75 0.381148\n", "3 80 0.261047\n", "4 90 0.104122" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# What lengths do we want to ask it about?\n", "unknown = pd.DataFrame([\n", " { 'length_in': 55 },\n", " { 'length_in': 65 },\n", " { 'length_in': 75 },\n", " { 'length_in': 80 },\n", " { 'length_in': 90 },\n", "])\n", "\n", "# Save the predictions into a new column\n", "unknown['predicted'] = results.predict(unknown)\n", "unknown" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "According to the regression, we have a **26% chance to finish a 80-inch scarf**. While I'd like to say we can do better.... like I said, we're lazy.\n", "\n", "On the other hand, we only have a 85% chance for our 55-inch scarves!! Even though we finished *every single 55-inch scarf we've ever made*, it turns out that not completing a handful of 60-, 70-, and 82-inch scarves has made our regression a little pessimistic." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Examining our results\n", "\n", "Usually when you're doing logistic regression you don't care too much about making predictions. You're more concerned with figuring out things like how good your regression was, and - in this case - **how much each inch adds to our chance of failure**.\n", "\n", "There are a few ways to look at the results, but we'll only use the fanciest since we like \u2728\ud83c\udf1f\ud83d\udc8e \ud835\udcbb\ud835\udcb6\ud835\udcc3\ud835\udcb8\ud835\udcce \ud835\udcc9\ud835\udcbd\ud835\udcbe\ud835\udcc3\ud835\udc54\ud835\udcc8 \ud83d\udc8e\ud83c\udf1f\u2728. To get the results we just look at `results.summary()`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<table class=\"simpletable\">\n", "<caption>Logit Regression Results</caption>\n", "<tr>\n", " <th>Dep. Variable:</th> <td>completed</td> <th> No. Observations: </th> <td> 11</td> \n", "</tr>\n", "<tr>\n", " <th>Model:</th> <td>Logit</td> <th> Df Residuals: </th> <td> 9</td> \n", "</tr>\n", "<tr>\n", " <th>Method:</th> <td>MLE</td> <th> Df Model: </th> <td> 1</td> \n", "</tr>\n", "<tr>\n", " <th>Date:</th> <td>Mon, 30 Dec 2019</td> <th> Pseudo R-squ.: </th> <td>0.2282</td> \n", "</tr>\n", "<tr>\n", " <th>Time:</th> <td>12:07:56</td> <th> Log-Likelihood: </th> <td> -5.8499</td>\n", "</tr>\n", "<tr>\n", " <th>converged:</th> <td>True</td> <th> LL-Null: </th> <td> -7.5791</td>\n", "</tr>\n", "<tr>\n", " <th>Covariance Type:</th> <td>nonrobust</td> <th> LLR p-value: </th> <td>0.06293</td>\n", "</tr>\n", "</table>\n", "<table class=\"simpletable\">\n", "<tr>\n", " <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n", "</tr>\n", "<tr>\n", " <th>Intercept</th> <td> 7.8531</td> <td> 4.736</td> <td> 1.658</td> <td> 0.097</td> <td> -1.429</td> <td> 17.135</td>\n", "</tr>\n", "<tr>\n", " <th>length_in</th> <td> -0.1112</td> <td> 0.067</td> <td> -1.649</td> <td> 0.099</td> <td> -0.243</td> <td> 0.021</td>\n", "</tr>\n", "</table>" ], "text/plain": [ "<class 'statsmodels.iolib.summary.Summary'>\n", "\"\"\"\n", " Logit Regression Results \n", "==============================================================================\n", "Dep. Variable: completed No. Observations: 11\n", "Model: Logit Df Residuals: 9\n", "Method: MLE Df Model: 1\n", "Date: Mon, 30 Dec 2019 Pseudo R-squ.: 0.2282\n", "Time: 12:07:56 Log-Likelihood: -5.8499\n", "converged: True LL-Null: -7.5791\n", "Covariance Type: nonrobust LLR p-value: 0.06293\n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 7.8531 4.736 1.658 0.097 -1.429 17.135\n", "length_in -0.1112 0.067 -1.649 0.099 -0.243 0.021\n", "==============================================================================\n", "\"\"\"" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The part we're interested in is down at the bottom, where it says `length_in` and `coef`. `coef` stands for **coefficient**, and it's (kind of) the \"how much\" in \"how much does length affect completion.\"\n", "\n", "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZ0AAACFCAIAAAAPT7K3AAAYZ2lDQ1BEaXNwbGF5AABYhZV5eThUf9j+c+acmTFjjH3f933fyb7v+y4Mxr41dlqIZKlEUkKLSpFosyTtSiJR0aZSWohKhRLl/YP6br/r917v88ec57o/93M/27nOuWYGQDCHlpQUh2MDiE9IYbjbmIv5+vmLEceBGQSBBXhAnBaWnGTm6uoIAPD7+k+bHwEEAOC+Mi0pKe6/5/9fYw+nJ4cBIK4AEBqeHBYPgJwDwPLDkhgpAAQ9AJBMT0lKASAEAgAXw9fPH4CQBABckSt+PgBwha74ewGAi+HpbgFAaAJgotBojEgAagcAiKWFRaYAUB8BEDgSwqMTALimAQjGYVG0cABBJQBQio9PDAcQ9AUAudC/6UT+QzP0jyaNFvnHX+kFAACYLKOTk+Jomf/HcfzvFh+X+juHDABQohi27gDABYA8ik10cAcACgAynRDq7AIAHADI9+hwgBUfR45KtfVa4eOEwpIt/AGABwCnFk6zdAAAIQCcdUKcs+MqHhoRbW0HAGwAuIzoFDvP1dhCerKVx6pmDSPR3eW3H8GwMFuNbaYxAFb53amxXmar+o+i6Ha/9Wezojx9VmpGyWnR3s4AQAVAeZJjPRxWOKhUVpSF828OI9XdCwCkAFADeoKN+Yo+GhTBsHZf5TPik3/3ixZGRds5r/r7UqI8bVd1msJoVh4AwAeAdtATzLx+69CTfR1/9xJOt7Ra6R0dpCd4rfaLjiWlmLuvxn5JinNd5WNkepyNOwBIAGBCyWkeq7GYcQrDc3VHmHNSiqvnSp1YaAzN3nWlHiwDHMECLEEMUkEMQiERYiB6YLp9GsRWT6yBBgyIBDooryK/I3yABgxIABp4QBZ8gASgQ/KfOHOgAQPokAYJ8PMPuvKpDBFAAwakAR2SIRbeAgPiwQHigA6pwAA6JPzJ5g2vgQHR/8keBokQB4nAgOj/B2YGFuC4iqT+1hVj/c0kWBEsCbYEa4I8JoAZY4aYI2aMmWLGmAamh+n/rvYvPv4tfgj/Cj+MH8M/Do7ezPhXP2LgBGOQujorOoT+vWdMBtPAtDFzzAgzxvRBDOPBBEAZ08L0MDPMBDPEtDF9sFitPBX+q/2PHv429VUeSY2EI/GSTEly/46kKlC1/6jQIeEfE1qpNfTPXC3+nPw7v8XfJh0OieDwbyZaiJ5Fe9CraC/ahbaDGHoZ7UD70Yto+9/uotfAgMg/2dyBDgkQC3EQ/Z98tNWcDKBDstoJtUm1pZWzFHpGCgCARWJSJiM6MipFzCwpKY4uZpcQpqIkpqGmrg/g6+cvtvKY+uoOCAAgPHf/wugTAGvYAEiDf2ExuwAabwDwFv+FyQQA8CsBnL4XlspIW8EwAAA8kIEVuIAfREAS5EAZNEAHDMEUrMAeXMAT/CAIwiAK4oEB6bAecqEASmAn7IZ9cAAOw3E4CWegHbrgKtyEPhiEYXgKY/AGpmAG5mERQRAiwoJwIvyIKCKNKCIaiB5ijFghjog74oeEIJFIApKKrEfykBKkHNmHHEIakNPIeeQq0osMIY+Rl8gk8gX5gUNxFBwXThgng1PF6eHMcA44T9xaXCRuHS4Ll4/bgduLq8M14dpwV3F9uGHcGG4KN4cCyozyoOKoMqqHWqAuqD8agTLQjWgxWonWoc1oJ9qD3kfH0Gl0ASNgnJgYpowZYraYFxaGrcM2YtuwfdhxrA3rxu5jL7EZ7BeeBS+EV8Qb4O3wvvhIfDq+AF+Jr8e34m/gh/Fv8PMEAoGHIEvQJdgS/AgxhGzCNkItoYVwhTBEGCfMEYlEfqIi0YjoQqQRU4gFxCpiE/Ey8R7xDfE7EzOTKJMGkzWTP1MC02amSqZGpktM95jeMS2S2EjSJAOSCymclEkqJR0hdZLukt6QFsnsZFmyEdmTHEPOJe8lN5NvkEfJX5mZmSWY9ZndmKOZc5j3Mp9ivsX8knmBwkFRoFhQAimplB2UY5QrlMeUrywsLDIspiz+LCksO1gaWK6zPGf5TuWkqlDtqOHUTdRqahv1HvUjK4lVmtWMNYg1i7WS9SzrXdZpNhKbDJsFG41tI1s123m2h2xz7Jzs6uwu7PHs29gb2XvZJziIHDIcVhzhHPkchzmuc4xzopySnBacYZx5nEc4b3C+4SJwyXLZccVwlXCd5BrgmuHm4Nbi9ubO4K7mvsg9xoPyyPDY8cTxlPKc4Rnh+cErzGvGS+ct4m3mvcf7jU+Qz5SPzlfM18I3zPeDX4zfij+Wv4y/nf+ZACagIOAmkC6wX+CGwLQgl6ChYJhgseAZwSdCOCEFIXehbKHDQv1Cc8IiwjbCScJVwteFp0V4RExFYkQqRC6JTIpyihqLRotWiF4WfS/GLWYmFie2V6xbbEZcSNxWPFX8kPiA+KKErISXxGaJFolnkmRJPckIyQrJa5IzUqJSTlLrpU5IPZEmSetJR0nvke6R/iYjK+Mjs1WmXWZClk/WTjZL9oTsqByLnIncOrk6uQfyBHk9+Vj5WvlBBZyCtkKUQrXCXUWcoo5itGKt4pASXklfKUGpTumhMkXZTDlN+YTySxUeFUeVzSrtKh9VpVT9VctUe1R/qWmrxakdUXuqzqFur75ZvVP9i4aCRphGtcYDTRZNa81Nmh2an7UUteha+7UeaXNqO2lv1b6m/VNHV4eh06wzqSulG6Jbo/tQj0vPVW+b3i19vL65/ib9Lv0FAx2DFIMzBp8MlQ1jDRsNJ9bIrqGvObJm3EjCiGZ0yGjMWMw4xPig8ZiJuAnNpM7klamkabhpvek7M3mzGLMms4/mauYM81bzbxYGFhssrliiljaWxZYDVhxWXlb7rJ5bS1hHWp+wnrHRtsm2uWKLt3WwLbN9aCdsF2bXYDdjr2u/wb7bgeLg4bDP4ZWjgiPDsdMJ52TvtMtp1FnaOcG53QVc7Fx2uTxzlXVd53rBjeDm6lbt9tZd3X29e48Hp0ewR6PHvKe5Z6nnUy85r1Sva96s3oHeDd7ffCx9yn3GfFV9N/j2+Qn4Rft1+BP9vf3r/ecCrAJ2B7wJ1A4sCBxZK7s2Y21vkEBQXNDFYNZgWvDZEHyIT0hjyBLNhVZHmwu1C60JnQmzCNsTNhVuGl4RPkk3opfT30UYRZRHTEQaRe6KnIwyiaqMmo62iN4X/TnGNuZAzLdYl9hjsctxPnEt8UzxIfHnEzgSYhO6E0USMxKHkhSTCpLG1hms271uhuHAqE9Gktcmd6RwpSSl9KfKpW5JfZlmnFad9j3dO/1sBntGQkZ/pkJmUea7LOuso9lYdlj2tfXi63PXv9xgtuHQRmRj6MZrmyQ35W96k2OTczyXnBube2ez2ubyzbN5Pnmd+cL5OfnjW2y2nCigFjAKHm413HqgECuMLhwo0iyqKvpVHF58u0StpLJkaVvYttvb1bfv3b68I2LHQKlO6f6dhJ0JO0fKTMqOl7OXZ5WP73La1VYhVlFcMbs7eHdvpVblgT3kPal7xvY67u2okqraWbW0L2rfcLV5dUuNUE1Rzbfa8Np7+033Nx8QPlBy4MfB6IOPDtkcaquTqas8TDicdvjtEe8jPUf1jjbUC9SX1P88lnBs7Lj78e4G3YaGRqHG0hO4E6knJpsCmwZPWp7saFZuPtTC01JyCk6lnnp/OuT0yBmHM9fO6p1tPid9rqaVs7W4DWnLbJtpj2of6/DrGDpvf/5ap2Fn6wWVC8e6xLuqL3JfLL1EvpR/afly1uW5K0lXpq9GXh2/Fnzt6XXf6w+63boHbjjcuHXT+ub1HrOey7eMbnX1GvSev613u71Pp6+tX7u/9Y72ndYBnYG2u7p3Owb1BzuH1gxdumdy7+p9y/s3H9g96Bt2Hh4a8Rp59DDw4dij8EcTj+Mef36S9mTxac4ofrT4GduzyudCz+teyL9oGdMZu/jS8mX/K49XT8fDxqdeJ79eepP/luVt5TvRdw0TGhNdk9aTg+8D3r+ZSppanC74wP6h5qPcx3OfTD/1z/jOvPnM+Lz8ZdtX/q/HZrVmr825zj2fj59f/Fb8nf/78QW9hZ4fPj/eLaYvEZf2/pT/2fnL4dfocvzychKNQQMAABQAcBERAF+OAbD4AXAOApADVr7nrRqKAOAAwBuxwpmhehgfnkxgIqox+ZHyyJcpBBYatZ2NzB7HcZtLm7uGF/hi+QcEdYR2Ck+JmoqVig9JkqX0pf1kYmXj5QLlzRWEFT4r3lSqUo5VMVJlUX2h1qKeo+GmKa75Qeu89hYdN10h3Td6zfoZBmaGZMP7a2qMwo2VjL+YtJuuNzM3p5i/sLhk2WhVa11ms9GWZmdiz2f/2aHfsdmp1vmQS5fruDveg99TwIvNG/Ve8ln0A39SADWQZS22di7oVfBgyBXa2dD6sKrwYnpmRGSkZ5R5tFaMQqx4HH88awKaMJv4Kmlw3QXGkeQdKZtSC9JaM7BMetaV9bBBZqPBJrucgNzUzTvydudnb9HaMl5QutW1ULqIuRhKcNvYt8vtMC513ulT5l/uv8u3wnu3Z6XbHue9DlU2+8yrjWv0azX3Kx9QOKh2yKEu7/DYUbv6pmNTDeyN0ifUmwxPWjY7tficCj4ddSbpbPq5ja2b27a0F3aUnC/t3H2hpqv+4rlLNy4/vDJ2deRay/WIbr7uWzcqb6b3RNxa2+tz263Pod/mju2A5911gweHHt9nfqA6bDFi99Dqkd5j6SfUJwtPJ0YfPbv6/PCLvLHIl16vnMedXru8cXlr/05/gndibLL4vdb7sanj01kfbD8yfWz4ZPNpfObw54wvQV9dZp3mYuavfd/6o/2n5fLy6v7VUQydxMbw44QZJpSkQ45irqGMURVY09lucvBzZnI94NHg3cz3TEBbsEBoUERA1FesTLxLYlRyTmpe+r3MHdnDcgx5YwUmhQeKB5RilLWVf6ncVN2h5qMuqv5Oo1kzTctIG9G+oVOs66LHqTeiX2UQYChsOLqmxijQmN/4ocke0wAzGbNF82GL05bbrOjWa2zYbd7adtnttk9zoDuGOkU5J7rEu4a6ubgbeih4CnpRvXHe8z7vfEf8rvs3B1QHFq/NCooO9g2xpKmG8oUhYe/Dh+ndEa2R9VGV0fkxibF+cabxsgksCbOJL5Oer5tNFk8JTq1Ku5r+KGM8czprYT3zBpGNcpvEcgg5L3JbN5fmMfKDtngV+G6NLswrqi0+WdK6rW37uR2nS0/ubCg7Wn5wV3XF7t2llUV7Nu/NrErcF1kdXZNTe/mA/MHjdbKHy4/cP7pwjHpcoEGyUeGEepPuSeNmyxanU36n484UnD187lLrUNvz9omOr53oBd4uxYuGl0wv614Rv4q7+upaz/XW7mM3qm/u7NlyK6uXcTulr6i/a4Dn7obBZ/cE7ps88ByOGMl5ePTR3cezTzlGlZ85Pk96sWfswst7r56Pv3o99Rb/Tn8id3Join1a7YP2R5lPrJ++z7z9/PDL7a/nZw/NbZr3/ib7bf5710LWD8NFypLlz8nV/asgU7haNAiTxxPxnwmTxPdMr0ifmckUaRYzqj9rLlsT+xDHMpc0txVPDO8WvgP85wRuCN4Suil8QeSQaIaYudgP8SMSDhJTkoVSslLXpIOkF2QqZNVkb8tFyhPljynYKrxTLFCSU7qhHKYCKrWqa1QfqaWqc6i3aDhqTGjmaYlodWi7a0/rbNEV1W3Xc9Wb0N9kwGNwwtDM8N6asDUfjbKNicbVJlomI6ZZZiJmHeYu5o8toiyWLeusXK1J1tdt1ttq2b63q7MPdOBzGHHc7eThzOrc65Lnaug669biHush6/Ha85DXWm9+7wc+pb62vst+rf5xAVIBzwIr1zqvnQ+qCJYOPhdiFvKElhEqEfoorCo8im4ToRupH2UXTYuJj6XFmcSzxY8mHE2MT9JOWlp3nVGc7JrCnfI09UBaeLpM+tuM/ZlWmaNZcdlc2ffXX9hwaWP3pus553MbNlfm5eUnbgkosNqqUIgvfFBUVexfIlWyuG1s+50d50sP7txYFlBusEtg10LFyO4zlXv2bN9bXnVo39nqmzWPat/vXzzIckisTvOw7ZHAo4n1G48VHd/WkNNIO6HbRG36cvJD88IpymmRMxpnXc9lt55r+96hfz6ps+rCqa6Oixcu9V6eu2pz7Xy3x425nspezdsP+rcPhAza3TN7YD4S95g6OvVq4P3c7MLyMsDK730AAAQdgF25AL4FAF4mAGXdADLDALxkAFcWAE99wMlEAI7SD4jByz/vDwRQIAAzsAMfiIIsqIEBWIIL+EMEJEMulMJ+aIZLcBdewixCRIQQdcQGCUbSkTKkCbmFvMURcHI4R1wyrhbXh1tG9dA09Dz6C7PBdmGv8Jr4QvwLggGhirBIpBFvM+kyHSMJksrIzOQiZjLzTooA5RiLFksX1YjayarHeoHNlu0pewoHG8dJTkvOIS5PriFuF+57PME833mr+Iz4nvNvEBAU6BQMEiIJdQmniWiJfBU9I8YQ1xZfkuiRrJSKkl4jQ5UZkz0rVygfqmCmKKNEVVpU/qjyWnVYrVU9W0Nd47lmoZa21iftDp1y3Uy9cH1HAzVD3jVUIxXjalNFs+3mvRafrJisuW34bYXspOy1HJwd1zntde52+eIm6e7jscOzxxvzsfQt8OsP4AkMXdsY9DqEQGMPJYTOhb0JH6W/j2SNcojeHfMubk18RcLHJPt1jcmUlHWpT9KtMzqylLPrN4htrM7hyS3LI+fnbpnbGlM4VVyyLX5Haxn7LoGKD5UNe4P38VQP1m4/YHNwrq70CNfRwvr547ENX07sPGnVwn7q85m35ybapjredY53fb7Me9XietCNkB6PXpM+1Tvyd3WGEu5/f4g9IY0eeMH58tIb6sT6KbMPLZ8Wv+jMWs+Tv23/fnth4sebxcdL537u/BW6rLa8vLp/IlCAA/hBHBRAE4zAFjwhBOIhG4qgChrgPPTBM5hB8IgAoobYIMFIJlKBnEIGkA84Vpwmzh+XhzuDe4OKosHoEXQa08HysWG8PD4XP0owIFQTgRhFHGayYuogqZIayfLkJmYt5ssUV8o4SwaVRN3NKs56is2C7Sl7OgcPRzunN+cHrg3cZO69PMo8t3kT+Xj5rvBHC3AJXBFMFJISGhWuEvEV5RN9LFYrHi6hJgmSD6ROSOfLBMpqybHIvZfvVzireECpVDlPZb1qilqYuqkGRWNAs1jLQZtX+7POY90evTb9OoNthllr0oyKjDtMvplpmodblFjWW7VZX7C5YHvRrtf+pSPOScHZ22WLa7vbtIeUZ7BXrfdzXwm/GP+2QOJan6B9wTdChmjXQhvCCsOj6e4RtpF+UZujr8SyxIXGdyUKJGWte5ZsntKQxpqelNGXJZ6dtn5wo/amI7mCmyvyyVuyC6YLaUWvSrK2q5Xidj4rP12RVqm150vV6erUWoP9Pw7W12kcrj3yrl72WNTxU428J2pOGjV/OFV1Rv/sQCutbbGjrtOtCy42XHa88vnage7Qmwa3xG9jfXfupN0lDBbfo9yvGw5+6PQ47umxZ+/GRF+5vs59e2mSd2rnR5mZO18r5rct2C1qLO3/+frX59X9Y0ACNuADcVAEHbAAVwiCeNgApXAIzsEteA6zCAWRQUyRtUg2Uo1cRF7iSDhNHA23GzeIcqN09CImhOVg7/F++DsEC8JFohHxKpMj0zNSMpmVfIrZm4JS2lnWUdWp31lvsFWxp3L4cdpx2XO78djz6vLJ82sLBAtmCqUIh4p4ijqLOYk7SThKOkm5SwfLJMtul2uUv6UwqcSirKsSobpPbURDQDNcq0V7UddV745B0Ro/Y7zJTtMlcweLPMt6q3brLptLtgN2iw4Ojm3OKi5NbirubZ4WXiM+8X5k/6ZA7yD2EObQ4PAA+utIw6iS6Lex7nH9CU6J99YFMCZSstNE0p9n3sy+sqF2k1fOj82H8r0LRLfOFF0s2bY9otSmjL+8ryJi9/yevCr2fXU1OrV3DkQcQupqjugdHT6W2iDYeKtpU7PNKdUz1uc2tdV1lHb6dfFefHi5+qrfdWL30ZtaPRd6rW4/7M8YUB1Eh2buTwwPPSx7LPuk9umvZ1bPi1/0vWR95TV+8PXkW/V3sRMHJ2+9fz+N/yD0Ue2T5YzPZ9qX8K+usxKzc3Pb54XmG7/pf9v3beG7z/e2BZ4FxkLbwuIP0x/5P3oXqYsei3sWB5eYlkyXMpZOL03+FP/p97P85+2fP3+p/wr/tedX369fy+rL9OW9y/3LywDJEZoaAACAUMwB8M+Xl7/KABDLAX6WLS8v1i0v/zwMgI4CXIlb+Q8JAIDABlAzAwDQ1/EjB/5l/wN+1dn/6wtEmQAAAAlwSFlzAAALEwAACxMBAJqcGAAABVFpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADw/eHBhY2tldCBiZWdpbj0i77u/IiBpZD0iVzVNME1wQ2VoaUh6cmVTek5UY3prYzlkIj8+IDx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IkFkb2JlIFhNUCBDb3JlIDUuNi1jMTQ1IDc5LjE2MzQ5OSwgMjAxOC8wOC8xMy0xNjo0MDoyMiAgICAgICAgIj4gPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4gPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iIHhtbG5zOnhtcD0iaHR0cDovL25zLmFkb2JlLmNvbS94YXAvMS4wLyIgeG1sbnM6ZGM9Imh0dHA6Ly9wdXJsLm9yZy9kYy9lbGVtZW50cy8xLjEvIiB4bWxuczpwaG90b3Nob3A9Imh0dHA6Ly9ucy5hZG9iZS5jb20vcGhvdG9zaG9wLzEuMC8iIHhtbG5zOnhtcE1NPSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvbW0vIiB4bWxuczpzdEV2dD0iaHR0cDovL25zLmFkb2JlLmNvbS94YXAvMS4wL3NUeXBlL1Jlc291cmNlRXZlbnQjIiBleGlmOlBpeGVsWERpbWVuc2lvbj0iNDEzIiBleGlmOlBpeGVsWURpbWVuc2lvbj0iMTMzIiB4bXA6Q3JlYXRlRGF0ZT0iMjAxOS0xMi0xMFQyMDoxMzoyMC0wNTowMCIgeG1wOk1vZGlmeURhdGU9IjIwMTktMTItMTBUMjA6MzI6MTItMDU6MDAiIHhtcDpNZXRhZGF0YURhdGU9IjIwMTktMTItMTBUMjA6MzI6MTItMDU6MDAiIGRjOmZvcm1hdD0iaW1hZ2UvcG5nIiBwaG90b3Nob3A6Q29sb3JNb2RlPSIzIiBwaG90b3Nob3A6SUNDUHJvZmlsZT0iRGlzcGxheSIgeG1wTU06SW5zdGFuY2VJRD0ieG1wLmlpZDo1ZTkyMTdiMi00NWQzLTRmYzEtOTVjZS02YzM4MzEyOTZlNmYiIHhtcE1NOkRvY3VtZW50SUQ9InhtcC5kaWQ6NWU5MjE3YjItNDVkMy00ZmMxLTk1Y2UtNmMzODMxMjk2ZTZmIiB4bXBNTTpPcmlnaW5hbERvY3VtZW50SUQ9InhtcC5kaWQ6NWU5MjE3YjItNDVkMy00ZmMxLTk1Y2UtNmMzODMxMjk2ZTZmIj4gPHhtcE1NOkhpc3Rvcnk+IDxyZGY6U2VxPiA8cmRmOmxpIHN0RXZ0OmFjdGlvbj0ic2F2ZWQiIHN0RXZ0Omluc3RhbmNlSUQ9InhtcC5paWQ6NWU5MjE3YjItNDVkMy00ZmMxLTk1Y2UtNmMzODMxMjk2ZTZmIiBzdEV2dDp3aGVuPSIyMDE5LTEyLTEwVDIwOjMyOjEyLTA1OjAwIiBzdEV2dDpzb2Z0d2FyZUFnZW50PSJBZG9iZSBQaG90b3Nob3AgQ0MgMjAxOSAoTWFjaW50b3NoKSIgc3RFdnQ6Y2hhbmdlZD0iLyIvPiA8L3JkZjpTZXE+IDwveG1wTU06SGlzdG9yeT4gPC9yZGY6RGVzY3JpcHRpb24+IDwvcmRmOlJERj4gPC94OnhtcG1ldGE+IDw/eHBhY2tldCBlbmQ9InIiPz6Nh7q2AAA4j0lEQVR4nO2dd1wTydvAnw1NNAFRA3igAopyioByqByinnI0z17OclasiOUUD1GxoOJhQQQOFEXFiqJiO2wgYsUChuJPFATpvScEQpJ9/5gzby4kITSN3Hw//BFmZneenZ199plnZucBsnWwxGFsbAwAL1++bJw1depUAHj69CmLxTpx4gQAbN26NSIiAgDWr1//8OFDANiwYUNOTs6RI0du376dl5fn6uoKAFu2bGEwGAAAAAwGg8ViLV26FJ3q9evXR44ciY+PT09PHzt2LABERkbOnj0bAEJCQlgs1uTJkwEgKirKz88PALy8vFgsVnR09C+//BISEmJpaQkAcXFxSUlJVlZWAJCQkCAidlhYGADMnDmTxWKdOnUKADw8PJ49ewYAK1asQIIZGxvX1NQwmUxzc3MkJACYmZklJSXFx8cDAJVKZbFYISEhBw4cyMnJEdt08kBxcTEAfPz4Ef1LpVKfPHmycuXKxYsXC8rY2tru3r37xo0bhoaGKAVdI4vFMjY2PnXqFIvFOnfunI6OjuAQLy8vGxubqqoqAHj79i1KjIqKQkcNGzbs3LlzKHHGjBne3t5FRUVLlixBd3zatGmfPn1isVh0Oj06OrpdL19GvLy8UO8VTkQd2MrKqnGiCCJ9LCgoCABcXFxYLFZ2djYA6Ovrs1iskJCQwMDAwMDAsLCwvLw81KMqKioAwNzcHBU4ceJEaWkp6qKWlpYsFktfXx+1Kp1Op9PpqArULXNyclrfLcvKyoQfcA8Pj4kTJ0rPRc+4oFNZWlp6e3ujxtHR0Rk9enR6ejrKSk9PnzFjho6OzuzZsxcvXrxmzRrhqj98+AAADx8+bCyVsF5SbNzircfAwCAlJSUrK2vQoEEoZc+ePbm5ue7u7hkZGQAwePBg1LjoUtetWwcAN2/epFKpAODo6Ni1a9fMzEwPD4+SkhKRk9vb2xsaGgqnaGlpPXjwwNXVlclkihRGGrZv374AwOFwUO0WFhYAMGLEiIsXLwLA2rVr0b+Co3JzcwcMGCDp6hwcHJC0SkpKADBt2jSUbm5uTqFQAMDS0jI+Ph51GgaDYWJiggowmUwmkzlr1ixZ2vCrQxCE8L/du3d//fq14N/s7OwePXpIOpZGo6FD8vLyGhoaUEMVFBTQ6XRUoLCwUE9PDwCKi4sFiXw+H/1AD0ZVVZWHh8e+ffueP3/u4uKyb9++/fv3t931fWn09fXXrFlTV1d38eJFBoMRGBh4+PBhkTJDhgwBgO7duxsZGaWmpvL5/LVr16JebWZmNmHCBBsbm6ioqOvXrwMAelsrKiru3r178eLFMoqBlGYru2WnTp2oVGpGRgZ6wNPS0nR1daXnogLoAQcAGo3W0NCQkpIyZswYb29vZ2dn9OwAwPv37w8fPqyurg4Ajo6OkyZNunTpUnBwMHoF6ujoGBkZpaeno6dYEhTZL0Z2kAUUEBDQ0NCALszLy+v06dN0Ol1TUxMAkEVQVFQEANra2qqqqjNnzkxLSwsKCqJSqRYWFtHR0d7e3ubm5gkJCRcuXBA+edeuXUWqO3z48KVLl1xcXN6/f79mzRrhLBUVFQAQNFn37t0BoLCwEEnl4uJy7dq1nj17AsCTJ08YDEZsbGxERISpqamUq6NSqdOmTWMwGEePHqXT6cOHD0fp6HUKAPn5+QDQq1cvADAyMmIwGAwG4+7du9euXevUqVNtbW1NTY3gGf5W+Pnnn+/cuYNMj7CwsNTU1NGjR0s/ZMiQIVQq9eDBgyRJvn379uTJk/b29ijL29sbtYO/v/+kSZMAQFNT8+7du1wu9+XLlw8ePACAoKCg5cuXA8CoUaOGDBnC5XIBQFlZuaysDP3+tvjuu++WLVu2Zs2akJAQALh9+3bjMqjn8Hi83NxcKpVKoVDCwsIuX758+fLlffv2wWddtnHjRgBwdHSsrq5esGBBbW1tTEzM+/fvxdZbUlJSX1/P4/HQixY9Aq3vlhMmTPD396+srExOTg4LC0OdITQ0NDIyUmzu999/r6Oj4+/vT5Lkmzdv7t27Z2VltWfPntWrVzs5OdXX17PZbKQuQkNDPT09ORzO48ePY2NjHRwcBgwY8Pz584iICB6PFx0dnZqaigwjKbSLXnNycqLT6bGxsSNHjnRyckJ2mbOzs6qq6vjx4wFg3bp1N27ccHd3BwA7OzvUEABQUlIye/ZsBQUFZGl/9913FArl3LlzIPQyF7EjAKC8vBwA+vTpU1JScuXKFeHCIvz8888A4Onpefny5dWrV588eXLQoEFoaHzkyBEGg7FkyZIpU6YoKytLv8CZM2cCQF5e3rx58xQUFFAiGucePXr06tWrOjo6w4YNGzZsWGpq6s2bN588eWJnZ3fw4EFFRcWffvpJW1sbmdPyCWphkXb+4YcffHx8xo8f36VLl7Vr14aGhvbu3VvwwhABpVOp1MuXL584cYJKpQ4bNmzZsmVTpkxBp6VQKHQ6XVtbW1lZeevWrQDg4uJy9uxZdXX1ZcuWIWfCypUrCwsLu3Xr1rNnz+TkZGdnZwCws7ObOXNmYGBgO7eBrCBljdi1axdKTExMFCROmTJF5BA0FBC8BYXx9va+dOmSm5sbk8n85ZdfAOCnn35ycHBwcHBAHhL0YigpKdHX1zcyMkKDeg0NDQ0NjTNnzjQ+Ye/evQFg8+bN69evRyl0Or1NuuX27dtZLJaOjs6IESMWL17s6OgIAMHBwbdu3RKbSxBEeHj4xYsXqVTqyJEjd+7caWFh8eDBA39//x6fWbVqFQD88ccfCQkJGhoa06dPDwwM7N27t6mp6Z49e3777Tc1NbWJEyf6+PigcZg02sO/xmKxUlNTbWxsBLWsWLGipKSExWJVVVWtXr1akH7w4EFUHtlQAHDjxg0Wi5WTk2NkZIRSUC+fP39+YmIiAMyePRsdgvxrz549i4mJEdi3yFQ8fvw48q8hTxx6v0VGRrJYrKNHj6LCOjo6e/fuZbFYRUVFqDAAGBkZXblypfHlCPvXWCxWaWkpKo8chQL/GhpS0en0e/fusVispKSkYcOGoZLW1tZJSUmsz87H+Ph4SU0nz9TU1GRlZTGZzGYdlZOTU1VVJZJYWFhYWFgonFJRUYGcaMLk5uZmZ2e3QNT2BvnXhJk6dWpsbKxIIpVKRYnCTjfUt9+/fy9IQf41gQlsbm6elpYmtl7Uwzdu3MhisZhMJnrFCno+8qYJ/Guo36KOjcyLsrKyNuyWubm5paWlzcrNzc2trq5u8sx5eXkixaqqqjIzM6UcK6yXCJIkpeq9JqitrZWSy2az8/PzdXV10XhQAIfDKS0t1dLSEhg7jeHz+UVFRTQaTaCzpFBfX19eXt69e/cmTS105tLSUjQiFhaVxWJJ8RkJSExMjIuLW79+vaGhIdJoaWlpZmZms2fPDg4OLikp0dTUFDZ2KisrFRUVZbkKzH+W06dPr1y5Migo6Ndff2UymWi0KCOlpaUKCgoaGhpic2traysrK7/77juR9I7XLTt37iz43S7zBgJUVVWRz14EZWXlxg0tAoVCQZ4vWVBRUZG9MIVCEVFqAKCqqqqqqirL4bdu3fLy8qJSqY1HQxQKRUtLSySxsUMQg5GEioqKiBHQJNJfxp07dxZ+4AV07G7ZvvZah6Sqqqq8vLxXr16Kiv+8FTgcTlZWFpVKlV23YjDCVFVVFRcXa2pqonlATAsQVt9Yr2EwmI6AsF5rl/lQDAaD+YpgvYbBYDoaWK9hMJiOBtZrGAymo4H1GgaD6Wi0dj4U0yT8BpJbTwKAggqhoCT6ERgGg2lz2nddLobD5Jem1hcl1ZMkqTmoE32giooatpExmPYF67V2hM8lKz81ZD2qrcjgAAC7lKegQmgNVqEoYqsNg2lHsO3QjpB8YFfw2BU8kg8kH+qq+LWlXP63t8UOBvONgfVaO0JQoJM6pZP6P9/2K9MoqhoKFIlf+mMwmLYBj0PbEYoCoaGv3NuaVFGjkHzQHKTSrZ8yHoRiMO0Nng9tZ0jgNZANbD6QoNiJoqhCAFZrGEw7g/UaBoPpaGD/GgaD6WjIvV4jScAWJQaDaQ7yOm9AksBmQ2UlVFcDAKipQdeu0KkTSAgUgsFgMALkUq/xeFBUBK9ewcuXgML26OiAhQVYWICWFijKpcySENibWCNjMF+K1uqItt8vlySJwkKFyEjFS5eIpCSCyQQAkkolBw/mTp/OGz+e/O47aBRqTx6pqyMKCigfPxLl5aSGBl9fn9TRAdlCKGAwHQmxARbaFfmzfdhsyqtXimFhlJcv4XP4W6K6moiLU+TzSU1Nnq0tfPFmai5ETQ3l3j2lY8cIBoNoaAAFBf7AgdylS3kODmS3bl9bOgymgyN3eo2orFR4/pySlAQiMb15PEpyssKzZ3wLC1LO9VpdHeXePWV3d6KwED5HaKYkJCht3gwkyZ0yBbp0+boCYjAdGzlz+pAkUVlJpKYCiyUml8Ui3r0jKipAQjh3OYEoKFA6cuQfpYb8ayQJPB5RVqZ45AglOxvP8GIw7Yqc6TUAqK4m8vPFP/kkSRQWQk3NF5epOZAk5cMHSlLSP0pNJCs5mXj7VtQUxWAwbYqc6TWCIBoaiPp6ifl1dURDw5eUqNnw+UReHtTViVfNPB4lJwd4vC8uFgbzH0LO9Jr0VbhoGlS+B6EAAHV1EjUXWpeHx6EYTHsiZ3oNgCQIiZoLab1vYpGHFCEJAus1DKZdkTu9BhSKdM1FfhMLXKVoLpLEa3TlltLS0lGjRr18+bJZWWIpKyt79+7du3fvMjMz21RGuYPJZL77zNeW5R++tQeMQvkGrDUeT6JqJgjgchtrPW9v7/79+9vY2LS4zqKiokOHDt2/f18kvbKyskuXLo6Oji0+87fOq1evughhbGy8f/9+ngRHAYfDiY+PZ4mbjpeSJZbQ0NAffvjhhx9+WLx4MZfL9fLy0tPT09LSWr16NZPJFC4pNrempmbt2rUocd26dSjRz89vuBBPnz5tVlNIF4PJZG7YsEFPT8/Y2HjLli1cLreysnL4v3F2dm4sxokTJ374TLPkaT/kbP0aGoRKMXak58oHhHRjrVHuhw8fPD09zc3N165d2+JKCwsLt27dumTJkp9//rnFJ+nA0On0pUuXstnsa9eu7dixo1u3bk5OTl+gXg8PDzs7u8DAwD179jg6Onbr1u3EiRMcDufo0aOCMmJzV6xYce3atWnTphUXFx87dozD4QQGBkZHR3/69GnWrFnoQHV19WYJI12MdevWXbhwYdasWQRB+Pr6crlcd3d3DQ0NlFtRUZGSktK7d28AEBHDwsLizp07Tk5OeeirRzlAzvQaSRLSPWgUinzOG4SHhx84cCAlJWWEhcVpCwtDySXzCgoWjR//+u1bU1PTffv2sdnsBQsWAEB2dvbTp08nTJggS3VBQUGHDh3Ky8szNDTctWuXvr7+3LlzASAsLKywsPDixYvh4eF79uyprKycN29e48Pr6up27doVHh4OALNmzdq2bZuiouKECRO6d+/ev3//w4cPnz17dseOHePHj09OTn7w4EFRUVGLWqW11NbWiqjpgwcPjhgxogWn6t+//5YtWwBg9OjRkydPjo2N1dDQqK+vnzx5sqq4j9uCgoLOnj0r+HfQoEE7duxoQb2DBw8eMmTIokWLqFTquXPnlJWV3717d/bsWV9fX0G9J06cEMndt2/ftWvXLC0tT58+XVtbS6fTIyMjASAuLs7W1tbZ2ZnD4Xz//feKzfxWunFFAjG4XO6FCxcMDQ1DQkIA4N27dwEBAdu2bbtz5w461tXV9dOnT/v27ZMkBo1Ga0H7tBNyNw4lKZQmnFPyx6NHjxYuXFhUVDR37txXr14dCwmRWJQgTp458/Hjx6lTpyYmJtrZ2SkqKg4fPhwAevfubWJiIkt1Hz58cHV11dbW9vHxqaysRK9N9LT37NnT0tLy7du3CxcuLCgomDBhwpkzZxqfwd3d3dfX19zcvHfv3gcPHty/fz8AJCcnI20IADwej8Fg7Nmz58aNG3Q6vSWN0kZoa2tra2sDAIPBYDAYDa1b5VNbWxsdHQ0AdDr906dPS5Ys0dPT27t3b0FBgUhJKpWKqkb1Zmdnt7hSkiTT0tIGDBigrKwMAMbGxgBQWFgoJbegoMDHx8fd3R0A0tLSAKB///4lJSVMJvPq1atDhw4dMWLE2LFjKysr20oMRG1tLZ/PBwA03M7Pz0fpb9++DQoK2r9/v76+fivF+DLInV4Dgmhi3kBB7gKfHD9+HADCw8ODg4P9/Pz69esnUf+SJAVg1+7drq6uy5YtYzKZKSkpixcvBgBbW9s5c+bIUh2ynsrLy+l0emho6LVr1/T09FavXg0Ao0ePXrdu3dWrVwHA19fX39//0qVLIoezWKzg4GAzM7Pdu3cHBASgkoLcEydOFBUV9e3bF/3LYDBSUlKa1RptSOfOna9cuXL+/HlkU8ycOXPkyJEtO9XTp0+7dOlCp9P9/f0BYN68eWvXro2IiBgxYsTu3bv79eu3bNmy5ORkQfl58+ZduXJlyZIl6F9vb+8WX0VdXR0AIG0CAOhaqqqqpOTW1dUtX7583Lhxjx8/trW1BQA3N7fU1FQAsLKyCgsLs7W1jY+PP3bsWFuJoaioaGNjk5eXt2DBAicnJ6RMBc7EjRs3UqnU6dOnA0ArxfgyyNk4tEnk0l57//49AAwdOhQAFi5cqMhiwebNkgpTAFauWMH+/G9WVpaenl6zqhs+fPisWbPCwsLQGHPq1KkiT3tOTg4AmJqaAsCgQYNEDkemB4PBEJiHTCZT4EKeMWOGoKS9vb2hoZQh9RfC1dX1+fPnZmZmf/31F9HSVT5UKhWZtP369ZszZ46ZmRkA2Nra2travnv3bvXq1efOnauqqjp06JDgkPfv36Mn+cyZM6ampgLjpbkgDVLz+TsZ9KNHjx5N5vr5+bm7u1Op1CtXrowbNw4AmEwmaoGBAweamJjIPjnbpBgA4O/vv3jx4qtXr+rr65ubm8fHx3fv3h0A0tLSYmNjFy9ejLblsLa2bo0YXwY5s9dkWNtFyJ9/DQ2UiouLAeD27dsRV69KKawAcOjQIQaD8eLFi4iICOQXaxZlZWXLly9nMBgnT57U0dG5evWqYBoUDSLQyDE3Nxc+G3fCoM5qZGSERlh37969du1ap06dAIBKpVKE1qB07dq1ubK1OcHBwcgrdPHixdZsd2Nqanr9+vXr168fPHjQ3NwcJfJ4vDt37vzxxx/Pnz8HAAsLC0H5ysrKadOmAcDGjRunTp3auosAIyOjlJSUmpoakiSRFtDU1GSxWDU1NXw+X2xuQECAu7u7kZHRy5cv7e3tASA0NNTBweHevXsAgN5DVCq1DcV4/fq1s7NzSUlJUlJSeXk5lUrt2bMnALx+/RoARo0ahU7SejG+AHJmr6F5AynI5aJWR0fHqKio9evXT5o0aee2bYv+PX0uAg8gLCysk4bG5cuXb926dffu3eZWl5iYOG3atAULFmzYsGHEiBFXrlxRV1dHvtuYmJizZ8/a2Nj4+vq6u7uz2ezTp0+LHK6pqTls2LCXL1/evHlTQ0PDxcXF2toauedFPOgtNo7aiv/973+///47ANja2t6+fRsAzMzMhLVPa3j48KGLi0tmZqa+vv7evXvnzp3bvXt3gVH2+++/o3VnPXv2ROMspFxaxvz58zdv3rxw4UItLa20tLTffvtNWVnZ0tIyNTWVwWA0zmWz2W5ubgCgpaWFhsAUCuWXX355/PhxRkbG8uXLz58/DwAyzjLJKMaVK1euXbu2bdu2kpKSzMxMNzc31Kni4uIAYPDgwegkWlparRTjCyBneg2ABKmDTT5fDr83cHJyysjICAgIuHHjhsXQodPs7eHkSUlXMWbcuD+jox++eEGn0z09PUeOHBkTE9Os6mxsbJYsWXL8+PHQ0FAAWLBggZWVFZfLtbS0fP78+Z9//pmcnLxhw4aDBw/OmTMHeWcU/u2UPH78+JIlSzw8PADA2tr6r7/+Eqniq2s0hMCdf/XqVeQ03Lp1a3P1mqRref369aBBgw4fPjxmzBiFRk7bjx8/oh/r169HP/73v/81q15hnJyc4uPjr1y5AgCWlpYiU6uNcwULXGNjY2NjY9FvNEHp6em5bds2APDw8GiuISldjD179qSnp3t6egLA/Pnz0awFADx58gQ+D0oAwN7evpVifAFaG2evjffLJUnKs2cqzs5Eerr4/L596//6i29lJYdL9uvr62tqanpoaCgePars5iZ+PQpBNLi7s9euLWOxtLS0WlMdi8UqKyvr1q2b8CiAzWZTKBQVFRX0m8ViCTtQRKisrFRUVJTDQcSXgcvlNl4nkZ+fb2hoeOvWrZ9++kn2LLH4+Ph4eHi4u7vb2Ngg115NTU1DQ0M3CRuLSs8V0NDQUFRUpKmpKZgBaC5SKiJJsri4WF1dHfklZBSjtLT0zZs3q1atysvLE7tuGe+X29TMgPwNQgWoqKioqKgAjwc8njQ5eTxFBYVWKjUAQKvnRRKFB5KqqqpiV2YJkAf32VekuYu/WsbevXujo6ORSS59hZeM67+UlJR0dXVbI5KUigiCkLFnCouRmJg4efLk1ojU5siZXiMIQvqWHnx+Ew44OeCfpcUStpCTz3XFGASVSt24cSNaVS97llgWLlyIHE9f3lr5wlhaWjIYDJAb9wXInV5rMlqo9FW78gEp5aMIggD5W3+HEaCmpibpuwIpWWLp1q1bk4PKjkHnzp3lYT2QMHLnpSKlq3y5V2oAUnfsIElQVv42rgKD+WaRO70GFIq0OQGCkPd9ikiSkOJfQ19TYJMNg2lP5FJHSDVn5GUELwkuF6qrpV2CoqIcTuZiMB0JOXvACIKQvhOR3I/gCA6HkPIZMLJG5ca9isF0SORMrzW50zdByHvQk9paqK6WmEuSJJWK7TUMpl1p7XxoG89hkyRIXXJFAKh06gSdO8u1apCypaqKinKPHsqqqtjFhsG0H/KnHaTsU4SyvshyyhZCklBbC+XlEguoqoKqKh6HYjDtilzqNUlIWe8qJxAEMJnSxqHKyqCkhPUaBtOuyJlea9J9xufL9Xp9koS6Oqirk1gA2WsYDKY9kTO91uT3BlKii8oJ9fUgZS8AVVVQUfmC0mAw/0XkTK9BUzusyX9cZA4HOBzxWQQBXboAlSrXQ2kM5ttH/vSagkIT6zzkeSYRjUPZbIm56urQqZO8q2YM5htH/vTaNw1JQnW1RHsNAKhUUFbGeg2DaVfkTK9JCIf+rwJf1r/2/v37kf/m1q1bwgXevn3r6OiopaU1cuTIRw8eQFWVFPnD79yxmzz54sWLwokHDhzYvn17G8p8+vTpxnGbNm3aJHwVEydOBIDCwsK5c+dqaWmZm5sLpKqtrV2/fv2AAQMcHR3R7tstpqamRktLq6KionFWVFSUo6Nj37593dzcUKC29PR0YQnT09Ph38378OHDZtV+4MCBAQMG9O3bd8OGDXWNJnMa5+bn5wta4+bNmyDD3W8WUloDIXzjysrKnJycevXqZWpqKtjMPS4ubty4cVpaWvb29oI9dWVHuKeJ7Q/SBRbbW65cuSI4CdqcWS4g2wF05vr6ekkFfHx8jh49+s8/fP4/f+j3gwdk374kgPg/Q0MyOprk8WQX5uXLl97e3ikpKS27FiaT+fAzKJ6eyKkMDAw2btyYlpa2e/fuLgB1mzeTFIpY4fkA5QsXnvL1BYAPHz6QJBkeHu7s7AwAv//+e8vEa0xmZiaNRhs4cKBIenJysuBCrKysnJycSJIcP3784sWLU1NTg4ODBVLNnDlz5syZiYmJZ86codFoZWVlLRAjNzfXy8sLhbwqKSkRyWUwGDQaLTw8PDExcezYsdu3bydJMjw8fMyYMa8+w2KxyH83LwBUVVXJKEBYWBiNRouJiUlISDAxMfHw8Ggy18TEZMWKFR8+fEDBCRMTE5u8+23SGgiRG+fg4DBjxozU1NQbN24AwKNHj1CYseDg4LS0NA8PD11dXdREstC4p4ntD9IFFttbXFxcNm7ciG5ZfHx8CxqnPWhHvVZXVyelgIGBAcnjkUwmmZ9PfvxIfvpElpWRdXVkVBTZp49EvTZgABkTQ3K5/+hBGfDz8wOACxcutPKiqqqqDAwM7t69K5zI5/NpNNr9+/dJkkxISOgKwF6xgiQI8cJTKOTOnZyqKgB4/vw5SZK3bt0KCQmxs7NrK73G5XLHjBnj5eXVWK8JOHLkiLW1NYfD4fP5Y8eOzczMJEmSw+EgqfLy8gCgvLycy+WSJPnp0yfZHx5h8vLyQkJCULDOxk/y6tWr3dzcSJJsaGioqanJyckhSXL79u07duzIzs4uLS1FxUSaFwBkV7IzZszYvHkz+n3+/HkDAwPpuejC8/LyUOL48eN3794tKC/27suO9NYgG904FP757du3KNfa2hqFghVcBQo89ujRIxkFkNLTBP1BusBiewtJklZWVnfv3k1PT29ZP2kn2levbd26dejQoX5+frq6ujQa7cCBAyRJorgbCgAOxsYVJ09Wb9hw18DgmJKSn7Z29MaNZFgYu3t3SXrtk5pa8YULZE5O4atXi375RZNGMzUxOXPmDEmSmZmZQ4cOXbdunbW1NY1GmzBhQn5+/smTJzU1NQFAV1d3z549rbmo3bt329nZiabyeL4+Pr11dXfu3Dlw4MDNixeTc+dKVMpKSglOTlMmTrSysmpoaBCcw83Nra30mo+Pz5o1a54+fSpJr5WXl9NotKdPnwonhoeHT5o0CUkVGxtrYGCAoizTaLS//vqrNfKgAWbjJ3nMmDHr16/v378/AIwdOzY3N5ckyUmTJgm2qHZycuLxeCRJHjp0SPdz87q4uMhetZWVVXBwMPqNIo9Iz0WiJiQkoMShQ4euWLFCUF783W8mklqDFHfjysvLUQsghXv//v3IyEgajYaGQeXl5S14WzfuaWL7g3SBhXsLj8cDoY3Fjx071ix52o/29a9lZWUlJCSsWbOmX79+NTU1rq6u+fn5Y8eOpQDoAqzq0YMaEKB06NDPGRlODQ3Lior67N9f5utLSJpPBKirrr7v6so/fPiOvT391q3NlpaKJSUL582LiYmpr69PSEjw9fXt1KnTjz/+ePPmzTlz5vTu3XvgwIEAMGjQoO+//77FF1JdXb1169ZNmzb98z9JQk0NpKRARIR2ZOQPubmF9+9XZmfXFhXxSkslnkVRMe5//0tOSencuXN5efk/a4zbbs3H27dvAwIC0HhNEoGBgWZmZj/++KNwYkRERHJyMpKqvLw8IyOjpKQkNTXV19d31apVb968aSsJBRQXF/v4+Pj4+CQnJ9fX1yM1qqio6OXlxWazHz9+HBISgsImZWVl5ebmxsXF5eTklJaWsiX3DRE4HI4ggkHjUAaNc9XV1ceOHbt8+fKoqKgtW7YkJCQIQpCI3n0ZYLFYggF1uZTv6gBAwo3T0NCgUCipqak//fTThAkTxo0bZ2lpCQArV66MiYlZsGABAMjeGpIQ2x+kI9xbSktLR4wYERERweFw/P39ly5diqJuf33aQ1miM9fV1aGA5MieQlFmnzx5QpIkDcCdTicHDeKLc0XxJPinSIKoU1ZmKyryCYIPUK+kxDY1TVy0qDfA3DlzUlNTAcDExITH4/H5fGQVZmdnt8k49MKFC5qammhoRvL5ZF4e6etLjhjB7dGjAqBBQ4M0N6/esmWaqmr5wIES7TUNDfLiRQ6bPWbo0JMbNpAREeT58+T9+3+uWOG6Zk3rm33o0KGenp7x8fHHjx/X1dUVmB7CGBgYiH2pcjgcCwuL/fv3R0ZGAkB+fj5Kt7KyOnToUItFkmShWFhYCIyv6OhoAOD922fq5OS0bt26rKwsAGAwGCRJFhQU0Gi0yMhIGau2srIKDAxEv9GEQ5O5eXl5s2bN0tXVnTdv3rJly9avX48K/Ovuy0ZiYqLgERP0PUmtIenGnThxAgB+//13gUvn+fPnY8aMMTAw8PT0NDExuX79uuwikeLsNUn9QbrAgt4ikq6rq3v16tVmidROfIlvyJH3Ee2AXl9fDwBaAKNqayEjQ2zwdoqkGU+SVBFaQqHc0MBPTFRKTPwZIP39e5RoYWGBIpZbWVm9evUKuVpbz+XLl2fOnPlPlMmyMjhzBnx9obhYgSS7AkBFBcTH01JT/1BVBcn2Gl9BgaKkpPT2rQeFYhwcDH5+wOMBlTqle/cngwdDTQ3IFpFIErm5uQEBAQEBAWw2u6amxt7eXiTYe1JSUkZGhqOjI/o3IyPDzMwsPT1dU1NTSUlp5MiR7969Q9F/BcH31NTUOFKWrbQUfX19weAFbQlTWVm5f/9+Dw8P9K+CggKHw/nw4QMAIItbW1vb3Nz8w4cPDg4OslTRp0+ftLQ09Ds9PR2NeaXnvnv3LigoCIXpGjdunCAs5r/uvmyYmJiQMlviYm/cqVOn1q5dGxUVNW7cOFSsurqaxWI9ePCAIIja2tpt27bp6enJLlJjRPqDdMT2lsTExEePHiFzG8HlclsjUpvRHsoSnVlgr6WmppIkuXnzZgCIjo4meTwrgCRlZYn+9eb88Tt1KrK0jL94EdlrAifIjBkzACAlJQXZa+fPn2/NFWlqaiKrk2xoIO/f5+jpiRWGB8Dr1EmSqJU0GmfVqrohQzgA/3/tBMGjUAq6dyevXyclzyA3C2E3zf379wMCAtDv06dP6+rqCorV1dXRaDQ3N7eamprMzExNTc3AwEA+n4/8WXw+//Xr1/DZPdwyRF74AmFOnz6tqamZnZ3d0NAwZ84cOzs7Ho9Ho9H27t3L5/NTU1NpNNrly5fRa8nPz6++vh4J8/DhQxmrDgkJodFoaWlpxcXFAwcORHaKQACxuXPmzHFxcamvr0dh8T59+oRO9f93v3VIag0BghtXX1+PWqD2Mw0NDfX19ZqamufPn+dyuV5eXgYGBnyZZ88QIvaaSH9oLJKwwGJ7S1JSEgDExsaSJIkmkbOzs5vXKO3D19FrDp07v2m1RhMsnkgECF2yJPVzOO4DBw4EBAQAgK6uLofDCQwMBICff/45KiqqZZfDZDIBIC4ujiRJsrqa3LOHr6DQAlErCaIEgC8uq15RkZwyhczKapP2F9Zrbm5ugr67detWBwcH4ZJ3795F8yoA4OTkhCbFEhISBNaNl5dXayQReZIFwvB4PBcXF1SFhYVFWloaSZJRUVECYdzd3dFDe/LkSYFlt3PnTtmr5nA4M2fORAdaWVkVFxcLCyA2NyUlBQUwptFox48fR+f5191vHZJaQ4DgxgkPYxEnT54kSTI8PNzAwAD17VevXjVXgE2bNgnrtcb9QUQkEYHF9paDBw+iFkMLd5orUjvxNfQan79v/vwogIY20mu5XbvW/f030msmJiao6TU1NdEseGpqKkoRnt5qOSUl5NKlLZRWin1KEGTfvuTDh7KvX2kr+Hx+UVFR40U5ZWVlzfIotQAOh1NZWdlYGOHJYpIkuVxuYWGhyEIEGWEymRUVFc3KFaxxkVtatqKwTRDbWxoaGtCLQX5oF73WNKWlvM2b+WpqbWCyKSqSo0aRb96kvnsHAPPmzePxeAUFBcImOp/PZ7FYIk9LCykuJp2d22QELfrXqxd5/Top308UBvNN8JW+o+ralTJlCjFyZGs3vyUI6NYNfvwRdHQEH11SKBRtbW3h0NMEQXTu3LnxZH9LUFUFQ0NQUmqDUwmD9gFWU8OfjmIwrecr7amtoACDB4OrK1RVwbNnLV/DpaYGo0fD1KmgoaGvro78zW0qaCNUVcHSEoyMICmpLU+rrAzGxqCvj/UaBtN6FHbs2PF1alZUBF1dsLCA0lJ4/77Zqo1CgZ49YcoUWLkSTE1BWVlBQaFHjx7trtcoFFBXB4KAV68k7kfUXAgCevcGNzcYMkSuozdgMN8IX0+vAYCCAmhpwejRoKAAyclQXy/TUQQBGhrw88+wfj0sXAj9+4OycjsL+m+UlUFfH/h8SEyUVWYpEAQYGsKOHeDgAG0b3AuD+a9CkPKwd2tlJZw+DV5e8O91pGJQVQUbG5gzB0aNAjq97f1cspOfD0eOgI+PtKh60iEIoNFg3DhwdoaRI6FTpzaVD4P57yIfeg0Aqqvh6lXYtQsyMiSW0dYGd3eYNAm+++5rajQEnw9FRXDoEPz1l7SABpJQUoI+feDXX2HuXOjX7+tfDgbTgZAbvQYATCbcvg3e3vDmjejmkQQBBgbg4QFTpoCa2leSrxEkCXl54O8P/v7N87VRqTByJCxdCqNHQ9eucr2zOQbzDSJPeg0Aqqvh5Us4eRLu3IGKCiBJUFCArl1hxAhYsQLGjpU7DxRJQk4OHDsGhw8Dk9n07IeCAhgYwPz5MGMGGBhgMw2DaQ9aq9dqWzAEk05dHZGfT3n3jvL+PVRUgLo6v39/vokJqaMjp1qAJInCQoXr15UOHyZyciSqNgqF7NaNP2oUd84c/ogRZNeueEkH5j9C5y9ujsifXgMAPh84HILDAT6fVFICZWVQVJR3LVBbS0lOVgwJUbh+nWAy/5VFoZCdO5ODB/N++YVnb8/X08NTBJj/FFivfdsQlZWUpCSFv/+mvHhBFBX9o9F69eJbWfHGjCENDUkqFShyFisHg2lnsF779iFJYLOJsjKiuJiorSVVVEg6Heh0UlUVzw9g/ptgvdaBEEzpYgMN89/my+s1/NVOu4HVGQbzlcDPHgaD6WhgvYbBYDoaWK81zf3790eNGqWnpzdv3rzc3FyxZQ4fPiwSJ43JZOrp6aGdlAXExcXZ2toK/k1JSZkyZYqenp6jo+PLly9bL+rhw4fNzMyMjY3d3d3r6upkyT169OioUaOsrKzQ5ulpaWk2/+b27dstlkdsIyBiYmKmTJlibGy8devWqqoq4SwXFxdPT0/0++XLl46Ojnp6epMmTXr/OTqPjDS3NQoKChYtWqSnp2dlZYXicrVSABmF4XK53t7eKHfbtm0igXKEWyMhIcHe3l5PT2/ChAkMBqPFwjQpktjc8vLylStX9u/ff/jw4efPnxcpLD3A45emlftSshqxYsUKAHj69GnjrOZy/vz53bt3FxYWslgsOp1Op9ObPARFvWp91QLevn0LAEePHmUwGLNnz7a0tBQpcPbs2WXLlgHA6tWrUUpaWtrOnTuNjY0BIDs7GyU+fPjQzc2NSqUaGRmhlIKCAiqV6unpmZSUtGvXLiqViq60xYSGhlKp1Dt37jx79szY2HjTpk1N5vr5+RkZGUVFRUVFRdHp9L///ru4uPjOZ1BoiNevX7dAGLGNICAuLo5KpZ49e/bFixejR4/evHmzIOvChQsA4OLiwmKxUDyqgICApKSkTZs26ejolJSUtF9rGBsbL1myJDEx8cyZMwDw4sWL1ggguzBeXl6GhoaPHj169OiRkZHR9u3bxbZGUVERlUr18fF58+bNsmXL6HQ6k8lsgTCyiCQ219bWdurUqQwGIzw8HADu3bvHEtf/G9NKJdMC2t5eQyGg24Tz589v3bq1urpa9kPWrl3r7e3dVgIAwP37983NzX/77TdDQ8Pt27c/f/68oKBAuICqquqQIUNsbGwEKQRBaGpqLlq0SLiYkpISsvgEKejlv379+r59+65atYrJZGZI+eZfBq5fv75y5Upra2tTU9MNGzZcvHixyVx/f/+9e/daWloOGzbs1atXJiYmXbp0sba2tra2NjMz279//40bN1oWT1psIwgIDQ1dtmzZlClTjIyMLl26tHDhQpReVFS0bdu2NWvWoH9fvHihr6+/aNGivn37bt26NS8vT/Ygzc1tjYKCgpSUlE2bNvXr12/q1Kn29vaRkZGtEUB2YR4+fLhmzRpzc3Nzc/Pp06cLqhBpjby8vJEjRy5fvrx///6//vprSUlJQ0NDC4SRRaTGuVwu9969e1u2bDE0NHR0dLSyskpOTgZx/V8eaN9xaFxcnI2NjZaWlq2tLTKbPT09raysgoKC+vfvr6WlhYLgAcDly5fNzMz09PS8vLysrKy2b9++cePGW7duAcDEiRPv378PAGw228vLS0tLq3///iiWbWPu3buH4oRLqqi5FBcX9+3bF/3+7rvvUIpwAXt7+/nz56MYqYiePXvOnz9/zpw5wsXMzMzmz58vPAgdOnRoWloa2q/8wYMHAKCvr98yIREFBQV9+vRBv3v16pWZmSk9t6GhIS0tLSEhQU9PT01Nbffu3V26dBGUDwoK6tu3ryB4ZXMR2wgCUlJSuFyumZmZurr6zJkzUSOQJOni4rJz5046nY6KUanUkpISNC5DY9X8/HwZBWhua6BrF9zcwsLCvLy81ggguzBnz5797bffAIDL5d6+fXvQoEEgrjUGDBhw5cqVmpqaM2fOrFq1auPGjcqt2Hmwue2jqKiYl5eHApUVFBQ8ffp0wIABIK7/ywPtqNeysrLGjRuXnp4+derUxMREOzu7ysrK7OxsBoPh6upqYGDAZDLd3d0LCgrS0tIWLFiQlpY2adKk4OBgBoORnZ09cOBA9JxbWlr27NkTAJhMZnBwcL9+/fLy8tatWye20uTkZOSoEluRjJInJCTEx8fHx8dnZWU1NDQoff4uFUXGbavIrwoKCmpqagBw+/bt6dOn//nnn2qt26qEw+EIYjg0DubQOBc9qGFhYdeuXbtz505YWNipU6dQgZqamp07d7q6uspeO4vFiv9MRUWF9MIlJSV+fn5//vnnq1evOBzOhg0bAODs2bOKioqTJk0SFBs+fDgArF279tGjR0uXLgWAxm4gSTS3NdTU1EaPHr169eqYmJgdO3YwGAwWi9ViAURaQ7owXbp0UVZWrqysnDdvXmFhIfLkNG4NRFlZ2a1bt3Jzc+vr61vTFZvbPgDQtWtXCoWCQlM7OjqOGTOmxbW3N+2o15CT4s8//3R1dV22bBmTyYyIiEBZISEhd+7cmTx5MgB8+vQJ2WWBgYGHDx8ODQ1FZRYtWoReXO7u7shNAwAPHz588uQJlUpNS0uT5aaKVCSj5A4ODqNGjRo1apSvr6+ysjL78x5EyOxvzUtSBA6H4+rqOn369CNHjggHzW4ZysrKgqeu8ePXOFdFRQUAtm/fbmZmZm1tjaKLowJ37tyh0+k//vij7LVnZmaO+ozgPJLo0qXLihUr7O3tBw4c6OHhcePGjby8vBUrVsybN4/BYOTk5BQUFHz8+FFdXf3mzZuZmZnOzs7m5ubGxsbdunWTUZ7mtgYAhISEGBgYLF++PDc3d/HixXQ6vcUCiLSGdGEA4NWrV8OGDauqqoqNje3Ro0dBQUHj1kAl9fT0Ll68mJCQ4Ofn17JBsZQWaDL3zJkzQ4YMsbe3P3v2LCHHn2y347pc5HN1cnISpGRlZaEfSE+h8R2HwykqKgIA5PIXRORtDJVK1dPTA4ABAwbEx8dzudwmQ0yJVCSj5EVC2/aePHkSDWzhs2bU0tKS8TxN4uLikpiYmJKS0soRKKJXr16C3p+RkYHaU0oulUqlUqmCsWfnzp0FTXTt2rVp06YpNOfDL2NjY5bMWwfr6elRqVT0W1VVFQBKS0vpdLqzszMAlJSUAACFQvH396+trb19+zZBELW1tbt27RIMjpqkua0BAO/fvz98+LC6ujoAODo6Tpo0qaampmUCiLRGZGSkFGFSUlLGjBnj7e3t7OxMoVAAoKysrHFrODo6BgcHo3eGjo6OkZFRenq6hYWFjA0iQgva5+zZs66urn///bc8W2qIdrTX0PMfEBDAYDBevHgRERExd+5clIUsBcrnFfna2toAEB8fDwCNX0H8zx8koQegWYhU1AJGjhyZmpp6/fp1Npt96NAhIyMjTU3N9PT0/fv3i12+IDvJycnnzp27dOmStrY2m81ms9lk675p++mnn06ePPnx48fS0lI/Pz97e3sAiImJOXr0qNhcgiCmTp0aFBTEZrPLyspOnTol6K9Pnz5t8QMjBYEwDg4OoaGhubm5XC43KCjIxsbG1NT002d27drl4uJy6tQpFRWVhQsXhoeH83i8wMBAfX19ZMK3R2sAQGhoqKenJ4fDefz4cWxsrIODQ2sEkF2YPXv2rF692snJqb6+ns1mNzQ0GBsbN26NAQMGPH/+PCIigsfjRUdHp6amDh48uAXCtKx9kLvg6NGjw4cPR921rRwy7UE76jXUUufOnWMwGLt27ZoyZUqRhPAF48ePB4BNmzZNmzZtxowZgnQ04vP19UWm31fB0NBw7969c+bM6dGjR2RkJLrrGRkZO3bsEF511QKb/Pnz5wAwcODAHp9p5Xzo7Nmz7ezsTExM+vTpo6GhgbxjDx48OHjwoKTcvXv31tfX9+jRo3fv3kOGDEET9iwWq6SkRDBb0oYIhPn111+nTZs2YMAAdXX1jx8/+vj4iC2vrKx86NAhT09PNTW14ODg06dPy97OLWiNP/74IyEhQUNDY/r06YGBgb17926NADIKQ5LkgwcP/P39Bd1g1apVYk9iamq6Z8+e3377TU1NbeLEiT4+PgL/TNuKJDb3w4cPTCZz7ty5AjnDwsIEZ5O3MWnbf/e+bt26Y8eOPXv2zNTU1NfXd8uWLQBAp9NXr169YcOGJUuWXLhwgcFgGBoa7tixY//+/ZGRkaNHj46JiQkODuZwOHZ2dr///vusWbNCQkKuXr2KFkacOnVq48aN8HkkOGrUqPj4+LKysk6NNjIzNjbOzMxksViSKmrZZTY0NFRVVfXo0aNlh39JWCwWl8tFgykZc5lMpqKiYuPGbG8aGhrYbLYssyUVFRUaGhotqKIFrVFZWUmj0UTG4C0WQHZhZITL5ZaXl3fv3r1ZXoKWidQmAkOH3M+Dy+WWlZVJ90nFxcVt3brVyspq06ZNx44dc3d337Ztm5ubGwDweLy6ujpVVdXGY8msrCyRKc6ePXvK7n/BYDBfhg6o12Q8ydSpUx8/foz+HTZs2JUrV5qceHry5ElCQoJwypAhQ6ytrVsvDwaDaUP+o3oNAEiSTEpKKi4u1tfX19fXbxMbG4PByAP/Xb2GwWA6Kl9er+H9PDAYTEcD6zUMBtPRwHoNg8F0NLBew2AwHQ2s1zAYTEejtfOhGAwGI29gew2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0sF7DYDAdDazXMBhMRwPrNQwG09HAeg2DwXQ0/g+vaD8F5R2o7gAAAABJRU5ErkJggg==\">\n", "\n", "In this case, `length_in` has a coefficient of around **-0.1**. With linear regression, we could use the coefficient to put things into words easily, saying something like \"every 1000 miles we get in one more accident.\" Sadly, **it's not so easy with logistic regression!**\n", "\n", "## Odds and odds ratios\n", "\n", "LOGistic regression is, yes, all about LOGarithms, and to turn this coefficient into something cool and normal we need to ask Python to undo the logarithm. The opposite of a logarithm is an exponent, so we're going to ask [np.exp](https://docs.scipy.org/doc/numpy/reference/generated/numpy.exp.html) to perform a little magic for us." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>coef</th>\n", " <th>odds ratio</th>\n", " <th>name</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>0</th>\n", " <td>7.853131</td>\n", " <td>2573.780516</td>\n", " <td>Intercept</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", " <td>-0.111171</td>\n", " <td>0.894786</td>\n", " <td>length_in</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " coef odds ratio name\n", "0 7.853131 2573.780516 Intercept\n", "1 -0.111171 0.894786 length_in" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coefs = pd.DataFrame({\n", " 'coef': results.params.values,\n", " 'odds ratio': np.exp(results.params.values),\n", " 'name': results.params.index\n", "})\n", "coefs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **odds ratio** column tells you what happens to the **odds** of completing a scarf each time you add one more of each variable.\n", "\n", "* If the odds ratio is 1, it means nothing changes as we add or remove inches.\n", "* If it's above 1, the odds increase. For example, odds ratio of 2.0 would mean our odds to finish are twice as good for each additional inch.\n", "* If it's below 1, the odds decrease. For example, an odds ratio of 0.5 would mean our odds drop by 0.5x (half) for each additional inch.\n", "\n", "In this case, `length_in` has an odds ratio of about 0.89. This means for each extra inch we add, we have 0.89x the odds we did before - a drop of 11%.\n", "\n", "## What are odds?\n", "\n", "Even though we keep talking about odds, let's be honest: the worst thing in history at the moment is that **you probably have no idea what odds are**. You're _definitely_ (probably maybe) thinking about **probability instead**, because in real life we never ever ever talk about odds.\n", "\n", "Odds and probability (or % chance) are closely related to each other, though. Let's compare the two:\n", "\n", "|words|probability|odds|\n", "|---|---|---|\n", "|---|`finished / attempted`|`finished / unfinished`|\n", "|I start on 100 scarves, I finish half|50/100 = **0.5 or 50% chance**|50:50 odds, 50/50 = **1.0 or 1:1 odds**|\n", "|I finish 20, ignore 80|20/100 = **0.2 or 20% chance**|20:80 odds, 20/80 = **0.25 or 1:4 odds**|\n", "|I finish 75, ignore 25|75/100 = **0.75 or 75% chance**|75:25 odds, 75/25 = **3.0 or 3:1 odds**|\n", "|I finish 996, ignore 4|996/1000 = **0.996 or 99.6% chance**|99.6:0.4 odds, 99.6/0.4 = **249.0 or 249:1 odds**|\n", "\n", "Yeah, it gets pretty weird. If you're at a 100% chance, the odds are 100:0, or 100/0, which is _infinity_.\n", "\n", "### Odds and probability conversion\n", "\n", "To try and make this up to you, I'll show you how to escape from odds and back to probability (...is that even a reward?). Here are a couple formulas to convert between the two measurements:\n", "\n", "```python\n", "probability = odds / (odds + 1)\n", "odds = probability / (1 - probability)\n", "```\n", "\n", "For example, up above our **probability** for finishing a 55-inch scarf was about **0.85**, or 85%. What is that probability as odds?" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.666666666666666" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "0.85 / (1 - 0.85)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "According to The Power Of Mathematics, an 85% probability is an odds of around 5.67." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using the odds ratio\n", "\n", "There's a reason we're spending so much time discussing this, and it's so we can talk about the **odds ratio**. Remember when we looked at the coefficient/odds ratio table before?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>coef</th>\n", " <th>odds ratio</th>\n", " <th>name</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>0</th>\n", " <td>7.853131</td>\n", " <td>2573.780516</td>\n", " <td>Intercept</td>\n", " </tr>\n", " <tr>\n", " <th>1</th>\n", " <td>-0.111171</td>\n", " <td>0.894786</td>\n", " <td>length_in</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " coef odds ratio name\n", "0 7.853131 2573.780516 Intercept\n", "1 -0.111171 0.894786 length_in" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coefs = pd.DataFrame({\n", " 'coef': results.params.values,\n", " 'odds ratio': np.exp(results.params.values),\n", " 'name': results.params.index\n", "})\n", "coefs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we saw above:\n", "\n", "> In this case, `length_in` has an **odds ratio** of about 0.89. This means for each extra inch we add, we have 0.89x the odds we did before (aka 11% less).\n", "\n", "The odds ratio is important in data-driven reporting because it's often used to **explore charges of bias**. For example, you might find that being African American increases the odds of being searched at a traffic stop, or living in a low-income neighborhood decreases your odds of a mortgage being approved. **Learning to use the odds ratio helps you successfully investigate and report on these topics.**\n", "\n", "Time to put our odds ratio to work!\n", "\n", "Let's say we're aiming for a 55-inch scarf, but we'd maybe like it a little longer. But if we're less likely to finish it if it's longer, maybe we should keep it short? Let's use our odds ratio to see how our chance of finishing a 55-inch scarf changes as we **add inches**.\n", "\n", "To start, we need the **odds** for completing a 55-inch scarf. When we did the prediction up above, it gave us a **probability of 0.85** for a 55-inch scarf. Let's convert that probability to odds using our `probability / (1 - probability)` formula." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.666666666666666" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Odds for finishing a 55-inch scarf\n", "odds_55 = 0.85 / (1 - 0.85)\n", "odds_55" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want to add an inch, we use the coefficient/odds ratio table above to find the **odds ratio** of 0.895. If we go one inch longer, we have to multiply the odds by the odds ratio, 0.895. Note that since the odds ratio is **less than 1**, that will make the odds go **down**." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.071666666666666" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Odds for a 56 inch scarf\n", "# Multiply the 55-inch scarf odds by the odds ratio\n", "odds_55 * 0.895" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll multiple it twice if we want to add two inches, moving up to a 57 inch scarf." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.539141666666667" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Odds for a 57 inch scarf\n", "# Multiply the 55-inch scarf odds by the odds ratio twice\n", "odds_55 * 0.895 * 0.895" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And one more for a 58-inch scarf..." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.062531791666666" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Odds for a 58 inch scarf\n", "odds_55 * 0.895 * 0.895 * 0.895" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we'll jump ahead to a 60-inch scarf. Five extra inches means five times multiplied by the odds ratio." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.2541895284197917" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "odds_55 * 0.895 * 0.895 * 0.895 * 0.895 * 0.895" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even though we did all this math, we're still not really hyped about a 3.25 odds ratio. **No one understands odds, what is this number in probability?** We can use the `probability = odds / (odds + 1)` formula." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7649376001422655" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "odds_60 = odds_55 * 0.895 * 0.895 * 0.895 * 0.895 * 0.895\n", "odds_60 / (odds_60 + 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There we go! Adding five extra inches dropped our probability of completion from 85% to 76%.\n", "\n", "We calculated this the long long confusing way, though. Remember our friend `results.predict`? It allows us to **make predictions** if we give it a dataframe. Let's feed it the details for a sixty-inch scarf and see what it thinks the probability should be." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 0.765466\n", "dtype: float64" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sample_df = pd.DataFrame([\n", " {'length_in': 60}\n", "])\n", "\n", "results.predict(sample_df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**It matches our calculation perfectly!** ...or close enough, since we rounded 0.894786 to 0.895 so we could save a little bit of typing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why not percent change?\n", "\n", "**Why don't we just talk about percent increase instead of odds ratio?** Even though that would make way way way more sense to us, we can't! Let's look at why.\n", "\n", "Say we have a really skilled friend who helps us out with knitting, and her help has an **odds ratio of 2.0**. She doubles the odds of completion!!! Incredible!!! But because **odds ratios are about odds and _not_ probability**, her help does **not** take us from 25% to 50%. Here's what a 2x increase in odds does:\n", "\n", "|probability|odds|odds are doubled|convert to %|new percentage|increase|\n", "|---|---|---|---|---|---|\n", "|50/100 or 50%|1:1 or **1.0**|2|`2 / (2 + 1)`|0.67 or **67%**|+17 percentage points|\n", "|20/100 or 20%|1:4 or **0.25**|0.5| `0.5 / (0.5 + 1)`|0.33 or **33%**|+13 percentage points|\n", "|75/100 or 75%|3:1 or **3.0**|6| `6 / (6 + 1)`|0.86 or **86%**|+11 percentage points|\n", "\n", "Even though it was all a 2x odds ratio, each one had a different percent increase! **An odds ratio above 1.0 will always increase your probability, but the change depends on the probability you started at.** For example, the probability change from 55 inches to 57 inches will be different than than 65 inches to 67 inches, even though they're both a 2-inch increase." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Review\n", "\n", "TODO" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.8" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }