01 - Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Session 01: Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the wine quality data for this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\hat{y} = \\hat{f}(X) + \\epsilon$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\epsilon_i = y_i - \\hat{y}_i$ is the ith residual, the difference between the ith response value and our predicted value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We defined the _residual sum of squares_ RSS as:\n",
    "\n",
    "$RSS = \\sum_{i} e_{i}^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The least squares approach chooses parameters to minimize the RSS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The accuracy of $\\hat{y}$ depends on two quantities, which we call as reducible error and irreducible error.\n",
    "\n",
    "The inability of $\\hat{f}$ to be a good estimate of $f$ is called reducible error, since it is possible to improve this by using the most appropriate statistical learning techniqie.\n",
    "\n",
    "However, even if $\\hat{f}$ is a perfect estimate, Y is also a function of $\\epsilon$, which cannot be predicted using X. This is called irreducible error, because no matter how well we estimate $f$, we cannot reduct this error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is more popularly seen in the following form:\n",
    "\n",
    "$$E(Y - \\hat{Y})^2 = E[f(X) + \\epsilon - \\hat{f}(X)]^2  \\\\\n",
    "  = [f(X) - \\hat{f}(X)]^2 + Var(\\epsilon)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Univariate Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is basically linear regression with 1 predictor variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model:   \n",
    "$y = \\theta_0 + \\theta_1 x$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways to fit a linear regression model, depending on the number of features. For a smaller number of features (say <10k), you can use the Normal equation form to solve for $\\theta$ directly.\n",
    "\n",
    "For more than 10k features, you probably want to use gradient descent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that we write the model in matrix form as:\n",
    "\n",
    "$Y = X \\Theta$,\n",
    "\n",
    "where $\\Theta \\in \\mathbb{R}^{n}, X \\in \\mathbb{R}^{m x n}, y \\in \\mathbb{R}^{m}$\n",
    "\n",
    "The Normal equation is given by:\n",
    "\n",
    "$\\Theta = (X^T X)^{-1} X^T y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be shown that by similfying the above Normal equation, for the single predictor case, we get:\n",
    "    \n",
    "$$\\theta_1 = \\frac{\\sum_{i = 1}^{n} (x_i - \\overline{x}) (y_i - \\overline{y})}{\\sum_{i = 1}^{n} (x_i - \\overline{x})^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\theta_0 = \\overline{y} - \\theta_1 \\overline{x}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "cars = np.genfromtxt('cars.csv', delimiter=',', names=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.ones(len(cars)) # Models the intercept\n",
    "x = cars['speed']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([x0, x]).T\n",
    "y = cars['dist']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-17.57909489,   3.93240876])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VdXV//HPMkKJUIlWjRhowaHUAW0gVC0OwaFotZriUKwDtlpaqz62P0Wh9bH2aS3BAdSqKIpTVYYighXUohgnCpWIGBFRCooEhLYYJAIKYf/+2DcaIMOdz7nnft+vly+Sk3vuWWyuK/uuu/Y+5pxDRESia6egAxARkcxSohcRiTglehGRiFOiFxGJOCV6EZGIU6IXEYk4JXoRkYhTohcRiTglehGRiNs56AAA9thjD9e9e/egw4jLp59+SseOHYMOI1AaA0/j4GkcvCDGobq6+j/OuT3belwoEn337t2ZN29e0GHEpaqqivLy8qDDCJTGwNM4eBoHL4hxMLMP4nmcSjciIhGnRC8iEnFK9CIiEadELyIScUr0IiIR12aiN7P7zWyNmb3V5NhNZvaOmb1pZk+YWVGTnw03syVmttjMBmQqcBGRdJs6v5Z+lbPoMWw6/SpnMXV+bdAhpUU8M/oHgZO2OzYTOMQ5dyjwLjAcwMwOAgYBB8fOucvMCtIWrYhIhkydX8vwKTXU1m3EAbV1Gxk+pSYSyb7NRO+cewlYu92xvzvntsS+nQN0jX19OjDBOfeZc24ZsAT4ThrjFRHJiJueXczGzQ3bHNu4uYGbnl0cUETpk44FUz8FJsa+LsEn/kYrYsd2YGZDgCEAxcXFVFVVpSGUzKuvr8+ZWDNFY+BpHLyojMOgbuuhW3M/WR/X3y/M45BSojez3wJbgEcbDzXzsGbvPu6cGwuMBSgrK3O5srJOqwA1Bo00Dl5UxuG3lbOordu4w/GSokIuP7e8zfPDPA5Jd92Y2WDgVOBc51xjMl/Btr8TuwIrkw9PRCQ7hg7oSWG7bT9SLGxXwNABPQOKKH2SSvRmdhJwDXCac25Dkx89CQwys6+YWQ/gAOCfqYcpIpJZFaUljBjYi5KiQgw/kx8xsBcVpc1Wn3NKm6UbMxsPlAN7mNkK4Hf4LpuvADPNDGCOc+4XzrmFZjYJeBtf0rnUOdfQ/DOLiIRLRWlJJBL79tpM9M65c5o5PK6Vx98A3JBKUCIikj5aGSsiEnFK9CIiEadELyIScUr0IiIRp0QvIhJxSvQiIhGnRC8iEnFK9CIiEadELyIScUr0IiIRp0QvIhJxSvQiIhGnRC8iEnFK9CIiEadELyIScUr0IiIRp0QvIhJxSvQiIhGnRC8iEnFK9CIiEadELyIScUr0IiIR12aiN7P7zWyNmb3V5NjuZjbTzN6L/blb7LiZ2e1mtsTM3jSz3pkMXkRE2hbPjP5B4KTtjg0DnnfOHQA8H/se4GTggNh/Q4Ax6QlTRESS1Waid869BKzd7vDpwEOxrx8CKpocf9h5c4AiM+uSrmBFRCRxydboi51zqwBif+4VO14CfNjkcStix0REJCA7p/n5rJljrtkHmg3Bl3coLi6mqqoqzaFkRn19fc7EmikaA0/j4GkcvDCPQ7KJfrWZdXHOrYqVZtbEjq8AujV5XFdgZXNP4JwbC4wFKCsrc+Xl5UmGkl1VVVXkSqyZojHwNA6exsEL8zgkW7p5Ehgc+3owMK3J8Qti3TdHAOsaSzwiIhKMNmf0ZjYeKAf2MLMVwO+ASmCSmV0ELAfOij18BvB9YAmwAfhJBmIWEZEEtJnonXPntPCj45t5rAMuTTUoERFJH62MFRGJOCV6EZGIU6IXEYm4dPfRi4jkrKnza7np2cWsrNvIPkWFDB3Qk4rS3F/zqUQvIoJP8kMnL2Bzg1/jWVu3kaGTFwDkfLJX6UZEBPj93xZ+keQbbW5w/P5vCwOKKH2U6EVEgI83bE7oeC5RohcRiTglehERoKiwXULHc4kSvYgIcP1pB9Nup2034G23k3H9aQcHFFH6qOtGRIQvO2vUXikiEmEVpSWRSOzbU+lGRCTilOhFRCJOiV5EJOKU6EVEUjB1fi39KmdRU7uOfpWzmDq/NuiQdqAPY0VEkjR1fi3Dp9SwcXMDdPP74wyfUgOEa38czehFIqZxhtlj2PTQzjCj4qZnF/sk38TGzQ3c9Ozi+J7ggw9gw4YMRLYtJXqRCGmcYdbWbcTx5QxTyT4zVtZtTOj4F956Cy64AHr3hgULMhDZtpToRSIk5RmmJGSfosKEjjN7Npx2Gpx4Ihx0ECxdCkcemcEIPSV6kQhJeoYpSRk6oCeF7Qq2OVbYroChA3p+ecA5mD4djjkGzj8fTjkFli2DYcOgc+esxKkPY0UiZJ+iQmqbSeotzjAlJU23TYD1lDTdNmHLFpg4EUaOhIICn9jPOAN2zn7aTWlGb2a/NrOFZvaWmY03sw5m1sPM5prZe2Y20czapytYEWldXDNMSauK0hJeHXYcvUo68+qw46j41u5w111wwAFw331w003w+uvwox8FkuQhhURvZiXA/wBlzrlDgAJgEDASGO2cOwD4GLgoHYGKZELUOlQqSksYMbAXJUWFGFBSVMiIgb1C1eoXVTuvXw833AA9esDMmTB+PLzwAgwYAGZtP0EmY0vD+YVmthnYBVgFHAf8OPbzh4DrgTEpXkck7bbpgSa8PdCJiurGXKG1ciWMHs3hY8fCwIE+uR94YNBRbSPpGb1zrha4GViOT/DrgGqgzjm3JfawFYBecRJK6lCRlLz7LvzsZ3DIIbBlC/PuvRceeCB0SR7AnHNtP6q5E812Ax4HfgTUAX+Nff8759z+scd0A2Y453o1c/4QYAhAcXFxnwkTJiQVR7bV19fTqVOnoMMIVFTGoKZ2XYs/61XSdjdEVMYhVfk2Dl9dvJhu48dTtGABtRUV1FZUsKVz50DGoX///tXOubK2HpdK6eYEYJlz7t8AZjYF+C5QZGY7x2b1XYGVzZ3snBsLjAUoKytz5eXlKYSSPVVVVeRKrJkSlTH4beWsZjtUSooKufzc8jbPj8o4pCovxsE5eP5530GzeDFceSVMn06Pjh3pEXtImMchla6b5cARZraLmRlwPPA28AJwZuwxg4FpqYUokhnqUJE2NTTA5MnQty9ccQWcdx4sWeK/7tgx6OjilvSM3jk318wmA68DW4D5+Bn6dGCCmf0xdmxcOgIVSbeo3jpu6vzayP2dsu6zz+CRR+DGG2H33eG66+DUU2Gn3FxjmlLXjXPud8Dvtju8FPhOKs8rki1R61CJaidR1qxfD/fcA6NHw2GHwb33wtFHB94emarc/PUkIs1SJ1GS1qyBa6/1PfCvvw4zZvj/jjkm55M8KNGLRIr2uknQsmVw2WXwrW/B2rXwz3/CY4/52XyEKNGLREjCuynmq5oa/8Fq376w666waJHftmDffYOOLCOU6EUiRJ1EbXjlFf+h6oABcOih8K9/wZ/+BMXFQUeWUdq9UiRCotpJlJKtW329vbISVq+GoUN9y2SHDkFHljVK9CIRE7VOoqRt3vzlNsHt2sHw4X4vmoKCts+NGCV6EYmWDRtg3Di45RbYbz8YNQpOOCES3TPJUqIXkWhYuxbuvBPuuAOOOgomTYLvaEkPKNGLSK5bscIvcHrwQaiogBdf9O2SSYjqqmJ13YhIbnrnHbjoIt89A7BggS/ZpJDkh0+pobZuI44vVxXn+s1oQIleRHLNa6/5e68eeyx07+43GbvlFujaNaWnjfKqYpVuRCT8JQvn4LnnfIvkkiVw1VXw8MNp3UEyyquKlehF8lyoN0JraIApU3yC/+wzuOYaGDTIt0um2T5Fhc3enyAKq4pVuhHJc6EsWWzaBGPH+nr7rbfC738Pb74J55+fkSQP0V5VrBm9SJ4LVcnik0/g7rt9cu/dGx54gKkde/iy0itPZ7SsFOVVxUr0InkuFCWL1avhttv8LP6kk+CZZ+DQQ7NeVorqqmKVbkTyXKAli6VL4Ze/hAMP9LP5117zd3aKtUyGsqyUg5ToRfJcRWkJIwb2oqSoEMPfHH3EwF5tzmynzq+lX+UsamrX0a9yVmL95gsWwI9/7Feu7r6774m/4w5/448mQlVWymEq3YhIwiWLbUoq3eIsqTgHL7/sO2gWLIBf/9rX43fdtcXrhKKsFAGa0YtIwhIqqWzdCk8+Cf36wcUX+x0kly71vfCtJHmIdidMNmlGLyIJi6uksnmzvy3fjTdCYaHfJriiIqFtgqPcCZNNSvQikrBWSyqffgr33ee3JejZE26/HY47LultgqPaCZNNKt2ISMKaK6l02fwp970/3X+g+sorfkXrzJlw/PF5vRd8GKQ0ozezIuA+4BDAAT8FFgMTge7A+8DZzrmPU4pSJENCv8dLSDUtqXT671JGvvIMP6x5jvZnn+WT/De/GXCE0lSqpZvbgGecc2eaWXtgF+A3wPPOuUozGwYMA65J8ToiaRfqPV5yQEWHT6hY/Bc2T5lCuyFDYOJCKNG4hVHSpRsz2xU4BhgH4Jz73DlXB5wOPBR72ENARapBimSCFuMkac4c+OEPobwc9t+fuY88AjfdpCQfYuacS+5Es28DY4G3gcOAauAKoNY5V9TkcR8753Zr5vwhwBCA4uLiPhMmTEgqjmyrr6+nU6dOQYcRqKiMQU3tuhZ/1qukc5vnR2Uc4uIcu732Gl8fP54OH33Eh2efzUcnn8zWDh3yaxxaEcQ49O/fv9o5V9bW41JJ9GXAHKCfc26umd0GfAJcHk+ib6qsrMzNmzcvqTiyraqqivLy8qDDCFRUxqBf5axmO0dKigp5ddhxbZ4flXFo1ZYtMHkyjBzptwweNgzOPht2/rLqmxfjEIcgxsHM4kr0qXTdrABWOOfmxr6fDPQGVptZl1gQXYA1KVxDJGO0GKcVmzb5Vas9e/obbv/xj19uW7CzurJzTdKJ3jn3EfChmTX+X3E8vozzJDA4dmwwMC2lCEUyJNk9XiJt3Tq/RUGPHjBjhr+L08svwymnqEUyh6X6q/ly4NFYx81S4Cf4Xx6TzOwiYDlwVorXEMmYKC7GSapl9KOP/B7w997rk/rMmXDIIXFdZ1C39fy2cpZaU0MspUTvnHsDaK4+dHwqzysiyUm4ZXTJErj5Zpg0Cc47D6qr/Q23E7lOvJuaSWC0MlYkQuJuGZ0/39979cgjYa+9YPFiv1VBHEk+oetIKCjRiyQopX3YM6zVzcacg6oqfwenH/zA7wW/dCn83//Bnnum7zoh0Phv1GPY9ND9GwVBH5+LJCDsJYvmNhszt5VBK+fDkX+Aujq4+mqYNg2+8pW0XqfxeNC04nlHmtGLJCDsJYumLaPtGjZz1pszeW7cpVw9f4pP8AsXwk9/mlKS3/46jcLSmhr2f6MgaEYvkoCwlywqSkso2PApy0bcypkv/ZUP9+7OqhtuZr8hZ6W1PbKitIR5H6xl/NwPASgw44w+4ehgCvu/URCU6EUSEOaSBf/5D/z5z/zgrrv8/u9Vz7JP794ZudTU+bU8Xl1LQ2xlfYNzPF5dS9k3dg882Yf63yggKt2IJCCUJYvly+GKK/zWwB99BLNnw8SJkKEkD+Euj4Ty3yhgmtGLJKDpPuywnpIM72Hf6uKnhQv9bfqeesrfi3XhQujSJSNxbC/M5RHdfnBHSvQiCWpcTVtVVcXl55Zn7DotdY/stqCaY58YB3Pn+pn8bbdBUVEbz5ZeYS+PRHHFcyqU6EVCapvyiHOUL53HJXMfp2v9f+GG/4UJE/xNtwMwdEDPbX4JgcojYaZELxJSK+s2UrC1gVPeeZlL5kwGYMwRZzHjW0ex5JLTAo0tzF03siMlepEw2riRSxf9nbNfnMiqXfdk5LEXUrVvHzCjJATlkTB33ciOlOhFwqSuDu66C26/nXO/dRj/r+Jq/rH3l+WQsJRHWuu6UaIPHyV6iYyktucNi5Ur/TbB48b5fWhmzaLLQQfxo/m1LM/C3ynRsQtz143sSIleIiFn9zd57z1/Y+3Jk+GCC/yukl//+hc/zkb3SDJjF/auG9mWFkxJJIR5AU+zqqv9vVe/+13YZx94910/o2+S5LMlmbHToqTcohm9REJOlBKcgxde8LfqW7QIrrwS7r8fOnUKNKxkxi7bC8ckNUr0EgmhLiU0NPhtgSsrYf16uOYaf5Pt9u2DjgxIfuyytXBMUqfSjURCKEsJn33mP1w96CAYORJ+8xu/TcGFF4YmyUNIx07SSjN6iYRQ7W+yfj2MHQujR/sbbN9zDxx7bFq3CU6nUI2dZIQSvURG4Pub/Pvf/r6rd98NJ5wAf/sblJZmPYxk2kwDHzvJKJVuRFL1/vtw+eXQs6ffE37OHBg/PrAkP3xKDbV1G3F82SqZ7/dMzXcpJ3ozKzCz+Wb2VOz7HmY218zeM7OJZhaeYqRIGnVcuhTOPx/69PGdM2+/DWPGwH77BRZTzrWZSlako3RzBbAI2DX2/UhgtHNugpndDVwEjEnDdUTC4dVXobKSQ2fP9vdhveMO6Ny51VOSXbWrFauSDinN6M2sK3AKcF/sewOOAybHHvIQUJHKNURCwTmYPh2OPtqvYD31VOaOH+9bJeNI8smUU5I5r6WWyFC0mUpgUi3d3ApcDWyNff81oM45tyX2/QpAn/BI7tqyBR59FA49FK69Fi67DBYvhp//nK1xtkgmW07RilVJF3OxbUYTPtHsVOD7zrlfmlk5cBXwE+Afzrn9Y4/pBsxwzvVq5vwhwBCA4uLiPhMmTEjub5Bl9fX1dAp4JWPQ8mEMdtq0iS5PP023iRPZ2KULy885h4/79t2mRTLecaipXdfiz3qVtPxuINnz6jZuZvW6TXzesJX2BTtR3LkDRYXt2owzWfnweohHEOPQv3//audcWVuPS6VG3w84zcy+D3TA1+hvBYrMbOfYrL4rsLK5k51zY4GxAGVlZa68vDyFULKnqqqKXIk1UyI9Bh9/DHfeCX/+M/TrB9Om0eHww9mtmYfGOw6/rZzV7MrTkqLCVleUJntetkX69ZCAMI9D0qUb59xw51xX51x3YBAwyzl3LvACcGbsYYOBaSlHKZJptbVw1VWw//6wdClUVcGUKXD44Sk/dbLlFJVhJF0ysWDqGmCCmf0RmA+My8A1RNLiuSdeYuOfKjn6rZd5ps8AdnvkGQac3Det10h25Wk2V6zm9F7+0qa0JHrnXBVQFft6KfCddDyvSMa89hq1w6+n9B+v8lDvU/nfIWOpK9yVwtn/ZcTetWlPcsmuPA3rfvSSW7QyVvKHc/Dcc357gjPPZFJhD476+Thu73cOdYV+GUg+Li7SIqvoU6KXuE2dX0u/ylnU1K6jX+Ws0C2rv3ZqDfsNn0H3YdPZb/gMrp3qZ6U0NPg7OPXtC1dc4fvglyzh9oNPZmP7Djs8T74tLtIiq+jTpmYSl23e3ncL39v7a6fW8Mic5V983+Ack175F2V/n0zFc+Pha1+D666DU0+Fnfz8JtR72GeRxiH6NKOXuIT97f34uR9+8XWnzzYwZO7jvHTPRRQ9OwPuuw9mz4bTTvsiyYO6WhppHKJPM3qJS9jf3jc4xx6ffsyF1X/jx288w0s9SvnJWdezaK99ef+YY5o9Jxe6WrLRDZPsODTGNqjben5bOUudOiGmRC9xCfXb+2XL+MPMMfzg7Rd58sBjOf2CUXxYtDcABW3c7CPMXS3Z7IZJdBzCXsqTbal0I3EJ5dv7N9+Ec8+Fvn3p1qMLJ1w8huu+d8kXSR7gnMO7BRdfTDb3usmWMMcmO1Kil7hUlJYwYmAvSmIz+JKiQkYM7JX92Ztz8PLLcMopcNJJ8O1vw9KllE++l5OOP+yLGXyBGecd8XX+WLHDNktZl2zZK8zlsjDHJjtS6Ubi1vj2vqqqKvt7rWzd6rcJrqyENWv8PvCPPw4dvmyP/GNFr1Ak9u0lW/YKc7kszLHJjjSjl3DbvBkefhh69YLrr4df/QreeQd+9rNtknyYRXGvmzDHJjvSjF7C6dNPYdw4uOUWv9HYrbf6Fa1tfLgaRrmw102imsYG6ykJUWyyIyV6ybiEWgTXroU77uCzW2/nH/scxOjyX/GfAw9j6B49qcjBJN8ozHvdJCvQUp4kRIleMiruFsEVK2DUKHjwQT445ntc8qMRvN15H/8zte6JpEQ1esmoNtvw3nkHfvpTf6u+nXaCN9/kx0cM+TLJN3eOiCREiV7ilsymZi212+319hswcCAceyzsuy8sWQI33wxduybdutcYX49h0+OOL5Vzwrq5m8j2VLqRuCS7EnKbNjznOPr9+VwyZzI9PlkNf7gWHnkEdtml5XO2Ox5XfMQXX8rnaEWo5AjN6CUuya6EHDqgJx0L4NRFL/HUQ7/i2ln3MbV0AP98dg5cdtkOSb7xnERb95KJL1vniARNM3qJS1LllE2bqPjnU5zw8Aj+VfBVRh91Lot7H81VJx/I6a3MfitKS5j3wVrGz/2QBucoMOOMPq13nyQTX7bOEQmaEr3EJaFyyrp1cPfdcNtt0Ls3ncY/wmFHHRX3zYOnzq/l8epaGpwD/M6Uj1fXUvaN3VtM9smUe7J1jkjQVLqRuMRVTvnoIxg+3H+4WlMDzzwDTz0FRx2V0LWSKY8kU+7J1jkiQdOMPmSysf94MlpdCfmvf/mOmYkT/W6S1dXQvXvS10qmPJLMKtJUz9GKUMkVSvQhks39x5Oxw0rIN96Ac66CmTPhkkt8T/xee6V8nWTLI8msIk3lHK0IlVyRdOnGzLqZ2QtmtsjMFprZFbHju5vZTDN7L/bnbukLN9pyoqPDOTq/8QacfLLfKrisDJYtgz/8IS1JHlQeEUm3VGr0W4ArnXMHAkcAl5rZQcAw4Hnn3AHA87HvJQ6h7ujYuhWmToXvfpeeo0bBmWfC0qVw5ZXw1a+m9VIVpSWc0adkm73l2+q6ySYtmJJck3Tpxjm3ClgV+3q9mS0CSoDTgfLYwx4CqoBrUooyT4Syo+Pzz+Gxx2DkSOjYEYYP559FRZQff3zGLplM1022aMGU5KK0dN2YWXegFJgLFMd+CTT+MkjP+/k8EKqSRX293xp4//3h0UfhjjvgtdfgjDOgoKDt81MQ5hJWmGMTaYm52Kwp6Scw6wS8CNzgnJtiZnXOuaImP//YObdDnd7MhgBDAIqLi/tMmDAhpTiypb6+nk6dOmXs+es2bmb1uk183rCV9gU7Udy5A0WF7TJ2ve21W7eOkieeYJ9p06g77DA+POcc1vfsuU1su7XfysefZy62mtp1Lf6sV0nntF8vEU1jKy6E1U3egAUdW1Ay/f9ErghiHPr371/tnCtr63EpJXozawc8BTzrnBsVO7YYKHfOrTKzLkCVc67VKWlZWZmbN29e0nFkU1VVFeXl5UGHkX7Ll/ttgh9+2M/ahw6Fb37zix83LVlc2WsLt9TsTGG7gozcN7Zf5axmS1glRYW8Ouy4tF4rUU1jaxwHCEdsQYns/xMJCmIczCyuRJ9K140B44BFjUk+5klgcOzrwcC0ZK8hWfD223DhhVBaCu3bw1tvwb33bpPkIbWSRaI7RIaqhLWdMMcm0pJU+uj7AecDNWb2RuzYb4BKYJKZXQQsB85KLURpS1KLrObM8TfanjMH/ud//DbBu7XcCZvK1sGJrg3QLfRE0iuVrptXgJbu7Za5lgzZRkKJ1Dl49lmf4D/4wJdnxo+Hwra7epLtCGrtnUBbK1DDmjy1YEpyjfa6yXFxlVS2bPEJvbQUrr4ahgyB996DX/4yriQPyZcsQr02QCRPaAuEHNdqIt24ER58EG66Cbp2hT/9ya9oTeIm28mWLEK5NkAkzyjR57jmEumum+r55TszYd+LoW9f+MtfoF+/lK+VTMli6ICe25SWIL53AmHd3E0kFynR57imiXTP+rVcNG8ag978O58c9z2/2dghhwQaXzIfrIZ9czeRXKNEn+MqSkvouHwZn/6pkvI3X2Rm6YnM/etMBpxyeNChfSHRD1aT/QBXRJqnRJ/LXn8dRo7kxFmz/AerTz3AWXvuGXRUKdMHuCLppa6bXOMcvPACDBgAp58ORxzhtwn+/e8hAkkeWv6gVh/giiRHiT5XbN0KU6b4xH7JJTBokL+z069/DRHbZ0SrT0XSK+9KNznXzfH55/DII3DjjbDrrv6erKefDjul9js6zOMQ5pWxIrkorxJ9TnVzrF/v95wZNQoOPhjGjIHy8qR64LeXC+MQ5pWxIrkmr0o3ObGX+L//DdddB/vuC3PnwpNP+m0L+vdPS5KHHBkHEUmbvEr0oe7m+OADv7lYz56wejXMng0TJ0Lv3mm/VKjHQUTSLq8SfSi7Od56Cy64wCf0XXaBhQvhnnvggAMydslQjoOIZExeJfpQdXPMng2nnQYnnggHHeQ7aCoroUuXjF86VOMgIhmXVx/GBt7N4RzMmOFvtF1b67cJnjgx7h0k0yXwcRCRrMqrRA8BdXNs2eIT+siRvi1y2DA480zYObjhV1eLSP7Iu0SfVRs2wAMPwM03wze+4XvhBwxIW/eMiEg8lOgz4eOP4a674M9/9itZH3sMjjwyrlPDvJBJRHKTEn06rVwJo0fD/ffDD34As2b5D1rjlAsLmUQk9+RV103GvPsu/Oxnfu/3LVtg/nx/Z6cEkjxoIZOIZEbezejTWhqZN89/wFpVBZde6u/D+rWvJR2bFjKJSCbkVaJPS2nEOXjuOZ/gFy+GK6/0H7imYQdJ3V9VRDIhr0o3KZVGGhpg8mT6/OIXfquC886DJUvgiivStk2wFjKJSCZkLNGb2UlmttjMlpjZsExdJxFJlUY++wzuuw8OPBBuvpn3zz/fb1sweDC0b5/W+CpKSxgxsBclRYUYUFJUyIiBvfRBrIikJCOlGzMrAO4ETgRWAK+Z2ZPOubczcb14JVQaWb/e7zkzejQceqjfMviYY/jviy+mvBd8a7SQSUTSLVMZ6zvAEufcUufc58AE4PQMXStucZVG1qyBa6+FHj2guhqmT4enn4Zjj9VCJxHJSeacS/+Tmp0JnOScuzj2/fnA4c65y5o8ZggwBKC4uLjPhAkT0h5Hc+o2bmb1uk3pwDNbAAAEkUlEQVR83rCV9gU7Udy5A0WF7eiwahXdJk1ir+efZ03//nx49tlsKtlxZl1fX0+niN26L1EaA0/j4GkcvCDGoX///tXOubK2Hpeprpvmpr7b/EZxzo0FxgKUlZW58vLyDIXShpoa30Hz9NMwZAiMGUPJ3nvTUvGkqqqKwGINCY2Bp3HwNA5emMchU6WbFUC3Jt93BVZm6FrJeeUVOPVUv/fMoYfC0qUwYgTsvXfQkYmIpFWmZvSvAQeYWQ+gFhgE/DhD14rf1q1+m+DKSvjoI7j6apg8GTp0CDoyEZGMyUiid85tMbPLgGeBAuB+59zCTFwrLps3w4QJvkTTvr3fJviMM6CgoO1zRURyXMZWxjrnZgAzMvX8cdmwAcaN89sE77cfjBrl7+ik7hkRySPR3AJh7Vq480644w7o1w8mTYLDDw86KhGRQERrC4QVK/zeM/vvD8uW+c3GpkxRkheRvBaNRP/OO3DRRb57xjlYsMDvCX/ggUFHJiISuNwu3axd6/eBf/lluOyylLcJFhGJotxO9J07w8knw8MPQ8eOQUcjIhJKuZ3oCwrg4ouDjkJEJNSiUaMXEZEWKdGLiEScEr2ISMQp0YuIRJwSvYhIxCnRi4hEnBK9iEjEKdGLiERcRu4Zm3AQZv8GPgg6jjjtAfwn6CACpjHwNA6exsELYhy+4Zzbs60HhSLR5xIzmxfPzXijTGPgaRw8jYMX5nFQ6UZEJOKU6EVEIk6JPnFjgw4gBDQGnsbB0zh4oR0H1ehFRCJOM3oRkYhToo+Tmb1vZjVm9oaZzQs6nmwxs/vNbI2ZvdXk2O5mNtPM3ov9uVuQMWZDC+NwvZnVxl4Tb5jZ94OMMdPMrJuZvWBmi8xsoZldETueV6+HVsYhtK8HlW7iZGbvA2XOubzqFzazY4B64GHn3CGxYzcCa51zlWY2DNjNOXdNkHFmWgvjcD1Q75y7OcjYssXMugBdnHOvm9lXgWqgAriQPHo9tDIOZxPS14Nm9NIq59xLwNrtDp8OPBT7+iH8izzSWhiHvOKcW+Wcez329XpgEVBCnr0eWhmH0FKij58D/m5m1WY2JOhgAlbsnFsF/kUP7BVwPEG6zMzejJV2Il2yaMrMugOlwFzy+PWw3ThASF8PSvTx6+ec6w2cDFwaeysv+W0MsB/wbWAVcEuw4WSHmXUCHgd+5Zz7JOh4gtLMOIT29aBEHyfn3MrYn2uAJ4DvBBtRoFbH6pSN9co1AccTCOfcaudcg3NuK3AvefCaMLN2+OT2qHNuSuxw3r0emhuHML8elOjjYGYdYx+6YGYdge8Bb7V+VqQ9CQyOfT0YmBZgLIFpTG4xPyTirwkzM2AcsMg5N6rJj/Lq9dDSOIT59aCumziY2b74WTzAzsBjzrkbAgwpa8xsPFCO35lvNfA7YCowCfg6sBw4yzkX6Q8qWxiHcvzbdAe8D/y8sVYdRWZ2FPAyUANsjR3+Db4+nTevh1bG4RxC+npQohcRiTiVbkREIk6JXkQk4pToRUQiToleRCTilOhFRCJOiV5EJOKU6EVEIk6JXkQk4v4/YdIMk2N9lqIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108c36e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(x, y)\n",
    "x1, x2 = ax.get_xlim()\n",
    "y1, y2 = np.dot(theta, [1, x1]), np.dot(theta, [1, x2])\n",
    "ax.plot([x1, x2], [y1, y2], linewidth=1, c='red')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-17.579094890510959, 3.9324087591240882)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sanity check\n",
    "theta_1 = np.dot((x - np.mean(x)), (y - np.mean(y))) / np.sum((x - np.mean(x))**2)\n",
    "theta_0 = np.mean(y) - theta_1 * np.mean(x)\n",
    "\n",
    "theta_0, theta_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is linear regression with multiple predictor variables.\n",
    "\n",
    "Model:\n",
    "    \n",
    "$y = \\theta_0 + \\theta_1 x_1 + \\theta_2 x_2 + \\theta_3 x_3 \\ldots + \\theta_n x_n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "wine = pd.read_csv('winequality-red.csv', sep=';')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1599, 12)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wine.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "wine_train, wine_test = wine[:1200], wine[1200:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.ones((len(wine_train), 1)) # models the intercept\n",
    "X = wine_train.drop('quality', axis=1).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.concatenate([x0, X], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = wine_train['quality'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  4.26148353e+01,   4.68226978e-02,  -1.06858847e+00,\n",
       "        -2.69078072e-01,   3.46373237e-02,  -1.60559854e+00,\n",
       "         3.48715209e-03,  -3.77404055e-03,  -3.91307833e+01,\n",
       "        -2.53835069e-01,   7.71874679e-01,   2.69457226e-01])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.ones((len(wine_test), 1)) # models the intercept\n",
    "x_test = wine_test.drop('quality', axis=1).values\n",
    "X = np.concatenate([x0, x_test], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = wine_test['quality']\n",
    "predictions = np.matmul(X, theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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>expected</th>\n",
       "      <th>predictions</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1205</th>\n",
       "      <td>7</td>\n",
       "      <td>6.123699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1206</th>\n",
       "      <td>7</td>\n",
       "      <td>6.123699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1207</th>\n",
       "      <td>5</td>\n",
       "      <td>5.425015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1208</th>\n",
       "      <td>7</td>\n",
       "      <td>6.123699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1209</th>\n",
       "      <td>7</td>\n",
       "      <td>6.192969</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1210</th>\n",
       "      <td>6</td>\n",
       "      <td>5.648277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1211</th>\n",
       "      <td>5</td>\n",
       "      <td>5.419992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1212</th>\n",
       "      <td>6</td>\n",
       "      <td>5.648277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1213</th>\n",
       "      <td>6</td>\n",
       "      <td>6.076723</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1214</th>\n",
       "      <td>6</td>\n",
       "      <td>5.939480</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      expected  predictions\n",
       "1205         7     6.123699\n",
       "1206         7     6.123699\n",
       "1207         5     5.425015\n",
       "1208         7     6.123699\n",
       "1209         7     6.192969\n",
       "1210         6     5.648277\n",
       "1211         5     5.419992\n",
       "1212         6     5.648277\n",
       "1213         6     6.076723\n",
       "1214         6     5.939480"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'expected': expected, 'predictions': predictions})[5:15]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the regression model doesn't do too badly on some of the records, though we would need a more formal way to measure this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split our data into train and test sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.ones((len(wine_train), 1)) # models the intercept\n",
    "X = wine_train.drop('quality', axis=1).values\n",
    "X = np.concatenate([x0, X], axis=1)\n",
    "\n",
    "train_y = wine_train['quality'].values\n",
    "train_y_hat = np.matmul(X, theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.ones((len(wine_test), 1)) # models the intercept\n",
    "X = wine_test.drop('quality', axis=1).values\n",
    "X = np.concatenate([x0, X], axis=1)\n",
    "\n",
    "test_y = wine_test['quality'].values\n",
    "test_y_hat = np.matmul(X, theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A good benchmark is to check the performance against the most naive model. In the regression case, we can use a model that always predicts the mean response as a benchmark:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "benchmark_y_hat = np.ones(len(test_y)) * np.mean(test_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Residual Sum of Squares (RSS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula:\n",
    "    \n",
    "$$RSS = \\sum_{i = 1}^{n} (y_i - \\hat{y_i})^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rss(y_true, y_pred):\n",
    "    return np.dot(y_true - y_pred, y_true - y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "488.86658602015399"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_rss = rss(train_y, train_y_hat)\n",
    "train_rss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "183.08042389268061"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_rss = rss(test_y, test_y_hat)\n",
    "test_rss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "252.79699248120298"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rss(test_y, benchmark_y_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Residual Standard Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RSE is an estimate of the standard deviation of $\\epsilon$. Roughly speaking, it is the average amount that the response will deviate from the true regression line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula:\n",
    "\n",
    "$$RSE = \\sqrt{\\frac{RSS}{n - 2}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rse(y_true, y_pred):\n",
    "    return np.sqrt(rss(y_true, y_pred) / (len(y_true) - 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.63880273682707245"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_rse = rse(train_y, train_y_hat)\n",
    "train_rse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6790874449701394"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_rse = rse(test_y, test_y_hat)\n",
    "test_rse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.79797759556625414"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rse(test_y, benchmark_y_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean Squared Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse(y_true, y_pred):\n",
    "    return rss(y_true, y_pred) / len(y_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.40738882168346163"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_mse = mse(train_y, train_y_hat)\n",
    "train_mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.45884818018215695"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_mse = mse(test_y, test_y_hat)\n",
    "test_mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.63357642225865407"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mse(test_y, benchmark_y_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### R-Squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula:\n",
    "\n",
    "$$R^2 = 1 - \\frac{RSS}{TSS}$$\n",
    "\n",
    "where \n",
    "\n",
    "$$TSS = \\sum_{i = 1}^{n} (y_i - \\overline{y})^2$$\n",
    "\n",
    "It is a measure of the proportion of variance of the residuals explained by the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tss(y_true):\n",
    "    t = y_true - np.mean(y_true)\n",
    "    return np.dot(t, t)\n",
    "\n",
    "def r_squared(y_true, y_pred):\n",
    "    t = tss(y_true)\n",
    "    r = rss(y_true, y_pred)\n",
    "    return 1 - (r/t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.37750170499006286"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_r_squared = r_squared(train_y, train_y_hat)\n",
    "train_r_squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adjusted R-Squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A drawback of $R^2$ is that by adding more parameters to the model, the $R^2$ value will always increase. A variant of this statistic called adjusted $R^2$ penalizes the model for over-parameterizing the data.\n",
    "\n",
    "The adjusted R-squared increases only if the new term improves the model more than would be expected by chance. It decreases when a predictor improves the model by less than expected by chance. The adjusted R-squared can be negative, but it’s usually not.  It is always lower than the R-squared."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formula:\n",
    "\n",
    "$$Adjusted\\_R^2 = 1 - \\frac{(1 - R^2)(N - 1)}{N - p - 1}$$\n",
    "\n",
    "where N is the total number of samples, and p is the number of parameters in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adjusted_r_squared(y_true, y_pred, p):\n",
    "    n = len(y_true)\n",
    "    r_sq = r_squared(y_true, y_pred)\n",
    "    return 1 - (1 - r_sq)*(n - 1)/(n - p - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.37173783188811904"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = wine_train.shape[1] - 1 # #predictors\n",
    "train_adj_r_squared = adjusted_r_squared(train_y, train_y_hat, p)\n",
    "train_adj_r_squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Bias-Variance Tradeoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For cases where we have more than ~10k features (extremely common if you're considering polynomial regression, or including interaction terms), we can no longer use the Normal equation form. For such cases, we can use gradient descent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In gradient descent, we move in the opposite direction of the gradient. For least-squares regression (what we've been doing till now), the gradient can be derived as follows:\n",
    "\n",
    "$J(\\theta) = \\frac{1}{n} \\sum_{i = 1}^{n} (\\hat{y_i} - y_i)^2 $\n",
    "\n",
    "$J(\\theta) = \\frac{1}{n} \\sum_{i = 1}^{n} (h_{\\theta}(x_i) - y_i)^2 $\n",
    "\n",
    "where \n",
    "\n",
    "$h_{\\theta}(x_i) = \\theta^T X$\n",
    "\n",
    "We have:\n",
    "\n",
    "$\\frac{dh}{d\\theta_j} = X_j$, where $X_j$ is the jth feature.\n",
    "\n",
    "Thus we get the partial derivative of J w.r.t. $\\theta_j$ as \n",
    "\n",
    "$\\frac{d J(\\theta)}{d\\theta_j} = \\frac{1}{n} \\sum_{i = 1}^{n} 2(h_{\\theta}(x_i) - y_i)X_{ij} $, \n",
    "\n",
    "where $X_{ij}$ is the jth feature in the ith example.\n",
    "\n",
    "The gradient is just the vector of all the partial derivatives. Also remember that our $X_0$ was just a dummy variable, and is always 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "cars = np.genfromtxt('cars.csv', delimiter=',', names=True)\n",
    "\n",
    "x0 = np.ones(len(cars)) # Models the intercept\n",
    "x = cars['speed']\n",
    "\n",
    "X = np.array([x0, x]).T\n",
    "y = cars['dist']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_gradient(X, y, y_pred):\n",
    "    n = len(y_pred)\n",
    "    return (2/n) * np.matmul((y_pred - y), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_descent(X, y, learning_rate=0.001, iters=20000):\n",
    "    # Initialize theta to all zeros\n",
    "    theta = np.zeros(X.shape[1])\n",
    "\n",
    "    for i in range(iters):\n",
    "        # Predict\n",
    "        y_pred = np.matmul(X, theta)\n",
    "        \n",
    "        # Calculate gradient\n",
    "        theta_grad = evaluate_gradient(X, y, y_pred)\n",
    "        \n",
    "        # Update weights\n",
    "        theta = theta - learning_rate*theta_grad\n",
    "\n",
    "    return theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "theta = gradient_descent(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-17.2934172 ,   3.91577301])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt4VdW19/HvEFGjVCNFUxtRKKB9kRxFouJBbcD2IF5KRG1Ba6mXom2l1VZKqPVSaw9BvFSs0lJv2AJRAVEriijGCxV6wkW5iaBUIKIoCiWACsl8/5g7GiAh+77WXvl9nseHZO21soeTzcjcY485lznnEBGR6Nor6ABERCSzlOhFRCJOiV5EJOKU6EVEIk6JXkQk4pToRUQiToleRCTilOhFRCJOiV5EJOL2DjoAgHbt2rkOHToEHUZctmzZwgEHHBB0GIHSGHgaB0/j4AUxDvPmzfvIOXdIc+eFItF36NCBqqqqoMOIS2VlJSUlJUGHESiNgadx8DQOXhDjYGbvxnOeSjciIhGnRC8iEnFK9CIiEadELyIScUr0IiIR12yiN7MHzGy9mS1ucGy0mb1pZm+Y2eNmlt/gsRFmttLMlptZ30wFLiKSbtMWVNOrfBYdy56mV/kspi2oDjqktIhnRv8QcMYux2YC3Zxz/wW8BYwAMLOuwEDgmNg195pZq7RFKyKSIdMWVDNi6iKqN27DAdUbtzFi6qJIJPtmE71z7mXg412OPeec2xH7dg5weOzr/kCFc+4z59wqYCVwYhrjFRHJiNEzlrNte+1Ox7Ztr2X0jOUBRZQ+6VgwdSnwSOzrQnzir7c2dmw3ZjYEGAJQUFBAZWVlGkLJvJqampyJNVM0Bp7GwYvKOAxsvxnaN/bI5rj+/8I8DiklejO7DtgBTKg/1Mhpjd593Dk3DhgHUFxc7HJlZZ1WAWoM6mkcvKiMw3Xls6jeuG2344X5eQy9qKTZ68M8Dkl33ZjZYOBs4CLnXH0yX8vOvxMPB95LPjwRkewY1vdo8lrv/JFiXutWDOt7dEARpU9Sid7MzgCGA991zm1t8NCTwEAz29fMOgJdgH+lHqaISGaVdi9k5IAiCvPzMPxMfuSAIkq7N1p9zinNlm7MbBJQArQzs7XAjfgum32BmWYGMMc5d6VzbomZPQosxZd0fuacq238J4uIhEtp98JIJPZdNZvonXODGjl8/x7O/wPwh1SCEhGR9NHKWBGRiFOiFxGJOCV6EZGIU6IXEYk4JXoRkYhTohcRiTglehGRiFOiFxGJOCV6EZGIU6IXEYk4JXoRkYhTohcRiTglehGRiFOiFxGJOCV6EZGIU6IXEYk4JXoRkYhTohcRiTglehGRiFOiFxGJOCV6EZGIU6IXEYm4ZhO9mT1gZuvNbHGDY23NbKaZrYj9eXDsuJnZGDNbaWZvmNnxmQxeRESaF8+M/iHgjF2OlQEvOOe6AC/EvgfoB3SJ/TcEGJueMEVEJFnNJnrn3MvAx7sc7g+Mj309HihtcPxh580B8s3ssHQFKyIiiUu2Rl/gnFsHEPvz0NjxQmBNg/PWxo6JiEhA9k7zz7NGjrlGTzQbgi/vUFBQQGVlZZpDyYyampqciTVTNAaexsHTOHhhHodkE/0HZnaYc25drDSzPnZ8LdC+wXmHA+819gOcc+OAcQDFxcWupKQkyVCyq7KyklyJNVM0Bp7GwdM4eGEeh2RLN08Cg2NfDwaeaHD8h7Hum57ApvoSj4iIBKPZGb2ZTQJKgHZmtha4ESgHHjWzy4DVwAWx06cDZwIrga3AJRmIWUREEtBsonfODWriodMbOdcBP0s1KBERSR+tjBURiTglehGRiFOiFxGJuHT30YuI5KxpC6oZPWM5723cxtfz8xjW92hKu+f+mk8lehERfJIfNvl1ttf6NZ7VG7cxbPLrADmf7FW6EREBfvfUki+SfL3ttY7fPbUkoIjSR4leRAT4ZOv2hI7nEiV6EZGIU6IXEQHy81ondDyXKNGLiAA3ffcYWu+18wa8rfcybvruMQFFlD7quhER4cvOGrVXiohEWGn3wkgk9l2pdCMiEnFK9CIiEadELyIScUr0IiIpmLagml7ls1hUvYle5bOYtqA66JB2ow9jRUSSNG1BNSOmLmLb9lpo7/fHGTF1ERCu/XGU6EUiJqo7MIbR6BnLfZJvYNv2WkbPWB7/mDsHZs2flwKVbkQipH6GWb1xG44vZ5hhLCdEwXsbtyV0/At1dfCPf8App8DLL2cgsp0p0YtEyJ5mmJJ+X8/PS+g427fD3/4G//VfcMMNMHQo9OqVwQg9JXqRCEl6hilJGdb3aPJat9rpWF7rVgzre/TOJ27ZAmPGQOfO8NBDcMcdMG8efP/7sHfmK+iq0YtEyNfz86huJKk3OcOUlDTcNgE2U7jrZyIbNsCf/gT33AOnngqPPQYnnpj1OFOa0ZvZNWa2xMwWm9kkM9vPzDqa2VwzW2Fmj5jZPukKVkT2LO4ZpqRNafdCZpf1oajwIGaX9fFJfs0auOYa6NLFf/3KKzBlSiBJHlJI9GZWCPwcKHbOdQNaAQOBUcCdzrkuwCfAZekIVCQT6nugO5Y9Hdoe6ESUdi9k5IAiCvPzMKAwP4+RA4rUdZMtS5fCj34Exx3nSzKLFsF998HRwf6iTbV0szeQZ2bbgf2BdUAf4MLY4+OBm4CxKT6PSNrt1ANNeHugExXVjblC7bXX6HbddfD22/4D1pUr4eCDg47qC0nP6J1z1cBtwGp8gt8EzAM2Oud2xE5bC+gVJ6GkDhVJiXMwfTp861tw4YV8fMIJsGoVXHddqJI8gDnnmj+rsQvNDgamAN8HNgKPxb6/0TnXOXZOe2C6c66okeuHAEMACgoKelRUVCQVR7bV1NTQpk2boMMIVFTGYFH1piYfKyo8qNnrozIOqWpp42C1tRzy4oscMWkSAKsvvJAPS0rYvG1b1sehd+/e85xzxc2dl0rp5tvAKufchwBmNhX4byDfzPaOzeoPB95r7GLn3DhgHEBxcbErKSlJIZTsqaysJFdizZSojMF15bMa7VApzM9j6EUlzV4flXFIVYsZh61b4cEH4bbb4MgjfTfNGWfQNbaqNczjkErXzWqgp5ntb2YGnA4sBV4Ezo+dMxh4IrUQRTJDHSoSl08+gVtugY4dYeZMmDgRKiuhX7+Mb12QLknP6J1zc81sMjAf2AEswM/QnwYqzOyW2LH70xGoSLpF9dZx2usmTdauhTvv9LP4/v3hxReha9ego0pKSl03zrkbgRt3OfwOEEyzqEiCotahEtVOoqx6800YPRoef9y3Sr7+OrRvH3RUKdEWCCIRok6iFMydCwMGwGmnQYcOvkXyjjtyPsmDtkAQiRTtdZMg5+C556C8HN55B6691m86dsABQUeWVkr0IhGivW7itGMHTJ4Mo0b5r4cP9xuMtW4ddGQZodKNSISok6gZ27bBn//styS45x7fTfPGG/CDH0Q2yYNm9CKREtVOopRt3Ahjx/qtgouLYfx4f9OPFkKJXiRiotZJlJL33oM//hHuvx/OOsv3wXfrFnRUWafSjYhEz4oVMGSIT+qffQbz58PDD7fIJA+a0YtIlFRV+Q9YKyvhZz+Dt96Cdu2CjipwSvQiktucgxde8C2Sb70Fv/ylX82axAZjUV1VrEQvIrmpthamTvUz+K1bfYvkoEGwT3I3tYvyqmIlehHJLZ9+6uvto0f7sswNN8DZZ8NeqX3kuKdVxUr0IpLzcqJk8Z//+B74P/4Runf3nTSnnpq2HSSjvKpYiV6khQt9yeL99+Guu+Cvf4W+feGZZ+DYY9P+NFFeVaz2SpEWLrQboa1cCVde6bcG3rwZ/u//YMKEjCR5iPaqYiV6kRYudCWLBQtg4EDo2RMOOYTpj71Ir8MH0PEvS+lVPotpC6oz8rSl3QsZOaCIwvw8DH+nsZEDisLxriZFKt2ItHChKFk453vfy8thyRK45hr461+ZtvI/WS0rRXVVsWb0Ii1coCWLujrfInnSSfCTn/gdJN9+G371K/jKV8JbVsoxSvQiLVyyJYtpC6rpVT6LRdWbEi+pfPaZ75rp2tX3wY8YAUuXwqWXwr77fnFa6MpKOUqlGxFJuGSxU6dO+wRKKps3w7hx/l6s3br5HSVLSppskQxFWSkCNKMXkYQlXFJZvx5++1vo2NHvR/PUU/Dss9C79x774KPcCZNNSvQikrC4SyqrVvnNxb75Tdiwwd+XddIkv+ApDlHuhMkmlW5EJGHNllRef93X3p97zm8XvGwZFBQk9VxR7YTJJs3oRSRhjZZU9t6L8nafwJlnQr9+ftb+zjvwv/+bdJKX9EhpRm9m+cB9QDfAAZcCy4FHgA7Av4HvOec+SSlKEQmVhrcspG4T339vPsMXTqPtp5th2DDfMrnffgFHKfVSLd3cBTzrnDvfzPYB9gd+A7zgnCs3szKgDBie4vOIZERObOYVUqXHHELp6zPZcv1NHNCuHdwwAs49F1q1av5iyaqkSzdmdiBwGnA/gHPuc+fcRqA/MD522nigNNUgRTKhvkWweuM2HF+2CGZqiX1k1NT4HSQ7dYIJE1gxdKjfh+b885XkQyqVGv03gA+BB81sgZndZ2YHAAXOuXUAsT8PTUOcImmnVZcJ+ugjuPFG3yL5z3/CtGkwcyYbe/RI21bBkhnmnEvuQrNiYA7Qyzk318zuAv4DDHXO5Tc47xPn3MGNXD8EGAJQUFDQo6KiIqk4sq2mpoY2SdyiLEqiMgaLqjc1+VhR4UHNXh+VcWjOvu+/T/tHH6Xg+ef58LTTWDNwINsOP/yLx1vKODQniHHo3bv3POdccXPnpZLovwbMcc51iH1/Kr4e3xkocc6tM7PDgErn3B5XNxQXF7uqqqqk4si2yspKSkpKgg4jUFEZg17lsxptESzMz2N2WZ9mr4/KODRp0SK49VaYPh0uvxyuvhoOO2y30yI/DnEKYhzMLK5En3Tpxjn3PrDGzOqT+OnAUuBJYHDs2GDgiWSfQySTtOqyCa++CuecA//zP3DMMX6TsVGjGk3ykhtS7boZCkyIddy8A1yC/+XxqJldBqwGLkjxOUQyomGLYJS6bpLqJKqr8zP38nJYtw5+/Wt47LE9tkjWP8/A9pu5rnxWJMYuqlJK9M65hUBjbxtOT+XnimRL1FZdJnxbwO3boaLCz9j32QeGD4fzzoO995wakt7UTAKhlbEiERJ3J9GWLTBmDHTuDA8+CHfcAfPm+f3gm0nyCT2PhIISvUiCUtqHPcOa3Wxswwa4+WbfIvnSS748M2uWr8cn0CIZ9n3i6/+OOpY9Hbq/oyAo0YskoOEiKwjfIqum9mk/js3+9nxdusDq1fDKKzBlCpx4YlqfJwz7xGsh3O6U6EUSEPaSxa6dRJ0/Ws2dz/yRR8Ze6VetLloE990HR6fWWRTmjqWw/x0FQdsUiyQg7CWL+g9Cp4+byvkvTKR43Zu8/8Mh7FP5KBy827rFlJ6n6t2PmTR3DQCtzDivRzg+2A7731EQlOhFEhDqW9s5B888Q+moUZSuXu13kbzkEtrmpT+2aQuqmTKvmtrYgsta55gyr5riI9sGnuxD/XcUEJVuRBIQypLFjh0wcSIcd5y/yfaVV8KKFfDTn0IGkjyEuzwSyr+jgGlGL5KAnfZhZzOFGV5ktcfFT1u3+tbI226DI47wi53OOCMrG4yFuTwS1YVwqVCiF0lQ/SKryspKhl5UkrHnaWrxU+v/bOSsVx6HP/0Jevb0s/mTT85YHI0Je3kkagvhUqXSjUhI7VoeKdj8EdfMGMep/U72+8/MmuW3Cs5ykgeVR3KNEr1ISNWXQTptWMOo6Xcx44GraOXqOONHY3zJpmvXwGIr7V7IeT0KaRUrE4Wp60Z2p9KNSEidvvldzn/+7/RYu4y/HX8WJUPGsTHvQApDUB4Jc9eN7E4zepEwcQ5mzIA+fRgz9X+p6ngsp11xH2N6DWJj3oGhKY+EuetGdqcZvUgY7NgBkyf7XSS3b4eyMvb//vfptng9bUPYPRLmrhvZnRK9REZS+7AH7dNP4aGHYPRo+PrX4fe/hzPPhL38m+1sdY8kOnZh77qRnal0I5GQcxtZbdwII0f6XSSffhrGj/cbjZ199hdJPluSGTt13eQWJXqJhJypGb/3nr97U6dOsGwZzJwJTz0Fp5wSWEjJjF1p90JGDij64oPhwvw8Rg4oCv87qBZKpRuJhNDXjFes8OWZyZPh4oth/nw48sigowKSH7tsLRyT1GlGL5EQ2v3Rq6rgggvgv//b31z7rbfgrrtCk+QhxGMnaaNEL5EQqpqxc/D88/Dtb8O550KvXrBqFfzud9CuXfbjaUaoxk4yQqUbiYRQbGRVWwtTp/oWya1b/Y22Bw3yN93OokQ7aEIxdpJRSvQSGYFtZPXZZ/Dww3DrrX7Gfv31cM45We+egaY3QgOaTfZK7NGV8ivRzFqZ2QIz+0fs+45mNtfMVpjZI2aW3emMSJa02rLFJ/eOHf3mYvffD//8J/TvH0iShxzqPpKsSseM/hfAMuDA2PejgDudcxVm9mfgMmBsGp5HJBzefx/uuoue997r+96feQaOPXaPlyS7mCvR60LffSSBSGnaYWaHA2cB98W+N6APMDl2ynigNJXnEAmNt9+Gn/zE7xq5eTPz/vxnmDAhriSfzGKuZK5TB400JtX3l38Efg3Uxb7/KrDRObcj9v1aQIU/yW0LFsDAgXDSSb4G/+ab8Kc/8elhh8V1ebLllGSuUweNNMZcbJvRhC80Oxs40zn3UzMrAa4FLgFec851jp3THpjunCtq5PohwBCAgoKCHhUVFcn9H2RZTU0Nbdq0CTqMQLWIMXCO/IULOWLiRA74979Zc8EFrDv7bGr33/+LU+Idh0XVm5p8rKjwoLRft3Hbdj7Y9Cmf19axT6u9KDhoP/LzWjcbZ7JaxOshDkGMQ+/evec554qbOy+VRD8SuBjYAeyHr9E/DvQFvuac22FmJwM3Oef67ulnFRcXu6qqqqTiyLbKykpKSkqCDiNQkR6Dujr/weqoUbBpk9+u4KKLYN99dzs13nHoVT6r0Q3ACvPzmF3WJ+3XZVukXw8JCGIczCyuRJ906cY5N8I5d7hzrgMwEJjlnLsIeBE4P3baYOCJZJ9DJGs++wweeMDX30eNgrIyWLoULr200SSfiGTLKSrDSLpkoo9+OFBhZrcAC4D7M/AcImnxj1eXs2rknZz/8mTePewb7Lj2Zk758QUQu0VeOiS7ICmbC5lycotniVtaEr1zrhKojH39DnBiOn6uSMasX8/y3/yBXhMfwh15LJeffwNLCjqRt6YVIxe+l/Ykl+yCpGwsZEp2kZXkDu11Iy3LqlVw1VXwzW+ydPEq+l98O0P7D2dJQSegZS4u0iKr6FOil7hNW1BNr/JZLKreRK/yWaG7qcdvpy2i04jpdCh7mk4jpvPbaYu+fPCNN/yHqiecAAceCEuX8suSK1h98O4tki1tcZEWWUWf9rqRuOz09r59+N7e/3baIv4+Z/UX39c6x99fe5cOS6q4fPZjsHAhXH013HsvHORbE3U7PE/jEH2a0Utcwv72ftLcNV98ba6O76yYw9S/X8u3b/+t3yr4nXd8q+RBX/afq6vF0zhEn2b0Epewv72vdY7Wtdvpv/Qlrpg7hW2t92XsSecz46iTeefH3230Gm3P62kcok+JXuIS6rf3NTVcVvUEl/5rGu+0LeTGb1/BP488Fsxo1UybZLa2583WpmbJSmYc6mMb2H4z15XP0i+HEFOil7gM63v0Ti14EIK39x99BHffDffeS+lRx3PFgOtY/LXOO50y6KT2AQX3pWTbF8Pc9hj2z2xkZ6rRS1xKuxcyckARhbEZfGF+HiMHFAXzj/rdd+HnP4ejjoJ162D2bIpmz+C40j5fzOBbmfGDnkdwS+lu2yxlXTY3NcuWMMcmu9OMXuJW//a+srKSoReVZD+AxYv9jT6efhouvxyWLPE33I65pbQoFIl9V8l+vhHmz0XCHJvsTjN6Cb9XX/W35vvOd/xeNG+/7fejiXOb4KAlu0d8mPeWD3Nssjslegmnujr4xz/glFNg8GA46yzfIllWBvn5QUeXkChuahbm2GR3Kt1IxiXUObJ9O1RUsOnmP/DBllrGnDCAN/p9m1+e1JXSvNycLebCpmaJahgbbKYwRLHJ7pToJaPi7hzZssXfXPv22/nw0MMpO+FiXmjvWyTZ/HnOd3SEeVOzZAX+mY3ETaUbyahmuzM2bICbb4aOHeGll+Cxxyg97/e8cMRxO20VrI4OkeQp0UvcktnUrKkuDLd6NVxzDXTp4tslX34ZpkyBE09MuqOjPr6OZU/HHV8q14R1czeRXal0I3FJdoHMritqO3+0mivnTuE7b/8LfjoEFi2CwsI9XtPweFzxEV98KV+jhUKSIzSjl7gku0Cmvjvj+OpljJt6C5Mm/Ya17Qp59ZnX4LbbdkvyDa9pqLmOjmTiy9Y1IkHTjF7iklQ5xTlK33+DU6b/js9Xvcufi0sp/8H1/PycYzlrD7Pf0u6FVL37MZPmrqHWOVqZcV6PPX8omUx82bpGJGhK9BKXhMopO3bAo4/6RU1Au+HD4Xvf4+a943u5TVtQzZR51dQ6B/idKafMq6b4yLZNJvtkyj3ZukYkaCrdSFziKqds3Qr33OM/YP3LX6C83N/w48ILIc4kD8mVR5Ip92TrGpGgaUYfMtnaljZRe1wg88kn/s5Nd98NPXvCxIlw8slJP1cy5ZFkFheleo0WCkmuUKIPkTBvS1sfw04LZKqr4dpr4YEHoH9/mDXL70WTomTLI8ksLkrlGi0Uklyh0k2I5EpHx/6rV8Nll0FREdTWwuuvw4MPpiXJg8ojIumWdKI3s/Zm9qKZLTOzJWb2i9jxtmY208xWxP48OH3hRlvoOzr+9S8YMIDjfvELOPJIWLEC7rwT2qf35h6l3Qs5r0fhTnvLN9d1k01aMCW5JpUZ/Q7gV865/wf0BH5mZl2BMuAF51wX4IXY9xKHUG796hw89xz06QMXXAC9ezNn4kS44Qb46lcz8pRNdd2EIaHWl9fqS0v15bUwxCbSlKQTvXNunXNufuzrzcAyoBDoD4yPnTYeKE01yJYiVCWLHTvgkUegRw/45S/hkktg5UoYOpS6DO8iGeYSVphjE2mKudisKaUfYtYBeBnoBqx2zuU3eOwT59xu5RszGwIMASgoKOhRUVGRchzZUFNTQ5s2bTL28zdu284Hmz7l89o69mm1FwUH7Ud+XuuMPd+u9vr8c7727LO0f+QRPm/bltUXXsiGk06Cvb6cE2R6DBZVb2rysaLCgzL2vPFoGFtBHnzQoKoWdGxByfTrIVcEMQ69e/ee55wrbu68lBO9mbUBXgL+4JybamYb40n0DRUXF7uqqqqU4siWyspKSkpKgg4j/TZtgrFj4a67oLgYhg/3N/1ooL71c2D7zVSs+UrG2gp7lc9qtOumMD+P2WV90v58iWgY26+KdnD7It+4FobYghLZfxMJCmIczCyuRJ9S142ZtQamABOcc1Njhz8ws8Nijx8GrE/lOSTD1q3zSf0b34ClS2HmTHjqqUaTfLK16UR3iAxVCWsXYY5NpCmpdN0YcD+wzDl3R4OHngQGx74eDDyRfHiSMStWwJAhcMwx8OmnMH8+PPwwdOvW6OnJ1qYb/oJwxPcLorR7ISMHFFGYn4fhZ8sjBxSFouumYWwQrthEmpLKgqlewMXAIjNbGDv2G6AceNTMLgNWAxekFqI0J6HVtFVVfg+aykr46U/hrbegXbtmnyPZ1s89/YJobgVqWJOnFkxJrkk60TvnXgWsiYdPT/bnSmLiWk3rHLzwgk/wb74Jv/qVX+CUwAdHya5WDf3aAJEWQCtjc9weSyq1tTB5MpxwAvz85/CDH8Dbb8PVVyeU5CH52nQo1waItDDa6ybHNTYz3mfHdk6rfBYeuNIvarr+ejjnnJ1aJBOV7GZew/oevdM7DojvF0RYN3cTyUVK9DmuYUmlzWdbuWjhdC6tepJ3vt4Z7rsPTj11p5tspyKZ2nQyO0SGfXM3kVyjRJ/jhvU9mtsefplBcx5n0OszeLljd64Y9Ht+9JP+EJKkmOgHq8l+gCsijVOiz2Vvv03puNs4c+IknupaQv8f3kFdh445X+bQB7gi6aVEn4sWLPAdNM8/D1deyT4r3uK8Qw/lvKDjShPdrk8kvdR1kyucgxdfhL594eyzfSfNqlVwyy1w6KFBR5dWWn0qkl4tbkafc90cdXUwbZqfwW/aBL/+NTz5JOy7b0o/NszjkMwHuCLStBaV6HOqm+Ozz2DCBLj1VjjwQBgxwt+uL4UWyXq5MA5hXhkrkmtaVOkmJ/YS37wZbr8dOnXy+8GPHQtz58K556YlyUOOjIOIpE2LmtGHuptj/XoYMwb+8hc4/XS/g2T37hl5qlCPg4ikXYua0YdyOf6qVXDVVfDNb8KGDTBnDlRUZCzJQ0jHQUQypkUl+lB1c7zxBlx0kb/Jx1e+4veCHzvWl2wyLFTjICIZ16JKN4F3czgHr7ziO2gWLPCbi917LxyU3VvQBT4OIpJVLSrRQ0DdHHV1vuY+ahR8+KFvkZwyBfbbL7txNKCuFpGWo8Ul+qz6/HOYONG3SO6/P5SV+e6ZVq2av1ZEJE2U6DOhpsbvHHnHHXDUUb6b5vTT07aLpIhIIpTo0+mjj+Duu/2Hqt/6Fkyd6j9sTUCYV6yKSG5qUV03GfPuu/4OTkcdBevWwauvwmOPJZXkE72RtohIc5ToU7F4Mfzwh3D88f6D1cWLYdw4n/CToBWrIpIJLa50k5bSyOzZUF4OVVV+Jj9mDOTnpxybVqyKSCa0qESf0mZedXUwfTrdR4yArVth2DB49FHIS99qUu3DLiKZkLHSjZmdYWbLzWylmZVl6nkSkVRpZPt2+Nvf4Nhj4frrWXvuubB8OVx5ZVqTPGjFqohkRkYSvZm1Au4B+gFdgUFm1jUTz5WIhEojW7f6DprOneHBB+G222D+fD7s0wf2zswbodLuhYwcUERhfh4GFObnMXJAkbpuRCQlmSrdnAisdM69A2BmFUB/YGmGni+Zny8SAAAE10lEQVQucZVGNmyAe+7x/51yiu+eOfHErMWoFasikm6ZKt0UAmsafL82dixQeyyNrFkD11wDXbr4dsmXXvLbFGQxyYuIZII559L/Q80uAPo65y6PfX8xcKJzbmiDc4YAQwAKCgp6VFRUpD2Oxmzctp0PNn3K57V17NNqLzptWs8x0ybTbvZs3u/Xj7Xnn89nhxzS5PU1NTW0adMmK7GGlcbA0zh4GgcviHHo3bv3POdcswt2MlW6WQu0b/D94cB7DU9wzo0DxgEUFxe7kpKSDIXShDlzfIvka6/B0KEwYQLt27bdKejGVFZWkvVYQ0Zj4GkcPI2DF+ZxyFSi/z+gi5l1BKqBgcCFGXqu+DkHzz7rE/zq1XDttX7Tsf33DzoyEZGMyUiid87tMLOrgBlAK+AB59ySTDxXXHbs8D3vo0b5ZF9WBt/7Xsa6Z0REwiRjmc45Nx2YnqmfH5dt23xr5OjRcMQRMHIk9OunXSRFpEWJ5pT2k0/8nZvuvht69vTlmZNPDjoqEZFARGtTs+pqX3fv1AlWrIBZs2DaNCV5EWnRopHoly+Hyy+HoiKorYWFC+Ghh6Br4ItxRUQCl9ulm48/hh//2N9w+6qr/Cz+q18NOioRkVDJ7UR/0EH+w9WHH4YDDgg6GhGRUMrtRN+qlS/ZiIhIk6JRoxcRkSYp0YuIRJwSvYhIxCnRi4hEnBK9iEjEKdGLiEScEr2ISMQp0YuIRFxGbiWYcBBmHwLvBh1HnNoBHwUdRMA0Bp7GwdM4eEGMw5HOuabvfRoTikSfS8ysKp57NEaZxsDTOHgaBy/M46DSjYhIxCnRi4hEnBJ94sYFHUAIaAw8jYOncfBCOw6q0YuIRJxm9CIiEadEHycz+7eZLTKzhWZWFXQ82WJmD5jZejNb3OBYWzObaWYrYn8eHGSM2dDEONxkZtWx18RCMzszyBgzzczam9mLZrbMzJaY2S9ix1vU62EP4xDa14NKN3Eys38Dxc65FtUvbGanATXAw865brFjtwIfO+fKzawMONg5NzzIODOtiXG4Cahxzt0WZGzZYmaHAYc55+ab2VeAeUAp8CNa0OthD+PwPUL6etCMXvbIOfcy8PEuh/sD42Nfj8e/yCOtiXFoUZxz65xz82NfbwaWAYW0sNfDHsYhtJTo4+eA58xsnpkNCTqYgBU459aBf9EDhwYcT5CuMrM3YqWdSJcsGjKzDkB3YC4t+PWwyzhASF8PSvTx6+WcOx7oB/ws9lZeWraxQCfgOGAdcHuw4WSHmbUBpgBXO+f+E3Q8QWlkHEL7elCij5Nz7r3Yn+uBx4ETg40oUB/E6pT19cr1AccTCOfcB865WudcHfBXWsBrwsxa45PbBOfc1NjhFvd6aGwcwvx6UKKPg5kdEPvQBTM7APgfYPGer4q0J4HBsa8HA08EGEtg6pNbzLlE/DVhZgbcDyxzzt3R4KEW9XpoahzC/HpQ100czOwb+Fk8wN7AROfcHwIMKWvMbBJQgt+Z7wPgRmAa8ChwBLAauMA5F+kPKpsYhxL823QH/Bu4or5WHUVmdgrwCrAIqIsd/g2+Pt1iXg97GIdBhPT1oEQvIhJxKt2IiEScEr2ISMQp0YuIRJwSvYhIxCnRi4hEnBK9iEjEKdGLiEScEr2ISMT9f7gOC3/zt3aSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111a5c2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(x, y)\n",
    "x1, x2 = ax.get_xlim()\n",
    "y1, y2 = np.dot(theta, [1, x1]), np.dot(theta, [1, x2])\n",
    "ax.plot([x1, x2], [y1, y2], linewidth=1, c='red')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we get the same regression coefficients as we got from the Normal equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactions and Polynomial Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Polynomial regression is when you model the relationship between the outcome and the predictor as an nth degree polynomial.\n",
    "\n",
    "Interaction terms basically account for interactions between variables. For example, instead of just using the `height` and `area` of a house to predict the price, you may also want to use `height*area`, i.e. people might prefer spacious houses, rather than just really tall or really vast houses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regularization prevents overfitting. In the context of regression, the two most common forms are:\n",
    "\n",
    "1. LASSO: Least Absolute Shrinkage Selector Operator) (L1)\n",
    "2. Ridge (L2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both cases, $\\lambda$ is a hyperparameter, which decides the extent to which you want to penalize large coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### LASSO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cost function:**\n",
    "\n",
    "$$J(\\theta) = \\frac{1}{n} \\left[ \\sum_{i = 1}^{n} (h_{\\theta}(x_i) - y_i)^2 + \\lambda \\sum_{j = 1}^{k} |\\theta_{j}| \\right]$$\n",
    "\n",
    "Note that we do not penalize the bias term ($\\theta_0$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Advantages**\n",
    "\n",
    "1. Feature Selection: The above form drives certain coefficients towards zero, and thus results in sparse coefficients.\n",
    "\n",
    "**Disadvantages**\n",
    "\n",
    "1. No analytical solution: Unlike the simple model, and the ridge model, there is no Normal equation to directly calculate the values of the weight vector $\\theta$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cost function:**\n",
    "\n",
    "\n",
    "$$J(\\theta) = \\frac{1}{n} \\left[ \\sum_{i = 1}^{n} (h_{\\theta}(x_i) - y_i)^2 + \\lambda \\sum_{j = 1}^{k} \\theta_{j}^2 \\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is also known as *Tikhonov Regularization*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For instance, imagine that you are trying to perform polynomial regression. Without regularization, it is plausible that your model will choose a more complex model, which fits the training data extremely well, but fails to generaize.\n",
    "\n",
    "In this case, you may want to introduce a regularization term (say L1), so your model is penalized for using too many coefficients. In this case, it will push the potentially non-essential coefficients to zero, resulting in a model that can generalize well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing coefficients as a function of $\\lambda$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LASSO**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](http://scikit-learn.org/stable/_images/sphx_glr_plot_lasso_lars_001.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ridge**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](http://scikit-learn.org/stable/_images/sphx_glr_plot_ridge_path_001.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The L1/L2 Norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Till now, we've been using the L2 norm, or least squares regression for our regression analysis, i.e. our loss function for a single example was defined as:\n",
    "\n",
    "$L(y_i, \\hat{y}_i) = (y_i - \\hat{y}_i)^2$\n",
    "\n",
    "However, in some cases, we may want to use the L1 norm, defined as\n",
    "\n",
    "$L(y_i, \\hat{y}_i) = |y_i - \\hat{y}_i|$\n",
    "\n",
    "One advantage of using the L1 norm is **robustness to outliers**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assumptions of the Linear Regression Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. A straight-line relationship between the predictors and the response\n",
    "2. Uncorrelated error terms\n",
    "3. Error terms have constant variance"
   ]
  }
 ],
 "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.6.4"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}