{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "skip"
},
"tags": [
"hide-cell"
]
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"import matplotlib as mp\n",
"import sklearn\n",
"from IPython.display import Image, HTML\n",
"import statsmodels.api as sm\n",
"from sklearn import model_selection\n",
"from sklearn import metrics\n",
"\n",
"import laUtilities as ut\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Logistic Regression"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"So far we have seen linear regression: a continuous valued observation is estimated as linear (or affine) function of the independent variables."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Today we will look at the following situation.\n",
"\n",
"Imagine that you are observing a binary variable -- a 0/1 value.\n",
"\n",
"That is, these could be pass/fail, admit/reject, Democrat/Republican, etc."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"You believe that there is some __probability__ of observing a 1, and that probability is a function of certain independent variables."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"So the key properties of a problem that make it appropriate for logistic regression are:\n",
" \n",
"* What you can observe is a __categorical__ variable\n",
"* What you want to estimate is a __probability__ of seeing a particular value of the categorical variable."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## What is the probability I will be admitted to Grad School?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"```{note}\n",
"The following example is based on http://www.ats.ucla.edu/stat/r/dae/logit.htm.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"A researcher is interested in how variables, such as _GRE_ (Graduate Record Exam scores), _GPA_ (grade point average) and prestige of the undergraduate institution affect admission into graduate school. \n",
"\n",
"The response variable, admit/don't admit, is a binary variable."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"There are three predictor variables: __gre,__ __gpa__ and __rank.__ \n",
"\n",
"* We will treat the variables _gre_ and _gpa_ as continuous. \n",
"* The variable _rank_ takes on the values 1 through 4. \n",
" * Institutions with a rank of 1 have the highest prestige, while those with a rank of 4 have the lowest. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"hide_input": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"df1 = df[df['rank']==1]\n",
"df2 = df[df['rank']==2]\n",
"df3 = df[df['rank']==3]\n",
"df4 = df[df['rank']==4]\n",
"#\n",
"fig = plt.figure(figsize = (10, 5))\n",
"ax1 = fig.add_subplot(221)\n",
"df1.plot.scatter('gre','admit', ax = ax1)\n",
"plt.title('Rank 1 Institutions')\n",
"ax2 = fig.add_subplot(222)\n",
"df2.plot.scatter('gre','admit', ax = ax2)\n",
"plt.title('Rank 2 Institutions')\n",
"ax3 = fig.add_subplot(223, sharex = ax1)\n",
"df3.plot.scatter('gre','admit', ax = ax3)\n",
"plt.title('Rank 3 Institutions')\n",
"ax4 = fig.add_subplot(224, sharex = ax2)\n",
"plt.title('Rank 4 Institutions')\n",
"df4.plot.scatter('gre','admit', ax = ax4);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Logistic Regression"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Logistic regression is concerned with estimating a __probability.__\n",
"\n",
"However, all that is available are categorical observations, which we will code as 0/1.\n",
"\n",
"That is, these could be pass/fail, admit/reject, Democrat/Republican, etc."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Now, a linear function like $\\alpha + \\beta x$ cannot be used to predict probability directly, because the linear function takes on all values (from -$\\infty$ to +$\\infty$), and probability only ranges over $(0, 1)$.\n",
"\n",
"However, there is a transformation of probability that works: it is called __log-odds__."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"For any probabilty $p$, the __odds__ is defined as $p/(1-p)$. Notice that odds vary from 0 to $\\infty$, and odds < 1 indicates that $p < 1/2$.\n",
"\n",
"Now, there is a good argument that to fit a linear function, instead of using odds, we should use log-odds. That is simply $\\log p/(1-p)$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"hide_input": false,
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEhCAYAAACUW2yNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt7UlEQVR4nO3deXhU5d3/8fcXAoQtbFnYd8K+iIi4gLKoqHWp1bZal9Zard0XH7v92tra1W7P01bb2qpYt1ZbtS61FUXEXUB2DGGHBMhCSCD7dv/+OCc6JgFmyMycyeTzuq65kjlz5pzvPZPMZ859n8Wcc4iIiITqFHQBIiKSeBQOIiLSgsJBRERaUDiIiEgLCgcREWlB4SAiIi0oHEREpAWFg4iItKBwaMfM7JNm5sxsbJzWd5uZuWbTLjWzr4X5/KZ6W7stik3VYdV11Da01uagxKsWM/udmT3dhucPMrNGMzvDvz/UX+YbZlbpv98jo1ZwG5nZV81svZnp8zCEXgyJxF+A05pNuxQIKxxCXOEvJ/T2dluLa4NLOXobWmtz0jKzMcBNwA/asJhLgSLgDf/+WOCjwCHglbbUFyN/BDKB64IuJJGkBF2AtB/OuTwgLwqLWuuc2xaF5cRcFNvcXnwFWOecW9WGZVwKPOWca/Tvr3DOZQGY2Q3AuW2qMMqcc1Vm9lfgFuC+oOtJFNpySHJmttjfnK8yszIze9LMxrcy35VmlmNm1Wa2wcwuNrPlZrY8ZJ4PdGuY2RK8b1tDQrqHdrWx3iWtLaN5LaH1mNk4M3vWzMrNbLeZfa95F4GZTTezJ8zsoP9abDGzbx2vDUfpSjvuaxpJbW0R5fe3G3A18HArz3/bzB41sx+a2XZ/OevNbGGz+dKAs4Enm6aFhETchVs38DdgkpmdHkSdiUjhkMTMbDHwLFAOfAy4GZgCvGpmQ0LmOwd4CMgBPgL8EvhfIPs4q7gd+DdeF0JT99CHwyits5mlhNw6R9KuVjwBLMP7xvokXpfIe10EZjYbr4tjDPBV4ELg18DQSNsQ7msaQW0j/RC5LbImx+T9nQP0pVnXj5mlAFPxXrcz8V7DK/F6Hh43s/SQ2S8EaoEXIm1PtEVY91rgMLA4zmUmLuecbu30BnwScMDYozy+CtgKpIRMGwXUAb8OmfY6sBGwkGkz/WUvD5l2m/cn84F1LAHyIqy3+e3VZsvb1cpzl4fWEloP8Klm0zcAz4fcXwHsBXocpa6jtqF5myN4TcOtbQRQD3wvjNfvRGsJ9/39BtAIdG223mn+vC8DnUOmn+VP/3DItL8Djx2jDTf4zxkZh/+PsOv2p78S+t509Ju2HJKUmfXE+wD4u3Ouvmm6c24n8BrePwj+t/ZZwD+d/x/iz/cOsDNG5X0YOCXk9uk2Lu/ZZvc3AsMBzKwHcAbwkHOusi0rCfc1Dbc2/7m7nXMpzrkfxqKWCN/fwcBh51xts+kn+z+/7ZxrCJme4/8c4K+rK9437ycjaUtrzGyRHX3PttDb8mMsJqy6QxThvQaCBqSTWT/AgP2tPHYA7xsrQDrQBShsZb6C2JTGRhfdAemSZvdrgFT/93543afRGFQO9zUNt7Z41BLJ+5vq19fcTGCfc+61ZtObPkibXtuFQHdaBuKJeB2YGMZ8xwr8cOtuUoVXv6BwSGaH8DadB7by2EDgoP97MV43RGYr82UBe2JS3dFVA11bmT6A92uOxCG8rpLWxgNOZFnhvKbxEIv39yBe6DQ3E8hvZfrH8D6cm8YoLgVeds6VHrv04/O38nKOO+OxhVt3k/54r5egAemk5ZyrAFYDV4QO+JrZCOB0vH5Y/M3tVcBHzMxC5jsZr//6eGqI7ret3UBW6GChefvet9gDJxz+h8yrwNVmdrQ6w2pDuK9pPMTo/c0BupjZ0JD5OgHTgRH+AG/T9MHA54DfO+cq/GVfRBS6lKIh3LqbPW0UsCV+VSY2hUNyWGxmlze7nQN8FxgHPGNmF5nZlcBSoAz4Vcjzvw9MBp4wswvM7FrgMbzuiePthrgZ6G9mN5vZKWY2tY1teQzvG/FDZnaemX0C+Bdt+0Z3C96Wxxtmdo2ZzTezT5vZ7/zHI2lDuK9pWMxshJnVm9n3In1uBLWE+/6u8H/ODpk2AeiJF6BLzGyhmV2H9607x182eHs6DcJ7r1pr5+VmdjnvjwOc709rbZwmGsKtu6m+vnh7b61APEGPiOt24jeOvvePw+vXB2+A8A28/tQyvH/e8a0s6yq8b001wCa8QeM1wBMh89xGy72VegKP8H43x64w6m1176qQ+S7FG7itAtbhHTS1nKPvrZTSbPqS5nUAJwFPA6X+cnOAbxyvDUdp83Ff03BrA0b6890Wxvt9QrWE+/76870F3Bdy/2q/vinAM3jdMQXA74DeIfP9HFh5jNqP9ne6/HjtPsH/jbDqDpn/E3hdmgNi+T/bnm7mvzAiH+B3LWwDfuycuz3oeiS6jvb+mtkngf8DBjnnKs3sN8Dlzrlhx1leDvCAc+7HMSw7bOHWHTL/c0Cxc+6a2FbWfigcBL8v/td4By4VA6OBW/EGLCc751rbI0baiUjeX3/8YgNwr3Pul2b2MnDIOXdp3Atvg0jqNrMZwJvAFNdOTusSD9pbSQAa8PZw+T1e33wFXt/sFQqGpBD2++ucazCz64GZ/iDzDLwjqtuNE6h7IN7BigqGENpyEBGRFrS3koiItJA03Urp6elu5MiRQZchItKurF69utg5l9F8etKEw8iRI1m1qi2noBcR6XjMbHdr09WtJCIiLSgcRESkBYWDiIi0oHAQEZEWFA4iItKCwkFERFpQOIiISAsKBxGRdmp7UTm/XppL4eHqqC9b4SAi0k69vv0gv31xK7UNx7smV+QUDiIi7dT6vaUM6NmVIX2jeaVej8JBRKSd2pBfxtShfQi5PHjUKBxERNqhytp6cguOMG1o35gsX+EgItIObd53mEYH04b0icnyFQ4iIu3QurwyAKYNVTiIiIhvQ14pA9NSyUxLjcnyFQ4iIu3Qen8wOlYUDiIi7czh6jp2FFUwXeEgIiJNNuZ74w1TY7SnEiR4OJjZODOrNrMHg65FRCRRrG8ajI7RnkqQ4OEA3AmsDLoIEZFEsiGvjGH9u9OvZ9eYrSNhw8HMPg6UAi8GXIqISEJZn18as4PfmiRkOJhZGvBD4OvHme9GM1tlZquKioriU5yISIBKKmrZW1IV0y4lSNBwAG4H7nHO7T3WTM65u51zs5xzszIyMuJUmohIcDbkNx381jem60mJ6dJPgJnNABYBJwVciohIwnltWzFdOhtThqTFdD0JFw7A2cBIYI9/psFeQGczm+ScmxlgXSIigWpsdDy9bh/zxmXQO7VLTNeViN1KdwNjgBn+7Y/As8B5wZUkIhK81XsOsb+smotnDI75uhJuy8E5VwlUNt03s3Kg2jmnEWcR6dCeWruP1C6dWDQxK+brSrhwaM45d1vQNYiIBK2+oZF/b9jPwolZ9OwW+4/uROxWEhGRZl7ffpCDFbVcNC32XUqgcBARaReeWreP3t1SOHt8fHbbVziIiCS4mvoG/rvxAOdOHkhql85xWafCQUQkwS3fUsSRmnoumj4obutUOIiIJLglr+0is3c3zhibHrd1KhxERBLY6t2HeGPHQW6cN5ouneP3ka1wEBFJYHe+tI1+Pbpw1anD47pehYOISILamF/GspxCrj9jFD26xvewNIWDiEiCumv5Nnp3S+Ha00fGfd0KBxGRBLSt8AjPbTzANaeNoE/32J5krzUKBxGRBPSr53PpltKJ688cFcj6FQ4iIglm+ZZCntt4gM+fPZb0Xt0CqUHhICKSQKrrGvjevzYxOqMnN541OrA6Ev6srCIiHcldL21jT0klD99wKt1S4nOqjNZoy0FEJEHsKCrnjy/v4JIZgzk9jkdDt0bhICKSAOobGrn1H+vp1qUT37lwYtDlqFtJRCQR/OaFXFbtPsT/fXwGmb1Tgy5HWw4iIkFbkVvEXcu38/FThnHJjCFBlwMoHEREAlV4uJqv/n0t4zJ78f2LJgddznvUrSQiEpDqugY+//A7VNTW87er5tC9a3B7JzWncBARCUBjo+N//rGelbsO8dsrT2JcVu+gS/oAdSuJiATgl89v4el1+7h18Xgunj446HJaUDiIiMTZw2/t4a7l27ly9jBuPmtM0OW0SuEgIhJH/1ydx3ee3MBZ2RncfskUzCzoklqlcBARiZMn1uRxyz/WcfqYAfzx6pNJieNlPyOVuJWJiCSRJ9fk8/VH13Ha6AH85dpTEmrPpNZobyURkRhb8tpOfvDMZuaMGsA91yV+MIDCQUQkZpxz3PHfLfxh+XbOm5zF/338JFK7JH4wgMJBRCQmauob+NbjG3j8nXyuOnU4t18yhc6dEnPwuTUKBxGRKCs8XM1ND65mzZ5SvnZONl9cMDZh90o6GoWDiEgUrd1byk0PrOJwVT13fWImF0wdFHRJJ0ThICISBc45HnhzNz965l0y07rxz5tPZ9LgtKDLOmEKBxGRNiqrquOb/1zPcxsPcPb4DH790Rn079k16LLaJOHCwcy6AXcBi4D+wDbg28655wItTESkFW/vLOFrj67lQFk1375gAjecOZpO7Wjg+WgSLhzwatoLnAXsAS4AHjWzqc65XUEWJiLSpLqugV8vzeXPr+xgWL8ePPrZ05g5vF/QZUVNwoWDc64CuC1k0jNmthM4GdgVRE0iIqHW7i3l1n+sI7egnE+cOpxvXzCRnt0S7uO0TRK+NWaWBWQDm1p57EbgRoDhw4fHuTIR6Wgqaur55fNbWPL6LrJ6p3Lfp05h/vjMoMuKiYQOBzPrAjwE3O+cy2n+uHPubuBugFmzZrk4lyciHYRzjqWbC/jB05vZV1bF1aeO4NbF4+md2iXo0mImYcPBzDoBDwC1wBcCLkdEOqhdxRXc9vQmlm8pIjurF4/ddBqzRvYPuqyYS8hwMO9QwnuALOAC51xdwCWJSAdzuLqOO1/axn2v7qJrSie++6FJXHvaCLok8Gm2oykhwwH4AzARWOScqwq6GBHpOOobGvnbyr38ZmkuBytquWzmEL65eAKZaalBlxZXCRcOZjYCuAmoAQ6EnI/kJufcQ4EVJiJJzTnHcxsP8Mv/bmFHcQWnjurPkgsnMXVon6BLC0TChYNzbjfQ/o8gEZF2wTnHK1uL+eXzW1ifV8a4zF7cfc3JnDMpq92dLC+aEi4cRETiwTnH69sP8uuluazefYghfbvzi8uncdnMoe3q1NqxonAQkQ7FOcfLuUXc+dI2Vu46xMC0VG6/dAofnTWUbint40I88aBwEJEOoaHRsXTzAe58aTsb8ssY1CeVH1w8mY+dMqzdXJ0tnhQOIpLUqusaeGJNPn9esYMdxRWMGNCDn39kKh8+aShdUzrGbqknQuEgIknpYHkND721h7++sZvi8hqmDEnj91edxPlTBmlMIQwKBxFJKjkHDnP/67t4/J18auobOXt8Bp+ZO5rTxwzo0HsfRUrhICLtXn1DIy+8W8CS13fx5o4SUrt04rKZQ/j0maMYm9k76PLapYjCwcwyAJxzRf79qcDHgE3OuUeiX56IyNEVHK7mkbf38Mjbeyg4XMOQvt351vkT+Ngpw+jbo31fiS1okW45PIp3Mrx7zSwdWAHsA75oZoOdc7+KdoEiIqEaGh0rthbxyFt7eDGnkIZGx1nZGfzo0hEsmJCp8YQoiTQcpgFv+r9fDmxzzp1iZpcAvwAUDiISE/mlVTy2ai+Prcojv7SKAT27csPcUVx5ynBGpvcMurykE2k4dAfK/d8XAU/5v78DDItWUSIi4O2GunRzAY+tzuOVrUUAnDk2nW9dMIFzJw3UrqgxFGk4bAUuM7N/AufibS2Ad2rt0ijWJSIdlHOONXtLefydPJ5au4/D1fUM6dudLy4YxxUnD2VY/x5Bl9ghRBoOPwAewes+etE595Y//TxgTTQLE5GOZW9JJU+uyefxNfnsLK6gW0onzp8ykCtmDeO00QPopLGEuIooHJxzj5vZcGAwsC7koReAf0azMBFJfiUVtTy7YT9Prsln9e5DAMwZ3Z+bzx7D+VMGJvVlOBNdxMc5OOcKgIJm0946yuwiIh9QXlPP0s0HeGrtPl7ZWkx9oyM7qxe3Lh7PxdMHM7Sfuo0SwXHDwczuDXdhzrnr21aOiCSjqtoGXtpSyNPr9rEsp5Ca+kaG9O3ODXNHc/H0wUwc1FtHLyeYcLYcMprdnwc0Ahv8+1OATnjHPIiIAF4gLN9SyLMb9vPiu4VU1TWQ3qsbHz9lGB+aPpiTh/fTOEICO244OOcuavrdzL4FVAGfcs5V+NN6AvfwfliISAdVUVPP8i1F/HvDfpbleIEwoGdXLps5hAunDuLU0QN0kFo7EemYw5eAhU3BAOCcqzCz24EXgR9HszgRSXxlVXW8+G4Bz208wIrcImrqG0nv1Y2PnDyE86cM4tRR/UnprOMR2ptIw6EX3p5Km5tNHwRoFEmkgyg8XM1/Nxfw/KYDvLH9IPWNjoFpqVw5eziLpwzklJH9tYXQzkUaDv8E7jOz/+H902jMAX4OPB7NwkQksWwrLGfp5gL+u+kAa/eWAjA6vSc3zB3NeZOzmD60r8YQkkik4XAz3gFwS4CmHZDr8cYcboleWSIStIZGx5o9h1i6uYClmwvYUez1Jk8b2odbzs3m3MkDGZfZS3sZJalID4KrAj7nbzmMAQzv5HsVx36miLQH5TX1vLq1iKWbC3lpSyElFbV06WzMGT2AT54xkkUTsxjct3vQZUocnNDFfvwwWB/lWkQkAHtLKlmWU8gL7xbw1o4Sahsa6dO9C/PHZ7BoUhbzsjNI05HKHU44B8G9BLhwFuacW9DmikQkpuobGnlnTynLcgpZllNAboF3ouXR6T257vQRLJyYxawR/bSHUQcXzpbDxpDfOwOfAA4ATafMmI23t9KD0S1NRKKlpKKWl3MLWZZTxIrcIsqq6kjpZMwe1Z+PzhrGggmZjM7oFXSZkkDCOQjui02/m9lvgPuBLzvnXMj0/8UbfxCRBNDY6Ni07zAvbfHGDtbuLcU5SO/VlUUTs1g4MZMzx6Wru0iOKtIxh2uB00KDwXcX3q6tX45KVSISsbLKOlZsLWL5liJezi2iuLwGM5g2pA9fXjiO+eMzmTqkj3Y3lbBEGg4GTAVym02fGp1yRCRcTVsHy7cU8nJuEe/sOUSjg749ujB3XAbzx2cwLzuD9F7dgi5V2qFIw+Fe4C9mNo4PHgR3K3BfNAsTkZZKKmp5ZWsRL28pYsXWIorLawHv2IPPzx/L2eMzmTGsr45OljaLNBxuBQrxuo9+4k/bD/wM7+A4EYmi+oZG1uWV8vKWIl7eWsz6PG/soF+PLszLzuCsbG0dSGxEehBcI3AHcIeZpfnTDseiMJGOan9ZFStyvXGDV7cWc7i6nk4GJw3vx1cXZTMvO4OpQ/po60BiKuKD4MxsKjAe79iHXHSqbpE2qa5r4O2dJbyc6+1murXQO+5gYFoqi6cM5KzsTM4cm06fHtqzSOIn7HAws5PxxhUm8/5uq87MNuJd3+GdaBVlZv3xztd0LlAMfMs593C0li8SJOcc2wrLvTDYWsxbOw5SU99I15ROzB7pHXcwLzuD7Cydt0iCE1Y4mNl44CVgK3AN3im7DZgEfB14ycxmO+e2RKmuO4FaIAuYATxrZuucc5uitHyRuCqtrOW1bQdZkesNJO8vqwZgTEZPrpw9nLPGZzBn1AC6d+0ccKUiHmt5yEIrM5k9AvQELml+jIN5X23+BVQ4565sc0HeleUOAVOcc7n+tAeAfOfcN4/2vFmzZrlVq1a1dfUiUfHeQHJuMStyi1ifV0qjg96pKZw5Np152RnMHZfO0H66DIoEy8xWO+dmNZ8ebrfSAuCiVg5+wznnzOwneAERDdlAQ1Mw+NYBZzWf0cxuBG4EGD58eJRWL3Ji9pVWvbdl0DSQbAbTh/blCwvGcVZ2OtOH9tU5i6RdCDcc+gL7jvF4HtCnzdV4egFlzaaVAb2bz+icuxu4G7wthyitXyQsVbUNvLXzICtyi1mxtYhtzQaS52VncObYdPr26BpwpSKRCzcc8oBp/s/WzADyo1EQUA6kNZuWBhyJ0vJFTohzjq2F5e8dgPbWzhJq/YHkU0f15+OneAPJugCOJINww+EJ4Bd+31RB6ANmNhDvILgnolRTLpBiZuOcc1v9adMBDUZL3JVV1vHa9uL3AqFpIHlsZi+umTOCedkZnDqqP6ldNJAsySXccLgduBDYZmYPAjl4xzlMBq7C22r4UTQKcs5VmNnjwA/N7Aa8rZJLgNOjsXyRY2lsdGzIL3vvmIM1e0tpaHTvDSR/eaF3RLKuhibJLqxwcM6VmdlpeKfM+BjQz3+oFHgA+I5zrjSKdX0O7zxOhcBB4GbtxiqxUlxewyv+2Uxf2VpMSUXte2cz/dzZYzh7fIYGkqXDCfsgOP/D/3Nm9nkgw59c1NoeTG3lnCsBLo32ckUAGhoda/ceYvkWLxA25Hv7Pwzo2ZWz/XMVzR2XzgCdr0g6sIhPn+GHQaGZDcU7EE57CUnCKy6v4eUtRSzPff9KaJ0MZg7vx9fPyeas8RlMGaxrHYg0iTgcQmzGGw/YEZ1SRKKnsdGxPr+Ml3IKWb6lkPX5Zf6V0LqxaGIW8ydkMHdshs5XJHIUbQkHfcWShFJWVccrW4tYllPIy1uKOOiPHcwY1pevLcpm/oRMJg1K09aBSBjaEg4igXLOsaO4gmXvFvJiTgErdx2iodHRt0cX5o3LYMGETOZlZ9C/pw5CE4lUW8LhJ0BJtAoRCUddQyMrd5bwwruFLMspYNfBSgAmDOzNjfNGs2BCJicN055FIm11wuHgnPtpNAsROZqyqjqWbynkhXe98YMj1fV0TenE6WMG8OkzRzF/QqZOYCcSZRGFg5nde5SHHFANbAP+7pw71nmYRI5rb0klL7xbwNLNBby9s4T6Rkd6r64snjyQRZOyOHNsOj27qVdUJFYi/e/KAOYCjcBGf9oUvMHp1cBleEc2z3XOrY1WkZL8nHO8u/8Iz28+wPObCti837v67LjMXnxm3mjOmZTFjKF9NZgsEieRhsNreCfG+7RzrhLAzHoAf8Y7rfYFwF+BXwELo1inJKGGRsc7ew7x340H+O/mA+wtqcIMZo3ox7cvmMA5kwYyKr1n0GWKdEiRhsOXgQVNwQDgnKs0sx8DLzrn7jCznwMvRLNISR51DY28taOE5zbu5/nNBRQdqaFr506cMXYAnz97LAsnZpHRW0cmiwQt0nDoBQwC3m02faD/GMDhE1iuJLHa+kZe217Mcxu8QCitrKNH187MH5/JeVMGMn98Br1TdTCaSCKJ9EP8CeAeM7sVWIk3ED0buAN43J9nNt5pt6UDq2to5LVtxTy7fj//3XSAw9X19O6WwqJJWSyeMpCzsjN0mmuRBBZpOHwW+DXwYMhz6/HOoHqLf/9d4DNRqU7alfqGRt7cUcIz6/fxn00HKK2so3dqCudMyuLCqYM4c1w63VIUCCLtQUTh4I81fNbMvg6MwdtLaZtzriJknrVRrVASmnOONXtLeWrtPp5Zv4/i8lp6du3MOZOy+NC0wczNViCItEcnOjbQgLc7q/N/lw5mW2E5/1qbz7/W7mNPSSVdUzqxcEImF00fzIIJmeoyEmnnIj0ILgX4KfAFoCvelkONmf0O74I/ddEvURLFwfIanl63j8fX5LM+r4xOBmeMTeeLC8Zy3pSBpGlQWSRpRLrlcAdwJd7Yw6v+tLl4gdGJ98cdJEnU1jeyLKeQf6zOY/mWQuobHZMGpfH/LpzIRdMHk5WWGnSJIhIDkYbDVcD1zrl/h0zbbmZFwF9QOCSNnAOHeXRlHk+uzaekopaM3t24/sxRXDZzCBMGpgVdnojEWKTh0AfY3sr07UDfNlcjgSqvqeeptfv4+6q9rNtbSpfOxjmTsrji5GHMHZeuM52KdCCRhsM64EvA55tN/7L/mLRDG/LKePjtPTy1Np+K2gbGZ/Xmux+axIdPGqJrIYh0UJGGw63Av83sHOANvL2VTgMGA+dHuTaJoeq6Bp5et48H39zNurwyUrt04kPTBnPl7OHMHN4XM53gTqQji/Q4hxVmlo235TABb2+lx4B/A1/h/UFqSVB7Syp58M3d/G3lXsqq6hib2YvbLprEh2cOpU937W0kIp6Ij3Pwr9XwndBpZjYd+Ei0ipLocs7x9s4S7n1tJ0s3F2BmnDspi2tPG8mc0f21lSAiLegEeUmsrqGRf2/Yz19e2cmG/DL69ujCZ88aw9VzRjC4b/egyxORBKZwSEIVNfX8beVe7nllB/vKqhmd0ZMff3gKl500lO5ddeSyiByfwiGJlFbWsuT1Xdz32i7KquqYPbI/t186hfnjM3UFNRGJSFjhYGZPHWcWHRUVoOLyGv7yyk4efHM35TX1nDMpi8+eNYaTR/QLujQRaafC3XI4GMbjO9tYi0SopKKWu1fs4P7Xd1FT38CF0wbz+fljdASziLRZWOHgnPtUrAuR8B2uruPPK3Zw76s7qaxr4OLpg/nSwnGMyeh1/CeLiIRBYw7tSHVdAw++uZs7X9rGoco6Lpw6iK8sGse4rN5BlyYiSUbh0A4453h6/X5+/lwO+aVVzB2XzjcWT2DKkD5BlyYiSUrhkODW7DnED5/ZzJo9pUwalMYdl0/jjLHpQZclIklO4ZCgDpbX8PP/5PDoqjwye3fjjsun8ZGZQ+msXVJFJA4SJhzMrBtwF7AI6A9sA77tnHsu0MLirLHR8dDbe/jFf3KorG3gpnmj+eLCcfTqljBvlYh0AIn0iZMC7AXOAvYAFwCPmtlU59yuIAuLl60FR/jm4xtYvfsQp48ZwA8vmczYTA02i0j8JUw4OOcqgNtCJj1jZjuBk4FdQdQUL3UNjdz50jbufGkbPbul8KsrpnPZzCE6IZ6IBCZhwqE5M8sCsoFNx5jnRuBGgOHDh8epsujaWnCErz66lo35h7l4+mC+d9Ek0nt1C7osEengEjIczKwL8BBwv3Mu52jzOefuBu4GmDVrlotTeVHR2Oi47/Vd/Pw/OfTqlsIfr57J4imDgi5LRASIYziY2XK88YTWvOacO9OfrxPwAFALfCE+1cVXSUUttzy2jmU5hSyamMlPL5tGRm9tLYhI4ohbODjnzj7ePOZ1st8DZAEXOOfqYl1XvL29s4QvPbKGkopafnjJZK6ZM0JjCyKScBKtW+kPwERgkXOuKuhiosk5x/2v7+L2Z99lWL/uPP6503WEs4gkrIQJBzMbAdwE1AAHQr5N3+SceyiwwqKgpr6B7z65kUdX5bFoYha/+dh0eqfqes0ikrgSJhycc7uBpOtfOVhew2f+uop39pTypQVj+cqibF14R0QSXsKEQzLaW1LJtfe+zb7SKu68aiYXTtPeSCLSPigcYmRjfhmfWrKS2vpGHv7MqZw8on/QJYmIhE3hEAOrdpXwyftWkpaawiM3n6ZTYIhIu6NwiLLVu0u47t63yUpL5aHPnMqgPt2DLklEJGIKhyh6Z88hrrt3JZlpqTxy4xyy0lKDLklE5IR0CrqAZLExv4zr7nmb9F5deeQzCgYRad8UDlGQX1rF9UtW0js1hUdunMPAPgoGEWnf1K3URoer67j+vpVU1Tbwj5tP1xiDiCQFhUMb1DU08rkH32F7UTn3Xz+b8QO1V5KIJAeFQxv87LkcXt1WzC8un8YZY9ODLkdEJGo05nCClm4u4J5Xd/LJ00dyxaxhQZcjIhJVCocTkF9axS2PrWPKkDS+dcGEoMsREYk6hUOE6hoa+dIja2hodPz+ypl0S+kcdEkiIlGnMYcI/enl7azefYjfXnkSI9N7Bl2OiEhMaMshAnsOVvK7Zdu4cOogLp4+OOhyRERiRuEQJucc33tqIymdjO9+aFLQ5YiIxJTCIUz/2XiA5VuK+Nq543UEtIgkPYVDGMpr6vnB05uZNCiN604bEXQ5IiIxpwHpMPzp5e0cOFzNXVfPJKWz8lREkp8+6Y6jrKqOJa/t4oKpA5k5vF/Q5YiIxIXC4TgeeGMXR2rq+dzZY4MuRUQkbhQOx1BZW889r+5k/vgMpgzpE3Q5IiJxo3A4hoff2sOhyjq+sGBc0KWIiMSVwuEoqusauHvFDk4bPYCTR2isQUQ6FoXDUTz+Tj6FR2r4wgKNNYhIx6NwOIrHVu9lwsDenD5mQNCliIjEncKhFXtLKlmzp5RLZgzBzIIuR0Qk7hQOrXhm/X4APjRtUMCViIgEQ+HQiqfW7WPm8L4M698j6FJERAKhcGhmW+ER3t1/mIt0Sm4R6cAUDs08tW4/nQwunKouJRHpuBQOIZxzPLNuH3NGDyAzTaflFpGOS+EQYtO+w+worlCXkoh0eAkZDmY2zsyqzezBeK73mfX7SelknD9lYDxXKyKScBIyHIA7gZXxXunbOw9y0vC+9O3RNd6rFhFJKAkXDmb2caAUeDGe661vaGTTvsNMG9o3nqsVEUlICRUOZpYG/BD4epjz32hmq8xsVVFRUZvWnVtQTk19I9OG6tTcIiIJFQ7A7cA9zrm94czsnLvbOTfLOTcrIyOjTSvekF8KoC0HERHiGA5mttzM3FFur5rZDGAR8Jt41RRqXV4ZvVNTGKGjokVESInXipxzZx/rcTP7CjAS2OOf7K4X0NnMJjnnZsa6vg15ZUwd0odOnXSiPRGRROpWuhsYA8zwb38EngXOi/WKa+obyDmgwWgRkSZx23I4HudcJVDZdN/MyoFq51zbRprDkLP/CHUNToPRIiK+hAmH5pxzt8VrXevzSgEUDiIivkTqVgrM+rwy+vfsypC+3YMuRUQkISgcgA353mC0rvomIuLp8OFQWVtPbsERpqtLSUTkPR0+HDbvO0yjg6naU0lE5D0dPhzW5ZUBGowWEQnV4cNhQ14pWWndyNLFfURE3pOwu7LGy7is3gzSXkoiIh/Q4cPh8/PHBl2CiEjC6fDdSiIi0pLCQUREWlA4iIhICwoHERFpQeEgIiItKBxERKQFhYOIiLSgcBARkRbMORd0DVFhZkXA7giekg4Ux6icRKZ2dyxqd8cTadtHOOcymk9MmnCIlJmtcs7NCrqOeFO7Oxa1u+OJVtvVrSQiIi0oHEREpIWOHA53B11AQNTujkXt7nii0vYOO+YgIiJH15G3HERE5CgUDiIi0oLCQUREWkjacDCz/mb2hJlVmNluM7vqGPN+1cwOmFmZmd1rZt3iWWs0hdtuM7vOzFab2WEzyzOzO8ys3V4ZMJL3O+Q5y8zMted2Q8R/66PN7BkzO2JmxWZ2RzxrjaYI/tbNzH5kZvn+//hyM5sc73qjwcy+YGarzKzGzJYcZ942fa4lbTgAdwK1QBbwCeAPrf1BmNl5wDeBhcBIYDTwg/iVGXVhtRvoAXwF72jKU/Haf0ucaoyFcNsNgJl9guS5TG64f+tdgaXAMmAgMBR4MI51Rlu47/kVwPXAXKA/8AbwQLyKjLJ9wI+Ae481U1Q+15xzSXcDeuL90WSHTHsA+Fkr8z4M/CTk/kLgQNBtiHW7W3nu14Cng25DPNoN9AFygTmAA1KCbkM82g7cCLwSdM0BtPsbwKMh9ycD1UG3oY3t/xGw5BiPt/lzLVm3HLKBBudcbsi0dXh/FM1N9h8LnS/LzAbEsL5YiaTdzc0DNsWkqtiLtN0/Af4AHIh1YXEQSdvnALvM7Dm/S2m5mU2NS5XRF0m7/waMNbNsM+sCXAf8Jw41BqnNn2vJGg69gLJm08qA3mHM2/R7a/Mmukja/R4z+xQwC/hljOqKtbDbbWazgDOA38WhrniI5D0fCnwc+C0wGHgW+Jff3dTeRNLu/cArwBagCq+b6asxrS54bf5cS9ZwKAfSmk1LA46EMW/T763Nm+giaTcAZnYp8DPgfOdcez2LZVjtNrNOwF3Al51z9XGqLdYiec+rgFedc88552rxvgwMACbGtsSYiKTd3wdOAYYBqXh978vMrEdMKwxWmz/XkjUccoEUMxsXMm06rXebbPIfC52vwDl3MIb1xUok7cbMFgN/Bi5yzm2IQ32xEm670/C2kP5uZgeAlf70PDObG/syYyKS93w93hhLMoik3dOBvzvn8pxz9c65JUA/YFLsywxM2z/Xgh5YieGAzd+AR/AGrs7A26ya3Mp8i/H6nifh/cEsI4wB3ES9RdDuBcBBYF7QNcer3YDh7aXTdDsF78NyCNA16DbE4T0fD1QCi4DOeF0r29tr2yNo9/eBV/H2auoEXANUAH2DbsMJtDkFb+vnp3gD8Km0skNFND7XAm9sDF/E/sCT/h/BHuAqf/pwvE2u4SHzfg0oAA4D9wHdgq4/1u0GXgLq/WlNt+eCrj8e73fIc0bSzvdWirTtwGXANv9vfXlrH6bt5RbB33oq3m6v+/12vwMsDrr+E2zzbf7fbOjttlh8runEeyIi0kKyjjmIiEgbKBxERKQFhYOIiLSgcBARkRYUDiIi0oLCQUREWlA4iITBzG4zs41RWI4zs8sjmed490ViQeEgScfMlvgfoM7M6sxsh5n90sx6Bl1bmAYBT4fzuJmN9Ns5Ky6VSYeRLBc7EWnuBbzTJHTBu8jLX/BOs3Bz6Ez+VeAaXAIdDeqcO+apxI/3uEg0aMtBklWNc+6Ac26vc+5h4CHg0qbuITP7pJltB2qAnmY23L/k5BH/9riZDW2+UDO7wcz2mFmVmT1pZukhj51iZs/710o4bGavmtlprdQ20MyeNbNK//KWVzdbxzG7jZo9vtP/udKfvtzM5vlbTAObPe/HZrY+rFdPOjyFg3QUVXhbEQCjgKvwzus/HS8gnsQ7MdsCYD7e9Q6eNDMLWcZI4GrgEryT143jg5dr7I13MrS5wGxgLfDv0ADx/QB4CpgB3A38tQ3dQrP9n4vxupsuc86twDuh3rVNM/mnK78WuOcE1yMdjLqVJOmZ2Wy8MHjRn9QVuMY5V+A/fg5eSIxxzu3yp12Fd4K6hXhdVADdgWudc3v8eW4CXjGzcc65rc65Zc3W+0XgI3gf3KHXan7cOfcn//cfm9l8vOt5f2ALIkxF/s+Dzbqb/gJ8GrjDv38ekEn7vma0xJG2HCRZLTazcjOrxrug/Argi/5jeU3B4JsI7GsKBgDn3A68i7mHnvM/vykYfG8Bjf7zMbNMM/uTmeWaWRnehVUy8c6YGeqNVu5H+9oC9wOjzex0//71wJOufV6nRAKgLQdJViuAG4E6vA/+OgC/l6ii2bzG0S+CE8lA9f14XVNfBXbhdVe9iLelElfOuSIzewq43sy2ABcDF8W7Dmm/tOUgyarSObfNObe7KRiOYTMwxMxGNk0ws9F44w6bQ+YbYmbDQu7Pxvsfete/fybwO+fcs865TXhbDoNaWd+cVu6/28p84aj1f3Zu5bE/Ax8FbsI7r/8Lrcwj0iptOYh4H5rrgIfM7Et4WxK/w7soTOg4QhVwv5l9DW/84Y/As865rf7jucDVZvYW3m6zd/D+h3eoy8xsJd7Fdi7HG9c49QRrL/TrOs/MdgHVzrmmi8kvxbva3/fxrgLWeILrkA5IWw7S4fnHOFyKN7i7HO8qeQeAS5sd/7AL79KUT+OFxg7gUyGPXw/0Alb7893rP6e52/AGqtfjHXfxKefcylbmC6f2euBLwA14YyT/atau+/D20rrvRJYvHZeuBCeSxMzsD8BY59w5Qdci7Yu6lUSSkJn1AU7GO7bhowGXI+2QwkEkOf0Lb8D8Hufcs0EXI+2PupVERKQFDUiLiEgLCgcREWlB4SAiIi0oHEREpAWFg4iItPD/AY3+JEXuynppAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"alphas = [-4, -8,-12,-20]\n",
"betas = [0.4,0.4,0.6,1]\n",
"x = np.arange(40)\n",
"fig = plt.figure(figsize=(8, 6)) \n",
"ax = plt.subplot(111)\n",
"\n",
"for i in range(len(alphas)):\n",
" a = alphas[i]\n",
" b = betas[i]\n",
" y = np.exp(a+b*x)/(1+np.exp(a+b*x))\n",
"# plt.plot(x,y,label=r\"$\\frac{e^{%d + %3.1fx}}{1+e^{%d + %3.1fx}}\\;\\alpha=%d, \\beta=%3.1f$\" % (a,b,a,b,a,b))\n",
" ax.plot(x,y,label=r\"$\\alpha=%d,$ $\\beta=%3.1f$\" % (a,b))\n",
"ax.tick_params(labelsize=12)\n",
"ax.set_xlabel('x', fontsize = 14)\n",
"ax.set_ylabel('$p(x)$', fontsize = 14)\n",
"ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), prop={'size': 16})\n",
"ax.set_title('Logistic Functions', fontsize = 16);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Parameter $\\beta$ controls how fast $p(x)$ raises from $0$ to $1$\n",
"\n",
"The value of -$\\alpha$/$\\beta$ shows the value of $x$ for which $p(x)=0.5$\n",
"\n",
"Another interpretation of $\\alpha$ is that it gives the __base rate__ -- the unconditional probability of a 1. That is, if you knew nothing about a particular data item, then $p(x) = 1/(1+e^{-\\alpha})$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The function $f(x) = \\log (x/(1-x))$ is called the __logit__ function.\n",
"\n",
"So a compact way to describe logistic regression is that it finds regression coefficients $\\alpha, \\beta$ to fit:\n",
"\n",
"$$\\text{logit}\\left(p(x)\\right)=\\log\\left(\\frac{p(x)}{1-p(x)} \\right) = \\alpha + \\beta x$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Note also that the __inverse__ logit function is:\n",
"\n",
"$$\\text{logit}^{-1}(x) = \\frac{e^x}{1 + e^x}$$\n",
"\n",
"Somewhat confusingly, this is called the __logistic__ function."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, the best way to think of logistic regression is that we compute a linear function:\n",
" \n",
"$$\\alpha + \\beta x$$\n",
" \n",
"and then \"map\" that to a probability using the $\\text{logit}^{-1}$ function:\n",
"\n",
"$$\\frac{e^{\\alpha+\\beta x}}{1+e^{\\alpha+\\beta x}}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Logistic vs Linear Regression"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Let's take a moment to compare linear and logistic regression."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In __Linear regression__ we fit \n",
"\n",
"$$y_i = \\alpha +\\beta x_i + \\epsilon_i.$$\n",
"\n",
"We do the fitting by minimizing the sum of squared error ($\\Vert\\epsilon\\Vert$). This can be done in closed form. \n",
"\n",
"(Recall that the closed form is found by geometric arguments, or by calculus)."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Now, if $\\epsilon_i$ comes from a normal distribution with mean zero and some fixed variance, \n",
"\n",
"then minimizing the sum of squared error is exactly the same as finding the maximum likelihood of the data with respect to the probability of the errors.\n",
"\n",
"So, in the case of linear regression, it is a lucky fact that the __MLE__ of $\\alpha$ and $\\beta$ can be found by a __closed-form__ calculation."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In __Logistic regression__ we fit \n",
"\n",
"$$\\text{logit}(p(x_i)) = \\alpha + \\beta x_i.$$\n",
"\n",
"\n",
"with $\\text{Pr}(y_i=1\\mid x_i)=p(x_i).$\n",
"\n",
"How should we choose parameters? "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Here too, we use Maximum Likelihood Estimation of the parameters.\n",
"\n",
"That is, we choose the parameter values that maximize the likelihood of the data given the model."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"$$ \\text{Pr}(y_i \\mid x_i) = \\left\\{\\begin{array}{lr}\\text{logit}^{-1}(\\alpha + \\beta x_i)& \\text{if } y_i = 1\\\\\n",
"1 - \\text{logit}^{-1}(\\alpha + \\beta x_i)& \\text{if } y_i = 0\\end{array}\\right.$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"We can write this as a single expression:\n",
"\n",
"$$\\text{Pr}(y_i \\mid x_i) = \\text{logit}^{-1}(\\alpha + \\beta x_i)^{y_i} (1-\\text{logit}^{-1}(\\alpha + \\beta x_i))^{1-y_i} $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"We then use this to compute the __likelihood__ of parameters $\\alpha$, $\\beta$:\n",
"\n",
"$$L(\\alpha, \\beta \\mid x_i, y_i) = \\text{logit}^{-1}(\\alpha + \\beta x_i)^{y_i} (1-\\text{logit}^{-1}(\\alpha + \\beta x_i))^{1-y_i}$$\n",
"\n",
"which is a function that we can maximize via various kinds of gradient descent."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Logistic Regression In Practice"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"So, in summary, we have:\n",
"\n",
"**Input** pairs $(x_i,y_i)$\n",
"\n",
"**Output** parameters $\\widehat{\\alpha}$ and $\\widehat{\\beta}$ that maximize the likelihood of the data given these parameters for the logistic regression model.\n",
"\n",
"**Method** Maximum likelihood estimation, obtained by gradient descent."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The standard package will give us a correlation coefficient (a $\\beta_i$) for each independent variable (feature).\n",
"\n",
"If we want to include a constant (ie, $\\alpha$) we need to add a column of 1s (just like in linear regression)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Index(['gre', 'gpa', 'rank', 'intercept'], dtype='object')"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['intercept'] = 1.0\n",
"train_cols = df.columns[1:]\n",
"train_cols"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimization terminated successfully.\n",
" Current function value: 0.574302\n",
" Iterations 6\n"
]
}
],
"source": [
"logit = sm.Logit(df['admit'], df[train_cols])\n",
" \n",
"# fit the model\n",
"result = logit.fit() "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
Logit Regression Results
\n",
"
\n",
"
Dep. Variable:
admit
No. Observations:
400
\n",
"
\n",
"
\n",
"
Model:
Logit
Df Residuals:
396
\n",
"
\n",
"
\n",
"
Method:
MLE
Df Model:
3
\n",
"
\n",
"
\n",
"
Date:
Mon, 01 Nov 2021
Pseudo R-squ.:
0.08107
\n",
"
\n",
"
\n",
"
Time:
11:56:41
Log-Likelihood:
-229.72
\n",
"
\n",
"
\n",
"
converged:
True
LL-Null:
-249.99
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
LLR p-value:
8.207e-09
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
z
P>|z|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
gre
0.0023
0.001
2.101
0.036
0.000
0.004
\n",
"
\n",
"
\n",
"
gpa
0.7770
0.327
2.373
0.018
0.135
1.419
\n",
"
\n",
"
\n",
"
rank
-0.5600
0.127
-4.405
0.000
-0.809
-0.311
\n",
"
\n",
"
\n",
"
intercept
-3.4495
1.133
-3.045
0.002
-5.670
-1.229
\n",
"
\n",
"
"
],
"text/plain": [
"\n",
"\"\"\"\n",
" Logit Regression Results \n",
"==============================================================================\n",
"Dep. Variable: admit No. Observations: 400\n",
"Model: Logit Df Residuals: 396\n",
"Method: MLE Df Model: 3\n",
"Date: Mon, 01 Nov 2021 Pseudo R-squ.: 0.08107\n",
"Time: 11:56:41 Log-Likelihood: -229.72\n",
"converged: True LL-Null: -249.99\n",
"Covariance Type: nonrobust LLR p-value: 8.207e-09\n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"gre 0.0023 0.001 2.101 0.036 0.000 0.004\n",
"gpa 0.7770 0.327 2.373 0.018 0.135 1.419\n",
"rank -0.5600 0.127 -4.405 0.000 -0.809 -0.311\n",
"intercept -3.4495 1.133 -3.045 0.002 -5.670 -1.229\n",
"==============================================================================\n",
"\"\"\""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Notice that all of our independent variables are considered significant (no confidence intervals contain zero)."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Using the Model"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Note that by fitting a model to the data, we can make predictions for inputs that were never seen in the data. \n",
"\n",
"Furthermore, we can make a prediction of a probability for cases where we don't have enough data to estimate the probability directly -- e.g, for specific parameter values."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Let's see how well the model fits the data.\n",
"\n",
"We have three independent variables, so in each case we'll use average values for the two that we aren't evaluating."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"GPA:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "-"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAFWCAYAAADtzQ3ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8RElEQVR4nO3deXxU1fnH8c8DhB2CYAAFwr4ICCgRd0TRuqKoFbe61KqttrXaatW2olKrtrX1p61LsVr3Bfd9RVHBBXAJyL6FJMi+hCSQ/fn9cSc6xAlkQjIzmXzfr9e8MnPuuXeeM5ObJ+fcc+81d0dERER21CTeAYiIiCQiJUgREZEIlCBFREQiUIIUERGJQAlSREQkAiVIERGRCJQgJSGZ2YVm5qHHgAjLx4QtPzpOMbqZ3VSP26/8DHrtot7DYZ9F1cdLVer+wcyyzazMzL4OlXU1s1fMbFNonSvruB1XmtlpdbnN0HZ/0Jad1E0xs8vM7GMz22xmpWa22sxeM7PzzKxZWN0Lq3yG+WaWaWa/Cq8XqtvKzPJC9YbXdRslvprtuopIXOUD5wE3VCk/P7SsXcwj+t7BQG4c3z/ceuDkCOWbKp+Y2SjgL8DfgZcIPj+AicARwIXAaiCrjmO7EpgOvFBXG9xJWyLVbQe8CYwEHgitswXoTvCZ/Q8oAZ6psuoZBN9v+9DzfwGdCT6vSqeFlkPwO/m7WjdKEo4SpCS6F4CfmNlED13VwsxaAacDzxP8Ua8TZtbC3YtrWt/dP6ur964DJTWIZ5/Qz/vdfXmV8kx3f7F+QqsX1bUlkn8BGcAR7v55lWVPmtl+QKsI633t7ktDz98xs34EyT48QV5A8E/IEuBcM/u9u5dH0Q5JYBpilUT3GNATOCys7FSgKUGC3IGZHWBmz5lZrpltN7NFZnZrKKmG15tmZtPNbJyZfWVmxcDloWX7h4biiswsJzSUd7OZeZVt7DDEamY3hcr6m9nrZlZgZivNbKKZNQmr19LM7jSzb0J11pjZq2Y2qA4+r4jMbBrwcOjlslCcD4faNAY4PGxIsVdond5m9oSZrTezYjP72sxOjbDt4Wb2opltDPvMrw8tyyL4/s4N2/7DVbdRZXujzOy90GdTaGZTQz3GnbXlpmq21Q34CfCfCMkRAHf/yt0/2VlMIbOAdmbWOWzbY4Gngf8CXYBja7AdaSDUg5REtxL4iGCY9eNQ2fnAi0BBhPrpwNcEf0DzgSEE//H3Ac6qUncAcDfwZ2A5sMnM9gSmAt+G3qcEuAroFUXMLxIM290JjANuBnJCZQAtCIaGbyEY0uxIkJw/M7NB7r4mivf6TtXjYyHloZ735QSJ4nqCYcHVocf9wH+A8lAdgNVm1gP4HFhH0P71wJnA82Y23t1fCb3nKGAasDRULxfoDwwLbetU4A0gE7gpVLZ+J20YBnwIzCcYHXDgOuBDMzvI3TOraUt1Q91jCP6Zeq2694xCb4LPqfL37jyCTsajwAKCnuoFBO2VZODueuiRcA++/+PYD7gI2Ay0BPYCyoBjCP74OXB0Ndswgn8CfwJUAJ3Clk0LlY2oss6tBEmxe1hZK2BtsLvsUNeBm8Je3xQq+2mVenOBd3bS1qZAa4KEflWEz6DXLj6rh0P1Ij2uDqt3caTtERwfnFal7EGCRNapSvm7BEOPla8/Ikj+rXcSXxbweA2/9+cIjg92CCtrTzCM+cKu2hJhe9eG6g2s5nej8tEkwuc+MLRsD+DnBMnxpbB684GFYa+fAorCY9ejYT80xCoNwbMEva5xwLnAGoJe3g+YWXsz+6uZLQOKgVKCYVoj6NmEy3L3r6uUHQR86u7f9UjcfTvwehTxVq37DUHPNjzOCWb2uZltIUj4hUBbgj/KtbEOOCDC47Fabu84gp5Qnpk1q3wAbwPDQ59za+BQ4Al331bL96lqNPCau2+pLHD3rcArBBOJomXVlF9L8LtR+Xg0Qp2FoWWbgHuBJwj+WavsOe/Djp/vIwS/pxNqEackIA2xSsJz93wLTlc4j2Co8wl3rzCL+Lfvf8DRBMOqXxMknlHAPQQ90HCrI6y/F0FCq2ptFCFvqvK6OPy9zWwcwYzJRwiGXzcQ9GbfiBBjTZW6++xarhtJZ4Ih5vOrWd6JoKfdhLqdyduRyN/LGoKeXLRyQj/TgUVh5Q8D74Wev1LNuqcStC0fWOnuRWHLLgj9fNXMOoSezyLodV8ATK5FrJJglCCloXiUoGfWBDg7UgUzawmcQjDseVdY+b7VbDPSvd5WEySHqrpEFe3OnQUsdfcLKwvMLIUgOSSKjQTHfP9azfJvCYaGK4Budfi+m4CuEcq78sN/PGpiGkGMJxEMDwPgwXHeNQBmVlLNut/497NYv2Nmzfn+eHZmhPXSzKxfpHWlYVGClIbiXWAKsMXd51VTpwXBH+3SKuUXRvE+nwFXm1n3ymHW0AzYE6MLd6daEwyrhjuPIPZE8RbBeZ7zQkPMEZnZdILTcCbtpF4xkU+jiORD4EQza+fu+aH3aEcwvD6tpsFXcvdVZvYE8HMze9KrmckapXEE/8zcHCGmLgSzWs9nx9NBpAFSgpQGwYNzyyL2HMPq5JnZZ8DvzGw1wdDlRUTXw/kncBnwtpndTPDH/behn3V1d/G3gPFmdifB7MqRwBUEk1Nqq7mZHRShfJu7z6nF9iYCM4GPzOzfBBNt9gCGAn3c/aJQvasJktqnZvYPgiHJPgSTn34dqjOf4DSSkwh6bRvcPaua9/0zQW9vqpn9leAzv5bgn4pJtWgHwK8Ijj9/YGYPEAytbiZIcqMJeqfVXmggggsIZrLe4e4/mEltZlcB55vZje6uO9I3YJqkI8nmbOALgmOODxP8Qf5NTVd29w0E57ZtJhjWvZfgD+qLQF4dxfgAwVVgzgReJeidjtvN7acBn0Z4PFmbjbl7NsHJ9ZkEM3vfBe4jmCjzfli9WQQTdXIITnN4A7iGHY9LXk9w/G8KwXG6m3byvnMIZidvJThG+xhBMjrCg1M8atOWraG4rwm16THgA4LTWwYBPyNIortkZmnA8cCUSMkx5EGCcz9rM6lIEojpHxyRnTOzpsCXBD2fsfGOR0RiQ0OsIlWY2Z8JTnxfSTBb82KCE99PiGdcIhJbSpAiP+QEx+D2Dj2fA4x39zfjGpWIxJSGWEVERCLQJB0REZEIlCBFREQiaFTHIPfcc0/v1atXvMMQEZEE8cUXX2xw97RIyxpVguzVqxezZ9fl5SpFRKQhM7OV1S3TEKuIiEgESpAiIiIRKEGKiIhEoAQpIiISgRKkiIhIBEqQIiIiEShBioiIRKAEKSIiDcucKXDnULipQ/BzzpR6eZtGdaEAERFp4OZMgVevgNLtweu8nOA1wLAJdfpW6kGKiMRDjHpBSWfqpO+TY6XS7UF5HVMPUkQk1mLYC0o6ebnRle8G9SBFRGIthr2gpJPaPbry3aAEKSISazHsBSWdsRMhpdWOZSmtgvI6pgQpIhJrMewFJZ1hE2Dc3ZDaA7Dg57i762VoWscgRURibezEHY9BQr31gpLSsAkxOVarHqSISKzFsBcktacepIhIPMSoFyS1px6kiIhIBDFNkGbW0cxeNLNCM1tpZufUYJ33zczNrFlY2TQzKzKzgtBjUf1GLiIijU2se5D3ACVAF+Bc4D4zG1JdZTM7l+qHgX/l7m1Dj4F1H6qIiDRmMUuQZtYGOB24wd0L3H068ApwXjX1U4Ebgd/HKkYREZFKsexBDgDK3X1xWFkmUF0P8lbgPmBNNctvM7MNZjbDzMZU96ZmdqmZzTaz2evXr69F2CIi0hjFMkG2BfKqlOUB7apWNLMM4FDgX9Vs61qgD9ANmAy8amZ9I1V098nunuHuGWlpabWNXUREGplYJsgCoH2VsvZAfniBmTUB7gV+4+5lkTbk7p+7e767F7v7I8AM4IR6iFlERBqpWCbIxUAzM+sfVjYcmFelXnsgA3jGzNYAs0LluWZ2eDXbdsDqMlgREWncYnahAHcvNLMXgElmdjEwAjgFOKRK1Txg77DXPYCZwEhgvZl1AA4EPgTKgDOB0cCV9Ri+iIg0MrG+ks7lwEPAOmAjcJm7zzOzdGA+MNjdswmbmGNmLUNP17p7WWh26y3AIKAcWAiMd3edCykiInUmpgnS3TcB4yOUZxNM4om0ThZhw6fuvh44oH4iFBERCehScyIiIhEoQYqIiESgBCkiIhKBEqSIiEgESpAiIiIRKEGKiIhEoAQpIiISgRKkiIhIBEqQIiIiEShBioiIRKAEKSIiEoESpIiISARKkCIiIhEoQYqIiESgBCkiIhKBEqSIiEgESpAiIiIRKEGKiIhEoAQpIiISgRKkiIhIBEqQIpKc5kyBO4fCTR2Cn3OmxDsiaWCaxTsAEZE6N2cKvHoFlG4PXuflBK8Bhk2IX1zSoKgHKSLJZ+qk75NjpdLtQblIDSlBikjyycuNrlwalLLyCu6btox1W4vq9X00xCoiySe1ezCsGqlcGrSsDYVcNeVrvsreQrMmxiWj+9Tbe6kHKSLJZ+xESGm1Y1lKq6BcGiR356mZ2Zxw98csW1fAXWeNqNfkCOpBikgyqpyIM3VSMKya2j1Ijpqg0yCtzy/m+hfm8N6CdRzarxN///Fw9u7Qatcr7iYlSBFJTsMmKCEmgXfnr+W65+eQX1zGxJMGc+EhvWjSxGLy3jEdYjWzjmb2opkVmtlKMzunBuu8b2ZuZs3CyqLejoiINBwFxWVc+9wcLnl0Nl3at+S1Xx/GRYf1jllyhNj3IO8BSoAuwAjgdTPLdPd5kSqb2blEjjGq7YiISMPxxcpNXPVMJjmbt3HZmL5cdfQAmjeL/ZSZmCVIM2sDnA4MdfcCYLqZvQKcB1wXoX4qcCNwPvBpbbcjIiINQ2l5BXe9t4R7py1l7w6teObSgxnVu2Pc4ollD3IAUO7ui8PKMoEjqql/K3AfsGY3tyMiIglu6bp8rnomk7mr8jhjZHcmjhtMu5YpcY0plgmyLZBXpSwPaFe1opllAIcCvwGqnrhU4+2EtnUpcClAenp61EGLiEj9qahwHvtsJbe+sYDWzZty/0/257ihe8U7LCC2CbIAaF+lrD2QH15gZk2Ae4HfuHuZ2Q8OyNZoO5XcfTIwGSAjI8NrFbmIiNS5NXlFXPNcJh8v2cCYgWn87cfD6NyuZbzD+k4sE+RioJmZ9Xf3JaGy4UDViTXtgQzgmVBybBoqzzWzM4Ava7gdERFJUK/PWc0fXpxLcVk5t4wfyrkHphOhQxRXMUuQ7l5oZi8Ak8zsYoLZp6cAh1SpmgfsHfa6BzATGAmsd/eSGm5HREQSzNaiUm58eR4vfrWK4d1TufPMEfRJaxvvsCKK9WkelwMPAeuAjcBl7j7PzNKB+cBgd88mbGKOmVX2t9e6e9nOthOjNoiISC18tnwjv5uSyZqtRfxmbH9+dVQ/Upom7hVPY5og3X0TMD5CeTbB5JtI62QBVqUs4nZERCTxFJeV8493FvPAx8vp1akNz/3iYPZL3yPeYe2SLjUnIiL1ZsHqrVz1zNcsXJPPuQem88cT96F184aRehpGlCIi0qBUVDj/nb6cO95eTPtWKTx0YQZHDeoS77CiogQpIiJ1KnfzNq5+NpPPlm/iR4O7cNtp+9KpbYt4hxU1JUgREdnRnCm1ulWYu/PS16uY+NI8Ktz524+HccbI7gl3+kZNKUGKiMj35kyBV6+A0u3B67yc4DXsNElu2VbCH1/8htfnriaj5x7ceeYIenRsHYOA648SpIiIfG/qpO+TY6XS7UF5NQnyo8XrufrZTDZvK+H3xw3k56P70jSGt6WqL0qQIiLyvbzcGpdvLynn9jcX8MinK+nXuS0PXXgAQ7ul1nOAsaMEKSIi30vtHgyrRioPMzc3jyuf+Ypl6wv56aG9uPa4QbRMafrD9RqwxL2EgYiIxN7YiZDSaseylFZBOVBWXsG/31/CqffOoLC4nMd/diA3jhuSdMkR1IMUEZFwlccZI8xiXbmxkKue+Zovs7dw0rC9uGX8UDq0bh7feOuREqSIiOxo2IQdJuS4O8/MzGbSa/Np2sS466wRnDKiWxwDjA0lSBERqdaGgmKue34u7y1Yy8F9OvGPCcPZu0OrXa+YBJQgRUQkovfmr+W6F+awtaiMP524Dxcd2psmSXD6Rk0pQYqIyA62FpXyl9cW8MzsHPbZqz1PXDyCgV3bxTusmFOCFBGR73y8ZD3XPjeHNVuL+MURfbnqmP60aJZ8M1RrQglSREQoKC7j1jcW8OTn2fRJa8Pzlx3SIO7ZWJ+UIEVEGrlPlm7gmufm8G3edi4d3YffHjMgKc9rjJYSpIhII1VYXMZf31rIo5+upPeebXjuFwczsmfHeIeVMJQgRUQaoc+Wb+Sa5zLJ3bydnx3Wm6t/NJBWzdVrDKcEKSLSiGwrKeNvby3i4U+y6NmpNc9cejCjeqvXGIkSpIhIIzEraxPXPJtJ1sZtXHhIL35/3EBaN1caqI4+GRGRJFdUWs4dby/iwRkr6NahFU9dchAH9+0U77ASnhKkiEgS+2LlZq55NpPlGwo576CeXHf8INq00J/+mtCnJCKShIpKy7nz3cU88PFy9kptxRMXH8ih/faMd1gNihKkiEiS+TpnC1c/m8nSdQWcPSqdP5wwiHYtU+IdVoOjBCkikiSKy8q5670l3P/hMrq0b8mjF41i9IC0eIfVYClBiogkgTm5Qa9x8doCJmR0508nDaa9eo27RQlSRKQBKymr4F/vL+HeacvYs21z/vfTAzhyYOd4h5UUmsTyzcyso5m9aGaFZrbSzM6ppt5ZZrbIzPLMbJ2ZPWJm7cOWTzOzIjMrCD0Wxa4VIiKJ4ZtVeZz87+n86/2ljB/RjXeuPELJsQ7FNEEC9wAlQBfgXOA+MxsSod4M4FB3TwX6EPR0b6lS51fu3jb0GFifQYuIJJLS8gr+773FjL9nBhsLS3jwggz+MWE4qa01pFqXYjbEamZtgNOBoe5eAEw3s1eA84Drwuu6e06V1cuBfjEJVEQkgS1YvZXfTclk/uqtnLpfN24cN5gOrZvHO6ykFMtjkAOAcndfHFaWCRwRqbKZHQa8DrQHtgGnVqlym5ndDiwC/uju0+o8YhGRBFFWXsH9Hy7jrqlLSG2Vwn/OG8mxQ7rGO6ykFssE2RbIq1KWB7SLVNndpwOpZtYNuATIClt8LTCfYLj2LOBVMxvh7suqbsfMLgUuBUhPT9/NJoiIxN6iNflc/Wwmc1flMW743tx88hA6tlGvsb7F8hhkAUFvMFx7IH9nK7n7KuAt4Omwss/dPd/di939EYJjlidUs/5kd89w94y0NJ0PJCINR1l5BfdOW8q4f03n2y3bue/c/fnX2fspOcZILHuQi4FmZtbf3ZeEyoYD82qwbjOg706WO2C7GZ+ISMJYui6f3z07h8ycLZywb1f+fMpQOrVtEe+wGpWYJUh3LzSzF4BJZnYxMAI4BTikal0zOxf4GMgB0oG/AFNDyzoABwIfAmXAmcBo4Mr6boOISH0rr3D++/Fy/vHuYto0b8q/z9mPk4btHe+wGqVYXyjgcuAhYB2wEbjM3eeZWTrBMcXB7p4NDAb+CuwBbAbeAK4PbSOF4JSPQQSzWxcC491d50KKSIO2bH0B1zybyZfZWzh2SBduGb8vae3Ua4wXc/d4xxAzGRkZPnv27HiHISKyg/IK538zVvD3txfRMqUpk04ZwsnD98ZMR47qm5l94e4ZkZbpUnMiInGUtaGQa57LZFbWZo7epzO3nrovndu3jHdYghKkiEhclFc4j36axV/fWkjzpk3454ThnLpfN/UaE4gSpIhIjC1em8+1z8/hq+wtHDkwjdtOG0bXVPUaE80uE6SZdYx2o+6+qXbhiIgkr+Kycu75YBn3TVtK2xbN+L8zR3DKCB1rTFQ16UFuIDjPsKbczAa4+/JaxpS45kyBqZMgLxdSu8PYiTBsQryjEpEGYHbWJq59fg7L1hdy6n7d+NOJ++i8xgRX0yHWHwM16RUawSkZyWfOFHj1CijdHrzOywleg5KkiFQrv6iUv761kMc/y6Zbh1Y8ctEojhigq3o1BDVJkCuBj9x9Y002aGbLgdLdiioRTZ30fXKsVLo9KFeCFJEI3pm3hokvz2NdfhE/O6w3vz1mAG1aaOpHQ7HLb8rde0ezQXcfWvtwElhebnTlItJorcsv4qZX5vHG3DUM6tqO/5w3kuE9OsQ7LIlSVP/KmNn5wDPuXlylvDlwlrs/WpfBJZTU7sGwaqRyERHA3ZkyO4e/vL6AorIKrjl2IJeO7kNK01jfm17qQrTf2v+A1Ajl7ULLktfYiZDSaseylFZBuYg0eis2FHLOA59z7fNz2Wev9rz1m8P55ZH9lBwbsGgHw43IM1rT+eG9HpNL5XFGzWIVkTCl5RU88PFy7npvCc2bNeG20/blzIweNGmiUzcauholSDObS5AYHfjQzMrCFjcFepKss1fDDZughCgi35mTu4Vrn5/LgtVbOX5oV24+eYguE5dEatqDfC70cyjwOsHNjyuVAFnA83UXlohI4tpWUsY/31nMQzNWkNauBf85byTHDuka77CkjtUoQbr7zQBmlkUwSaeoPoMSEUlUHy1ezx9enEvu5u2ce2A61x4/iPYtU+IdltSDqI5Buvsj9RWIiEgi21RYwi2vzeeFr1bRJ60NU35+MKN6R30lTmlAanIt1q1AH3ffYGb57OSyc+7evi6DExGJN3fnlcxvufnV+WzdXsqvj+rHL4/sR8uUpvEOTepZTXqQvwbyw543njssi0ijlrt5G3966RumLVrPiB4duP30fRnUVf2AxqImV9J5JOz5w/UajYhIAiivcB75JIs73lkEwI3jBnP+wb1oqlM3GpVaXxTQzFpS5UID7r5ttyMSEYmjhWu2cu3zc8nM2cKYgWncMn4o3fdoHe+wJA6ivdRcT+Bu4EigTYQqGpQXkQapqLScf7+/lPs/XEb7VincddYITh6uezU2ZtH2IB8HWhIci1yLjkeKSBL4fPlGrn9hLss3FHLa/t3404mD6dimebzDkjiLNkHuBxzg7gvqIxgRkVjaWlTK7W8u5MnPs+nRsRWP/WwUh/fXvRolEG2CzATSACVIEWnQ3vpmDRNf/oYNBcVccnhvrjpmAK2b616N8r1ofxsuBe42s7uBb6hyY2R3z66rwERE6sParUXc+PI83pq3hn32as9/L8hgWPcO8Q5LElC0CbIJ0Bl4kR2PP1be5UOTdEQkIVVUOE/PyuG2NxdQUlbBtccN4uLDe+t2VFKtaBPkI8B6YByapCMiDcSy9QVc/8JcZq7YxMF9OnHrafvSe89IE/FFvhdtghwEjHD3xfURjIhIXSotr2DyR8u5a+oSWjZrwt9OH8YZGd116obUSLQJcibQG1CCFJGE9smyDUx8eR5L1xVw4r57cePJg+ncTvdqlJqLNkHeB/yfmf0DmMsPJ+l8WVeBiYjUxrqtRdzy+gJeyfyWHh1b8eAFGYzdp0u8w5IGKNoE+VTo5+QIy3Y5ScfMOgIPAj8CNgDXu/uTEeqdBdwMdAWKgTeBX7v71mi2IyKNR1l5BY98upI7311MSXkFV4ztz+Vj+uquG1Jr0SbI3rv5fvcAJUAXYATwuplluvu8KvVmAIeGbrHVFvgPcAtwRZTbEZFGYFbWJm546RsWrsnniAFp3HzyEHppEo7spmhvmLyytm9kZm2A04Gh7l4ATDezV4DzgOuqvE9OldXLgX7RbkdEktv6/GJue3MBL3y5im4dWnH/T0Zy7JAumoQjdaImN0w+v6Ybc/dHd7J4AFBeZQZsJnBENe97GPA60B7YBpxam+2ISPIpK6/gic+zueOdRRSVlnP5mL786qh+uhKO1Kma/DbdU+V1cyAFqAi9bkIwWacY2FmCbAvkVSnLA9pFquzu04FUM+sGXAJk1WY7ZnYpwRWASE9P30l4ItIQfLFyMze89A3zV2/lsH57cvMpQ+ib1jbeYUkS2uUlJNy9XeUDOAuYAxxOcFePlqHnXwPn7GJTBQS9wXDtgfxdvP8q4C3g6dpsx90nu3uGu2ekpekixCIN1caCYn7/XCan3/cJmwpLuOec/XnsZ6OUHKXeRDsecQdwkbt/GlY2w8yuBB4GXtvJuouBZmbW392XhMqGAzWZWNMM6FsH2xGRBqa8wnl6VjZ/e2sRhcVl/Hx0H64Y2582LTScKvUr2t+wXkBhhPJtwE7HL9290MxeACaZ2cUEs09PAQ6pWtfMzgU+BnJC2/0LMDXa7YhIw5aZs4UbXv6GObl5HNSnI5NOGcqALhGPpojUuWiv0vs5wd08ulUWhJ7fCXxWg/UvB1oB6wjOqbzM3eeZWbqZFZhZZZIdDHxCMJw6A1hEcBxyp9uJsi0ikqC2bCvhDy/OZfy9M1idV8RdZ43gqUsOUnKUmDL3ml9v3Mz6Ai8RXJN1Vai4G0ECG+/uS+s6wLqUkZHhs2fPjncYIlKNigrn2S9yuP3NhWwtKuOCg3tx1TH9adcyJd6hSZIysy/cPSPSsmjPg1xmZsOAYwiSpAHzgfc8mkwrIlLFN6vymPjyN3yZvYUDeu3BpFOGss9eVefjicRO1Ee5Q4nwndBDRGS35G0v5Z/vLOKxz1ayR+vm3HHGcE7fv5tO9pe4izpBhq6DehzB5Jnm4cvcfVIdxSUiSc7deeHLVdz25gI2FZbwk4N68rsfDSS1lYZTJTFElSDN7CCCq9sUA2kExyH3Cr3OApQgRWSXFq7Zyg0vfcOsrM2M6NGBh386iqHdUuMdlsgOou1B/h14AvgNsBU4iuC0j6cI7q4hIlKt/KJS7nx3CY98mkX7ls346+n7csbIHjRpouFUSTzRJshhwM/c3c2sHGjh7svN7FrgSYLkKSKyA3fnlcxvueX1BWwoKOasA9L5/bED2aNN812vLBIn0SbIkrDna4GewAKC8xX3rqugRCR5LFmbz8SX5/Hp8o3s2y2VB87PYESPDvEOS2SXok2QXwIHEFzubRpwi5l1AX5CcI1WEREACovLuHvqEh6cvoI2LZpxy/ihnD0qnaYaTpUGItoE+Ue+v2vGnwju3vEvgoT50zqMS0QaKHfnjblr+PNr81mztYgJGd259rhBdGrbIt6hiUQl2gsFzA57vh44vs4jEom1OVNg6iTIy4XU7jB2IgybEO+oGqRl6wu46ZV5fLxkA4P3as895+7PyJ57xDsskVrR5fClcZszBV69Akq3B6/zcoLXoCQZhW0lZfz7/aU88PFyWqY05eaTh3Dugek0axrt5Z5FEscuf3vN7A0zq/EJSmb2fOi4pEjimzrp++RYqXR7UC675O689c0ajvnnR9w7bRnjhu3N+78bwwWH9FJylAavJj3IY4GuZta0hts8BmhT+5BEYigvN7py+c43q/K45fX5fLZ8EwO7tOOZSw/iwD6d4h2WSJ2pSYKsvCC5SPJJ7R4Mq0Yql4jWbi3i728v4vkvc+nQKoU/nzKEs0alk6IeoySZmiTII2ux3VW7riKSAMZO3PEYJEBKq6BcdrC9pJzJHy3n/g+XUVZRwSWH9+GXR/bTtVMlae0yQbr7h2bWmuAyc+OBFOA94Ap331C/4YnUs8qJOJrFWq2KCuelr1fxt7cWsWZrEccP7cp1xw+iZycdSZHkVtNZrDcDFxJcSm47cA5wH3BG/YQlEkPDJighVmPmik3c8vp85uTmsW+3VO4+ez9G9e4Y77BEYqKmCfI0gmuwPg1gZk8AM8ysqbuX11t0IhIX2Ru3cdubC3jzmzV0bd+Sf04YzvgR3XRRcWlUapogewAfV75w95lmVkZw/dUIMxxEpCHaWlTKPe8v5X8zsmjaxLjq6AFcMro3rZvrlGlpfGr6W9+UHS9UDlAWxfoiksDKyit4alYOd767mM3bSjh9/+5c/aOBdE1tGe/QROKmpgnOgMfNrDisrCXwgJltqyxw95PrMjgRqX/TFq3jL68vYMm6Ag7s3ZEbThqsmxeLUPME+UiEssfrMhARia1Fa/L5yxsL+Gjxenp2as39PxnJsUO6YKbjjCJQwwTp7rpTh0iS2FBQzJ3vLuapmdm0adGMP524D+cf3IvmzXSiv0g4HUMUaSSKSst5+JMs7nl/KdtKyznvoJ785ugBdGzTPN6hiSQkJUiRJFd5f8bb31pAzqbtHDWoM384YR/6dW4b79BEEpoSpEgS+zpnC7e8Np/ZKzczqGs7HvvZKA7vnxbvsEQaBCVIkST07Zbt/O2thbz09bfs2bY5t522LxMyetBUJ/qL1JgSpEgSKSwu4/4PlzH5o+U4cPmYvlw2pi/tWuqC4iLRimmCNLOOwIPAj4ANwPXu/mSEehcAVwD9ga3Ak8Af3L0stHwacBDBxQoAVrn7wHpvgEiCKq9wnv8il7+/s4j1+cWMG7431x43kO57tI53aCINVqx7kPcQXJGnCzACeN3MMt19XpV6rYErgc+BNOAV4Grg9rA6v3L3/9Z3wCKJ7pNlG7jltQXMX72VET06cP9PRjKy5x7xDkukwYtZgjSzNsDpwFB3LwCmm9krwHnAdeF13f2+sJerQhdHr819KUWS1vL1Bdz6xkLeW7CWbh1acffZ+zFu2F460V+kjsSyBzkAKHf3xWFlmcARNVh3NFC1l3mbmd0OLAL+6O7T6iRKkQS3ZVsJd09dyqOfZtGiWROuOXYgPzusNy1TmsY7NJGkEssE2RbIq1KWB7Tb2Upm9lMgA7g4rPhaYD7BcO1ZwKtmNsLdl0VY/1LgUoD09PRaBy8Sb6XlFTz26UrumrqE/KJSzjygB1cdM4DO7XRBcZH6EMsEWQC0r1LWHsivbgUzG09w3PFod99QWe7un4dVe8TMzgZOAP5VdRvuPhmYDJCRkeG1DV4kXtyd9xas47Y3FrB8QyGH9uvEn04czD57Vd2dRKQuxTJBLgaamVl/d18SKhvOD4dOATCz44AHgBPdfe4utu0EdxwRSSozV2zi728vZFbWZvqkteHBCzI4alBnHWcUiYGYJUh3LzSzF4BJZnYxwSzWU4BDqtY1s6OAJ4BT3X1mlWUdgAOBDwlO8ziT4BjllfUYvkhMzc3N4453FvHh4vWktWvBn08Zwlmj0klpqguKi8RKrE/zuBx4CFgHbAQuc/d5ZpZOcExxsLtnAzcAqcAbYf8pf+zuxwMpwC3AIKAcWAiMd/dFMW2JSD1Yui6ff767mDfmrqFD6xT+u99yjlp1P03eXgWfdYexE2HYhHiHKdIoxDRBuvsmYHyE8myCSTyVr6s9pcPd1wMH1Ed8IvGSs2kbd01dwgtf5tIqpSlXjO3PLzp+Qeu3/gKl24NKeTnw6hXBcyVJkXqnS82JxNG6/CLueX8pT87Mxsy46NDeXDamL53atoA7T/s+OVYq3Q5TJylBisSAEqRIHGzZVsL9Hy7n4U9WUFruTMjowRVj+7FXaqvvK+XlRl65unIRqVNKkCIxVFBcxv+mr2DyR8spKCnj5OF7c9XRA+i1Z5sfVk7tHgyrRioXkXqnBCkSA0Wl5TzxeTb3frCUjYUlHL1PF373owE7P5dx7MTgmGP4MGtKq6BcROqdEqRIPSotr+D5L3K5a+oSVucVcUjfTlx97ED2T6/BxcQrjzNOnRQMq6ZqFqtILClBitSDigrntbmrufPdxazYUMiIHh2444zhHNpvz+g2NGyCEqJInChBitQhd2fqgnXc8c4iFq7JZ2CXdjxwfgZH76Or34g0NEqQInXkk2Ub+Pvbi/gqews9O7XmrrNGcNKwvWnaRIlRpCFSghTZTV/nbOGOtxcxfekGurZvya2n7ssZGd11WTiRBk4JUqSWFq3J5x/vLOKd+WvZo3UKfzpxH35yUE/dl1EkSShBikRp5cZC/u+9Jbz09SraNm/GVUcP4KLDetGuZUq8QxOROqQEKVJDa/KKuPv9JUyZlUPTJsalo/vwi9F92aNN83iHJiL1QAlSZBc2FZZw37SlPPrpSsornLNHpfOro/rRpX3LeIcmIvVICVKkGvlFpfz34xU8OH0FhSVlnLpfN64cO4D0Tq3jHZqIxIASpEgVRaXlPPppFvdOW8aWbaUcN6Qrv/3RAAZ0aRfv0EQkhpQgRUJKyip4ZnYO/35/CWu3FjN6QBpX/2gAw7p3iHdoIhIHSpDS6JWWV/Dy199y99QlZG/axsiee3DXWftxUJ9O8Q5NROJICVIaraLScp77Ipf7P1xG7ubtDN6rPf+78ADGDEzTZeFERAlSGp9tJWU8+Xk2kz9azrr8Yob36MCN44YwdlBnmuiycCISogQpjUbe9lIe+zSLh2ZksamwhIP6dOSfE0ZwaL9O6jGKyA8oQUrS21hQzEMzVvDoJyvJLy7jyIFp/Oqofozs2THeoYlIAlOClKS1Jq+IyR8t56mZ2RSVlXP80K5cPqYfQ7ulxjs0EWkAlCAl6eRs2sZ9Hy7judm5lLtzyvC9ufzIvvTrrPMYRaTmlCAlaSxdl8+9Hyzj5cxvaWrGjzO684vRfXXlGxGpFSVIafC+WZXHPR8s5a15a2jZrCkXHtKLSw7vQ9dUXStVRGpPCVIarNlZm/j3B0uZtmg97Vo045dj+vHTQ3vRqW2LeIcmIklACVIaFHdnxtKN/PuDJXy2fBMd2zTnmmMHct7BPWmv+zGKSB1SgpQGoaLCmbpwHf/+YCmZOVvo0r4FN5w0mLNH9aB1c/0ai0jd018WSWjlFc7rc1dz7wdLWbgmnx4dW3Hrqfty+shutGjWNN7hiUgSaxLLNzOzjmb2opkVmtlKMzunmnoXmNkXZrbVzHLN7G9m1iza7UjDVVJWwZRZORz9zw+54qmvKKtw/jlhOB/8bgznHJiu5Cgi9S7WPch7gBKgCzACeN3MMt19XpV6rYErgc+BNOAV4Grg9ii3Iw1MUWk5U2bn8J8Pl7Nqy3aG7N2e+87dn2OHdNV1UkUkpmKWIM2sDXA6MNTdC4DpZvYKcB5wXXhdd78v7OUqM3sCODLa7UjDUVBcxhOfreSBj1ewoaCYjJ57cMupQxkzQHfWEJH4iGUPcgBQ7u6Lw8oygSNqsO5ooLJ3uDvbkQSzZVsJD3+Sxf9mZJG3vZTD++/JL4/cjwN7d1RiFJG4imWCbAvkVSnLA3Z6/S8z+ymQAVxcm+2Y2aXApQDp6enRRSz1Zn1+Mf+dvpzHP11JYUk5xwzuwi+P7MeIHh3iHZqICBDbBFkAtK9S1h7Ir24FMxtPcNzxaHffUJvtuPtkYDJARkaGRx211Klvt2z/7gLipeUVnDQsuE7qoK5Vv1IRkfiKZYJcDDQzs/7uviRUNpzvh053YGbHAQ8AJ7r73NpuRxLD0nX5PPDRCl74Khd3OG3/blw2ph+992wT79BERCKKWYJ090IzewGYZGYXE8w+PQU4pGpdMzsKeAI41d1n1nY7El/uzkdLNvDg9BV8tHg9LZo14ZxR6Vx6RF+6dWgV7/BERHYq1qd5XA48BKwDNgKXufs8M0sH5gOD3T0buAFIBd4Im6jxsbsfv7PtxK4ZsjPbS8p54atc/jcji6XrCjiv9efMSX2GdsVrseXdofdE6DAh3mGKiOxUTBOku28CxkcozyaYfFP5+sjabEfia01eEY9+msWTM7PZsq2Ufbul8uwhuWTMnYwVbw8q5eXAq1cEz4cpSYpI4tKl5mS3ZeZs4aEZK3h9zmoq3Dl2SFcuOqw3GT33wP7vF1C6fccVSrfD1ElKkCKS0JQgpVbKyit4Z/5aHpy+gi9WbqZdi2ZceEgvLjikFz06ht2gOC838gaqKxcRSRBKkBKVvO2lPDMrm0c+WcmqLdtJ79iaG8cN5oyMHrRtEeHXKbV7MKwaqVxEJIEpQUqNrNhQyMMzVvDsF7lsKynnoD4duXHcYMbu04WmO7tG6tiJwTHH8GHWlFZBuYhIAlOClGq5O58u28iD01fw/qJ1pDRpwrjhe/PTQ3sxtFtqzTZSeZxx6qRgWDW1e5AcdfxRRBKcEqT8QFFpOa9kfstD01ewcE0+ndo054qj+nPuQel0btcy+g0Om6CEKCINjhKkfGddfhGPf5bNE5+tZGNhCYO6tuNvPx7GycP3pmWK7r8oIo2LEqQw79s8HpqexauZ31JaUcHYQZ256NDeHNy3k+6oISKNlhJkI1Ve4UxdsJaHZqzgs+WbaN28KWeP6sGFh/bW9VFFRFCCbHQKist4dnYOD3+SxcqN2+jWoRV/OGEQZ2akk9o6Jd7hiYgkDCXIRiJn0zYe/iSLKbNyyC8uY2TPPfj9sYM4dkgXmjVtEu/wREQSjhJkEnN3ZmVt5qHpK3hn/hqamHHCvntx0WG9dWNiEZFdUIJMQiVlFbw+91semp7F3FV5pLZK4edH9OX8g3uyV6puMyUiUhNKkElkY0ExT36ezaOfrWR9fjF909rwl1OHctp+3WnVXKdpiIhEQwmygauocD5bvpEnZ2bz9rw1lJY7owek8fcf92J0/zSa7OwycCIiUi0lyAZqQ0Exz32Ry9Mzs8nauI3UVin85KCenDMqnf5d2sU7PBGRBk8JsgGpqHBmLNvAUzOzeWfeWsoqnFG9OvKbo/tz/NC9dLUbib05U3SdXUlaSpANwLr8Ip6dncvTs7LJ2bSdDq1TuPCQXpw1qgf9Oqu3KHEyZ8qOd2rJywleg5KkJAUlyARVXuF8vGQ9T83MZuqCdZRVOAf16cjVPxrIsUO6qrco8Td10o63MYPg9dRJSpCSFJQgE8zarUVMmZXD07NyWLVlOx3bNOdnh/XmzAN60CetbbzDE/leXm505SINjBJkAiivcD5cvI6nZubw/sJ1lFc4h/brxPUnDOKYwV1o0Uy9RUlAqd2DYdVI5SJJQAkyjr7dsp0ps3OYMiuHb/OK2LNtcy45vA9nHdCDXrpguCS6sRN3PAYJkNIqKBdJAkqQMVZWXsG0RcGxxQ8WraPC4fD+e3LDSYMZu08XmjfTdVGlgag8zqhZrJKklCBjJHfzNqbMymHK7FzWbC0irV0LLh/TjzMP6EGPjq3jHZ5I7QyboIQoSUsJsh6Vllfw/sJ1PDUzmw8XrwfgiAFp3HzKEI4a1JkU3UVDRCRhKUHWg5xN23hmVg5TZuewLr+Yru1b8usj+zHhgB5030O9RRGRhkAJso6Ullfw3vy1PDkzm+lLN2DAkQM7c/aodMYMTNM9F0VEGhglyN20cmMhT8/K4dnZuWwoKGbv1JZcOXYAEw7orltLiYg0YDFNkGbWEXgQ+BGwAbje3Z+MUG8o8A9gJNDJ3a3K8mnAQUBZqGiVuw+sx9B3UFJWwTvz1/D0zBymL91A0ybGUYM6c86odEYPSKOp7qAhItLgxboHeQ9QAnQBRgCvm1mmu8+rUq8UmALcC7xUzbZ+5e7/rac4I1q3tYgHp6/guS9y2VhYQrcOrfjdMQM4I6MHXVNbxjIUERGpZzFLkGbWBjgdGOruBcB0M3sFOA+4Lryuuy8CFplZv1jFVxNFpRU8NGMFYwd14ewD0zms357qLYqIJKlY9iAHAOXuvjisLBM4opbbu83MbgcWAX9092m7Gd8upXdqzew/HkNq65T6fisREYmzWE6tbAvkVSnLA2pzv6ZrgT5AN2Ay8KqZ9Y1U0cwuNbPZZjZ7/fr1tXirHSk5iog0DrFMkAVA+ypl7YH8aDfk7p+7e767F7v7I8AM4IRq6k529wx3z0hLS4s6aBERaZximSAXA83MrH9Y2XCg6gSd2nCgcR8MnDMF7hwKN3UIfs6ZEu+IREQatJglSHcvBF4AJplZGzM7FDgFeKxqXQu0BJqHXrc0sxah5x3M7NhQWTMzOxcYDbwdq7YknMo7u+flAP79nd2VJEVEai3Wl3e5HGgFrAOeAi5z93lmlm5mBWaWHqrXE9jO973L7QSTcQBSgFuA9QTnUv4aGB+a+do47ezO7iIiUisxPQ/S3TcB4yOUZxNM4ql8nUU1Q6buvh44oH4ibKB0Z3cRkTqnC4Qmg+ru4K47u4uI1JoSZDIYOzG4k3s43dldRGS3KEEmg2ETYNzdkNoDsODnuLt1I1sRkd2gu3kkC93ZXUSkTqkHKSIiEoESpIiISARKkCIiIhEoQYqIiESgBCkiIhKBEqSIiEgESpAiIiIRKEGKiIhEYO4e7xhixszWAyt3YxN7EtxBJBmoLYlJbUlMaktiqou29HT3tEgLGlWC3F1mNtvdM+IdR11QWxKT2pKY1JbEVN9t0RCriIhIBEqQIiIiEShBRmdyvAOoQ2pLYlJbEpPakpjqtS06BikiIhKBepAiIiIRKEGKiIhE0GgTpJm1MLMHzWylmeWb2Vdmdnw1dS8wsy/MbKuZ5ZrZ38ysWdjyaWZWZGYFocei2LUk6rZcaGblYbEWmNmYsOUdzexFMysMbe+cWLUj9P7RtOX+Ku0oNrP8sOVx/V5CMTxuZqtDvzuLzezindS9yszWmFmemT1kZi3ClsX1ewnFUKO2JPr+Eoqhpm1J6P0lFENN25Lw+0tYLP1DsTy+kzr1v7+4e6N8AG2Am4BeBP8onATkA70i1L0MOBxoDnQDvgCuC1s+Dbi4gbTlQmD6Trb1FPAM0BY4DMgDhiRiWyKs+zDwUKJ8L6EYhgAtQs8HAWuAkRHqHQusDdXfIxT77YnyvUTZloTeX6JsS0LvL9G0JcJ6Cbe/hMXyDvAx8Hg1y2OyvzTaHqS7F7r7Te6e5e4V7v4asAIYGaHufe7+sbuXuPsq4Ang0FjHXJ1o2rIzZtYGOB24wd0L3H068ApwXt1HHVlt2xIW+yOxiLOm3H2euxdXvgw9+kaoegHwYKj+ZuDPBH+cE+J7gZq3JdH3F4jqe6lWQ/tewiXq/gJgZmcBW4CpO6kWk/2l0SbIqsysCzAAmFeD6qMj1LvNzDaY2YzwIZh4qEFb9gvFutjMbggb/hoAlLv74rC6mQT/pcVFFN/L6cB64KMq5XH/XszsXjPbBiwEVgNvRKg2hOCzrpQJdDGzTiTQ91LDtlSVkPtLFG1J+P2lFt9LQu4vZtYemAT8bhdVY7K/KEECZpZC8F/uI+6+cBd1fwpkAHeEFV8L9CEYTpoMvGpmUf03Wldq0JaPgKFAZ4Kd5GzgmtCytgRDEeHygHb1E+3ORfO9EPxH+aiHxldCEuJ7cffLCT7Dw4EXgOII1ap+9pXP20VYVrk85t9LDdvynUTeX2rYlgaxv0T7vZC4+8ufCXqGObuoF5P9pdEnSDNrAjwGlAC/2kXd8cDtwPHu/t0Fct39c3fPd/did38EmAGcUH9RVxvfLtvi7svdfUVo+HIuwX9rPw4tLgDaV1mlPcExwJiK8nvpARwBPBpenijfSyiW8tBQT3eCY3RVVf3sK5/nR1hWuTzm3wvUqC1A4u8voVh22paGsr9AVN9LQu4vZjYCOBq4swbVY7K/NOoEaWYGPAh0AU5399Kd1D0OeAAYF9pRdsYBq7NAayCatlQRHutioJmZ9Q9bPpyaDTvXmVq05XzgE3dfvot6Mf9eImhG5OND8wg+60rDgbXuvpEE+V4iqK4tCb+/RFBtW6pIuP0lgl21JVH3lzEEk/OyzWwNcDVwupl9GaFubPaXaGf1JNMDuB/4DGi7i3pHARuB0RGWdSCYUdWS4BfzXKAQGJigbTke6BJ6Pgj4BrgxbPnTBDPA2hBMrIjHrLwatSWs/iLgokT7XgiG5c4iGPJpGoqnEDglQt3jCGYfDiaYlfc+O87Ki+v3EmVbEnp/ibItCb2/RNOWsHUSdX9pDXQNe9wBPAekRagbk/0lJg1PxAfQk+A/pCKCLnnl41wgPfQ8PVT3A6CsSr03Q8vSgFkE3fctBH/Yj0ngttxBMD26EFhOMGSUEratjsBLoeXZwDmJ2pZQ/YNDsbarsp1E+F7SgA9D778VmAtcEloWqS2/DX03W4H/EZq6nyDfS43b0gD2l2jakuj7S7S/Ywm7v0Ro202ETvOI1/6ia7GKiIhE0KiPQYqIiFRHCVJERCQCJUgREZEIlCBFREQiUIIUERGJQAlSREQkAiVIERGRCJQgRUREIlCCFBERiUAJUiRJmFkbM3vUzArMbK2ZXW9mr5nZw6HlWWZ2k5k9HqqzxsyurrKN35rZHDMrNLNVZvZfM+sQj/aIxJsSpEjy+AfBbYxOJbhg+HCC+wOG+y2wANgfuBG41cxOC1teAVxJcHPZc4BRwL/qNWqRBKVrsYokATNrC2wCznf3p0NlbYBc4GV3v9DMsoAl7n5M2Hr/BQa5+2HVbPc44GWglbtX1HMzRBKKepAiyaEvkALMrCxw90KC2zOF+zTC68GVL8zsKDN718xyzSyf4O70zQluPyTSqChBiiSHyhvb1npIyMx6Aq8TDMGeAYwELgotbr5b0Yk0QEqQIslhKVBKcMwQADNrDQytUu+gCK8XhJ5nECTCq9z9U3dfDOxdP+GKJL5m8Q5ARHafuxeY2UPAX81sA7Aa+BPBP8HhvcqDzOx6gju1jwHOJ7gZNcCSUP0rzewFguR5ZUwaIJKA1IMUSR5XAx8DrwAfAHOA2UBRWJ1/AsOAr4BbgInu/hyAu88BfkMw03U+cHFomyKNkmaxiiQpM2sBrAT+7u7/CM1i/be73xHfyEQaBg2xiiQJM9sP2IdgJms74NrQz2fiGZdIQ6UEKZJcfgsMBMqAr4HR7p4b14hEGigNsYqIiESgSToiIiIRKEGKiIhEoAQpIiISgRKkiIhIBEqQIiIiEShBioiIRPD/+mQBioMDs+sAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"bins = np.linspace(df.gpa.min(), df.gpa.max(), 10)\n",
"groups = df.groupby(np.digitize(df.gpa, bins))\n",
"prob = [result.predict([600, b, 2.5, 1.0]) for b in bins]\n",
"ax = plt.figure(figsize = (7, 5)).add_subplot()\n",
"ax.plot(bins, prob)\n",
"ax.plot(bins,groups.admit.mean(),'o')\n",
"ax.tick_params(labelsize=12)\n",
"ax.set_xlabel('gpa', fontsize = 14)\n",
"ax.set_ylabel('P[admit]', fontsize = 14)\n",
"ax.set_title('Marginal Effect of GPA', fontsize = 16);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"GRE Score:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "-"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFWCAYAAAAR2CYTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0IklEQVR4nO3deXyc1X3v8c/Pu6zNliWvsrybzdgGBGGHBtImTQkkJJSGQEhCaKAtWcq9SdqGEJK2l6QNtySUhNwkkJCNNtCSQNI2tCxmF4uNDXjf8SJZlizJiyzpd/84j+TxMLI18uiRZp7v+/WalzTnOfPMOTMj/eac5yzm7oiIiCTRsMEugIiIyGBREBQRkcRSEBQRkcRSEBQRkcRSEBQRkcRSEBQRkcRSEJTYmdm1ZubRbX6G4xemHL94kMroZnbrAJ6/+zWYeZR896a8Fum3f0vL+1dmtsnMOszs1Shtspk9bGaN0WM+k+N6fMbMPpDLc0bnfVtdjpB3pJndYGZPmdluMztoZtvM7NdmdrWZjUjJe23aa9huZmvN7O/MbEzaeS+03l97N7Nxua63xG/E0bOIDJgW4GrgS2np10THSmMv0SFnAVsG8flT1QPvy5De2P2LmZ0B/C3wDeDfCK8fwC3ABcC1wDZgQ47L9hlgCfBgrk54hLpkylsK/AY4Dfhe9JgmoJrwmv0QaAd+kfbQDxHe31Lg/cAXo9//IsPT3AS8mCG913JJ/lAQlMH0IPARM7vFo1UbzKwIuBz4JeEfd06Y2Wh3P9DX/O7+XK6eOwfa+1CeE6Kf33H3dWnpS939oYEp2oDorS6ZfAuoBS5w9+fTjv3UzE4BijI87lV3XxP9/l9mNg/4hJl92t270vK+McQ+D5JD6g6VwfRjYAZwbkra+4HhhCB4GDM73cz+1cy2mNk+M1sZdWMVpeV73MyWmNklZvaKmR0AboyOnRp1m+03s81Rt9tXzMzTznFYd6iZ3RqlzTOzR8ys1cw2mtktZjYsJd8YM7vDzJZHebab2a/M7PgcvF4ZmdnjwL3R3bVROe+N6nQhcF5KF97M6DGzzOwnZlZvZgfM7FUze3+Gcy8ys4fMbFfKa/7F6NgGwvt3Vcr5700/R9r5zjCz30WvTZuZPRa1/I5Ul1t7Odc04CPAdzMEQADc/RV3f+ZIZYq8TAiWlX3IKwVELUEZTBuBJwldok9FadcADwGtGfLXAK8S/km2ACcRuvtmA1em5Z0P3Al8FVgHNJpZJfAY8Fb0PO3AZ4GZWZT5IUIX2x3AJcBXgM1RGsBoQrfa1wjdjxWEAPycmR3v7tuzeK4eqde1UnRGLegbCcHgi8AHoufdBnwH+C7QGeUB2GZm04HngZ2E+tcDfwz80swuc/eHo+c8A3gcWBPl2wLMAxZG53o/8CiwFLg1Sqs/Qh0WAk8ArxNa+Q58AXjCzM5096W91KW3bukLCV+Yft3bc2ZhJtAM7MpwbFiG19/dvTMHzyuDzd110y3WG4f+Ac4FPg7sBsYAU4AO4F2Ef3AOXNzLOYzwJe4jQBcwIeXY41Ha4rTH/B0h8FWnpBUBO8KfwmF5Hbg15f6tUdrH0vK9BvznEeo6HBhLCNqfzfAazDzKa3VvlC/T7eaUfNdlOh/het3jaWnfJwSrCWnp/0XoJuy+/yQhwI89Qvk2APf38X3/V8L1unEpaWWEa5sPHq0uGc73+Sjfcb18NrpvwzK87sdFx8ZHn8EO4M/TztP9Gcx0Wz7Yf0e65eamlqAMtn8Bvk1oVc0AthNaa+enZzSzMuCvgQ8C04GRKYfncfi3+A3u/mraKc4EnnX3npaFu+8zs0eAj/WxvI+k3V8OnJJWziuAvyT8oy1POXRcH58j3U7gvRnSN/fzfO8mtOCa01o4/wF8I3qdO4BzgG+4+95+Pk+684Ffu3tTd4K77zGzhwnvf7asl/TPA3+fcv8nhC9Lqd5Mu//P7v7tXs73Z8ALaWn7+lRCGfIUBGVQuXuLhaH+VxO6pH7i7l1mGf+//RC4mNAF+irQBpwB3EVoSabaluHxUwhBK92OLIrcmHb/QOpzm9klhJGI9xG6ShsIrdJHM5Sxrw66e10/H5vJREJ38DW9HJ9AaDEPI7cjZCvI/L5sJ7TIstX9JaAGWJmSfi/wu+j3h3t57PsJdasCPgfcaGbPu/uPMuRdlePXX4YQBUEZCn5EaGENA/4kU4ZoDtelhC7Kf0pJP7mXc2baI2wbIQCkm5RVaY/sSmCNu1/bnWBmIwkBYKjYRbgGe3svx98idON2AdNy+LyNwOQM6ZN5+5eLvnicUMY/InTlAuDhuut2ADNr7+Wxyz0aHWpm/w0sI7SCf+nubf0oi+QpjQ6VoeC/gAcIQ+JX9JJnNOEf88G09GuzeJ7ngLPMrLo7IRpZmqmrsb/GEroSU11NKPtQ8VvC4JYV7l6X4XYg6gJdQpjCkmmKQbcDZJ6CkMkTwHujuX1Azzy/S6JjWXH3rYSuzj81s3dk+/iU8xwA/hfhC9KNR8kuBUYtQRl0HkbZZWwBpuRpNrPngL80s22EbsaPk11L5ZvADcB/mNlXCP/APxf9zNXu0r8FLjOzOwijFk8jTLZuOoZzjjKzMzOk73X3Zf043y2Ea1xPmtm3CYNbxgMLgNnu/vEo382E4PSsmf0joftwNmHAUfek8tcJUzD+iND6anD3Db0871cJrbbHzOx2wmv+ecIXh9v6UQ+APydcD/4fM/seoRt0N6HlfT6hlXnUSe3u/rCZvQjcbGbfdvfUa34nmFmm0cqvqdWY/xQEJZ/8CXA34RrgPkLr8dP0cYi8uzeY2UWEqRM/InQLfocwN6y362PZ+h5h0M7HgT8lrDRyCWFqRX9VAc9mSF9BCFxZcfdNZlZLGPH6d9H5dxGul96Xku9FMzuHEKC+RWiNb+TQdBAIUxm+R3gviqLHX9vL8y4zswsJq8HcRxjY8hxhovvSbOsRnXOPmV0AfBL4MPBRoJjwJekl4BPAz/t4ur8hDA76FGEKTLc7e8l/OqBrhXnO3HP1BVgk/5jZcMJE6QZ3v2iwyyMi8VJLUBLFzL5KmPy9kTAK8jrC9bE/HMxyicjgUBCUpHHCNbGp0e/LgMvc/TeDWioRGRTqDhURkcTSFAkREUksBUEREUmsgromWFlZ6TNnzhzsYoiIyBDy0ksvNbh7VaZjBRUEZ86cSV2dpu2IiMghZraxt2PqDhURkcRSEBQRkcRSEBQRkcRSEBQRkcRSEBQRkcRSEBQRkcRSEBQRkcSKNQiaWYWZPWRmbWa20cw+3Es+M7OvmdlWM2s2s8fN7KQ4yyoiIr1Y9gDcsQBuHRd+Lnsgt/ljFHdL8C6gHZgEXAXc3Utw+xBhU9LzCDtEPwv8OK5CiohIL5Y9AL+6CZo3Ax5+/uqm3gNbtvljFlsQNLNi4HLgS+7e6u5LgIeBqzNknwUscfd17t4J3A+cGFdZRUSkF4/dBgf3HZ52cF9Iz0X+mMXZEpwPdLr7qpS0pUCmluDPgblmNt/MRgIfBX6b6aRmdr2Z1ZlZXX19fc4LLSIiKZq3DGx6zOIMgiVAc1paM1CaIe824ClgJbCP0D362Uwndfd73L3W3WurqjKujyoiIrlSXj2w6TGLMwi2AmVpaWVAS4a8XwZOB6YDY4CvAP9tZmMHtIQiInJkF90CI4sOTxtZFNJzkT9mcQbBVcAIM5uXkrYIWJEh7yLgF+6+xd073P1eYDy6LigiMrgWXgGX3Anl0wELPy+5M6TnIn/MzN3jezKznwMOXAcsBh4Fznb3FWn5vgy8izCQpp4wkvQ7wDR3b+rt/LW1ta6tlEREJJWZveTutZmOxb2f4I3AD4CdwC7gBndfYWY1wOvAie6+CbgdmAi8ChQDa4DLjxQARUREshVrEHT3RuCyDOmbCANnuu/vB/4suomIiAwILZsmIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJFWsQNLMKM3vIzNrMbKOZfbiXfN8xs9aU2wEza4mzrCIiUvhGxPx8dwHtwCRgMfCImS119xWpmdz9U8Cnuu+b2b1AV3zFFBGRJIitJWhmxcDlwJfcvdXdlwAPA1f38XH3DXwpRUQkSeLsDp0PdLr7qpS0pcBJR3nc5UA98ORAFUxERJIpziBYAjSnpTUDpUd53EeBH7m7ZzpoZtebWZ2Z1dXX1+egmCIikhRxBsFWoCwtrQzodcCLmU0HLgB+1Fsed7/H3WvdvbaqqionBRURkWSIMwiuAkaY2byUtEXAil7yA1wDPOPu6wa0ZCIikkixBUF3bwMeBG4zs2IzOwe4FPjxER52DXBvDMUTEZEEinuy/I1AEbAT+Blwg7uvMLOaaD5gTXdGMzsLqAb+JeYyiohIQsQ6T9DdG4HLMqRvIgycSU17FiiOp2QiIpJEWjZNREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSK9YgaGYVZvaQmbWZ2UYz+/AR8s42s1+bWYuZNZjZ1+Msq4iIFL64W4J3Ae3AJOAq4G4zOyk9k5mNAv4L+G9gMlAN3B9jOUVEJAFiC4JmVgxcDnzJ3VvdfQnwMHB1huzXAm+5+zfdvc3d97v7srjKKiIiyRBnS3A+0Onuq1LSlgJvawkCZwIbzOw3UVfo42Z2ciylFBGRxIgzCJYAzWlpzUBphrzVwJXAncBU4BHg36Nu0sOY2fVmVmdmdfX19TkusoiIFLI4g2ArUJaWVga0ZMi7D1ji7r9x93bgH4AJwAnpGd39HnevdffaqqqqXJdZREQKWJxBcBUwwszmpaQtAlZkyLsM8FhKJSIiiRVbEHT3NuBB4DYzKzazc4BLgR9nyH4/cKaZXWxmw4HPAA3AG3GVV0RECl/cUyRuBIqAncDPgBvcfYWZ1ZhZq5nVALj7SuAjwHeA3YRg+b6oa1RERCQnRsT5ZO7eCFyWIX0TYeBMatqDhJajiIjIgNCyaSIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklgKgiIiklixBkEzqzCzh8yszcw2mtmHe8l3rZl1mllryu3COMsqIiKFb0TMz3cX0A5MAhYDj5jZUndfkSHvs+5+bpyFExHJyrIH4LHboHkLlFfDRbfAwisGu1SShdhagmZWDFwOfMndW919CfAwcHVcZRARyZllD8CvboLmzYCHn7+6KaRL3oizO3Q+0Onuq1LSlgIn9ZL/FDNrMLNVZvYlM4u71Soi0rvHboOD+w5PO7gvpEveiDOwlADNaWnNQGmGvE8CC4CNhCD5C6AD+Pv0jGZ2PXA9QE1NTQ6LKyJyBM1bskuXISnOlmArUJaWVga0pGd093Xuvt7du9z9NeA24IOZTuru97h7rbvXVlVV5bzQIiIZlVdnly5DUpxBcBUwwszmpaQtAjINiknngA1IqURE+uOiW2Bk0eFpI4tCuuSN2IKgu7cBDwK3mVmxmZ0DXAr8OD2vmb3HzCZFvx8PfAn497jKKiJyVAuvgEvuhPLpgIWfl9yp0aF55qjXBM2sItuTuntjL4duBH4A7AR2ATe4+wozqwFeB050903ARcC9ZlYC7ADuB/4u23KIiAyohVco6OW5vgyMaSB0R/aVm9l8d1/3tgMhOF6WIX0TYeBM9/2bgZuzeE4REZGs9XV06AeB3lp3qQx4tP/FERERiU9fguBG4El339WXE5rZOuDgMZVKREQkBkcNgu4+K5sTuvuC/hdHREQkPlmNDjWza8xsdIb0UWZ2Te6KJSIiMvCynSLxQ6A8Q3ppdExERCRvZBsEjcwjRWt4+5JoIiIiQ1qfRoea2WuE4OfAE2bWkXJ4ODADjQoVETk6bb80pPR1isS/Rj8XAI8Q1gHt1g5sAH6Zu2KJiBSg7u2Xunef6N5+CRQIB0mfgqC7fwXAzDYAv3D3/QNZKBGRgnSk7ZcUBAdFVlspuft9A1UQEZGCp+2Xhpy+rB26B5jt7g1m1sIRllBz9/StkkREpFt5dbQTfYZ0OSTG66Z9aQn+BYf2/PsLsltHVEREul10y+HXBEHbL6WL+bppX1aMuS/l93tzXgIRkaTo/ieu0aG9i/m6aVbXBFOZ2RjS5hm6+95jLpGISCHT9ktHFvN102yXTZthZv8eXSdsI3STpt5ERET6r7frowN03TTbluD9wBjCtcEd6PqgiIjkUszXTbMNgqcAp7v7GwNRGBERSbiYr5tmGwSXAlWAgqBIvtFyXZIvYrxumm0QvB6408zuBJaTtnmuu2/KVcFEJIe0XJdIRtkGwWHAROAhDr8e2L27xPAclUtEcknLdYlklG0QvA+oBy5BA2NE8oeW6xLJKNsgeDyw2N1X9efJzKwC+D7w+0AD8EV3/+lRHvPfwO8BI92940h5RaQXWq5L8oi7s2pHK41t7Zw1Z8KAPle2QfAFYBbQryAI3EXYemkSsBh4xMyWuvuKTJnN7Kp+lFFE0iV1uS4NBsobTXvbeWp1A0+uquep1Q1s37OfuRNL+N3nLhjQ5802wNwN/F8z+0fgNd4+MObl3h5oZsXA5cACd28FlpjZw8DVwBcy5C8HvgxcAzybZTlFJFUSl+vSYKAhraOzi6VbmnhiVQh8y7Y00eVQXjSSc+dWcsH8Ks6bXzng5TD3vl/WM7OuIxx2d+91YIyZnQI84+5FKWk3Axe4+yUZ8t8FrCEMwllPH7pDa2trva6u7ii1EJFEuGNBL13A0+Gzy+Mvj/BW0z6eXFXPk6vrWbK6gT37OxhmsGj6OC6YX8X586tYVD2O4cMsp89rZi+5e22mY9m2BGcdQzlKgOa0tGagND2jmdUC5wCfBo540cLMridM3aCmpuYYiiciQ1423ZsaDDTo9rV38tz6XTy1qoEnV9ezZmcrAJPLxvCeBVM4f34V58ydwLixowatjNluqrvxGJ6rFUjfb7CMtDVHzWwY8M/Ap929w+zI3wjc/R7gHggtwWMon4gMZdl2b2owUOzcnTe3t/S09l5cv5v2zi5GjxjGGbMquPL06Zw/v4p5E0s42v/2uPRlU91r+noyd//REQ6vAkaY2Tx3Xx2lLQLSB8WUAbXAL6IXqbuLdYuZfcjdn+preUSkgGQ71zGpg4Fi1tB6gCWrQ0vvqdUN1LccAGD+pBKuOWsG58+v4oxZFYwZOTSnkfelJXhX2v1RwEig+/rgMMIAmQNAr0HQ3dvM7EHgNjO7jjA69FLg7LSszcDUlPvTCaNSTyPMURSRJMq2ezOJg4FicKCjk5c27uap1Q08tbqe5Vv3ADBubBjQcv78Ks6fV8Xk8jGDXNK+6cumuj3X7MzsvcCtwGeA56PkdwDfBL7ah+e7EfgBsBPYBdzg7ivMrAZ4HTgxWnpte8pzdr+SOzRPUCTB+tO9qb37jpm7s7a+lSdXhaD33LpG9h3sZMQw49Sa8dz8+/M5b14VC6aV53xASxyyHRjzD8DH3T11ysLTZvYZ4F7g10d6sLs3ApdlSN9EGDiT6TEbCMuyiRQuzWc7OnVvxqaxrZ0laxp4alU9S9Y0sK15PwCzKov5UG01582r4szZFZSOGTnIJT122QbBmYTNdNPtBTQ0U6Q/NJ+tb9S9OWAOdHTy8sYmnoqu6y1/qxl3KBszgnPnVXLTvCrOnVvJ9Iqxg13UnMt2nuDj0a9XufvWKG0a8OPoXL+X8xJmQfMEJS9pPpvErHtZsu6g98L6w7s4z51XyXnzKlk4AHP2BkMu5wl+Avg3YIOZbY3SpgErydDNKSJ9oPlsEoOde/azZE0DS1Y3sGRNAzujUZzXldfxT2N+QvnwnXjZNIad/WVYeNYglzY+2c4TXGtmC4F3ERbTNsKAlt95Nk1KETlE89lkAOxt7+D59Y0h6K1uYOWOMCW7ongU58yt5Ly5lVzc+QQVj93d0xVve7Ykris+68Wpo2D3n9FNRI6VBnxIDnR2Ocu2NPH0mtDSe3ljE+2dXYwaMYwzZlbw/lOnce7cSk6cUsaw7i7OO/5P4veZzDoIRtshvZswEOawtW7c/bYclUskOTTgQ/rB3Vnf0NYT9J5Zu4uW/WEW2UlTy7j2nJmcO7fyyBPV1RWfXRA0szOBRwgT46uArcCU6P4GQEFQpD80n036oL7lAM+sDd2bT69p4K1o6sK0cUW89+QpnDO3krPnTGBCyei+nVBd8Vm3BL8B/ISwsPUe4J2EKRM/I2yWKyIiOdJ2oIMX1jf2tPbe3B6u65UXjeTsORP4s3dWcu7cSmoqxvZvLU51xWcdBBcCn3B3N7NOYLS7rzOzzwM/JQRIERHph/aOLl7dHK7rPbO2gVc2NdHR5YwaMYzaGeP53+8+jnPnVnLS1BytzqKu+KyDYHvK7zuAGcAbhB0ipmZ8hIiIZNTVFXZdeHpNA0+vDfP19rZ3YgYLp5XzyfNnc86cSmpnjh+4BagT3hWfbRB8GTidsCPE48DXzGwS8BFgWW6LJiJSWNydTY17eWbtLpasaeDZtbtobAttizlVxXzwtGrOnlPJWbMnUD42/5ckywfZBsG/5tAmuH9D2DXiW4Sg+LEclkv6S2tQSq7os5QTO/bs55m1DTyzZhfPrN3F1qZw/W1S2WguPK6Kc+ZUcs7cyrzZdaHQZDtZvi7l93rgPTkvkfSf1qCUXNFnqd+a9rbz7NoQ8J5Z28Da+rDc8rixIzlr9gQ+dcFszppTyZyq4iGzsWySZT1PUIawbDcdFemNPkt91nqggxc3NPLs2l08vaaB17ftwR3Gjhoe7aZew1lzJhw+SV2GjL7sLP8o8Cfu3tyXE5rZL4Eb3X3HsRZOsqSJr5Ir+iz1al972FT22XXhmt7SLc10djmjhg/j1Bnj+OzF8zln7gQWVo9j5PBhg11cOYq+tAT/AJhsZn0dmvQuoLj/RZJ+08RXyRV9lnoc6OjklU1NPLt2F8+u28Wrm8JyZMOHGQury0P35uxKTpsxnqJRAzSCUwZMX4Jg9yLZMtRp4qvkSoI/Swc7u1i25VDQq9uwmwMdXQwzWDCtnI+dM5Mz50zg9JkVlIzWFaV815d3sD97BG49ehbJOU18lVxJ0GfpYGcXr21t5rl1u3huXSN1G8JcPYATppRx1TtmcNacCZwxq4LyIk1bKDR92lTXzMYSlky7DBgJ/A64yd0bBrR0WdKmuiJyNB09Qa+R59btom5DI21R0DtuUilnzq7gzNkTeMfsCVQUjzrK2SQf5GJT3a8A1xKWRdsHfBi4G/hQLgooIjJQOjq7WP7Wnqilt4sX1x8KevMnlXD5adWcOTu09Cr7uvC0FIy+BsEPENYM/TmAmf0EeNrMhrt754CVTkQkS93dm8+va+T59eGaXuuBsMXQ3IklfODU6qilp6AnfQ+C04Gnuu+4+wtm1kFYLzTDELLMor0Ivw/8PtAAfNHdf5oh35WE1udkwjZNvwH+wt339PW5RCQZDnR0smxLM8+v28Xz6xt5aePunmt68yaWcNkpU3nHrAmcOXsCVaUKenK4vgbB4Ry+eDZARxaP73ZXdJ5JwGLgETNb6u4r0vI9DZzj7g1mVgJ8F/gacFOWzyciBWb/wU5e3dzU09J7edNu9h/sAuD4yaV86LRq3qHuTemjvgYxA+43swMpaWOA75nZ3u4Ed39frycwKwYuBxa4eyuwxMweBq4GvpCa193TW5edwNw+llVECkjbgQ5e3rSbF9Y38vz6Rl7d3ER7RxdmcMLkMv7kjJpwTW9mBeM1kEWy1NcgeF+GtPuzfK75QKe7r0pJWwpckCmzmZ1L2MW+DNgLvD/L5xORPNS89yAvbmjkhQ0h6C3fGlZkGT7MWDC1jI+eNYN3zArz9LTTghyrPgVBd8/FDhElQPrSa80c2pUi/TmXAOVmNg34JLAhUz4zux64HqCmpiYHxRSRONW3HOCF9Y28sD5c01u5owV3GDV8GIunj+OGC+ZwxqwKTp0xXpPTJefi/ES1Elp1qcqAliM9yN23mtlvgZ8Dp2Y4fg9wD4R5grkpqsgx0jZEGXXvp/fC+kZe3NDIixt2s74h7LIwdtRwTpsxnveePIUzZlWwaPq4gdtIViQSZxBcBYwws3nuvjpKWwSkD4rJZAQwZ8BKJpJL2oaoR2eX8+b2Pby4vpEXN+7mxfWN7GwJQwvKi0Zy+szxXHn6dN4xewInTS3TgtMSu9iCoLu3mdmDwG1mdh1hdOilwNnpec3sKsKUjM1ADfC3wGNxlVXkmCR4G6Lu6QrdLb2XNu6mZX+Yoze1fAxnRWtunjGrgrlVJdpaSAZd3B3sNwI/AHYCu4Ab3H2FmdUQFuk+0d03AScCtwPjgd3Ao8AXYy6rSP8kaBui3W3tvLRxN3Ubd1O3oZFlW5pp7wzTFeZNLOGSRVM5Y2YFp8+qYNq4okEurcjbxRoE3b2RsP5oevomwsCZ7vt/Dfx1fCUTyaEC3YbI3dncuI8XNzRSt7GRug27Wb2zFYCRw42Tox0WTpsxntqZFVp3c7DoenRWNNRKJNcKZBuig51dvLFtD3UbdlO3MQxiqY+u55WNGcFpM8Zz2SnTOH1mBQuryzWIZSjQ9eisKQiK5FqebkPUtLedVzY1UbcxXMtburmZfQfD8mPV44s4d27YOPb0mRXMm6jreUNSgq9H95eCoMhAWHjFkP6n4+6sb2jjpY27e27dXZvDhxknTS3jj0+fTu3M8Zw2YzxTynU9Ly8k6Hp0rigIiiTA/oNh1ObLm3ZTt2E3L2/aTWNbWA64u2vz0sVTOW1GBYumlzN2lP415KUCvR49kPRJFykw7s6W3ft4edNuXtnUxMubdvP6W3vo6AprScyqLOadx08MA1hmjGeOpioUjgK5Hh0nBUGRPLf/YCfLt4ZW3ssbQ9DrnpBeNHI4C6vL+eT5szm1Zjyn1IzTzgqFLE+vRw8mBUHJjoZfD6ruVt4rm5t4JWrprXirmYOdoZU3vaKIs+dM4NQZ4zm1ZjzHTy5lRFyrsOizMTQM8evRQ42CoPSdhl/Hru1AB8u2NPPK5hDwXtnURENraOWNGTmMhdPG8fFzZ3FqTQh6g7ZprD4bkqcUBKXvNPx6QHV1Oesa2kILb3MIeCu37yG6lMesymLOn1fJKTXjOKVmPMdNLh06a23qsyF5SkFQ+k7Dr3OqofUAr25qYumWJl7d3MTSzU3sidbZLB09gsU143jX783llJrxLJ4+bmhvGKvPhuQpBUHpOw2/7rfuwSuvbm7ilSjgbdkdWk7DDI6bXMZ7F05l8fRyTq3JwxGbhfLZ0HXNxFEQlL7T8Os+6exy1ta3snRzaOG9urmJldtbeqYoTBtXxKLp5Vxz1gwWTx/Pgmll+T8vrxA+G7qumUh5/pcnsdLw67fpHq25dEsTy7Y0s3RzE8u3NtPWHpYbKx09goXTy/nTC2azqHoci6ePY2LZmEEu9QAohM+GrmsmkoKgZCfhw68bWg+wbEsTSzc39wS+7pVXRg0fxglTy7j8tGoWVY9j0fRyZlfmWbfmscj3z8ZQva6pLtoBpSAo0ovmvQd5bWszy7Y28dqWZpZtaWZrU2gpmIX98t55/EQWTR/Houpyjp9cxqgRQ2S0pmRvKF7XVBftgFMQFAFaD3SwfGtzCHZbm1m2pYmNu/b2HK+pGMvi6eP46NkzWFg9jgXTyikZrT+fgjIUr2uqi3bA6a9YEqftQAevb9vDa1uaQ0tvSxPrGtrwaD7etHFFnDytnCtqp7OwupyTp5UzbuwQnp4guTEUr2sO1S7aAqIgKAWt9UAHK7aGYLfirT28trWZtfWtPQGvqnQ0i6rLed+iaSHgVZdrbc0kG2rXNYdiF22BURCUgtGy/yDLt+5h+dZmlr8VAt/6lBbepLLRLJhazntPnsLJ00LAm1SIIzWlcAzFLtoCoyAoeWlX6wFWvLWHFW/tYflbzbz+1h7WN7T1HJ9SPoaTppZz6aJpnFxdxoJp5UwsVcCTPDMUu2gLjIKgDGnuzrbm/SHYRV2aK95qZlvz/p481eOLWDC1nA+cMo0F1eUsmFo+eAtJS99p6H/fDLUu2gITaxA0swrg+8DvAw3AF939pxnyfRS4CZgH7AF+CvyVu3fEWFyJWWeXs76hjde3hUD3etTS656HN8xgdlUJZ8yqYMHUck6aWsaJU8s0aCUfaei/DBFxtwTvAtqBScBi4BEzW+ruK9LyjQU+AzwPVAEPAzcD/ye2ksqA2tvewZvbW3j9rT28vm0Pr7+1hze372H/wS4ARg435k8q5V0nTGLBtDJOnFrOCVNKB295MbVacktD/2WIiO0/ipkVA5cDC9y9FVhiZg8DVwNfSM3r7nen3N1qZj8Bfi+uskruuDv1LQd4fdse3tjW0tPKSx2wUjZmBCdOLePDZ8zoad3NqSoZOhPP1WrJPQ39lyEizq/V84FOd1+VkrYUuKAPjz0fSG8tyhDT3tHF2vpW3ti2J7q18Ma2PeyKujMhzME7cWoZ71s0lROnhIA3bVwRZkN4aTG1WnJPQ/9liIgzCJYAzWlpzUDpkR5kZh8DaoHrejl+PXA9QE1NzbGXUvpkV+uBniD3xrbQpbm2vpWDnaF5N2rEMOZPCsuKnTCljBOmlHHilDLKx44c5JL3g1otuaeh/zJExBkEW4GytLQyoKW3B5jZZYTrgBe7e0OmPO5+D3APQG1treekpNJjX3snq3e28Ob2FlZGtze3t9DQeqAnz8TS0ZwwpYwLj5vICVNKOXFKGbMqixkxVHY9P1ZqteSehv7LEBFnEFwFjDCzee6+OkpbRC/dnGb2buB7wHvd/bWYyphYnV3Oxl1tPUFu5fYWVu5oYcOuQ9fuRo8YxvxJpVx4XBXHTy7l+MllnDCllAmDvcLKQA9aUatlYGjovwwBsQVBd28zsweB28zsOsLo0EuBs9Pzmtk7gZ8A73f3F+IqYxJ0D1R5c3sLq3YcCnird7b0jMw0g5kTijluUinvWzSV4yeXctzkUmZMKGb4UNsWKI5BK2q1iBQsc4+vBzGaJ/gD4F3ALuAL7v5TM6sBXgdOdPdNZvY/wHnA/pSHP+Xu7znS+Wtra72urm6ASp9/2g50sHJHajfmHlZub2H33oM9eSpLRvcEueMml3L85FLmTSylaNTwQSx5Fu5Y0EtX5XT47PL4yyMiQ46ZveTutZmOxTrpyt0bgcsypG8iDJzpvq/pEFlo2X+QNTtbw62+lbU7W1m5o4XNjYe678aOGs78SaX8wUmTewLecZOGQFfmsdKgFRE5Blo2LU+4Ow2t7T2Bbs2OlvBzZys79hwapDJyuDGrspiF1eO44rTpUeuujOrxRYW5w7kGrYjIMVAQHGK6upytTfsOtey6g97OVpr3HerGLB41nLkTSzhnbiVzJ5Ywt6qEuRNLqKkYWzijMvtCg1ZE5BgoCA6S9o4uNu5qe1ugW1vf2jNABWBC8SjmTCzhvQun9AS6uRNLmFI+ZmhPMI+LBq2IyDFQEBxge9s7WLuzjTX1LT0Bb/XOVjbt2ktH16FBSdPGFTFnYgnvmDWBeZNKelp344u1OPRRaai9iPSTgmCO7G5r72nNrd5xaIDK1qZD3XTDhxkzJ4xl3sQS3rNgchToSpldVUzxaL0VIiJx03/eLHTvbfe2LsydrYetjzlm5DDmVJVQO3M8V1ZN7+nCnDGheOgsCi0iIgqCmXR0drGpce/bAt3a+jZaDxza0rC8aCRzJ5Zw8QmTmDephDlRF+a0cQU6ElNEpMAoCKa59ocv8MyaXbR3HhqcMrlsDHMnlvDB06p7At3ciSVUlozS4BQRkTymIJjmlOnjOW5yaU+gmzOxhLIxebjzgYiIHJWCYJpPXzxvsIsgIiIx0SgNERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJLAVBERFJrFiDoJlVmNlDZtZmZhvN7MO95FtgZv9hZg1m5pnyiIiIHKu4W4J3Ae3AJOAq4G4zOylDvoPAA8AnYiybiIgkTGy7SJhZMXA5sMDdW4ElZvYwcDXwhdS87r4SWGlmc+Mqn0heWfYAPHYbNG+B8mq46BZYeMVgl0ok78S5ldJ8oNPdV6WkLQUuiLEMIvlv2QPwq5vg4L5wv3lzuA8KhCJZirM7tARoTktrBkqP5aRmdr2Z1ZlZXX19/bGcSiQ/PHbboQDY7eC+kC4iWYkzCLYCZWlpZUDLsZzU3e9x91p3r62qqjqWU4nkh+Yt2aWLSK/iDIKrgBFmlrp1+yJgRYxlEMl/5dXZpYtIr2ILgu7eBjwI3GZmxWZ2DnAp8OP0vBaMAUZF98eY2ei4yioypF10C4wsOjxtZFFIF5GsxD1F4kagCNgJ/Ay4wd1XmFmNmbWaWU2Ubwawj0OtxH3AypjLKjI0LbwCLrkTyqcDFn5ecqcGxYj0g7kXzlz02tpar6urG+xiiIjIEGJmL7l7baZjWjZNREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSS0FQREQSK9YgaGYVZvaQmbWZ2UYz+/AR8n7WzLabWbOZ/cDMRg94AZc9AHcsgFvHhZ/LHhjwp0yEbF/XgX4f9D6LSCTuluBdQDswCbgKuNvMTkrPZGZ/AHwBuAiYCcwGvjKgJVv2APzqJmjeDHj4+aub9A/yWGX7ug70+6D3WURSxBYEzawYuBz4kru3uvsS4GHg6gzZPwp8391XuPtu4KvAtQNawMdug4P7Dk87uC+kS/9l+7oO9Pug91lEUsTZEpwPdLr7qpS0pcDbWoJR2tK0fJPMbEJ6RjO73szqzKyuvr6+/6Vr3pJduvRNtq/rQL8Pep9FJEWcQbAEaE5LawZK+5C3+/e35XX3e9y91t1rq6qq+l+68urs0qVvsn1dB/p90PssIiniDIKtQFlaWhnQ0oe83b9nypsbF90CI4sOTxtZFNKl/7J9XQf6fdD7LCIp4gyCq4ARZjYvJW0RsCJD3hXRsdR8O9x914CVbuEVcMmdUD4dsPDzkjtDuvRftq/rQL8Pep9FJIW5e3xPZvZzwIHrgMXAo8DZ7r4iLd+7gXuBdwLbgF8CL7j7F450/traWq+rq8t9wUVEJG+Z2UvuXpvpWNxTJG4EioCdwM+AG9x9hZnVmFmrmdUAuPtvga8D/wNsjG5fjrmsIiJS4EbE+WTu3ghcliF9E2EwTGraN4FvxlMyERFJIi2bJiIiiaUgKCIiiaUgKCIiiaUgKCIiiaUgKCIiiaUgKCIiiRXrZPmBZmb1hDmF6SqBhpiLM9hU52RIYp0hmfVWnftvhrtnXFy6oIJgb8ysrrfVAgqV6pwMSawzJLPeqvPAUHeoiIgkloKgiIgkVlKC4D2DXYBBoDonQxLrDMmst+o8ABJxTVBERCSTpLQERURE3kZBUEREEitvg6CZjTaz75vZRjNrMbNXzOw9KccvMrM3zWyvmf2Pmc1IOWZmdruZ7YpuXzczG5yaZMfM7jezbWa2x8xWmdl1KccKss7dzGyeme03s/tT0gq2zmb2eFTf1ui2MuVYIdf7SjN7w8zazGytmZ0XpRdcnVPe2+5bp5l9K+V4wdUZwMxmmtmjZrbbzLab2bfNbER0LN46u3te3oBi4FZgJiGY/xHQEt2vBJqBDwFjgG8Az6U89k+BlUA1MA14HfjUYNepj/U+CRgd/X48sB04rZDrnFKH/wSeAu6P7hd0nYHHgesypBdsvYF3ERa8ODP6u54W3Qq2zil1KAZagfMT8D4/Ctwb1Wsy8Bpw02DUedBfjBy/sMuAy4HrgWfSPlz7gOOj+88A16cc/0TqC50vN+A4YBtwRaHXGbgSeIDwxac7CBZ6nXsLggVb76jsn0hSnVPK/FFgHYcGLBZsnYE3gD9Muf8N4LuDUee87Q5NZ2aTgPnACkJraWn3MXdvA9ZG6aQfj34/iTxhZv9sZnuBNwlB8FEKuM5mVgbcBvxl2qGCrXOKvzezBjN72swujNIKst5mNhyoBarMbI2ZbYm6yYoo0Dqn+SjwI4/+u1PYdf4n4EozG2tm04D3AL9lEOpcEEHQzEYCPwHuc/c3gRJCkzpVM1Aa/Z5+vBkoyZf+dHe/kVCX84AHgQMUdp2/Cnzf3TenpRdynQE+D8wmdPvcA/zKzOZQuPWeBIwEPkj4bC8GTgH+hsKtMwBmVgNcANyXklzIdX6CELz2AFuAOuDfGIQ6530QNLNhwI+BduDPo+RWoCwtaxnhmmGm42VAa8o3sCHP3TvdfQmhb/wGCrTOZrYYuBi4I8PhgqxzN3d/3t1b3P2Au98HPA38IYVb733Rz2+5+zZ3bwC+SWHXuds1wBJ3X5+SVpB1jv5n/wfhC3wx4TrgeOB2BqHOeR0Eo+j/fcI3yMvd/WB0aAWwKCVfMTAnSn/b8ej3FeSnERyqWyHW+ULCYKdNZrYduBm43MxepnDr3BsHjAKtt7vvJrQKMv1DK8g6p7iGw1uBULh1rgCmA9+OvuDtAn5I+LITf50H+wLpMV5c/Q7wHFCSll5FaCZfThhhdDuHjzD6FOHC7DRgavQiDvlRVcBEwgCREmA48AdAG3BpAdd5LGH0WPftH4B/jepbkHWOyj4uen/HEL7oXBW918cVeL1vA16MPuvjCaOBv1rgdT47em9L09ILuc7rgC9En+1xwEOES1qx13nQX4xjeBFnEL4x7ic0kbtvV0XHLyYMHNlHGGU3M+WxBnwdaIxuXycakTWUb9EH5AmgidCX/hrwyZTjBVfnDK/BrUSjQwu5ztF7/SKhG6iJ8GXvXQmo90jgn6M6bwfuBMYUeJ2/C/y4l2OFWufFUX12E/YL/Bdg4mDUWWuHiohIYuX1NUEREZFjoSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoIiKJpSAoUuDMbNRgl0FkqFIQFMkzZlZsZj8ys1Yz22FmXzSzX5vZvdHxDWZ2q5n9wMyaCGsyYmZnm9kTZrbXzLaa2d3RXo0iiaUgKJJ//pGw99z7gXcSVtI/Ly3P5wjrL9YCf2VmJwP/CTwc5f8AYf3GH8RTZJGhSWuHiuQRMyshLBx8jbv/PEorJmxB9O/ufq2ZbQBec/dLUh73I+Cgu38iJW0x8Aowyd13xlcLkaFjxGAXQESyMoew08IL3Qnu3mZmy9Py1aXdPw2Ya2Z/nJLWvRv3HEBBUBJJQVAkv3QHrqN14bSl3R8G/D/gjgx5tx5roUTylYKgSH5ZAxwEzgDWA5jZWGABsPYIj3sZOMnd1wx4CUXyiAbGiOQRd28lDGa53cwuMrMTCS28YRy5dXg7cIaZfcfMTjGzuWb2R2b23RiKLTJkqSUokn9uBooJIz1bCV2ck4D9vT3A3ZeZ2fnA14AngOHAOuChAS+tyBCm0aEiec7MRgMbgW+4+z8OdnlE8olagiJ5xsxOAU4gjBAtBT4f/fzFYJZLJB8pCIrkp88BxwEdwKvA+e6+ZVBLJJKH1B0qIiKJpdGhIiKSWAqCIiKSWAqCIiKSWAqCIiKSWAqCIiKSWAqCIiKSWP8fft8JGtClGNEAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"prob = [result.predict([b, 3.4, 2.5, 1.0]) for b in sorted(df.gre.unique())]\n",
"ax = plt.figure(figsize = (7, 5)).add_subplot()\n",
"ax.plot(sorted(df.gre.unique()), prob)\n",
"ax.plot(df.groupby('gre').mean()['admit'],'o')\n",
"ax.tick_params(labelsize=12)\n",
"ax.set_xlabel('gre', fontsize = 14)\n",
"ax.set_ylabel('P[admit]', fontsize = 14)\n",
"ax.set_title('Marginal Effect of GRE', fontsize = 16);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Institution Rank:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "-"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAFWCAYAAADdQRz4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA69klEQVR4nO3deXxU9fX/8dfJAmELiwQQMKKIWkBAiIjgVm1rtSBULFURV0TBta2t3bBq/WqtVX8FkcWlWMUFEXettSgqqwYFFBFUFgERgmAgYQ/n98ed6DgmJJNJZknez8fjPjL5zGc+c+5cyJnPXc41d0dERESqJi3RAYiIiKQyJVIREZEYKJGKiIjEQIlUREQkBkqkIiIiMVAiFRERiYESqSQFM7vIzDy0HF7G8yeHPf+jBMXoZnZTDY5f+hl0qKDf5LDPInJ5NqLvH83sczPba2YLQ21tzOx5M9sces111bwe15nZWdU5Zmjc761LOf1mRnwm28xstpmdWd0xlfHebma31vT7SHLJSHQAIhG2AcOA0RHtF4SeaxL3iL51HLA2ge8frgAoKzFsLn1gZr2B/wPuBJ4l+PwAbgROAi4C1gOrqjm264BZwPTqGnA/61KexcDloce5wB+B6WbWz93nV1dcIqBEKslnOnC+md3ooWohZtYAGAw8TfDHv1qYWX1331XZ/u4+r7reuxrsrkQ8Pwj9nODuKyLaF7n7MzUTWo0ob13Ksy3s85lnZnOAz4GLASVSqVbatSvJ5hHgYOD4sLafA+kEifQ7zOwYM5tmZmvNbIeZLTOz20LJN7zfTDObZWYDzOx9M9sFjAo919PM3jaznWa2JrQL8WYz84gxvrNr18xuCrV1MrOXzKzIzFab2Y1mlhbWL8vM7jGzD0N9vjSzF8zsyGr4vMpkZjOByaFfPwvFOTm0TicDJ4Tt+uwQes0hZjbFzArMbJeZLTSzn5cxdncze8bMvgr7zP8Qem4VwfYbGjb+5MgxIsbrbWb/C302xWY2IzQD3d+63BTN5+Huawlm8bkR732zmb1nZoVmtsnMXjezPhF9Sg8rnGlm94b6FZjZo2bWrIJ1axja1uvNrHs0MUvq0IxUks1q4C2C3btvh9ouAJ4BisronwssJPhDuw3oQrDr8lDgnIi+hwNjgL8CK4DNZtYSmAF8EXqf3cCvgA5RxPwM8C/gHmAAcDOwJtQGUJ9gl/StBLtSWxAk8XlmdqS7fxnFe33DzMr6/1sSmsmPAs4H/gCcFXrf9cAEYCJQEuoDsN7MDiKYqW0kWP8C4JfA02Y2yN2fD71nb2Am8Gmo31qgE9AtNNbPgZeBRcBNobaC/axDN+BN4COCvQ0O/B5408z6uPuictYlql3sZtYEOAD4LOKpdgTbbS3QKPQ+b5lZnrsvjuj7T+BF4DzgCODvBJ/jheW8Z4tQ/5ZAX3dfGU3MkkLcXYuWhC98+0f0MOASYAuQBRwI7AV+TDCTcuBH5YxhBF8Ozwf2AQeEPTcz1NYj4jW3ESTP9mFtDYANwX+P7/R14Kaw328KtV0c0e8D4L/7Wdd0oCFB4v9VGZ9Bhwo+q8mhfmUt14f1G17WeATHL2dGtD1IkPAOiGh/DVgY9vtbBF8SGu4nvlXAo5Xc7tOAr4FmYW3ZBMd6p1e0LuWMOTO0jhmh5RDgKYIvCYdWsF0ygGXAP8PaS//dPRzR/15gJ2AR/0ZuJfiCtxR4F8hJ9P8vLTW7aNeuJKOnCGZxA4ChwJcEs8bvMbNsM7vDzD4DdgF7CHYPG8FMKdwqd18Y0dYHmOvBrj8A3H0H8FIU8Ub2/ZDv70IcYmbzzexrgi8GxUBjgplNVWwEjiljeaSK4/2UYCZZaGYZpQvwKtA99Dk3BPoBU9x9exXfJ9KJwIvu/nVpg7tvBZ4nOCGqqvoR/FvYQ7D3YQAw2COOr5rZj8zsDTP7imC77CHYc1HWdonczh8Q/DttHdHeGZhD8IXjh+5e7oxcagft2pWk4+7bLLiMYxjBLtYp7r7PzMrq/i/gRwS7cxcSJKjewDiCGW249WW8/kCCxBdpQxQhb474fVf4e5vZAOBJ4GGC3b6bCGbHL5cRY2Xtcff8Kr62LK0Idm1fUM7zBxDM3NOo3jOXW1D2dvkSaB7DuIsIZrHpBLv77wCeMrOjShObmfUk2AavApeG4igBHqDs7VLWdqaMvicSfF6/cfeyDkdILaNEKsnq3wQzgDTg3LI6mFkWMJBgd+s/w9qPKmfMsu4ZuJ4giUSKnGXE4hzgU3e/qLTBzDIJkkiy+IrgmPQd5Tz/BUFS2kdwXLG6bAbalNHehu8nrmgUhX3RmG9mK4HXCXbHXxlqH0wwCz3L3feUvtDMmhPsbq6qiUBT4FEz2+vu3ztJTmoX7dqVZPUaMJXgcocl5fSpT/DHfU9E+0VRvM884Dgza1/aEDrj92dRjFGRhgR/sMMNI4g9WfyH4IShJe6eX8ayK7Q7dxbB5UkN9jPWLoLjzJXxJvCz0MlAwDcnBg0IPVct3P0NgpPChodt64YEM9BvvmCZ2SlE7Jav2tv5VQR7RZ4ws1/EOJ4kOc1IJSm5ewnlzETD+hSa2TzgN2a2nmCX6SVEN2O6GxgJvGpmNxMkgV+HflbXXe//Awwys3sIzuLsBVxDbLOeepGXaYRs9++fbVoZNwLvEJyxei/BCUPNga4EJ+hcEup3PUGCm2tmdxHs5j2U4CSuq0N9PiK4vKY/wS7aTe6+qpz3/SvQH5hhZncQfOY3ECS5W6qwHvtzI8EejBuAqwm2y3XAZDP7F8Gx0dHAuup4M3e/zsxKgMfMLM3dn6yOcSX5aEYqqe5cYAHBt//JBH+4r63si919E3AqwVnC/wbuA/5HMHsprKYY7yeoyvNL4AWC2e6AGMfPAeaWsTxWlcHc/XMgj+DY4m0EewTGE5zw83pYv3cJTuRZA4wlOMb4W7573PQPBGe+TiU4a/Wm/bzvYoKzYrcSHEN+hOAyp5M8uPSl2rj7hwTHqoeb2YHu/irBF5p+BF9wLiE4RvxpNb7nbwi+rE0xs/1+MZTUZe7V9aVbpHYws3TgPYKZ1KmJjkdEkpt27UqdZ2Z/JZiFrCY423I4wfHCMxIZl4ikBiVSkeC43I1A29DjxcAgd38loVGJSErQrl0REZEY6GQjERGRGCiRioiIxKBOHSNt2bKld+jQIdFhiIhIkliwYMEmd8+JZYw6lUg7dOhAfn51licVEZFUZmarYx1Du3ZFRERiENdEamYtzOwZMys2s9Vmdl45/S4ysxIzKwpbTg57fqaZ7Qx7blm81kFERCRcvHftjiO4FVNroAfwkpktKqco+Vx3P34/Y13l7g/UQIwiIiKVFrcZqZk1Irht0Wh3L3L3WQQ37x0WrxhERESqWzx37R4OlLj78rC2RQQ33S3L0Wa2ycyWm9loM4ucPd8een52+G5fERGReIpnIm3M9+92UQg0KaPvWwS3b2pFMIs9l+AOE6VuILh1UztgEvCCmXUs603NbISZ5ZtZfkFBQWxrICIiEiGeibQIyI5oywa2RXZ09xXuvtLd97n7BwT3JTw77Pn57r4tdLPhh4HZlFNg3N0nuXueu+fl5MR0qZCIiMj3xDORLgcyzKxTWFt3oKwTjSI5YDE8n/oWT4V7usJNzYKfi6cmOiIRESGOidTdi4HpwC1m1sjM+hHcrf6RyL5mdrqZtQ49PpLgrvXPhX5vZmanmVmWmWWY2VDgRODVeK1L3C2eCi9cA4VrAA9+vnCNkqmISBKId0GGUUADYCPwODDS3ZeYWW7oetDcUL9TgcVmVgy8TJCAbws9lwncChQAm4CrCW55VXuvJZ1xC+zZ8d22PTuCdhERSai4Xkfq7puBQWW0f05wMlLp79cD15czRgFwTA2FmJwK10bXLiIicaMSgamgafvo2kVEJG6USFPBqTdCZoPvtmU2CNpFRCShlEhTQbchMGAMND0IsODngDFBu4iIJFSduo1aSus2RIlTRCQJaUYqIiISAyVSERGRGCiRioiIxECJVEREJAZKpCIiIjFQIhUREYmBEqmIiEgMlEhFRERioEQqIiISAyVSERGRGCiRioiIxECJVEREJAZKpCIiIjFQIhUREYmBEqmIiEgMlEhFRERioEQqIiISAyVSERGRGMQ1kZpZCzN7xsyKzWy1mZ1XTr+LzKzEzIrClpOjHUdERKSmxXtGOg7YDbQGhgLjzaxLOX3nunvjsGVmFcepNZ7KX8P23XsTHYaIiISJWyI1s0bAYGC0uxe5+yzgeWBYIsZJNcu+3MYNTy/m3Pvns7l4d6LDERGRkHjOSA8HStx9eVjbIqC8meTRZrbJzJab2Wgzy6jiOLXCEW2aMP78Xny8fitnj5/Dms3bEx2SiIgQ30TaGCiMaCsEmpTR9y2gK9CKYPZ5LvDbKoyDmY0ws3wzyy8oKKhi6MnhtC5teHT4sWwq2sXg8XNYun5rokMSEanz4plIi4DsiLZsYFtkR3df4e4r3X2fu38A3AKcHe04obEmuXueu+fl5OTEtALJ4JgOLXjqir6kmTFk4lzmrfgq0SGJiNRp8Uyky4EMM+sU1tYdWFKJ1zpg1TBOrXBEmyY8PaovrZrU54KH3uE/H65PdEgiInVW3BKpuxcD04FbzKyRmfUDBgKPRPY1s9PNrHXo8ZHAaOC5aMepzdo1a8C0K/rSpW02I6e8xyPzVic6JBGROinel7+MAhoAG4HHgZHuvsTMckPXiuaG+p0KLDazYuBlgsR5W0XjxGslkkXzRvV4bHgfTjmiFaOf/ZC7X1uOuyc6LBGROsXq0h/evLw8z8/PT3QY1W5vyT7+MP0DnlqwlnN75/LXgV3ISFfRKhGRipjZAnfPi2WMjIq7SLLLSE/j72d3o1V2fca98RmbinYx9tyjycpMT3RoIiK1nqYttYSZ8dvTjuSmAZ3539INDHtwPoXb9yQ6LBGRWk+JtJa5qN8hjD33aBau+ZohE+eyvnBHokMSEanVlEhrof7d2vLwxb1Z9/UOBt83h083lnmJrYiIVAMl0lqq72EteWJEH3aXOGdPmMuC1VsSHZKISK2kRFqLdW3XlOkj+9KsQSZDH5jHjKUbEh2SiEito0Ray+Ue0JBpI/vSqVUTRjyygKn5axIdkohIraJEWge0bFyfx0f0oW/HA/jdtMWMe+NTFW4QEakmSqR1ROP6GTx44TEM7NGWO19dxs0vfMS+fUqmIiKxUkGGOqReRhr3DOlBy8b1eXDWSgqKdnH3kO7Uz1DhBhGRqlIirWPS0ozR/TvTOrs+t738MVuKdzNxWC+aZGUmOjQRkZSkXbt11IgTO3LXL7ozf+Vmzpk0j43bdiY6JBGRlKREWocN7tWeBy7MY0VBMYPHz2HVpuJEhyQiknKUSOu4Hx7RiscuO5ainXsZPH4Oi9d+neiQRERSihKpcHRuc6aN7EtWZjrnTJrH258UJDokEZGUoUQqAHTMacz0UX3JbdGQSya/y3ML1yU6JBGRlKBEKt9onZ3Fk5cfR8/c5lz7xEIeeHtFokMSEUl6SqTyHU0bZPLwJb05vWsbbn1pKbe/slRVkERE9kOJVL4nKzOde8/ryfl9cpn45gp+89Qi9pTsS3RYIiJJSQUZpEzpacZfB3alVZMs7n5tOZuLd3Pf0J40rKd/MiIi4TQjlXKZGdec2onbzzqKt5YXcO7989lcvDvRYYmIJBUlUqnQub1zmXB+Lz5ev5Wzx89hzebtiQ5JRCRpxDWRmlkLM3vGzIrNbLWZnVeJ17xuZm5mGWFtM81sp5kVhZZlNRu5/KRLGx4dfiybinYxePwclq7fmuiQRESSQrxnpOOA3UBrYCgw3sy6lNfZzIZS/nHcq9y9cWg5ovpDlUjHdGjBU1f0Jc2MIRPmMm/FV4kOSUQk4eKWSM2sETAYGO3uRe4+C3geGFZO/6bAX4DfxStGqdgRbZrw9Ki+tG6axQUPvcN/Plyf6JBERBIqnjPSw4ESd18e1rYIKG9GehswHviynOdvN7NNZjbbzE6utiilQu2aNeCpy4+ja9tsRk55j0fmrU50SCIiCRPPRNoYKIxoKwSaRHY0szygHzC2nLFuAA4F2gGTgBfMrGNZHc1shJnlm1l+QYFqyFaX5o3qMWV4H045ohWjn/2Qu19brsINIlInxTORFgHZEW3ZwLbwBjNLA+4DrnX3vWUN5O7z3X2bu+9y94eB2cAZ5fSd5O557p6Xk5MT80rItxrUS2fisF4MyWvPmBmf8MdnPmSvCjeISB0Tz6vrlwMZZtbJ3T8JtXUHlkT0ywbygCfNDCA91L7WzH7h7m+XMbYDVgMxSwUy0tO4Y3A3cprUZ9wbn7GpaBdjzz2arMz0il8sIlILxG1G6u7FwHTgFjNrZGb9gIHAIxFdC4G2QI/QUjrT7AXMN7NmZnaamWWZWUbozN4TgVfjsBpSBjPjt6cdyU0DOvO/pRsY9uB8CrfvSXRYIiJxEe/LX0YBDYCNwOPASHdfYma5oetBcz3wZekClB7Y3ODuu4FM4NZQ+ybgamCQu+ta0gS7qN8hjD33aBatKeQXE+ewvnBHokMSEalxVpdOEMnLy/P8/PxEh1Hrzfl0EyMeWUB2Vgb/vrQ3h7X63vlkIiJJwcwWuHteLGOoRKBUu76HteSJEX3YXeKcPWEuC1ZvSXRIIiI1RolUakTXdk2ZPrIvzRpkMvSBecxYuiHRIYmI1AglUqkxuQc0ZNrIvnRq1YQRjyxgav6aRIckIlLtlEilRrVsXJ/HR/Shb8cD+N20xYx741MVbhCRWkWJVGpc4/oZPHjhMQzs0ZY7X13GzS98xL59SqYiUjvEsyCD1GH1MtK4Z0gPchrX54FZKyko2sXdQ7pTP0OFG0QktSmRStykpRl/7t+ZVtn1ue3lj9lSvJuJw3rRJCsz0aGJiFSZdu1K3I04sSN3D+nOOys388uJ89i4bWeiQxIRqTIlUkmIs3q254EL81i5qZjB4+ewalNxokMSEakSJVJJmJOPaMXjI/pQvKuEwePnsHjt14kOSUQkakqkklA9DmrGtCuOIysznXMmzePtT3TPWBFJLUqkknCH5jRm+qi+5LZoyCWT3+W5hesSHZKISKUpkUpSaJ2dxdQrjqNnbnOufWIhD7y9ItEhiYhUihKpJI3srEwevqQ3p3dtw60vLeX2l5eqcIOIJD0lUkkqWZnp3HteT4b1OZiJb63g+qcWsadkX6LDEhEplwoySNJJTzNuGdiFVk3qc9dry9m8fTf3De1Jw3r65yoiyUczUklKZsbVp3bib2cdxVvLCzj3/vlsLt6d6LBERL5HiVSS2jm9c5lwfi8+Xr+Vs8fPYc3m7YkOSUTkO5RIJen9pEsbHh1+LJuKdjF4/ByWrt+a6JBERL6hRCop4ZgOLZg2si9pZgyZMJd5K75KdEgiIoASqaSQw1s34elRfWndNIsLHnqHVz5Yn+iQRESUSCW1tGvWgGlXHEfXttmMeuw9Hpm3OtEhiUgdp0QqKadZw3pMGd6HU45oxehnP+Tu/y7DXYUbRCQx4ppIzayFmT1jZsVmttrMzqvEa143MzezjLC2qMeR2qVBvXQmDuvFkLz2jHn9U/74zAfsVeEGEUmAeF/hPg7YDbQGegAvmdkid19SVmczG0rZMUY1jtROGelp3DG4G62aZHHvG5+yqWg3Y889mqzM9ESHJiJ1SNxmpGbWCBgMjHb3InefBTwPDCunf1PgL8DvYhlHajcz4/rTjuDmM7vwv6UbGPbgfAq370l0WCJSh8Rz1+7hQIm7Lw9rWwR0Kaf/bcB44MtYxjGzEWaWb2b5BQW612VtdWHfDow992gWrSnkFxPnsL5wR6JDEpE6Ip6JtDFQGNFWCDSJ7GhmeUA/YGws4wC4+yR3z3P3vJycnKiDltTRv1tbJl9yDF98vZPB983h043bEh2SiNQB8UykRUB2RFs28J2/dmaWBtwHXOvue6s6jtRNfTu25IkRfdhd4pw9YS4LVm9JdEgiUsvFM5EuBzLMrFNYW3cg8gShbCAPeNLMvgTeDbWvNbMTohhH6qiu7ZoyfWRfmjXIZOgD85ixdEOiQxKRWixuidTdi4HpwC1m1sjM+gEDgUciuhYCbQnOxu0BnBFq7wXMj2IcqcNyD2jItJF9Obx1E0Y8soCp+WsSHZKI1FLxLsgwCmgAbAQeB0a6+xIzyzWzIjPL9cCXpQtQeobQBnffvb9x4rwukuRaNq7P45f1oW/HA/jdtMWMe+NTFW4QkWpnFf1hMbMW0Q7q7purHFENysvL8/z8/ESHIXG2e+8+fjdtEc8u/IKL+nbgxv6dSUuzRIclIknAzBa4e14sY1SmIMMmIJqv8W5mh7v7iirGJFKt6mWkcfeQHrRsXJ8HZq2koGgXdw/pTv0MFW4QkdhVtrLR2UBlZpkGvFz1cERqRlqa8ef+nWmVXZ/bXv6YLcW7mTisF02yMhMdmoikuMok0tXAW+5eqRtAmtkKQKVlJCmNOLEjLRvX53fTFvPLifOYfMkxtGqSleiwRCSFVXiykbsfUtkkGurf1d11iqQkrbN6tueBC/NY9VUxg8fPYeWm4kSHJCIpLKqzds3sAjOrX0Z7PTO7oPrCEqlZJx/Riscu60PxrhLOHj+HxWu/TnRIIpKior385V9A0zLam4SeE0kZPQ5qxrQrjqNBvXTOmTSPt5arFrOIRC/aRGqUfQZvLt+vfyuS9A7Nacz0kX3JbdGQSya/y3ML1yU6JBFJMZU6a9fMPiBIoA68aWbhNXDTgYPR2bqSolplZzH1iuMY8e98rn1iIQXbdjH8hEMTHZaIpIjKXv4yLfSzK/ASQeH4UruBVcDT1ReWSHxlZ2Uy+eLe/HrqQm59aSkF23Zxw0+PVOEGEalQpRKpu98MYGargCfdfWdNBiWSCFmZ6Yw9tycHNFrCxLdWULBtF3ec3Y3M9HhX0hSRVFLZGSkA7v5wTQUikgzS04xbBnahVZP63PXacr4q3s19Q3vSqH5U/1VEpA6p8Ku2mW01s5ahx9tCv5e51Hy4IjXPzLj61E787ayjePuTAs57YD6bi3dX/EIRqZMq8zX7ar69afbVRFd3VyRlndM7lwMa1+eqx97j7PFzePiS3hzUomGiwxKRJFPh3V9qE939Raoif9VmLpn8LlmZ6Tx8SW9+cGB2okMSkWpSHXd/qfJZFGaWZWYNw5dYAhFJVnkdWjBtZF/S04whE+Yyb0WlK2aKSB0QbYnAg83sudDx0GKCXb7hi0itdHjrJjw9si+tm2ZxwYPvMGX+at0kXESAKM/aBR4FsgiOlW5Ax0ulDmnbrAHTrjiOa59YyJ+e+ZD3P/+aWwd1JStT9zUVqcuiTaRHA8e4+9KaCEYk2TVrWI+HLjqGMTM+Yczrn/DRF1uZcH4vcg/QkQ2RuiraY6SLgJyaCEQk6S2eCvd0Jf2W5vzqw7N48aT1rN2ynf5j3+aNjzcmOjoRSZBoE+kI4C9mNtDMOppZbvhSEwGKJIXFU+GFa6BwDeBQuIYuC/7M6z8poH3zhlw8+V3ufm05Jft0tEOkrok2kaYBrYBngOXAytCyKvRTpHaacQvs2fHdtj07aDn/b0wf1Zeze7VnzIxPuHjyu2xR8QaROiXaY6QPAwXAAHSykdQlhWvLbc/KTOfOs7vRM7c5Nz2/hP5jZzHh/F4c1b6sW/eKSG0T7Yz0SOAKd3/J3fPdfUH4UtGLzayFmT1jZsVmttrMziun3zlmtszMCs1so5k9bGbZYc/PNLOdZlYUWpZFuR4i0Wnafr/tZsZ5x+by1BXH4e4MnjCHJ9/9PI4BikiiRJtI3wEOieH9xhHcdq01MBQYb2Zdyug3G+jn7k2BQwlmzrdG9LnK3RuHliNiiEmkYqfeCJkNvtuW2SBoD9P9oGa8eM0J9O7Qghue/oDfP72YnXtK4hioiMRbtIl0PPD/zGy4mR1rZj3Dl/290MwaAYOB0e5e5O6zgOeBYZF93X2Nu28KayoBDosyVpHq020IDBgDTQ8CLPg5YEzQHqFFo3o8fElvrvrhYTzx7hp+MWEua7dsj3/MIhIXUdXaNbN9+3na3b3cK9PN7Ghgjrs3CGu7HjjJ3QeU0f94gpuIZwPbgZ+7+39Dz80EugAGLAP+5O4zK4pftXYl3l77aAO/fnIh6enGP885mpMO19VjIskkEbV2D9nPcmgFr20MFEa0FQJNyurs7rNCu3bbA3cSnBlc6obQ+7UDJgEvmFnHssYxsxFmlm9m+QUFBRWEKFK9fty5NS9cfTxtsrO46F/vMGbGJ+zTJTIitUpUidTdV+9vqeDlRQSzy3DZVFCj193XAf8Bnghrm+/u29x9V+hm47OBM8p5/SR3z3P3vJwczQYk/jq0bMQzo/oxqEc77n5tOcP/nU/h9j2JDktEqkmFl7+Y2QWVHczd/72fp5cDGWbWyd0/CbV1B5ZUYugMoMwZZ+lbE+zmFUlKDeqlc/eQ7hyd24y/vvgRA+6dxfjze9KlrS6REUl1FR4jNbPIGWM9IBMoPV6aBuwBdrn7fm/UaGZPECS94UAP4GWgr7svieg3FHgbWAPkAv8GvnL3s8ysGXAs8CawF/glwe7dnu6+38tgdIxUksGC1Vu4csp7bNm+m9t+fhSDe5VzaY2I1Li4HCN19yalC3AOsBg4geAuMFmhxwuBMq8JjTAKaABsBB4HRrr7klCJwaKwMoOdgTkEu4NnE5xQdFnouUyCS2EKgE0Ed6IZVFESFUkWvQ5uzovXHE/P3Ob85qlF/OmZD9i1V5fIiKSqaM/aXQpc4u5zI9qPAyYn+/WcmpFKMtlbso87/7uMiW+uoPtBzRg/tCdtmzWo+IUiUm0ScdZuB4IbekfaTrALVkQqKSM9jT+c/gMmnN+TzzYW0X/sLGZ/uqniF4pIUok2kc4HxphZu9KG0ON7gHnVGZhIXfHTrgfy3FX9OKBRPYY9OJ9xb3yqS2REUki0ifRS4ABglZmtMrNVBNd3tuLbY5giEqWOOY159sp+nHHUgdz56jIuf3QBW3fqEhmRVBDVMVIAMzPgxwQF7A34CPifRztQAugYqSQ7d+dfs1dx28tLad+8AROG9eLINvs9GV5EYlAdx0ijTqSpTIlUUsW7qzYzasp7FO3cy+1nHcWgo9tV/CIRiVp1JNJo70eKmbUAfkpwclG98Ofc/ZZYghGRwDEdWvDS1cdz1WPvc92TC3n/8y386WedqZcR7dEYEalpUSVSM+tDUEh+F5ADrAMODP2+ClAiFakmrbKzmHLZsfztlY95cNZKPlhXyH1De9GmaVaiQxORMNF+vb0TmEJQLH4ncArBzDQfuKN6QxORzPQ0RvfvzL3nHc3HX26j/9i3mfvZV4kOS0TCRJtIuwH3hk4sKgHqu/sGgrux3FTNsYlISP9ubXnuyn5kN8jk/AfnM+mtz6hL5zeIJLNoE+nusMcbgINDj4uAttUSkYiUqVPrJjx3ZT9+0rk1t738MaOmvMc2XSIjknDRJtL3gGNCj2cCt5rZhcAYghq8IlKDmmRlct/QnvzpjB/w3482MHDcbD7ZsN87EYpIDYs2kf4J+CL0+M8EhePHAs2BEdUYl4iUw8y47MRDmTL8WLbu2MPAcbN5cfEXFb9QRGqEriMVSWFfFu7kysfeY8HqLVzS7xD+cMaRZKbrEhmRykpE0XoRSSJtmmbx+GV9uKhvBx6avZLz7p/Hxq07Ex2WSJ1SYSI1s5fNrGllBzSzp82sdWxhiUhl1ctI46Yzu/DPc3rw4bqt/GzsLN5dtTnRYYnUGZWZkZ4GtDGzFpVZCOrwNqrZsEUk0sAe7Xj2yn40rp/BuZPm8eCslbpERiQOKlPZqLQwvYgkuSPaNOG5q/rxm6mL+OuLH/H+51u4Y3A3GtWPuhqoiFRSZf53/bAK466rwmtEpBpkZ2Uy8fxeTHjrM/7x6jKWfbmNCcN60TGncaJDE6mVKnXWrpk1JCgPOAjIBP4HXOPum2o0umqms3alrpn96Saufvx9du/dx51nd+P0ow5MdEgiSSWeZ+3eDFxEULD+cYLjoONjeWMRqXn9DmvJi1cfT8dWjRk55T1uf3kpe0v2JToskVqlson0LOBSdx/h7tcCPwMGmVl6zYUmItWhbbMGTL28D+f3yWXiWys4/8H5FGzbleiwRGqNyibSg4C3S39x93eAvai+rkhKqJ+Rzq2DjuKuX3Tn/c+/pv/Yt1mwekuiwxKpFSqbSNP5bsF6CBJptPczbWFmz5hZsZmtNrPzyul3jpktM7NCM9toZg+bWXa044jIdw3u1Z7po/pSPyOdcybN5eE5q3SJjEiMKpsIDXjUzML3B2UB95vZ9tIGdz+zgnHGESTk1kAP4CUzW+TuSyL6zQb6ufsmM2sMTARuBa6JchwRidClbVNeuOp4fj11IX95fgnvf76F2846iob1dImMSFVUdkb6MEGx+q/ClkeBNRFt5TKzRsBgYLS7F7n7LOB5YFhkX3dfE3FGcAlwWLTjiEjZmjbM5P4L8vjNjw/nuUVf8PNxc1i5qTjRYYmkpEp9BXX3i6vhvQ4HStx9eVjbIuCksjqb2fEEZwlnA9uBn1dlHBEpW1qacfWpneh2UDOufeJ9zhw7i7uGdOcnXdokOjSRlBLPovWNgcKItkKgSVmd3X2WuzcF2hNcw7qqKuOY2Qgzyzez/IKCgiqGLlJ7nXR4Di9cdTwdWjZixCML+Pt/PqZkn46bilRWPBNpEcHsMlw2sN+7Erv7OuA/wBNVGcfdJ7l7nrvn5eTkRB20SF1wUIuGPHXFcZzb+yDum/kZFz70Dl8V6RIZkcqIZyJdDmSYWaewtu5AZU4QygA6VsM4IlKOrMx0bj+rG3cMPop3Vm1mwNhZLFzzdaLDEkl6cUuk7l4MTAduMbNGZtYPGAg8EtnXzIaaWa4FDgb+D5gR7TgiEr1fHpPL01f0JS3NGDJhLo/OW61LZET2I9439h4FNAA2EpQaHOnuS0JJs8jMckP9OgNzCHbjzgaWAZdVNE6c1kGk1juqfXCJzHEdD+DPz37I9U8tZueekkSHJZKUKlW0vrZQ0XqR6JTsc8bM+IR/zviEzgdmM+H8XuQe0DDRYYlUm3gWrReROig9zfjVjw/nXxcdw9ot2+k/9m1e/3hDosMSSSpKpCJSoR8e2YoXrz6B9s0bcsnkfO5+bbkukREJUSIVkUrJPaAh00f1ZXDP9oyZ8QkXT36XLcWRJbhF6h4lUhGptKzMdP7xi27838+7Mu+zr+g/dhYfrI2sjyJStyiRikhUzIyhxx7MU1cch7szeMIcnnz380SHJZIwSqQiUiXdD2rGi9ecQO8OLbjh6Q+4YZoukZG6SYlURKqsRaN6PHxJb678YUeezF/DLybMZc3m7RW/UKQWUSIVkZikpxm/Pe1I7r8gj1Wbihlw7yxmLtuY6LBE4kaJVESqxY87t+b5q4+nTXYWF09+lzEzPmGfLpGROkCJVESqzSEtGzF9VF8G9WjH3a8tZ/i/8yncvifRYYnUKCVSEalWDetlcPeQ7twysAtvf1LAgHtnseQLXSIjtZcSqYhUOzPjguM68MSI49i1t4Sz7pvDtAVrEx2WSI1QIhWRGtPr4Oa8ePUJHJ3bjOufWsSfnvmAXXt1iYzULkqkIlKjcprU59FLj+Xykw5lyvzPGTJxHuu+3pHosESqjRKpiNS4jPQ0/nD6Dxg/tCefbSxiwNhZzPpkU6LDEqkWSqQiEjenH3Ugz13VjwMa1eOCh+Yz7o1PdYmMpDwlUhGJq445jXn2yn6ccdSB3PnqMkY8soDCHbpERlKXEqmIxF2j+hmMPfdobuzfmZnLNjLw3ll8/OXWRIclUiVKpCKSEGbGJccfwuMj+lC8u4RB42bz7PvrEh2WSNSUSEUkoY7p0IKXrj6ebu2acd2TC/nLcx+ye+++RIclUmlKpCKScK2ys5hy2bFcevwhPDx3NedMmsuXhTsTHZZIpSiRikhSyExPY3T/ztx73tF8/OU2+o99mzmf6RIZSX5xTaRm1sLMnjGzYjNbbWbnldPvQjNbYGZbzWytmf3dzDLCnp9pZjvNrCi0LIvfWohITerfrS3PXdmP7AaZnP/AfCa++RnuukRGkle8Z6TjgN1Aa2AoMN7MupTRryFwHdASOBY4Fbg+os9V7t44tBxRcyGLSLx1at2E567sx2ld2nD7Kx8z8tH32LZTl8hIcopbIjWzRsBgYLS7F7n7LOB5YFhkX3cf7+5vu/tud18HTAH6xStWEUm8JlmZ3De0J3864we8tnQDA8fNZvmGbYkOS+R74jkjPRwocfflYW2LgLJmpJFOBJZEtN1uZpvMbLaZnVw9IYpIMjEzLjvxUB699Fi27tjDoHGzeWHRF4kOS+Q74plIGwORNyUsBJrs70VmdjGQB/wjrPkG4FCgHTAJeMHMOpbz+hFmlm9m+QUFBVWNXUQS6LiOB/Di1SfwgwOzufrx97nlhY/YU6JLZCQ5xDORFgHZEW3ZQLn7asxsEPA34HR3/+b0PXef7+7b3H2Xuz8MzAbOKGsMd5/k7nnunpeTkxPrOohIgrRpmsXjl/Xhor4deGj2Ss66bw7vfb4l0WGJxDWRLgcyzKxTWFt3vr/LFgAz+ylwPzDA3T+oYGwHrFqiFJGkVS8jjZvO7MK483qyYetOzrpvDtc/tYiCbbsSHZrUYXFLpO5eDEwHbjGzRmbWDxgIPBLZ18xOITjBaLC7vxPxXDMzO83Msswsw8yGEhxDfbXm10JEksHPuh3I69efzOUnHcpzC9dxyj9m8sDbK7S7VxIi3pe/jAIaABuBx4GR7r7EzHJD14PmhvqNBpoCL4ddK/pK6LlM4FagANgEXA0McnddSypShzSun8EfTv8B/7nuRI4+uDm3vrSUM/75NnM+VREHiS+rSxc65+XleX5+fqLDEJFq5u689tEG/vrSR6zZvIMzjmrDn37WmXbNGiQ6NElyZrbA3fNiGUMlAkUk5ZkZP+nShtd+dRK//vHhzFi6kVPvmsnYGZ+wc09JosOTWk6JVERqjazMdK45tRMzfnMSPzyiFXe9tpyf3PMWr320QWUGpcYokYpIrdO+eUPGn9+LRy89lnoZaVz273wunvwuKwqKEh2a1EJKpCKSWhZPhXu6wk3Ngp+Lp5bb9fhOLXnl2hP4889+QP6qLZz2/97ib698TPGuvfGLV2o9JVIRSR2Lp8IL10DhGsCDny9cs99kmpmexvATDuX160/izO7tmPDmZ5xy10yeW7hOu3ulWiiRikjqmHEL7Nnx3bY9O4L2CrRqksVdQ7rz9Mi+5DSpz7VPLOSXk+axdP3WGgpW6golUhFJHYVro2svQ6+Dm/Pclcdz28+P4pMN2/jZmLe58bkP+Xr77moKUuoaJVIRSR1N20fXXo70NOO8Y3N54/qTGXrswTw6bzU//MdMHn/nc0r2aXevREeJVERSx6k3QmZEkYXMBkF7FTRrWI+/DurKi1efwGGtGvOH6R8waNxsFcOXqCiRikjq6DYEBoyBpgcBFvwcMCZoj0HnttlMvfw4/nlODxXDl6ipRKCISJiiXXsZ+/onPDRrJVkZ6Vz7o05c2LcDmemad9RGKhEoIlLNwovh9wwrhj9bxfClHEqkIiJl6JjTmMkXH8P9F+Sxc28JQx+Yz6gpC1j39Y6KXyx1ihKpiEg5zIwfd26tYviyX0qkIiIVUDF82R8lUhGRSiothj9l+LfF8C/6l4rh13VKpCIiUep32LfF8BesVjH8uk6JVESkClQMX0opkYqIxCC8GH6rJlnfFMP/6AsVw68rlEhFRKpBr4Ob8+yV/b4pht9/rIrh1xVKpCIi1SS8GP75fVQMv65QIhURqWbNGtbjloEqhl9XKJGKiNSQ8GL4G7cFxfB/M3URG7ftTHRoUo3imkjNrIWZPWNmxWa22szOK6ffhWa2wMy2mtlaM/u7mWVEO46ISKKZGQN7tGPGb07mipM68vyidZz6jzd54O0V7CnZl+jwpBrEe0Y6DtgNtAaGAuPNrEsZ/RoC1wEtgWOBU4HrqzCOiEhSaFw/g9+ffqSK4ddCcbuNmpk1ArYAXd19eajtEWCdu/++gtf+Gvihuw+IZRzdRk1EkoG787+lG7nlxSWs2byDM45qw59+1pl2zRpU/GKpVql2G7XDgZLS5BeyCKjMTPJEYElVxjGzEWaWb2b5BQUFVQhbRKR6RRbDf/3joBj+GBXDT0nxTKSNgcKItkKgyf5eZGYXA3nAP6oyjrtPcvc8d8/LycmJOmgRkZpSWgz/f78+iVOObMXdry3nx/e8qWL4KSaeibQIyI5oywa2lfcCMxsE/A043d1LDyREPY6ISDJr37wh9w0NiuHXz0hXMfwUE89EuhzIMLNOYW3d+XaX7XeY2U+B+4EB7v5BVccREUkVZRXDv/2VpRSpGH5Si9vJRgBm9gTgwHCgB/Ay0Nfdl0T0OwV4Cvi5u79V1XEi6WQjEUkVG7ft5I5XlvH0e2tpnV2fP57xA87s3hYzS3RotUqqnWwEMApoAGwEHgdGuvsSM8s1syIzyw31Gw00BV4OtReZ2SsVjRO/1RARqVllFsOfqGL4ySiuM9JE04xURFJRyT7nyXfXcOerH1O4Yw/n9zmYX//4cJo1rJfo0FJeKs5IRUQkSuUVw39svorhJwMlUhGRFBFeDL9Tqyb88ZmgGP6C1SqGn0hKpCIiKaZz22yevLzPN8XwB49XMfxEUiIVEUlBKoafPJRIRURSWMKK4S+eCvd0hZuaBT8XT63Z90tiSqQiIrVAx5zGTL74GO6/II+de0sY+sB8Rj66gLVbtlf/my2eCi9cA4VrAA9+vnBNnU2mSqQiIrVEZDH8N5Zt5Ed3v1n9xfBn3AJ7dny3bc+OoL0OUiIVEallarwYfuHa6NprOSVSEZFaqsaK4TdtH117LadEKiJSy4UXw3+vOorhn3ojZEbchDyzQdBeBymRiojUAZnpaQw/4VBmXH8SA3u0Y+KbKzj1rpk8t3Bd9Lt7uw2BAWOg6UGABT8HjAna6yDV2hURqYMWrN7CTc8v4YN1hfTu0IKbzuxC57aRt3qu/VRrV0REqqTXwc159sp+3H7WUXyycRv9x77N6Gc/5OvtuxMdWspRIhURqaPS04xze39bDH/KfBXDrwolUhGROk7F8GOjRCoiIoCK4VeVEqmIiHyjrGL4p6gY/n4pkYqIyPeUFsN/9boT6RUqhn96PIrhpyAlUhERKdehYcXwd9V0MfwUpUQqIiL7Fbdi+ClKiVRERCqlrGL4T79XNwvVh8tIdAAiIpJaSovhv7NyM0fnNkt0OAkX1xmpmbUws2fMrNjMVpvZeeX062pmr5rZJjP73lXBZjbTzHaaWVFoWVbz0YuISLjeh7QgM107NuP9CYwDdgOtgaHAeDPrUka/PcBU4NL9jHWVuzcOLUdUf6giIiIVi9uuXTNrBAwGurp7ETDLzJ4HhgG/D+/r7suAZWZ2WLziExERqYp4zkgPB0rcfXlY2yKgrBlpZdwe2vU728xOjjU4ERGRqohnIm0MFEa0FQJNqjDWDcChQDtgEvCCmXUsq6OZjTCzfDPLLygoqMJbiYiIlC+eibQIiLzZXTawLdqB3H2+u29z913u/jAwGzijnL6T3D3P3fNycnKiDlpERGR/4plIlwMZZtYprK07sKQaxnbAqmEcERGRqMQtkbp7MTAduMXMGplZP2Ag8EhkXwtkAfVCv2eZWf3Q42ZmdlqoLcPMhgInAq/Ga11ERERKxfvyl1FAA2Aj8Dgw0t2XmFlu6HrQ3FC/g4EdfDtb3QGUXiuaCdwKFACbgKuBQaEzfUVEROIqrpWN3H0zMKiM9s8JTkYq/X0V5eyqdfcC4JiaiVBERCQ6KkkhIiISA3P/XgW+WsvMtvHtLuJU1ZJgl3aqSvX4QeuQLLQOySHV1+EId6/KZZjfqGtF65e5e16ig4iFmeWn8jqkevygdUgWWofkkOrrYGb5sY6hXbsiIiIxUCIVERGJQV1LpJMSHUA1SPV1SPX4QeuQLLQOySHV1yHm+OvUyUYiIiLVra7NSEVERKqVEqmIiEgMalUiNbMWZvaMmRWb2WozO6+cfheZWUmoLGHpcnJ8oy0zrqtCt3zbZWaTK+j7KzP70swKzeyh0lrEiVbZdUjWbQBgZvXN7MHQv6FtZva+mZ2+n/5Jty2iWYdk3RZm9qiZrTezrWa23MyG76dv0m0DqPw6JOs2CGdmncxsp5k9up8+SbkdoOL4Y9kGtSqRAuOA3UBrYCgw3szKu3H4XHdvHLbMjFeQ+/EFQR3hh/bXycxOA34PnAp0ILg36801HVwlVWodQpJxG0BwffUa4CSgKTAamGpmHSI7JvG2qPQ6hCTjtrgd6ODu2cCZwK1m1iuyUxJvA6jkOoQk4zYINw54t7wnk3w7QAXxh1RpG9SaRGpmjYDBwGh3L3L3WcDzwLDERlZ57j7d3Z8Fvqqg64XAg+6+xN23AH8FLqrh8ColinVIWu5e7O43ufsqd9/n7i8CK4Gy/gAm5baIch2SUugz3VX6a2jpWEbXpNwGENU6JDUzOwf4Gpixn25Jux0qGX+V1ZpEChwOlLj78rC2RUB5M9KjzWxTaHfLaDNLpSpPXQjWrdQioLWZHZCgeKoqJbaBmbUm+PdV1r1zU2JbVLAOkKTbwszuM7PtwMfAeuDlMrol9Tao5DpA8m6DbOAW4DcVdE3K7RBF/FDFbVCbEmljoDCirRAoq4biW0BXoBXBLPZc4Lc1Gl31ilzX0scx1YuMs5TYBmaWCUwBHnb3j8vokvTbohLrkLTbwt1HEXyWJxDcz3hXGd2SehtUch2SdhsQzCwfdPc1FfRL1u1Q2firvA1qUyItArIj2rKBbZEd3X2Fu68M7fL6gODbytlxiLG6RK5r6ePvrWuySoVtYGZpBDee3w1cVU63pN4WlVmHZN8W7l4SOlTTHhhZRpek3gZQ8Tok6zYwsx7Aj4B7KtE96bZDNPHHsg1qUyJdDmSYWaewtu6UvysrnFPO/U+T1BKCdSvVHdjg7il7XJIk2wZmZsCDBCeuDXb3PeV0TdptEcU6REqqbREmg7KPLybtNihDeesQKVm2wckEJw59bmZfAtcDg83svTL6JuN2OJnKxx+p8tvA3WvNAjwBPA40AvoR7FroUka/04HWocdHAh8Cf0mC+DOALIIz/R4JPc4oo99PgS+BzkBz4HXgb4mOP8p1SMptEBbfBGAe0LiCfsm8LSq7Dkm3LQh2r51DsLswHTgNKAYGpso2iHIdkm4bhGJpCLQJW/4BTANyUmE7RBl/lbdBQjdSDXxoLYBnQ/9YPwfOC7XnEux2yA39/g9gQ6jfCoIpfGYSxH8T357ZV7rcFBl/qO+vQ+uwFfgXUD/R8UezDsm6DUKxHRyKe2co5tJlaKpsi2jWIRm3BZADvElwpuVW4APgstBzqbINKr0OybgNylmnm4BHU2k7VDb+WLaBau2KiIjEoDYdIxUREYk7JVIREZEYKJGKiIjEQIlUREQkBkqkIiIiMVAiFRERiYESqYiUy8xmmtm9iY5DJJkpkYqkIDObbGYeWvaa2edmNt7Mmic6NpG6RolUJHX9DziQoJbocGAAcF8iAxKpi5RIRVLXLnf/0t3Xuvt/gSeBnwCYWbqZPWhmK81sh5l9Yma/C90NhlCfyWb2oplda2brzGyLmf3LzBqW94ZmdqqZfW1ml9f86omkhqS4cayIxMbMDiUoGl56h5c0YB0wBCgAegOTgK8I7ghT6gSCm03/CDgImEpwJ6Xby3iPwQT1U4e7+9QaWRGRFKREKpK6fmpmRQR3FskKtf0awINbpt0Y1neVmfUkuFlxeCLdCox0973AUjN7CjiViERqZiOAO4GzQ7NfEQlRIhVJXW8BI4AGwGUE97kcU/qkmV1BcOz04FCfTGB1xBgfhZJoqS+AYyP6DAQuB05097nVuQIitYGOkYqkru3u/qm7f+Du1xDce3E0gJn9Evh/wGSC+2D2IDgRqV7EGJE3+3a+/3dhMcHu30tDNwsXkTBKpCK1x83ADWbWFjgemO/u97r7e+7+KcGMtSpWAicTnMg0SclU5LuUSEVqCXefCSwB/kxwwlBPMzvdzDqZ2WjgpBjGXgH8kOCEJiVTkTBKpCK1y93ApcDzBGfgPga8S3Ct6V2xDOzunxHMTH8KTFQyFQmYuyc6BhERkZSlGamIiEgMlEhFRERioEQqIiISAyVSERGRGCiRioiIxECJVEREJAZKpCIiIjFQIhUREYmBEqmIiEgM/j9XLe5fQTurqwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"prob = [result.predict([600, 3.4, b, 1.0]) for b in range(1,5)]\n",
"ax = plt.figure(figsize = (7, 5)).add_subplot()\n",
"ax.plot(range(1,5), prob)\n",
"ax.plot(df.groupby('rank').mean()['admit'],'o')\n",
"ax.tick_params(labelsize=12)\n",
"ax.set_xlabel('Rank', fontsize = 14)\n",
"ax.set_xlim([0.5,4.5])\n",
"ax.set_ylabel('P[admit]', fontsize = 14)\n",
"ax.set_title('Marginal Effect of Rank', fontsize = 16);"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Logistic Regression in Perspective"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"At the start of lecture I emphasized that logistic regression is concerned with estimating a __probability__ model from __discrete__ (0/1) data. \n",
"\n",
"However, it may well be the case that we want to do something with the probability that amounts to __classification.__"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"For example, we may classify data items using a rule such as \"Assign item $x_i$ to Class 1 if $p(x_i) > 0.5$\".\n",
"\n",
"For this reason, logistic regression could be considered a classification method."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In fact, that is what we did with Naive Bayes -- we used it to estimate something like a probability, and then chose the class with the maximum value to create a classifier."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Let's use our logistic regression as a classifier.\n",
"\n",
"We want to ask whether we can correctly predict whether a student gets admitted to graduate school.\n",
"\n",
"Let's separate our training and test data:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
" df[train_cols], df['admit'],\n",
" test_size=0.4, random_state=1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Now, there are some standard metrics used when evaluating a binary classifier.\n",
"\n",
"Let's say our classifier is outputting \"yes\" when it thinks the student will be admitted."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"There are four cases:\n",
"* Classifier says \"yes\", and student __is__ admitted: __True Positive.__\n",
"* Classifier says \"yes\", and student __is not__ admitted: __False Positive.__\n",
"* Classifier says \"no\", and student __is__ admitted: __False Negative.__\n",
"* Classifier says \"no\", and student __is not__ admitted: __True Negative.__"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"__Precision__ is the fraction of \"yes\"es that are correct:\n",
" $$\\mbox{Precision} = \\frac{\\mbox{True Positives}}{\\mbox{True Positives + False Positives}}$$\n",
" \n",
"__Recall__ is the fraction of admits that we say \"yes\" to:\n",
" $$\\mbox{Recall} = \\frac{\\mbox{True Positives}}{\\mbox{True Positives + False Negatives}}$$"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Precision: 0.586, Recall: 0.340\n"
]
}
],
"source": [
"def evaluate(y_train, X_train, y_test, X_test, threshold):\n",
"\n",
" # learn model on training data\n",
" logit = sm.Logit(y_train, X_train)\n",
" result = logit.fit(disp=False)\n",
" \n",
" # make probability predictions on test data\n",
" y_pred = result.predict(X_test)\n",
" \n",
" # threshold probabilities to create classifications\n",
" y_pred = y_pred > threshold\n",
" \n",
" # report metrics\n",
" precision = metrics.precision_score(y_test, y_pred)\n",
" recall = metrics.recall_score(y_test, y_pred)\n",
" return precision, recall\n",
"\n",
"precision, recall = evaluate(y_train, X_train, y_test, X_test, 0.5)\n",
"\n",
"print(f'Precision: {precision:0.3f}, Recall: {recall:0.3f}')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Now, let's get a sense of average accuracy:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"PR = []\n",
"for i in range(20):\n",
" X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
" df[train_cols], df['admit'],\n",
" test_size=0.4)\n",
" PR.append(evaluate(y_train, X_train, y_test, X_test, 0.5))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Average Precision: 0.581, Average Recall: 0.227\n"
]
}
],
"source": [
"avgPrec = np.mean([f[0] for f in PR])\n",
"avgRec = np.mean([f[1] for f in PR])\n",
"print(f'Average Precision: {avgPrec:0.3f}, Average Recall: {avgRec:0.3f}')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Sometimes we would like a single value that describes the overall performance of the classifier.\n",
"\n",
"For this, we take the harmonic mean of precision and recall, called __F1 Score__:\n",
"\n",
"$$\\mbox{F1 Score} = 2 \\;\\;\\frac{\\mbox{Precision} \\cdot \\mbox{Recall}}{\\mbox{Precision} + \\mbox{Recall}}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Using this, we can evaluate other settings for the threshold."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"import warnings\n",
"warnings.filterwarnings(\"ignore\")\n",
"def evalThresh(df, thresh):\n",
" PR = []\n",
" for i in range(20):\n",
" X_train, X_test, y_train, y_test = model_selection.train_test_split(\n",
" df[train_cols], df['admit'],\n",
" test_size=0.4)\n",
" PR.append(evaluate(y_train, X_train, y_test, X_test, thresh))\n",
" avgPrec = np.mean([f[0] for f in PR])\n",
" avgRec = np.mean([f[1] for f in PR])\n",
" return 2 * (avgPrec * avgRec) / (avgPrec + avgRec), avgPrec, avgRec\n",
"\n",
"tvals = np.linspace(0.05, 0.8, 50)\n",
"f1vals = [evalThresh(df, tval)[0] for tval in tvals]"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "fragment"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzV0lEQVR4nO3dd3xUdfb/8ddJJwkhQAKEFDpSFBBDEVGxI66CKyrYsfe27qrrfnd1Xf2t2+xl7e66Kqio2EUUG6IEpNfQQw2dJJB6fn/ciw4hZQiZuTOZ83w85sHM3M/c+547Yc7cz733c0VVMcYYE7mivA5gjDHGW1YIjDEmwlkhMMaYCGeFwBhjIpwVAmOMiXBWCIwxJsJZITBhRRwvich2EfkxyMv+WEQuDeYy3eX+RUS2iMjGBrx2lYicHIhc1ZYzVUSubOBra80oIsNEpODQ0pn6WCGIQO5/vD0iUuRza+9Oe1ZElohIlYhc5nHUmgwFTgGyVHVgoBYiIveKyKu+z6nq6ar6SqCWWUuObOA3QC9VbVdt2oU+n98e9zP7+TMNZk4T3qwQRK4zVTXZ57befX4OcD0wy8NsdekArFLVYq+DBEkHYKuqbq4+QVX/t+/zA04H1vt+pge7IBGJaYS8JgxZITD7UdUnVXUKsLe+tiJyhoj8JCK7RGStiNzrMy1BRF4Vka0iskNEZohI21rmc5eILBeR3SKyUETOrqXdFcDzwNHur977ROQyEfm2WjsVka7u/ZdF5EkR+dCd/w8i0sWnbW8RmSwi20Rkk4j8XkSGA78HzneXM8dt+3P3h4hEicgfRGS1iGwWkf+ISAt3Wkc3w6Uissbt1rmnjvXYwn19oTu/P7jzPxmYDLR3c7xc32dSi34iMldEdorIeBFJcJc7TEQKROROt9vpJXe5+z6PrSIyQURaue3r+0w7iMh37nr+TETSfN7jWSKywH3dVBHpWcu6aOZ+ZttFZCEwoIHv2RwEKwTmUBQDlwCpwBnAdSIyyp12KdACyAZaA9cCe2qZz3LgWLf9fcCrIpJRvZGqvuDO53v3V++f/Mw51p1vSyAfeABARJoDnwOfAO2BrsAUVf0EeBAY7y6nbw3zvMy9nQB0BpKBJ6q1GQocBpwE/LG2Lz/gcZz33hk4HmedjlPVz9n/l/5lfr7f6s4DhgOdgD5u7n3aAa1wtjyuBm4GRrk52gPbgSfdtvV9phcA44A2QBxwB4CIdAdeB24F0oGPgPdFJK6GrH8Curi309xlmgCzQhC53nV/ne0QkXcbMgNVnaqq81S1SlXn4vxnP96dXI7zZdFVVStVdaaq7qplPm+q6np3PuOBZUBj9v9PVNUfVbUC+B/Qz33+V8BGVf2nqu5V1d2q+oOf87wQ+JeqrlDVIuBuYEy17pX7VHWPqs7B6XI7oKCISDRwPnC3u/xVwD+BixvwPmvzmLt+twHv88v7B6gC/qSqpaq6B7gGuEdVC1S1FLgXGO2+r/o+05dUdak7nwk+yzkf+FBVJ6tqOfAPoBkwpIas5wEPqOo2VV0LPNYoa8DUyQpB5BqlqqnubVRDZiAig0TkS7dLYyfOL8R93QH/BT4F3hCR9SLyNxGJrWU+l4jI7H2FCTjcZz6NwfdomxKcX+/g/LJd3sB5tgdW+zxeDcQAvl0ltS3XVxrOr+fq88psYK6a1JWjUFV9uwE7AO/4fBaLgEqc91XfZ1rbcvZbV6paBayl5vfY3p22z+oa2phGZoXAHIrXgElAtqq2AJ4BBEBVy1X1PlXthfPL71c4XR77EZEOwHPAjUBrVU0F5u+bjx+KgUSf+bWro211a3G6IGpS37C863G+NPfJASqATQexfIAtOL+0q89r3UHOp6Gqv8+1wOk+PxJSVTVBVdf5+5nWYL91JSKCU4Rreo8b3Gn75BzMmzENY4XA7EdE4tydiQLEujsIa/s7aQ5sU9W9IjIQp49433xOEJEj3K6PXThfdpU1zCMJ58uo0H3dOJwtAn/NAXqLSD83970H8doPgHYicquIxItIcxEZ5E7bBHSs472/DtwmIp1EJJlf9ilUHMTyUdVKnG6UB9zldwBuB16t+5UB84ybpQOAiKSLyEj3vr+faXUTgDNE5CR3C+I3QCkwrZa2d4tISxHJAm469Ldk6mOFwFT3Gc4OwCHAs+7942ppez3wZxHZDfwR5z/xPu2At3C+MBYBX1HDl5uqLsTpE/8e58v3COA7f8Oq6lLgzzg7fZcB39b9iv1euxvnnIQzcbo1luHs/AV40/13q4jUdCjtizhdJV8DK3GOsmrol9ZNOFs2K3Dyv+bO3wuP4mzlfeZ+rtOBfcXRr8+0OlVdAlyEs1N8C876PlNVy2pofh9Od9BKnL/F/x7KmzH+EbswjTHGRDbbIjDGmAhnhcAYYyKcFQJjjIlwVgiMMSbChd0gU2lpadqxY0evYxhjTFiZOXPmFlVNr2la2BWCjh07kpeX53UMY4wJKyJS61na1jVkjDERzgqBMcZEOCsExhgT4awQGGNMhLNCYIwxEc4KgTHGRDgrBMYYE+HC7jwCE95UleWFRXy9dAsJsdH0zGhOj3YpNIuL9jqaMRHLCoHxi6ryZl4BXy0tpEViLK2T4mhV7dY6KZ6WSbHEx+z/pb63vJLvl2/li8Wb+XLJZgq2738NexHo1DqJnhkp9MxoTnarRFo0i93vltIsltho24A1JhCsEISgqirluW9W0CYlnlH9MnGu7OedrUWl3DVxHpMXbqJ9iwRKK6rYXlJGVS2XskiKi6ZVchytEuOIj41mztodlFZU0Sw2mmO6pnH9sK4cf1g6lZXKwg27WLRhF4s37mLeup18OG9DrTkyU5tx3bAunJebTVyMFQVjGkvYXZgmNzdXm/IQE3vKKrlt/Gw+WeBcB/yc/lncP6o3iXG11+yqKuWtmQV8k7+Fu0/vQfvUZo2W56ulhdzx5hx2lpTzu+GHcfkxnYiKEqqqlB17ytlWXMrWojK2FZexvaSc7SVlbC0qY3uJ89zuveX0yUrlxB5tGNipFQmxdXcB7d5bzqZde9m5p5yde8rZtafi5/tTl2xm1podZLVsxi0ndePsIzOJsa0EY/wiIjNVNbfGaVYIQkfh7lKu/E8ecwt2cM+InhSVVvDolGV0a5PMUxf2p2ub5ge8ZvHGXfzhnfnkrd6OCLRKjOPJC/szuHPrepe3s6Sc6GghOf7AIrO3vJK/fryYl6et4rC2zXlkTD96ZqQ0yvtsKFVl6tJC/vXZUuat20nntCRuObkbZ/ZpD8D6nXtYtrmI/E1FLNu8mw0793LdsC4M6ZLmaW5jQoEVgjCwbNNuxr08gy1FpTw65khO690OgG+WFXLrG7PZU17Jg2cfwagjMwEodovEC9+uJCUhhrtH9KR/TipX/3cmq7eWcM+Inow7pmON3UqFu0t5/ItlvPbDGiqqlNZJceS0TiSnVSIdWiXSrkUzXp62kqWbihh3TEfuHN6j3l/ywaSqfLZwEw9PXsrijbtpl5LArr3llJT9ch31tOR4RKCktII3rj6aI7JaeJjYGO9ZIQhx0/K3cM2rM4mPiebFy3Lpk5W63/SNO/dy0+uzmLFqO2MH5nBstzT+8sFC1u/cy5gB2dw5vActk+IAp2vl9glzmLxwE2cfmcmDZx/x8xE5u/eW89w3K3n+mxWUVlRxXm422a2asXZbCWu2lbB6awnrd+yhSiG9eTz/OLcvx3evcdTakFBVpXw4bwMfzt1ARmoC3do0p1vbZLqmJ9MyKY5Nu/ZyztPT2FNWyVvXDaFTWpLXkY3xjBWCEKWqvDmzgN9PnEfn9CRevGwAWS0Ta2xbUVnFPz5byjNfLQegR7vm/GXU4eR2bHVA26oq5Ykv83n486X0ykjh8bFHMnVJIU98mc+24jLOOCKD35zanc7pyQe8tryyivU79pCWHE9SDV1G4WZFYRGjn/mepPho3r5uCG2aJ3gdyRhPeFYIRGQ48CgQDTyvqn+tNn0Y8B6w0n1qoqr+ua55NpVCMHP1Nh76eAk/rtrG0K5pPHVRf1ISYut93VdLC1m7rYTzB2TXezjlF4s3ccsbs9m9twKAIV1ac+fwHvTNTm2MtxA2Zq/dwQXPTadD6yTGXzPYr/VsTFPjSSEQkWhgKXAKUADMAMaq6kKfNsOAO1T1V/7ONxQKwYxV29haVMaww9IPuu986abd/O2TJXy+aBNpyfHccnI3xvjxpd5QK7cU89SX+ZzZtz3Hdkvz/FBUr3y1tJArXp5BbseWvDxuYEjt8zAmGOoqBIHc9h8I5KvqCjfEG8BIYGGdrwphywuL+H8fLeLzRZsBSEmI4Vd923NO/yz656TW+SW7bsceHp68lImzCkiKi+GOU7tz+dBOdR4W2hg6pSXx93P7BnQZ4eD47un849y+3Dp+NrdPmM3jY/sTHRWZRdGY6gL5LZQJrPV5XAAMqqHd0SIyB1iPs3WwoHoDEbkauBogJycnAFHrtr24jEenLOPV6atJiI3mzuE96N0+hXd+WsfEWQW89sMaOqclcc5RWRzTNY1Nu/ay1t35umZbyc87Y6OihCuGduL6YV1/3rlrgmfUkZlsKSrlLx8uok3zhdx7Vm+vIxkTEgLZNXQucJqqXuk+vhgYqKo3+bRJAapUtUhERgCPqmq3uubb2F1DxaUVvD9nPQAtk+L2GzqhWVw0//1+NY9NWUZRaQVjB+Zw2yndSUuO//n1u/eW8/G8jbw1q4AfV27bb97NE2Lo0DqRDq2S6JSWxNhBOWQ24slepmHu/2AhL3y7kvtH9ubiozt6HceYoPCqa6gAyPZ5nIXzq/9nqrrL5/5HIvKUiKSp6pYA5gKco3DG563l4cnL2FJUWmfb47qnc8+InhzW7sATuponxHLegGzOG5DNmq0lLFi/k8yWzejQKokWibZTMhT9fkRPVm0p5t73F9IxLYlju4XuIbLGBEMgtwhicHYWnwSsw9lZfIFv14+ItAM2qaqKyEDgLaCD1hHqULcI9p2M9NAni1lRWMyAji353fAeZKY2Y1tx2X63HSVlHNWxVUgfS28apqi0gtFPT2Pdjj28c/2QGs/aNqYp8WSLQFUrRORG4FOcw0dfVNUFInKtO/0ZYDRwnYhUAHuAMXUVgUM1c/U2HvxoMTNXb6dLehLPXZLLyT3b/LyTtzHH6DGhLTk+hucvzWXUk99x+ct5vHvDMbSy/TYmQkXMCWVvzyzgN2/OIb15PLef0p1zj8qyAcsMs9ZsZ8yz0+mXncqrVwyyUU1Nk1XXFkHE/NWf3Kstvxt+GF/9dhhjB+ZYETAA9M9pyd9H9+HHldu45515hNsPI2MaQ/iPIeCnFs1iuX5YV69jmBA0sl8mywuLeWzKMjqlJ9nfiYk4EVMIjKnLbSd3Y+WWYv72yRKaJ8Ry8eAOXkcyJmisEBgDiAj/PLcve8oq+L935xMXLZw/IPgnLxrjBesoN8YVFxPFkxf257ju6dw1cR4TZxV4HcmYoLBCYIyP+Jhonr34KI7u3Jo73pzz81nnxjRlVgiMqSYhNprnL80lt0Mrbh0/m0/mb/Q6kjEBZYXAmBokxsXw4rgB9MlqwU2vz+KLxZu8jmRMwFghMKYWyfExvDxuID3apXDtq7OYuXq715GMCQgrBMbUoUWzWF65fCAZLRK4+j95rN1W4nUkYxqdFQJj6tEqKY4XLxtAeWUVl788g117y72OZEyjskJgjB+6pCfzzEVHsXJLMTf8bxYVlVVeRzKm0VghMMZPQ7qm8cDZh/PNsi3c9/5CG5fINBl2ZrExB+H8ATms2FLMv79aQef0JMYd08nrSMYcMisExhykO0/rwaotxdz/wUI6tE7kxB5tvY5kzCGxriFjDlJUlPDw+f3o1T6Fm177iTlrd3gdyZhDYoXAmAZIjIvhhUsH0Co5jotf+IH563Z6HcmYBrNCYEwDtU1J4PWrBpPSLJYLn7diYMKXFQJjDkFWy0Rev2owyfExXPTCDyxcv8vrSMYcNCsExhyi7FaJvHbVIJrFRnPRCz+wZONuryMZc1CsEBjTCDq0TuL1qwYTGy1c8Nx0lm6yYmDChxUCYxpJxzSnGERFOcVgbsEOryMZ4xcrBMY0os7pybx+1WCio4RRT37HAx8upKSswutYxtTJCoExjaxrm2Q+u/V4zh+QzXPfrOS0R77m66WFXscyplZWCIwJgBaJsfy/X/dh/NWDiY2O4pIXf+S28bPZWlTqdTRjDmCFwJgAGtS5NR/dfCw3n9iVD+au5+R/fcXkhXa1MxNarBAYE2AJsdHcfuphfHjzsWS1TOSG/81i+oqtXscy5mdWCIwJku5tm/PqFYPIaZ3I1f/JY5kdYmpChBUCY4KoRWIsL102gPjYaC57aQabd+31OpIxgS0EIjJcRJaISL6I3FVHuwEiUikiowOZx5hQkN0qkRcvHcD2kjIuf2UGxaV2eKnxVsAKgYhEA08CpwO9gLEi0quWdg8BnwYqizGh5oisFjx5QX8WbdjNDa/ZpS+NtwK5RTAQyFfVFapaBrwBjKyh3U3A28DmAGYxJuSc0KMN9488nKlLCvm/9+bbpS+NZwJZCDKBtT6PC9znfiYimcDZwDN1zUhErhaRPBHJKyy0E3NM03HBoBxuOKELr/+4lqe/Wu51HBOhAlkIpIbnqv/keQS4U1Ur65qRqj6rqrmqmpuent5Y+YwJCXecehjDe7fj0c+X2f4C44lAFoICINvncRawvlqbXOANEVkFjAaeEpFRAcxkTMgRES4f2onSiiqmLLYeUhN8gSwEM4BuItJJROKAMcAk3waq2klVO6pqR+At4HpVfTeAmYwJSbkdWtKmeTwfzKn+W8mYwAtYIVDVCuBGnKOBFgETVHWBiFwrItcGarnGhKOoKGHEERlMXVrI7r3lXscxESag5xGo6keq2l1Vu6jqA+5zz6jqATuHVfUyVX0rkHmMCWW/6pNBWUUVUxZZ95AJLjuz2JgQ0T+nJRktEvhg7gavo5gIY4XAmBCxr3vo66WF7Nxj3UMmeKwQGBNCzuiTQVllFZ/bUNUmiKwQGBNCjsxOJTO1GR/Os+4hEzxWCIwJISLCGX0y+GZZITtLrHvIBIcVAmNCzBlHZFBeqXy2cKPXUUyEsEJgTIjpk9WCrJbN7OghEzRWCIwJMfu6h77L38L24jKv45gIYIXAmBB0Zp/2VFRZ95AJDisExoSg3u1T6NA60bqHTFBYITAmBIkIZxyRwbTlW9laVLrftMoqZfLCTfz+nXls3m3XPDaHLsbrAMaYmp3RJ4Onpi7n0wWbuGBQDluLShmft5b/TV/Duh17AGiZGMtvT+vhcVIT7qwQGBOiemWk0Dktidd/XEPeqm18MHcDZZVVHN25NX84oyev/biGt2eu4/ZTDiM6qqbrQBnjHysExoSofUcPPf5FPisKixgzMJuLB3egW9vmAFQp3PDaLL7N38Lx3e3KfabhrBAYE8KuOb4LPdqlcFz3NJonxO437eRebUhNjGVC3lorBOaQ2M5iY0JYcnwMZ/TJOKAIAMTHRDOqXyaTF2xiR4mdb2AazgqBMWHs3NwsyiqreG+2XeLSNJwVAmPCWO/2LeiVkcKbM9d6HcWEMSsExoS583KzmL9uFwvX7/I6iglTVgiMCXMj+2USFx1lWwWmwawQGBPmWibFcXKvNrz70zrKKqq8jmPCkBUCY5qAc3Oz2V5SzpRFdolLc/CsEBjTBBzXLZ22KfG8ObPA6ygmDFkhMKYJiI4SzumfxdQlm9m0ywaiMwfHCoExTcToo7KoUpg4a53XUUyY8asQiMhQERnn3k8XkU6BjWWMOVid05MZ0LElb85ci6p6HceEkXoLgYj8CbgTuNt9KhZ4NZChjDENc+5R2awoLGbWmu1eRzFhxJ8tgrOBs4BiAFVdDzQPZChjTMOM6JNBUlw0r0xb7XUUE0b8KQRl6mxnKoCIJPk7cxEZLiJLRCRfRO6qYfpIEZkrIrNFJE9Ehvof3RhTXXJ8DBcd3YEP5q5n5ZZir+OYMOFPIZggIv8GUkXkKuBz4Ln6XiQi0cCTwOlAL2CsiPSq1mwK0FdV+wGXA88fRHZjTA2uHNqZ2Ogonvoy3+soJkzUWQhERIDxwFvA28BhwB9V9XE/5j0QyFfVFapaBrwBjPRtoKpF+sterSTcrQ5jTMOlN49n7MAc3vlpHWu3lXgdx4SBOguB+yX9rqpOVtXfquodqjrZz3lnAr6DnxS4z+1HRM4WkcXAhzhbBQcQkavdrqO8wsJCPxdvTOS65vjOiMC/v17udRQTBvzpGpouIgMaMO+aLqJ6wC9+VX1HVXsAo4D7a5qRqj6rqrmqmpuebldiMqY+GS2aMfqobCbMKLATzEy9/CkEJ+AUg+Xujt15IjLXj9cVANk+j7OAWq+eoapfA11EJM2PeRtj6nHd8V2oVOXZr1d4HcWEOH+uWXx6A+c9A+jmnny2DhgDXODbQES6AstVVUWkPxAHbG3g8owxPnJaJzKyX3v+98Nqrh/WhdbJ8V5HMiGq3i0CVV0NpAJnurdU97n6XlcB3Ah8CiwCJqjqAhG5VkSudZudA8wXkdk4Rxidr3ZKpDGN5vphXSmtqOKFb1d6HcWEMKnve1dEbgGuAia6T50NPOvnkUONLjc3V/Py8rxYtDFh6YbXZvHVkkK+u/NEWiTGeh3HeEREZqpqbk3T/NlHcAUwSFX/qKp/BAbjFAZjTBi48YSuFJVW8PK0VV5HMSHKn0IgQKXP40pqPiLIGBOCemakcHLPtrz43UqKSiu8jmNCkD+F4CXgBxG5V0TuBaYDLwQ0lTGmUd10Yld27innyS/zbWRScwB/dhb/CxgHbAO2A+NU9ZEA5zLGNKK+2amccUQGT09dzuhnvmduwQ6vI5kQ4s8w1IOBZar6mKo+CuSLyKDARzPGNKbHxx7J387pw+qtxZz1xHfc8eYcNtvJZgb/uoaeBop8Hhe7zxljwkhUlHDegGy+vGMY1xzXmfdmr+OEf0zl6anLKa2orH8Gpsnya2ex77H9qlqFfyeiGWNCUPOEWO4e0ZPPbjueo7uk8dAnixn99PdUVdm+g0jlTyFYISI3i0ise7sFsHPWjQlzndKSeP7SXP48sjfz1u1k+go7qT9S+VMIrgWG4AwTsQ4YBFwdyFDGmOA5LzeblIQYxuetrb+xaZLq7eJR1c044wQZY5qghNhoRvbLZHzeWv5cUm5nH0egWrcIROQqEenm3hcReVFEdrojkPYPXkRjTKCdPyCbsooqJs2tdYBg04TV1TV0C7DKvT8W6At0Bm4HHg1sLGNMMPVun0LPjBQmzLDuoUhUVyGoUNVy9/6vgP+o6lZV/RznspLGmCZCRDg/N4t563aycP0ur+OYIKurEFSJSIaIJAAn4Vy0fp9mgY1ljAm2kf0yiYuOYoLtNI44dRWCPwJ5ON1Dk1R1AYCIHI8dPmpMk9MyKY5Te7fl3dnr7ASzCFNrIVDVD4AOQE9V9R12Og84P9DBjDHBd15uNjtKypm8cJPXUUwQ1XkegapWqOr2as8Vq2pRba8xxoSvY7qm0b5FAhPyCryOYoLInxPKjDERIjpKGJ2bzTfLClm3Y4/XcUyQWCEwxuzn3KOyUIW3Z9pWQaRoUCEQkR6NHcQYExqyWyVyTNfWTMhbawPRRYiGbhF81qgpjDEh5bzcbAq277GB6CJErWMNichjtU0CUgOSxhgTEk7r3e7ngeiGdE3zOo4JsLq2CMYB84GZ1W55QFngoxljvLJvILqP529kS1Gp13FMgNVVCGYA81X1leo3YHeQ8hljPHLpkI6oKvd/sNDrKCbA6ioEo4HZNU1Q1U4BSWOMCRld2yRzwwldeW/2er5YbCeYNWV1FYJkVS0JWhJjTMi5flhXurdN5p535rN7b3n9LzBhqa5C8O6+OyLyduCjGGNCTVxMFH89pw8bd+3lb58s8TqOCZC6CoH43O8c6CDGmNDUP6cl44Z04r/TVzNj1Tav45gAqKsQaC33jTER5o7TupPVshl3vj2XveU2MmlTU1ch6Csiu0RkN9DHvb9LRHaLiF9XrhCR4SKyRETyReSuGqZf6F76cq6ITBORvg19I8aYwEmMi+HBs49gRWExj3+xzOs4ppHVNQx1tKqmqGpzVY1x7+97nFLfjEUkGngSOB3oBYwVkV7Vmq0EjlfVPsD9wLMNfyvGmEA6rns6o4/K4t9frbCrmDUxgRx0biCQr6orVLUMeAMY6dtAVaf5DHM9HcgKYB5jzCH6wxk9SU2M486351JRWeV1HNNIAlkIMgHfa94VuM/V5grg45omiMjVIpInInmFhYWNGNEYczBSE+O476zezFu3k79/akcRNRWBLARSw3M17nQWkRNwCsGdNU1X1WdVNVdVc9PT0xsxojHmYI04oh0XDc7h31+v4PUf13gdxzSCWgedawQFQLbP4yxgffVGItIHeB44XVVtqENjQpyIcO+ZvVm7bQ9/eHc+WS2bcWw3+4EWzgK5RTAD6CYinUQkDhgDTPJtICI5wETgYlVdGsAsxphGFBMdxRMXHEm3Nslc/+oslm6y4cfCWcAKgapWADcCnwKLgAmqukBErhWRa91mfwRaA0+JyGwRyQtUHmNM42qeEMsLlw0gIS6acS/NoHC3jVIarkQ1vM4Vy83N1bw8qxfGhIp5BTs579/f071dc964ajDN4qK9jmRqICIzVTW3pml2zWJjzCE5IqsFj47px9yCHdw+YbZd3jIMWSEwxhyyU3u3454RPfl4/kYe+dx294UbKwTGmEZxxdBOjD4qi8e/zLdrHYcZKwTGmEYhItx3Vm86tk7itvGz2VFiV7QNF1YIjDGNJik+hsfGHMmWolLuense4XYwSqSyQmCMaVRHZLXgt6cdxicLNvLGjLX1v8B4zgqBMabRXTm0M8d2S+O+9xeQv9lONgt1VgiMMY0uKkr457l9SYyL4abXZ1NaYRezCWVWCIwxAdEmJYG/j+7Dog27eOhjG6k0lFkhMMYEzEk923Lp0R148buVfLlks9dxTC2sEBhjAuruET05rG1zfj9xnl3vOERZITDGBFRCbDR/HtmbDTv38sK3K72OY2pghcAYE3CDOrfmlF5teXrqcrYU2SilocYKgTEmKO46vQd7yit59PNlXkcx1VghMMYERZf0ZC4clMNrP64hf3OR13GMDysExpigueWkbjSLjeavHy/2OorxYYXAGBM0rZPjuW5YFz5ftMlGKA0hVgiMMUF1xdBOZLRI4MGPFtlFbEKEFQJjTFAlxEbz29MOY27BTt6fu97rOAYrBMYYD4zql0nv9in87ZMldpJZCLBCYIwJuqgo4Z4RPVm3Yw8vT1vldZyIZ4XAGOOJIV3TOLFHG578It+uZuYxKwTGGM/8bvhhFJVV8MxXK7yOEtGsEBhjPNOjXQpn9W3Py9NWsnn3Xq/jRCwrBMYYT912cnfKK5Unv8j3OkrEskJgjPFUx7QkzsvN4rUf11CwvcTrOBHJCoExxnM3ndgNEbEB6TxihcAY47n2qc24aFAH3p5VwPJCG5Au2KwQGGNCwvUndCEhNpp/TV5aa5vv8rdw2sNf8/XSwiAma/oCWghEZLiILBGRfBG5q4bpPUTkexEpFZE7ApnFGBPa0pLjufyYTnw4dwML1u/cb1pVlfL4lGVc/MIPLNm0m+e+scNNG1PACoGIRANPAqcDvYCxItKrWrNtwM3APwKVwxgTPq46rjMpCTH887Nftgq2F5dx+Ssz+OfkpZzVtz1XDu3Et/lb2LBzj4dJm5ZAbhEMBPJVdYWqlgFvACN9G6jqZlWdAZQHMIcxJky0aBbLNcd34YvFm5m5ehuz1+7gV49/y7T8rfxl1OE8fH4/Lj66A6owcdY6r+M2GYEsBJnAWp/HBe5zB01ErhaRPBHJKyy0vkFjmrJxx3QkLTmO28bP4dxnpgHw1nVHc9HgDogIHVonMaBjS96eVYCqDWPdGAJZCKSG5xr0qanqs6qaq6q56enphxjLGBPKEuNiuPGErqzZVsKx3dL58Oah9MlK3a/NOf2zWFFYzOy1OzzJ2NTEBHDeBUC2z+MswAYfN8bU69IhHTkypyVHZLYgKurA35Qj+mTwp0kLeHtWAUfmtPQgYdMSyC2CGUA3EekkInHAGGBSAJdnjGkiRIS+2ak1FgGAlIRYhh/ejkmz19v1DBpBwAqBqlYANwKfAouACaq6QESuFZFrAUSknYgUALcDfxCRAhFJCVQmY0zTcU7/LHbtrWDKos1eRwl7gewaQlU/Aj6q9twzPvc34nQZGWPMQTmmaxrtUhJ4e1YBZ/TJ8DpOWLMzi40xYSk6Sji7fyZfLS20IawPkRUCY0zYOqd/FpVVyns/2XEoh8IKgTEmbHVtk0zf7FQ7p+AQWSEwxoS10f0zWbxxNwvW7/I6StiyQmCMCWtn9m1PXHQUb88q2O/5qipl6pLNXPnKDO6eOM+jdOEhoEcNGWNMoKUmxnFyrza8N3s9d5/ek70VlbyVV8B/p69m5ZZiYqOF8krlhhO6kNUy0eu4Icm2CIwxYe+c/llsKy7j6v/mMfjBKfz5g4WkJsby6Jh+fHzLcQBMmmM7lGtjWwTGmLB3XPd02qUkMG35Vs7s055Lh3TYb3yi/jmpTJq9nuuHdfUuZAizQmCMCXux0VFMuvEYYqOjaJkUd8D0kf0y+dOkBSzeuIse7eofvGBnSTktEmMDETUkWdeQMaZJaJOSUGMRADijTwbRUcK7fpxv8F3+Fo68/zMmL9zU2BFDlhUCY0yTl5Ycz7Hd0pg0ex1VVbWfb6CqPDx5KVUKD360iPLKqiCm9I4VAmNMRBjVL5P1O/cyY9W2WttMX7GNvNXbOaVXW1ZuKea1H9YEMaF3rBAYYyLCKb3a0iw2mndn19499PgXy0hvHs/jY4/k6M6teeTzpeza2/SvpGuFwBgTEZLiYzi1d1s+mreBsooDu3xmrt7GtOVbuea4ziTERnPPGT3ZXlLO01OXe5A2uKwQGGMixqh+mezcU87UJQdew+DxL/JplRTHBYNyADg8swW/PjKTF75dybode4IdNaisEBhjIsbQbmm0SorjvWrdQ3MLdjB1SSFXDO1EYtwvR9X/5rTDEOAfny4JctLgskJgjIkYsdFR/KpPBp8v2sRun77/J77IJyUhhkuO7rBf+8zUZlwxtBPv/LSOeQU7gx03aKwQGGMiysh+mZRWVPHJ/I0ALNqwi88WbmLcMZ1onnDgSWTXDutCq6Q4HvhoYZMd6toKgTEmovTPSSWnVeLP3UNPfJlPcnwMlx/Tqcb2KQmx3HpyN6av2MYXi5vm9ZGtEBhjIoqIMLJfe6Yt38L3y7fy0bwNXHJ0hzqHlBg7MIfOaUk8+NEiKprgSWZWCIwxEWdkv0yqFK7730wSYqK5YmjNWwP7xEZHcdfpPVheWMxjX+QHKWXwWCEwxkScrm2SOTwzhR0l5Vw4KIfWyfH1vuaUXm359ZGZPDZlGY9PWRaElMFjo48aYyLSmAE5PLR1MVcf19mv9iLC38/tC8A/Jy9FgZtP6hbAhMFjhcAYE5EuHJTD6KOySIiN9vs10VFuMRD41+SlVKly68ndA5gyOKwQGGMikogcVBHYJzpK+PvovgjCI58vQxVuOyW8i4EVAmOMOUjRUcLfRvchSuDRKctQ4LaTuyEiXkdrECsExhjTANFRwkPn9EEEHpuyjILtJVxzXBcOa9fc62gHzQqBMcY0UFSU8Ndf96FlYhwvTVvFxFnrOLpzay4d0pGTe7YhJjo8DsyUQJ4yLSLDgUeBaOB5Vf1rteniTh8BlACXqeqsuuaZm5ureXl5AUpsjDENs624jPEz1vLq9NWs27GHzNRmXDg4h1N6tmVveRW7S8sp2ltBUWkFxaUV7C2vIj42ioTYaJrtu8VFkxgXTfe2zUmKb9zf6SIyU1Vza5wWqEIgItHAUuAUoACYAYxV1YU+bUYAN+EUgkHAo6o6qK75WiEwxoSyisoqPl+0mVemreL7FVsbNI8ogR7tUujfIZX+OS3pn9OSDq0TD2kfRF2FIJBdQwOBfFVd4YZ4AxgJLPRpMxL4jzrVaLqIpIpIhqpuCGAuY4wJmJjoKIYf3o7hh7dj6abdLFy/i6T4GJLjY2ie4PybFB9DQmwUpRVV7CmrZG95JXvKK9lTVsnOPeXMW7eTWWu2886sdbw63blcZuukOK4b1oUrj/XvvIeDytzoc/xFJrDW53EBzq/++tpkAlYIjDFhr3vb5nRvW/vO49qmnNq7HQCVVcrSTbuZtWY7P63ZQXrz+s+AbohAFoKatmGq90P50wYRuRq4GiAnJ+fQkxljTBiIjhJ6ZqTQMyOFCwd1qP8FDRTIXdoFQLbP4yyg+lWj/WmDqj6rqrmqmpuent7oQY0xJpIFshDMALqJSCcRiQPGAJOqtZkEXCKOwcBO2z9gjDHBFbCuIVWtEJEbgU9xDh99UVUXiMi17vRngI9wjhjKxzl8dFyg8hhjjKlZQE8oU9WPcL7sfZ97xue+AjcEMoMxxpi6hcdpb8YYYwLGCoExxkQ4KwTGGBPhrBAYY0yEC+igc4EgIoXAaq9z+EgDtngdog6hng9CP2Oo5wPL2BhCPR8cWsYOqlrjiVhhVwhCjYjk1TaQUygI9XwQ+hlDPR9YxsYQ6vkgcBmta8gYYyKcFQJjjIlwVggO3bNeB6hHqOeD0M8Y6vnAMjaGUM8HAcpo+wiMMSbC2RaBMcZEOCsExhgT4awQ+EFEhovIEhHJF5G7apjeQ0S+F5FSEbkjRDNeKCJz3ds0EekbYvlGutlmi0ieiAwNZj5/Mvq0GyAilSIyOpj53GXXtx6HichOdz3OFpE/hlI+n4yzRWSBiHwVzHz+ZBSR3/qsv/nuZ90qhPK1EJH3RWSOuw4PfdRmVbVbHTecIbSXA52BOGAO0KtamzbAAOAB4I4QzTgEaOnePx34IcTyJfPLPqs+wOJQW4c+7b7AGVV3dKhlBIYBHwT7b/Ag8qXiXLc8x33cJtQyVmt/JvBFKOUDfg885N5PB7YBcYeyXNsiqN9AIF9VV6hqGfAGMNK3gapuVtUZQLkXAfEv4zRV3e4+nI5zNbhQylek7l82kEQNlyz1OqPrJuBtYHMww7n8zegVf/JdAExU1TXg/N8JwYy+xgKvByWZw598CjQXEcH5AbUNqDiUhVohqF8msNbncYH7XCg52IxXAB8HNNH+/MonImeLyGLgQ+DyIGXbp96MIpIJnA08gzf8/ZyPdrsNPhaR3sGJBviXrzvQUkSmishMEbkkaOkcfv9fEZFEYDhO4Q8Wf/I9AfTEuazvPOAWVa06lIUG9MI0TYTU8FyoHXPrd0YROQGnEASzD96vfKr6DvCOiBwH3A+cHOhgPvzJ+Ahwp6pWOj/Ggs6fjLNwxpQpEpERwLtAt0AHc/mTLwY4CjgJaAZ8LyLTVXVpoMO5Dub/85nAd6q6LYB5qvMn32nAbOBEoAswWUS+UdVdDV2obRHUrwDI9nmchVOJQ4lfGUWkD/A8MFJVtwYpGxzkOlTVr4EuIpIW6GA+/MmYC7whIquA0cBTIjIqKOkc9WZU1V2qWuTe/wiIDeJ69GcdFgCfqGqxqm4BvgaCeeDCwfwtjiG43ULgX75xON1rqqr5wEqgxyEtNZg7asLxhvMLZgXQiV923vSupe29eLOzuN6MQA7OtaGHhGi+rvyys7g/sG7f41DJWK39ywR/Z7E/67Gdz3ocCKwJ1nr0M19PYIrbNhGYDxweSuvQbdcCp+89KQQ/46eBe937bd3/K2mHslzrGqqHqlaIyI3Apzh79F9U1QUicq07/RkRaQfkASlAlYjcirOnv8Gbao2dEfgj0BrnVyxAhQZppEU/850DXCIi5cAe4Hx1/9JDKKOn/Mw4GrhORCpw1uOYYK1Hf/Kp6iIR+QSYC1QBz6vq/GDk8zej2/Rs4DNVLQ5WtoPIdz/wsojMw+lKulOdrasGsyEmjDEmwtk+AmOMiXBWCIwxJsJZITDGmAhnhcAYYyKcFQJjjIlwVgjMIRGR1j4jNW4UkXXu/R0isjAAy7v3YEd4FZGiWp5/uaYRRMUZTXa2iPwkIl0amDNWRP4qIsvcESx/FJHT3WmrGuskLxE5a98IlSKSLiI/uLmPFZGPRCS1AfMcJiJDfB5f68FQECaI7DwCc0jUOUO5Hzhf0kCRqv5DRDoCH9T3ehGJUdVDGjArAEYB76nqn/xp7A7+Jbr/eC/3Axk4J0uVikhb4PjGDqqqk4BJ7sOTcEZtvdR9/E0DZzsMKAKmucvw/BwKE1i2RWACKVpEnnPHTP9MRJoBuAOOPeiORX+LiBwlIl+5g5B9KiIZbrubRWShONcpeMNnvr3ceawQkZv3PSkit7u/vue7J/XtRxxPuPP8EGf48OptRgC3AleKyJe1zVdEOorIIhF5Cmd8n2yfeSQCVwE3qWopgKpuUtUJNSzvXfd9LxCRq93not2tlfkiMk9EbqttfYjIZe576gf8DRjhbs00893yEJFL3NfNEZH/us+d6bMF8bmItHUL+LXAbe58jvXdChORfiIy3Z3XOyLS0uczfcjd8lkqIsfW/mdhQk4wT5+2W9O+4TPEBtARZ2jcfu7jCcBF7v2pwFPu/VicX57p7uPzcc6mBGeMlXj3fqrPMqYB8UAasNWdx1E4IzEm4QzNuwA40n1Nkfvvr4HJOGdstgd2UMMwEdXeR43zdd9fFTC4htf3AX6qYz2twh0SAGjl/tsMZ7iF1u4yJ/u0T61jfVwGPFH9vu9ygN7AkhqW2ZJfTiq9Evhn9fdfw/qYCxzv3v8z8IjPZ7rv9SOAz73+e7Sb/zfrGjKBtFJVZ7v3Z+J8ee4z3v33MOBwnBEUwfmS3uBOmwv8T0TexRlFc58P1fmlXSoim3HGWxkKvKPukAAiMhE4FvjJ53XHAa+raiWwXkS+8OM91DbfScBqVZ3uxzzqcrOInO3ez8YZKXQJ0FlEHscZkvszd3pt66M+JwJvqTsMgf4ymmYWMN7dAovDGbysViLSAqcA7buq2CvAmz5NJrr/Vv+sTYizriETSKU+9yvZf5/UvjFcBFigqv3c2xGqeqo77QzgSZxfyDNFZN/ra5qvv+NCH+yYKnXNt7ZxaPKBHBFpXueMRYbhDLV9tKr2xSlaCepcQKgvzq/sG3BGjIXa10d9hJrf9+M4WxBHANcACX7Orzb7Ppfqn7UJcVYIjNeWAOkicjT8fLRNbxGJArJV9UvgdziXOEyuYz5fA6NEJFFEknAGDau+s/RrYIzbB58BnOBHPn/mux9VLQFeAB4TkTj3fWWIyEXVmrYAtqtqiYj0AAa7bdOAKFV9G/g/oH8D1oevKcB5ItLanf++6++2wBm5EuBSn/a7gQOKmKruBLb79P9fDAT9msOm8VnVNp5S1TJxDuF8zO16iMG5AMxS4FX3OQEeVtUdUssFYVR1loi8DPzoPvW8qv5Urdk7ON0k89z51/slVtt83Z2qdfkD8BdgoYjsxdl6qH4h+U+Aa0VkLk5B3NfNlAm85H75A9yN02Xm9/qo9h4WiMgDwFciUomz5XEZTt//myKyzl12J/cl7wNvichInEtz+roUeMbdIb4CZ2x8E+Zs9FFjjIlw1jVkjDERzgqBMcZEOCsExhgT4awQGGNMhLNCYIwxEc4KgTHGRDgrBMYYE+H+P3WpuQ2mcuMoAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(tvals,f1vals)\n",
"plt.ylabel('F1 Score')\n",
"plt.xlabel('Threshold for Classification')\n",
"plt.title('F1 as a function of Threshold');"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Based on this plot, we can say that the best classification threshold appears to be around 0.3, where precision and recall are:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "-"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Best Precision: 0.416, Best Recall: 0.695\n"
]
}
],
"source": [
"F1, Prec, Rec = evalThresh(df, 0.3)\n",
"print('Best Precision: {:0.3f}, Best Recall: {:0.3f}'.format(Prec, Rec))"
]
},
{
"cell_type": "markdown",
"metadata": {
"hide_input": true,
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The example here is based on\n",
"http://blog.yhathq.com/posts/logistic-regression-and-python.html\n",
"where you can find additional details."
]
}
],
"metadata": {
"anaconda-cloud": {},
"celltoolbar": "Tags",
"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.7.10"
},
"rise": {
"scroll": true,
"theme": "beige",
"transition": "fade"
}
},
"nbformat": 4,
"nbformat_minor": 1
}