{
"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": 28,
"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": true,
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEhCAYAAACUW2yNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt7UlEQVR4nO3deXhU5d3/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:
Tue, 07 Nov 2023
Pseudo R-squ.:
0.08107
\n",
"
\n",
"
\n",
"
Time:
09:21:06
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: Tue, 07 Nov 2023 Pseudo R-squ.: 0.08107\n",
"Time: 09:21:06 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,
"scrolled": true,
"slideshow": {
"slide_type": "-"
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAAFWCAYAAADtzQ3ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8RElEQVR4nO3deXxU1fnH8c8DhB2CYAAFwr4ICCgRd0TRuqKoFbe61KqttrXaatW2olKrtrX1p61LsVr3Bfd9RVHBBXAJyL6FJMi+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": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAFWCAYAAAAR2CYTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0IklEQVR4nO3deXyc1X3v8c/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": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAFWCAYAAADdQRz4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA69klEQVR4nO3deXxU9fX/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.584, Average Recall: 0.221\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": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0UUlEQVR4nO3dd3hUZfr/8fedngBJgARCQkloAqGpEZSiYFtAFFR2xbr2Rdeyrv5WXde2frfYdkXFde2uunZF1gIqKkUQCdJC7yT0mpCQnvv3xxxwCClDyLTM/bquuZiZ88xzPmeGzD2nPUdUFWOMMaErzN8BjDHG+JcVAmOMCXFWCIwxJsRZITDGmBBnhcAYY0KcFQJjjAlxVghMUBGXV0Vkn4j86ON5fyEiv/blPJ35/p+I7BaR7Q147UYROdsbuarN5zsRub6Br601o4gME5G840tn6mOFIAQ5f3jFIlLodkt1pr0gIqtEpEpErvZz1JoMAc4B2qvqAG/NREQeEpE33Z9T1ZGq+rq35llLjg7AnUAvVU2pNu1yt8+v2PnMDn+mvsxpgpsVgtB1vqo2d7ttdZ5fDNwM/OTHbHXpBGxU1SJ/B/GRTsAeVd1ZfYKqvnXo8wNGAlvdP9NjnZGIRDRCXhOErBCYI6jqJFWdDpTU11ZEzhORhSJSICK5IvKQ27QYEXlTRPaIyH4RmS8ibWvp5x4RWSciB0RkuYhcWEu764CXgNOcX70Pi8jVIjK7WjsVka7O/ddEZJKIfOb0P09Euri1zRSRr0Rkr4jsEJE/isgI4I/AJc58FjttD2/+EJEwEfmTiGwSkZ0i8h8RSXCmpTsZfi0im53NOvfV8T4mOK/f5fT3J6f/s4GvgFQnx2v1fSa16C8iS0QkX0TeFZEYZ77DRCRPRO52Nju96sz30OexR0TeE5FWTvv6PtNOIvK98z5/KSJJbst4gYgsc173nYj0rOW9iHU+s30ishw4pYHLbI6BFQJzPIqAq4BE4DzgJhEZ60z7NZAAdABaAxOA4lr6WQcMddo/DLwpIu2qN1LVl51+5jq/eh/0MOelTr8tgbXAXwBEpAXwNTAVSAW6AtNVdSrwV+BdZz79aujzauc2HOgMNAeerdZmCHACcBbwQG1ffsAzuJa9M3AGrvf0GlX9miN/6V/t4fJW9ytgBJAB9HVyH5ICtMK15nEjcBsw1smRCuwDJjlt6/tMLwOuAdoAUcBdACLSHXgb+B2QDHwO/E9EomrI+iDQxbn9wpmn8TIrBKFrsvPrbL+ITG5IB6r6naouVdUqVV2C64/9DGdyOa4vi66qWqmqC1S1oJZ+3lfVrU4/7wJrgMbc/v+Rqv6oqhXAW0B/5/nRwHZVfVJVS1T1gKrO87DPy4F/qOp6VS0E7gXGV9u88rCqFqvqYlyb3I4qKCISDlwC3OvMfyPwJHBlA5azNk877+9e4H/8vPwAVcCDqlqqqsXAb4D7VDVPVUuBh4BxznLV95m+qqqrnX7ec5vPJcBnqvqVqpYDTwCxwKAasv4K+Iuq7lXVXODpRnkHTJ2sEISusaqa6NzGNqQDERkoIt86mzTycf1CPLQ54A1gGvCOiGwVkcdEJLKWfq4SkUWHChPQ262fxuB+tM1BXL/ewfXLdl0D+0wFNrk93gREAO6bSmqbr7skXL+eq/eV1sBcNakrxy5Vdd8M2An42O2zWAFU4lqu+j7T2uZzxHulqlVALjUvY6oz7ZBNNbQxjcwKgTke/wWmAB1UNQF4HhAAVS1X1YdVtReuX36jcW3yOIKIdAJeBG4BWqtqIpBzqB8PFAFxbv2l1NG2ulxcmyBqUt+wvFtxfWke0hGoAHYcw/wBduP6pV29ry3H2E9DVV/OXGCk24+ERFWNUdUtnn6mNTjivRIRwVWEa1rGbc60Qzoey8KYhrFCYI4gIlHOzkQBIp0dhLX9P2kB7FXVEhEZgGsb8aF+hotIH2fTRwGuL7vKGvpohuvLaJfzumtwrRF4ajGQKSL9ndwPHcNrPwVSROR3IhItIi1EZKAzbQeQXseyvw3cISIZItKcn/cpVBzD/FHVSlybUf7izL8T8Hvgzbpf6TXPO1k6AYhIsoiMce57+plW9x5wnoic5axB3AmUAnNqaXuviLQUkfbArce/SKY+VghMdV/i2gE4CHjBuX96LW1vBv4sIgeAB3D9ER+SAnyA6wtjBTCDGr7cVHU5rm3ic3F9+fYBvvc0rKquBv6Ma6fvGmB23a844rUHcJ2TcD6uzRprcO38BXjf+XePiNR0KO0ruDaVzAQ24DrKqqFfWrfiWrNZjyv/f53+/WEirrW8L53P9QfgUHH06DOtTlVXAVfg2im+G9f7fb6qltXQ/GFcm4M24Pq/+MbxLIzxjNiFaYwxJrTZGoExxoQ4KwTGGBPirBAYY0yIs0JgjDEhLugGmUpKStL09HR/xzDGmKCyYMGC3aqaXNO0oCsE6enpZGdn+zuGMcYEFRGp9Sxt2zRkjDEhzgqBMcaEOCsExhgT4qwQGGNMiLNCYIwxIc4KgTHGhDgrBMYYE+KC7jwC0zSUVVSxbGs+i3L3MyCjFZmpCf6OZEzIskJgjkv+wXKWbsl3bvvZXVhGu4QYUhNjST30b2IsiXGRLN9aQPamfSzYuI/FefspragCoH3LWL664wxio8Lrnd/6XYXc/eESHhidSZ/2VjyMaQxWCIJMfnE5b/6wiZLySu44uzthYZ5e0bFxVFYp72fnMmvNbpZuyWfz3oOHp3VsFUfb+Gh+2ryPz5duo7zy6GtdRIQJmWkJXHFqJ07u1JIwESa8uYBnvlnDH0b0qHPeVVXK3R8uYf7Gffxp8lI+vnmwz5ffmKbICkGQ2F1YyiuzN/DG3E0cKP35aoh3nntCg/ssKa9k1prdFJaWM6pPO6Ij6v5FvnZnIX/4YDE/bd5PWmIs/TokcOmAjvRJS6B3WjyJcVGH21ZVKbsLS9myv5it+0vYW1RKt7Yt6Nc+8ahf/uNObs8LM9dz4YlpdGvbotb5v/HDJuZv3Mc5vdry1fIdfLxwCxef3L7By2+McQm6K5RlZWVpKI01tHV/MS/MXM878zdTWlHFqN7tuGlYF96Yu4l3s3N58pf9junLsKyiillrdvHZkm18tXzH4aKSEh/DTcO6cMkpHYiJPPKLuqKyihdmreepr9cQGxnOQxf0Ymz/NFzXID9+ewpLOfPJGZyQ0oJ3bzy1xn5z9x7kF0/N5JT0Vrx69Slc+Nz3bMsv4du7htEs2n7PGFMfEVmgqlk1TbO/oABQUl7JrgOlbC8oYUdBCTsKStlRUMLmPQeZvnIHqjD2xDRuGtaFLsnNAXhkbG827z3IPR8toUOrOAZktKpzHvPW7+H9BXl8uWw7BSUVxMdEMLJPCuf1TSVM4Jnpa3lwyjImfbuWCWd04bKBHYmJDGfl9gL+3/tLWLolnxGZKfx5bCZtWsQ06vK3bh7NPSN7cO9HS/nwpy2Mq1bYVJV7P1qKAH+9qA9hYcID5/fi4n/N5fkZ645rrcgYY2sEfpV/sJx/fr2aN3/YREXVkZ9DVEQYKfExDD8hmRtO70z7lnE1vv7C575n38EyPr55MOlJzY5qs6+ojEc+Xc5HC7fQIjqCczNTGN23HYO7JhEV8fPRw6rK3PV7eHr6Gn5Yv5ek5tEMPyGZyYu2EB8TySNjezOqT7vGfxMcVVXKuOfnsHHPQb6584wjNjO9Nz+XP3y4hEfG9ubKUzsdfv62txcybdl2pt95Ro3vjzHmZ3WtEVgh8IOqKuW97Fwem7aK/QfL+OXJHTi5U0vaJsSQEh9D2/hoEmIjPdr0snF3EWOf+55WzaL4+KbBJMRFAq4v9s+XbufBKTnsP1jOzcO6cPPwrkdt9qnJvPV7mDh9DXPW7WFs/1QeOD+TVs2i6n3d8VqxrYDRz8zmV1nt+dtFfQHYUVDC2f+YQa928bx9w6lH7Bzeur+YM5/8jrN6tmXSZSd5PZ8xwcw2DQWQhZv38eCUZSzJy+eU9JY8dMGA4zqGPj2pGf++4mSueHkeN721gNevHcC+ojL+NDmHL5fvoE9aAm9cN5Ce7eI97nNg59b8t3NrDpZVEBflu/8iPdvFc+3gdF6ctYFxJ7fnpI4tue/jHMorq3j04r5HHSGUmhjLb07vwsTpa7h60F5OSa9589iKbQVs2F3k1TUaY4KZV9cIRGQEMBEIB15S1b9Xmz4M+ATY4Dz1kar+ua4+g3WNYE9hKX/7YiUfLMijTYto/jiqJ2P6pzbaDtcPFuRx1/uLGdotiUW5+ymrqOL353TnuiEZRIQHzwnkRaUVnP2PGSTERvKbMzpzx7uLuW9UT244vXON7YvLKjnzye9Iah7NJ7898nDS3YWlPPnlat6dv5kqhT+PyeSq09J9tCTGBBa/rBGISDgwCTgHyAPmi8gUVV1ereksVR3trRyBYE9hKb/891xy9x7kN2d05tYzu9G8kY90GXdyezbsLmTSt+sYkNGKRy/uS0YN+wwCXbPoCB48P5MJby7gzvcW069DItcOyai1fWxUOPeM7MHt7yzig5/y+FVWB0orKnl9zkaemb6W4vJKrh6UwcY9RTw0ZRkdW8Ux7IQ2PlwiYwKfN9f7BwBrVXU9gIi8A4wBqheCJq2wtIJrXpvPln3FvHndQAZ2bu21ed117gmc1yeVHiktgvpEq19ktuWsHm2YtWY3j4/rS3g9y3JBv1Rem7ORx6etIiYynCe/XMWmPQcZfkIy953Xi65tmlNUWsEvn5/LLf9dyAc3nUaPFM83lRnT1Hlzm0EakOv2OM95rrrTRGSxiHwhIpk1dSQiN4pItohk79q1yxtZj8nanQf43TsLWbY1v852pRWV3PifbJZtLeC5y0/yahEAEBF6pcYHdREA13JMuvwkvvr96XSv4wQz9/YPnp/JrgOl3Pb2QiLDw3jtmlN49ZoBdG3jOty2WXQEL1+dRbPocK57LZudB0q8vRjGBA2v7SMQkV8Cv1DV653HVwIDVPVWtzbxQJWqForIKGCiqnarq19/7yPYdaCUsZO+Z8v+YiLChN+d3Y0JZ3Q5ajt8ZZXy27d+Yuqy7fzjV/246CQ7A9bb3pi7EUQYf0oHImvZL5KzJZ9fPj+X7ikteOeGUz0a38iYpqCufQTeXCPIAzq4PW4PbHVvoKoFqlro3P8ciBSRJC9mOi4l5ZXc+EY2e4pK+c+1AxjVpx1PfLmacc/PZd2uwsPtVJU/TV7K1GXbeWB0LysCPnLlaelceWqnWosAQO+0BCaO78+SvP3c+f4iqqqC6/BpY7zBm4VgPtBNRDJEJAoYD0xxbyAiKeIcNiMiA5w8e7yYqcGqqpS73l/Mws37eeqS/pzePZmnLz2RZy87kY17ijjv6Vm8+v0GqqqUx6et4u0fc7n1zK517ug0/nFuZgr3jerJ50u38/iXq/wdxxi/89rOYlWtEJFbgGm4Dh99RVWXicgEZ/rzwDjgJhGpAIqB8erjM9wqKqv4Imc7nVrH0bd9Yq3tnvp6NZ8u2cY9I3swovfPx6OP7pvKgPRW3P3hEh7+33LemreZtTsLuXxgR35/TncfLIFpiOuGZLB+dxH/+m4dzaMjuHlYl0Y7lNeYYBPSZxYv31rAHz5cTM6WAgDO6tGGO87pTu+0I0/w+nBBHne+v5hLsjrw94v71PiFoaq8Oz+XRz5dzpk92/LUJf3rPdrF+Fd5ZRV3vb+YTxZt5dIBHXlkTGZQnXNhzLGwISaqKa2o5Nlv1vKv79aRGBfJ/aN7kbv3IC/O2kB+cTnn9GrL787uRmZqAvPW7+GKl+eR1akVr1874IjxeWpSUl5JdESY/boMElVVyhNfruK579Yx/IRknr3sJBvN1DRJVgjc/LR5H3d/sIQ1Owu56KQ07j+vFy2dcXQKSsp5dfZGXpq9ngMlFfwisy3zNuw9ahwf0/S8NW8T90/OoVdqPK9cfUqjj7BqjL9ZIQAOllXw5JereeX7DaTEx/DXi/owvJYzTPOLy3l59gZenb2BiHCpdWRP07R8s3IHv31rIa2aRfH6tafQtU395zAYEyysEADvZefyhw+WcMWpHbl7RA9axNT/676gpJyyiiqSmkc3JKoJQkvy9nPta/MPj9XUrW0L0pOa0S4+JuhP1DOhzQoBrm3BOVvz6zwyyBhwXQ3tutfns3rHz+eGREWE0alVHOlJzRjcpTW/HpRu+4FMULFhqIGwMLEiYDzSoVUcU28/ne0FJWzcU8TG3QfZtKeIDbuLWLOzkK+W76Bnu3ivDxlijK+ETCEw5liEhQmpibGkJsYyqMvPz5eUVzLk0W+Y9N06KwSmybCDpo05BjGR4VwzOIOZq3eRs6XuQQeNCRZWCIw5Rlee1okW0RE8991af0cxplFYITDmGMXHRHLVoE58kbOdtTsL63+BMQHOCoExDXDN4AyiI8J4fsY6f0cx5rhZITCmAZKaRzP+lI5MXriFLfuL/R3HmONihcCYBrrh9M4AvDhzvZ+TGHN8rBAY00BpibFceGIab/+4md2Fpf6OY0yDWSEw5jhMGNaFssoqXv1+g7+jGNNgVgiMOQ5dkpszsncK/5mziYKScn/HMaZBrBAYc5xuHtaVA6UVvDF30+HnCksr+GblDh75dDmjJs7in1+t9mNCY+pmQ0wYc5x6pyVwevdkXpm9gbKKKuas283CzfupqFKiIsJoFRfFCzPXc93QDOI9GPXWGF+zNQJjGsEtw7uyp6iMp79ZQ1lFFTec3pm3rh/IkgfP5d9XnkxxeSWfLNzi75jG1MjWCIxpBAMyWvHF7UNJTYg96kp2fdsnkJkaz1vzNnPFqZ1s+GoTcGyNwJhG0rNdfI2XMxURLhvYkZXbD7Aod7/vgxlTDysExvjAmP5pxEWF8995m/0dxZijWCEwxgeaR0cwpn8q/1uylfxiO8zUBBYrBMb4yGUDOlFSXsVk22lsAowVAmN8pE/7BPqkJfD2j5sJtmuFm6bNCoExPnRop/FPm/f7O4oxh1khMMaHLuiXSjPbaWwCjBUCY3yoWXQEY05M49MlW8k/aDuNTWCwQmCMj102oCOlFVV8vDDP31GMAbxcCERkhIisEpG1InJPHe1OEZFKERnnzTzGBILeaQn0a5/Af22nsQkQXisEIhIOTAJGAr2AS0WkVy3tHgWmeSuLMYHmsoEdWb2jkAWb9vk7ijFeXSMYAKxV1fWqWga8A4ypod2twIfATi9mMSagjO6bSvPoCNtpbAKCNwtBGpDr9jjPee4wEUkDLgSer6sjEblRRLJFJHvXrl2NHtQYX2sWHcGFJ6YxZfFW/vHVakrKK/0dyYQwbxaCmoZYrL5B9CngblWt869AVV9Q1SxVzUpOTm6sfMb41V3nnsB5fdvx9PQ1nPvPmXyzcoe/I5kQ5c1CkAd0cHvcHtharU0W8I6IbATGAc+JyFgvZjImYCTERTJx/In894aBREWEce1r2dzwn2zy9h30dzQTYrxZCOYD3UQkQ0SigPHAFPcGqpqhqumqmg58ANysqpO9mMmYgDOoSxKf3zaUu0f0YPaa3Zz9jxlM+nYt5ZVV/o5mQoTXCoGqVgC34DoaaAXwnqouE5EJIjLBW/M1JhhFRYRx07AufH3nGQzr3obHp63iH3adY+MjEmzHMWdlZWl2dra/YxjjVXe9v5iPF27h01uH0LNdvL/jmCZARBaoalZN0+zMYmMC0H2jepIQG8k9Hy2lsiq4fqyZ4GOFwJgA1LJZFA+M7sXi3P28MXejv+OYJs4KgTEBakz/VIZ2S+LxaavYur/Y33FME2aFwJgAJSL8ZWwfKlV54JMcG5fIeI0VAmMCWMfWcdxxdne+XrGTqTnb/R3HNFFWCIwJcNcNyaBXu3genLLMLnxvvMIKgTEBLiI8jL9d1IfdhaU8NnXl4ecrKqtYsGkfz0xfw/gX5nL+M7MpLrMxi8yxi/B3AGNM/fp1SOTqQRm88v0GWsZFsXL7Aeat38OB0gpEIL11MzbsLmLG6l2M6J3i77gmyFghMCZI3Hlud6Yt286z366lU+s4RvdLZUjXJE7r0poWMRGc8pevmZqzzQqBOWZWCIwJEs2iI/jklsEUl1XSoVXcUdPP6dmWqTnbKa2oJDoi3A8JTbCyfQTGBJGk5tE1FgGAkX1SOFBawZy1e3ycygQ7KwTGNBGDuybRIjqCL3K2+TuKCTJWCIxpIqIjwjmzZxu+Wr6DChvC2hwDKwTGNCEje6ew72A58zbs9XcUE0SsEBjThJzRvQ2xkeG2ecgcEysExjQhsVHhDO+RzLRlO6iy4auNh6wQGNPEjOjdjl0HSlmweZ+/o5ggYYXAmCbmzB5tiIoI44ulNkid8YwVAmOamObREZzeLYmpOdts6GrjESsExjRBI3q3Y2t+CYvz8v0dxQQBKwTGNEHn9GxLRJjY0UPGI1YIjGmCEuIiOa1La6bmbLfNQ6ZeVgiMaaJG9m7Hpj0HWbHtgL+jmABnhcCYJurczLaECUy1zUOmHh4VAhEZIiLXOPeTRSTDu7GMMccrqXk0AzJa8YVd69jUo95CICIPAncD9zpPRQJvejOUMaZxjOzdjjU7C1m7s9DfUUwA82SN4ELgAqAIQFW3Ai28GcoY0zhG9E4hIkx4adZ6f0cxAcyTQlCmrsMOFEBEmnk3kjGmsbSNj+HqQem8m53Lotz9/o5jApQnheA9Efk3kCgiNwBfAy960rmIjBCRVSKyVkTuqWH6GBFZIiKLRCRbRIYcW3xjTH1uP7sbSc2jefCTHBuIztSozkIgIgK8C3wAfAicADygqs/U17GIhAOTgJFAL+BSEelVrdl0oJ+q9geuBV461gUwxtStRUwk943qyeK8fN7LzvV3HBOA6rx4vaqqiExW1ZOBr46x7wHAWlVdDyAi7wBjgOVu/bvvwWqGs/nJGNO4xvRP5b/zNvPo1JWM6J1CYlyUvyOZAOLJpqEfROSUBvSdBrj//MhznjuCiFwoIiuBz3CtFRxFRG50Nh1l79q1qwFRjAltIsLDYzIpKKngiS9X+TuOCTCeFILhuIrBOmd7/lIRWeLB66SG5476xa+qH6tqD2As8EhNHanqC6qapapZycnJHszaGFNdz3bxXHlqJ96at5mcLTYYnfmZJ4VgJNAZOBM4Hxjt/FufPKCD2+P2wNbaGqvqTKCLiCR50LcxpgHuOKc7rZtF8YDtODZu6i0EqroJSMT15X8+kOg8V5/5QDcRyRCRKGA8MMW9gYh0dXZIIyInAVHAnmNaAmOMxxJiI7l7RA9+2ryfD3/K83ccEyA8ObP4duAtoI1ze1NEbq3vdapaAdwCTANWAO+p6jIRmSAiE5xmFwM5IrII1xFGl6gNlWiMV118UntO6pjIo1NXkl9c7u84JgBIfd+7zv6A01S1yHncDJirqn19kO8oWVlZmp2d7Y9ZG9Nk5GzJ5/xnZzOqTzseH9eXuKg6DyA0TYCILFDVrJqmebKPQIBKt8eV1Lwj2BgTJHqnJfD7s7vz2ZJtjHhqFnPX2RbZUOZJIXgVmCciD4nIQ8APwMteTWWM8bpbz+rGuzeeighc+uIP3D85h6LSCn/HMn5Q76YhOLwjdwiuNYGZqrrQ28FqY5uGjGlcxWWVPD5tFa/O2UBaYiyPXdyXQV3t4L2m5rg2DYnIqcAaVX1aVScCa0VkYGOHNMb4R2xUOA+c34v3f3MakeFhXPbSPO6fnEOlHV4aMjzZNPQvwH0oiCLnOWNME5KV3orPbxvKtYMzeOOHTTw2baW/Ixkf8eRQAXE/pFNVq0TEDjEwpgk6tHZQWlHJv2esp3dqAuf3S/V3LONlnqwRrBeR20Qk0rndDthVLoxpwh48P5OsTi35wwdLWL61wN9xjJd5UggmAIOALc5tIHCjN0MZY/wrKiKM5644ifjYCG58I5t9RWX+jmS8yJMhJnaq6nhVbePcLlPVnb4IZ4zxnzYtYnj+ipPZWVDKrW8vpKKyyt+RjJfUWghE5AYR6ebcFxF5RUTynRFIT/JdRGOMv5zYsSX/N7Y3s9fu5tGptvO4qaprjeB2YKNz/1KgH65RSH8PTPRuLGNMoPjVKR246rROvDhrA58s2uLvOMYL6ioEFap6aESq0cB/VHWPqn6N62pixpgQcf/oXgxIb8XdHy5h5XbbedzU1FUIqkSknYjEAGfhumj9IbHejWWMCSSR4WFMuvwkYiLDefQL20TU1NRVCB4AsnFtHpqiqssAROQM7PBRY0JOcotobhjamW9X7WJx7n5/xzGNqNZCoKqfAp2Anqp6g9ukbOASbwczxgSeq07rRGJcJBOnr/F3FNOI6jx8VFUrVHVfteeKVLWwttcYY5quFjGRXD8kg29W7mRJ3n5/xzGNxJMTyowx5rBfD0onITaSp22toMmwQmCMOSaH1gq+XrGTpXn5/o5jGkGDCoGI9GjsIMaY4PHrwenEx0TYvoImoqFrBF82agpjTFCJj4nk+qGd+XrFDnK22FpBsKt1OGkRebq2SUCiV9IYY4LG1YPTeWnWeiZOX8OLV9V44SsTJOpaI7gGyAEWVLtlAzYUoTEhLj4mkuuGdOar5bZWEOzqKgTzgRxVfb36DTjgo3zGmAB29eB0WsRE2BFEQa6uQjAOWFTTBFXN8EoaY0xQSYiN5LohGXy5fIddwCaI1VUImqvqQZ8lMcYEpWsGZ9haQZCrqxBMPnRHRD70fhRjTDBKiI3kmkHpTF22nTU7bKtxMKqrEIjb/c7eDmKMCV7XDM4gLiqc575b5+8opgHqKgRay31jjDlCy2ZRXD6wI1MWb2XzHtuiHGzqKgT9RKRARA4AfZ37BSJyQERsr5Ax5gg3DO1MeJjwrxm2VhBs6hqGOlxV41W1hapGOPcPPY73pHMRGSEiq0RkrYjcU8P0y51rIC8RkTki0u94FsYY4z9t4mP4VVZ7PlyQx/b8En/HMcfAa4POiUg4MAkYCfQCLhWRXtWabQDOUNW+wCPAC97KY4zxvt+c3oVKVV6YadeuCibeHH10ALBWVderahnwDjDGvYGqznG73sEPQHsv5jHGeFmHVnGM6Z/Kf3/cxJ7CUn/HMR7yZiFIA3LdHuc5z9XmOuCLmiaIyI0iki0i2bt27WrEiMaYxnbzsK6UVlTxyvcb/B3FeMibhUBqeK7Go49EZDiuQnB3TdNV9QVVzVLVrOTk5EaMaIxpbF3bNGdk7xT+M2cTBSXl/o5jPODNQpAHdHB73B7YWr2RiPQFXgLGqOoeL+YxxvjIzcO6cqC0gjfmbvJ3FOMBbxaC+UA3EckQkShgPDDFvYGIdAQ+Aq5U1dVezGKM8aHeaQkMPyGZl2dv4GBZhb/jmHp4rRCoagVwCzANWAG8p6rLRGSCiExwmj0AtAaeE5FFIpLtrTzGGN+65cyu7C0q4+0fc+tvbPxKVIPrpOGsrCzNzrZ6YUwwGP/CXDbsLmLG/xtOTGS4v+OENBFZoKo1XkHILl5vjPGa353dnR0FpUz6dq2/o5g6WCEwxnjNqZ1bc+GJaTw/Yx1rd9rIpIHKCoExxqvuO68ncVER3PdxDsG2KTpUWCEwxnhVUvNo7hnZg3kb9vLhT1v8HcfUwAqBMcbrLsnqwMmdWvLXz1ewr6jM33FMNVYIjDFeFxYm/OXC3hQUl/O3L1bU2TZnSz6r7UpnPmWFwBjjEz1S4rluaAbvZefx44a9R03ff7CMez9ayuhnZjPhjQV+SBi6rBAYY3zm9rO6kZYYy30fL6WsogoAVeXjhXmc9eQM3svOpV+HRNbvLiJ3r13pzFesEBhjfCYuKoI/j8lkzc5CXpy1nvW7Crni5Xnc8e5iOrSK43+3DOGJcX0BmL12t5/Tho4IfwcwxoSWs3q2ZURmChOnr2Hi9DVEh4fxyJhMLhvYifAwQVVJiY9h9prdXDqgo7/jhgQrBMYYn3vwgl4smrSfrPSWPDC6F23iYw5PExGGdEviq+U7qKxSwsNqGtHeNCYrBMYYn2uXEMvce89EpOYv+aHdkvhgQR45W/Lp1yHRt+FCkO0jMMb4RW1FAGBw1yTA9hP4ihUCY0zASWoeTa928cxcbZem9QUrBMaYgDS0exI/bd5HUald2MbbrBAYYwLS0K7JlFcq8zbYFWy9zQqBMSYgZaW3JDoijFlrbD+Bt1khMMYEpJjIcAZktLJC4ANWCIwxAev0bsms3VnItvxif0dp0qwQGGMC1pBursNIba3Au6wQGGMCVo+UFiQ1j2a2FQKvskJgjAlYIsLQbknMXrubqiq7zKW3WCEwxgS0IV2T2FtUxvJtBf6O0mRZITDGBLSh3Wy4CW+zQmCMCWht4mM4oW0LZq2x4Sa8xQqBMSbgDe2WxPyN+yguq/R3lCbJCoExJuAN6ZZEWUUVP248+lrH5vhZITDGBLyBGa2JCg9jtm0e8gqvFgIRGSEiq0RkrYjcU8P0HiIyV0RKReQub2YxxgSv2KhwstJb2ollXuK1QiAi4cAkYCTQC7hURHpVa7YXuA14wls5jDFNw9BuyazcfoAHPskhd+9Bf8dpUry5RjAAWKuq61W1DHgHGOPeQFV3qup8oNyLOYwxTcCVp3XiV1ntefvHzQx74jvueHcRq3cc8HesJsGbhSANyHV7nOc8d8xE5EYRyRaR7F27bBuhMaGoeXQEj43rx8w/DOfqQelMzdnOuf+cyfWvZ/PT5n3+jhfUvFkIarogaYPOEVfVF1Q1S1WzkpOTjzOWMSaYtUuI5f7RvZhzz5ncflY3sjft5aLn5vDy7A3+jha0vFkI8oAObo/bA1u9OD9jTAhp2SyKO87pzvd3n8mwE5J5Ytoqtu634aobwpuFYD7QTUQyRCQKGA9M8eL8jDEhqFl0BI+M6U2VKv/32XJ/xwlKXisEqloB3AJMA1YA76nqMhGZICITAEQkRUTygN8DfxKRPBGJ91YmY0zT1KFVHLcM78rnS7czc7XtRzxWohpcQ7tmZWVpdna2v2MYYwJMaUUlI56aBcDU3w0lOiLcz4kCi4gsUNWsmqbZmcXGmCYhOiKchy7IZMPuIl6cud7fcYKKFQJjTJNxRvdkRvZO4dlv19pJZ8fACoExpkm5f3QvBOHh/9mOY09ZITDGNCmpibHcdlY3vl6xg+krdvg7TlCwQmCMaXKuG5JB1zbNeeh/yygpt2sY1McKgTGmyYmKCOPPF2SSu7eY575b5+84Ac8KgTGmSRrUNYnRfdvx7xnr2GJnHNfJCoExpsm6d1RPAB6futLPSQKbFQJjTJOVlhjL9UMzmLxoK4ty9/s7TsCyQmCMadJuGtaVpObRPPLpcoJtJAVfsUJgjGnSmkdHcNe53VmwaR+fL93u7zgByQqBMabJ+2VWB3qktOBvX6yww0lrYIXAGNPkhYcJ94/uRd6+Yl6bs9HfcQKOFQJjTEgY3DWJs3u24dlv1rK7sLTGNht3F/GXz5bzxdJtIbU/wQqBMSZk3DuqJyXllfzzq9VHPL+zoIT7Pl7K2f+YwYuzNnDTWz9x4xsL2J5f4qekvmWFwBgTMrokN+eKUzvx9o+bWb3jAPnF5Tw2dSWnP/4t787P5bKBHZl775n8cVQPZq3ZxTn/mMFb8zZRVdW01w7swjTGmJCyr6iMYU98R5sW0ew8UEp+cTlj+qfy+3O606l1s8PtNu0p4t6PljJn3R4GZLTibxf1oUtycz8mPz52YRpjjHG0bBbFHWd3Y83OQk7smMhntw1h4vgTjygCAJ1aN+Ot6wfy2Li+rNxWwMiJs3hp1vomue/A1giMMSFHVdleUEK7hFiP2u88UML9k3OYtmwHN57emXtH9kBEvJyycdkagTHGuBERj4sAQJsWMfzr8pO56rROvDBzPX+anNOk9htE+DuAMcYEg7Aw4eELMomLiuD5GesoLqvksXF9iQgP/t/TVgiMMcZDIsLdI06geXQ4T3y5muLySiaOP5GoiOAuBsGd3hhjfExEuOXMbtw/uhdf5Gznxjeyg37YCisExhjTANcNyeDvF/VhxupdXPnyPGav2U15ZZW/YzWIbRoyxpgGGj+gI7FR4dz70VKueHkeiXGRnN2zLSN7pzC4axIxkeH+jugRKwTGGHMcxvRP49xeKcxcs4upOduZtmw7HyzIo1lUOMN7tOHik9tzRrdkwsIC93BTKwTGGHOcYqPC+UVmCr/ITKGsooo563Yzbdl2pi3bwadLttGpdRxXntqJX57cgYS4SH/HPYqdUGaMMV5SVlHF1GXb+c+cjWRv2kdMZBhj+6dx5WmdyExNoLiskm35xWzLL2Hrfte/RaUVjO6bSp/2CY2apa4TyqwQGGOMDyzbms8bczcxedEWSsqriI+JoKCk4qh2EWFCRZWS1akl1wzO4BeZbRvlXAW/FQIRGQFMBMKBl1T179WmizN9FHAQuFpVf6qrTysExphgln+wnPcX5LJxTxHtEmJplxBDu4RYUhNjaBsfQ1llFe/Nz+X1uRvJ3VtMakIMVw1KZ/wpHUiMi2rwfP1SCEQkHFgNnAPkAfOBS1V1uVubUcCtuArBQGCiqg6sq18rBMaYUFBZpUxfsYNXv9/I3PV7iIkM465zT+D6oZ0b1F9dhcCbO4sHAGtVdb0T4h1gDLDcrc0Y4D/qqkY/iEiiiLRT1W1ezGWMMQEvPEw4NzOFczNTWLGtgNe+30hqoufjIx0LbxaCNCDX7XEerl/99bVJA44oBCJyI3AjQMeOHRs9qDHGBLKe7eJ5dFxfr/XvzTOLazpotvp2KE/aoKovqGqWqmYlJyc3SjhjjDEu3iwEeUAHt8ftga0NaGOMMcaLvFkI5gPdRCRDRKKA8cCUam2mAFeJy6lAvu0fMMYY3/LaPgJVrRCRW4BpuA4ffUVVl4nIBGf688DnuI4YWovr8NFrvJXHGGNMzbw6xISqfo7ry979uefd7ivwW29mMMYYUzcbhtoYY0KcFQJjjAlxVgiMMSbEBd2gcyKyC9jk7xxukoDd/g7hgWDIaRkbTzDkDIaMEBw5PcnYSVVrPBEr6ApBoBGR7NrG7wgkwZDTMjaeYMgZDBkhOHIeb0bbNGSMMSHOCoExxoQ4KwTH7wV/B/BQMOS0jI0nGHIGQ0YIjpzHldH2ERhjTIizNQJjjAlxVgiMMSbEWSHwkIiMEJFVIrJWRO6pYXoPEZkrIqUicleAZrxcRJY4tzki0i9Ac45xMi4SkWwRGRJoGd3anSIilSIyzpf5nHnX9z4OE5F8531cJCIP+DqjJzmdNsOcjMtEZEagZRSR/+f2PuY4n3mrAMyZICL/E5HFznvp2UCeqmq3em64Rk9dB3QGooDFQK9qbdoApwB/Ae4K0IyDgJbO/ZHAvADN2Zyf91/1BVYGWka3dt/gGlhxXKBlBIYBn/r6M25AzkRcl7Dt6DxuE2gZq7U/H/gmQN/LPwKPOveTgb1AVH192xqBZw5ff1lVy4BD118+TFV3qup8oNwfAfEs4xxV3ec8/AHXhYB8zZOcher8TwaaUcNV6/yd0XEr8CGw05fhHJ5m9DdPcl4GfKSqm8H1txSAGd1dCrztk2RH8iSnAi1ERHD9oNoLVNTXsRUCz9R2beVAcqwZrwO+8GqimnmUU0QuFJGVwGfAtT7Kdki9GUUkDbgQeB7/8PTzPs3ZTPCFiGT6JtoRPMnZHWgpIt+JyAIRucpn6Vw8/tsRkThgBK4fAL7mSc5ngZ64rvS4FLhdVavq69ir1yNoQjy6trKfeZxRRIbjKgQ+3/aO59ep/hj4WEROBx4BzvZ2MDeeZHwKuFtVK10/vnzOk4w/4RpfplBERgGTgW7eDlaNJzkjgJOBs4BYYK6I/KCqq70dznEsf9/nA9+r6l4v5qmNJzl/ASwCzgS6AF+JyCxVLairY1sj8EwwXFvZo4wi0hd4CRijqnt8lM3dMb2XqjoT6CIiSd4O5saTjFnAOyKyERgHPCciY32SzqXejKpaoKqFzv3PgUgfv4/g+bXLp6pqkaruBmYCvjyQ4Vj+T47HP5uFwLOc1+DazKaquhbYAPSot2df7/AIxhuuXyzrgQx+3kmTWUvbh/DPzuJ6MwIdcV0WdFAgv5dAV37eWXwSsOXQ40DJWK39a/h+Z7En72OK2/s4ANjsy/fxGHL2BKY7beOAHKB3IGV02iXg2ubezJfv4TG+l/8CHnLut3X+dpLq69s2DXlAPbj+soikANlAPFAlIr/DtUe/zlUyX2YEHgBa4/r1ClChPh5V0cOcFwNXiUg5UAxcos7/7ADK6FceZhwH3CQiFbjex/G+fB89zamqK0RkKrAEqAJeUtWcQMroNL0Q+FJVi3yVrQE5HwFeE5GluDYl3a2utaw62RATxhgT4mwfgTHGhDgrBMYYE+KsEBhjTIizQmCMMSHOCoExxoQ4KwTmuIhIa7dRGbeLyBbn/n4RWe6F+T10rKO7ikhhLc+/VtOooeIaSXaRiCwUkS4NzBkpIn8XkTXOaJU/ishIZ9rGxjqxS0QuODQKpYgki8g8J/dQEflcRBIb0OcwERnk9niCH4Z9MD5k5xGY46Kus5P7g+tLGihU1SdEJB34tL7Xi0iEqtY7KJaPjQU+UdUHPWnsDPAleuSYLo8A7XCdGFUqIm2BMxo7qKpOAaY4D8/CNVLrr53HsxrY7TCgEJjjzMPv500Y77I1AuNN4SLyojMu+pciEgvgDC72V2fc+dtF5GQRmeEMODZNRNo57W4TkeXiujbBO2799nL6WC8itx16UkR+7/z6znFO6DuCuDzr9PkZrqHDq7cZBfwOuF5Evq2tXxFJF5EVIvIcrjF9Orj1EQfcANyqqqUAqrpDVd+rYX6TneVeJiI3Os+FO2srOSKyVETuqO39EJGrnWXqDzwGjHLWZmLd1zxE5CrndYtF5A3nufPd1iC+FpG2TgGfANzh9DPUfS1MRPqLyA9OXx+LSEu3z/RRZ81ntYgMrf2/hQk4/jhV2m5N84bb8BpAOq7hb/s7j98DrnDufwc859yPxPXLM9l5fAmuMybBNY5KtHM/0W0ec4BoIAnY4/RxMq7RFpvhGn53GXCi85pC59+LgK9wnZWZCuynhqEhqi1Hjf06y1cFnFrD6/sCC+t4nzbinPYPtHL+jcU1tEJrZ55fubVPrOP9uBp4tvp99/kAmcCqGubZkp9PKr0eeLL68tfwfiwBznDu/xl4yu0zPfT6UcDX/v7/aDfPb7ZpyHjTBlVd5NxfgOvL85B3nX9PAHrjGiURXF/S25xpS4C3RGQyrpEzD/lMXb+0S0VkJ64xVYYAH6tz+r+IfAQMBRa6ve504G1VrQS2isg3HixDbf1OATap6g8e9FGX20TkQud+B1yjg64COovIM7iG4f7SmV7b+1GfM4EP1BlqQH8eObM98K6zBhaFa4CyWolIAq4CdOgKYq8D77s1+cj5t/pnbQKcbRoy3lTqdr+SI/dJHRqvRYBlqtrfufVR1XOdaecBk3D9Ql4gIodeX1O/no4FfaxjqtTVb21jzqwFOopIizo7FhmGa3jt01S1H66iFaOuiwf1w/Ur+7e4RouF2t+P+gg1L/czuNYg+gC/AWI87K82hz6X6p+1CXBWCIy/rQKSReQ0OHy0TaaIhAEdVPVb4A+4LmfYvI5+ZgJjRSRORJrhGiCs+s7SmcB4Zxt8O2C4B/k86fcIqnoQeBl4WkSinOVqJyJXVGuaAOxT1YMi0gM41WmbBISp6ofA/cBJDXg/3E0HfiUirZ3+D11rNwHX6JQAv3ZrfwA4qoipaj6wz237/5WAz68vbBqfVW3jV6paJq5DOJ92Nj1E4Lroy2rgTec5Af6pqvullovAqOpPIvIa8KPz1EuqurBas49xbSZZ6vRf75dYbf06O1Xr8ifg/4DlIlKCa+2h+sXjpwITRGQJroJ4aDNTGvCq8+UPcC+uTWYevx/VlmGZiPwFmCEilbjWPK7Gte3/fRHZ4sw7w3nJ/4APRGQMrstxuvs18LyzQ3w9rvHvTZCz0UeNMSbE2aYhY4wJcVYIjDEmxFkhMMaYEGeFwBhjQpwVAmOMCXFWCIwxJsRZITDGmBD3/wFfU2g7OUAhQQAAAABJRU5ErkJggg==\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.462, Best Recall: 0.663\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 (ipykernel)",
"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.9.12"
},
"rise": {
"scroll": true,
"theme": "beige",
"transition": "fade"
}
},
"nbformat": 4,
"nbformat_minor": 1
}