011_MML_Logistic_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Logistic Regression\n",
    "***\n",
    "\n",
    "- Logistic regression is the most fundamental algorithm for handling binary classification problems\n",
    "- It's a commonly used algorithm for constructing a baseline model, and is generally a good first model to build, because it's highly interpretable\n",
    "- Logistic regression is similar to linear regression\n",
    "    - We're still dealing with the line equation for making predictions, but this time the result is passed through the sigmoid function to turn ouptut into the probability\n",
    "    - The probability tells you the chance of this instance belonging to the positive classs\n",
    "        - For example, if you have two classes (bad, good), \"good\" is the positive class, and output of 0.8 means the model is 80% condifdent the instance belongs to the positive class\n",
    "- In a nutshell, logistic regression model outputs probabilities instead of raw classes\n",
    "    - We then get the classes by applying a certain threshold\n",
    "    - For example, if the probability is greater than 0.5, we assign the positive class, and negtive class in the other case\n",
    "    - This threshold can (and should) be altered depending on the problem and the type of metric you're optimizing for\n",
    "    - More on that later\n",
    "    \n",
    "**Assumptions** (source: https://www.statisticssolutions.com/assumptions-of-logistic-regression/)\n",
    "- Logistic regression **DOES NOT** assume:\n",
    "    - Linear relationship between the dependent and independent variables\n",
    "    - Normal distribution of the residuals\n",
    "    - Homoscedasticity - constant variance in the residuals\n",
    "- Logistic regression **DOES** assume:\n",
    "    - Independent observations\n",
    "    - Little or no multicollinearity among independent variables\n",
    "    - Linearity of indenpendent variables and log odds\n",
    "    - Large sample size - minimum of 10 cases with the least frequent outcome for each independent variable in your model\n",
    "        - E.g., if you have 5 independent variables and expected probability of the least frequenct outcome is 0.1, then you need a minimum sample size of 500 (10 * 5 / 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Math behind\n",
    "\n",
    "- In a way, we're still dealing with a line equation\n",
    "- Formula of a line equation\n",
    "\n",
    "$$ \\large \\hat{y} = wx + b $$\n",
    "\n",
    "- The **Sigmoid/Logistic** function is applied to the result\n",
    "- The formula for this function is as follows:\n",
    "\n",
    "$$ \\large S(x) = \\frac{1}{1 + e^{-x}} $$\n",
    "\n",
    "- Let's visualize the function to see what happens behind the scenes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0QAAAHCCAYAAADPbwnCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABJcklEQVR4nO3deZhcVZ3/8fe3t6Szr5BANnYSlgTIwg6yyOKCyqCog6Io4Mg4LjPiqKPjNm7jMg4qguOgP1R0FARZZN8hkBBCWBMSCNlDQkL2pLfz+6MqRaXpTjpJd93urvfreeqpOufeqvrWzU13ffrce26klJAkSZKkclSRdQGSJEmSlBUDkSRJkqSyZSCSJEmSVLYMRJIkSZLKloFIkiRJUtkyEEmSJEkqWwYiSeqiIuLCiEgRcWHWteyOiDg5/zn+fSee8+/555y8E885ICJuiIhl+ee+vvPVZi8ixuTrvybrWiSpO6jKugBJUk5EVAIfBf4eOAzoC6wGlgGPAzellG7KrsKuK79t/wLsD/w/YBGwOcuaticiEnB/SunkrGuRpO7OQCRJnUD+C/vNwJnA68At5L60DwL2Az4AHAwUB6IbgKnA0lLW2gEeB8YCKzvwPfYBxgFXp5Qu7sD3KYXF5LbXmqwLkaTuwEAkSZ3D+8mFoaeAk1JK23zZjYhewJTivvw6Xf5LcUppI/BCB7/NXvn7JR38Ph0upVRPx28vSSobnkMkSZ3Dsfn7a5qHIciFhpTSvcV92zuHKCLOiIiHI2JDRKyKiL9ExMERcU3+OWOK1i2ckxIR+0XEnyLitYhYFxF3RMSh+fWGRsRVEbE0IjZHxLSIeEtLHyYi+kfEtyNidn7d1RFxe0Sc1sK6rZ5DFBFHRcTf8rWsjYi7IuKY7W/KN71GAu7PN7+af6/C+7W0TXZUW0Tcl++viogvRsSLEbElIhZGxHcjoqaVWg6OiF9FxPz8+q9GxIMR8Yn88gvz9QKcVFRrcb2tnkMUEcMj4qf516+LiBURcX1EHNXCuoX9JyLekv9MW7fzLRExtm1bWJK6NkeIJKlzeC1/f+DuvlBEvA/4HbAF+CO5Q+qOBR4lNwLVmjHAY8DzwDX59ruB+/Ih5G/AWuAP5A7lOx+4LSIOTCktKHr/AcDD5A5Rmwb8GBgCvBe4IyI+kVL6RRs+x7HAXUANcD0wF5gA3Afcs6PnF/la/rN8mFwwui/ff1/Lq++U3wEnALeR2zZnA58H9gA+UrxiRLwN+D+gB7lt+XtgADA+/5yfAzPz9X4VeIXcv8NW2603IvYBHiI3GnZP/vVHAucBb4uIc1NKN7fw1LcD5+Q/w5Xk/t3OBiZFxLiUUkceyihJ2UspefPmzZu3jG/AEUAd0ETupP/3AKN38JwLgQRcWNS3dSKGLcD4Zut/J79+AsYU9Y8p6v9Ss+f8W75/FbkvyxVFyy7IL/tRs+f8It//CyCK+g8gd4jflmbvf3J+/X8v6gtyh4Ul4Jxmr/9PRfWe3Mbt+6b3KFp2TfNtsqPnkQsnCXgCGFTU35tccGsEhhX1D8l/9jpyh0Q2f58RzdoJuK+Vz7L13+uaZv23t/JveCzQQC5092lh/2kATm32nG/nl30+6/8b3rx589bRNw+Zk6ROIKX0JLnZ5Zbn7/8MzM8funZDRLyjjS91DrlRh9+mlJqPBn2T3IQNrZlPLjQV+3X+vgfwLymlpqJlvyP3ZXrC1o6IqM7Xvx7415TS1sO/SCm9CPyE3IjPh3bwOY4FDgIeSCnd2GzZFcC8HTy/VC5PKa3a2kgpbQB+S+6Q9IlF630Y6Af8PKV0P82klBbtThERMQJ4K7AA+F6z136E3GjRIHJBu7nrUkp3N+u7Kn8/eXfqkqSuwEAkSZ1ESumPwCjgDOAb5GadqwDeBdwUEb+OiNjByxyRv3+ohddfT+6QrNbMTCk1NuvbOgnBnJTSumav10guwI0o6j4Y6AU8VRwUimw91O2IFpYVOzJ/31J4aKSFz5eR6S30LczfDyzqOzp/f1sH1bF1ez6YcpMuNLe97d7WzyBJ3ZKBSJI6kZRSfUrpjpTSV1JK7yB3qNX7gA3kRlXO2cFL9M/fL29leWv90MKMdSmlhtaW5TUA1S28f2tTgW/tH7CdOopfp7V6l+3g+SWRUnq9he6t26yyqG9A/n5xB5WyO9v99eYdRf/ulc2XSVJ3YyCSpE4spdSYHzn6Ub7rlB08ZW3+fs9WlrfW3162BqdhrSwf3my9Hb1Oa/W29vq7YuthgC1NNDSgnd7j9fz93u30es2113aXpLJjIJKkrmHr4Wo7OmTuyfz98c0XREQfis736SCzgY3AhIho6XCrrdN0z9jB62xdflLzBfmL2L7p8+2G1fn7kS0sm9hC366Ymr8/q43rN7FzozOFf/eIaCnYtXW7S1LZMRBJUicQEe+PiNMj4k0/lyNiGPDxfPOBHbzUjeRGAT4YEeObLfsy7Tfi0aKUUh25SQX6AF8vXhYR+wGfAurJzaS3PY+QC1cnRkTzwwQvA/Zrl4JzHs/ff7y4MyIOIzejXXv4NbnRu09ExInNF+YnRSj2Gi0HtBblJ2W4k9wMdJ9u9tpTgA+QC3437EzRklQOvA6RJHUOU8h9+V4WEQ8BL+f79wHeBtSSCzt/2t6LpJTWRsQ/ANcCj0RE8XWIxpObpOAk3jhMrCN8gdy1eS6LiEnAvbxxHaK+wGUppZe383xSSikiLiL3Jf/PEbH1OkTjgdPIXcfnzHaq90bgReD9+WDyGLnJLc7JL3vv7r5BSmllRHyA3L/fvRFxGzCL3Mxzh5MLP/sUPeVu4PyI+Cu5qb0byM24t71AfCm56z99PyLeSm6yhK3XIWoCPtJ8YgxJkoFIkjqLH5D7Un4auS/IZwA9yY0U3EduiuvfFU9j3ZqU0u8iYjW5awi9j9x1fx4AjgH+M7/a2laevttSSqvyF3L9V3LTPH8W2ERuJOb7KaU72vg6D0fECcC3eONQs8fIXRvoDNopEKWUNkfEqeS2zenAJOAZcqMqq2iHQJR/n1siYiJwOXAquWmyV5O73tK3m62+9VpLp5K7SGoFuQu2thqIUkov5V//y/nnnEzu3/lvwLdSStPa43NIUncTbfjdKknqBvLn3rwE9EgpteekBJIkdVmeQyRJ3UxEDIiIXs36gtzIwSjg+kwKkySpE3KESJK6mYg4E/gDcAcwn9wEB0eTm2FuITAxpfRqVvVJktSZGIgkqZuJiH2AbwLHAUPJnS+6CLgZ+I+U0vYuzipJUlkxEEmSJEkqW11+lrkzzzwz/e1vf8u6DEmSJEmdV6sXNu/ykyqsXLky6xIkSZIkdVFdPhBJkiRJ0q4yEEmSJEkqWwYiSZIkSWXLQCRJkiSpbBmIJEmSJJUtA5EkSZKksmUgkiRJklS2DESSJEmSypaBSJIkSVLZMhBJkiRJKlsGIkmSJElly0AkSZIkqWyVLBBFxK8i4tWIeKaV5RERP4mIuRExKyKOLFVtkiRJkspTKUeIrgHO3M7ys4AD8reLgZ+XoCZJkiRJZayqVG+UUnogIsZsZ5VzgN+klBIwNSIGRMTwlNLS0lQoSZKkljQ1JeoaG6lrbGRL/r6uYdv2loZGGpuaaGhKNDQ10dDURGN643FDUxNNTdu2G5pS/jlNNOaXNaZmr9GUaEq5W0qJRKIpUegrXtZif9H6qU3rv7kPIOXvc4/z96RCu/jx1mUtrffmZVtfMxXaxY+3LmtpvTcva/29aWFZc2+8UlFfC+u28vRtttFW3z/rLXxg/CGtPKNzKFkgaoO9gYVF7UX5PgORJElSMyklNtTVs2bLFtZtqWPtljo21tWzsT5321TfkH/cwKb6ejbUvfG4+H5jfT1bGhrZ0thAXWMTWxoaqG9s2iboNDQ1Zf1x1UVtaWjMuoQd6kyBKFroazGARsTF5A6rY9SoUR1ZkyRJUodpakq8vnkzqzZt5rWNm1i1aROrNuYe59qbeX3zZtZtqWPN5q3BZwtrt9SxbktdYfRC0q7rTIFoETCyqD0CWNLSiimlq4CrACZOnOhPAkmS1Kms31LHsvUbWLpuPUvWrWfp2vWF9tJ161m+fiOr8oGnq4SaHlWV1FRW0qOykurKymbtCmoqc+3KiqCqooKqigoqKyqojDfaVRVBZbPHlZF7/MZzYpvnbr2viCACKvKPi2/F/VHcz47Xr2hh/W1eI98PEAFbW5HvLG63tOyNx1tfIwrt4sdtWVZcR/F7sxt1NddSd7SwcitPf9O6vaurW1mz8+hMgegm4LKIuA6YAqzx/CFJktTZNDUllq5fz/zVa3jl9TW8snot81/PPV66Lhd61m2pK0kttdVV9O/Rg749aujbo4Y+NTX0qq6itrqaXtVV9Kqupra6it411fm+6m2Wb+3vWZULNjWVldQUBZ2afPCpqqho8Uux1B2ULBBFxO+Bk4EhEbEI+CpQDZBSuhK4FTgbmAtsBD5SqtokSZKae23jJmavXMWcFauYvfI15qxcxYuvreaV19dS19h+50X061HDoF61DKrtyeBetQzuVcug2loG9+rJoNpaBtT2pH/PGvr26EG/HjXbBKDqysp2q0MqV6WcZe79O1iegE+WqBxJkiQANtbV8+yrK5m17FWeWvoqz766kjkrV7Fy46Zdfs0eVZUM69Ob4X37FG7D+vZmr/zjPfv0YnCvXgyq7UlNlaFGylJnOmROkiSpQ22oq2P64mU8sXgZs5at4Klly5mzcvVOn8czuFctowf0Y/SA/owZ2J8xA/ozemA/RvTry/C+fRhY29NDzKQuwkAkSZK6rYVr1jJ14RIeXbCYxxYu4allr9LY1Lbw06u6igOHDMrdBg/iwKG5+30G9qdfzx4dXLmkUjEQSZKkbmPpuvXcPW8+d819hYcXLGLRmnU7fE4EHDB4EOOH7cHhw4Zy+LA9OHjoYEb060tFhaM8UndnIJIkSV3W5voGHlmwmLvmzeeuufN5evmK7a4fAeOGDmHyiOFMGL4n44fvwaF7DqF3TU2JKpbU2RiIJElSl7J602b++sKL/OW5F7nv5QVsqm9odd3eNdVM2ns4R4/ci2NG7c3kEcMZUNuzhNVK6uwMRJIkqdNbtXETN8+ey5+fmc09Ly2goampxfWqKys4ZuTenLb/GE7ZdzTjh+1BVWVFiauV1JUYiCRJUqe0dvMW/vL8nB2GoP0HD+S0/cZw+v5jOHHMSPr08PA3SW1nIJIkSZ1GSomHXlnEr2c8zQ3PzWFjK4fDTdx7GO855CDeNfYA9hk0oLRFSupWDESSJClzqzdt5v89+QxXT3+Kua+tbnGdrSHoPeMOZPTA/iWuUFJ3ZSCSJEmZef7VlVwxdQbXzXquxdGgQ/YYwgfGj+PcQw4yBEnqEAYiSZJUUlsPi/vRw9O4bc5Lb1rev2cP3nfYwXzoiMM4cq89ifBaQJI6joFIkiSVREqJO+a+zLfue5Rpi5a+aflhew7l0slH8L7DD/a6QJJKxkAkSZI61NYg9M17H2H64mXbLIuAtx+0P586diLHjdrb0SBJJWcgkiRJHebRBYv54h33M3Xhkm36e1RV8vfjD+Gfjp3IAUMGZVSdJBmIJElSB3hx5Sr+7a4HufH5F7fp71FVyUVHjedzx09mr359MqpOkt5gIJIkSe1m3ZY6vnXfI/x06oxtLqRaXVnBRUeN519OmGIQktSpGIgkSdJuSynxx2de4Au338eydRu2WXbeoQfztVOP9wKqkjolA5EkSdotL616nU/edAf3vbxgm/5jR+3Nd888mYl7D8+oMknaMQORJEnaJY1NTVwxdQZfu+chNhVdVHVYn978xxkncf5hY501TlKnZyCSJEk7bfaK1/j4X/62zfWEKiuCf5hyJF8++Vj69eyRYXWS1HYGIkmS1GYpJX45/Skuv/2+bUaFDttzKFe+6wyO3GtYhtVJ0s4zEEmSpDZZsWEjn7jxdm6ZPa/QV1NZyb+edDSfPW4yNVWVGVYnSbvGQCRJknboofkLueD/bmbZ+jdmkDtkjyH877lv47BhQzOsTJJ2j4FIkiS1KqXEjx6exlfufpDGplTo/+TRR/LN006kZ7VfJSR1bf4UkyRJLVqzeQsfv+E2/vrC3ELf0N61/PLdZ/PWA/bJsDJJaj8GIkmS9CbzXlvNub+7gdkrVxX6pozci2vPewcj+vfNsDJJal8GIkmStI0HXl7A+/9wE6s2bS70ffLoI/mP009y4gRJ3Y6BSJIkFfzvE7P41M130dDUBECPqkquPOdMzj98bMaVSVLHMBBJkiRSSnz9nof5zgNTC33D+vTmD+9/F5NHDM+wMknqWAYiSZLKXENjE5+6+U7+d8bThb7xw/bg/z7wLkb275dhZZLU8QxEkiSVsU319XzoT7dwc9FMcmcesA/Xvvcd9K6pybAySSoNA5EkSWVq3ZY6zv3d9Tw4f1Gh74Pjx/Hzc86gutLJEySVBwORJEllaO3mLZxz7Z+ZunBJoe8zx03iW6efSERkWJkklZaBSJKkMvP6ps2889o/M23R0kLfN08/kc8dPznDqiQpGwYiSZLKyOpNm3n7b/6PGUuWF/q+f+ZbuOyYozKsSpKyYyCSJKlMrNtSxznX/nmbMPTjt53KJZOPyLAqScqWgUiSpDKwub6B835/wzaHyV3xjtO5aOL4DKuSpOxVZF2AJEnqWPWNjXzwjzdx/8sLC30/OvtUw5AkYSCSJKlba2pKfOyG27h1zkuFvq+fdgKXTvEwOUkCA5EkSd3aF++8nz8+/UKh/S8nTOFfTpiSYUWS1LkYiCRJ6qZ+NnUG//XI9EL7kkkT+Nqpx2dYkSR1PgYiSZK6oRuff5F//ts9hfY7Dt6fH5x9ihddlaRmDESSJHUzjy1cwoV/uoWUcu3JI4Zzzblvo7LCX/uS1Jw/GSVJ6kYWrVnH+677C5sbGgDYb9AA/vSBd9OrpjrjyiSpczIQSZLUTWyqr+d91/2F5es3AjC4Vy03/v25DO3dK+PKJKnzMhBJktQNpJS49MbbmbFkOQCVFcFv3/sO9hs8MOPKJKlzMxBJktQN/OChx7eZXvsHZ53CSfuMyrAiSeoaDESSJHVxd8+bz1fufrDQvuiow7l40oTsCpKkLsRAJElSF7ZozbptZpQ7bvQIfnj2qU6vLUltZCCSJKmLqm9s5IL/+ysrN24CYFif3lx73juoqarMuDJJ6joMRJIkdVFfvvMBpi5cAkBFBL857+0M69s746okqWsxEEmS1AXd+PyL/OTRJwrtfz/1eE4YMzLDiiSpazIQSZLUxSxcs5ZP3Hh7oX32gfvyueMmZ1iRJHVdBiJJkrqQxqYmPnb9bazetBmAkf37cvW7z6KiwkkUJGlXGIgkSepCfvjwNB6YvxDInTf0v+e+jUG9ajOuSpK6LgORJEldxBOLl/H1ex4utC8/cQrHjR6RYUWS1PUZiCRJ6gI21NVx4Z9voaGpCYDJI4bzxZOOzbgqSer6DESSJHUBX7rzAea+thqAvj1quObct1FV6a9xSdpd/iSVJKmTu++lBfzi8ZmF9g/PPoV9Bg3IrB5J6k4MRJIkdWLrttRx6Y1/K7TfftB+fHD8IRlWJEndi4FIkqRO7Et33M8rr68FYGBtT37yjtOJcIptSWovJQ1EEXFmRMyOiLkR8YUWlvePiL9GxFMR8WxEfKSU9UmS1JncM+8Vrp7+VKH9w7NPZXjfPhlWJEndT8kCUURUAj8FzgLGAe+PiHHNVvsk8FxKaTxwMvCDiKgpVY2SJHUW67fU8Ymbbi+033Hw/rzvsIMzrEiSuqdSjhBNBuamlF5KKdUB1wHnNFsnAX0jdyxAH2AV0FDCGiVJ6hS+ce/DLMgfKjeotic/ebuHyklSRyhlINobWFjUXpTvK3YFMBZYAjwN/FNKqan5C0XExRExPSKmr1ixoqPqlSQpE08uWc4VU2cU2t8982SG9e2dYUWS1H2VMhC19Get1Kx9BjAT2AuYAFwREf3e9KSUrkopTUwpTRw6dGh71ylJUmYaGpv45E130JRyvyJP3meUs8pJUgcqZSBaBIwsao8gNxJU7CPA9SlnLvAy4AHTkqSy8bPHZvDk0uUA9Kiq5L+dVU6SOlQpA9E04ICI2Cc/UcL5wE3N1lkAnAoQEXsCBwEvlbBGSZIy88rra/j6vQ8X2v964jHsP3hghhVJUvdXVao3Sik1RMRlwO1AJfCrlNKzEXFpfvmVwDeAayLiaXKH2F2eUlpZqholScrS5/92Hxvq6gEYt8dgPnPcpIwrkqTur2SBCCCldCtwa7O+K4seLwHeWsqaJEnqDO6c+zI3Pf9iof3f73grNVWVGVYkSeWhpBdmlSRJb7aloYHP3npPof33Ew7h2FHNJ2KVJHUEA5EkSRn7yaNPMPe11QD079mDb55+YsYVSVL5MBBJkpShRWvW8Z37Hy20/+0tx7FnH685JEmlYiCSJClDX7j9PjbWNwBw6J5DuGTShGwLkqQyYyCSJCkjD81fyJ+fnV1o/+jsU6mq9FezJJWSP3UlScpAU1Pi8tvvK7TPO/Rgjh8zsvUnSJI6hIFIkqQM/H7Wc8xYshyAnlVVTqQgSRkxEEmSVGIb6ur4yt0PFtr/dOxERg3ol2FFklS+DESSJJXYjx+ezpK16wEY1qc3nzt+csYVSVL5MhBJklRCi9eu44cPP15of/WU4+nboybDiiSpvBmIJEkqoa/f83Bhmu3Dhw3lgiMOybgiSSpvBiJJkkrkuVdXcu3MZwvt75xxMpUV/iqWpCz5U1iSpBL56l0P0pQSAG/dfwxv2Xd0xhVJkgxEkiSVwMOvLOLm2fMK7a+f5jTbktQZGIgkSepgKSX+7a43ptk+//CxjB++R4YVSZK2MhBJktTBbp3zEo8uWAxAdWUFXz3luIwrkiRtZSCSJKkDNTY18ZW7Hii0Pz5xAmMGDsiuIEnSNgxEkiR1oD88/TzPvfoaAH1qqrn8xCkZVyRJKmYgkiSpg9Q3NvKt+x4ttD91zET26NM7w4okSc0ZiCRJ6iDXznyWl1a9DsDA2p586tiJ2RYkSXoTA5EkSR1gS0MD377/jdGhzxw7if49e2RYkSSpJQYiSZI6wDUznmbhmnUADO1dyyemHJFxRZKklhiIJElqZ5vq6/nOA1ML7c8dP4U+PWoyrEiS1BoDkSRJ7ezqaU+xbN0GAIb17c3Fk8ZnXJEkqTUGIkmS2tHGunp++PDjhfblJxxNbXV1hhVJkrbHQCRJUjv61ROzWL5+IwB79+vLR446LOOKJEnbYyCSJKmdbK5v4AdFo0P/fPxkelRVZViRJGlHDESSJLWTa2Y8vc25Qxce6eiQJHV2BiJJktrBloYG/vOhN0aHPnfcZHpWOzokSZ2dgUiSpHZw7cxnWbw2d92hPXr34qNHHZ5xRZKktjAQSZK0m+obG/neg48V2p8+bhK9apxZTpK6AgORJEm76bpZz7Pg9bUADOlVy8cnet0hSeoqDESSJO2GpqbED4rOHfrHY46iT4+aDCuSJO0MA5EkSbvhljnzmL1yFQB9e9Rw8aQJ2RYkSdopBiJJknZRStuODn1s4ngG1PbMsCJJ0s4yEEmStIseXrCYxxYuAaCmspLLjj4q44okSTvLQCRJ0i76YdHo0AfGj2Ovfn0yrEaStCsMRJIk7YJnlq/gtjkvARABnzluUsYVSZJ2hYFIkqRd8MOHphUev/PgAzhwyKAMq5Ek7SoDkSRJO+mV19fwx2eeL7Q/d/zkDKuRJO0OA5EkSTvpJ488QWNTAuDEMSOZNGJ4xhVJknaVgUiSpJ3w2sZNXDNjVqHt6JAkdW0GIkmSdsKVjz3JxvoGAA4fNpTT9x+TbUGSpN1iIJIkqY021NXxs8dmFNqfPX4yEZFhRZKk3WUgkiSpja6Z8QyrNm0GYMzA/pw77qCMK5Ik7S4DkSRJbVDf2Mh/PfLGVNv/dMxEqir9NSpJXZ0/ySVJaoM/PTubhWvWATCkVy0fOuLQjCuSJLUHA5EkSTuQUuK/H3mi0P6HKUfSq6Y6w4okSe3FQCRJ0g48vGAxTy5dDkDPqio+Pml8xhVJktqLgUiSpB244tE3Roc+MH4cQ3r3yrAaSVJ7MhBJkrQd81e/zl9fmFtoX3b0kRlWI0lqb1U7WiEiKoCTgZOAMUAtsAKYAdyRUlrYgfVJkpSpnz32JE0pAXDafmMYu8eQjCuSJLWnVkeIIqI2Ir4ELARuAd4K9AHqgH2ArwIvR8StEXF0KYqVJKmU1m7ewjUzni60//GYozKsRpLUEbY3QvQi8ChwMbmRoPrmK0TEaOADwB8i4psppas7pkxJkkrvN08+w7otdQAcNGQQp+03JtuCJEntbnuB6MyU0jPbe3JK6RXg2xHxA2B0u1YmSVKGGpua+OljMwrtTx59JBUVkWFFkqSO0Oohc8VhKCK2e65RSqkupfRiexYmSVKWbpk9j/mr1wAwsLYnHxx/SMYVSZI6QltnmVsaEf8ZEWM7tBpJkjqJK6a+MdX2RUcd7oVYJambamsg+iJwLPBMRDwaERdFRJ8OrEuSpMzMXLqcB+cvAqCqooJLJh+RcUWSpI7SpkCUUro6pXQscCjwEPBNcqNGv4qI4zqyQEmSSq34QqzvOeRARvTvm2E1kqSOtFMXZk0pPZ9S+hdgBLlRow8AD0TECxFxaf6aRZIkdVlL163nj8+8UGhfdrRTbUtSd7ZTASYiaiLifOA24EfAVOBC4H+BfwN+t4PnnxkRsyNibkR8oZV1To6ImRHxbETcvzP1SZK0u66eNpP6xiYAjh65F5NGDM+4IklSR9ru7HFbRcSRwEeB9wP1wG+ATxbPLBcRdwMPbuc1KoGfAqcDi4BpEXFTSum5onUGAD8jN+X3gojYY6c/kSRJu2hTfT1XT3uq0PZCrJLU/bUpEAGPA3eSu0jrjSmlhhbWeR64bjuvMRmYm1J6CSAirgPOAZ4rWucDwPUppQUAKaVX21ifJEm77Q+zXmDlxk0AjOzfl3cefEDGFUmSOlpbA9F++YuwtiqltAH4yHZW2RtYWNReBExpts6BQHVE3Af0Bf4rpfSb5i8UEReTC2eMGjVqh8VLkrQjKSV+/vgbF2L9xJQjqar01FhJ6u7a+pP+3ogY3LwzIgZExEttfI2WLu+dmrWrgKOAtwFnAP8WEQe+6UkpXZVSmphSmjh06NA2vr0kSa2bunAJs5atAKC2uooLjzws44okSaXQ1hGiMUBlC/09yI38tMUiYGRRewSwpIV1VuZHmzZExAPAeGBOG99DkqRd8ovHZxYev++wsQys7ZldMZKkktluIIqI9xQ13xYRa4ralcCpwPw2vtc04ICI2AdYDJxP7pyhYjcCV0REFVBD7pC6H7Xx9SVJ2iXL12/g+udmF9qXTJ6QXTGSpJLa0QjRn/L3CfifZsvqyYWhz7XljVJKDRFxGXA7uTD1q5TSsxFxaX75lSml5yPib8AsoAn4ZUrpmTZ9EkmSdtGvZzxdmGp7ysi9mDB8z4wrkiSVynYDUUqpAiAiXgYmpZRW7s6bpZRuBW5t1ndls/b3ge/vzvtIktRWDY1NXD39jam2L5k0IbtiJEkl16ZziFJK+3R0IZIkZeHWOfNYtGYdAEN71/KeQ940l48kqRtrNRBFxGeBn6WUNucftyql9MN2r0ySpBK4atrMwuMLjzycHlVtnW9IktQdbO+n/j8CvwY25x+3JgEGIklSl/PiylXcPS93mb2KCC6aeHjGFUmSSq3VQFR8mJyHzEmSuqNfFI0OnX3gvowe0D+7YiRJmfAS3JKksrShro5rZz5baF8y+YgMq5EkZWVH5xC1iecQSZK6mj/MeoE1m7cAsN+gAZyy7+iMK5IkZWFH5xC1hecQSZK6lJQSv5j2ZKF98aQJVFREhhVJkrLSpnOIJEnqTh5duIRZy1YAUFtdxQVHHJpxRZKkrHgOkSSp7Fz1+MzC4/MPG8vA2p7ZFSNJypTXIZIklZXl6zdw/XOzC+2LJ0/IrhhJUua8DpEkqaxc88TT1Dc2ATBl5F5MGL5nxhVJkrLkdYgkSWWjobGJXz7xVKF9iaNDklT2PIdIklQ2bp0zj0Vr1gEwtHct7xl3YMYVSZKy1uZAFBHviogHImJl/vZgRLy7I4uTJKk9/aJoMoULjzycHlXbO3JcklQO2hSIIuJzwB+A2cDn87cXgN9FxD93XHmSJLWPOStXcc9LrwBQEcFFEw/PuCJJUmfQ1j+N/TNwWUrp6qK+X0XE48DXgf9s98okSWpHV02bWXh89oH7MnpA/+yKkSR1Gm09ZK4PcG8L/ffml0mS1GltqKvj2pnPFtqXTD4iw2okSZ1JWwPRX4C/a6H/XOCmdqtGkqQOcN2s51mzeQsA+w8eyCn7js64IklSZ7GjC7NuNRf4QkS8BXg033d0/uY1iCRJnVZKaZvJFD4+cTwVFZFdQZKkTmVHF2Yttho4MH8r7ruQ3HlEkiR1Oo8uXMLTy1cAUFtdxQVHHJpxRZKkzqRNF2aVJKmr+sXjTxYen3/YWAbW9sywGklSZ+OFWSVJ3dby9Ru44bk5hfbFkydkV4wkqVNq8xXpIuJAchMrjAJqipellD7aznVJkrTbrnniaeobmwCYMnIvJgzfM+OKJEmdTZsCUUS8Dfgz8CRwFDAN2A/oATzYYdVJkrSLGhqbuHr6U4X2JY4OSZJa0NZD5r4OfC2ldAywBbgAGAPcBdzXIZVJkrQbbp0zj8Vr1wEwtHct7xl34A6eIUkqR20NRAcBf8g/rgd6pZQ2kwtKn+6AuiRJ2i3FU21feOTh9Khq81HikqQy0tZAtA7YOi3PUmD//OMqYGB7FyVJ0u6Ys3IV97z0CgAVEVw08fCMK5IkdVZt/XPZY8DxwHPALcAPImI88G7euFCrJEmdwlXTZhYen33gvowe0D+7YiRJnVpbA9FngT75x/8O9AXOBebkl0mS1ClsqKvj2pnPFtqXTD4iw2okSZ1dmwJRSumloscbgU90WEWSJO2G62Y9z5rNWwDYf/BATtl3dMYVSZI6s506wzQiTgHG5ZvPpZTuaf+SJEnaNSmlbSZT+PjE8VRURHYFSZI6vbZeh2gf4HrgMGBJvnuviHgaOLd4BEmSpKw8unAJTy9fAUBtdRUXHHFoxhVJkjq7ts4y9z/AWmDflNKolNIoYF/gdeCXHVSbJEk75RePP1l4fP5hYxlY23M7a0uS1PZD5o4Bjk4pLdjakVJaEBGfwVnmJEmdwPL1G7jhuTmF9sWTJ2RXjCSpy2jrCNECoLaF/p7AwvYrR5KkXXPNE09T39gEwJSRezFh+J4ZVyRJ6graGog+B/wkIo6OiMr87Wjgx/llkiRlpqGxiaunP1VoX+LokCSpjVo9ZC4i1gGpqKsn8DDQlG9XAI3Ab4F+HVWgJEk7csuceSxeuw6Aob1rec+4AzOuSJLUVWzvHKLLSlaFJEm74aqiqbYvPPJwelTt1FUlJEllrNXfGCmlX5eyEEmSdsWclau456VXAKiI4KKJh2dckSSpK2nzn9AiogfwQXIXZk3As8DvU0pbOqg2SZJ2qPhCrGcfuC+jB/TPrhhJUpfTpkkVImIc8CLwQ2AKsHVChTkRMbbDqpMkaTs21NVx7cxnCu1LJh+RYTWSpK6orbPM/RfwJDAqpXRCSukEYBTwFLlgJElSyV0363nWbqkDYP/BAzll39EZVyRJ6mraesjcccCklNLarR0ppbUR8SVgaodUJknSdqSUtjlc7uMTx1NREdkVJEnqkto6QrQZGNBCf//8MkmSSuqRBYt5evkKAGqrq7jgiEMzrkiS1BW1NRD9Fbg6Io4rujDr8cAvgJs6rjxJklp21bSZhcfnHzaWgbU9sytGktRltTUQ/RO5SRUeJDcitBm4H5gDfLpDKpMkqRXL1m3ghufmFNoXT56QXTGSpC5th+cQRUQFsBe5KbeHAWOBAJ5LKc3t2PIkSXqz/31iFvWNTQBMGbkXE4bvmXFFkqSuqi2TKiRgJjAuH4AMQZKkzNQ3NnL19KcK7U9McaptSdKu2+EhcymlBMwGhnZ8OZIkbd9NL8xl6br1AOzZpxfvHntgxhVJkrqytp5D9Hng+xExISKc01SSlJkrH3uy8Piio8ZTU1WZYTWSpK6urdch+iPQE3gCaIiILcULU0r92rswSZKae3rZCh56ZREAVRUVXDRxfMYVSZK6urYGoss6tApJktrgF4+/MTp0ztgD2KtfnwyrkSR1B20KRCmlX3d0IZIkbc/qTZv5/aznCu1LnUxBktQO2jpCRET0BD4AjMt3PQf8PqW0qSMKkySp2G+efIaN9Q0AHD5sKMeN2jvjiiRJ3UGbJlWIiCOBl4EfAJPzt/8EXsovkySpwzQ1pW0Ol7t08hE4x48kqT20dZa5q4AHgREppRNTSicCI4EH8sskSeowd8x9mZdXrwFgQM8evO+wsRlXJEnqLtp6yNwhwIdSShu2dqSUNkTE14HpHVKZJEl5VxaNDn34yMPoVVOdYTWSpO6krSNELwB7tdA/HJjTfuVIkrStea+t5vYXXwYgAi6eNCHbgiRJ3UpbR4i+DPwkPyI0Nd93dL7/CxExaOuKKaVV7VuiJKmc/WLazMLjMw/Yl30HDcisFklS99PWQPTX/P3vgJR/vPVs1huL2gnwkuGSpHaxoa6O3zz5TKF96WSn2pYkta+2BqK3tMebRcSZwH+RC02/TCl9p5X1JpEbiXpfSulP7fHekqSu57pZz7Nm8xYA9h88kNP2G5NtQZKkbqetF2a9f3ffKCIqgZ8CpwOLgGkRcVNK6bkW1vsucPvuvqckqetKKW0zmcIlkyZQUeFU25Kk9tXqpAoRsU9bXyRyRu5gtcnA3JTSSymlOuA64JwW1vtH4M/Aq219f0lS9/Pg/IU8s3wlAL1rqvn7CYdkXJEkqTva3ixzj0bE/0TEMa2tEBEDI+ITwHO0HG6K7Q0sLGovyvcVv97ewLuBK7f3QhFxcURMj4jpK1as2MHbSpK6ov9+9InC4w8cPo4BtT0zrEaS1F1t75C5g4EvAbdERCPwBLAU2AwMBMYBY4HHgU+nlHZ0iFtLxzmkZu0fA5enlBq3dwXylNJV5C8IO3HixOavIUnq4ua9tppb5swrtC875qgMq5EkdWetBqKU0uvAv0TEV4C3AccDo4FaYCXwa+D2lNIzrb1GM4uA4sPqRgBLmq0zEbguH4aGAGdHRENK6S9tfA9JUjfw08dmkPJ/7jrzgH04cMig7T9BkqRdtMNJFVJKm4A/5W+7YxpwQP7cpMXA+cAHmr1X4byliLgGuNkwJEnl5fVNm7eZatvRIUlSR2rrtNu7LaXUEBGXkZs9rhL4VUrp2Yi4NL98u+cNSZLKwzUznmZDXT0A4/YYzCn7js64IklSd1ayQASQUroVuLVZX4tBKKV0YSlqkiR1Hg2NTfy8aKrty44+iu2dUypJ0u7a3ixzkiSV1E0vvMiC19cCMKRXLecfPjbjiiRJ3Z2BSJLUaVwxdUbh8ccmjqe2ujrDaiRJ5WC7gSgiji5VIZKk8jZ98VIeXbAYgOrKCi6ePCHbgiRJZWFHI0QPRMQ3IqKk5xpJksrPFY++MTp03qEHM7xvnwyrkSSVix0ForOAC4DHI2JcCeqRJJWhxWvX8ednZxfalx3tVNuSpNLYbiBKKd0NHAY8CUyPiM+WpCpJUln5xeMzaWhqAuD40SM4Yq89M65IklQudjipQkppXUrpInIjRd+LiPURsbb41vFlSpK6q4119fzP9KcKbS/EKkkqpTadGxQRE4FvAi8C/wk0dGRRkqTy8bunnmPVps0AjBnYn7cftF/GFUmSysl2A1F+MoWvApcDPwO+kFLaXIrCJEndX1NT4oqpTxTa/zDlSCorvCKEJKl0djRCNA0YBJyVP59IkqR2c+e8l5m9chUAfXvU8OEjDs24IklSudnRn+GeAQ43DEmSOsIPHnq88PjDRxxGv549MqxGklSOtjtClFK6oFSFSJLKy2MLl/Dg/EUAVFVU8KljnUxBklR6HqgtScpE8ejQ+w47mJH9+2VYjSSpXBmIJEklN3vFa9w8e26h/dnjJ2dYjSSpnBmIJEkl96OHp5FS7vHZB+7LuD2GZFuQJKlsGYgkSSW1aM06fjfruULb0SFJUpYMRJKkkrpi6hPUNzYBcMyovTlu9IiMK5IklTMDkSSpZFZv2sz/TH+q0P6co0OSpIwZiCRJJXP1tJmsr6sHYOzQwZx1wL4ZVyRJKncGIklSSWyqr+enU2cU2p89fhIVFZFhRZIkGYgkSSXy/558llc3bARg7359ee+hYzOuSJIkA5EkqQQaGpv48SPTCu1PHXMUNVWVGVYkSVKOgUiS1OFueH4OL69eA8DA2p589KjDM65IkqQcA5EkqUOllPjhQ48X2pdOnkCfHjUZViRJ0hsMRJKkDnX3vFeYufRVAGqrq/jElCMzrkiSpDcYiCRJHSalxLfvf7TQ/tARhzK0d68MK5IkaVsGIklSh3lw/kIeWbAYgOrKCj57nBdilSR1LgYiSVKH+Y+i0aELJhzKqAH9MqxGkqQ3MxBJkjrEw68s4v6XFwJQWRH8ywmODkmSOh8DkSSpQ3zn/qmFxx8cfwhjBg7IrhhJklphIJIktbvHFi7hrnnzAaiI4PMnTMm2IEmSWmEgkiS1u+8UnTv0vsPGst/ggRlWI0lS6wxEkqR29djCJfztxZcBiIDLT3R0SJLUeRmIJEnt6hv3Plx4fN6hB3PQ0MEZViNJ0vYZiCRJ7eahVxZx97xXgNy5Q18++diMK5IkafsMRJKkdpFS4mt3P1Rof3D8OA4YMijDiiRJ2jEDkSSpXdz70gIeemURAFUVFXzx5GMyrkiSpB0zEEmSdltKia/d88bo0IePONTrDkmSugQDkSRpt93+4ss8vmgpADWVlVx+0tEZVyRJUtsYiCRJu6WpKfGVux4stC+aeDgj+/fLsCJJktrOQCRJ2i1/ePp5nl6+AoBe1VV8/gRHhyRJXYeBSJK0y7Y0NGxz7tCnjpnIsL69M6xIkqSdYyCSJO2yq6c9xSuvrwVgcK9aPnPcpIwrkiRp5xiIJEm7ZO3mLXz3gamF9uUnHk2/nj0yrEiSpJ1nIJIk7ZIfPzKNlRs3ATBqQD8unjQ+44okSdp5BiJJ0k5bum49P3n0iUL7q6ccR4+qqgwrkiRp1xiIJEk77Wt3P8SGunoADttzKO87bGzGFUmStGsMRJKknfLU0lf5zcxnCu1vn3ESlRX+OpEkdU3+BpMktVlKictvv5eUcu2zDtyXU/cbk2lNkiTtDgORJKnNbp49j/tfXghAZUXwH289KeOKJEnaPQYiSVKb1DU08q+331dof2zieA4eOji7giRJagcGIklSm1z5+JPMW/U6AP179uDLJx+bbUGSJLUDA5EkaYeWr9/At+57pND+wolHM6R3rwwrkiSpfRiIJEk79OU7H2DtljoA9h88kE9MOSLjiiRJah8GIknSdj26YDHXzny20P7h2ad4EVZJUrdhIJIktaqxqYnP3HJ3of3OsQdw+v77ZFiRJEnty0AkSWrVL6c/xVPLXgWgZ1UV3zvj5GwLkiSpnRmIJEktWrFhI1+75+FC+/MnTGH0wP4ZViRJUvszEEmSWnT53+5j9abNAOwzsD+fOW5SxhVJktT+DESSpDe5e958fj/ruUL7R287lZ7VTqQgSep+ShqIIuLMiJgdEXMj4gstLP9gRMzK3x6JiPGlrE+SBBvr6vnHv95ZaP/doQdxxgH7ZliRJEkdp2SBKCIqgZ8CZwHjgPdHxLhmq70MnJRSOhz4BnBVqeqTJOV8+/5HeXn1GgAG9OzB9888JeOKJEnqOKUcIZoMzE0pvZRSqgOuA84pXiGl9EhKaXW+ORUYUcL6JKnsPb1sBT9+ZHqh/a3TT2JY394ZViRJUscqZSDaG1hY1F6U72vNRcBtLS2IiIsjYnpETF+xYkU7lihJ5auhsYl/uOl2GpqaADhu9AguPPKwjKuSJKljlTIQRQt9qcUVI95CLhBd3tLylNJVKaWJKaWJQ4cObccSJal8/fiRaUxfvAyA6soKrnjH6VRUtPSjW5Kk7qOUUwYtAkYWtUcAS5qvFBGHA78EzkopvVai2iSprD336kq+ce8jhfaXTz6Wg4cOzrAiSZJKo5QjRNOAAyJin4ioAc4HbipeISJGAdcDF6SU5pSwNkkqWw2NTXz8htuoa2wE4Mi99uSzx03OuCpJkkqjZCNEKaWGiLgMuB2oBH6VUno2Ii7NL78S+AowGPhZRAA0pJQmlqpGSSpHP3z4cWYsWQ5ATWUlv3z3WVRVepk6SVJ5KOlV9lJKtwK3Nuu7sujxx4CPlbImSSpnzyxfwTfve+NQuX97y7GM3WNIhhVJklRa/glQksrUpvp6Pvynm6lvzM0qN2nEcD597KSMq5IkqbQMRJJUpr54xwM892pu7pra6iqufteZHionSSo7/uaTpDJ025x5XPn4k4X29854Cwc5q5wkqQwZiCSpzCxbt4FL/vK3QvvtB+/PRRMPz7AiSZKyYyCSpDLS1JS4+C+3sWLDJgCG9e3Nz9/5VvIze0qSVHYMRJJURr7/4GPcOXd+of0/7z6bIb17ZVeQJEkZMxBJUpm4/+UFfP3ehwvtzx43iVP2G51hRZIkZc9AJEllYMna9Xzo/26mKSUAjhs9gq+dekLGVUmSlD0DkSR1cw2NTXz4Tzfz6oaNAOzRuxe/+bu3O8W2JEkYiCSp2/vSnffz0CuLAKiI4Nd/93b26tcn46okSeocDESS1I1dO/MZfvLoE4X2v73lWE7ed1SGFUmS1LkYiCSpm3ps4RI+edOdhfbbD96fz59wdIYVSZLU+RiIJKkbWrx2HedfdyN1jY0AjNtjML96z9lUVHi9IUmSihmIJKmb2VBXx3t//xeWrd8AwKDanvzp/e+mb4+ajCuTJKnzMRBJUjfS2NTEhX+6hRlLlgNQWRH89r3vZJ9BA7ItTJKkTspAJEndREqJz916DzfPnlfo+8FZpziJgiRJ22EgkqRu4ocPT+MX02YW2p85bhKXTD4iu4IkSeoCDESS1A1cN+t5vnznA4X2eYcezDdPOzHDiiRJ6hoMRJLUxd38wlw+dsOthfbxo0dw9bvPdEY5SZLawEAkSV3YPfNe4YN//CuNTQmAsUMH88f3v4seVVUZVyZJUtdgIJKkLurRBYs57/c3FK41tO+gAdzy4fMYWNsz48okSeo6DESS1AVNW7SUd//2ejbWNwCwV78+3Pqh8xjet0/GlUmS1LV4TIUkdTFTFy7hnP/3J9ZuqQNgaO9abv3QeYwe2D/jyiRJ6noMRJLUhTz0yiLefe2fWV9XD8DgXrX89YLzOGjo4IwrkySpazIQSVIXcf/LC3hP0WFyQ3vXcuuH38uhew7NuDJJkrouA5EkdQE3PDeHC/90S2EChWF9enPrh89j7B5DMq5MkqSuzUAkSZ3c1dNm8k+33EXKzazN8L59+NuF7+XAIYOyLUySpG7AQCRJnVRKiW/d9wjfuu/RQt8Bgwfy1wv+zgkUJElqJwYiSeqENtc38A833cHvZz1X6Ju49zCu/+B7GNq7V4aVSZLUvRiIJKmTWb5+A++77kYeW7ik0Hf6/mP43XvfSZ8eNRlWJklS92MgkqRO5OllKzj3d9ezcM26Qt9HjjyMH7/tNGqqKjOsTJKk7slAJEmdxO+eepbL/nonm/LTaldE8O23nsQ/HnMUEZFxdZIkdU8GIknK2JaGBv7ltnu5evpThb6+PWr49d+9jbMO3C/DyiRJ6v4MRJKUoZdXvc4Ff7qZJxYvK/QdNGQQv3/fO73GkCRJJWAgkqQMpJT47VPP8plb7mZ9XX2h/9xDDuLn55xBXydPkCSpJAxEklRiqzdt5rK/3sH1z84p9FVVVPDtt57EJ48+0vOFJEkqIQORJJXQLbPn8Y9/vZOl69YX+vYfPJD/PfdsJu49PMPKJEkqTwYiSSqBFRs28rlb7+H/nnlhm/6PHnU43zvzZHrXeIicJElZMBBJUgdqasqdK/TFO+5n5cZNhf49evfiv99xOu8ce0CG1UmSJAORJHWQJ5cs5zO33s1jC5ds0//3Ew7hu2eczKBetRlVJkmStjIQSVI7W75+A9+69xH+54lZNKVU6B/Zvy9XvOOtvPWAfTKsTpIkFTMQSVI72VBXx08eeYIfPvz4NlNp11RW8uljJ/L5E6d4rpAkSZ2MgUiSdtOWhgaumfE033lgKsvWbdhm2en7j+EHZ53CAUMGZVSdJEnaHgORJO2iLQ0N/HrGM3zvwcdYvHbdNsvG7TGYb5x2ImcduK/XFZIkqRMzEEnSTlq7eQu/emIWV0yd8aYgNKxvb776luP5+wmHUFVZkVGFkiSprQxEktRGi9eu4+ePPckvpz/Fms1btlm2Z59efPa4yXxs4nh61VRnVKEkSdpZBiJJ2o6UEg/MX8iVjz/JX1+YS2NT2mb5nn168ZnjJvNxg5AkSV2SgUiSWvDq+g38ftbz/HrG0zy/4rU3LT9g8EA+fewkPjB+HD2r/VEqSVJX5W9xScrbVF/P7S++zO9nPc+ts+fR0NT0pnVOHDOSTx59JG8/aH8qKpwsQZKkrs5AJKmsbQ1Bf352NrfNeYkNRdcP2qp3TTUfHD+OSyYfwbg9hmRQpSRJ6igGIkllZ2NdPbfPfZnrtxOCAI4ZtTcfPuJQ3nPIQfTt4QVVJUnqjgxEkrq9lBLPLF/JnXNf5u55r/DwgkVsaWhscd0DhwziPeMO5P3jx3GgF1OVJKnbMxBJ6pZeXb+Be15awN3z5nPXvPksW7eh1XW3hqBzDz2IQ/YY4oVUJUkqIwYiSV1eU1Pi+RUrmbpwCY8uWMzUhUuYt+r17T5n7NDBvGvcAbznEEOQJEnlzEAkqUtJKbFo7TpmLVvBU0uX89jCJTy2aOmbLpTa3MDanpyy72hO338Mp+w3mpH9+5WoYkmS1JkZiCR1WlsaGnjxtdU8vXwFTy19lVnLXmXWshW8tnHTDp9bU1nJxL2Hcdp+Yzh1/zEctdeeVFZUlKBqSZLUlRiIJGUqpcSKDRuZ89pq5qxYxeyVr/Hia6uZvXIV81evoSmlNr3OHr17cfSovTl65F4cM2pvjhi+Bz2q/BEnSZK2z28LkjpUSonXNm5i8dr1vPL6GuavXsP819fwytb719e2Ou11a/r37MHhw4Zy+LA9OGL4nhwzam/2Gdjf84AkSdJOMxBJ2iVbGhpYsWETKzduZOWGTazYuJFX129kydp1LFm3niVr17Nk3XqWrlvf6hTXOxIBYwb0Z+zQwYwfvgfjh+/J+GFDGT3A8CNJktqHgUgqYyklNjc08PrmLazZvCV/v5k1m+tYs3kzr2/ewuubNvPaxk2s3LiJFRs2Fu7Xbalrtzr61FRzwOBBHDhkEAcNGcSBQ3P3+w8aSM9qf0xJkqSOU9JvGhFxJvBfQCXwy5TSd5otj/zys4GNwIUppRmlrFHqLFJK1DU2sqWhkc0NDWxpbKSuoZEt+b7cfQMb6+vZUFfPxvp61m+pZ0N9PRvrtr3fUJe/FfWt27KFNZvrqGvctdGbndGvRw179evDqP79GDOwP6MH9N/mflBtT0d8JElSJkoWiCKiEvgpcDqwCJgWETellJ4rWu0s4ID8bQrw8/x9l5PyJ4IXnw+e2LZva7ulvp1bp4X3baW96+u0T81sd3u88d4JaEqpxVtKLS9Lhcfbee52X7eo3bTteg1NTYVbY1Oivqkx3952WUNj0ePmy5rdtjQ0FgLP1pBT19jA5obGwrLOrLIiGNKrliG9ejGkdy1De/diSK9e7NW3D3v168Pe/fqyV78+7NW3D3161GRdriRJUotKOUI0GZibUnoJICKuA84BigPROcBvUu7b+NSIGBARw1NKS0tY5y45/qpreWLxsqzLkHZaTWUlA3r2oH/+NqBnz/x9D/r17EH/Hj0Y3KuWoX16MaTX1uBTy4CePamocFRHkiR1baUMRHsDC4vai3jz6E9L6+wNdPpAtHVERWpP1ZUV9KyqokdlJTVVlfSorKRHVSU9q6qoyT/uVV1Nr+pqetdU07u6ml411fSpKeorerz1vk9NDb1rqhnYs6fn6EiSpLJWym9CLf0puXmKaMs6RMTFwMUAo0aN2v3KOljxqRGR/4hb+6LoIzfv2/46W9uxTbulvl1fZ8f1sBs1b+/9KyLedIugWfvN61REK88tes2Wn9f68yOCqoqK/K34cQXVlRVUFrWLl795WW55dUUllRVBdWUlPQsh542AUwhAlZWOwEiSJHWwUgaiRcDIovYIYMkurENK6SrgKoCJEyd2iqGZBz/+92988ffkcEmSJKlLqCjhe00DDoiIfSKiBjgfuKnZOjcBH4qco4E1XeH8IYCKitxIgmFIkiRJ6jpKNkKUUmqIiMuA28lNu/2rlNKzEXFpfvmVwK3kptyeS27a7Y+Uqj5JkiRJ5Se6+mQAEydOTNOnT8+6DEmSJEmdV6uHcZXykDlJkiRJ6lQMRJIkSZLKloFIkiRJUtkyEEmSJEkqWwYiSZIkSWXLQCRJkiSpbBmIJEmSJJUtA5EkSZKksmUgkiRJklS2DESSJEmSypaBSJIkSVLZMhBJkiRJKluRUsq6ht0SESuAV7Kuo8gQYGXWRXRzbuOO5zbueG7j0nA7dzy3ccdzG3c8t3HHy3obr0wpndnSgi4fiDqbiJieUpqYdR3dmdu447mNO57buDTczh3Pbdzx3MYdz23c8TrzNvaQOUmSJElly0AkSZIkqWwZiNrfVVkXUAbcxh3Pbdzx3Mal4XbueG7jjuc27nhu447Xabex5xBJkiRJKluOEEmSJEkqWwYiSZIkSWXLQLSTIuK8iHg2IpoiYmKzZf8aEXMjYnZEnNHK8wdFxJ0R8WL+fmBpKu+6IuIPETEzf5sfETNbWW9+RDydX296icvs0iLi3yNicdF2PruV9c7M799zI+ILpa6zK4uI70fECxExKyJuiIgBraznfryTdrRfRs5P8stnRcSRWdTZlUXEyIi4NyKez/8O/KcW1jk5ItYU/Rz5Sha1dmU7+v/vvrx7IuKgov1zZkSsjYhPN1vH/XgnRcSvIuLViHimqK9N33c7y/cKzyHaSRExFmgCfgH8c0pper5/HPB7YDKwF3AXcGBKqbHZ878HrEopfSf/Dz8wpXR5KT9DVxYRPwDWpJS+3sKy+cDElJIXVttJEfHvwPqU0n9uZ51KYA5wOrAImAa8P6X0XEmK7OIi4q3APSmlhoj4LkBL//fdj3dOW/bLfMD/R+BsYArwXymlKRmU22VFxHBgeEppRkT0BZ4A3tVsO59M7vfi27Opsuvb0f9/9+X2k//ZsRiYklJ6paj/ZNyPd0pEnAisB36TUjo037fD77ud6XuFI0Q7KaX0fEppdguLzgGuSyltSSm9DMwlF45aWu/X+ce/Bt7VIYV2QxERwHvJBU+V3mRgbkrppZRSHXAduf1ZbZBSuiOl1JBvTgVGZFlPN9KW/fIccr+oU0ppKjAg/wVfbZRSWppSmpF/vA54Htg726rKkvty+zkVmFcchrRrUkoPAKuadbfl+26n+V5hIGo/ewMLi9qLaPmXxZ4ppaWQ+wUD7FGC2rqLE4DlKaUXW1megDsi4omIuLiEdXUXl+UPwfhVK0Pbbd3HtWMfBW5rZZn78c5py37pvtuOImIMcATwWAuLj4mIpyLitog4pLSVdQs7+v/vvtx+zqf1P7C6H+++tnzf7TT7c1UWb9rZRcRdwLAWFn0ppXRja09roc/jEduojdv8/Wx/dOi4lNKSiNgDuDMiXsj/1UJsfxsDPwe+QW6f/QbwA3Jf2rd5iRae6z5epC37cUR8CWgAftvKy7gf75y27Jfuu+0kIvoAfwY+nVJa22zxDGB0Sml9/tCuvwAHlLjErm5H///dl9tBRNQA7wT+tYXF7sel02n2ZwNRC1JKp+3C0xYBI4vaI4AlLay3PCKGp5SW5oe5X92VGrubHW3ziKgC3gMctZ3XWJK/fzUibiA3FOsXyby27tcRcTVwcwuL2rqPl6027McfBt4OnJpaOYHT/XintWW/dN9tBxFRTS4M/TaldH3z5cUBKaV0a0T8LCKGeD5c27Xh/7/7cvs4C5iRUlrefIH7cbtpy/fdTrM/e8hc+7kJOD8iekTEPuT+mvB4K+t9OP/4w0BrI07a1mnACymlRS0tjIje+RN9iYjewFuBZ1paV2/W7Bj0d9PytpsGHBAR++T/unY+uf1ZbRARZwKXA+9MKW1sZR33453Xlv3yJuBD+Rm6jiY3McvSUhfaleXP4fwf4PmU0g9bWWdYfj0iYjK57xivla7Krq2N///dl9tHq0ecuB+3m7Z83+003yscIdpJEfFu4L+BocAtETEzpXRGSunZiPgj8By5w2E+uXWGuYj4JXBlfka67wB/jIiLgAXAeZl8kK7nTcf6RsRewC9TSmcDewI35H+GVQG/Syn9reRVdl3fi4gJ5Iaq5wOXwLbbOD872mXA7UAl8KuU0rMZ1dsVXQH0IHcYDMDUlNKl7se7p7X9MiIuzS+/EriV3Kxcc4GNwEeyqrcLOw64AHg63rj0wReBUVDYzn8HfCIiGoBNwPmtjYSqRS3+/3dfbl8R0YvcrGaXFPUVb2P3450UEb8HTgaGRMQi4Ku08n23s36vcNptSZIkSWXLQ+YkSZIklS0DkSRJkqSyZSCSJEmSVLYMRJIkSZLKloFIkiRJUtkyEEmSJEkqWwYiSVKXFhEVEfFARNzUrL9XRMyOiJ9nVZskqfMzEEmSurSUUhNwIXBKRHy0aNF3yV3g8p+zqEuS1DV4YVZJUreQv9r894DDgP3JXf385JTSQ5kWJknq1AxEkqRuIyJuB2qBMcB1KaXPZ1uRJKmzMxBJkrqNiNgHmJe/HZpS2pJxSZKkTs5ziCRJ3clHgU3ACGDfjGuRJHUBjhBJkrqFiJgEPAK8E/gEsCdwbEqpMdPCJEmdmiNEkqQuLyJ6Ar8Brkkp3QZcTG5iBc8hkiRtlyNEkqQuLyJ+BLwLODyltC7fdz7wa+ColNIzGZYnSerEDESSpC4tIk4E7gFOSynd12zZH8mdS3R0Sqkhg/IkSZ2cgUiSJElS2fIcIkmSJElly0AkSZIkqWwZiCRJkiSVLQORJEmSpLJlIJIkSZJUtgxEkiRJksqWgUiSJElS2TIQSZIkSSpb/x/nn6nyj9IAYgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "rcParams['figure.figsize'] = (14, 7)\n",
    "rcParams['axes.spines.top'] = False\n",
    "rcParams['axes.spines.right'] = False\n",
    "\n",
    "\n",
    "def sigmoid(x):\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "\n",
    "xs = np.arange(-10, 10, 0.1)\n",
    "ys = [sigmoid(x) for x in xs]\n",
    "\n",
    "plt.plot(xs, ys, c='#087E8B', lw=3)\n",
    "plt.title('Sigmoid function', size=20)\n",
    "plt.xlabel('X', size=14)\n",
    "plt.ylabel('Y (probability)', size=14)\n",
    "plt.savefig('images/003_LogisticRegression_sigmoid_function.png', dpi=300, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Essentially, sigmoid function returns a value between 0 and 1\n",
    "- This value is interpreted as a probability of the positive class\n",
    "- We can make classification from probabilities\n",
    "    - If the probability is greater than some threshold (commonly 0.5), assign the positive class\n",
    "    - If the probability is less than some threshold (commonly 0.5), assign the negative class\n",
    "    \n",
    "- As with linear regression, there are two parameters we need to optimize for - weights and bias\n",
    "- To perform the optimization, we'll need to declare **Cost function**\n",
    "    - We can't use Mean Squared Error as with linear regression\n",
    "    - We can, in theory, but it's not a good idea, because MSE was designed to evaluate regression models\n",
    "    - Instead, we can use the **Cross Entropy** function\n",
    "    \n",
    "- Formula of the **Binary Cross Entropy** function (2 classes):\n",
    "\n",
    "$$ \\large BCE = - \\frac{1}{n} \\sum_{i}^{n} y_i \\cdot log \\hat{y}_i + (1 - y_i) \\cdot log (1 - {\\hat{y}_i}) $$\n",
    "\n",
    "- The $ \\hat{y} $ in the formula can be re-written as a sigmoid function for a given input\n",
    "- A brief introduction to BCE is made later\n",
    "- Next you'll need to use this cost function in gradient descent to iteratively update weights and bias\n",
    "- To do so, you'll have to calculate partial derivatives of the cross entropy function wrt. weights and bias parameters:\n",
    "\n",
    "$$ \\large \\partial_w = \\frac{1}{N}\\sum 2x_i(\\hat{y} - y_i)) $$\n",
    "<br>\n",
    "$$ \\large \\partial_b = \\frac{1}{N}\\sum 2(\\hat{y} - y_i)) $$\n",
    "\n",
    "- This $ 2 $ can be ommited, or you can leave it - it's not important\n",
    "- Next, we're updating the existing weights and bias according to the following formulas:\n",
    "\n",
    "$$ \\large w = w - \\alpha \\cdot \\partial_w $$\n",
    "<br>\n",
    "$$ \\large b = b - \\alpha \\cdot \\partial_b $$\n",
    "\n",
    "- Where $ \\alpha $ is the learning rate\n",
    "- This process is then repeated for a predefined number of iterations\n",
    "- Before we see this in action, let's make a brief introduction to the Binary Cross Entropy function\n",
    "\n",
    "<br>\n",
    "\n",
    "## A Brief Introduction to Binary Cross Entropy\n",
    "- Binary cross entropy function is a common cost function for evaluating binary classification machine learning algorithms\n",
    "- Commonly refered to as *log loss*\n",
    "- It's value can be calculated from the previously mentioned formula\n",
    "- This cost function \"punishes\" wrong predictions much more than it \"rewards\" good predictions\n",
    "\n",
    "**Example 1: Calculating BCE for a correct prediction**\n",
    "- Let's say your model predicts the positive class with 90% probability\n",
    "- This means the model is only 10% confident the negative class should be predicted\n",
    "- What's the cross entropy value? \n",
    "\n",
    "![](images/003_logreg_crossentropy_small.jpg)\n",
    "\n",
    "- As you can see, the value of cross entropy (or *loss*) is rather small\n",
    "\n",
    "**Example 2: Calculating BCE for a incorrect prediction**\n",
    "- Let's say your model predicts the positive class with 10% probability\n",
    "- This means the model is 90% confident the negative class should be predicted\n",
    "- What's the cross entropy value? \n",
    "\n",
    "![](images/003_logreg_crossentropy_big.jpg)\n",
    "\n",
    "- As you can see, the loss here is quite big\n",
    "- We can calculate both with Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_cross_entropy(y, y_hat):\n",
    "    def safe_log(x): return 0 if x == 0 else np.log(x)\n",
    "    \n",
    "    total = 0\n",
    "    for curr_y, curr_y_hat in zip(y, y_hat):\n",
    "        total += (curr_y * safe_log(curr_y_hat) + (1 - curr_y) * safe_log(1 - curr_y_hat))\n",
    "    return - total / len(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We need the `safe_log()` function inside because `log(0)` returns infinity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-3-f6e7c0610b57>:1: RuntimeWarning: divide by zero encountered in log\n",
      "  np.log(0)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-inf"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's evalutate for **Example 1**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10536051565782628"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_cross_entropy(y=[1, 0], y_hat=[0.9, 0.1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- And now for **Example 2**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.302585092994046"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_cross_entropy(y=[1, 0], y_hat=[0.1, 0.9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- You now know everything needed to implement Logistic regression algorithm from scratch\n",
    "\n",
    "<br><br>\n",
    "\n",
    "## Implementation\n",
    "- The `LogisticRegression` class is written to follow the familiar Scikit-Learn syntax\n",
    "- The coefficients are set to None at the start - `__init__()` method\n",
    "- The `fit()` method calculates the coefficients\n",
    "- The `predict()` method essentially implements the line equation passed through a sigmoid function\n",
    "- The `_binary_cross_entropy()` private function is used to calculate loss at every iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression:\n",
    "    '''\n",
    "    A class which implements logistic regression model with gradient descent.\n",
    "    '''\n",
    "    def __init__(self, learning_rate=0.1, n_iterations=1000):\n",
    "        self.learning_rate = learning_rate\n",
    "        self.n_iterations = n_iterations\n",
    "        self.weights, self.bias = None, None\n",
    "        \n",
    "    @staticmethod\n",
    "    def _sigmoid(x):\n",
    "        '''\n",
    "        Private method, used to pass results of the line equation through the sigmoid function.\n",
    "        \n",
    "        :param x: float, prediction made by the line equation\n",
    "        :return: float\n",
    "        '''\n",
    "        return 1 / (1 + np.exp(-x))\n",
    "    \n",
    "    @staticmethod\n",
    "    def _binary_cross_entropy(y, y_hat):\n",
    "        '''\n",
    "        Private method, used to calculate binary cross entropy value between actual classes \n",
    "        and predicted probabilities.\n",
    "        \n",
    "        :param y: array, true class labels\n",
    "        :param y_hat: array, predicted probabilities\n",
    "        :return: float\n",
    "        '''\n",
    "        def safe_log(x): \n",
    "            return 0 if x == 0 else np.log(x)\n",
    "        total = 0\n",
    "        for curr_y, curr_y_hat in zip(y, y_hat):\n",
    "            total += (curr_y * safe_log(curr_y_hat) + (1 - curr_y) * safe_log(1 - curr_y_hat))\n",
    "        return - total / len(y)\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        '''\n",
    "        Used to calculate the coefficient of the logistic regression model.\n",
    "        \n",
    "        :param X: array, features\n",
    "        :param y: array, true values\n",
    "        :return: None\n",
    "        '''\n",
    "        # 1. Initialize coefficients\n",
    "        self.weights = np.zeros(X.shape[1])\n",
    "        self.bias = 0\n",
    "        \n",
    "        # 2. Perform gradient descent\n",
    "        for i in range(self.n_iterations):\n",
    "            linear_pred = np.dot(X, self.weights) + self.bias\n",
    "            probability = self._sigmoid(linear_pred)\n",
    "            \n",
    "            # Calculate derivatives\n",
    "            partial_w = (1 / X.shape[0]) * (2 * np.dot(X.T, (probability - y)))\n",
    "            partial_d = (1 / X.shape[0]) * (2 * np.sum(probability - y))\n",
    "            \n",
    "            # Update the coefficients\n",
    "            self.weights -= self.learning_rate * partial_w\n",
    "            self.bias -= self.learning_rate * partial_d\n",
    "            \n",
    "    def predict_proba(self, X):\n",
    "        '''\n",
    "        Calculates prediction probabilities for a given threshold using the line equation \n",
    "        passed through the sigmoid function.\n",
    "        \n",
    "        :param X: array, features\n",
    "        :return: array, prediction probabilities\n",
    "        '''\n",
    "        linear_pred = np.dot(X, self.weights) + self.bias\n",
    "        return self._sigmoid(linear_pred)\n",
    "    \n",
    "    def predict(self, X, threshold=0.5):\n",
    "        '''\n",
    "        Makes predictions using the line equation passed through the sigmoid function.\n",
    "        \n",
    "        :param X: array, features\n",
    "        :param threshold: float, classification threshold\n",
    "        :return: array, predictions\n",
    "        '''\n",
    "        probabilities = self.predict_proba(X)\n",
    "        return [1 if i > threshold else 0 for i in probabilities]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "## Testing\n",
    "- We'll use the *breast cancer* dataset from Scikit-Learn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_breast_cancer\n",
    "\n",
    "data = load_breast_cancer()\n",
    "X = data.data\n",
    "y = data.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The below code cell applies train/test split to the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- You can now initialize and train the model, and afterwards make predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-6-822322b64438>:18: RuntimeWarning: overflow encountered in exp\n",
      "  return 1 / (1 + np.exp(-x))\n"
     ]
    }
   ],
   "source": [
    "model = LogisticRegression()\n",
    "model.fit(X_train, y_train)\n",
    "preds = model.predict(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- These are the \"optimal\" weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 6.48458109e+01,  8.46348283e+01,  3.76652184e+02,  1.54054003e+02,\n",
       "        5.83159961e-01, -2.99918435e-01, -1.18719824e+00, -5.08264671e-01,\n",
       "        1.09780244e+00,  4.70543004e-01,  2.55378613e-01,  6.20056326e+00,\n",
       "       -1.54946182e+00, -1.60920898e+02,  3.11666417e-02, -6.77737261e-02,\n",
       "       -1.32537252e-01, -1.64448048e-02,  9.92337789e-02,  1.02283101e-02,\n",
       "        6.80812543e+01,  1.05097039e+02,  3.79967070e+02, -2.23235180e+02,\n",
       "        7.12011129e-01, -1.31516705e+00, -2.55931308e+00, -5.99329656e-01,\n",
       "        1.40482786e+00,  4.28992752e-01])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- And this is the \"optimal\" bias:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.486537143743831"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.bias"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's evaluate the model with two metrics - accuracy score and the confusion matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score, confusion_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9473684210526315"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(y_test, preds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[43  0]\n",
      " [ 6 65]]\n"
     ]
    }
   ],
   "source": [
    "print(confusion_matrix(y_test, preds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The model seems to be performing good\n",
    "- Has more false negatives than false positives\n",
    "\n",
    "<br><br>\n",
    "\n",
    "## Threshold Optimization\n",
    "- There's no guarantee that 0.5 is the best classification threshold for every classification problem\n",
    "- We can change the threshold by altering the `threshold` parameter of the `predict()` method:\n",
    "- Let's optimize the threshold for accuracy - you can choose any metric you want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-6-822322b64438>:18: RuntimeWarning: overflow encountered in exp\n",
      "  return 1 / (1 + np.exp(-x))\n"
     ]
    }
   ],
   "source": [
    "evals = []\n",
    "\n",
    "for thresh in np.arange(0, 1.01, 0.01):\n",
    "    preds = model.predict(X_test, threshold=thresh)\n",
    "    acc = accuracy_score(y_test, preds)\n",
    "    evals.append({'Threshold': thresh, 'Accuracy': acc})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0QAAAHCCAYAAADPbwnCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+/klEQVR4nO3deZhcVbWw8Xd1d7XMAZkMBCXKDBmAEECmYCQGQZAZRPlAuYiCInzKcPUqH15ERRS8gIiIiBcSHBlkMMosyBAkzEIiBNOAEAIBMpHuZH9/VHXnVKU6qe70VDnv73nq6a5zTp2zuvoEavdae+1IKSFJkiRJedTQ3wFIkiRJUn9xQCRJkiQptxwQSZIkScotB0SSJEmScssBkSRJkqTcckAkSZIkKbccEElSN0XEXRFRF2sXRMSxEZEi4tiBdI2IuKr0mk17K67e1Efv65jSNc7urWt0VSmeu/o7DknqCQ6IJOVe6cNdVx7H9nfMeRcR60TENyPioYh4MyLejYgZEfHriNinB68z4AYjfSEipkfE9P6OQ5L6QlN/ByBJA8D/q7LtK8Ag4CJgdsW+Kb0bjpYlIvYEfgesBzwDXAO8A2wO7AccFhH/CxyfUnq3l8P5A/AA8EovXuMhYGvg9V68RldtDczr7yAkqSc4IJKUeymlsyu3lbJAg4ALU0rT+zgkdSIitgFuAVYFvgxcnFJKmf2bANcDnwYWAp/rzXhSSm8Bb/XyNeYB/+jNa3RVSmlAxSNJK8KSOUlaQRHRFBH/GRFTM6Vb34uI5irHptLco/dFxBUR8VJELMqW4UXEzhHx24j4d0QsLJ3vpxGxUZXzfTAiLo+IaRExPyLeiIgnIuKyiFi3k3j3LsXwTkS8HRE3R8TWnRw7OCIuKZVQLYyImRHx+4jYsYvv0Ucj4t6ImFuK8fqI2Kor5yj5MbA68P2U0v9kB0MAKaUZwP7Am8BnI2K3ijja3/+NIuJXEfFa6X17JCI+VXHsVcCdpaffqiibHFM6puocovaSs4hYIyJ+VPodzo+IKRHxydIx2ftmQUT8MyJOrvLeLVW2FxFnL6+0M3Nsc0ScHBG3RMSLpXv0jYj4S0TsW+1awAeAD1Sc86rK97FKrIMi4ryIeLb0M70ZEX+KiI8u6+eKiJGl+3B2RMyLiLsj4sOVr5Gk3mCGSJJW3LXAHsCtwNvAx4HTgQ2A46oc/16KZVZzgN8Di4FXASLiOOBnwLvAjcAMiqVgxwOfiIhdUkr/Kh07GHgYWIti1uR3wCrAUOAzwMXArIpr7w8cWIr1MmCbUrw7RcQ2KaWOsqyIGAr8FdgIuAOYAGwCHAbsFxGHpJT+uLw3JyIOBa6jmLG5jmJ52e7A34DHl/f6injGlt6b73d2XErplYi4AvgacAJwX8Uh6wD3UyyF/AWwNnA4cE1EbJxSOr903PWlr/8HuBu4K3OO6TWEXAD+TPH3fQPQDBwF/C4ixgFfBHam+Lt4l+L7+j8RMTOldN1yzn1XJ9s3AT4LzM9sey/F0s/7S/HMBAYDnwBuiYj/SCldkfm5/h/FklGACzPnmbKsgCJibYrv9TYU78sLKZY1Hg5MiogvpJR+WuWloyj+e/kbcAXwfuAQ4PaIGJlSenZZ15WkFZZS8uHDhw8fFQ+KHwwTsOkyjrmrdMwjwHsz21cHpgGLgPdVvCaVHlcDTRX7tqA4aJgGbFyx7yOl8/0hs+1LpXOdUiW21YFVM8+PLR3bBoytOPa80r7TK7b/qbT96xXbP1w6zyxgjSrXODazbY3Sca3AqIrz/CjzfnT6PmeO/0zp2PtqOHaf0rHTOnn/fw00ZLYPBd4ovf8fzGwfUzr+7E6us9TPXHH/3AS8J7N9j9L2NygOGtbO7Ptg6fqPVpxrmTFkjluL4gBzEXBwZvt7gCFVjh8EPFmKZdUq8U9fxrUScFfFtp+Wtv8UiMz2zSmWFb6b/T1nfq5q79/nS9sv7e6/YR8+fPio9WHJnCStuDNSSm+0P0kpzaU40b+B4l+/Ky0EvppSaqvY/gWKWYVTUkovZXeklO6gmDH6RESsWfG6+RXPSSnNTSkttR2YmFK6vWLb5aWvo9s3RMQQYBzwLyqyMSml+ylmi94LHFzlGlkHlo67NqU0uWLf2XRt/s3g0tcZNRzbfsxSZYYUBwxnpJQWt29IKb1AsRyvQHHg1VO+kjKNHVJK9wIvUMxSnZFSmp3Z9zzFDMuwiGjsykUiogn4DTAM+FpK6feZ876bUmqpfE0qzn+6shTLTl25XpXrFyjO25oDnJVS6ijZSylNpfjeNgPHVHn5fSmlqyq2XUlx0D166cMlqWc5IJKkFVf5QR+WfCBfp8q+6Sml16ps37X0da/SvIqyB8USvEaKmSQoDpDmAJdExO8i4oSI2DYiogdi3b709d6UUmuV19xRcVxndih9vbtyR+kD+ZTlvD6r/eeqZe2nZR37r9IAqNJdpa/L+5lqNTul9M8q218ufX2kyr6XKP6O39fFa/2E4gD20pTSDyt3lu6LqyLi+dJcpvZ5RheUDtm4i9ertBWwGvBY9o8DGcu6X5a6J0v33KtU//cjST3KOUSStIKyf+XPaM/+VPtL/787OVV7E4SvLeeSa5Su+2JEjKaYaRnPkmzNjIj4QUrpx1Veu1SsKaW20hgqG+ug0tfO2km3b197ObG2n+fVTvZ39l4s65rvr+HYIRWvyVpeLIM62d9VnWW/2qBjQFh1H8VMVU0i4iyKc8xupth5r3L/LhQHJE3A7RQH0m9TnLs2kmIW7z21Xq8TK3K/zO7kNW1U//cjST3KAZEk9b3OMhztH5AHpZTerulEKT0DHFEqmRoBfJTi3KKLImJuSunn3YyxPZbOMhWDK45b3nk27GR/VzIhfy193TEi1u5kINquvatZZUOFWmLp1TbaPSkijgDOBR4FjkwpLapy2DcotinfO6V0V8Xrz6I4IFpRPXW/SFKfs2ROkgaOB0pf9+jqC1NKbSmlR1JK36PYyQzgkysQy6Olr7uXBluV9i59/ftyztO+f6/KHRExiGKGoialOTZ3UsxmdJpFi4gNKWZMYMn8qKz3R8SmVbaPKX19NLOtfYAx4DIVpbbUV1Ess9s/pTSnk0M3A96oHAyVLPV7KVlE137mZyku1DoyIqqVudV6v0hSn3NAJEkDx8UUu7H9KCK2qNxZWk9mj8zz0aUP/5Xat83rbiClSfh/BjZlSQvm9uvuDHyK4lo/f1jOqW4oHfepiKhsMHE2XS9PO4Xiz3VGRHyhcmdEbAz8kWIjh1+klKpliBqB70VEQ+Z1QymWm7UB/5s5tr1teS1len0mIjaj+N62AvullF5exuHTgfdGxPCKc3wO+Fgnr5kFrB8Rq9YST0ppIcVGImsA51Rc50MU39tW4Fe1nE+S+pIlc5I0QKSU/hERn6XYYeupiLgNeI7ifJL3U8wczaQ4gR2Kg5KTIuJuiq263wQ+RHF9mXcpX0OmO06kWHJ2fmndnMksWYdoMXBcSumd5fxMcyLiBIrrD90bEdl1iLYD7gH2rDWglNITEbE/8Fvg0og4iWLW6B2KmZD9KE7uv4Zi175qHqe4/s8jETGJ4qDsCIrzW06vaITwLMUMzJERsZBi170E/Cql9GKtcfeCH1Nc4+cO4OCIWKrbX0rp7NK3F1Ic+Pw1In5NsWxtFMXfwW+BQ6uc/3aKnedui4h7KN5Pj6WUblpGTGdSvEdPjoidKP5e2tchWhM4uZNmFpLUrxwQSdIAklL634h4DPi/FMuMxgFzKXYm+y3FgUW7CRTLxz5MsZvbqhQ/vE8ELkgpPbmCsTxfyup8g+LirWMoTsa/DTg3pfRwjef5bUSMB75F8cPxuxQHQrtS/BBd84CodL47Sxm0L1FcaPYYigvSzqS4QO3PUkqTlnGKN4F9KbYTP47i+j1PAz9IKV1bca1FEXEQ8F2WfLAPivOZ+nNAtFrp60dKj2rOBkgp3RYRn6D4ezyCYjncQxTvrw9SfUD03xQHiJ8AdqOYVfslxXWVqkopvRERuwJnUWzwcRrFlvAPAecv53ciSf0mMksFSJK0Uiu1mr47pTSmv2ORJA0MziGSJEmSlFsOiCRJkiTllgMiSZIkSbnlHCJJkiRJuVX3XebGjx+fbrvttv4OQ5IkSdLAFZ3tqPuSuddff72/Q5AkSZJUp+p+QCRJkiRJ3eWASJIkSVJuOSCSJEmSlFsOiCRJkiTllgMiSZIkSbnlgEiSJElSbjkgkiRJkpRbDogkSZIk5VZTfwdQ71JK3Dt9BhfdP5lHX3mVea1trFZoYvvBG/KV3XZi9w8MIaLThXElSZIk9SMHRCtg0tQXOOnGScxesIC5C1tJpe1vLXiXf7/zPPdMn8E6q67CJQeMY5/NhvZrrJIkSZKWZslcN1399yc4cuINtLz9DnMyg6F2CZizsJUZb73DERNu4Oq/P9EfYUqSJElaBgdE3TBp6gt85ebbmd/WVtPx89va+MrNtzNp6gu9HJkkSZKkrnBA1EUpJU66cVLNg6F289vaOPmmSaRUmUuSJEmS1F8cEHXRvdNnMHvBgm699o35C/jriy09HJEkSZKk7nJA1EUX3T+ZuQtbu/XaeQtbufC+yT0ckSRJkqTuckDURY++8upSDRRqlYApr7zak+FIkiRJWgEOiLpoXmvX5g4t/fruZZckSZIk9TwHRF20WmHFlm5arVDooUgkSZIkrSgHRF20/eANiW6+NoCRgzfsyXAkSZIkrYAVS3fk0CkfHsU902cwpxuNFRoimL1gAYdPuL7nA5MkSZL6UEME+235IT6z/Xb9HcoKcUDURXtsuglrr7JKtwZEi1Li/n+91AtRSZIkSX3vhmemstOQwWy1/rr9HUq3WTLXRRHBpQeOY9Umx5KSJEnStFlv9ncIK8RP9d2wz2ZDuXC/sXzl5tuZ37b8rnOrNDZy2u6jGT54gz6ITpIkSepdP/zrQzzU8goArYsX93M0K8YBUTcds8Mw3rfmGpx80yTemL+AeQtby9YnCmC15gLvXXUVLjlgHPtsNrS/QpUkSZJ61G+e+MeSAdEiB0S5NW7zoTx76gn89cUWLrzvYaa88hrzWltZrVBg5OANOXW3Uez2gSFEdLcvnSRJkjTwNDc2dnzfumhRP0ay4hwQraCIYI9NN2GPTTfp71AkSZKkPlFoXNKKoN4zRDZVkCRJktQlhYbMgGhxfWeIHBBJkiRJ6pJCpmRuoRkiSZIkSXlSXjJnhkiSJElSjpSXzJkhkiRJkpQjK1OXOQdEkiRJkrrELnOSJEmScivbVMGSOUmSJEm5kp1DtNCSOUmSJEl5UpYhsmROkiRJUp40Z+YQtVkyJ0mSJClPCg3ZhVktmZMkSZKUI3aZkyRJkpRb5V3mzBBJkiRJypFshmihGSJJkiRJeZJtu93mHCJJkiRJedLswqySJEmS8qq8ZM4MkSRJkqQcyZbM2WVOkiRJUq4ULJmTJEmSlFfl6xBZMleziBgfEc9GxLSIOLPK/nUi4g8R8XhEPBQR2/VlfJIkSZKWr9CQyRBZMlebiGgELgH2BbYBjoqIbSoO+09gSkppOHAMcFFfxSdJkiSpNs3ZDJElczUbDUxLKT2fUloITAQOrDhmG+B2gJTSP4BNI2LDPoxRkiRJ0nJk5xDZZa52GwMzMs9bStuyHgMOBoiI0cAHgCGVJ4qIEyJickRMnjlzZi+FK0mSJKmaJrvMdUtU2ZYqnn8XWCcipgBfAh4F2pZ6UUqXp5RGpZRGrb/++j0eqCRJkqTOZRdmbVtc3xmipj68VguwSeb5EODl7AEppbeB4wAiIoAXSg9JkiRJA0T5wqxmiGr1MLB5RAyNiGbgSODG7AERsXZpH8DxwD2lQZIkSZKkAaJ8YVYzRDVJKbVFxMnAn4BG4MqU0lMRcWJp/2XA1sDVEbEIeBr4XF/FJ0mSJKk2zSvRwqx9WTJHSukW4JaKbZdlvv8bsHlfxiRJkiSpa8pL5uo7Q9SnC7NKkiRJqn+NDQ1EqWVaSrCojrNEDogkSZIkdVlZ2VwdN1ZwQCRJkiSpy7KNFeq5bM4BkSRJkqQuK6wkjRUcEEmSJEnqsubGbOttB0SSJEmScqSpYUmGyJI5SZIkSbmSbb3dZsmcJEmSpDwp7zJnhkiSJElSjpR3mTNDJEmSJClHsiVzrYvNEEmSJEnKkYILs0qSJEnKq2zJnHOIJEmSJOVKecmcGSJJkiRJOdJsyZwkSZKkvCrvMmfJnCRJkqQcacpmiCyZkyRJkpQnzdk5RJbMSZIkScqTQkN2DpElc5IkSZJypGCGSJIkSVJelXWZW2yGSJIkSVKONJV1mTNDJEmSJClHykvmzBBJkiRJypFm225LkiRJyqvswqxmiCRJkiTlil3mJEmSJOVWwZI5SZIkSXlVKOsyZ8mcJEmSpBwpyxBZMidJkiQpT5ozc4jaLJmTJEmSlCeFhiUZIkvmJEmSJOWKXeYkSZIk5VZ2DpEZIkmSJEm5UpYhcg6RJEmSpDzJtt1uM0MkSZIkKU+ay0rmzBBJkiRJypHykjkzRJIkSZJyJFsyZ5c5SZIkSblilzlJkiRJuZUtmWuzy5wkSZKkPCk0LMkQWTInSZIkKVeaMxkiS+YkSZIk5Up2DpELs0qSJEnKlSa7zEmSJEnKq+zCrK2WzEmSJEnKk/KFWc0QSZIkScqR8oVZzRBJkiRJypHmsoVZzRBJkiRJypHykjkzRJIkSZJypLGhgYji9ynBojqdR+SASJIkSVK3lJfN1WeWyAGRJEmSpG4prARrEfXpgCgixkfEsxExLSLOrLJ/UETcFBGPRcRTEXFcX8YnSZIkqXaF7FpElswtW0Q0ApcA+wLbAEdFxDYVh50EPJ1SGgGMAS6IiOa+ilGSJElS7ZozjRUsmVu+0cC0lNLzKaWFwETgwIpjErBmRASwBvAG0NaHMUqSJEmqUVNDJkNkydxybQzMyDxvKW3LuhjYGngZeAI4JaW01DsbESdExOSImDxz5szeileSJEnSMmRbb7dZMrdcUWVbqnj+MWAKsBEwErg4ItZa6kUpXZ5SGpVSGrX++uv3dJySJEmSamCXua5pATbJPB9CMROUdRzw+1Q0DXgB2KqP4pMkSZLUBXaZ65qHgc0jYmipUcKRwI0Vx/wLGAsQERsCWwLP92GMkiRJkmqULZlrXVyfGaKmvrpQSqktIk4G/gQ0AlemlJ6KiBNL+y8Dvg1cFRFPUCyxOyOl9HpfxShJkiSpdoWykrn6zBD12YAIIKV0C3BLxbbLMt+/DIzry5gkSZIkdU+2ZK7NOUSSJEmS8qS8ZK4+M0QOiCRJkiR1i13mJEmSJOWWXeYkSZIk5VZTJkNkyZwkSZKkXGnOzCGyZE6SJElSrhQalmSI2iyZkyRJkpQnZV3mHBBJkiRJyhO7zEmSJEnKraYG1yGSJEmSlFPlJXNmiCRJkiTlSLZkzjlEkiRJknKlbGHWxWaIJEmSJOWIXeYkSZIk5VbBLnOSJEmS8qpglzlJkiRJeVWwqYIkSZKkvGq27bYkSZKkvCo0ZDJElsxJkiRJyhO7zEmSJEnKLbvMSZIkScqtsgyRJXOSJEmS8iTbdrvNDJEkSZKkPGkuK5kzQyRJkiQpR8pL5swQSZIkScqRbMmcXeYkSZIk5Ypd5iRJkiTlVrZkrs0uc5IkSZLypNCwJENkyZwkSZKkXGnOZIgsmZMkSZKUK9k5RC7MKkmSJClXmuwyJ0mSJCmvsguztloyJ0mSJClPyhdmNUMkSZIkKUfKF2Y1QyRJkiQpR5rLFmY1QyRJkiQpR8pL5swQSZIkScqRxoYGIorfpwSL6nAekQMiSZIkSd1WXjZXf1kiB0SSJEmSuq1Q52sROSCSJEmS1G2F7FpElsxJkiRJypPmTGMFS+YkSZIk5UpTQyZDZMmcJEmSpDzJtt5us2ROkiRJUp7kostcRHwyIhqXf6QkSZKkPMlLl7lrgJci4nsRsWVvBiRJkiSpfmRL5loXr6QZIuB9wLeAvYCnI+KvEXFcRKzee6FJkiRJGugKZSVzK2mGKKX0TkrppymlXYBhwIPAecArEfGziNilN4OUJEmSNDBlS+baVtY5RFkppaeBHwGXA83AEcC9EfFgRAzv4fgkSZIkDWDlJXMraYYIICIKEXF4RNwGvAB8BDgR2BD4APAccF2vRClJkiRpQKr3LnNNtRwUEf8DHAUk4FfAaaVMUbv5EfF1YHqPRyhJkiRpwMpLl7ltgJOBjVNKlYOhdi8Dey/rJBExPiKejYhpEXFmlf1fi4gppceTEbEoIt5bY4ySJEmS+lhTJkNUjyVzNWWIUkpjazimDbi7s/2ldYwuAfYBWoCHI+LG7OAqpXQ+cH7p+E8Ap6aU3qglRkmSJEl9rzkzh6geS+ZqXZj13Ig4scr2EyPi2zVeazQwLaX0fEppITAROHAZxx8FTKjx3JIkSZL6QaFhSYaobSUumfsM8GiV7Y8Ax9R4jo2BGZnnLaVtS4mI1YDxwO862X9CREyOiMkzZ86s8fKSJEmSelo2Q7QyzyHaAKg28phFsctcLaLKttTJsZ8A7uusXC6ldHlKaVRKadT6669f4+UlSZIk9bRCnXeZq3VA9C9gjyrb96SY6alFC7BJ5vkQio0YqjkSy+UkSZKkAa+pob7XIaqpqQLwU+BHEdEM3FHaNhY4D/hejed4GNg8IoYCL1Ec9Hyq8qCIGATsBXy6xvNKkiRJ6iflJXP1lyGqtcvcBRGxHvBjoLm0eSFwUUrp+zWeoy0iTgb+BDQCV6aUnmpv1pBSuqx06EHApJTS3C78HJIkSZL6QbZkrh7nENWaISKldFZE/DfFNYkCeDqlNKcrF0sp3QLcUrHtsornVwFXdeW8kiRJkvpH2cKsi1fSDFG7Utbm4V6KRZIkSVKdac5Lhigi9qa4NtD7WVI2B0BK6SM9HJckSZKkOtCUk4VZjwVuBdYExlBswb0OsAPwdC/FJkmSJGmAK9R5l7la225/FTg5pXQU0AqclVLaHvhfoEvziCRJkiStPOq9ZK7WAdEHgb+Uvn8XWKP0/cXAsT0ckyRJkqQ6Uajzttu1DohmUSyXg+IaQtuVvl8XWLWng5IkSZJUHwoNmQxRHZbM1dpU4V5gHPAE8GvgxxGxD8XFWf/cS7FJkiRJGuCam7IZopV3QHQysErp+/OANmA3ioOj/+6FuCRJkiTVgaZMhqgeu8wtd0AUEU3AkcD1ACmlxcD3ejcsSZIkSfWgbA5RHZbMLXcOUUqpDTgfKPR+OJIkSZLqSbbLXFsdZohqbarwALBjbwYiSZIkqf5k1yFauBLPIfoZ8IOIeD/wCDA3uzOl9PeeDkySJEnSwFdeMld/GaJaB0TXlr7+sMq+BDRW2S5JkiRpJVeo84VZax0QDe3VKCRJkiTVpfKSuZU0Q5RSerG3A5EkSZJUf7Ilc2112GWupgFRRBy8rP0ppd/3TDiSJEmS6klzTkrmftvJ9lT66hwiSZIkKYfqvWSuprbbKaWG7ANoBnYG7gX27M0AJUmSJA1cZU0V6rBkrtZ1iMqklNpSSg8D/wlc2rMhSZIkSaoXZW2367BkrlsDoozZwId6IA5JkiRJdajQkJ1DVH8lc7U2VdihchMwGDgDeLSng5IkSZJUH8oXZq2/DFGtTRUmU2ygEBXbHwCO69GIJEmSJNWN8i5zK2mGiKUXZl0MzEwpLejheCRJkiTVkfIucytphsiFWSVJkiRVU14yV38ZopqaKkTEuRFxYpXtJ0bEt3s+LEmSJEn1oLGhgYYozqxJCRbV2TyiWrvMfYbqzRMeAY7puXAkSZIk1ZtslqjeFmetdUC0ATCzyvZZwIY9F44kSZKkepOdR1RvaxHVOiD6F7BHle17Ai09F44kSZKkepPtNFdvGaJau8z9FPhRRDQDd5S2jQXOA77XG4FJkiRJqg/1vBZRrV3mLoiI9YAfA82lzQuBi1JK3++t4CRJkiQNfE0N2bWIVsIBEUBK6ayI+G9gG4oLtD6dUprTa5FJkiRJqgvN2QzRylgyFxHvA5pSSi3Aw5ntQ4DWlNKrvRSfJEmSpAGukJlDVG8lc7U2VfgVsG+V7R8r7ZMkSZKUU3noMrcTcE+V7fcCo3ouHEmSJEn1pp67zNU6IGoC3lNl+yqdbJckSZKUE0113GWu1gHRg8AXqmw/icycIkmSJEn5ky2Za6uzDFGtXea+DtwRESOA20vbPgLsQHE9IkmSJEk5VV4ytxJmiFJKDwC7Ai8ABwOHAM+Xtq3Wa9FJkiRJGvDKF2ZdOTNEpJQeA46GjnbbxwF/AN4PNC7jpZIkSZJWYnnoMkdENEbEQRFxM8VM0SeBnwCb9VJskiRJkupAoY67zC03QxQRWwLHA8cAc4FrKa4/9JmU0tO9G54kSZKkga6wsnaZi4h7gQeAtYHDU0ofTCl9A0h9EJskSZKkOlBoWJIhaquzkrnlZYh2BS4BfpZSerIP4pEkSZJUZ5ozGaJ6K5lb3hyiURQHTfdGxKMRcWpEvK8P4pIkSZJUJ7JziFaqpgoppSkppZOAwcAPgQOBGaXX7RcR6/R+iJIkSZIGsqaGlXQOUbuU0oKU0q9SSmOArYHzgVOBf0fErb0YnyRJkqQBbmUumVtKSmlaSulMYBPgcGBhj0clSZIkqW5kS+ZWtqYKnUopLQJuKD0kSZIk5VTZwqyLV/IMkSRJkiRlNZctzFpfGSIHRJIkSZJWSFN2YdaVfQ7RioiI8RHxbERMi4gzOzlmTERMiYinIuLuvoxPkiRJUtcV6rjLXLfnEHVVRDRSXOR1H6AFeDgibkwpPZ05Zm3gUmB8SulfEbFBX8UnSZIkqXvKS+bMEHVmNDAtpfR8SmkhMJHiukZZnwJ+n1L6F0BK6bU+jE+SJElSNxQyJXP11mWuLwdEG1Nc1LVdS2lb1hbAOhFxV0Q8EhHHVDtRRJwQEZMjYvLMmTN7KVxJkiRJtSg0LMkQ1VvJXF8OiKLKtlTxvAnYEdgP+BjwXxGxxVIvSunylNKolNKo9ddfv+cjlSRJklSz5qb6XZi1z+YQUcwIbZJ5PgR4ucoxr6eU5gJzI+IeYATwXN+EKEmSJKmrmrIZIkvmOvUwsHlEDI2IZuBI4MaKY24A9oiIpohYDdgZeKYPY5QkSZLURdk5RPVWMtdnGaKUUltEnAz8CWgErkwpPRURJ5b2X5ZSeiYibgMeBxYDV6SUnuyrGCVJkiR1XbbLXL2tQ9SXJXOklG4BbqnYdlnF8/OB8/syLkmSJEndV7YOkSVzkiRJkvKkvGSuvjJEDogkSZIkrZBC2cKsZogkSZIk5Uh5yZwZIkmSJEk5ki2Za6uzLnMOiCRJkiStkOaykjkzRJIkSZJyxC5zkiRJknIr21Sh3hZmdUAkSZIkaYVk5xBZMidJkiQpVwoNSzJEbZbMSZIkScqT8oVZHRBJkiRJyhG7zEmSJEnKLbvMSZIkScqt8pK5RaSU+jGarnFAJEmSJGmFNDY00BABQEqwaLEDIkmSJEk5UpklqhcOiCRJkiStsHqdR+SASJIkSdIKq9dOcw6IJEmSJK2wel2LyAGRJEmSpBXW1LAkQ2TJnCRJkqRcac5miCyZkyRJkpQnhcwcIkvmJEmSJOWKXeYkSZIk5ZZd5iRJkiTlVpNd5iRJkiTlVbZkrs0MkSRJkqQ8KS+ZM0MkSZIkKUfKF2Y1QyRJkiQpR+wyJ0mSJCm3CnaZkyRJkpRXBbvMSZIkScqrQsOSDFGbJXOSJEmS8qQ5kyGyZE6SJElSrmTnENlUQZIkSVKuNDU4h0iSJElSTlkyJ0mSJCm3siVzNlWQJEmSlCtlC7MuNkMkSZIkKUeayxZmNUMkSZIkKUeasguzOodIkiRJUp4U7DInSZIkKa/KS+bMEEmSJEnKkUKmZM4uc5IkSZJypdCwJENkyZwkSZKkXGlucmFWSZIkSTnVlM0QWTInSZIkKU+yc4gsmZMkSZKUK9kuc65DJEmSJClXytYhsmROkiRJUp6Ul8yZIZIkSZKUI4WyhVnNEEmSJEnKkfKSOTNEVUXE+Ih4NiKmRcSZVfaPiYi3ImJK6fHNvoxPkiRJUvdkS+ba6qjLXFNfXSgiGoFLgH2AFuDhiLgxpfR0xaH3ppT276u4JEmSJK245rKSOTNE1YwGpqWUnk8pLQQmAgf24fUlSZIk9RK7zC3fxsCMzPOW0rZKu0bEYxFxa0RsW+1EEXFCREyOiMkzZ87sjVglSZIkdUG2qYILs1YXVbaliud/Bz6QUhoB/A9wfbUTpZQuTymNSimNWn/99Xs2SkmSJEldlp1DZMlcdS3AJpnnQ4CXsweklN5OKc0pfX8LUIiI9fouREmSJEndUWhYkiFqs2SuqoeBzSNiaEQ0A0cCN2YPiIj3RUSUvh9dim9WH8YoSZIkqRvKF2atnwFRn3WZSym1RcTJwJ+ARuDKlNJTEXFiaf9lwKHAFyKiDZgPHJlSqiyrkyRJkjTA1GuXuT4bEEFHGdwtFdsuy3x/MXBxX8YkSZIkacXZZU6SJElSbpWXzC2iXgq9HBBJkiRJWmGNDQ00FNsBkBIsWuyASJIkSVKOVGaJ6oEDIkmSJEk9oh7nETkgkiRJktQj6rHTnAMiSZIkST2iHtcickAkSZIkqUc0NSzJEFkyJ0mSJClXmrMZIkvmJEmSJOVJITOHyJI5SZIkSblSniFyQCRJkiQpRwoNdpmTJEmSlFNNdpmTJEmSlFfZdYjazBBJkiRJypNCw5LhxULnEEmSJEnKk/KFWc0QSZIkScqRbMmcXeYkSZIk5UpTWcmcGSJJkiRJOVKwy5wkSZKkvCrvMueASJIkSVKOFCyZkyRJkpRXBZsqSJIkScqr7BwiM0SSJEmSciVbMmdTBUmSJEm5UrCpgiRJkqS8arZkTpIkSVJeFRoyTRUsmZMkSZKUJ03ZhVnNEEmSJEnKk+zCrAudQyRJkiQpT8q7zJkhkiRJkpQj2XWI7DInSZIkKVfKS+bMEEmSJEnKkWyGyC5zkiRJknKlKdt225I5SZIkSXniwqySJEmScquQmUPUZsmcJEmSpDwpa7ttyZwkSZKkPLHLnCRJkqTcarLLnCRJkqS8Ki+ZM0MkSZIkKUeyJXPOIZIkSZKUK+ULs5ohkiRJkpQjdpmTJEmSlFt2mZMkSZKUW3aZkyRJkpRbhYYlGaI2S+YkSZIk5UlzJkNkyZwkSZKkXClk225bMidJkiQpT+wyJ0mSJCm3ChUlcymlfoymNg6IJEmSJPWIxoYGGiI6ni9a7IBIkiRJUo4UylpvD/zGCn06IIqI8RHxbERMi4gzl3HcThGxKCIO7cv4JEmSJK2Y8sVZB/48oj4bEEVEI3AJsC+wDXBURGzTyXHfA/7UV7FJkiRJ6hnljRXMEGWNBqallJ5PKS0EJgIHVjnuS8DvgNf6MDZJkiRJPaC8ZM4MUdbGwIzM85bStg4RsTFwEHDZsk4UESdExOSImDxz5sweD1SSJElS9xTKSubMEGVFlW2VbScuBM5IKS3znUspXZ5SGpVSGrX++uv3VHySJEmSVlC2ZK6tDuYQNfXhtVqATTLPhwAvVxwzCpgYxVZ96wEfj4i2lNL1fRKhJEmSpBWSzRDVQ8lcXw6IHgY2j4ihwEvAkcCnsgeklIa2fx8RVwF/dDAkSZIk1Y/misVZB7o+GxCllNoi4mSK3eMagStTSk9FxIml/cucNyRJkiRp4Cs0ZDJElsyVSyndAtxSsa3qQCildGxfxCRJkiSp5zTZZU6SJElSXmUXZnUdIkmSJEm5Ur4w68DPEPVpyVxfaW1tpaWlhQULFvR3KJJWwCqrrMKQIUMoFAr9HYokSapR+cKsAz9DtFIOiFpaWlhzzTXZdNNNKbXwriqlxL3TZ3DR/ZN59JVXmdfaxmqFJrYfvCFf2W0ndv/AkGW+XlLvSSkxa9YsWlpaGDp06PJfIEmSBoTmsoVZzRD1iwULFix3MDRp6gucdOMkZi9YwNyFrR0rxL614F3+/c7z3DN9BuusugqXHDCOfTbzw5jU1yKCddddl5kzZ/Z3KJIkqQuaykrmBn6GaKWdQ7SswdDVf3+CIyfeQMvb7zAnMxhql4A5C1uZ8dY7HDHhBq7++xO9Gquk6szQSpJUfwp2mRvYJk19ga/cfDvz29pqOn5+Wxtfufl2Jk19oZcjkyRJkuqfXeYGsJQSJ904qebBULv5bW2cfNMkUqrMJXWusbGRkSNHMmLECHbYYQfuv//+roYLwIUXXsi8efNq2rfGGmt06xrLctVVV3HyySd36TWbbropr7/++lLbzz77bH7wgx+scEy33XYbW265JZttthnf/e53qx6TUuLLX/4ym222GcOHD+fvf/97TeeeOXMmhUKBn/70pysc50D35ptvctBBBzF8+HBGjx7Nk08+2bFv0003ZdiwYYwcOZJRo0Z1eo677rqLkSNHsu2227LXXnt1bP/sZz/LBhtswHbbbVd2/BlnnMHw4cM55phjOrb96le/4qKLLurBn0ySJPWneusyl6sB0b3TZzC7m53n3pi/gL++2FLz8auuuipTpkzhscce47zzzuOss87q1nW7MiCqRVsXB4MDzaJFizjppJO49dZbefrpp5kwYQJPP/30UsfdeuutTJ06lalTp3L55ZfzhS98oabz/+Y3v2GXXXZhwoQJPR16mYHwe/jOd77DyJEjefzxx7n66qs55ZRTyvbfeeedTJkyhcmTJ1d9/ezZs/niF7/IjTfeyFNPPcVvfvObjn3HHnsst912W9nxb731Fvfffz+PP/44ixYt4oknnmD+/PlcddVVfPGLX+z5H1CSJPWLQilDFAGLFteeUOgvuRoQXXT/ZOYubO3Wa+ctbOXC+6p/MFyet99+m3XWWafj+fnnn89OO+3E8OHD+da3vgXA3Llz2W+//RgxYgTbbbcd1113HT/+8Y95+eWX2Xvvvdl7773LztnZvq9//euMGDGCXXbZhVdffRUofjg97bTT2HvvvTnjjDP45z//yfjx49lxxx3ZY489+Mc//gEUBwPbbbcdI0aMYM899+w458svv8z48ePZfPPNOf300zu2T5gwgWHDhrHddttxxhlnVP3Zzz33XLbccks++tGP8uyzz3br/ct66KGH2GyzzfjgBz9Ic3MzRx55JDfccMNSx91www0cc8wxRAS77LILs2fP5pVXXlnu+SdMmMAFF1xAS0sLL730Usf2q6++muHDhzNixAg+85nPAPDqq69y0EEHMWLECEaMGMH999/P9OnTy7IiP/jBDzj77LMBGDNmDP/5n//JXnvtxUUXXcRNN93EzjvvzPbbb89HP/rRjt/XnDlzOO644xg2bBjDhw/nd7/7HT//+c859dRTO877s5/9jNNOO61b72G7p59+mrFjxwKw1VZbMX369I4YanHttddy8MEH8/73vx+ADTbYoGPfnnvuyXvf+96y4xsaGli4cCEpJebPn0+hUOD888/ny1/+sm21JUlaifzw4x9hzrdOY97ZX+X4nUb0dzjLtVJ2mevMo6+8ulQDhVolYMortX9YnD9/PiNHjmTBggW88sor3HHHHQBMmjSJqVOn8tBDD5FS4oADDuCee+5h5syZbLTRRtx8881A8a/pgwYN4oc//CF33nkn6623Xtn5v/zlLy+1b+7cueyyyy6ce+65nH766fzsZz/jG9/4BgDPPfccf/nLX2hsbGTs2LFcdtllbL755jz44IN88Ytf5I477uCcc87hT3/6ExtvvDGzZ8/uuNaUKVN49NFHec973sOWW27Jl770JRobGznjjDN45JFHWGeddRg3bhzXX389n/zkJzte98gjjzBx4kQeffRR2tra2GGHHdhxxx2Xeq+uueYazj///KW2b7bZZvz2t78t2/bSSy+xySabdDwfMmQIDz744FKvrXbcSy+9xODBg6v+vgBmzJjBv//9b0aPHs3hhx/Oddddx2mnncZTTz3Fueeey3333cd6663HG2+80fE72GuvvfjDH/7AokWLmDNnDm+++Wan54diVuXuu+8GiiVrDzzwABHBFVdcwfe//30uuOACvv3tbzNo0CCeeOKJjuOam5sZPnw43//+9ykUCvziF7+oWtZ3xBFHVB14nnbaaWVlagAjRozg97//PbvvvjsPPfQQL774Ii0tLWy44YZEBOPGjSMi+PznP88JJ5yw1Dmfe+45WltbGTNmDO+88w6nnHLKUtfIWnPNNTnkkEPYfvvtGTt2LIMGDeLhhx/mm9/85jLfM0mSVF8KmTlE9SBXA6J5rStWpjSvtfbsUnvJHMDf/vY3jjnmGJ588kkmTZrEpEmT2H777YFiNmDq1KnssccefPWrX+WMM85g//33Z4899uhyfM3Nzey///4A7Ljjjvz5z3/u2HfYYYfR2NjInDlzuP/++znssMM69r377rsA7Lbbbhx77LEcfvjhHHzwwR372z+8AmyzzTa8+OKLzJo1izFjxrD++usDcPTRR3PPPfeUDYjuvfdeDjroIFZbbTUADjjggKpxH3300Rx99NE1/YzV5nFV60RW63FZEydO5PDDDwfgyCOP5HOf+xynnXYad9xxB4ceemjHwLM983HHHXdw9dVXA8U5Y4MGDVrugOiII47o+L6lpYUjjjiCV155hYULF3astfOXv/yFiRMndhzXnl38yEc+wh//+Ee23nprWltbGTZs2FLnv+6665Z5/awzzzyTU045hZEjRzJs2DC23357mpqK/0m477772GijjXjttdfYZ5992GqrrcqyhlAs+3vkkUe4/fbbmT9/Prvuuiu77LILW2yxRafXPP300zuyjMcffzznnHMOV1xxBZMmTWL48OEdA3hJkqS+kqsB0WqFJt5a8O4KvL57ZT277rorr7/+OjNnziSlxFlnncXnP//5pY575JFHuOWWWzjrrLMYN25cl/9yXigUOj70NzY2ls1TWX311QFYvHgxa6+9dsdgLeuyyy7jwQcf5Oabb2bkyJEdx7znPe/pOKb9vLU2mKilbXJXMkRDhgxhxowZHc9bWlrYaKONlnptrcdlTZgwgVdffZVrrrkGKJYKTp06lZRSze2fm5qaWJxpL7mgYs5a++8B4Etf+hKnnXYaBxxwAHfddVdHaV1n1zv++OP5zne+w1ZbbcVxxx1X9fpdyRCttdZa/OIXv+i45tChQzsGZe3v1QYbbMBBBx3EQw89tNSAaMiQIay33nqsvvrqrL766uy555489thjyxwQtXv00UcB2GKLLTjllFO45557OPLII5k6dSqbb775cl8vSZLUU3I1h2j7wRvS3VVNAhg5eMNuvfYf//gHixYtYt111+VjH/sYV155JXPmzAGKpV2vvfYaL7/8Mqutthqf/vSn+epXv9rRFW3NNdfknXfeqXreZe3rzFprrcXQoUM7JsCnlHjssccA+Oc//8nOO+/MOeecw3rrrVc2oKi08847c/fdd/P666+zaNEiJkyYUNZlDIrzSP7whz8wf/583nnnHW666aaq5zr66KOZMmXKUo/KwRDATjvtxNSpU3nhhRdYuHAhEydOrJp5OuCAA7j66qtJKfHAAw8waNCgjnK5sWPHls0PAnj22WeZO3cuL730EtOnT2f69OmcddZZTJw4kbFjx/LrX/+aWbNmAXSUzI0dO5af/OQnQLHZw9tvv82GG27Ia6+9xqxZs3j33Xf54x//2Ol7+NZbb7HxxhsD8Mtf/rJj+7hx47j44os7nrdnnXbeeWdmzJjBtddey1FHHVX1nNddd13V97JaKdvs2bNZuHAhAFdccQV77rkna621FnPnzu24r+bOncukSZOW6hYHcOCBB3LvvffS1tbGvHnzePDBB9l66607/Xmz/uu//otzzjmH1tZWFpXacTY0NHS5SYgkSdKKytWA6JQPj2L15u5leVZrLvCV3TpvP1ypfQ7RyJEjOeKII/jlL39JY2Mj48aN41Of+hS77rorw4YN49BDD+Wdd97hiSeeYPTo0YwcOZJzzz23o3TohBNOYN99912qqcLy9i3LNddcw89//nNGjBjBtttu29GU4Gtf+1pHk4Q999yTESM6nwQ3ePBgzjvvPPbee++O1uIHHnhg2TE77LADRxxxBCNHjuSQQw7pVhlgpaamJi6++GI+9rGPsfXWW3P44Yez7bbbAsUM12WXXQbAxz/+cT74wQ+y2Wab8R//8R9ceumlQDFDNm3atKUm/E+YMIGDDjqobNshhxzChAkT2Hbbbfn617/OXnvtxYgRIzqaGVx00UXceeedDBs2jB133JGnnnqKQqHAN7/5TXbeeWf2339/ttpqq05/lrPPPpvDDjuMPfbYo2yO2De+8Q3efPPNjgYXd955Z8e+ww8/nN12262sSUd3PfPMM2y77bZstdVW3HrrrR2tr1999VV23313RowYwejRo9lvv/0YP348UP4eb7311owfP76jbffxxx/fMXA66qij2HXXXXn22WcZMmQIP//5zzuue/3117PTTjux0UYbsfbaa3f8W4iIZd5zkiRJvSG6srbOQDRq1KhU2Rb4mWeeqfqX6pQSW/zwclre7lpWBWCTQWvy7Kkn1Fw6pYHpySef5Morr+SHP/xhf4fSLfvvvz+nnnpqR3e4POjs37MkSVIXdPohPlcZoojg0gPHsWpT16ZOrdrUxCUHjHMwtBLYbrvt6nIwNHv2bLbYYgtWXXXVXA2GJEmSeluumioA7LPZUC7cbyxfufl25tewOOaqTU1ctP9H2WezoX0QnVTd2muvzXPPPdffYUiSJK10VtoB0bI6gx2zwzDet+YanHzTJN6Yv4B5C1vL1icKinOG3rvqKlxywDgHQ1I/qfeSXkmSNPCtlAOiVVZZhVmzZrHuuut2Oigat/lQnj31BP76YgsX3vcwU155jXmtraxWKDBy8IacutsodvvAEMvkpH6SUmLWrFmsssoq/R2KJElaia2UTRVaW1tpaWlZag0YSfVllVVWYciQIRS6uQaYJElSSadZjpUyQ1QoFDoWmJQkSZKkzuSqy5wkSZIkZTkgkiRJkpRbDogkSZIk5VbdN1WIiJnAi/0dR8Z6wOv9HYTqjveNusP7Rt3hfaPu8L5Rdwyk++b1lNL4ajvqfkA00ETE5JTSqP6OQ/XF+0bd4X2j7vC+UXd436g76uW+sWROkiRJUm45IJIkSZKUWw6Iet7l/R2A6pL3jbrD+0bd4X2j7vC+UXfUxX3jHCJJkiRJuWWGSJIkSVJuOSCSJEmSlFsOiLohIsZHxLMRMS0izqyyPyLix6X9j0fEDv0RpwaWGu6bo0v3y+MRcX9EjOiPODWwLO++yRy3U0QsiohD+zI+DUy13DcRMSYipkTEUxFxd1/HqIGnhv9PDYqImyLisdJ9c1x/xKmBJSKujIjXIuLJTvYP+M/FDoi6KCIagUuAfYFtgKMiYpuKw/YFNi89TgB+0qdBasCp8b55AdgrpTQc+DZ1MhFRvafG+6b9uO8Bf+rbCDUQ1XLfRMTawKXAASmlbYHD+jpODSw1/vfmJODplNIIYAxwQUQ092mgGoiuAqoueFoy4D8XOyDqutHAtJTS8ymlhcBE4MCKYw4Erk5FDwBrR8Tgvg5UA8py75uU0v0ppTdLTx8AhvRxjBp4avnvDcCXgN8Br/VlcBqwarlvPgX8PqX0L4CUkveOarlvErBmRASwBvAG0Na3YWqgSSndQ/Fe6MyA/1zsgKjrNgZmZJ63lLZ19RjlS1fvic8Bt/ZqRKoHy71vImJj4CDgsj6MSwNbLf+92QJYJyLuiohHIuKYPotOA1Ut983FwNbAy8ATwCkppcV9E57q2ID/XNzU3wHUoaiyrbJ3eS3HKF9qviciYm+KA6LdezUi1YNa7psLgTNSSouKf7SVarpvmoAdgbHAqsDfIuKBlNJzvR2cBqxa7puPAVOAjwAfAv4cEfemlN7u5dhU3wb852IHRF3XAmySeT6E4l9KunqM8qWmeyIihgNXAPumlGb1UWwauGq5b0YBE0uDofWAj0dEW0rp+j6JUANRrf+fej2lNBeYGxH3ACMAB0T5Vct9cxzw3VRcxHJaRLwAbAU81Dchqk4N+M/Flsx13cPA5hExtDSR8EjgxopjbgSOKXXV2AV4K6X0Sl8HqgFlufdNRLwf+D3wGf9Kq5Ll3jcppaEppU1TSpsCvwW+6GAo92r5/9QNwB4R0RQRqwE7A8/0cZwaWGq5b/5FMatIRGwIbAk836dRqh4N+M/FZoi6KKXUFhEnU+zm1AhcmVJ6KiJOLO2/DLgF+DgwDZhH8S8qyrEa75tvAusCl5b+2t+WUhrVXzGr/9V430hlarlvUkrPRMRtwOPAYuCKlFLVlrnKhxr/e/Nt4KqIeIJiGdQZKaXX+y1oDQgRMYFi18H1IqIF+BZQgPr5XBzFrKckSZIk5Y8lc5IkSZJyywGRJEmSpNxyQCRJkiQptxwQSZIkScotB0SSJEmScssBkSSpz0TEphGRIqLPW8pHxF0RcfEKnmNMKf71lnHMoRFhC1dJqhMOiCRJPaI0UFjW46r+jlGSpEouzCpJ6imDM9/vD/ysYtt8YJ3unDgiCiml1hWITZKkqswQSZJ6RErp3+0PYHbltpTSW5nDPxARf46IeRHxdETs074jU5b28Yh4KCIWAh+LotMj4p8RMT8inoiIT2djiIhvRsSLEfFuRPw7Iq6uCLMhIr4TEa9HxGsR8YOIaMi8fp2I+GVEvFm6xl8iYttl/dwRcUzpmvMi4o/Aht17ByVJ/cEBkSSpP5wL/BgYATwMTIyINSqO+R7wDWAr4EHgv4HPAScB2wDnAT+NiP0AIuIQ4KvAF4HNKWapHqo459FAG/Bh4GTgK8ARmf1XATsDBwKjgXnAbRGxarUfIiJ2Lr3mcmAkcBNwTm1vgSRpIIiUnPcpSepZEXEo8JuUUlRs3xR4ATgxpfTT0raNgRZgj5TSXyNiDHAncGhK6XelY1YHXgfGpZTuzZzvQmCLlNLHI+I04PPAdtXK6yLiLuA9KaVdM9v+DLyYUjo+IjYHngP2SindU9o/CPgX8H9TSldkYls/pfR6RFxb+j6b4boC+Fzlzy5JGpjMEEmS+sPjme9fLn3doOKYyZnvtwFWoZitmdP+AL4AfKh0zG9Kx7wQET+PiMMi4j3LuG77tduvuzWwGPhb+85Smd8TpetXs3X2+JLK55KkAcymCpKk/tCRwUkppYiApf9INzfzffu+T1DM2Cx1rpTSjIjYEhgLfBS4APhWROycUpqbPTYjZc69rIxOZ+UUZoEkqc6ZIZIk1YOngXeBD6SUplU8Xmw/KKW0IKV0c0rpVGAnYFtgty5cowHIltStBQwr7evsNbtUbKt8LkkawMwQSZIGvJTSOxHxA+AHUUwn3QOsQXHwsTildHlEHEvx/2sPAnMoNktoBabWeI2pEXEDxUYNJ1DslHcu8DZwbScv+zFwf0ScBfwWGAMc1J2fUZLUP8wQSZLqxX8BZ1PsJPcU8GfgEIpNGqA4gPkccC/wZGnfwSmlFypPtAzHUexMd2Pp62rA+JTS/GoHp5QeKF3zCxTnJx1cilGSVCfsMidJkiQpt8wQSZIkScotB0SSJEmScssBkSRJkqTcckAkSZIkKbccEEmSJEnKLQdEkiRJknLLAZEkSZKk3HJAJEmSJCm3/j+V8qyQerWsowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "evals_df = pd.DataFrame(evals)\n",
    "\n",
    "best_thresh = evals_df.sort_values(by='Accuracy', ascending=False).iloc[0]\n",
    "\n",
    "plt.plot(evals_df['Threshold'], evals_df['Accuracy'], lw=3, c='#087E8B')\n",
    "plt.scatter(best_thresh['Threshold'], best_thresh['Accuracy'], label=f\"Best threshold = {best_thresh['Threshold']}, Accuracy = {(best_thresh['Accuracy'] * 100):.2f}%\", s=250, c='#087E8B')\n",
    "plt.title('Threshold Optimization', size=20)\n",
    "plt.xlabel('Threshold', size=14)\n",
    "plt.ylabel('Accuracy', size=14)\n",
    "plt.legend()\n",
    "plt.savefig('images/003_LogisticRegression_threshold_optimization.png', dpi=300, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- You can now retrain with the best threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-6-822322b64438>:18: RuntimeWarning: overflow encountered in exp\n",
      "  return 1 / (1 + np.exp(-x))\n"
     ]
    }
   ],
   "source": [
    "model = LogisticRegression()\n",
    "model.fit(X_train, y_train)\n",
    "preds = model.predict(X_test, threshold=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's evaluate again for accuracy and confusion matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.956140350877193"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(y_test, preds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[43  0]\n",
      " [ 5 66]]\n"
     ]
    }
   ],
   "source": [
    "print(confusion_matrix(y_test, preds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "## Comparison with Scikit-Learn\n",
    "- We want to know if our model is good, so let's compare it with `LogisticRegression` model from Scikit-Learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dradecic/opt/anaconda3/lib/python3.8/site-packages/sklearn/linear_model/_logistic.py:762: ConvergenceWarning: lbfgs failed to converge (status=1):\n",
      "STOP: TOTAL NO. of ITERATIONS REACHED LIMIT.\n",
      "\n",
      "Increase the number of iterations (max_iter) or scale the data as shown in:\n",
      "    https://scikit-learn.org/stable/modules/preprocessing.html\n",
      "Please also refer to the documentation for alternative solver options:\n",
      "    https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression\n",
      "  n_iter_i = _check_optimize_result(\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "lr_model = LogisticRegression()\n",
    "lr_model.fit(X_train, y_train)\n",
    "lr_preds = lr_model.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.956140350877193"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy_score(y_test, lr_preds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[39  4]\n",
      " [ 1 70]]\n"
     ]
    }
   ],
   "source": [
    "print(confusion_matrix(y_test, lr_preds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- They are more or less the same, accuracy-wise"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}