05_linear/05_linear.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression\n",
    "\n",
    "## Princeton Salary Datasets\n",
    "- We've learned the theory. Let's get to a real example\n",
    "- In this task, we will study the impact of various factors (education, experience and gender) on salaries in academia.\n",
    "- The dataset is a historic dataset of salaries in an academic department and is easy to study.\n",
    "- According to the website\n",
    "```\n",
    "These are the salary data used in Weisberg's book, consisting of observations on six variables for 52 tenure-track professors in a small college.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sx</th>\n",
       "      <th>rk</th>\n",
       "      <th>yr</th>\n",
       "      <th>dg</th>\n",
       "      <th>yd</th>\n",
       "      <th>sl</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>male</td>\n",
       "      <td>full</td>\n",
       "      <td>25</td>\n",
       "      <td>doctorate</td>\n",
       "      <td>35</td>\n",
       "      <td>36350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>male</td>\n",
       "      <td>full</td>\n",
       "      <td>13</td>\n",
       "      <td>doctorate</td>\n",
       "      <td>22</td>\n",
       "      <td>35350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>male</td>\n",
       "      <td>full</td>\n",
       "      <td>10</td>\n",
       "      <td>doctorate</td>\n",
       "      <td>23</td>\n",
       "      <td>28200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>female</td>\n",
       "      <td>full</td>\n",
       "      <td>7</td>\n",
       "      <td>doctorate</td>\n",
       "      <td>27</td>\n",
       "      <td>26775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>male</td>\n",
       "      <td>full</td>\n",
       "      <td>19</td>\n",
       "      <td>masters</td>\n",
       "      <td>30</td>\n",
       "      <td>33696</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       sx    rk  yr         dg  yd     sl\n",
       "0    male  full  25  doctorate  35  36350\n",
       "1    male  full  13  doctorate  22  35350\n",
       "2    male  full  10  doctorate  23  28200\n",
       "3  female  full   7  doctorate  27  26775\n",
       "4    male  full  19    masters  30  33696"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sb\n",
    "\n",
    "data = pd.read_csv(\"http://data.princeton.edu/wws509/datasets/salary.dat\", sep='\\s+')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, we start with exploring the dataset for a little."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 52 entries, 0 to 51\n",
      "Data columns (total 6 columns):\n",
      "sx    52 non-null object\n",
      "rk    52 non-null object\n",
      "yr    52 non-null int64\n",
      "dg    52 non-null object\n",
      "yd    52 non-null int64\n",
      "sl    52 non-null int64\n",
      "dtypes: int64(3), object(3)\n",
      "memory usage: 2.8+ KB\n",
      "None \n",
      "\n",
      "       sx    rk  yr         dg  yd     sl\n",
      "0    male  full  25  doctorate  35  36350\n",
      "1    male  full  13  doctorate  22  35350\n",
      "2    male  full  10  doctorate  23  28200\n",
      "3  female  full   7  doctorate  27  26775\n",
      "4    male  full  19    masters  30  33696 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "print data.info(), '\\n'\n",
    "print data.head(), '\\n'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First note we only have a very limited set of data points (52 rows).\n",
    "\n",
    "As we can read on the [Princeton website](http://data.princeton.edu/wws509/datasets/#salary), these columns mean the following:\n",
    "- sx = Sex: male or female,\n",
    "- rk = Rank: full, associate or assistant professor,\n",
    "- yr = Number of years in current rank,\n",
    "- dg = Highest degree: doctorate or masters,\n",
    "- yd = Number of years since highest degree was earned,\n",
    "- sl = Academic year salary, in dollars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Salary Distribution\n",
    "\n",
    "Let's view the distribution of salaries first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAECCAYAAAD0JMwBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADu9JREFUeJzt3H2MZXddx/H3bHenMMx00q2zKNtNCX/s1/4h8iCgFSkN\nCFQhEGJiLK2AD1BAszEq0kJqgmBJmqIiD5JiaakgpKQISKT9w6Zg/yjy4HP9Doi0u6TpDjvLstNt\n2dmd8Y97xw4wuz1z7j0zs9/7fiWTzj17vvf3Pb/76+eeOfdhbHl5GUlSPds2uwFJUjcMeEkqyoCX\npKIMeEkqyoCXpKIMeEkqanvTHSPiKuBlwDjw/sy8sbOuJEkDa3QGHxHPB34uMy8CLgb2dNmUJGlw\nTc/gXwT8e0T8HXAO8IfdtSRJGoamAT9D76z9pcBTgM8AP9lVU5KkwTV9kfU7wB2ZeSIzZ4FHIuLH\nOuxLkjSgpmfw/wTsA94dEU8CngAcWmvH5eXl5bGxsVbNzM7OcsVVH2Niete6a48dOcgt117G3r17\nW42t9fPxkoaqXXCeRqOAz8zPRcTzIuJL9M7635iZa35L2djYGHNzR1s1Mz+/wMT0LibP3d26vu3Y\nXZiZmdpS/Qxbtcdro1RfF+vhXDxqZmZq6PfZ+G2SmflHQx9dktQZP+gkSUUZ8JJUlAEvSUUZ8JJU\nlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEv\nSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUUZ8JJUlAEvSUVt\nb7pjRHwVONK/+c3M/M1uWpIkDUOjgI+IxwFk5iXdtiNJGpamZ/A/DUxExO39mqsz857u2pIkDarp\nNfiHgOsy88XAlcBHI8Lr95K0hTUN6VngowCZ+XXgEPATXTUlSRpc00s0vwH8FPCmiHgScA7wwKl2\nnpmZatXM4cOTrepW7Nw52Xrsrmy1foap4uO1UUb1uNfiXHSnacD/NXBTRHwRWAZem5lLp9p5bu5o\nq2bm5xda1a2ubzt2F2ZmprZUP8NW7fHaKNXXxXo4F4/q4omuUcBn5iLwqqGPLknqjC+USlJRBrwk\nFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXA\nS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JR\nBrwkFWXAS1JRBrwkFbW96Y4RsQv4CvCCzJztriVJ0jA0OoOPiB3AB4GHum1HkjQsTS/RXAd8AHig\nw14kSUP0mJdoIuI1wFxm3hERVwFjp9t///79fPvbc62auf/++1rVASydPDFQPcCePRcwPj4+0H2c\nSY4fP87+/e3n7Ex9vAY97kHGljbS2PLy8ml3iIi7gOX+z9OABF6emQ+utf8Vv/MnyweOndeqmYX9\nX2Jyz7OZPHf3umsPfuurwBgT07tajX3syEFuufYy9u7d26r+TDQ7O8sVV32s9ZwdOnAv551/4Rn3\neA163KO4VrQhTnvy3MZjnsFn5sUrv0fEncDrTxXuABNPOIfJs9f/PzzA4qFzWtX9/9jTu1qFzYr5\n+QXm5o4O1MNqMzNTQ72/YZufXxhozo4dOeUyaGSzHq9Bj3uQsWHrr4uN5Fw8amZmauj36dskJamo\nxm+TBMjMS7pqRJI0XJ7BS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwk\nFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXA\nS1JRBrwkFWXAS1JRBrwkFWXAS1JRBrwkFWXAS1JR25vsFBFnATcAe4Fl4MrM/M8uG5MkDabpGfxL\ngaXMfC7wNuCd3bUkSRqGRgGfmZ8GXt+/+WTgcFcNSZKGo9ElGoDMPBkRNwOvAH6lu5bOPMePH2f/\n/vt+ZPvhw5PMzy88Zv3i4iIAO3bsaDX+nj0XMD4+3qpWUl2NAx4gM18dEU8E7omICzPz4WE2c9a2\nzX3Nd+fOSWZmptZdNzs7y77rPsPE9K5W4x46cC+PnzqvVf2xIwe55drL2L1777prDx+eXHfNVtL2\n8RrGcbcde8UgtdU4F91p+iLr5cD5mfku4GFgqf8zVCeXljb1bT3z8wvMzR1tVTcxvYvJc3e3GvfY\nkQcHqh+k7zPZZh5327GhF2hta6txLh7VxRNd0zP424APR8RdwA5gX2Z+f+jdSJKGplHAZ+Yx4Fc7\n7kWSNER+0EmSijLgJakoA16SijLgJakoA16SijLgJakoA16SijLgJakoA16SijLgJakoA16SijLg\nJakoA16SijLgJakoA16SijLgJakoA16SijLgJakoA16SijLgJakoA16SijLgJakoA16SijLgJako\nA16SijLgJakoA16SijLgJakoA16SitreZKeI2AHcCFwAnA28IzM/22VjkqTBND2DfxUwl5nPA14C\nvLe7liRJw9DoDB64Ffhk//dtwIlu2pEkDUujgM/MhwAiYope2L+1y6bU3NLJE9x//32tatvWbQWj\neNzHjx9n//7Bet+z5wLGx8eH1JG2uqZn8ETEHuA24H2Z+fEumjlr2+a+5rtz5yQzM1Prrjt8eLKD\nbpp5ZOEQ139inonpB9Zde+jAvZx3/oUddNW9zT7utmtlRZva2dlZ9l33GSamd7Ua89iRg9xy7WXs\n3r23VX1XBplHnV7TF1mfCNwBvDEz7+yqmZNLS5v6tp75+QXm5o62qttME9O7mDx397rrjh15sINu\nNs5mHnfbtQK9QGu7ztoe8+r7aNt3F9rORUVdPNE1PYO/GpgGromIa/rbLs3MR4bekSRpKJpeg98H\n7Ou4F0nSEPlBJ0kqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkq\nyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqyoCX\npKIMeEkqyoCXpKIMeEkqyoCXpKIMeEkqat0BHxHPiYg7u2hGkjQ829ezc0S8GbgcWOimHUnSsKz3\nDP4bwCuBsQ56kSQN0boCPjNvA0501IskaYjWdYmmsqWTJ7j//vta1bat05lpkLWyuLjIgw8+gYWF\n4+uuHXSdDdo3wI4dO4Zaf/jwJPPzp7/i29XYTe3ZcwHj4+Otajfblgr4s7Zt3pt6Hlk4xPWfmGdi\n+oF11x46cC/nnX9hB11pKxp0rTx+6jwmpne1qh1knW1W34PWb+bYx44c5JZrL2P37r2txt5sbQN+\neahd9J1cWtrU921OTO9i8tzd6647duTBDrrRVjbIWtnMdbYZfQ9av5ljA8zPLzA3d7RV7XrMzEwN\n/T7XHfCZ+S3goqF3IkkaKj/oJElFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lF\nGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCS\nVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFGfCSVJQBL0lFbW+yU0RsA94PPBX4PvBbmfk/XTYmSRpM\n0zP4VwDjmXkR8Bbg+u5akiQNQ9OA/3ng8wCZeQ/wM511JEkaiqYBfw7wvVW3T/Yv20iStqhG1+Dp\nhfvUqtvbMnNprR0XHzrE0sLRVs0sHn2Axe07W9U+fHQeGGtVO2j9mTr2mdr3qI59pvZ9Jo997MjB\n1uNuBU0D/m7gZcCtEfGzwL+dascb3/v29o+EJGlomgb8p4BfjIi7+7df21E/kqQhGVteXt7sHiRJ\nHfCFUkkqyoCXpKIMeEkqyoCXpKIe8100EfEc4F2ZeUlEPB34LPD1/j+/PzNvjYjfBl4HnADekZmf\ni4jHA38DzABHgVdn5nf6b7P88/6+d2Tm24d/WMMXETuAG4ELgLOBdwD3AjcBS8B/AG/KzOXq83GK\nuTgA/D0w299tJNZGRJwF3ADsBZaBK+l9X9NNjN66WGsuxhnBdbEiInYBXwFeQG893MQGrovTnsFH\nxJvpPWBn9zc9E3h3Zl7S/7k1In4c+F3gIuDFwLURMQ68AfjXzHwe8BHgbf37+Cvg1zLzucBzIuJp\njWZq870KmOsfz0uA99H7Tp6r+9vGgJePyHysNRfPAK4fwbXxUmCp3/PbgD9ldNfFD8/FOxnddbFy\nIvRB4CF66+DdbPC6eKxLNN8AXsmjHwN7JvDLEXFXRHwoIiaBZwN3Z+ZiZn6vX/NUVn1/Tf+/L4yI\nKXpfWva//e23Ay98jB62iluBa/q/bwMWgWdk5hf62/6B3rE8i/rzsdZcjOTayMxPA6/v33wycBh4\n5iiuizXm4ruM6Lrouw74APBA//aG58VpAz4zb6P3p8CKe4A/yMyLgW8Cf0zvKwyOrNrnKDDND35/\nzVrbVm/f8jLzocxc6E/yrfSeUVfP3+pjLD0fa8zFW4EvMbpr42RE3Az8BfBRfvBz8SOzLmDNuRjJ\ndRERr6H3V+4d/U1jbMK6WO+LrJ/KzK+t/A48nR/9npopes/cq7evtY1+w99dZw+bJiL2AP8IfCQz\n/5betbQVK8cyEvPxQ3PxcUZ8bWTmq4EAPgQ8btU/jdS6gB+YixvoXScexXXxWnqf/r8TeBpwM73r\n6Ss2ZF2sN+Bvj4hn9X9/IfBles/QvxARZ0fENHAhvRcQ7gZ+qb/vpcAXMvMocDwinhIRY8CLgC9w\nBoiIJwJ3AG/OzJv6m78WERf3f7+U3rGUn49TzMVIro2IuDwi3tK/+TBwEvjyiK6LH56LJeC2UVwX\nmXlxZj4/My8B/gX4deDzG70umn4Xzcr3GbwB+MuIWKR3Xel1/T/V3wN8kd4TxtWZ+f2I+ABwc0R8\nkd67Ci7r38eV9P50Owu4PTP/uWEPm+1qen8OXRMRK9ef9wHv6b8o8l/AJ/uvilefj7Xm4veAPxvB\ntXEb8OGIuAvYQW9N/Ddwwwiui7Xm4gCjmxmrLQO/zwavC7+LRpKK8oNOklSUAS9JRRnwklSUAS9J\nRRnwklSUAS9JRRnwklSUAS9JRf0fNBZ/qy6fizsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10acb5510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "f = data.sl.hist(bins=20)\n",
    "# f = plt.hist(data.sl, bins=20) #Same thing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What correlates with salary?\n",
    "\n",
    "Next, we would want to what factors may be correlated with salary.  We would expect that as **experience** or **education** increases, the salary should increases as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAADhCAYAAACN6NC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVPWZ7/FPg0D3hAbzmrRb4piQjI+4geCCS0QjEZQY\nGJWIxkSzuMWIMeGaBElwwagTF8QxOpEx4OgdblDEhUvDjCGAZFRQxBiGJ5oQEydmgrlBmoTCVvr+\ncU611U11V3VTZ6v6vl8vXnafOnXqqWPVc57+nd9S19bWhoiIiIiIxKtP0gGIiIiIiNQiFeIiIiIi\nIglQIS4iIiIikgAV4iIiIiIiCVAhLiIiIiKSABXiIiIiIiIJUCEuiTCzfzKz/zSzPgXb+prZz8zs\n+iRjC2O5zsw+F/Nrnm1my+N8TRGRYtKSo5PIxZViZl82s8siPP5vzGxEVMeXeOyRdABSs74BrAW+\nDdwYbvs28DYwI6mg8tw98RhERBKUihyd8Vx8AvDzCI+vhWCqQJ0W9JGkmNmhwGrgRIK7M08CR7r7\nG2b2JeCycPufgK+6u5vZgcDdwPuA/YAXgXPcfYeZ7QAWAcOAzwKfBiYSXDj+BFzo7n/oFMMJwG1A\nX4KkdpO7LzSzucDP3f02M8sBNwGfDF/zTne/M3z+t4HPA+8Ar4SvsbWr+Iucg+uB88J9XgE+6O4n\nm1l/4Jbw3PQF1gFT3L3FzI4GfgD0A34FHAB8HagD7gS2AX8DHAOMBa4B+gN/Baa6+zPha18DnBnG\n+BvgK+7+Rsn/cSJSE1KSo+dSRi7u9JxjgNkEefBtgry33Mx2Ah9w9/8X7rcT+ABwOO/lzvcBVwPf\np4xcambXAh8G9iHIxZuBc4BRwBxgO3Cju99TEN+HgVXAhvC5o4EvAhOA+jCGqe6+qKvjh/8PNgFn\nAb8E/i+w2t2/vcv/SEk1dU2RxLj7ywSJ7V/Cf18Mk8toguL24+4+giAhLgyf9mXgR+5+HPAx4CPA\n6eFj/YDH3f0g4I/AlQQXjaOAZcDRRcK4Frjd3Y8kSIQnh9vbeK+1oT+w2d1PAM4Gbjaz/mb2aeAC\nYJS7HwZsAr5aIv52ZjaBoBAeBhwHDC54zW8Bre4+0t2HA2+Er9sXeAS4xt2HEVxshhc87xBgsrsf\nQZC0bwROC+O4BFhoZn9jZp8HDgWODvddQnDREBEBUpOjS+biwp3NrB9BsX9tmJcvAmaZWV2Jt5vP\nncMJiveycmn43BOAs919KPBn4BJ3fxR4nOD6cg+7+iBwvbsbMAD4BHBimNenA4Xdf3Y5fsFjexKc\nuydUhGeTuqZIotz9n8zsbOAX7r403DyeIIH/zMzyu77fzPYEvgmcamb/CzCCVpGBBYdcFf73dWA9\nsM7MlgBL3P0nRUL4MXC3mZ0B/AfBRaeYx8L/riNImgOBMcCP3f2t8L18A8DM/rGr+N19S8ExxwCP\nuPtfwufdD0wJH/sUMNjMPhn+3h/4H+AwoC1/rtz9p2b2csExf+fuvwt//iSwL/CTgjjeDWP7FHAU\nsDZ8rC/Q0MV7F5EalYIc3VnnXPw+gsI57zDgHXdfEsb/AkFjBwWxFlOYOzv/3l0ubQOWu/u2grje\nX3Ccrv4AeAf4zzDG18zsQuBzZvZRgtb09xXs29Xx64B/BVqBu7p7c5JeahGXNNhE0MUirw/wr+5+\nRNgaMQI4Kixi5xO0cPwGuB14gY6JbhuAu7e5+2iCFus/AXeY2azOL+zuPyRI3P9OcOvxJTMbVCTG\n7fnjhr/XESS/dmY2OLzl2F38hXbS8Tv4bqdzMKXgGEcDk8J9Oif2wudtK/i5D/BU/hjhcY4FXg4f\nu7lg+5EErS4iIp0llqOLKJaLC71Dp77TZnaome1RuH/nlnQ65s7Ov3eXSwFyBfu2dYqpq/6/O9x9\nZxjLCIKifCCwlKBbYuG1oavjtwE3AM8B/9jF60jKqRCXtChMXMuAc81sn/D3y4Cnwp9PJbidtyD8\n/RiC1twOzGxY2FK80d1vBmYR9APsvN9q4Ah3n0dwu29PgtaGUrcx2wha0M80s8Zw23XAVSXiL9QM\nTAoL+D5A4cwAS4ErzKxf+Ngc4HsEfQp3mNnYMP6jCVvJixz/JwQtUxbuOx54iaAP4lLgooLYbyBo\nWRERKSaRHN0LDrSZ2ZjwdUaEsdUR9K8+KtzvzB4cs7tc2vlaUVew7R2Cu5mlfBxY4+6zCO4Y/ANF\nzlmR40NQhH+F4FryyeJPkTRTIS5p0V5IuvsyghaBfzez9cBkgsQEMA141MzWAPcAPyW4Pdj5GOsJ\nup2sDfe9kKBI7uxq4Hoze4Eg2V7r7q/RsbDtXOS2ha+xBPgRsNrMXgL2AqaViL9d+Pz7CWYmeAbY\nUvBaNxC0KK0DfhFu/4a7v0swOOfaMOavA38gGDxU2JcSd98AXAzMN7MXCfrDn+HufyUo7J8Engkv\nhocStEyJiBSTVI4uGkMXv+PuOwiK7Blmti6M4R/cvZWg69/dZvY8wdia35dz7BK5tI1drxf535cA\nU8zsmyXey78BHzCzXxBcD1oIuvoMLHH8fHxvEhTj95vZ4CKvJSmmWVNEMibsg36ru//RzPYnmJXg\nI+6+NeHQREREpAciH6xpZnsBzwOnEPSJnRv+92XgcndvM7OLCP7afAeY6e6LzawBeBBoIvjr8AJ3\nf9PMRhHcwnoHWObuiS/+IhKz14CnzKyV4Bbll1SEi4iIZE+kLeLhNEI/BoYSzI/5fYKWvJVmdg9B\nP9VnCPqbjSSYteFpgoFjXwUGuvv1ZnYOcKy7fy28LfQP7r7JzBYTTOP2YmRvQkREREQkAlH3Ef8+\nQf+s/CIhI9x9ZfjzEoLp244imIS+NWzVe5VgwMbxBIPZCP87JhxY1t/dN4Xbl4bHEBERERHJlMi6\npoRzYm5292Xh6oOdR/q2ECxgMgh4q4vtW7vZlt8+pFQsbW1tbXV1pSbBEBFJpZpLXsrZIpJhPUpe\nUfYR/wLvTSE0HJhH0N87bxDBLBFbgcaC7Y1FthfbVniMbtXV1bF5c0vv3kWFNDU1Jh6D4lAcWYgj\nDTGkLY5ak4acXa60fE7KlaV4sxQrZCveLMUK2Yq3pzk7sq4p7j7a3U9y95MJZnX4PNAcLo0LcBqw\nkmAOzI+b2YBw2p2hBAM5V/PesrinASvdvQV428yGhMvVnhoeQ0REREQkU+Jc4r4N+AZwX7ii1Qbg\n4XDWlNkEk9j3IZiHeUc4mHOema0CdgDnhce5FHiIYLL7pe6+Jsb3ICIiIiJSEbEU4mGreN5JRR6f\nQ7DASOG27cBniuz7LMHSsiIiIiIimaWVNUVEREREEqBCXEREREQkASrERUREREQSEOdgTRFJUC6X\nY9HipQBMHD+W+vr6hCMSkWqjPCPSM2oRF6kBuVyOqdfNpnljA80bG5h63WxyuVzSYYlIFVGeEek5\nFeIiNWDR4qVsGzicPn33oE/fPdg2cFh7q5WISCUoz4j0nApxEREREZEEqBAXqQETx49l4LYX2flu\nKzvfbWXgtvVMHD826bBEpIooz4j0nAZritSA+vp6bp0xpWAQ1RQNohKRilKeEek5FeIiNaK+vp7J\nZ01IOgwRqWLKMyI9o64pIiIiIiIJUCEuIiIiIpIAFeIiIiIiIglQH3ERkQL5lQEbBw7glNEnabCZ\niEhCaiEfq0VcRCRUuDLggrV9tDKgiEhCaiUfqxAXEQlpZUARkXSolXysQlxEREREJAEqxEVEQloZ\nUEQkHWolH2uwpkjG5AevQJCoqnHwSlIKVwYMBgdpZUCRaqVcmm61ko9ViItkSH7wyraBwwF4eu1s\nbp1RnckpKfmVAZuaGtm8uSXpcEQkAsql2VAL+TjSQtzM+gL3AQcCbcClQH/gSeCX4W4/cPcFZnYR\ncDHwDjDT3RebWQPwINAEtAAXuPubZjYKmBXuu8zdr4/yfYikReHgFaB98IqWlBYRKZ9yqaRF1H3E\nPwXsdPcTgOnAjcAI4DZ3Pzn8t8DM9gGuAI4DxgI3mVl/4DJgvbufCDwQHgPgXuDc8LjHmNnwiN+H\niIiIiEhFRVqIu/tjwCXhrx8GtgAjgfFmtsLM5pjZQOBoYLW7t7r7VuBV4HDgeKA5fH4zMMbMGoH+\n7r4p3L4UGBPl+xBJi1oZvCIiEiXlUkmLura2tshfxMzmAROAScAHCVq515nZNOD9wIvAYe7+rYL9\nHwC+BVzh7hvNrA/wGjAKeMTdR4X7fgEY4u7f6SaE6N+kSExyuRzzH34CgMlnn6E+jdWvLukAEqCc\nLZFTLpWI9ChnxzJY090vMLO9gWeB49z99+FDjwJ3ASuBxoKnNBK0nm8t2F5sG8CgcHu3ku7kn5aB\nBoqjOuIYP3YcAC0trbS0tCYWRxTSEEPa4qhFaTj35UjL56RcWYo3jlgrmUt1bqOTpXh7mrMj7Zpi\nZueb2bfCX7cDO4GFZnZUuG0MsBZ4Dvi4mQ0ws8HAUOBlYDVwerjvacBKd28B3jazIWZWB5xKUMiL\niIiIiGRG1C3iC4EfmdkKoB9wJfA6cJeZtQJvABe7+zYzmw2sIvjjYJq77zCze4B5ZrYK2AGcFx73\nUuAhoC+w1N3XRPw+REREREQqKtJC3N3/CpxT5KETiuw7B5jTadt24DNF9n0WOLZCYYqIiIiIxE4L\n+ohIKuRXuQtWUDtJA6dERHpIq4VmT9TziIuIlJRf5a55YwML1vZh6nWzyeVySYclIpIZhXm0eWOD\n8mhGqBAXkcQVrnLXp+8e7avciYhIeZRHs0mFuIiIiIhIAlSIi0jiSq1yl8vlmP/IY8x/5DHdahWR\niqi2vKLVQrNJgzVFJHH19fXcOmNKwWDNKe2DjPL9HrcNHA7A02tnc+uMKRqEJCK9Vo15pTCPAkwc\nn+33UyvUIi4iqVBfX8/ksyZw4fmTOlw81O9RRCqtWvNKPo9OPmuCivCMUCEuIiIiIpIAFeIikmrq\n9ygilaa8ImmhPuIikmrq9ygilaa8ImmhQlykArSa2e7rbmXNfL9HEZGs0fVBuqOuKSK7SauZ7T6t\nrCkicYorb+v6IKWoEBfZTdU6+j5OOociEqe4co5ym5SiQlxEREREJAEqxEV20+6Mvs+v7Db3wQU1\nfbtSMxiISJwmjh9Lw5a1/OHVZ/jDq8/QsOX5SHKOcpuUosGaIrupt6PvO6/stuSn2V/Zrbe6W1lT\nRCQSffqw10dGBD9vfSGSl9DsLFKKCnGRCujNrB6FfQeB9r6DtTo7SP4cNjU1snlzS9LhiEgVW7R4\nKdsHjWjPv9sHHRFZ/tWsT9IddU0REREREUmACnGRhKjvoIhIMpR/JS0i7ZpiZn2B+4ADgTbgUmAH\nMBfYCbwMXO7ubWZ2EXAx8A4w090Xm1kD8CDQBLQAF7j7m2Y2CpgV7rvM3a+P8n2IREH9okVEkqG+\n25IWUbeIfwrY6e4nANOB7wG3AdPc/USgDphgZvsAVwDHAWOBm8ysP3AZsD7c94HwGAD3AueGxz3G\nzIZH/D5EIpHvO3jh+ZN0ERARiVE+/04+a4LyryQm0kLc3R8DLgl//TDwZ2Cku68Mty0BxgBHAavd\nvdXdtwKvAocDxwPN4b7NwBgzawT6u/umcPvS8BgiIiIiIpkReR9xd3/XzOYBdwIPEbSC57UAg4FB\nwFtdbN/azbbC7SIiIiIimRHL9IXufoGZ7Q08BxTe/xkEbCEorBsLtjcW2V5sW+ExutXU1Fhql8il\nIQZQHJ0pjo7SEEcaYoD0xFGLsnTusxQrZCveLMUK2Yo3S7FC9uItV9SDNc8HPuTuNwPbgXeBtWY2\n2t1XAKcBTxEU6Dea2QCCQn0owUDO1cDpwJpw35Xu3mJmb5vZEGATcCpwbalYkp6XOC1zIyuOQC6X\nKxgkeVLi/QOTPh9piiMNMaQtjlqUhnNfjrR8Tsq1O/Hm8yYEs45EnTdr6dzGLUuxQrbi7WnOjrpr\nykLgCDNbQdDH+0rgq8B1ZvYzgj8EHnb3/wFmA6sICvNp7r4DuAc4xMxWAV8GrguPeylBN5dngRfc\nfU3E70OqSH5Fy+aNDSxY24ep182u6eXlRURKKcybzRsblDdFKiTSFnF3/ytwTpGHTiqy7xxgTqdt\n24HPFNn3WeDYykQptUYrWoqI9Izypkg0tKCPiIiIiEgCVIhLzcnCimq5XI75jzzG/Ece0+1fEUlc\nFvJmWimfS3dimTVFJE3SvqJlvi/mtoHBOlVPr53NrTPSFaOI1BatRNk7yudSigpxqUn5FdV6MhI7\nrhkDdqcvZtyzGohI7cjnzd1VKk+lbVar3aG+9VKKuqaIlCELMwZkIUYRqW2l8pRmtZJao0JcMi2u\nvneFrRp9+u7R3qoRhd72xYwzRhGpPZXIt6XyVLXlMfWtl1LUNUUyq1r73qkvpoikTbXm26gpn0sp\nahGXzMpCK3Vv5ftiTj5rQtlJWy0vIhKVSuXbUnmqGvNYb/K51A61iIuUIQutGlmIUURqW6k8lfZZ\nrUQqTYW4ZNbE8WN5eu1stg0cBhC2nEyJ7PUqNWNAlLIQo4hkTyXzbak81ZtZrUSySoW4ZFattgBr\nikIRiVsW861ypWSBCnHJtFprAU7TgCld5EQkrdKUK0W6o8GaIhmSlqm9NGe5SG3J2nc+LblSpBQV\n4iLSY7rIidQWfedFoqFCXCRDqnFqLxGRSlOulKxQIS7dimvlSilPfsDUuINyjDsol1ifR13kROKR\nlhycte98WnKlSCkarCld0mCXdErDANUszqAgkjVpysFZ/M6nIVeKlKJCXLpU2CcQaO8TqMSWTZWe\n5UQXOZFopS0HV+o7rxmXRN6jQlykk2q8SKSpZU1E4pfPa8FqlScl9t1XLhLpSH3EpUtZ6xNYCb2d\noist/Ti7ohkPRLKnUjk4l8tx1YxZNG9sYMHaPlw1Y1ZieUq5SKSjyFrEzawfcD9wADAAmAm8DjwJ\n/DLc7QfuvsDMLgIuBt4BZrr7YjNrAB4EmoAW4AJ3f9PMRgGzwn2Xufv1Ub2HWpfFPoG7qze3gtXC\nIyJRqFQOXrDoSbYPGtGe17YPOoIFi57kc5PPrmi8ItJzUbaIfxbY7O4nAuOAu4ERwG3ufnL4b4GZ\n7QNcARwHjAVuMrP+wGXA+vD5DwDTw+PeC5zr7icAx5jZ8AjfQ83L9wmcfNYEFZZdyEILTy3e3RCp\nBpXIwetf3lDWtjgoF4l01GWLuJn9XXdPdPffljj2AuDh8Oc+QCswMji0TQBeAb4GHA2sdvdWoNXM\nXgUOB44Hbgmf3wx8x8wagf7uvincvhQYA7xYIhbJuN702+7NcyaOH8vTa2ezbeAwgPAiMWU3Ik9H\nn/NqvruRhvMrHVXg+pF6aelzXY5hhw5l3hPL+dDBJwHw+oblnHLGyERiqa+vZ+Y3L2bmrbMBmP7N\n3uWiNH3vs/RZkPTprmvKSqCti8fagCHdHdjd/wIQFs8LgGuAeuA+d19nZtOAGQRF9FsFT20BBgOD\ngK3dbMtv7zYOyb7edP3I94ncPmgEACuem8Ud132tZILsTcHaXfGepm4r1TjLSZrOr3SwW9ePtOv8\nuVvy03R/7iZNPIOfPf8Kr216HoC/328QkyaekUgsuVyO6bf8kG0DRwEw/ZYf9vjc9Ta/RyFrnwVJ\nn7q2tq5yZSDsk308QdeSJ4AjgMvcfUGpg5vZ/sBC4G53n2tmg939rfCxocBdwJ3AOHe/PNy+ELgR\nmAbc7O5rzGww8DRB95Vn3P2QcN8rgT3c/bYSoXT/JiXV5j64gAVr+7T3b9z5biuTjmzjwvMndfmc\n+370v1m0fkCH50wc9jYXfeG8SGLM5XLMf/gJACaffUZ7Eu5N7N0dTzrq7fnNmLqkA+it3bh+pDpn\nZ/FzVyqnxJVzKnHu4s7v3cniZ0Ei16OcXc5gzTuBq4Ezgb8S9PNeSNDK3SUz2xtYBnzF3ZeHm5ea\n2RXuvoagS8la4DngRjMbQNBiPhR4GVgNnA6sAU4DVrp7i5m9bWZDgE3AqcC15bzRzZtbytktMk1N\njYnHkNU4WrbtABp22dbd859+Zh00jNpl28RPdWwFquT5GD92XBBbSystLa3tcZYTe2Ecu7aw3Bhb\nC0saPh9RfzaiiCNKTU2NSYewO3p1/YDkc3Z3ovzcRWn82HHtn+t8joJ4c05Pzl1X38Fy83scsvpZ\nSEt+K1eW4u1pzi5nsGYfd18BjAceCfv29S3jedMIupN818yWm9ly4CrgjvDnYwlmSPkfYDawCngK\nmObuO4B7gEPMbBXwZeC68LiXAg8BzwIvhEW9VLHeDO4ZduhQXt+wvP05r29YzrBDh8YU8Xt6E3sW\nBn+mhQZ+pV5vrx+pVm2fuzhzTiXOXVryO1TfZ0HiV07XlBUEtxSnAgcDnwPOCmczyYq2pP+SSstf\nc1mNo6cDc3K5HFd953Ze+3Pw+wHvhztu+Pouz4vjfJQTe2Ec8x95jOaNDR1udY47KBdL/+40fD6i\n/mxEFUdUmpoas9w1pbfXj8RzdilZHaBX7HMdd84p9zvb1Xew3Pwelyx+FtKS38qVpXh7mrPL6Zry\nWeCLwJnu/v/MbD8g/o5YUtN6OtCwvr6em675ynsj86fu/m3W3hZ8PY09iplbopCWWQuqcRBqFana\n60f+cxd1gRDH9yzunJPL5Vi7bj0A404Z3eP3VF9fzx03fD0V+ScfTxyfBalOJVvEq0TirStp+YLW\nShyd+zwO3PZi0T6P5cZR7vF6q3McSRW5aTgftfIZ7UEcmW0R3w2J5+xyRfk5ieJ71l0rcxw5Z8uW\nLZxzyXT2PfzTALzx0uP8n3+eyZ577ll2rGmVpXizFCtkK96e5mwtcS9VqdJ9HuPut532hZTUj10k\nenF+z+LKOTNvnc2+h3+6/T3te/gZ7XcuRWpRZEvci4hUSlq6wYhIOignSLVQi7hUpUqPZNfI+I7i\nPB/52/PNGxto3tjA1Otmk8vlInktkTSpxrwzfeoU3njp8fb39MZLTzB9as/6oysnSDVRH/GYpKV/\nUy3F0V2LSW9GuUfZApPF/y9xzVaS1CwyKfp/oj7iKZa1wZpp+Fxv2bKlw0D6zv3DS+XnJGeW6k4a\nzm25shQrZCveKGZNEcmkrmbT6O2SxJqdoyOdD5HoVeP3bM899+TWmd8t+piWjJdao64pUnM00DBZ\nuVyO+Y88xvxHHivrdnI13p4XSZuefi+jUk5+Vk6QaqIWcRGJTefWrqfXlm7tqq+v59YZUwpuz6t1\nTKSSevO9TJJyglQTtYhLzVFrSnJ6ezci7dM5imRZmu4SlpuflROkWqhFXGpOYWtKMBhIrSkiImmg\n/Cy1Ri3iUpPyrSkXnj9JST5Guhshkj5p+14qP0stUYu41KTeTF8ou099O0XSJ23fS+VnqSUqxKXm\naHqsZFXjdGwiWZeW76Xys9QadU2RmpOmgUkiIvIe5WepNWoRl25FuZqkiEhWqfuEiFSCWsSlS/lb\nhM0bG2je2MDU62aXtdBDWhaG6EraBiaJSLYU5sYFa/t0mRvTngvTSPlZao0KcelSb24R9rZ4j1N+\nYNK4g3JMOrJN/Q+rmAohiUI5uTHOXFhNn3PlZ6k1KsSlorLSv0/TY2VLbwqNLPxRKNUrrlyYts95\nJf4oUH6WWqJCPKXS0MKhW4SSBr0tNLLyR6Fkz8TxY2nYspY/vPoMf3j1GRq2PJ9YbkzT5zxtfxSI\nZEFkgzXNrB9wP3AAMACYCfwXMBfYCbwMXO7ubWZ2EXAx8A4w090Xm1kD8CDQBLQAF7j7m2Y2CpgV\n7rvM3a+P6j0kpfP0TU+vTWb6pt7MLTtx/FieXjubbQOHAYTF+5TIY5XqVVhoAO2FRhqmWpMa1qcP\ne31kRPDz1hd2ebgWc6G+qyI9F2WL+GeBze5+IjAOuBu4DZgWbqsDJpjZPsAVwHHAWOAmM+sPXAas\nD/d9AJgeHvde4Fx3PwE4xsyGR/geEpGmFo78LcLJZ00o6w+Bwv594w7KqX+fJEZ3dCQqixYvZfug\nEe05evugI3bJ0XHlQn3ORbItyukLFwAPhz/3AVqBEe6+Mty2BDgVeBdY7e6tQKuZvQocDhwP3BLu\n2wx8x8wagf7uvincvhQYA7wY4fuQHkrLwhBSHXrbspi21QKl9sSRC9P0Oa/FuwAiuyuyQtzd/wIQ\nFs8LCFq0by3YpQUYDAwC3upi+9ZutuW3D4kg/EQpmYm8Z3cKDf1RKFFIW45Oy+c8TX8UiGRFXVtb\nW2QHN7P9gYXA3e4+18x+5+77h49NIGjNXgaMc/fLw+0LgRuBacDN7r7GzAYDTxN0X3nG3Q8J970S\n2MPdbysRSnRvMiK5XI75Dz8BwOSzz1AyE6lddUkHkIDU52zlaBHpQo9ydpSDNfcmKLK/4u7Lw83r\nzGy0u68ATgOeAp4DbjSzAUA9MJRgIOdq4HRgTbjvSndvMbO3zWwIsImga8u15cSzeXNLxd5bbzQ1\nNfY4hvFjxwHQ0tJKS0trYnFEQXEojjTHkLY4alEazn0p48eOa/+cVCpHRy0tn+tyZClWyFa8WYoV\nshVvT3N2lH3EpxF0J/mumX033HYlMDscjLkBeDicNWU2sIqgL/k0d99hZvcA88xsFbADOC88xqXA\nQ0BfYKm7r4nwPYhkSn7ZbQhun6uVTkQkWcrL0p0o+4hfSVB4d3ZSkX3nAHM6bdsOfKbIvs8Cx1Ym\nSpHqkZZpL0VEJKC8LKVoQR+RKpGmaS9FRER5WUpTIS4iIiIikgAV4lK1crkc8x95jPmPPFYTyyxr\nYQ8RSVqt5d1SlJellCgHa4okphb75WkOXxFJUi3m3VKUl6UUFeJVJOsjsysZf2G/PKC9X14aFr2I\nUloW9hCRytiyZQszb50NwPSpU9hzzz0TjqhrtZp3S1Felu6oa0qVyLdENG9soHljA1Ovm52p24JZ\nj19EpNK2bNnCOZdM582GUbzZMIpzLpnOli1bkg5LRCpIhXiVyPrI7ErHr355IpJ1M2+dzb6Hf7o9\nL+57+BntreOF0tIvW3lXpOfUNaVKtLa2Ag1FttUm9csTkVqQpn7ZyrsiPadCPIOK96Vu4/UNy/nQ\nwScB8Po3nZijAAAOgklEQVSG5fDRkckF2UMTx4/l6bWz2TZwGEDYkjJlt46pfnkikmXTp05h0pe/\nTb899wegdcvvmDXnpg77pK1ftvKuSM+oEM+Yrlo/+vXrz74HHs/m19YDsO+Bx9Ov384kQ+0RtaSI\niHRUX1/PkCH7s2PwkQAMeAvlRZEqoz7iGdNVX+pxp4zmjxuaaTpgGE0HDOOPG4JtWZJvSZl81gRd\nbESk5i1avJQdg49sz/c7Bo/cZezMuFNG88ZLj7f3y37jpScyl/tFapkK8QRVcoBN81Mr2OvgcWx+\nbT2bX1vPXgePpfmpFbu81twHF2g2EhGRGMQxiLL5qRX87YGf4FfPP86vnn+cvz3w5A65X0TSTYV4\nQno7XV93o9L77tGfvYccyd5DjqTvHv2LvtaCtX00NaCISMQqMSVrObOQtLa+zR9/vYaPjpzAR0dO\n4I+/XkNr69uVfCsiEiEV4gnp7XR99fX1zPzmxXxg+7N8YPuzzPzmxdTX13ebsLM+taGISNaUm3e7\nazXPj50Zd1COcQflupgNpY4PHXxy++t86OCTgbro3piIVJQGa2ZMLpdj+i0/ZNvAUQBMv+WH7cm5\nt4Mds74ip4hIFuVyOa76zu289uegcF7xzC+444avd8jBpWYh6devX1nbRCSd1CKekHGnjOb3BQNs\nfl/mAJvuWlm6GuzYXWu5VrQUEam8cgZRLlj0BP7fb7HXR0ay10dG4v/9FgsWPdGj19EiOiLZpkI8\nIU80/wd7Fwyu3PvgsTzR/B+RvFbh7c1JR7Z1uL2pbisiIpVXagA9wAsv/py/O+SU9vz7d4d8ghde\n/HmPXqe87isiklbqmpKQ9S9voG/DKPYeEswPu/PdVta/vI7PlXhebxe+ybeWNzU1snlzy+6GLyIi\nJeQH0EOQ4yGau41aREcku9QinpBhhw7l9Q3L228nvr5hOcMOHVryeZVu/dBtTRGRyisnt44Yftgu\n14ERww9LKGIRSYJaxBMyaeIZ/Oz5V3ht0/MA/P1+g5g08Yz2x7sbQFnJ1g+taCkiUnnl5NZJE89g\n1XP/xYY1jwJw8JD9OlwHRKT6RV6Im9kxwM3ufrKZHQE8AbwSPvwDd19gZhcBFwPvADPdfbGZNQAP\nAk1AC3CBu79pZqOAWeG+y9z9+qjfQxTq6+u544avFy22u1rGPqoCWbc1RUQqr5zc2rdffz521JnB\nz1tfiCMsEUmRSLummNnVwH3AgHDTSOB2dz85/LfAzPYBrgCOA8YCN5lZf+AyYL27nwg8AEwPj3Ev\ncK67nwAcY2bDo3wPSdAAShGR7Cu1suaixUvZPmhEe67fPugI5XqRGhN1H/FXgTN5b3WBkcB4M1th\nZnPMbCBwNLDa3VvdfWv4nMOB44Hm8HnNwBgzawT6u/umcPtSYEzE7yESmjZQRKR6KceLSDkiLcTd\nfSFBF5K8Z4Gp7j4a+DUwA2gE3irYpwUYDAwCtnazrXB75nTX6q0BlCIi2VbOnU3lehGJe7Dmo+6e\nL7ofBe4CVhIU43mNwBaCgruxm20QFOZbynnhpqbG0jtFrDCGxoEDdnm8ceCAcJ9G5v3TNcx/OFjY\nYfLZ11S0f3gazgUojs4UR7pigPTEUYuydO6Lxdp9jm/fEmmu70rWz22aZSneLMUK2Yu3XHEX4kvN\n7Ap3X0PQpWQt8Bxwo5kNAOqBocDLwGrgdGANcBqw0t1bzOxtMxsCbAJOBa4t54WTnju78/zdp4w+\niSU/7Tgf+Cmjp3TYZ/zYcQC0tLTS0tIaSRxJURyKI80xpC2OWpSGc1+Orj4n5eT4vChyfU/jTaMs\nxQrZijdLsUK24u1pzo6rEG8L/3sZcJeZtQJvABe7+zYzmw2sIugqM83dd5jZPcA8M1sF7ADOC49x\nKfAQ0BdYGhb1maNpA0VEqpdyvIiUI/JC3N1/QzAjCu6+DjihyD5zgDmdtm0HPlNk32eBY6OINW6a\nNlBEpHopx4tIKVpZU0REREQkASrERUREREQSoEJcRERERCQBKsRFRERERBKgQlxEREREJAEqxEVE\nREREEqBCXEREREQkASrERUREREQSoEJcRERERCQBKsRFRERERBKgQlxEREREJAEqxEVEREREEqBC\nXEREREQkASrERUREREQSoEJcRERERCQBKsRFRERERBKgQlxEREREJAEqxEVEREREEqBCXEREREQk\nASrERUREREQSsEfUL2BmxwA3u/vJZvYxYC6wE3gZuNzd28zsIuBi4B1gprsvNrMG4EGgCWgBLnD3\nN81sFDAr3HeZu18f9XsQEREREam0SFvEzexq4D5gQLjpdmCau58I1AETzGwf4ArgOGAscJOZ9Qcu\nA9aH+z4ATA+PcS9wrrufABxjZsOjfA8iIiIiIlGIumvKq8CZBEU3wAh3Xxn+vAQYAxwFrHb3Vnff\nGj7ncOB4oDnctxkYY2aNQH933xRuXxoeQ0REREQkUyItxN19IUEXkry6gp9bgMHAIOCtLrZv7WZb\n4XYRERERkUyJvI94JzsLfh4EbCEorBsLtjcW2V5sW+ExSqlramosvVfE0hADKI7OFEdHaYgjDTFA\neuKoQanI2eXKUqyQrXizFCtkK94sxQrZi7dccc+ass7MRoc/nwasBJ4DPm5mA8xsMDCUYCDnauD0\nwn3dvQV428yGmFkdcGp4DBERERGRTImrRbwt/O83gPvCwZgbgIfDWVNmA6sI/jCY5u47zOweYJ6Z\nrQJ2AOeFx7gUeAjoCyx19zUxvQcRERERkYqpa2trK72XiIiIiIhUlBb0ERERERFJgApxEREREZEE\nqBAXEREREUmACnERERERkQTEPY947MJpDl8Hfhlu+k93nxbj6/cBfkCwWugO4Mvu/qu4Xr9TLC/w\n3uJJv3b3L8X42scAN7v7yWb2MWAuwbzyLwOXu3sso4Y7xXEE8ATwSvjwPe7+4xhi6AfcDxwADABm\nAv9FjOekixheB57kve9K5OfDzPoC9wEHEsyudCnB92QuMX4+uoijPzGfj4J49gKeB04hOA9zSeD7\nkpSk83Y50pTby5Fk/u+JtFwrypGG60k50nDNKVdark3lqsQ1rOoLceCjwPPu/umEXn8i0N/djwu/\ntLeF22JlZvUA7n5yAq99NXA+sC3cdDvBNJUrw2kqJwCLEohjJHC7u98e9Wt38llgs7t/zszeD6wH\n1hHvOSkWw3XAbTGfj08BO939hHCNge+F2+P+fHSO40aCi2rc5yN/Ifpn4C8EqxEn8n1JWNJ5uxyp\nyO3lSDL/90RarhXlSNH1pBxpuOaUKy3XpnLt9jWsFrqmjAQ+aGY/MbPFZnZgzK9/PNAM4O7PAkfG\n/Pp5w4C/MbOlZvZUeOGIy6vAmQRFBcAId88vxLQEGJNQHCOB8Wa2wszmmNnAmOJYAHw3/LkP0Er8\n56RYDLGfD3d/DLgk/PXDwJ+BkXF/PorEsYXkPh/fB+4B3gh/T+r7kqSk83Y50pLby5Fk/u+JtFwr\nypGW60k50nDNKVcqrk3lqsQ1rKoKcTP7kpn9vPAf8Hvge+7+CYK/VB6MOaxBwNaC398Nb2nG7S/A\n9919LOGiSHHF4e4LgXcKNtUV/LwNGJxQHM8CU919NPBrYEZMcfzF3beZWSNB0plOx+9i5OekSAzX\nEKxym8T5eNfM5gF3EizWldTno3McsZ8PM7uQoDVoWbipjoTOR1xSmrfLkZbcXo7E8n9PpOVaUY60\nXE/KkYZrTrnSdG0q1+5ew1L3Rdwd7v4v7n5Y4T9gLfB4+PhqYL+Yw9oKNBb83sfdd8YcAwR9qx4C\ncPdXgD8B+yYQBwT9pvIaCVofk/Cou68Lf14EHBHXC5vZ/sBPgAfc/d9I4Jx0imE+CZ4Pd78AMGAO\nUF/wUKyfj4I47gOWJXA+vgB80syWA8OBeUBTweNJfl8ikdK8XY605PZypCn/90RarhXlSCx/liMN\n15xypenaVK7duYZVVSHehRnA1wDMbBjw25hffzVwevj6o4CXYn79vC8S9GHEzPYjaM15o9tnRGdd\n2JcK4DRgZXc7R2ipmR0V/nwKwcU/cma2N7AMuNrd54abYz0nXcQQ+/kws/PN7Fvhr9uBd4G1cX8+\nisSxE1gY9/lw99HuflLYl/dF4PNAc0q+L3FKOm+XIy25vRxpyv89kZZrRTkSuZ6UIw3XnHKl5dpU\nrkpcw2phsObNwINmNp6gr9GFMb/+owQtXKvD378Q8+vn/Qsw18xWEYzs/UICrTf5UcPfAO4zs/7A\nBuDhhOK4DLjLzFoJLkoXx/T60whuVX3XzPJ94a4EZsd4TorFcBVwR8znYyHwIzNbAfQjOA8bif/z\nUSyO10nm81GojeS/L0lIOm+XIy25vRxpyP89kZZrRTmSvp6UIw3XnHKl5dpUrt2+htW1tSU+W42I\niIiISM2pha4pIiIiIiKpo0JcRERERCQBKsRFRERERBKgQlxEREREJAEqxEVEREREEqBCXEREREQk\nASrERSrEzD5sZpuSjkNERESyQYW4iIiIiEgCamFlTZGKMLMHgFXufl/4+3Lgm8C9BKurrU8wPBER\n6aSLvH0MsAQ4BPiMu7+UYIhS49QiLlK++4HzAczsAKAJmANMdfeRwK8TjE1ERHZVLG8/C6x394NU\nhEvSVIiLlG8FsF+YzD8PPADs5+4/CR+fm1RgIiJSVLG8DUExLpI4FeIiZXL3NmAecB4wCXgQqCvY\n5d0k4hIRkeKK5O1/DR/anlhQIgVUiIv0zFzgUuC37v574DUzOz187LzEohIRka7M5b28/UbCsYh0\noEJcpAfc/XXgt7zXDeV8YIaZvQAMIRi0KSIiKVEkb4ukRl1bm+oGkXKZ2X7AT4FD3L014XBERKQE\n5W1JM7WIi5TJzM4GXgS+pWQuIpJ+ytuSdmoRFxERERFJgFrERUREREQSoEJcRERERCQBKsRFRERE\nRBKgQlxEREREJAEqxEVEREREEvD/ARBeIX2CCl4EAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d7ad410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(ncols=2,sharey=True)  #Force both subplots to use the same y axis\n",
    "f = data.plot(kind='scatter', x='yd', y='sl', ax=axes[0], title=\"Years since degree\", figsize=(12,3))\n",
    "f = data.plot(kind='scatter', x='yr', y='sl', ax=axes[1], title=\"Years in current rank\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Most importantly, while using **linear** regression, we should see if the relationship is **linear**.  \n",
    "- Does every increase in X lead to a constant increase in Y?  We may expect for the first 5-10 year of experience this true, but perhaps after 10 year the salaries stop increasing or jump dramatically--Linear regression would be unable to capture these nonlinearities (Unless we transform the data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### sklearn\n",
    "\n",
    "We will first use the **linear_model** module within sklearn to solve the linear regression coefficients for us.\n",
    "\n",
    "Let's first start with a simple model only trying to predict a professor's salary based on `yr` and `yd` (as well as the intercept)\n",
    "\n",
    "Our model will conceptually be this:\n",
    "\n",
    "$y = intercept + \\beta_1*yr + \\beta_2*yd$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import package and create the LinearRegression object.\n",
    "# At creation of object, note the arguments for \"fit_intercept\" (i.e. include an intercept term or not) \n",
    "from sklearn.linear_model import LinearRegression\n",
    "model = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   yr  yd\n",
      "0  25  35\n",
      "1  13  22\n",
      "2  10  23\n",
      "3   7  27\n",
      "4  19  30 \n",
      "\n",
      "Salary \n",
      "0    36350\n",
      "1    35350\n",
      "2    28200\n",
      "3    26775\n",
      "4    33696\n",
      "Name: sl, dtype: int64 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Create our X matrix of observations, and y vector of responses\n",
    "X = data[ ['yr','yd'] ]\n",
    "y = data[ 'sl' ]\n",
    "\n",
    "\n",
    "print X.head(), '\\n'\n",
    "print 'Salary \\n', y[0:5], '\\n'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5782019226431887"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit the model, and score\n",
    "model.fit(X, y)    # Fitting takes the X, y, and estimates for \"beta\". Beta = (X'X)^(-1) * X'y\n",
    "model.score(X, y)  # for regression, this is R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16555.689673113389"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# intercept of model\n",
    "model.intercept_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 489.29445905,  222.25130361])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# beta/slope per dimension\n",
    "model.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(489.29445904610566, 'yr'), (222.25130361475976, 'yd')]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# these go with each dimension\n",
    "coefs = zip(model.coef_, X.columns)\n",
    "coefs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sl = 16555.7 + 489.3 yr + 222.3 yd\n"
     ]
    }
   ],
   "source": [
    "# so the linear model would be the following\n",
    "print \"sl = %.1f + \" % model.intercept_ + \\\n",
    "    \" + \".join(\"%.1f %s\" % coef for coef in coefs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an intuition in the equation!\n",
    "\n",
    "* The base salary is almost \\$17K\n",
    "* Each year in the current position adds another \\$500\n",
    "* Each year after your degree adds another \\$200.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction\n",
    "\n",
    "Predicting new values is simple once we have fit our model. Pass in a matrix **X** of data points that you would like to predict to the model.predict function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 36566.84677578  27806.04632024  26560.41424671  25981.53608403\n",
      "  32519.82350343  29051.67839376  23667.73138879  28384.92448292\n",
      "  29584.05674916  29806.30805277  27316.75186119  28117.88132749\n",
      "  24737.61196598  26960.12500213  26293.37109128  23314.52044066\n",
      "  27361.54371301  25049.44697323  24782.4038178   24158.7338033\n",
      "  29496.18100099  27359.83575754  22870.01783343  25804.07663224\n",
      "  23626.35544791  23002.68543341  25049.44697323  22558.18282618\n",
      "  19579.33217556  21801.8452117   24782.4038178   28827.71913468\n",
      "  24070.85805513  25847.16052858  25936.74423222  19045.24586469\n",
      "  23581.56359609  19401.87272376  19624.12402737  19401.87272376\n",
      "  18912.57826471  20468.33739001  18111.44879842  18690.3269611\n",
      "  18201.03250205  17756.52989482  18867.78641289  17978.78119844\n",
      "  17267.23543577  17267.23543577  17267.23543577  17000.19228034] \n",
      "\n",
      "   yr  yd\n",
      "0  25  35\n",
      "1  13  22\n",
      "2  10  23\n",
      "3   7  27\n",
      "4  19  30\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X) # This is the same as y_pred = X*beta = X * (X'X)^(-1)*X'y\n",
    "print y_pred, '\\n'\n",
    "\n",
    "print X.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can always manually calculate a prediction based on the input values and beta coefficients.\n",
    "\n",
    "If you have had your degree for 6 years, and you've been in your current position for 2 years, what is your predicted salary?\n",
    "\n",
    "$salary = 16555.7 + 489.3 * yr + 222.3 * yd$\n",
    "\n",
    "$salary = $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36550"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "16555 + 489*25 + 222*35"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train/ Test\n",
    "Supervised problems, like linear regression, always require a training and a testset, to avoid having overfitted your model. (i.e. if we fit on the entire data set, it's like giving all the answers to a final exam)\n",
    "\n",
    "Similar to our KNN example, lets split up our data into train and test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=.8)\n",
    "\n",
    "model.fit(X_train,y_train) #Fit the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.593815495439\n",
      "0.310531863204\n"
     ]
    }
   ],
   "source": [
    "print model.score(X_train,y_train) #Print how high the R^2 is on the training data\n",
    "print model.score(X_test,y_test)   #Print how high the R^2 is on the test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are not so great numbers, and not that consistent either. This is not that surprising if you realize how little samples we have (52) and how little structure we saw in the plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "- We've seen how sklearn handles **fitting**, **predicting**, and **scoring** the models/data. \n",
    "- However, these are actually all things that are simple enough for us to do.\n",
    "\n",
    "We will first have a chance to use the sklearn features to do this for us automatically. Next, we will manually calculate these values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------\n",
    "#### SKLearn Implementation\n",
    "\n",
    "**1) First, set up your X and y matrix and vector. For now, use the entire dataset to create X as the Pandas dataframe containing yr and yd, and y as a numpy array containing salary.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create our X matrix of observations, and y vector of responses\n",
    "X = data[ ['yr','yd'] ]\n",
    "y = data[ 'sl' ]  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2) Next, fit a linear model using X and y. Print out the `intercept`, `yr`, and `yd` coefficient values.**\n",
    "\n",
    "(model.fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = LinearRegression()\n",
    "model.fit(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3) Generate the predicted y values given the fitted model. Save these as y_pred**\n",
    "\n",
    "(model.predict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_pred = model.predict(X) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4) Generate the $R^2$ value of the model for X and y.**\n",
    "\n",
    "(model.score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5782019226431887"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.score(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $R^2$ should roughly be around 57% as seen above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------\n",
    "#### StudentLearn Implementation\n",
    "\n",
    "**5) Set up your X and y matrix and vector.**\n",
    "\n",
    "For now, use the entire dataset to create X as a Pandas dataframe, and y as a numpy array. HOWEVER, when creating X, set the first column as all values of \"1\", second column as the yr values, and third column as the yd values. Adding an artificial column of 1's to our matrix is essentially adding an intercept term into model.\n",
    "\n",
    "**IF you initially use a dictionary to create the dataframe like I did, be aware that Pandas will resort your columns alphabetically!!!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>int</th>\n",
       "      <th>yr</th>\n",
       "      <th>yd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>19</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   int  yr  yd\n",
       "0    1  25  35\n",
       "1    1  13  22\n",
       "2    1  10  23\n",
       "3    1   7  27\n",
       "4    1  19  30"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = data[ ['yr','yd'] ].copy() \n",
    "y = data[ 'sl' ].copy()\n",
    "\n",
    "X['int'] = [1 for x in range(52)] \n",
    "X = X[ ['int','yr','yd' ] ]\n",
    "X.head()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6) Now your X should have 52 rows and 3 columns, y should be a vector of 52 elements. Now that we have X and y, let's solve for the beta coefficients.**\n",
    "\n",
    "**Recall that $ \\widehat{beta} = {(X^{T}X)}^{-1}X^{T}y $**\n",
    "\n",
    "**Calculate the coefficients now and save as `betas`. Confirm that the coefficients match what you generated in #2 especially the order.**\n",
    "\n",
    "Some helpful Pandas operations: data.transpose(), data.dot()\n",
    "\n",
    "Some helpful Numpy operations: np.dot(), np.transpose(), np.linalg.inv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 16555.68967311,    489.29445905,    222.25130361])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "betas = np.dot( np.dot( np.linalg.inv(X.transpose().dot(X)),X.transpose()) , y)\n",
    "betas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**7) Predict values of y for the X matrix you have created. $y_{pred} = X * betas$. Save the prediction as y_pred.**\n",
    "\n",
    "Hint: You may need to reshape the betas into an array of (3,1) in order to properly use np.dot or the Pandas dot function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_pred = X.dot(betas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**8) Calculate $R^2$. For OLS with an intercept, $R^2 = 1.0 - SSE/SST = \\frac{\\sum{(y_{pred} - y_{true})^2}}{\\sum{(y_{true} - E[y_{true}])^2}}$, where the $E[*]$ indicates taking the mean of the quantity**\n",
    "\n",
    "Flatten your y_pred first before doing calculations. This will prevent issues with numpy broadcasting array shapes and other funky errors -- y_pred.flatten()\n",
    "\n",
    "**Confirm that your $R^2$ value matches the one from #4**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5782019226431887"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sse = np.sum(np.power(y_pred - y,2))\n",
    "sst = np.sum(np.power(y - np.mean(y),2))\n",
    "\n",
    "\n",
    "1.0 - sse/sst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We'd like to somehow incorporate the categorical features as well, even though they're not numerical values\n",
    "- One way to do this is adding columns for each category, with a boolean flag\n",
    "- An easy off-the-shelves solution is `patsy`, a library for describing statistical models using symbolic formulas and building design matrices. Patsy is very similar to `R`.\n",
    "\n",
    "http://patsy.readthedocs.org/en/latest/overview.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "full         20\n",
       "assistant    18\n",
       "associate    14\n",
       "dtype: int64"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from patsy import dmatrices\n",
    "data.rk.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'd like to try a linear model predicting the salary, with features sex, years in rank, years since graduation, and rank."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y, X = dmatrices('sl ~ sx + yr + yd + rk', data=data, return_type='dataframe')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sl</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>36350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>35350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>28200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>26775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>33696</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      sl\n",
       "0  36350\n",
       "1  35350\n",
       "2  28200\n",
       "3  26775\n",
       "4  33696"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X is now our feature matrix\n",
    "- `Intercept` always contains a one, which offers a convenient way of computing the intercept\n",
    "- `sx[T.male]` is 1 if `sx == 'male'`, and 0 otherwise\n",
    "- `rk[T.associate]` is 1 if `rk == 'associate'`, and 0 otherwise\n",
    "- `rk[T.full]` is 1 if `rk == 'full'`, and 0 otherwise\n",
    "- `yr` is just `yr`\n",
    "- `yd` is just `yd`\n",
    "\n",
    "\n",
    "Note that the two `rk` columns can never be both one. If they're both zero, then `rk == assistant`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note that `y` is now a one-column dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Intercept</th>\n",
       "      <th>sx[T.male]</th>\n",
       "      <th>rk[T.associate]</th>\n",
       "      <th>rk[T.full]</th>\n",
       "      <th>yr</th>\n",
       "      <th>yd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>25</td>\n",
       "      <td>35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>19</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Intercept  sx[T.male]  rk[T.associate]  rk[T.full]  yr  yd\n",
       "0          1           1                0           1  25  35\n",
       "1          1           1                0           1  13  22\n",
       "2          1           1                0           1  10  23\n",
       "3          1           0                0           1   7  27\n",
       "4          1           1                0           1  19  30"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.84905558181749508"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = LinearRegression()\n",
    "model.fit(X, y).score(X, y) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Better score for sure!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 16754.56213785]\n",
      "[[     0.           -845.51782754   5017.30953105  10162.07521799\n",
      "     440.18151428    -55.30310331]]\n"
     ]
    }
   ],
   "source": [
    "print model.intercept_\n",
    "print model.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the first coefficient is zero, since we already fitted an intercept.\n",
    "\n",
    "(**Advanced: Look up 'rank deficient least squares'. This is what solvers do when there are duplicate columns**)\n",
    "\n",
    "Alternatively, we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.849055581817\n",
      "0.0\n",
      "[[ 16754.56213785   -845.51782754   5017.30953105  10162.07521799\n",
      "     440.18151428    -55.30310331]]\n"
     ]
    }
   ],
   "source": [
    "model = LinearRegression(fit_intercept=False)\n",
    "print model.fit(X, y).score(X, y)\n",
    "print model.intercept_\n",
    "print model.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note that `model.coef_` is now a matrix, since we gave the model a dataframe `y` instead of a series. It is possible to fit more `y`s at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sl = 16754.6 Intercept + -845.5 sx[T.male] + 5017.3 rk[T.associate] + 10162.1 rk[T.full] + 440.2 yr + -55.3 yd\n"
     ]
    }
   ],
   "source": [
    "coefs = zip(model.coef_[0], X.columns)\n",
    "print \"sl = \" + \" + \".join(\"%.1f %s\" % coef for coef in coefs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "We've so far done the following\n",
    "* Explored using sklearn's functions to build and predict models\n",
    "* Developed our own implementation\n",
    "* Added categorical features to our model features\n",
    "\n",
    "Let's see how the score changes as we try different transformations of our data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1) Using patsy, create a model using the entire data set of X and y on the following formula:**\n",
    "**`sl ~ sx + yr + yd`**\n",
    "\n",
    "Create the model, fit, and report the score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2) Let's plot out a scatter of yd on the x-axis and sl on the y-axis. We had done this earlier, but let's revisit the scatter**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3) Hm...it's somewhat linear up until 20 years since graduating, and seems more to flatten out afterwards. Let's try a transformation that will help capture this relationship. Scatter np.power(yd,1.0/3.0) on the x axis, and sl on the y-axis. Compare with the original scatter**\n",
    "\n",
    "By taking the cube root of yd, the relationship between yd and sl might be more linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4) The relationship between the cube root of yd and sl looks more linear than the original. Add this transformation to your dataframe. Call it `cubedyd`. Let's try building a new model of `sl ~ sx + yr + cubedyd`.**\n",
    "\n",
    "Build the model and report the resulting score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Your score should have improved from 58% to 60%.\n",
    "\n",
    "We want to try one more thing. Try fitting and building a new model:\n",
    "\n",
    "`sl ~ sx + yr + cubedyd + np.random.randn(52) + np.random.random(52)`.\n",
    "\n",
    "**5) What is your score now?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, your score should have increased from 60% to a higher $R^2$. Even though we added noise terms to our model, we fit non-zero coefficients to them. This is what is meant by 'gaming' $R^2$--just by adding terms, we were able to increase the score of our model.\n",
    "\n",
    "As a thought exercise, when there are 52 rows, all we need are 52 columns to get a square matrix and a perfect solution. i.e. $R^2$ = 1. \n",
    "\n",
    "** So we have a problem. We're sure that using the transformed `yd` factor should have increased our $R^2$. But, we also just saw that adding pure noise also increases our $R^2$. How do we get a measure of $R^2$ that isn't gamed just by adding additional factors?**\n",
    "\n",
    "**Solution** Adjusted $R^2$. The Adjusted $R^2$ is similar, but also penalizes for adding extraneous factors.\n",
    "\n",
    "(More here https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2) \n",
    "\n",
    "To get Adjusted R2, lets use `statsmodel`, which provides much more information for regression analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### statsmodels\n",
    "\n",
    "- Statsmodels is a relatively new package that provides convenient utilities for investigating the results of a model. \n",
    "- It uses `patsy` to provide R formula syntax\n",
    "\n",
    "A formula allows you to write a functional relationship between variables.  \n",
    "Example:\n",
    "```R\n",
    "Y ~ X1 + X2 + X3\n",
    "```\n",
    "\n",
    "It automatically assumes there  is an intercept term. You can make this explicit by using \n",
    "\n",
    "```R\n",
    "Y ~ 1 + X1 + X2 + X3\n",
    "```\n",
    "\n",
    "http://statsmodels.sourceforge.net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'data' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-e730cbbcf30b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mols\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformula\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"sl ~ yr\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msummary\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'data' is not defined"
     ]
    }
   ],
   "source": [
    "model = sm.ols(formula=\"sl ~ yr\", data=data).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use categorical features as well, and `statsmodels` will use `patsy` for that just as we have been before. i.e. dont need to manually create your own design matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = sm.ols(formula=\"sl ~ sx + yr\", data=data).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = sm.ols(formula=\"sl ~ sx + yr + rk\", data=data).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see `+` is not acting as an addition operator but as a separator between other variables.\n",
    "\n",
    "There are other operators that lose their algebraic meaning in a formula. \n",
    "\n",
    "- `+` adds a new variable.\n",
    "- `:` adds the _interaction_ of two variables. \n",
    "- `*` adds the original terms as well as their interaction effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Adding an interaction effect\n",
    "model = sm.ols(\" sl ~ sx*rk\", data).fit() #This is like sl ~ sx + rk + sx*rk\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- You can also convert a numerical field to a categorical variable, if you think that might be a better representation.\n",
    "- Please be aware though that changing numerical values into categories will increase the number of features in your dataset, as each category will have its own column.\n",
    "- This will increase the complexity of the model, and hence your chances of **overfitting**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Converting a field to a categorical variable\n",
    "model = sm.ols(\" sl ~ sx*rk + C(yr)\", data).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seaborn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The visualization library `seaborn` includes some linear regression functionality as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In scatter plots, seaborn can add the trendline and confidence bounds automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = sb.lmplot(x = 'yd', y='sl', data=data)\n",
    "f = plt.title(\"Salary vs years since degree\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = sb.lmplot(x = 'yd', y='sl', col='sx', data=data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = sb.lmplot(x = 'yr', y='sl', col='rk', row='sx', data=data)\n",
    "f = plt.ylim(0, 50000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously, we're a little short on data to make these last plots convincing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overfitting and regularization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that we haven't looked at cross-validation yet, while supervised problems, like linear regression, always require a training and a testset.\n",
    "- However, if you have a lot of data, and not so many features, trying to fit all these samples on a single line (or hyperplane in more dimensions) is a strong example of regularization already.\n",
    "- With many features, few samples, or a complex model (e.g., polynomial regression), overfitting is very likely, though, and not only need you to cross-validate your model, also you need to introduce regularization to make your model simpler.\n",
    "\n",
    "We will explore regularization and more into cross-validation after the break."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Further reading\n",
    "\n",
    "- [Linear Regression with Python](http://connor-johnson.com/2014/02/18/linear-regression-with-python/)\n",
    "- [Statsmodels Documentation](http://statsmodels.sourceforge.net/stable/index.html)\n",
    "- [Python 538 Model](https://github.com/jseabold/538model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise, we'll use one of sklearn's standard datasets to analyze Boston house prices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "boston = load_boston()\n",
    "X, y = boston.data, boston.target\n",
    "features = boston.feature_names\n",
    "# print boston.DESCR  # print this to get information on the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Let's make a histogram of all house prices in the dataset. \n",
    "- It's a little unclear in what units the prices are, let's assume it's in $K."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Find a linear model that fits the data, and with which you feel comfortable.\n",
    "- Compute the average $R^2$ of the model using cross validation sklearn.cross_validation cross_val_score with cv=10.\n",
    "- (*) Compute the MSE of the model as well. Remember that $MSE = \\frac{1}{N}\\sum{{(y_{pred}-y)}^2}$\n",
    "- Prices often behave expontential, rather than linear (people often say +2%, rather than +\\$200). Let's try modeling the logarithm of the price and see if that improves your model: `y = np.log(boston.target)` \n",
    "\n",
    "- (*) How much more worth would a house be in a crime-free community, compared to the same house in a community with a crime rate of 6 per resident? (Note that this might not be a fair analysis as features might correlate with themselves  as well.)\n",
    "- (**) What would be a good way of analyzing which feature has the biggest influence on the house price?\n",
    "\n",
    "- (\\*\\*\\*) Let's also examine the effect of outliers. Taking the entire X and entire y dataset, try fitting a model and report its score and coefficients. Then, multiply one of the values by 100000000 (or some commensurately high value). Now refit and report its new score and coefficients. How much did the coefficients change?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "# What's Next?\n",
    "\n",
    "Things we didn't cover, but will next time:\n",
    "* Having linearly dependent columns. What happens when theres multiple very similar columns? How do the coefficients vary? \n",
    "* Using linear regression for weight analysis. \n",
    "* Regularization of beta coefficients for better out of sample prediction"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}