project2logistics_regressionday17.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Logistic Regression\n",
    "\n",
    "\n",
    "Researchers are often interested in setting up a model to analyze the relationship between predictors (i.e., independent variables) and it's corresponsing response (i.e., dependent variable). Linear regression is commonly used when the response variable is continuous.  One assumption of linear models is that the residual errors follow a normal distribution. This assumption fails when the response variable is categorical, so an ordinary linear model is not appropriate. This newsletter presents a regression model for response variable that is dichotomous–having two categories. Examples are common: whether a plant lives or dies, whether a survey respondent agrees or disagrees with a statement, or whether an at-risk child graduates or drops out from high school.\n",
    "\n",
    "In ordinary linear regression, the response variable (Y) is a linear function of the coefficients (B0, B1, etc.) that correspond to the predictor variables (X1, X2, etc.,). A typical model would look like:\n",
    "\n",
    "    Y = B0 + B1*X1 + B2*X2 + B3*X3 + … + E\n",
    "\n",
    "For a dichotomous response variable, we could set up a similar linear model to predict individual category memberships if numerical values are used to represent the two categories. Arbitrary values of 1 and 0 are chosen for mathematical convenience. Using the first example, we would assign Y = 1 if a plant lives and Y = 0 if a plant dies.\n",
    "\n",
    "This linear model does not work well for a few reasons. First, the response values, 0 and 1, are arbitrary, so modeling the actual values of Y is not exactly of interest. Second, it is the probability that each individual in the population responds with 0 or 1 that we are interested in modeling. For example, we may find that plants with a high level of a fungal infection (X1) fall into the category “the plant lives” (Y) less often than those plants with low level of infection. Thus, as the level of infection rises, the probability of plant living decreases.\n",
    "\n",
    "Thus, we might consider modeling P, the probability, as the response variable. Again, there are problems. Although the general decrease in probability is accompanied by a general increase in infection level, we know that P, like all probabilities, can only fall within the boundaries of 0 and 1. Consequently, it is better to assume that the relationship between X1 and P is sigmoidal (S-shaped), rather than a straight line.\n",
    "\n",
    "It is possible, however, to find a linear relationship between X1 and function of P. Although a number of functions work, one of the most useful is the logit function. It is the natural log of the odds that Y is equal to 1, which is simply the ratio of the probability that Y is 1 divided by the probability that Y is 0. The relationship between the logit of P and P itself is sigmoidal in shape. The regression equation that results is:\n",
    "\n",
    "    ln[P/(1-P)] = B0 + B1*X1 + B2*X2 + …\n",
    "\n",
    "Although the left side of this equation looks intimidating, this way of expressing the probability results in the right side of the equation being linear and looking familiar to us. This helps us understand the meaning of the regression coefficients. The coefficients can easily be transformed so that their interpretation makes sense.\n",
    "\n",
    "The logistic regression equation can be extended beyond the case of a dichotomous response variable to the cases of ordered categories and polytymous categories (more than two categories)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mathematics behind Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem structure is the classic classification problem. Our data set $\\mathcal{D}$ is composed of $N$ samples. Each sample is a tuple containing a feature vector and a label. For any sample $n$ the feature vector is a $d+1$ dimensional column vector denoted by ${\\bf x}_n$ with $d$ real-valued components known as features. Samples are represented in homogeneous form with the first component equal to $1$: $x_0=1$. Vectors are bold-faced. The associated label is denoted $y_n$ and can take only two values: $+1$ or $-1$.\n",
    "\n",
    "$$\n",
    "\\mathcal{D} = \\lbrace ({\\bf x}_1, y_1), ({\\bf x}_2, y_2), ..., ({\\bf x}_N, y_N) \\rbrace \\\\\n",
    "{\\bf x}_n = \\begin{bmatrix} 1 & x_1 & ... & x_d \\end{bmatrix}^T \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The learning algorithm is how we search the set of possible hypotheses (hypothesis space $\\mathcal{H}$) for the best parameterization (in this case the weight vector ${\\bf w}$). This search is an optimization problem looking for the hypothesis that optimizes an error measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no sophisticted, closed-form solution like least-squares linear, so we will use gradient descent instead. Specifically we will use batch gradient descent which calculates the gradient from all data points in the data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Luckily, our \"cross-entropy\" error measure is convex so there is only one minimum. Thus the minimum we arrive at is the global minimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gradient descent is a general method and requires twice differentiability for smoothness. It updates the parameters using a first-order approximation of the error surface.\n",
    "\n",
    "$$\n",
    "{\\bf w}_{i+1} = {\\bf w}_i + \\nabla E_\\text{in}({\\bf w}_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To learn we're going to minimize the following error measure using batch gradient descent.\n",
    "\n",
    "$$\n",
    "e(h({\\bf x}_n), y_n) = \\ln \\left( 1+e^{-y_n \\; {\\bf w}^T {\\bf x}_n} \\right) \\\\\n",
    "E_\\text{in}({\\bf w}) = \\frac{1}{N} \\sum_{n=1}^{N} e(h({\\bf x}_n), y_n) = \\frac{1}{N} \\sum_{n=1}^{N} \\ln \\left( 1+e^{-y_n \\; {\\bf w}^T {\\bf x}_n} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll need the derivative of the point loss function and possibly some abuse of notation.\n",
    "\n",
    "$$\n",
    "\\frac{d}{d{\\bf w}} e(h({\\bf x}_n), y_n)\n",
    "= \\frac{-y_n \\; {\\bf x}_n \\; e^{-y_n {\\bf w}^T {\\bf x}_n}}{1 + e^{-y_n {\\bf w}^T {\\bf x}_n}}\n",
    "= -\\frac{y_n \\; {\\bf x}_n}{1 + e^{y_n {\\bf w}^T {\\bf x}_n}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the point loss derivative we can determine the gradient of the in-sample error:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla E_\\text{in}({\\bf w})\n",
    "&= \\frac{d}{d{\\bf w}} \\left[ \\frac{1}{N} \\sum_{n=1}^N e(h({\\bf x}_n), y_n) \\right] \\\\\n",
    "&= \\frac{1}{N} \\sum_{n=1}^N \\frac{d}{d{\\bf w}} e(h({\\bf x}_n), y_n) \\\\\n",
    "&= \\frac{1}{N} \\sum_{n=1}^N \\left( - \\frac{y_n \\; {\\bf x}_n}{1 + e^{y_n {\\bf w}^T {\\bf x}_n}} \\right) \\\\\n",
    "&= - \\frac{1}{N} \\sum_{n=1}^N \\frac{y_n \\; {\\bf x}_n}{1 + e^{y_n {\\bf w}^T {\\bf x}_n}} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our weight update rule per batch gradient descent becomes\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "{\\bf w}_{i+1} &= {\\bf w}_i - \\eta \\; \\nabla E_\\text{in}({\\bf w}_i) \\\\\n",
    "&= {\\bf w}_i - \\eta \\; \\left( - \\frac{1}{N} \\sum_{n=1}^N \\frac{y_n \\; {\\bf x}_n}{1 + e^{y_n {\\bf w}_i^T {\\bf x}_n}} \\right) \\\\\n",
    "&= {\\bf w}_i + \\eta \\; \\left( \\frac{1}{N} \\sum_{n=1}^N \\frac{y_n \\; {\\bf x}_n}{1 + e^{y_n {\\bf w}_i^T {\\bf x}_n}} \\right) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\eta$ is our learning rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Enough with the theory, now jump to the implimentation. We will look at 2 libraries for the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression with statsmodel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using the same dataset as UCLA's Logit Regression tutorial to explore logistic regression in Python. Our goal will be to identify the various factors that may influence admission into graduate school.\n",
    "\n",
    "The dataset contains several columns which we can use as predictor variables:\n",
    "\n",
    "   * gpa\n",
    "   * gre score\n",
    "   * rank or prestige of an applicant's undergraduate alma mater\n",
    "   * The fourth column, admit, is our binary target variable. It indicates whether or not a candidate was admitted our not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pylab as pl\n",
    "import statsmodels.api as sm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"binary.csv\")\n",
    "#df = pd.read_csv(\"https://stats.idre.ucla.edu/stat/data/binary.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>admit</th>\n",
       "      <th>gre</th>\n",
       "      <th>gpa</th>\n",
       "      <th>rank</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>380</td>\n",
       "      <td>3.61</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>660</td>\n",
       "      <td>3.67</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>800</td>\n",
       "      <td>4.00</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>640</td>\n",
       "      <td>3.19</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>520</td>\n",
       "      <td>2.93</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   admit  gre   gpa  rank\n",
       "0      0  380  3.61     3\n",
       "1      1  660  3.67     3\n",
       "2      1  800  4.00     1\n",
       "3      1  640  3.19     4\n",
       "4      0  520  2.93     4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>admit</th>\n",
       "      <th>gre</th>\n",
       "      <th>gpa</th>\n",
       "      <th>prestige</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>380</td>\n",
       "      <td>3.61</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>660</td>\n",
       "      <td>3.67</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>800</td>\n",
       "      <td>4.00</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>640</td>\n",
       "      <td>3.19</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>520</td>\n",
       "      <td>2.93</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   admit  gre   gpa  prestige\n",
       "0      0  380  3.61         3\n",
       "1      1  660  3.67         3\n",
       "2      1  800  4.00         1\n",
       "3      1  640  3.19         4\n",
       "4      0  520  2.93         4"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#rename the rank column b'oz there is also a dataframe method called \"rank\"\n",
    "\n",
    "\n",
    "df.columns = [\"admit\",\"gre\",\"gpa\",\"prestige\"]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(400, 4)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary Statistics & Looking at the data\n",
    "Now that we've got everything loaded into Python and named appropriately let's take a look at the data. We can use the pandas function which describes a summarized view of everything. There's also function for calculating the standard deviation, std.\n",
    "\n",
    "A feature I really like in pandas is the pivot_table/crosstab aggregations. crosstab makes it really easy to do multidimensional frequency tables. You might want to play around with this to look at different cuts of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>prestige</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>admit</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28</td>\n",
       "      <td>97</td>\n",
       "      <td>93</td>\n",
       "      <td>55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>33</td>\n",
       "      <td>54</td>\n",
       "      <td>28</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "prestige   1   2   3   4\n",
       "admit                   \n",
       "0         28  97  93  55\n",
       "1         33  54  28  12"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(df[\"admit\"],df[\"prestige\"],rownames=[\"admit\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<AxesSubplot:title={'center':'admit'}>,\n",
       "        <AxesSubplot:title={'center':'gre'}>],\n",
       "       [<AxesSubplot:title={'center':'gpa'}>,\n",
       "        <AxesSubplot:title={'center':'prestige'}>]], dtype=object)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAggElEQVR4nO3de7RcZZnn8e/PcDVBIQQDJpEDLY0GI7cMYDPdRGklXMZAD7rCyM1Low5pcTprScAewUHWoMsb4q2jQLCJQBRQGlChkaPNjIAE0QRChgBpOSQSlGtQwROe+WO/R4qTqnN2Xfeund9nrVpV9e69Tz21663n7Hr3u99XEYGZmVXXK4oOwMzMusuJ3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6EtG0oCkkLRVi9tvlLRnp+Mys/7VUjKx8oqISSOPJS0BhiLin4qLyMyK5iN6Myu1Vn/d2kuc6HtE0iJJD0p6VtJ9ko5L5RMkfVbSbyU9BBw9artBSZ+S9H9Ts8y/StpZ0lJJz0j6uaSBmvVD0uslnQa8B/jYyHa9fL9m45F0gKRfpO/EdyRdler6HElDks6U9BvgUkmvqPkO/U7SMkmTi34P/cKJvnceBP4aeDXwSeBySbsBfw8cA+wPzAaOr7PtfOAkYBrwF8DPgEuBycAq4JzRG0TEYmAp8JmImBQR/6XTb8isVZK2Aa4FlpDV4yuA42pW2TWV7w6cBnwEOBY4DHgt8CTwlZ4F3Oec6HskIr4TEesi4sWIuAp4ADgIeDfwxYh4JCKeAP53nc0vjYgHI+Jp4AfAgxHxbxExDHyH7J+EWT85hOwc4Zci4k8RcQ1wZ83yF4FzIuL5iPgD8EHg4xExFBHPA+cCx7tZJx/vpB6RdDLwj8BAKpoETCE7OnmkZtX/qLP5YzWP/1Dn+STM+strgUfj5aMq1n4PHo+IP9Y83x24VtKLNWWbgKnAo90Lsxp8RN8DknYHvgEsAHaOiB2BlYCA9cCMmtVf18GX9tCkVlbrgWmSVFNW+z0YXXcfAY6MiB1rbttFhJN8Dk70vTGRrOI+DiDpvcCb0rJlwEckTZe0E7Cog6/7GOA+9VZGPyM7Il8gaStJ88iaMhv5OnB+OmhC0i5pG8vBib4HIuI+4HNklfsxYBbwf9LibwA/An4J3A1c08GXvhiYKekpSd/r4N81a0tEvAD8HfB+4CngROB64PkGm1wIXAfcJOlZ4Hbg4O5HWg3yxCNmVgaS7gC+HhGXFh1L1fiI3swKIekwSbumpptTgDcDPyw6ripyrxszK8reZOeoJpFdZ3J8RKwvNqRqctONmVnFuenGzKziStF0M2XKlBgYGKi77LnnnmPixIm9DaiEvB8yY+2H5cuX/zYidulxSC0bqff9+Nn2Y8xQvbhz1/mIKPx24IEHRiO33nprw2VbEu+HzFj7AbgrSlCf895G6n0/frb9GHNE9eLOW+fddGNmVnFO9GZmFedEb2ZWcaU4GTuWFY8+zamLbmhqm7UXHD3+SmZWSQNj5IuFs4br5pOq5wwf0ZuZVZwTvZlZxTnRm5lVXOnb6M3KRNLewFU1RXsCnwB2JJv/9/FUfnZE3Njb6Mzqc6I3a0JErAb2A5A0gWwau2uB9wJfiIjPFhedWX1uujFr3eFkE7XXm+fXrDSc6M1aNx+4oub5Akm/knRJmhbSrBTcdGPWAknbAO8EzkpFXwPOI5sb+DyyqSPfV2e704DTAKZOncrg4CAbN25kcHCwF2F3TC9iXvHo0y1tt3BW42VTt8/60o9W9v3f7v52ojdrzZHA3RHxGMDIPYCkb5DNf7qZiFgMLAaYPXt2zJkzh8HBQebMmdP9iDuoFzE3e6FkHgtnDfO5FZunvbXvmdPx1+qkdve3m27MWnMCNc02knarWXYcsLLnEZk14CN6syZJeiXwduCDNcWfkbQfWdPN2lHLzArlRG/WpIj4PbDzqLKTCgrHbFxuujEzqzgnejOzinOiNzOrOCd6M7OKc6I3M6s4J3ozs4pz90oz2+KNNf1gI/00/aCP6M3MKm7cRC9phqRbJa2SdK+kM1L5ZEk3S3og3e9Us81ZktZIWi3piG6+ATMzG1ueI/phYGFEvBE4BDhd0kxgEXBLROwF3JKek5bNB/YB5gJfTRM0mJlZAcZN9BGxPiLuTo+fBVYB04B5wGVptcuAY9PjecCVEfF8RDwMrAEO6nDcZmaWU1MnYyUNAPsDdwBTI2I9ZP8MJL0mrTYNuL1ms6FUNvpvbTYudz2Nxo8eS9nHlm5FP45Z3g3eD2bNy53oJU0CrgY+GhHPSGq4ap2y2Kygzrjc9Vy09Pt1x48eS9nHlm5FP45Z3g3eD2bNy9XrRtLWZEl+aURck4ofGxmDO91vSOVDwIyazacD6zoTrpmZNStPrxsBFwOrIuLzNYuuA05Jj08Bvl9TPl/StpL2APYC7uxcyGZm1ow8bSKHAicBKyTdk8rOBi4Alkl6P/Br4F0AEXGvpGXAfWQ9dk6PiE2dDtzMzPIZN9FHxG3Ub3cHOLzBNucD57cRl5mZdYivjDUzqzgnejOzivOgZmZNkrQWeBbYBAxHxGxJk4GrgAGyycHfHRFPFhWjWS0f0Zu15q0RsV9EzE7P6w4JYlYGPqI364x5wJz0+DJgEDizqGDKppVhgK1znOjNmhfATZIC+Od0lXejIUFept7QH/04rEOzMTc7jEm3tDKkSiO9/MzarSNO9GbNOzQi1qVkfrOk+/NuWG/oj34c1qHZmE8tyRH9wlnDTQ+p0kgvh1ppt464jd6sSRGxLt1vAK4lG5210ZAgZoVzojdrgqSJknYYeQy8A1hJ4yFBzArnphuz5kwFrk2jt24FfDsifijp59QZEsSsDJzozZoQEQ8B+9Yp/x0NhgQxK5qbbszMKs5H9GZmLWjl2oC1FxzdhUjG5yN6M7OKc6I3M6s4J3ozs4pzojczqzgnejOzinOvGytMK70Wlsyd2IVIzKrNR/RmZhXnRG9mVnFuujEz65FWJ2Bpt8nSR/RmZhXnRG9mVnHjJnpJl0jaIGllTdlkSTdLeiDd71Sz7CxJayStlnREtwI3M7N88hzRLwHmjiqrO+O9pJnAfGCftM1XJU3oWLRmZta0cRN9RPwUeGJU8Tyyme5J98fWlF8ZEc9HxMPAGrJp1szMrCCt9rppNOP9NOD2mvWGUtlmJJ0GnAYwderUhjOctzJrey9nZ++VdmeBL6NmP1eo5n4w67ZOd69UnbKot2JELAYWA8yePTsazXB+0dLvNz1rey9nZ++VdmeBL6NTW7wytmr7wazbWk30j0naLR3N1854PwTMqFlvOrCunQDNykTSDOBbwK7Ai8DiiLhQ0rnA3wOPp1XPjogbi4myuwYW3cDCWcMt/aO2YrTavbLRjPfXAfMlbStpD2Av4M72QjQrlWFgYUS8ETgEOD11QgD4QkTsl26VTPLWn8Y9opd0BTAHmCJpCDgHuIA6M95HxL2SlgH3kX0hTo+ITV2K3azn0rmpkfNTz0paRYPzUGZlMW6ij4gTGiyqO+N9RJwPnN9OUGb9QNIAsD9wB3AosEDSycBdZEf9T9bZZrNOCP12gnnhrOGWOkmUQb/G3W4d8Vg3Zi2QNAm4GvhoRDwj6WvAeWSdD84DPge8b/R29Toh9NuJ9lNTG32znSTKoF/jbrcTgodAMGuSpK3JkvzSiLgGICIei4hNEfEi8A18/YiViBO9WRMkCbgYWBURn68p361mteOAlaO3NStK//2GMSvWocBJwApJ96Sys4ETJO1H1nSzFvhgEcGZ1eNEb9aEiLiN+hcGujullZabbszMKs6J3sys4pzozcwqzonezKzifDLWbAvW6mTV1l98RG9mVnFO9GZmFedEb2ZWcU70ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFedEb2ZWcU70ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFde1YYolzQUuBCYA34yIC7r1WmZlUGSd93DDNpauHNFLmgB8BTgSmAmcIGlmN17LrAxc563MunVEfxCwJiIeApB0JTAPuK9Lr2dWtI7VeR+dW6d1K9FPAx6peT4EHFy7gqTTgNPS042SVjf4W1OA3zbz4vp0M2v3jab3QxW99dNj7ofdexnLKOPWeWhY7/vus/1IH8YM/Rv3GPU+V53vVqJXnbJ42ZOIxcDicf+QdFdEzO5UYP3K+yFT4v0wbp2H+vW+xO+poX6MGbbcuLvV62YImFHzfDqwrkuvZVYGrvNWWt1K9D8H9pK0h6RtgPnAdV16LbMycJ230upK001EDEtaAPyIrKvZJRFxb4t/btzmnS2E90OmlPuhzTpfyvc0jn6MGbbQuBWxWTOimZlViK+MNTOrOCd6M7OKK02ilzRX0mpJayQtqrNckr6Ulv9K0gFFxNltOfbDHElPS7on3T5RRJzdJOkSSRskrWywvG/qgqQZkm6VtErSvZLOSOWTJd0s6YF0v1PNNmel97Za0hEFxj5B0i8kXd9HMe8o6buS7k/7/C19Evf/SPVjpaQrJG3X0bgjovAb2cmrB4E9gW2AXwIzR61zFPADsv7KhwB3FB13QfthDnB90bF2eT/8DXAAsLLB8r6pC8BuwAHp8Q7A/yMbIuEzwKJUvgj4dHo8M33u2wJ7pPowoaDY/xH49kh965OYLwM+kB5vA+xY9rjJLrZ7GNg+PV8GnNrJuMtyRP/ny8cj4gVg5PLxWvOAb0XmdmBHSbv1OtAuy7MfKi8ifgo8McYqfVMXImJ9RNydHj8LrCL7Ys8jS0qk+2PT43nAlRHxfEQ8DKwhqxc9JWk6cDTwzZrissf8KrKDhIsBIuKFiHiKksedbAVsL2kr4JVk12B0LO6yJPp6l49Pa2Gdfpf3Pb5F0i8l/UDSPr0JrVT6si5IGgD2B+4ApkbEesj+GQCvSauV5b19EfgY8GJNWdlj3hN4HLg0NTl9U9JESh53RDwKfBb4NbAeeDoibqKDcZcl0ee5fDzXJeZ9Ls97vBvYPSL2BS4CvtftoEqo7+qCpEnA1cBHI+KZsVatU9bT9ybpGGBDRCzPu0mdsiI+j63Imvy+FhH7A8+RNXk0Uoq4U9v7PLJmmNcCEyWdONYmdcrGjLssiT7P5eNbwiXm477HiHgmIjamxzcCW0ua0rsQS6Gv6oKkrcmS/NKIuCYVPzbS3JTuN6TyMry3Q4F3SlpL1nz4NkmXU6KYJW2UtOeo4iFgKCLuSM+/S5b4SxN3A38LPBwRj0fEn4BrgL+ig3GXJdHnuXz8OuDk1OPiELKfN+t7HWiXjbsfJO0qSenxQWSf4e96Hmmx+qYupM/qYmBVRHy+ZtF1wCnp8SnA92vK50vaVtIewF7Anb2KFyAizoqI6RExQFYHfxwRJxYVs6RBSR8YFeOkSENC15T9BnhE0t6p6HCyYaJLu6+TXwOHSHplqi+Hk53L6VzcvT7DPMaZ56PIeiQ8CHw8lX0I+FB6LLKJHR4EVgCzi465oP2wALiX7Kz77cBfFR1zF/bBFWRtlX8iO3p5f7/WBeA/k/2s/hVwT7odBewM3AI8kO4n12zz8fTeVgNHFhz/HF7qddNyzMBWbcQwSOpJk2Pd/YC70v7+HrBTP+xr4JPA/cBK4F/IetR0LO7Cvwhb+o3sp+UvgGeB7wBXAZ9KX7Ah4GyycajXAu+p2e7otN0zZCdmzi36vfi25d1SvTyL7Mj5SeBSYLua+nsm8JuUvF5B1mb+INmv0GUjySttc3kqf4rs1+1U4HxgE/BHYCPw5bR+AK9Pj3cG/jV9F36evj+31cT4BuBmsp5cq4F3F73fen0rS9PNFik1z1wLLAEmkx3JHlezyq5kEyVMI/vptrjmZ+lzwMlk/YSPBj4s6dhexG02ynuAI4C/AP4S+KdUvitZvd6dbLKVj5B1ETyM7KTjk2S/zCCr368ma3vemewX3B8i4uPAvwMLImuuWVDn9b9C9n3YNf2dkeYOUq+bm8muB3gNcALw1S2tt5oTfbEOIesp8KWI+FNkJ+pGt7X9z8j6y/4EuAF4N0BEDEbEioh4MSJ+RfZP4rBeBm+WfDkiHomIJ8iOwE9I5S8C56T6+wfgg2TNkUMR8TxwLnB86jv+J7IE//qI2BQRy2Ps3knAn+fq/a/pdX4fEffxUt9zgGOAtRFxaUQMR3ZNw9XA8R15532iWzNMWT6vBR6N9Psyqe0f+2REPFfz/D/SNkg6GLgAeBPZFYDbkjX9mPVabZ39cx0FHo+IP9Ys2x24VlJt3/xNZE00/0J2NH+lpB3JmnE+HlkvlLHsQpbHamOofbw7cLCkp2rKtkqvt8XwEX2x1gPTRnrRJLXdpnZKPz1HvI6XulF9m+zs+4yIeDXwder3rzXrtto6W1tHR/ftfoTsxOGONbftIuLR9Iv2kxExk6xr4TFkTZP1/k6tx4Fhsi6G9eJ5BPjJqNecFBEfbvI99jUn+mL9jOyIZoGkrSTNY/NLmT8paRtJf01W+UeO2ncAnoiIP6Zulv+tZ1GbvdzpkqZLmkzWeeCqBut9HThf0u4AknZJdR5Jb5U0KzXFPEPWlLMpbfcY2VWvm4mITWT9zs9N3RPfwEv/IACuB/5S0kmStk63/yTpje295f7iRF+gyMaz+Tuy7oNPASeSVczn0yq/ITthtQ5YSta98P607L8D/0vSs8AnyHowmBXh28BNwEPp9qkG611I9iv0plRvbwcOTst2JbvA6RmyPuQ/IWu+GdnueElPSvpSnb+7gOxE7kjvnitI36HIxhd6B9n1AOvSOp8ma+rcYniGqZKRdAfZkc/DwOURMX2cTcwKk66e/UBE/FvRsYyQ9Glg14g4ZdyVtxA+oi+YpMPS1a5bSToFeDPww6LjMusXkt4g6c3pSumDyH4hX1t0XGXiXjfF25us2WUS2YUkx0fE+pr+8mY2th3ImmteSzYezOd4abgAw003ZmaV56YbM7OKK0XTzZQpU2JgYKCrr/Hcc88xceLE8VfsMcfVnLHiWr58+W8jYpceh9SyRvW+rPu+EcfbfY1izl3nix5sJyI48MADo9tuvfXWrr9GKxxXc8aKC7grSlCf894a1fuy7vtGHG/3NYo5b513041ZHZIukbRB0sqasnMlPSrpnnQ7qmbZWZLWSFot6Yhiojarz4nerL4lwNw65V+IiP3S7UYASTPJLsjZJ23z1XSFp1kpONGb1RERPyUbvzyPecCVkY3S+DCwhs2HsjArTClOxlp3DCy6oaXt1l5wdIcjqZQFkk4mm8VoYUQ8STZfwO016wylss1IOo1sbHamTp3K4ODgZuts3LixbnlZrHj06Zc9n7o9XLR07G7rs6a9upshNaXs+7eedmN2ojfL72vAeWSjKZ5HdmHO+6g/amjdC1QiYjGwGGD27NkxZ86czdYZHBykXnlZnDrqAGLhrGE+t2LsVLL2PXO6GFFzyr5/62k3ZjfdmOUUEY9FNinGi8A3eKl5ZoiXD407nZeG6jUrnBO9WU6Sdqt5ehzZRM6Qjcg4X9K2kvYA9mLzmcLMCuOmG7M6JF1BNsH1FElDwDnAHEn7kTXLrCWbGo+IuFfSMrIJsoeB0yMbJ92sFJzozeqIiBPqFF88xvrnk82XalY6broxM6s4J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKc/dK28zAohtYOGt4s0vdx+LxcczKy0f0ZmYV50RvZlZxbrqxwrQyjPKSuf0116dZGbR1RC9pR0nflXS/pFWS3iJpsqSbJT2Q7nfqVLBmZta8dptuLgR+GBFvAPYFVgGLgFsiYi/glvTczMwK0nLTjaRXAX8DnAoQES8AL0iaRzbqH8BlwCBwZjtBmpnlNV6TYL0eZVXvNdZOG/2ewOPApZL2BZYDZwBTI2I9QESsl/SaehvnmVKtk8o6fVg341o4a7jlbadu39z2rbyHVuIr6+doVmbtJPqtgAOAf4iIOyRdSBPNNHmmVOuksk4f1s24mukHP1qe6eFqtTJVXCvxLZk7sZSfo1mZtdNGPwQMRcQd6fl3yRL/YyMz8aT7De2FaGZm7Wj5iD4ifiPpEUl7R8Rq4HCyGXbuA04BLkj3Y08Pb7m00hXRzAza70f/D8BSSdsADwHvJfuVsEzS+4FfA+9q8zXMzKwNbSX6iLgHmF1n0eHt/F0zM+scD4FgZlZxTvRmZhXnRG9Wh6RLJG2QtLKmrOHwHpLOkrRG0mpJRxQTtVl9TvRm9S0B5o4qqzu8h6SZwHxgn7TNVyVN6F2oZmNzojerIyJ+Cjwxqnge2bAepPtja8qvjIjnI+JhYA1wUC/iNMvDwxSb5ddoeI9pwO016w2lss3kGfqj7MM8jB66Is9wGb18P+PFUi/eMu9vaL9OONGbtU91yqLeinmG/ijrcB0jRg9dkWe4jFaGyGjVeENr1Iu3l/G1ot064aYbs/waDe8xBMyoWW86sK7HsZk15ERvlt91ZMN6wMuH97gOmC9pW0l7AHsBdxYQn1ldbroxq0PSFWTzKkyRNAScQzZ+02bDe0TEvZKWkY3zNAycHhGbCgncrA4nerM6IuKEBovqDu8REecD53cvIrPWuenGzKzinOjNzCrOid7MrOLaTvSSJkj6haTr0/OG44GYmVnvdeJk7BnAKuBV6fnIeCAXSFqUnp/ZgdepjNrZourNSG9m1kltHdFLmg4cDXyzprjReCBmZlaAdo/ovwh8DNihpqzReCAvk2fMj04q0/ghteNs5BknpAjNxtXKvm3lfZfpczTrFy0neknHABsiYrmkOc1un2fMj04q0/ghp45quhlvnJAiNBtXK2OFtNJktWTuxNJ8jmb9op0McyjwTklHAdsBr5J0OWk8kHQ0XzseiFXYgM8zmJVWy230EXFWREyPiAGySRd+HBEn0ng8EDMzK0A3+tFfALxd0gPA29NzMzMrSEcahyNiEBhMj39Hg/FAzMys93xlrJlZxTnRm5lVXPn69fUZ9zYxs7xazRdL5k5s63V9RG9mVnFO9GZmFeemG7OSWfHo001fNbz2gqO7FI1VgY/ozcwqzonezKzinOjNzCrObfRmTZK0FngW2AQMR8RsSZOBq4ABYC3w7oh4sqgYzWr5iN6sNW+NiP0iYnZ6PjKz2l7ALem5WSk40Zt1hmdWs9Jy041Z8wK4SVIA/5wm0enYzGqtzDrWy1m3RseWJ94i4xutXry9iq/V2eTanVnNid6seYdGxLqUzG+WdH/eDfPMrHbR0u83PetYKzN8tWp0H/88s5EVGd9o9eLtVXytzKoG7c+s5qYbsyZFxLp0vwG4FjiINLMagGdWs7JpOdFLmiHpVkmrJN0r6YxUPlnSzZIeSPc7dS5cs2JJmihph5HHwDuAlXhmNSuxdo7oh4GFEfFG4BDgdEkzce8Dq7apwG2SfgncCdwQET/EM6tZibXcRp9OPI2cfHpW0ipgGlnvgzlptcvIZp46s60ozUoiIh4C9q1T7pnVrLQ6cjJW0gCwP3AHHex90EntnrVupNWz6CNa6WHRC2WNq1ufo1mVtZ3oJU0CrgY+GhHPSMq1XZ7eB500ODjY1lnrRlo9iz4iT4+FIpQ1rnZ7H5htidrqdSNpa7IkvzQirknF7n1gZlYi7fS6EXAxsCoiPl+zyL0PzMxKpJ3f5ocCJwErJN2Tys4m622wTNL7gV8D72orQjMza0s7vW5uAxo1yPdd7wNP8m1mVeUrY83MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKq58Uwh1QL2RKBfOGm57Nigzs35U+kTv4YPNzNrTtaYbSXMlrZa0RtKibr2OWVm4zltZdSXRS5oAfAU4EpgJnCBpZjdey6wMXOetzLp1RH8QsCYiHoqIF4ArgXldei2zMnCdt9JSRHT+j0rHA3Mj4gPp+UnAwRGxoGad04DT0tO9gdUdD+TlpgC/7fJrtMJxNWesuHaPiF16GcyIPHU+leep92Xd94043u5rFHOuOt+tk7H15pJ92X+UiFgMLO7S629G0l0RMbtXr5eX42pOWeMiR52HfPW+xO+xLsfbfe3G3K2mmyFgRs3z6cC6Lr2WWRm4zltpdSvR/xzYS9IekrYB5gPXdem1zMrAdd5KqytNNxExLGkB8CNgAnBJRNzbjddqQs+aiZrkuJpTyrg6XOdL+R7H4Hi7r62Yu3Iy1szMysNj3ZiZVZwTvZlZxVUq0UuaIelWSask3SvpjDrrzJH0tKR70u0TPYhrO0l3SvpliuuTddaRpC+ly+d/JemAksTV8/1V89oTJP1C0vV1lvV8f3WbpEskbZC0suhY8sjzfSuTPPW9jMb6HuRV+kHNmjQMLIyIuyXtACyXdHNE3DdqvX+PiGN6GNfzwNsiYqOkrYHbJP0gIm6vWedIYK90Oxj4WrovOi7o/f4acQawCnhVnWVF7K9uWwJ8GfhWwXHklff7VhZ563vZjPU9yKVSR/QRsT4i7k6PnyXbOdOKjQoiszE93TrdRp8Fnwd8K617O7CjpN1KEFchJE0Hjga+2WCVnu+vbouInwJPFB1HXmX9vjVS5vreSI7vQS6VSvS1JA0A+wN31Fn8lvTz7QeS9ulRPBMk3QNsAG6OiNFxTQMeqXk+RA++NDniggL2F/BF4GPAiw2WF7K/rL5xvm+lkbO+l8kXGft7kEslE72kScDVwEcj4plRi+8mGx9iX+Ai4Hu9iCkiNkXEfmRXTB4k6U2jVsl1CX0BcfV8f0k6BtgQEcvHWq1OWamPzqpqnO9bqeSo76WR83uQS+USfWp7uxpYGhHXjF4eEc+M/HyLiBuBrSVN6VV8EfEUMAjMHbWo0EvoG8VV0P46FHinpLVko0C+TdLlo9bxkAMlMN73razG+B6WSZ7vQS6VSvSSBFwMrIqIzzdYZ9e0HpIOItsHv+tyXLtI2jE93h74W+D+UatdB5ycepMcAjwdEeuLjquI/RURZ0XE9IgYIBtK4McRceKo1Xq+v+zl8nzfyiTn97A0cn4Pcqlar5tDgZOAFakdDuBs4HUAEfF14Hjgw5KGgT8A86P7lwfvBlymbHKKVwDLIuJ6SR+qietG4ChgDfB74L1djilvXEXsr7pKsL+6StIVwBxgiqQh4JyIuLjYqMZU9/uWfvmVUd36XnBMPeEhEMzMKq5STTdmZrY5J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKc6I3M6u4/w9OeDjGqgF6ZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAggElEQVR4nO3de7RcZZnn8e/PcDVBIQQDJpEDLY0GI7cMYDPdRGklXMZAD7rCyM1Low5pcTprScAewUHWoMsb4q2jQLCJQBRQGlChkaPNjIAE0QRChgBpOSQSlGtQwROe+WO/R4qTqnN2Xfeund9nrVpV9e69Tz21663n7Hr3u99XEYGZmVXXK4oOwMzMusuJ3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6EtG0oCkkLRVi9tvlLRnp+Mys/7VUjKx8oqISSOPJS0BhiLin4qLyMyK5iN6Myu1Vn/d2kuc6HtE0iJJD0p6VtJ9ko5L5RMkfVbSbyU9BBw9artBSZ+S9H9Ts8y/StpZ0lJJz0j6uaSBmvVD0uslnQa8B/jYyHa9fL9m45F0gKRfpO/EdyRdler6HElDks6U9BvgUkmvqPkO/U7SMkmTi34P/cKJvnceBP4aeDXwSeBySbsBfw8cA+wPzAaOr7PtfOAkYBrwF8DPgEuBycAq4JzRG0TEYmAp8JmImBQR/6XTb8isVZK2Aa4FlpDV4yuA42pW2TWV7w6cBnwEOBY4DHgt8CTwlZ4F3Oec6HskIr4TEesi4sWIuAp4ADgIeDfwxYh4JCKeAP53nc0vjYgHI+Jp4AfAgxHxbxExDHyH7J+EWT85hOwc4Zci4k8RcQ1wZ83yF4FzIuL5iPgD8EHg4xExFBHPA+cCx7tZJx/vpB6RdDLwj8BAKpoETCE7OnmkZtX/qLP5YzWP/1Dn+STM+strgUfj5aMq1n4PHo+IP9Y83x24VtKLNWWbgKnAo90Lsxp8RN8DknYHvgEsAHaOiB2BlYCA9cCMmtVf18GX9tCkVlbrgWmSVFNW+z0YXXcfAY6MiB1rbttFhJN8Dk70vTGRrOI+DiDpvcCb0rJlwEckTZe0E7Cog6/7GOA+9VZGPyM7Il8gaStJ88iaMhv5OnB+OmhC0i5pG8vBib4HIuI+4HNklfsxYBbwf9LibwA/An4J3A1c08GXvhiYKekpSd/r4N81a0tEvAD8HfB+4CngROB64PkGm1wIXAfcJOlZ4Hbg4O5HWg3yxCNmVgaS7gC+HhGXFh1L1fiI3swKIekwSbumpptTgDcDPyw6ripyrxszK8reZOeoJpFdZ3J8RKwvNqRqctONmVnFuenGzKziStF0M2XKlBgYGKi77LnnnmPixIm9DaiEvB8yY+2H5cuX/zYidulxSC0bqff9+Nn2Y8xQvbhz1/mIKPx24IEHRiO33nprw2VbEu+HzFj7AbgrSlCf895G6n0/frb9GHNE9eLOW+fddGNmVnFO9GZmFedEb2ZWcaU4GTuWFY8+zamLbmhqm7UXHD3+SmZWSQNj5IuFs4br5pOq5wwf0ZuZVZwTvZlZxTnRm5lVXOnb6M3KRNLewFU1RXsCnwB2JJv/9/FUfnZE3Njb6Mzqc6I3a0JErAb2A5A0gWwau2uB9wJfiIjPFhedWX1uujFr3eFkE7XXm+fXrDSc6M1aNx+4oub5Akm/knRJmhbSrBTcdGPWAknbAO8EzkpFXwPOI5sb+DyyqSPfV2e704DTAKZOncrg4CAbN25kcHCwF2F3TC9iXvHo0y1tt3BW42VTt8/60o9W9v3f7v52ojdrzZHA3RHxGMDIPYCkb5DNf7qZiFgMLAaYPXt2zJkzh8HBQebMmdP9iDuoFzE3e6FkHgtnDfO5FZunvbXvmdPx1+qkdve3m27MWnMCNc02knarWXYcsLLnEZk14CN6syZJeiXwduCDNcWfkbQfWdPN2lHLzArlRG/WpIj4PbDzqLKTCgrHbFxuujEzqzgnejOzinOiNzOrOCd6M7OKc6I3M6s4J3ozs4pz90oz2+KNNf1gI/00/aCP6M3MKm7cRC9phqRbJa2SdK+kM1L5ZEk3S3og3e9Us81ZktZIWi3piG6+ATMzG1ueI/phYGFEvBE4BDhd0kxgEXBLROwF3JKek5bNB/YB5gJfTRM0mJlZAcZN9BGxPiLuTo+fBVYB04B5wGVptcuAY9PjecCVEfF8RDwMrAEO6nDcZmaWU1MnYyUNAPsDdwBTI2I9ZP8MJL0mrTYNuL1ms6FUNvpvbTYudz2Nxo8eS9nHlm5FP45Z3g3eD2bNy53oJU0CrgY+GhHPSGq4ap2y2Kygzrjc9Vy09Pt1x48eS9nHlm5FP45Z3g3eD2bNy9XrRtLWZEl+aURck4ofGxmDO91vSOVDwIyazacD6zoTrpmZNStPrxsBFwOrIuLzNYuuA05Jj08Bvl9TPl/StpL2APYC7uxcyGZm1ow8bSKHAicBKyTdk8rOBi4Alkl6P/Br4F0AEXGvpGXAfWQ9dk6PiE2dDtzMzPIZN9FHxG3Ub3cHOLzBNucD57cRl5mZdYivjDUzqzgnejOzivOgZmZNkrQWeBbYBAxHxGxJk4GrgAGyycHfHRFPFhWjWS0f0Zu15q0RsV9EzE7P6w4JYlYGPqI364x5wJz0+DJgEDizqGDKppVhgK1znOjNmhfATZIC+Od0lXejIUFept7QH/04rEOzMTc7jEm3tDKkSiO9/MzarSNO9GbNOzQi1qVkfrOk+/NuWG/oj34c1qHZmE8tyRH9wlnDTQ+p0kgvh1ppt464jd6sSRGxLt1vAK4lG5210ZAgZoVzojdrgqSJknYYeQy8A1hJ4yFBzArnphuz5kwFrk2jt24FfDsifijp59QZEsSsDJzozZoQEQ8B+9Yp/x0NhgQxK5qbbszMKs5H9GZmLWjl2oC1FxzdhUjG5yN6M7OKc6I3M6s4J3ozs4pzojczqzgnejOzinOvGytMK70Wlsyd2IVIzKrNR/RmZhXnRG9mVnFuujEz65FWJ2Bpt8nSR/RmZhXnRG9mVnHjJnpJl0jaIGllTdlkSTdLeiDd71Sz7CxJayStlnREtwI3M7N88hzRLwHmjiqrO+O9pJnAfGCftM1XJU3oWLRmZta0cRN9RPwUeGJU8Tyyme5J98fWlF8ZEc9HxMPAGrJp1szMrCCt9rppNOP9NOD2mvWGUtlmJJ0GnAYwderUhjOctzJrey9nZ++VdmeBL6NmP1eo5n4w67ZOd69UnbKot2JELAYWA8yePTsazXB+0dLvNz1rey9nZ++VdmeBL6NTW7wytmr7wazbWk30j0naLR3N1854PwTMqFlvOrCunQDNykTSDOBbwK7Ai8DiiLhQ0rnA3wOPp1XPjogbi4myuwYW3cDCWcMt/aO2YrTavbLRjPfXAfMlbStpD2Av4M72QjQrlWFgYUS8ETgEOD11QgD4QkTsl26VTPLWn8Y9opd0BTAHmCJpCDgHuIA6M95HxL2SlgH3kX0hTo+ITV2K3azn0rmpkfNTz0paRYPzUGZlMW6ij4gTGiyqO+N9RJwPnN9OUGb9QNIAsD9wB3AosEDSycBdZEf9T9bZZrNOCP12gnnhrOGWOkmUQb/G3W4d8Vg3Zi2QNAm4GvhoRDwj6WvAeWSdD84DPge8b/R29Toh9NuJ9lNTG32znSTKoF/jbrcTgodAMGuSpK3JkvzSiLgGICIei4hNEfEi8A18/YiViBO9WRMkCbgYWBURn68p361mteOAlaO3NStK//2GMSvWocBJwApJ96Sys4ETJO1H1nSzFvhgEcGZ1eNEb9aEiLiN+hcGujullZabbszMKs6J3sys4pzozcwqzonezKzifDLWbAvW6mTV1l98RG9mVnFO9GZmFedEb2ZWcU70ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFedEb2ZWcU70ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFde1YYolzQUuBCYA34yIC7r1WmZlUGSd93DDNpauHNFLmgB8BTgSmAmcIGlmN17LrAxc563MunVEfxCwJiIeApB0JTAPuK9Lr2dWtI7VeR+dW6d1K9FPAx6peT4EHFy7gqTTgNPS042SVjf4W1OA3zbz4vp0M2v3jab3QxW99dNj7ofdexnLKOPWeWhY7/vus/1IH8YM/Rv3GPU+V53vVqJXnbJ42ZOIxcDicf+QdFdEzO5UYP3K+yFT4v0wbp2H+vW+xO+poX6MGbbcuLvV62YImFHzfDqwrkuvZVYGrvNWWt1K9D8H9pK0h6RtgPnAdV16LbMycJ230upK001EDEtaAPyIrKvZJRFxb4t/btzmnS2E90OmlPuhzTpfyvc0jn6MGbbQuBWxWTOimZlViK+MNTOrOCd6M7OKK02ilzRX0mpJayQtqrNckr6Ulv9K0gFFxNltOfbDHElPS7on3T5RRJzdJOkSSRskrWywvG/qgqQZkm6VtErSvZLOSOWTJd0s6YF0v1PNNmel97Za0hEFxj5B0i8kXd9HMe8o6buS7k/7/C19Evf/SPVjpaQrJG3X0bgjovAb2cmrB4E9gW2AXwIzR61zFPADsv7KhwB3FB13QfthDnB90bF2eT/8DXAAsLLB8r6pC8BuwAHp8Q7A/yMbIuEzwKJUvgj4dHo8M33u2wJ7pPowoaDY/xH49kh965OYLwM+kB5vA+xY9rjJLrZ7GNg+PV8GnNrJuMtyRP/ny8cj4gVg5PLxWvOAb0XmdmBHSbv1OtAuy7MfKi8ifgo8McYqfVMXImJ9RNydHj8LrCL7Ys8jS0qk+2PT43nAlRHxfEQ8DKwhqxc9JWk6cDTwzZrissf8KrKDhIsBIuKFiHiKksedbAVsL2kr4JVk12B0LO6yJPp6l49Pa2Gdfpf3Pb5F0i8l/UDSPr0JrVT6si5IGgD2B+4ApkbEesj+GQCvSauV5b19EfgY8GJNWdlj3hN4HLg0NTl9U9JESh53RDwKfBb4NbAeeDoibqKDcZcl0ee5fDzXJeZ9Ls97vBvYPSL2BS4CvtftoEqo7+qCpEnA1cBHI+KZsVatU9bT9ybpGGBDRCzPu0mdsiI+j63Imvy+FhH7A8+RNXk0Uoq4U9v7PLJmmNcCEyWdONYmdcrGjLssiT7P5eNbwiXm477HiHgmIjamxzcCW0ua0rsQS6Gv6oKkrcmS/NKIuCYVPzbS3JTuN6TyMry3Q4F3SlpL1nz4NkmXU6KYJW2UtOeo4iFgKCLuSM+/S5b4SxN3A38LPBwRj0fEn4BrgL+ig3GXJdHnuXz8OuDk1OPiELKfN+t7HWiXjbsfJO0qSenxQWSf4e96Hmmx+qYupM/qYmBVRHy+ZtF1wCnp8SnA92vK50vaVtIewF7Anb2KFyAizoqI6RExQFYHfxwRJxYVs6RBSR8YFeOkSENC15T9BnhE0t6p6HCyYaJLu6+TXwOHSHplqi+Hk53L6VzcvT7DPMaZ56PIeiQ8CHw8lX0I+FB6LLKJHR4EVgCzi465oP2wALiX7Kz77cBfFR1zF/bBFWRtlX8iO3p5f7/WBeA/k/2s/hVwT7odBewM3AI8kO4n12zz8fTeVgNHFhz/HF7qddNyzMBWbcQwSOpJk2Pd/YC70v7+HrBTP+xr4JPA/cBK4F/IetR0LO7Cvwhb+o3sp+UvgGeB7wBXAZ9KX7Ah4GyycajXAu+p2e7otN0zZCdmzi36vfi25d1SvTyL7Mj5SeBSYLua+nsm8JuUvF5B1mb+INmv0GUjySttc3kqf4rs1+1U4HxgE/BHYCPw5bR+AK9Pj3cG/jV9F36evj+31cT4BuBmsp5cq4F3F73fen0rS9PNFik1z1wLLAEmkx3JHlezyq5kEyVMI/vptrjmZ+lzwMlk/YSPBj4s6dhexG02ynuAI4C/AP4S+KdUvitZvd6dbLKVj5B1ETyM7KTjk2S/zCCr368ma3vemewX3B8i4uPAvwMLImuuWVDn9b9C9n3YNf2dkeYOUq+bm8muB3gNcALw1S2tt5oTfbEOIesp8KWI+FNkJ+pGt7X9z8j6y/4EuAF4N0BEDEbEioh4MSJ+RfZP4rBeBm+WfDkiHomIJ8iOwE9I5S8C56T6+wfgg2TNkUMR8TxwLnB86jv+J7IE//qI2BQRy2Ps3knAn+fq/a/pdX4fEffxUt9zgGOAtRFxaUQMR3ZNw9XA8R15532iWzNMWT6vBR6N9Psyqe0f+2REPFfz/D/SNkg6GLgAeBPZFYDbkjX9mPVabZ39cx0FHo+IP9Ys2x24VlJt3/xNZE00/0J2NH+lpB3JmnE+HlkvlLHsQpbHamOofbw7cLCkp2rKtkqvt8XwEX2x1gPTRnrRJLXdpnZKPz1HvI6XulF9m+zs+4yIeDXwder3rzXrtto6W1tHR/ftfoTsxOGONbftIuLR9Iv2kxExk6xr4TFkTZP1/k6tx4Fhsi6G9eJ5BPjJqNecFBEfbvI99jUn+mL9jOyIZoGkrSTNY/NLmT8paRtJf01W+UeO2ncAnoiIP6Zulv+tZ1GbvdzpkqZLmkzWeeCqBut9HThf0u4AknZJdR5Jb5U0KzXFPEPWlLMpbfcY2VWvm4mITWT9zs9N3RPfwEv/IACuB/5S0kmStk63/yTpje295f7iRF+gyMaz+Tuy7oNPASeSVczn0yq/ITthtQ5YSta98P607L8D/0vSs8AnyHowmBXh28BNwEPp9qkG611I9iv0plRvbwcOTst2JbvA6RmyPuQ/IWu+GdnueElPSvpSnb+7gOxE7kjvnitI36HIxhd6B9n1AOvSOp8ma+rcYniGqZKRdAfZkc/DwOURMX2cTcwKk66e/UBE/FvRsYyQ9Glg14g4ZdyVtxA+oi+YpMPS1a5bSToFeDPww6LjMusXkt4g6c3pSumDyH4hX1t0XGXiXjfF25us2WUS2YUkx0fE+pr+8mY2th3ImmteSzYezOd4abgAw003ZmaV56YbM7OKK0XTzZQpU2JgYKCrr/Hcc88xceLE8VfsMcfVnLHiWr58+W8jYpceh9SyRvW+rPu+EcfbfY1izl3nix5sJyI48MADo9tuvfXWrr9GKxxXc8aKC7grSlCf894a1fuy7vtGHG/3NYo5b513041ZHZIukbRB0sqasnMlPSrpnnQ7qmbZWZLWSFot6Yhiojarz4nerL4lwNw65V+IiP3S7UYASTPJLsjZJ23z1XSFp1kpONGb1RERPyUbvzyPecCVkY3S+DCwhs2HsjArTClOxlp3DCy6oaXt1l5wdIcjqZQFkk4mm8VoYUQ8STZfwO016wylss1IOo1sbHamTp3K4ODgZuts3LixbnlZrHj06Zc9n7o9XLR07G7rs6a9upshNaXs+7eedmN2ojfL72vAeWSjKZ5HdmHO+6g/amjdC1QiYjGwGGD27NkxZ86czdYZHBykXnlZnDrqAGLhrGE+t2LsVLL2PXO6GFFzyr5/62k3ZjfdmOUUEY9FNinGi8A3eKl5ZoiXD407nZeG6jUrnBO9WU6Sdqt5ehzZRM6Qjcg4X9K2kvYA9mLzmcLMCuOmG7M6JF1BNsH1FElDwDnAHEn7kTXLrCWbGo+IuFfSMrIJsoeB0yMbJ92sFJzozeqIiBPqFF88xvrnk82XalY6broxM6s4J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKc/dK28zAohtYOGt4s0vdx+LxcczKy0f0ZmYV50RvZlZxbrqxwrQyjPKSuf0116dZGbR1RC9pR0nflXS/pFWS3iJpsqSbJT2Q7nfqVLBmZta8dptuLgR+GBFvAPYFVgGLgFsiYi/glvTczMwK0nLTjaRXAX8DnAoQES8AL0iaRzbqH8BlwCBwZjtBmpnlNV6TYL0eZVXvNdZOG/2ewOPApZL2BZYDZwBTI2I9QESsl/SaehvnmVKtk8o6fVg341o4a7jlbadu39z2rbyHVuIr6+doVmbtJPqtgAOAf4iIOyRdSBPNNHmmVOuksk4f1s24mukHP1qe6eFqtTJVXCvxLZk7sZSfo1mZtdNGPwQMRcQd6fl3yRL/YyMz8aT7De2FaGZm7Wj5iD4ifiPpEUl7R8Rq4HCyGXbuA04BLkj3Y08Pb7m00hXRzAza70f/D8BSSdsADwHvJfuVsEzS+4FfA+9q8zXMzKwNbSX6iLgHmF1n0eHt/F0zM+scD4FgZlZxTvRmZhXnRG9Wh6RLJG2QtLKmrOHwHpLOkrRG0mpJRxQTtVl9TvRm9S0B5o4qqzu8h6SZwHxgn7TNVyVN6F2oZmNzojerIyJ+Cjwxqnge2bAepPtja8qvjIjnI+JhYA1wUC/iNMvDwxSb5ddoeI9pwO016w2lss3kGfqj7MM8jB66Is9wGb18P+PFUi/eMu9vaL9OONGbtU91yqLeinmG/ijrcB0jRg9dkWe4jFaGyGjVeENr1Iu3l/G1ot064aYbs/waDe8xBMyoWW86sK7HsZk15ERvlt91ZMN6wMuH97gOmC9pW0l7AHsBdxYQn1ldbroxq0PSFWTzKkyRNAScQzZ+02bDe0TEvZKWkY3zNAycHhGbCgncrA4nerM6IuKEBovqDu8REecD53cvIrPWuenGzKzinOjNzCrOid7MrOLaTvSSJkj6haTr0/OG44GYmVnvdeJk7BnAKuBV6fnIeCAXSFqUnp/ZgdepjNrZourNSG9m1kltHdFLmg4cDXyzprjReCBmZlaAdo/ovwh8DNihpqzReCAvk2fMj04q0/ghteNs5BknpAjNxtXKvm3lfZfpczTrFy0neknHABsiYrmkOc1un2fMj04q0/ghp45quhlvnJAiNBtXK2OFtNJktWTuxNJ8jmb9op0McyjwTklHAdsBr5J0OWk8kHQ0XzseiFXYgM8zmJVWy230EXFWREyPiAGySRd+HBEn0ng8EDMzK0A3+tFfALxd0gPA29NzMzMrSEcahyNiEBhMj39Hg/FAzMys93xlrJlZxTnRm5lVXPn69fUZ9zYxs7xazRdL5k5s63V9RG9mVnFO9GZmFeemG7OSWfHo001fNbz2gqO7FI1VgY/ozcwqzonezKzinOjNzCrObfRmTZK0FngW2AQMR8RsSZOBq4ABYC3w7oh4sqgYzWr5iN6sNW+NiP0iYnZ6PjKz2l7ALem5WSk40Zt1hmdWs9Jy041Z8wK4SVIA/5wm0enYzGqtzDrWy1m3RseWJ94i4xutXry9iq/V2eTanVnNid6seYdGxLqUzG+WdH/eDfPMrHbR0u83PetYKzN8tWp0H/88s5EVGd9o9eLtVXytzKoG7c+s5qYbsyZFxLp0vwG4FjiINLMagGdWs7JpOdFLmiHpVkmrJN0r6YxUPlnSzZIeSPc7dS5cs2JJmihph5HHwDuAlXhmNSuxdo7oh4GFEfFG4BDgdEkzce8Dq7apwG2SfgncCdwQET/EM6tZibXcRp9OPI2cfHpW0ipgGlnvgzlptcvIZp46s60ozUoiIh4C9q1T7pnVrLQ6cjJW0gCwP3AHHex90EntnrVupNWz6CNa6WHRC2WNq1ufo1mVtZ3oJU0CrgY+GhHPSMq1XZ7eB500ODjY1lnrRlo9iz4iT4+FIpQ1rnZ7H5htidrqdSNpa7IkvzQirknF7n1gZlYi7fS6EXAxsCoiPl+zyL0PzMxKpJ3f5ocCJwErJN2Tys4m622wTNL7gV8D72orQjMza0s7vW5uAxo1yPdd7wNP8m1mVeUrY83MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKq58Uwh1QL2RKBfOGm57Nigzs35U+kTv4YPNzNrTtaYbSXMlrZa0RtKibr2OWVm4zltZdSXRS5oAfAU4EpgJnCBpZjdey6wMXOetzLp1RH8QsCYiHoqIF4ArgXldei2zMnCdt9JSRHT+j0rHA3Mj4gPp+UnAwRGxoGad04DT0tO9gdUdD+TlpgC/7fJrtMJxNWesuHaPiF16GcyIPHU+leep92Xd94043u5rFHOuOt+tk7H15pJ92X+UiFgMLO7S629G0l0RMbtXr5eX42pOWeMiR52HfPW+xO+xLsfbfe3G3K2mmyFgRs3z6cC6Lr2WWRm4zltpdSvR/xzYS9IekrYB5gPXdem1zMrAdd5KqytNNxExLGkB8CNgAnBJRNzbjddqQs+aiZrkuJpTyrg6XOdL+R7H4Hi7r62Yu3Iy1szMysNj3ZiZVZwTvZlZxVUq0UuaIelWSask3SvpjDrrzJH0tKR70u0TPYhrO0l3SvpliuuTddaRpC+ly+d/JemAksTV8/1V89oTJP1C0vV1lvV8f3WbpEskbZC0suhY8sjzfSuTPPW9jMb6HuRV+kHNmjQMLIyIuyXtACyXdHNE3DdqvX+PiGN6GNfzwNsiYqOkrYHbJP0gIm6vWedIYK90Oxj4WrovOi7o/f4acQawCnhVnWVF7K9uWwJ8GfhWwXHklff7VhZ563vZjPU9yKVSR/QRsT4i7k6PnyXbOdOKjQoiszE93TrdRp8Fnwd8K617O7CjpN1KEFchJE0Hjga+2WCVnu+vbouInwJPFB1HXmX9vjVS5vreSI7vQS6VSvS1JA0A+wN31Fn8lvTz7QeS9ulRPBMk3QNsAG6OiNFxTQMeqXk+RA++NDniggL2F/BF4GPAiw2WF7K/rL5xvm+lkbO+l8kXGft7kEslE72kScDVwEcj4plRi+8mGx9iX+Ai4Hu9iCkiNkXEfmRXTB4k6U2jVsl1CX0BcfV8f0k6BtgQEcvHWq1OWamPzqpqnO9bqeSo76WR83uQS+USfWp7uxpYGhHXjF4eEc+M/HyLiBuBrSVN6VV8EfEUMAjMHbWo0EvoG8VV0P46FHinpLVko0C+TdLlo9bxkAMlMN73razG+B6WSZ7vQS6VSvSSBFwMrIqIzzdYZ9e0HpIOItsHv+tyXLtI2jE93h74W+D+UatdB5ycepMcAjwdEeuLjquI/RURZ0XE9IgYIBtK4McRceKo1Xq+v+zl8nzfyiTn97A0cn4Pcqlar5tDgZOAFakdDuBs4HUAEfF14Hjgw5KGgT8A86P7lwfvBlymbHKKVwDLIuJ6SR+qietG4ChgDfB74L1djilvXEXsr7pKsL+6StIVwBxgiqQh4JyIuLjYqMZU9/uWfvmVUd36XnBMPeEhEMzMKq5STTdmZrY5J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKc6I3M6u4/w9OeDjGqgF6ZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.hist()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### dummy variables\n",
    "pandas gives you a great deal of control over how categorical variables can be represented. We're going dummify the \"prestige\" column using get_dummies.\n",
    "\n",
    "get_dummies creates a new DataFrame with binary indicator variables for each category/option in the column specified. In this case, prestige has four levels: 1, 2, 3 and 4 (1 being most prestigious). When we call get_dummies, we get a dataframe with four columns, each of which describes one of those levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "d_ranks = pd.get_dummies(df[\"prestige\"],prefix=\"prestige\")\n",
    "#dummy ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prestige_1</th>\n",
       "      <th>prestige_2</th>\n",
       "      <th>prestige_3</th>\n",
       "      <th>prestige_4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   prestige_1  prestige_2  prestige_3  prestige_4\n",
       "0           0           0           1           0\n",
       "1           0           0           1           0\n",
       "2           1           0           0           0\n",
       "3           0           0           0           1\n",
       "4           0           0           0           1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_ranks.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>admit</th>\n",
       "      <th>gre</th>\n",
       "      <th>gpa</th>\n",
       "      <th>prestige_2</th>\n",
       "      <th>prestige_3</th>\n",
       "      <th>prestige_4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>380</td>\n",
       "      <td>3.61</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>660</td>\n",
       "      <td>3.67</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>800</td>\n",
       "      <td>4.00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>640</td>\n",
       "      <td>3.19</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>520</td>\n",
       "      <td>2.93</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   admit  gre   gpa  prestige_2  prestige_3  prestige_4\n",
       "0      0  380  3.61           0           1           0\n",
       "1      1  660  3.67           0           1           0\n",
       "2      1  800  4.00           0           0           0\n",
       "3      1  640  3.19           0           0           1\n",
       "4      0  520  2.93           0           0           1"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# CREATING A CLEAN DATA FRAME\n",
    "cols_to_keep = [\"admit\",\"gre\",\"gpa\"]\n",
    "data = df[cols_to_keep].join(d_ranks.loc[:,\"prestige_2\":])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once that's done, we merge the new dummy columns with the original dataset and get rid of the prestige column which we no longer need.\n",
    "\n",
    "Lastly we're going to add a constant term for our logistic regression. The statsmodels function we would use requires intercepts/constants to be specified explicitly.\n",
    "\n",
    "### Performing the regression\n",
    "Actually doing the logistic regression is quite simple. Specify the column containing the variable you're trying to predict followed by the columns that the model should use to make the prediction.\n",
    "\n",
    "In our case we'll be predicting the admit column using gre, gpa, and the prestige dummy variables prestige_2, prestige_3 and prestige_4. We're going to treat prestige_1 as our baseline and exclude it from our fit. This is done to prevent multicollinearity, or the dummy variable trap caused by including a dummy variable for every single category."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>admit</th>\n",
       "      <th>gre</th>\n",
       "      <th>gpa</th>\n",
       "      <th>prestige_2</th>\n",
       "      <th>prestige_3</th>\n",
       "      <th>prestige_4</th>\n",
       "      <th>intercept</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>380</td>\n",
       "      <td>3.61</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>660</td>\n",
       "      <td>3.67</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>800</td>\n",
       "      <td>4.00</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>640</td>\n",
       "      <td>3.19</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>520</td>\n",
       "      <td>2.93</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   admit  gre   gpa  prestige_2  prestige_3  prestige_4  intercept\n",
       "0      0  380  3.61           0           1           0        1.0\n",
       "1      1  660  3.67           0           1           0        1.0\n",
       "2      1  800  4.00           0           0           0        1.0\n",
       "3      1  640  3.19           0           0           1        1.0\n",
       "4      0  520  2.93           0           0           1        1.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#ADDING THE INTERCEPT MANUALLY\n",
    "data[\"intercept\"] = 1.0\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_cols = data.columns[1:]\n",
    "\n",
    "logit = sm.Logit(data[\"admit\"],data[train_cols])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.573147\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "results  = logit.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we're doing a logistic regression, we're going to use the statsmodels Logit function. For details on other models available in statsmodels, check out their docs here.\n",
    "\n",
    "### Interpreting the results\n",
    "One of my favorite parts about statsmodels is the summary output it gives. If you're coming from R, I think you'll like the output and find it very familiar too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "ir = results.predict([800,4,0,0,0,1.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.73840825]\n"
     ]
    }
   ],
   "source": [
    "print(ir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>admit</td>      <th>  No. Observations:  </th>  <td>   400</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>Logit</td>      <th>  Df Residuals:      </th>  <td>   394</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                 <td>MLE</td>       <th>  Df Model:          </th>  <td>     5</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Sat, 02 Jan 2021</td> <th>  Pseudo R-squ.:     </th>  <td>0.08292</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>10:18:59</td>     <th>  Log-Likelihood:    </th> <td> -229.26</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>             <td>True</td>       <th>  LL-Null:           </th> <td> -249.99</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>  LLR p-value:       </th> <td>7.578e-08</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>gre</th>        <td>    0.0023</td> <td>    0.001</td> <td>    2.070</td> <td> 0.038</td> <td>    0.000</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>gpa</th>        <td>    0.8040</td> <td>    0.332</td> <td>    2.423</td> <td> 0.015</td> <td>    0.154</td> <td>    1.454</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prestige_2</th> <td>   -0.6754</td> <td>    0.316</td> <td>   -2.134</td> <td> 0.033</td> <td>   -1.296</td> <td>   -0.055</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prestige_3</th> <td>   -1.3402</td> <td>    0.345</td> <td>   -3.881</td> <td> 0.000</td> <td>   -2.017</td> <td>   -0.663</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>prestige_4</th> <td>   -1.5515</td> <td>    0.418</td> <td>   -3.713</td> <td> 0.000</td> <td>   -2.370</td> <td>   -0.733</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intercept</th>  <td>   -3.9900</td> <td>    1.140</td> <td>   -3.500</td> <td> 0.000</td> <td>   -6.224</td> <td>   -1.756</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                  admit   No. Observations:                  400\n",
       "Model:                          Logit   Df Residuals:                      394\n",
       "Method:                           MLE   Df Model:                            5\n",
       "Date:                Sat, 02 Jan 2021   Pseudo R-squ.:                 0.08292\n",
       "Time:                        10:18:59   Log-Likelihood:                -229.26\n",
       "converged:                       True   LL-Null:                       -249.99\n",
       "Covariance Type:            nonrobust   LLR p-value:                 7.578e-08\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "gre            0.0023      0.001      2.070      0.038       0.000       0.004\n",
       "gpa            0.8040      0.332      2.423      0.015       0.154       1.454\n",
       "prestige_2    -0.6754      0.316     -2.134      0.033      -1.296      -0.055\n",
       "prestige_3    -1.3402      0.345     -3.881      0.000      -2.017      -0.663\n",
       "prestige_4    -1.5515      0.418     -3.713      0.000      -2.370      -0.733\n",
       "intercept     -3.9900      1.140     -3.500      0.000      -6.224      -1.756\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}