worksheet02a_wip.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COMP90051 Workshop 2\n",
    "## Part A: Linear regression\n",
    "\n",
    "***\n",
    "\n",
    "Our aim for this part of the workshop is to fit a linear model from scratch—relying only on the `numpy` library. We'll experiment with two implementations: one based on iterative updates (coordinate descent) and another based on linear algebra. Finally, to check the correctness of our implementation, we'll compare its output to the output of `sklearn`.\n",
    "\n",
    "Firstly we will import the relevant libraries (`numpy`, `matplotlib`, etc.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check what a command does simply type `object?`. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Review\n",
    "In lectures, we saw that a linear model can be expressed as:\n",
    "$$y = w_0 + \\sum_{j = 1}^{m} w_j x_j = \\mathbf{w} \\cdot \\mathbf{x} $$\n",
    "where \n",
    "\n",
    "* $y$ is the *target variable*;\n",
    "* $\\mathbf{x} = [x_1, \\ldots, x_m]$ is a vector of *features* (we define $x_0 = 1$); and\n",
    "* $\\mathbf{w} = [w_0, \\ldots, w_m]$ are the *weights*.\n",
    "\n",
    "To fit the model, we *minimise* the empirical risk with respect to $\\vec{w}$. In this simplest case (square loss), this amounts to minimising the sum of squared residuals:\n",
    "\n",
    "$$SSR(\\mathbf{w}) = \\sum_{i=1}^{n}(y_i - \\mathbf{w} \\cdot \\mathbf{x}_i)^2$$\n",
    "\n",
    "**Note:** For simplicity, we'll consider the case $m = 1$ (i.e. only one feature excluding the intercept)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Data set\n",
    "We'll be working with some data from the Olympics—the gold medal race times for marathon winners from 1896 to 2012. The code block below reads the data into a numpy array of floats, and prints the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1896.            4.47083333]\n",
      " [1900.            4.46472926]\n",
      " [1904.            5.22208333]\n",
      " [1908.            4.15467867]\n",
      " [1912.            3.90331675]\n",
      " [1920.            3.56951267]\n",
      " [1924.            3.82454477]\n",
      " [1928.            3.62483707]\n",
      " [1932.            3.59284275]\n",
      " [1936.            3.53880792]\n",
      " [1948.            3.67010309]\n",
      " [1952.            3.39029111]\n",
      " [1956.            3.43642612]\n",
      " [1960.            3.20583007]\n",
      " [1964.            3.13275665]\n",
      " [1968.            3.32819844]\n",
      " [1972.            3.13583758]\n",
      " [1976.            3.0789588 ]\n",
      " [1980.            3.10581822]\n",
      " [1984.            3.06552909]\n",
      " [1988.            3.09357349]\n",
      " [1992.            3.16111704]\n",
      " [1996.            3.14255244]\n",
      " [2000.            3.08527867]\n",
      " [2004.            3.10265829]\n",
      " [2008.            2.99877553]\n",
      " [2012.            3.03392977]]\n"
     ]
    }
   ],
   "source": [
    "# CSV file with variables YEAR,TIME\n",
    "csv = \"\"\"1896,4.47083333333333\n",
    "1900,4.46472925981123\n",
    "1904,5.22208333333333\n",
    "1908,4.1546786744085\n",
    "1912,3.90331674958541\n",
    "1920,3.5695126705653\n",
    "1924,3.8245447722874\n",
    "1928,3.62483706600308\n",
    "1932,3.59284275388079\n",
    "1936,3.53880791562981\n",
    "1948,3.6701030927835\n",
    "1952,3.39029110874116\n",
    "1956,3.43642611683849\n",
    "1960,3.2058300746534\n",
    "1964,3.13275664573212\n",
    "1968,3.32819844373346\n",
    "1972,3.13583757949204\n",
    "1976,3.07895880238575\n",
    "1980,3.10581822490816\n",
    "1984,3.06552909112454\n",
    "1988,3.09357348817\n",
    "1992,3.16111703598373\n",
    "1996,3.14255243512264\n",
    "2000,3.08527866650867\n",
    "2004,3.1026582928467\n",
    "2008,2.99877552632618\n",
    "2012,3.03392977050993\"\"\"\n",
    "\n",
    "# Read into a numpy array (as floats)\n",
    "olympics = np.genfromtxt(io.BytesIO(csv.encode()), delimiter=\",\")\n",
    "print(olympics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll take the race time as the *target variable* $y$ and the year of the race as the only non-trivial *feature* $x = x_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = olympics[:, 0:1]\n",
    "y = olympics[:, 1:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting $y$ vs $x$, we see that a linear model could be a decent fit for this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGoxJREFUeJzt3X+4XVV95/H3h4SQEFEiCRYJNIo6/hqscKOonencYCmiDY4weWxVoODDmEcL1moGpt6MJEiH+/gUH8YhiNJOLFZNaWmBqojmFMcfwL0xEBKhEiyMETSJgPwov/nOH3vdnXNPzj13n5u7zz7n3M/rec5zzll73b3XujvZ37vW2mttRQRmZmYA+1VdADMz6x4OCmZmlnNQMDOznIOCmZnlHBTMzCznoGBmZjkHBTMzyzkomJlZzkHBzMxys6suQLsWLlwYS5YsqboYZmY9ZdOmTbsjYtFk+XouKCxZsoTR0dGqi2Fm1lMk3Vckn7uPzMws56BgZmY5BwUzM8s5KJiZWc5BwczMcg4KUzE8DLXa+LRaLUs3M+thDgpTsXQprFixJzDUatn3pUurLZeZ2T7quXkKXWFwEDZsyALBypWwbl32fXCw6pKZme0TtxSmanAwCwhr12bvDghm1gccFKaqVstaCEND2XvjGIOZWQ9yUJiKsTGEDRtgzZo9XUkODGbW4xwUpmJkZPwYwtgYw8hIteUyM9tHioiqy9CWgYGB8IJ4ZmbtkbQpIgYmy+eWgpmZ5RwUzMws56BgZmY5BwUzM8s5KJiZWc5BwczMcg4KZmaWc1AwM7Ocg4KZmeUcFMzMLOegYGZmuVKDgqR7Jd0h6TZJey1YpMylkrZL2iLpmDLLY2ZmrXXiyWuDEbF7gm3vAF6ZXm8G1qV3MzOrQNXdRycDX4rMzcDBkg6ruExmZjNW2UEhgG9J2iTp7CbbDwd+Vvd9R0ozM7MKlN199LaIuF/SocCNku6KiO/WbVeTn9nrAQ8poJwNcOSRR5ZTUjMzK7elEBH3p/edwDXAmxqy7ACOqPu+GLi/yX6uiIiBiBhYtGhRWcU1M5vxSgsKkuZLOmjsM3ACsLUh27XAaekupOOAX0fEA2WVyczMWiuz++glwDWSxo7zNxHxTUkfAoiIy4GvAycB24F/A/6oxPKYmdkkSgsKEfFT4A1N0i+v+xzAh8sqg5mZtafqW1LNzKyLOCiYmVnOQcHMzHIOCmZmlnNQMDOznIOCmZnlHBTMzCznoGBmZjkHBTMzyzkomJlZzkHBzMxyDgpmZpZzUDAzs5yDgpmZ5RwUzMws56BgZmY5BwUzM8s5KJiZWc5BwczMcg4KZmaWc1AwM7Ocg4KZmeUcFMzMLOegYGZmOQcFMzPLOSiYmVnOQcHMzHIOCmZmlnNQMDOznIOCmZnlHBTMzCznoGBmZjkHBTMzyzkomJlZzkHBzMxyDgpmZpYrPShImiVps6Trm2w7Q9IuSbel1wfLLo+ZmU1sdgeOcS5wJ/DCCbZ/LSI+0oFymJnZJEptKUhaDLwT+GKZxzEzs+lRdvfRZ4FVwPMt8pwiaYukqyUd0SyDpLMljUoa3bVrV3slGB6GWm18Wq2WpZuZ2TilBQVJ7wJ2RsSmFtmuA5ZExNHAt4H1zTJFxBURMRARA4sWLWqvIEuXwooVewJDrZZ9X7p077wOIGY2w5XZUngbsFzSvcBXgWWSrqrPEBG/ioin0tcvAMdOeykGB2HDhiwQrF6dvW/YkKU3aieAmJn1oUJBQdICSa+T9HJJhX4mIs6PiMURsQR4L7AxIt7fsN/D6r4uJxuQnn6Dg7ByJaxdm703Cwhj+YoGEDOzPjThBV7SiyT9d0l3ADcDnwc2APdJ+ltJU7pSSlojaXn6eo6kbZJuB84BzpjKPidVq8G6dTA0lL03dhHVKxpAzMz6kCKi+QbpRuBLwHUR8XDDtmOBDwB3RMSVpZeyzsDAQIyOjhb/gbEuoLG/+Bu/T5R/5cosgLilYGZ9QNKmiBiYLN+E8xQi4ndbbNsEtBpA7h4jI+Mv7GNdRCMje1/sGwPG4KC7kMxsRpmwpZBnkAS8D3h5RKyRdCTwGxFxaycK2KjtlkI7hoezQeX6AFCrZQFk1apyjmlm1gFFWwpFgsI6snkGyyLiNZIWAN+KiEpuySk1KJiZ9al97j6q8+aIOEbSZoCIeEjSnH0uoZmZdZ0it5c+I2kWEACSFtF6hrKZmfWoIkHhUuAa4FBJnwa+B1xUaqnMzKwSk3YfRcSXJW0CjgcEvDsiyplkZmZmlSq6dPYvgf+b8s+TdExE/Ki8YpmZWRUmDQqS1pLNNL6HNK6Q3peVVywzM6tCkZbCCuCoiHi67MKYmVm1igw0bwUOLrsgZmZWvSIthT8HNkvaCowtc01ELJ/4R8zMrBcVCQrrgYuBO/D8BDOzvlYkKOyOiEtLL4mZmVWuSFDYJOnPgWsZ333kW1LNzPpMkaDwxvR+XF2ab0k1M+tDRWY0+0ECZmYzxIRBQdL7I+IqSR9rtj0i/qK8YpmZWRVatRTmp/eDmmxr/RAGMzPrSa0ex/n59PHbEfH9+m2S3lZqqczMrBJFZjT/r4JpZmbW41qNKbwFeCuwqGFc4YXArLILZmZmndeqpTAHeAFZ4Dio7vUIcGr5ResTw8NQq41Pq9WydDOzLtNqTOEm4CZJ/yci7utgmfrL0qWwYgVs2ACDg1lAGPtuZtZlisxTcEDYF4ODWQBYsQJWroR16/YECDOzLlNkoNn21eBgFhDWrs3eHRDMrEs5KHRCrZa1EIaGsvfGMQYzsy4xaVCQ9CpJ30nPU0DS0ZI+WX7R+kT9GMKaNXu6khwYzKwLFWkpfAE4H3gGICK2AO8ts1B9ZWRk/BjC2BjDyEi15TIza6LIKqkHRsStkurTni2pPP1n1aq90wYHPa5gZl2pSEtht6SjSOsdSToVeKDUUpmZWSWKtBQ+DFwBvFrSz4F/Bd5faqnMzKwSReYp/BR4u6T5wH4R8Wj5xTIzsyoUufvoIkkHR8TjEfGopAWSLuxE4czMrLOKjCm8IyIeHvsSEQ8BJ5VXJDMzq0qRoDBL0gFjXyTNAw5okX8cSbMkbZZ0fZNtB0j6mqTtkm6RtKTofs3MbPoVCQpXAd+RdJakM4EbgfVtHONc4M4Jtp0FPBQRrwAuAS5uY79mZjbNJg0KETEMfBp4DfA6YG1Km5SkxcA7gS9OkOVk9gSYq4Hj1TAhwszMOqfILalExDeAb0xh/58FVtH8Oc8AhwM/S8d4VtKvgUOA3VM4lpmZ7aMidx8dJ2lE0mOSnpb0nKRHCvzcu4CdEbGpVbYmadFkX2dLGpU0umvXrskObWZmU1RkTOFzwB8AdwPzgA9S7BnNbwOWS7oX+CqwTNJVDXl2AEcASJoNvAh4sHFHEXFFRAxExMCiRYsKHNrMzKai0NLZEbEdmBURz0XEXwGTLtwTEedHxOKIWEK2gN7GiGicCX0tcHr6fGrKs1dLwczMOqPImMK/SZoD3CZpmGzdo/lTPaCkNcBoRFwLXAn8taTtZC0Er75qZlYhTfaHuaTfBH4JzAH+hKyL57LUeui4gYGBGB0dreLQZmY9S9KmiBiYLF+RW1Lvi4gnI+KRiLggIj4GvGRaSml7DA/v/eCdWi1LNzPrkAmDQpqJ/AeSPi7p9SntXZJ+QDb4bNNp6dLxT2Qbe2Lb0qXVlsvMZpRWYwpXkt0ZdCtwqaT7gLcA50XEP3SicDPK2BPZVqyAlSuzZznXP7HNzKwDWgWFAeDoiHhe0lyyCWWviIhfdKZoM9DgYBYQ1q6FoSEHBDPruFZjCk9HxPMAEfEk8BMHhJLValkLYWgoe28cYzAzK1mrlsKrJW1JnwUclb4LiIg4uvTSzSRjYwhjXUaDg+O/m5l1QKug8JqOlcJgZGR8ABgbYxgZcVAws46ZcJ6CJE02u7hInunmeQolGR7O7nSqD0C1WhaUVq2qrlxmNi2mY55CTdIfSzqyYcdzJC2TtJ49S1RYr/MtsWZG6+6jE4Ezga9IehnwMDAXmAV8C7gkIm4rv4jWEb4l1sxoERTSHUeXAZdJ2h9YCDxR/7xm6zO+JdZsxiu6SuozEfGAA0Kf8y2xZjNeoaBgM0D9LbFr1uzpSnJgMJtRHBQs0+qWWDObMYosnf0R4MsR8VBnitSab0k1M2vftC2dDfwGMCJpg6QTJTV7rrKZmfWBIs9T+CTwSrJVU88A7pZ0kaSjSi6bmZl1WNG7jwL4RXo9CywArk6P5zQzsz4x6TOaJZ1DNnN5N/BF4BMR8Yyk/YC7Aa+BYGbWJyYNCmST1t4TEffVJ6bnLLyrnGKZmVkVJg0KEbG6xbY7p7c4ZmZWJc9TMBse3nuSXq2WpZvNMA4KZl4h1ixXZEzBrL95hViznFsKZjB+hdiVKx0QbMZyULD29Gv/u1eINQMcFKxd/dj/7hVizXIOCtae+v731av3XEx7ubvFK8Sa5SZdJbXbeJXULrF69Z4ntK1Z0zzP8HDWgqgPGLVadrFd5YnwZp00naukmo1XtP+9H7uazPqcg0I/K2NQuJ3+937sajLrcw4K/ayMv9Tb7X/3rZ5mPcVjCv1uLBBUNSmr6uObGeAxBRtT5V/qRbua+nXug1kPclDod1VOyira1eQBabOu4e6jflb/l/rg4N7fu4m7mcxKVXn3kaS5km6VdLukbZIuaJLnDEm7JN2WXh8sqzwzUi9NyvKAtFlXKK2lIEnA/Ih4TNL+wPeAcyPi5ro8ZwADEfGRovt1S6FPuaVgVqrKWwqReSx93T+9equvyjrDaw+ZdY1SB5olzZJ0G7ATuDEibmmS7RRJWyRdLemICfZztqRRSaO7du0qs8hWhV7q5jLrcx0ZaJZ0MHAN8McRsbUu/RDgsYh4StKHgBURsazVvtx9ZGbWvsq7j+pFxMPAPwMnNqT/KiKeSl+/ABzbifKYmVlzZd59tCi1EJA0D3g7cFdDnsPqvi4H7iyrPGYd5Ql51qPKbCkcBtQkbQFGyMYUrpe0RtLylOecdLvq7cA5wBkllsesczwhz3qUJ6+ZlcW32VoX6aoxBbMZyRPyrAc5KJgV1e44QZXrTplNkYOCWVHtjBN4Qp71KAcFs6LaeZKcJ+RZj/JAs1m7Vq/OxgmGhrJWgFkP8ECzWRk8TmB9zkHBrCiPE9gM4KBgVlSvjBN4NrXtAwcFs6JWrdp7UHlwMEvvJp5NbftgdtUFMLNpVn+XlGdTW5vcUjDrR55NbVPkoGDWj3yXlE2Rg4JZv/FdUrYPHBTM+k2v3CVlXckzms3MZgDPaDYzs7Y5KJiZWc5BwczMcg4KZlXykhTWZRwUzKrkJSmsy3iZC7MqeUkK6zJuKZhVzUtSTB93x+0zBwWzqnlJiunj7rh95qBgViUvSTG92nmOtlsVTTkomFWpnSUpyriI9eOFsWh3nFsVzUVET72OPfbYMJuRNm6MWLgwe2/2vVv2WbWxOgwNTV6XdvL2OGA0ClxjK7/It/tyULAZrYyL2HTv8+KL997Hxo1ZetmmEuSGhrJL4dBQ+eWrUNGg4O4js15Sxp1K073PMrplinZztbtC7HQP8vdDd1yRyNFNL7cUbEbrhZZCGfussuusnZZPF3fH4e4jsz7TKxfGMdPdLVNVN1e7v/cuHadwUDDrN2X01Vd5YeyGQFNUuxf6LhyncFAws+lV9MJYNID02l/gRS/0VZdzAg4KZjb9ilwYp9IHP12Bpiy9Us4WHBTMbHqV9RfwdAea6dbOhb7KLr5JOCiY2fQp6y/gLu1qGafKgDR2rGn43VceFIC5wK3A7cA24IImeQ4AvgZsB24Blky2XwcFswqUcWHs4q6WrjMNwbNoUChz8tpTwLKIeAPwW8CJko5ryHMW8FBEvAK4BLi4xPKY2VStWrX3pLbBwSx9qtqdaNZP2p3k1sHl1UsLCik4PZa+7p9e0ZDtZGB9+nw1cLwklVUmM+siZQSaXtHurO8OLq9e6jIXkmZJug3YCdwYEbc0ZDkc+BlARDwL/Bo4pMwymZlVrp0lvju8vHqpQSEinouI3wIWA2+S9PqGLM1aBY2tCSSdLWlU0uiuXbvKKKqZWWcV7RLqcDebsvGH8kn6H8DjEfGZurQbgE9FxA8lzQZ+ASyKFoUaGBiI0dHR8gtsZlamsRZAh57NLWlTRAxMlq+0loKkRZIOTp/nAW8H7mrIdi1wevp8KrCxVUAwM+sLXfzEvTK7jw4DapK2ACNkYwrXS1ojaXnKcyVwiKTtwMeA80osj5lZd+jiO6861n00Xdx9ZGbWvsq7j8zMrPc4KJiZWc5BwczMcg4KZmaWc1AwM7Ncz919JGkXcF8HDrUQ2N2B43RKv9UH+q9O/VYf6L869XJ9fjMiFk2WqeeCQqdIGi1y+1av6Lf6QP/Vqd/qA/1Xp36rTzPuPjIzs5yDgpmZ5RwUJnZF1QWYZv1WH+i/OvVbfaD/6tRv9dmLxxTMzCznloKZmeVmTFCQ9JeSdkraWpf2Bkk/lHSHpOskvbBu2/mStkv6F0m/V5d+YkrbLqnSVV3bqZOk35W0KaVvkrSs7meOTenbJV1a1SNR2z1HafuRkh6T9PG6tJ48R2nb0WnbtrR9bkrvuXMkaX9J61P6nZLOr/uZrjhHko6QVEvl2ybp3JT+Ykk3Sro7vS9I6Uq//+2Stkg6pm5fp6f8d0s6faJjdr2ImBEv4D8CxwBb69JGgN9Jn88E1qbPrwVuBw4AXgbcA8xKr3uAlwNzUp7X9kid3gi8NH1+PfDzup+5FXgL2ZPwvgG8o9vrU7f974C/BT6evvfyOZoNbAHekL4fAszq1XME/CHw1fT5QOBeYEk3nSOyJf6PSZ8PAn6S/v8PA+el9POAi9Pnk9LvX8BxwC0p/cXAT9P7gvR5QVX/7vblNWNaChHxXeDBhuR/B3w3fb4ROCV9PpnsH/NTEfGvwHbgTem1PSJ+GhFPA19NeSvRTp0iYnNE3J/StwFzJR0g6TDghRHxw8j+dX8JeHf5pd9bm+cISe8m+8+3rS5/z54j4ARgS0Tcnn72VxHxXA+fowDmK3uq4jzgaeARuugcRcQDEfGj9PlR4E6yZ8efDKxP2daz5/d9MvClyNwMHJzOz++RPTPmwYh4iOz3cGIHqzJtZkxQmMBWYOyBP/8FOCJ9Phz4WV2+HSltovRuMlGd6p0CbI6Ip8jKv6NuW7fVqWl9JM0H/htwQUP+Xj5HrwJC0g2SfiRpVUrvyXMEXA08DjwA/D/gMxHxIF16jiQtIWtR3wK8JCIegCxwAIembL18bShkpgeFM4EPS9pE1nR8OqU366+NFundZKI6ASDpdcDFwH8dS2qyj26q00T1uQC4JCIea8jf7fWBies0G/ht4H3p/T9LOp7ur9NE9XkT8BzwUrJu2D+V9HK6sD6SXkDWFfnRiHikVdYmab1ybShkdtUFqFJE3EXWZEfSq4B3pk07GP8X9mJgrOtlovSu0KJOSFoMXAOcFhH3pOQdZPUY01V1alGfNwOnShoGDgael/QksInePUc7gJsiYnfa9nWy/vur6M1z9IfANyPiGWCnpO8DA2R/UXfNOZK0P1lA+HJE/H1K/qWkwyLigdQ9tDOlT3Rt2AH8p4b0fy6z3GWZ0S0FSYem9/2ATwKXp03XAu9Nfe4vA15JNtA3ArxS0sskzQHem/J2jYnqJOlg4J+A8yPi+2P5U9P4UUnHpTtaTgP+seMFn8BE9YmI/xARSyJiCfBZ4KKI+Bw9fI6AG4CjJR2Y+uF/B/hxr54jsi6jZemOnflkA7N30UXnKP0+rwTujIi/qNt0LTB2B9Hp7Pl9Xwuclup0HPDrdH5uAE6QtCDdqXRCSus9VY90d+oFfIWsb/MZsqh+FnAu2d0GPwH+J2kyX8r/Z2R3SPwLdXd6kN198JO07c96pU5k/1kfB26rex2atg2Q9QvfA3yu/vfQrfVp+LlPke4+6uVzlPK/n2zgfCswXJfec+cIeAHZnWHbgB8Dn+i2c0TWTRdkd32N/b84iezOr+8Ad6f3F6f8Av53KvcdwEDdvs4kuyllO/BHVf6725eXZzSbmVluRncfmZnZeA4KZmaWc1AwM7Ocg4KZmeUcFMzMLOegYGZmOQcF61uS5km6SdJ8SXdJ+vd121ZJurzVz09zWV4r6XZJm9MaO2Ufb26q+6yyj2X9xUHB+tmZwN9HxOPAR4HL0kzUw8nWfjq/5U8X0MZF9z3A1RHxxoi4dxr211JEPEm2cump07E/mzkcFKznSFqaHnAyN7UCtkl6fZOs7yMtTxAR3ySbiXsacAnwqciWOEbSeZJuTftcXXec65Q9kGibpA+mtNmSHpZ0oaRbyRZ9qy/bMZJuSfv6O0kvkrQc+AjwIUnfbsi/1/4kXSBpRNJWSZenpRiQ9CpJG1OL40djLY6Jyg/8Q/odmBVX9ZRqv/yaygu4EPgM2ZID5zfZPgf4RUPaS8mWZqjVpZ0EXEa2fMF+wDeBt6ZtY0sbHEi2TMMCskUkA3jPBOX6MfDb6fNFZMtFj5X3o03y77U/xi+p8BXSMitki/39fvo8N5WrVflnAzurPld+9dZrRq+Saj1tDdnCak8C5zTZvhB4uD4hIu6XtBG4vi75BOAdwOb0/QVkzzX4AfAn6a98yFa9PIpsbZynyVabHUfSIcDciPheSloP/HWBujTu73hJnyC78C8ENkm6GVgYEdelujyZjjlh+SPiWUkhaV5EPFGgHGYOCtazXkx2Adyf7OL5eMP2J1J6o+fTa4yACyPiyvpMkt5O9ujJ4yLiCUnfq9vfExHRbNGwqT43Od+fpAPJFrw7JiJ+LunCuuNOdMy9yl9nDvDUFMtlM5DHFKxXXQEMAV8me2jQOJGNF8xSevB9CzcAZ6WlnZG0WNJC4EXAgykgvA5YOlmBInsOwhOS3pqSPgDcVLRCyTyyoLVb0kHseZzqQynt91M556YAMlH5kfQSsmdxP9/kOGZNuaVgPUfSacCzEfE36W6dH0haFhEbG7J+i2xp5G/vtZMkIr4u6dXAzWk891Gyh8P8E3C2pNvJngFwS8HifQBYJ2keaQnlNqpGRPxK0nqyZbLvazju+4DPS/o0WZfTKS3KvxsYTPUwK8xLZ1vfkvRG4GMR8YGqy1IFSf8I/GlEbK+6LNY73H1kfSsiNgO1mTiBS9IBZPMiHBCsLW4pmJlZzi0FMzPLOSiYmVnOQcHMzHIOCmZmlnNQMDOz3P8HqiPQn9yUwZwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, y, 'rx')\n",
    "plt.ylabel(\"y (Race time)\")\n",
    "plt.xlabel(\"x (Year of race)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Iterative solution (coordinate descent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expanding out the sum of square residuals for this simple case (where $\\mathbf{w}=[w_0, w_1]$) we have:\n",
    "$$SSR(w_0, w_1) = \\sum_{i=1}^{n}(y_i - w_0 - w_1 x_i)^2$$\n",
    "Let's start with an initial guess for the slope $w_1$ (which is clearly negative from the plot)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1 = -0.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then using the maximum likelihood update, we get the following estimate for the intercept $w_0$:\n",
    "$$w_0 = \\frac{\\sum_{i=1}^{n}(y_i-w_1 x_i)}{n}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "786.0197711453593\n"
     ]
    }
   ],
   "source": [
    "def update_w0(x, y, w1):\n",
    "    return (y - w1*x).sum()/len(x) #can also use x.size\n",
    "\n",
    "w0 = update_w0(x, y, w1)\n",
    "print(w0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can update $w_1$ based on this new estimate of $w_0$:\n",
    "$$w_1 = \\frac{\\sum_{i=1}^{n} (y_i - w_0) \\times x_i}{\\sum_{i=1}^{n} x_i^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.3998724072997095\n"
     ]
    }
   ],
   "source": [
    "def update_w1(x, y, w0):\n",
    "    return ((y - w0)*x).sum()/np.square(x).sum()\n",
    "\n",
    "w1 = update_w1(x, y, w0)\n",
    "print(w1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's examine the quality of fit for these values for the weights $w_0$ and $w_1$. We create a vector of \"test\" values `x_test` and a function to compute the predictions according to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_test = np.arange(1890, 2020)[:, None]\n",
    "\n",
    "def predict(x_test, w0, w1): \n",
    "    return w0 + x_test*w1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the test predictions with a blue line on the same plot as the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYFFXWx/HvYQDBhCQjrqirLgwgYSSYwQBiwrhmVFwT5gSYWUxgdo0oKIZVMSCYRcCIAoOkAVTQVwVFQTEHFDjvH7dYWpzp6YHuru6e3+d56unu29VVZwqmz9y6VeeauyMiIlKRGnEHICIiuU2JQkREklKiEBGRpJQoREQkKSUKERFJSolCRESSUqIQEZGklChERCQpJQoREUmqZtwBpEOjRo28adOmcYchIpJXJk+e/LW7N65svYJIFE2bNqW0tDTuMERE8oqZfZrKejr1JCIiSSlRiIhIUkoUIiKSlBKFiIgkpUQhIiJJKVGIiEhSsSUKM6tjZhPNbJqZzTSz/lH7lmY2wczmmNnjZlY7rhhFRCTeHsUSoIu7bw+0BrqZWUdgIHCzu28DfAv0ylQAixbBuefC999nag8iIvkvtkThwU/Ry1rR4kAX4MmofRjQI1MxjBkDt90GzZvDs89mai8iIvkt1jEKMysys6nAQmA08BHwnbsvjVaZD2xWwWdPNrNSMytdtGjRau3/iCPg3XehYUM44AA46qjQyxARkZViTRTuvszdWwNNgPZAs/JWq+Czg929xN1LGjeutFRJhXbYAUpLoX9/ePJJaNYM/vtf8HL3KiJS/eTEVU/u/h3wGtAR2MDMVtSgagJ8ken9164Nl18OU6bA3/8ORx8N++8P8+Zles8iIrkvzqueGpvZBtHzusCewGxgHHBotFpPYGS2YiouhrffhptugrFjw+t77oHly7MVgYhI7omzR7EJMM7MpgOTgNHu/hzQBzjPzOYCDYEh2QyqqChcCVVWFk5LnXoqdOkCc+ZkMwoRkdxhXgAn40tKSjwTZcbdYehQOP98WLIEBgyAc86BmgVRnF1Eqjszm+zuJZWtlxNjFLnKDHr1glmzoGtXuPBC6NQJpk+POzIRkexRokjBppvCiBHw+OPw6afQrl0Y/F6yJO7IREQyT4kiRWZw+OEwe3a4/2LAAGjbNtyHISJSyJQoqqhhQ3joIXj+efjhB9hxxzD4/fPPcUcmIpIZShSrqXt3mDkzXBV1yy3QsmUoCSIiUmiUKNbA+uvDnXfC66+HK6H23BNOOgm++y7uyERE0keJIg123RWmTYM+feCBB0KRwWeeiTsqEZH0UKJIk7p14brrYMIE2HBDOOigMPj91VdxRyYismaUKNKsXTuYNAmuugpGjgy9i4ceUpFBEclfShQZUKsWXHIJTJ0K220Hxx0H++4Ln30Wd2QiIlWnRJFBzZrBm2/CrbeGAe/i4jD4rSKDIpJPlCgyrKgIzjorFBns1Al694bdd4cPP4w7MhGR1ChRZMmWW8LLL8P998OMGdCqFQwcCEuXVv5ZEZE4KVFkkRkcf3woMti9O/TtCx06hLEMEZFcpUQRg002gaefDlOvfv45lJSEwe/ffos7MhGRv1KiiNEhh4TexTHHwDXXQJs2MH583FGJiPyZEkXMGjQId3O/9BL88gvsvHMY/P7pp7gjExEJlChyRNeu4cqo3r3h9tuhRQt45ZW4oxIRUaLIKeutB//5D7zxBtSpE5LHCSfA4sVxRyYi1ZkSRQ7aeedwJVS/fqH8R/Pm8NRTcUclItWVEkWOqlMnDHBPmhSukjr00LB8+WXckYlIdaNEkePatIGJE0PSeO650Lt44AEVGRSR7FGiyAO1aoXTUFOnhkRxwgnQrRt88knckYlIdaBEkUf+8Y8w0H377eF+ixYtwuC3igyKSCYpUeSZGjXCJbRlZSvvudh1V3j//bgjE5FCpUSRp7bYAl58EYYNC3d3b799GMf444+4IxORQqNEkcfMwqRIs2fDAQeEelHt28N778UdmYgUEiWKArDRRvDEE6HQ4JdfhmTRrx/8+mvckYlIIVCiKCAHHRROQ/XsCdddB61bw1tvxR2ViOS72BKFmW1uZuPMbLaZzTSzs6P2BmY22szmRI/144oxH9WvD0OGwOjR8PvvsMsucMYZ8OOPcUcmIvkqzh7FUuB8d28GdAR6m1lzoC8wxt23AcZEr6WK9twzzKR39tlhnu7i4jD4LSJSVbElCndf4O7vRc9/BGYDmwEHAsOi1YYBPeKJMP+tuy7ccgu8/XZ43r17GPz+5pu4IxORfJITYxRm1hRoA0wANnL3BRCSCbBhfJEVhk6dYMoUuPRSePTRcHf3E0+oDIiIpCb2RGFm6wJPAee4+w9V+NzJZlZqZqWLFi3KXIAFYq21YMAAKC2FzTeHww+Hgw+GBQvijkxEcl2sicLMahGSxCPu/nTU/JWZbRK9vwmwsLzPuvtgdy9x95LGjRtnJ+ACsP328O67MGhQmFWvWTMYOlS9CxGpWJxXPRkwBJjt7jclvDUK6Bk97wmMzHZsha5mTbjwQpg2LSSOXr1g773h44/jjkxEclGcPYqdgGOBLmY2NVq6A9cBe5nZHGCv6LVkwLbbwrhxcNddMGECtGwZBr+XLYs7MhHJJeYFcM6hpKTES0tL4w4jr82bB6ecEi6h7dgx3IvRvHncUYlIJpnZZHcvqWy92AezJTdsvjk8/zw8/DDMmRMmTBowINy0JyLVmxKF/I8ZHH10KANy8MFw+eWwww7hSikRqb6UKOQvNtww3G8xciR8/TV06AAXXaQigyLVlRKFVOiAA2DmzHBV1PXXQ6tW8PrrcUclItmmRCFJbbABDB4MY8aEKVd33x1OOw1+SPnWSBHJd0oUkpIuXWD6dDjvvJA4iovD4LeIFD4lCknZOuvAjTfC+PFQrx7stx8cc0wYxxCRwqVEIVXWoUOYbvWKK2D48FAG5LHHVAZEpFApUchqqV0brrwSJk+GLbeEI4+EHj3g88/jjkxE0k2JQtZIy5bwzjtwww1hVr3mzeHee9W7ECkkShSyxoqK4Pzzw2B327Zw8smwxx7w0UdxRyYi6aBEIWnz97+Hy2jvuSeckmrZEm66SUUGRfKdEoWkVY0aoUcxc2boVZx/Puy4I5SVxR2ZiKwuJQrJiCZNYNSoUArk44/DKan+/VVkUCQfKVFIxpjBEUfA7Nlw2GHhKql27WDixLgjE5GqUKKQjGvUCB55BJ59Fr79Fjp1Cqekfvkl7shEJBVKFJI1++0Xxi7+9a8wyN2yZZhhT0RymxKFZFW9enD33SFB1KgRakidfDJ8/33ckYlIRZQoJBa77w7TpsGFF66cdvXZZ+OOSkTKo0QhsVl7bRg0CCZMgIYNw/wXRx4JixbFHZmIJFKikNiVlITpVv/9b3jqqVBk8JFHVAZEJFcoUUhOqF0bLrsMpkwJd3gfcwzsvz/Mmxd3ZCKSUqIws/pmVmxmW5mZkotkTHExvP023HxzGPAuLg6D38uXxx2ZSPVV4Ze+mdUzs4vNbAbwLnAPMBz41MyeMLPO2QpSqpeiIjjnHJgxA9q3D1OvdukCc+bEHZlI9ZSsd/AkMA/Yxd23c/ed3b3E3TcHrgMONLNeWYlSqqWttgqly4cMgalToVUruP56WLo07shEqhfzAhgxLCkp8dLS0rjDkAz64gs4/XQYOTIMfg8ZEhKHiKw+M5vs7iWVrVfpeIMFx5jZ5dHrv5lZ+3QEKZKqTTeFESPC1KuffRZqRl1+OSxZEndkIoUvlYHpO4FOwJHR6x+BOzIWkUgFzEJxwVmzwv0WAwZAmzZhhj0RyZxUEkUHd+8N/Abg7t8CtTMalUgSDRvCgw/CCy/ATz/BTjuFwe+ff447MpHClEqi+MPMigAHMLPGQFouVjSzoWa20MzKEtoamNloM5sTPdZPx76k8OyzTygyePrpcOut0KIFvPpq3FGJFJ5UEsVtwAhgQzO7GngLuCZN+38A6LZKW19gjLtvA4yJXouUa7314Pbb4Y03oFYt2Gsv6NULvvsu7shECkelicLdHwEuAq4FFgA93P2JdOzc3d8AFq/SfCAwLHo+DOiRjn1JYdtll1BksG9fGDYsFBl85pm4oxIpDKneZf0V8CYwHqhrZm0zFxIbufsCgOhxwwzuSwpI3bpw7bWhyOCGG8JBB8Hhh8NXX8UdmUh+S+Xy2AHAdMIpqBuj5YYMx1UpMzvZzErNrHSRyo1KgnbtYNIkuPrqcN9Fs2Zh8LsAbhkSiUUqPYrDga3dfXd37xwtXTIY01dmtglA9LiwvJXcfXB0p3hJ48aNMxiO5KNateDii8Md3c2aQc+e0L17uAdDRKomlURRBmyQ6UASjAJ6Rs97AiOzuG8pMM2awZtvwm23hcfiYrjjDhUZFKmKVBLFtcAUM3vZzEatWNKxczN7FHgH2M7M5ke1o64D9jKzOcBe0WuR1VajBpx5JpSVQadOcMYZsNtu8MEHcUcmkh8qrfVkZjMJlWNnkHD/hLu/ntnQUqdaT5Iq93BV1Lnnwq+/wpVXwvnnh1NVItVN2mo9AV+7+23uPs7dX1+xpCFGkawzg+OPh9mzYd99oV8/6NAhTJgkIuVLJVFMNrNrzayTmbVdsWQ8MpEM2njjMO3qk0+GyrQ77ACXXAK//RZ3ZCK5J5VTT+PKafYMX/lUJTr1JGti8eJw+umBB2C77UIJ8512ijsqkcxL26mnhEtiO2fp8liRrGrQAO6/H15+OfQodtkFzjorFBwUEahZ0Rtmdoy7P2xm55X3vrvflLmwRLJv773DlVEXXxzqR40aBYMHh3aR6ixZj2Kd6HG9cpZ1MxyXSCzWXXflPRd16kDXrnDCCeH0lEh1VWGPwt3viZ6+6u5vJ75nZjqDKwVtp53CXd0DBsDAgfDii+FGvUMOiTsykexL5aqn/6TYJlJQ6tQJ9aJKS8NUrIceGhLFggVxRyaSXcnGKDoBOwKNVxmnWB8oynRgIrmidWuYOBFuvBGuuALGjoWbbw71o8zijk4k85L1KGoTxiJq8ufxiR+AQzMfmkjuqFkT+vQJc160aBHGLbp2hU8+iTsykcxL5T6KLdz90yzFs1ry8j6KQYPCXV6dO69sGzcu1Me+6KL44pJKLV8Od98dEod7mAOjd+9QU0okn6TzPoqcThJ5a4cdwqw646L7GceNC6932CHeuKRSNWqEebrLylbec7HLLqEsiEgh0t9AcencGYYPD8nh8svD4/Dhf+5hSE7bYgt44YUwKdL774exjKuvhj/+iDsykfRSoqiqQYNW9gJWGDcutFdlHQhJ4bTTwjWYp52mJJGHzODYY2HWLOjRAy69NHQK33sv7shE0ieVqVC3NbMxZlYWvW5lZpdmPrQclcopo1RPK40bB3fdBZddFh5XTS6QetKRWG20ETz+OIwYEebobt8e+vYNpcxF8p67J12A14H2wJSEtrLKPpfNpV27dp5VY8e6N2rkftll4XHs2Kqvs+L9Fe2rvq7KegMHlv+5gQPX7OeU1bJ4sXuvXu7gvu227m+8EXdEIuUDSj2F79hUEsWk6DExUUxNZePZWrKeKNxDAoDwuDrrVOXLPV1JJ51SiT/VnzEXE10aYho92r1p0/Bf4PTT3b//Ps0xiqyhdCaKF4Gtgfei14cCL6ay8WwtedmjqKrKElO69pfqF2QqySmfe02pxl5JXD/95D5q54HembG++ebuL7yQhdhzMfFKTkpnotgKeBX4BfgceAtomsrGs7VkNVGk8wuyqvusLAmk0stJdV+pxJ7OhJmLvaaq/HyV/H/4fYNGfvwWYx3cr9lrrC9rWPWkk/I6aUpyKa9THbZVoMk3bYnifyuGarLrpbp+NpesJops/6eq6l/l6ejBVGVba3oKrirrpbuXlopUYk8xoSxv1MjH7XKZL6SR96g31h9/3H358nK2k85eWhqSXNp6hfm8rTiSbxaks0dxDbBBwuv6wFWpbDxbSyxjFNmSzr8gqyKNX5Bp6VFUJa50yVDC/PLky7xdu7B6jx7un39exX2m81hl+98wn7eV7eSbBelMFFPKaXsvlY1naynoRJGKdP+Fkot/jaYaV7pU5Zd5Nb5klo4e64MGudep416vnvt99yX0LtLRS8vVXmG+byvbyTfD0pkopgNrJbyuC8xMZePZWqp9okinOLrYcfWaksnSoP6HH7rvumv4TdxjD/fPH07Dl0yGk1xe9QLSua24km8GpTNRXBQNYPcCToyeX5TKxrO1KFGkUQ6dP/2TfI6rknWWLXO/6y737nXH+kIa+ROnj/WlS331e2m5euVaPm8rruSbYWkdzAb2AW4AbgS6pvKZbC5KFFIIvu030Pt2CFdGdejgXlbmq99LS0UuXl2Uq9uKI/lmQaqJotIy4/kgL8uMi5TDHR59NFSk/eGHUN2lTx+oXTvuyCQlqUwfkENTDKRaZjyV+Sg6EqY+bUaYzKgI+Nnd109HoOmgRCGFZtGikCweewxatoQhQ1SBXtIvbfNRALcDRwJzCAPZJ6E5s0UyqnHj0LMYORK++QY6dgx/bP7yS9yRSXWUUplxd58LFLn7Mne/H1A9bJEsOOCAUMK8Vy+4/nrYfnt47bW4o5LqJpVE8YuZ1QammtkgMzuXcJe2iGRBvXoweDCMGROmYe3cGU49Fb7/Pu7IpLpIJVEcG613BvAzsDlwSCaDAjCzbmb2gZnNNbO+md6fSK7r0gVmzIDzz4d774XiYnj++bijkuogpTmz3f03d//B3fu7+3nARpkMysyKgDsIl+U2B440s+aZ3KdIPlh7bbjhBnjnHahfH/bbD44+Ogx+i2RKhYnCzIrM7Egzu8DMWkRt+5nZeMIAdya1B+a6+8fu/jvwGHBghvcpkjfat4fJk+HKK+GJJ6B583CFVAFc7S45KFmPYgjhCqeGwG1mdj/hprtB7t4mw3FtBsxLeD0/avsfMzvZzErNrHSR/pySaqh2bbjiijA/91ZbwZFHwoEHwuefxx2ZFJpkiaIE2Mvd+wHdgcOA3d39mSzEZeW0/elvJXcf7O4l7l7SuHHjLIQkkptatIDx4+HGG+HVV0PvYvDgMPAtkg7JEsXv7r4cwN1/Az509y+zExbzCYPmKzQBvsjSvkXyTlERnHdeGOxu1w5OOQX22APmzo07MikEyRLFP8xserTMSHg9w8ymZziuScA2ZrZldGnuEcCoDO9TJO9tvXW4jPbee8MpqVatQk9j2bK4I5N8VjPJe82yFsUq3H2pmZ0BvEwoGTLU3WfGFY9IPjGDk06CffaB006DCy6Axx+HoUPDaSqRqqqw1pOZmVdSCCqVdbJBtZ5EyucOw4fDmWfCd9/BxReHRUUGBdJT62mcmZ1pZn9bZcO1zayLmQ0Deq5poCKSOWbwz3+GMiCHHw79+0PbtjBhQtyRST5Jlii6AcuAR83sCzObZWYfE4oDHgnc7O4PZCFGEVlDjRrBww/Dc8+F0h+dOoXB759/jjsyyQcpzUdhZrWARsCv7v5dxqOqIp16EkndDz9A375w113h/ot77w3lQaT6SWeZcdz9D3dfkItJQkSqZv314c47QxXaGjXCZbT/+lcYwxApT0qJQkQKz267wfTpYZ6LoUNDkcFRughdyqFEIVKN1a0LAweGwe2GDUMJkCOOgIUL445MckmlicLMzjCz+tkIRkTiUVICpaUwYACMGBHKgDzyiIoMSpBKj2JjYJKZDY/miCivDpOI5LnateHSS2HKFNhmGzjmmFDGfN68yj8rhS2V+SguBbYhVJM9HphjZteY2dYZjk1EYtC8Obz1FtxySxjwLi4OV0ipyGD1lepVTw58GS1LgfrAk2Y2KIOxiUhMiorg7LOhrAw6dIDTTw9TsM6ZE3dkEodUxijOMrPJwCDgbaClu58GtCMLU6KKSHy23BJeeQWGDIFp00KRwUGDYOnSuCOTbEqlR9EIONjdu7r7E+7+B0BUgny/jEYnIrEzgxNPDGVAunWDPn2gY8eQOKR6SGWM4nJ3/7SC92anPyQRyUWbbgpPPx2KDM6bF66UuuwyWLIk7sgk03QfhYikzAwOOyz0Lo46Cq66Ctq0gXfeiTsyySQlChGpsoYNYdgwePHFUFhwp53gnHPgp5/ijkwyQYlCRFZbt27hyqjTT4dbb4WWLWH06LijknRTohCRNbLeenD77fDGG+Gmvb33hl694Ntv445M0kWJQkTSYpddwpVQffuG01LNm4dyIJL/lChEJG3q1IFrr4WJE2HjjeHgg8PMel99FXdksiaUKEQk7dq2DcnimmtC6fJmzeDBB1VkMF8pUYhIRtSqBf36wdSpIVH07An77AOflntXluQyJQoRyah//APefBP+859QbLBFC7jjDhUZzCdKFCKScTVqwBlnhEtpd9wxPN9tN/jgg7gjk1QoUYhI1jRtCi+9BA88ADNnwvbbw3XXwR9/xB2ZJKNEISJZZRbGK2bNgv33D+MYHTqECZMkNylRiEgsNt4YnngCnnoKvvgCdtgBLr4Yfvst7shkVUoUIhKrgw+G2bPhuOPCPRitW8Pbb8cdlSRSohCR2NWvD0OHwssvhx7FLrvAmWfCjz/GHZmAEoWI5JC99w5XRp15ZriEtkWLkDwkXrEkCjM7zMxmmtlyMytZ5b1+ZjbXzD4ws65xxCci8Vl33VCJ9q23YO21Q4Xa44+HxYvjjqz6iqtHUQYcDLyR2GhmzYEjgGKgG3CnmRVlPzwRiduOO4YroS65BB55JNzd/eSTcUdVPcWSKNx9truXd6vNgcBj7r7E3f8PmAu0z250IpIr6tQJs+hNmgRNmoTZ9Q45BBYsiDuy6iXXxig2A+YlvJ4ftYlINda6NUyYEG7Oe/75UML8/vtVZDBbMpYozOxVMysrZzkw2cfKaSv3v4KZnWxmpWZWumjRovQELSI5q2ZN6NMHpk8PM+mdeCJ07QqffBJ3ZIUvY4nC3fd09xblLCOTfGw+sHnC6ybAFxVsf7C7l7h7SePGjdMZuojksG23hddeC1dFvfNOuDLqtttg2bK4IytcuXbqaRRwhJmtZWZbAtsAE2OOSURyTI0aYZ7umTNh113h7LPDvRezZ8cdWWGK6/LYg8xsPtAJeN7MXgZw95nAcGAW8BLQ2931d4KIlOtvfwtjFg89FCrRtm4NV1+tIoPpZl4Ao0ElJSVeWloadxgiEqOFC8ONesOHQ6tW4U7vdu3ijiq3mdlkdy+pbL1cO/UkIrJaNtwQHn8cRoyARYtCRdq+feHXX+OOLP8pUYhIQenRI5QwP/54GDgwzHnxxhuVfkySUKIQkYKzwQZw333w6quwdGmYTa93b/jhh7gjy09KFCJSsPbYA2bMgHPPhbvuCpfSvvBC3FHlHyUKESlo66wDN90E48fDeuvBvvvCscfC11/HHVn+UKIQkWqhY0d47z24/HJ47LFQBmT4cJUBSYUShYhUG2utBf37w+TJsMUW8M9/wkEHhalYpWJKFCJS7bRqFcp/XH99mBipeXMYMkS9i4ooUYhItVSzJlxwQRjsbt0aTjoJ9twTPv447shyjxKFiFRrf/87jB0L99wT5r1o0QJuvllFBhMpUYhItVejBpx8crhRr0sXOO882GmnUHRQlChERP6nSRN49ln473/ho4+gTRv497/h99/jjixeShQiIgnM4MgjQ+/i0EPhiiugpCSclqqulChERMrRuHHoWYwaBYsXh/swLrwQfvkl7siyT4lCRCSJ/fcPYxX/+hfccEO4tPa11+KOKruUKEREKlGvHtx9d7g6CqBzZzjlFPj++3jjyhYlChGRFHXuDNOnh/sv7rsPiovhuefijirzlChERKpg7bXDHd3vvAP164dTU0cdFSZLKlRKFCIiq6F9+1Azqn9/ePLJUAbk0UcLswyIEoWIyGqqXTtUo50yBbbeOvQsDjgA5s+PO7L0UqIQEVlDxcXw9tth3osxY8LrwYNh+fK4I0sPJQoRkTQoKgoz6ZWVhRv0TjklzLA3d27cka05JQoRkTTaaqswV/e994aJklq2DPdfLF0ad2SrT4lCRCTNzELZ8lmzYO+9wx3dO+4YSprnIyUKEZEM2WwzeOaZMPXqJ59A27ahdtSSJXFHVjVKFCIiGWQWplydNQuOOCJUo23XDiZMiDuy1ClRiIhkQaNG8NBD8PzzofRHp05h3ouff447ssopUYiIZFH37qHI4Kmnhpn0WrYMl9TmMiUKEZEsW399uPNOeP31MHf3nnuG6rTffRd3ZOWLJVGY2fVm9r6ZTTezEWa2QcJ7/cxsrpl9YGZd44hPRCQbdt0Vpk2Diy6CoUNDGZCRI+OO6q/i6lGMBlq4eyvgQ6AfgJk1B44AioFuwJ1mVhRTjCIiGVe3LgwcGAa3GzeGHj3CoPfChXFHtlIsicLdX3H3FbefvAs0iZ4fCDzm7kvc/f+AuUD7OGIUEcmmkhIoLYWrroIRI6BZM3j44dwoMpgLYxQnAi9GzzcD5iW8Nz9qExEpeLVqwSWXwNSpsN12cOyxsO++8Nln8caVsURhZq+aWVk5y4EJ61wCLAUeWdFUzqbKzadmdrKZlZpZ6aJCLgQvItVOs2bw5ptw661hwLu4GO66K74igxlLFO6+p7u3KGcZCWBmPYH9gKPd/9e5mg9snrCZJsAXFWx/sLuXuHtJ48aNM/VjiIjEoqgIzjorFBns2BFOPx123x0+/DD7scR11VM3oA9wgLv/kvDWKOAIM1vLzLYEtgEmxhGjiEgu2HJLeOWVcFXUjBmw/fYwaFB2iwzGNUZxO7AeMNrMpprZ3QDuPhMYDswCXgJ6u/uymGIUEckJZnDCCaEMyD77QJ8+0KFDuLQ2K/v3XBhSX0MlJSVeWloadxgiIlnx1FPQuzd8802Yv/ucc1ZvO2Y22d1LKlsvF656EhGRKjjkkNC7OProMAVrptXM/C5ERCTdGjSABx7Izr7UoxARkaSUKEREJCklChERSUqJQkREklKiEBGRpJQoREQkKSUKERFJSolCRESSKogSHma2CPg07jgq0Aj4Ou4gVpNij08+x6/Y47E6sW/h7pWW3y6IRJHLzKw0lVoquUixxyef41fs8ciUP5DXAAAIY0lEQVRk7Dr1JCIiSSlRiIhIUkoUmTc47gDWgGKPTz7Hr9jjkbHYNUYhIiJJqUchIiJJKVGsBjMbamYLzawsoW17M3vHzGaY2bNmtn7Ce/3MbK6ZfWBmXRPau0Vtc82sb67FbmZ7mdnkqH2ymXVJ+Ey7qH2umd1mZpZLsSe8/zcz+8nMLkhoy+njHr3XKnpvZvR+nag9p4+7mdUys2FR+2wz65fwmTiO++ZmNi6KZaaZnR21NzCz0WY2J3qsH7VbdFznmtl0M2ubsK2e0fpzzKxnDsZ+dBTzdDMbb2bbJ2xrzY69u2up4gLsCrQFyhLaJgG7Rc9PBAZEz5sD04C1gC2Bj4CiaPkI2AqoHa3TPMdibwNsGj1vAXye8JmJQCfAgBeBfXIp9oT3nwKeAC6IXufDca8JTAe2j143BIry4bgDRwGPRc/XBj4BmsZ43DcB2kbP1wM+jH4nBwF9o/a+wMDoeffouBrQEZgQtTcAPo4e60fP6+dY7DuuiAnYJyH2NT726lGsBnd/A1i8SvN2wBvR89HAIdHzAwm/OEvc/f+AuUD7aJnr7h+7++/AY9G6ORO7u09x9y+i9plAHTNby8w2AdZ393c8/E98EOiRS7EDmFkPwi/0zIT1c/64A3sD0919WvTZb9x9WZ4cdwfWMbOaQF3gd+AH4jvuC9z9vej5j8BsYLNo38Oi1Yax8jgeCDzowbvABtFx7wqMdvfF7v4t4Wfulkuxu/v4KDaAd4Em0fM1PvZKFOlTBhwQPT8M2Dx6vhkwL2G9+VFbRe1xqCj2RIcAU9x9CSHO+Qnv5VzsZrYO0Afov8r6+XDctwXczF42s/fM7KKoPeePO/Ak8DOwAPgMuMHdF5MDx93MmhJ6yROAjdx9AYQvZGDDaLWc/H1NMfZEvQg9I0hD7EoU6XMi0NvMJhO6ib9H7eWdQ/Yk7XGoKHYAzKwYGAicsqKpnG3kWuz9gZvd/adV1s+H2GsCOwNHR48Hmdke5Efs7YFlwKaEU63nm9lWxBy7ma1LOA15jrv/kGzVctpi/X2tQuwr1u9MSBR9VjSVs1qVYq9ZlZWlYu7+PuGUAWa2LbBv9NZ8/vwXehNgxemcitqzKknsmFkTYARwnLt/FDXPZ2W3FnIz9g7AoWY2CNgAWG5mvwGTyf3jPh943d2/jt57gTBG8DC5f9yPAl5y9z+AhWb2NlBC+Is2luNuZrUIX7SPuPvTUfNXZraJuy+ITi0tjNor+n2dD+y+SvtrmYwbqhw7ZtYKuI8wdvVN1JzsOyg1mRyMKeSFMECXOLi3YfRYg3Du+MTodTF/Hsz+mDC4VDN6viUrB5iKcyz2DaK4DilnG5MIg30rBlW751Lsq3zmSlYOZufDca8PvEcYDK4JvArsmw/HnfBX7P1RfOsAs4BWcR33KI4HgVtWab+ePw8ID4qe78ufB7MnRu0NgP+L/m3qR88b5FjsfyOMge64yvprfOwz/h+sEBfgUcI52D8I2boXcDbhqoQPgeuIbmaM1r+EcNXBByRcpUK4wuLD6L1Lci124FLC+eapCcuKL4gSwnnqj4DbE3/eXIh9lc9dSZQo8uG4R+sfQxiEL1vxRZAPxx1Yl3CV2UxCkrgw5uO+M+E0y/SE/8PdCVeSjQHmRI8NovUNuCOKcQZQkrCtEwlfxHOBE3Iw9vuAbxPWLU3Xsded2SIikpQGs0VEJCklChERSUqJQkREklKiEBGRpJQoREQkKSUKERFJSolCqhUzq2tmr5vZOmb2vpm1THjvIjO7O4uxNDezaWY2Jarlk+n91Yl+9qJM70sKixKFVDcnAk+7+8/AOcCd0RwEmxFqWfVL+ukUVOGL+GDgSXdv4+6fpGF7Sbn7b4SKr4emY3tSfShRSEEwsx2iCVvqRL2FmWbWopxVjwZGArj7S4Q7jo8Dbgau9KhMs5n1NbOJ0TYvT9jPsxYmcZppZidFbTXN7Dszu8rMJhIK4yXG1tbMJkTbesrM6pnZAcAZwKlm9uoq6/9le2bW38wmmVmZmd1tFiYsMrNtzWxs1DN5b0XPpKL4gWeiYyCSumzcRq9FSzYW4CrgBkIJhn7lvF8b+HKVtk0JZSnGJbR1B+4klHOoAbxEVD+HleUS1iaUqKhPqKXjwMEVxDUL2Dl6fg2h9PaKeM8pZ/2/bI8/l5h4lKgUDKHI4f7R8zpRXMnirwksjPvfSkt+LaoeK4Xk34Sieb8BZ5XzfiPgu8QGd//CzMYCzyU0702YIWxK9HpdwhwR44Fzo94AhCqcWxPq6vxOqLL7J2bWEKjj7m9FTcOAh1L4WVbd3h5mdiEhGTQCJpvZu0Ajd382+ll+i/ZZYfzuvtTM3MzquvuvKcQhokQhBaUB4UuxFuEL9edV3v81al/V8mhZwYCr3H1I4kpmtidhWtCO7v6rmb2VsL1f3b28wmmrO6f1/7ZnZmsTCgC2dffPzeyqhP1WtM+/xJ+gNrBkNeOSakhjFFJIBgOXAY8QJlr6Ew/jD0VmVl6ySPQy0CuaJQ8za2JmjYB6wOIoSRQDO1QWkIc5JX41sx2jpmOB11P9gSJ1CYnsazNbj5VT1X4bte0fxVknSioVxY+ZbUSY+3x5OfsRKZd6FFIQzOw4YKm7/ze6Smi8mXVx97GrrPoKoXzzq3/ZSMTdXzCzfwDvRmPGPxIm5HkeONnMpgHvE6alTMWxwF1mVpeoRHUVfjTc/RszG0YoL/7pKvs9GrjHzK4mnK46JEn8XwOdo59DJGUqMy7Vipm1Ac5z92PjjiUOZjYSON/d58Ydi+QPnXqSasXdpwDjquNNZ2a2FuG+DSUJqRL1KEREJCn1KEREJCklChERSUqJQkREklKiEBGRpJQoREQkqf8HPGwblm605AMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_fit(x_test, y_test, x, y): \n",
    "    plt.plot(x_test, y_test, 'b-')\n",
    "    plt.plot(x, y, 'rx')\n",
    "    plt.ylabel(\"y (Race time)\")\n",
    "    plt.xlabel(\"x (Year of race)\")\n",
    "    plt.show()\n",
    "\n",
    "plot_fit(x_test, predict(x_test, w0, w1), x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll compute the sum of square residuals $SSR(w_0,w_1)$ on the training set to measure the goodness of fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.960618520409062\n"
     ]
    }
   ],
   "source": [
    "def compute_SSR(x, y, w0, w1): \n",
    "    return np.square(y - (w0 + x*w1)).sum() #had to divide by 2 to make this possible\n",
    "\n",
    "print(compute_SSR(x, y, w0, w1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's obvious from the plot that the fit isn't very good. \n",
    "We must repeat the alternating parameter updates many times before the algorithm converges to the optimal weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration #0: SSR = 8.956016522018313\n",
      "Iteration #500: SSR = 6.996549644796752\n",
      "Iteration #1000: SSR = 5.587468721715798\n",
      "Iteration #1500: SSR = 4.574178284979233\n",
      "Iteration #2000: SSR = 3.8455065015753167\n",
      "Iteration #2500: SSR = 3.3215081012361094\n",
      "Iteration #3000: SSR = 2.944693285825805\n",
      "Iteration #3500: SSR = 2.6737203114411288\n",
      "Iteration #4000: SSR = 2.47885973313044\n",
      "Iteration #4500: SSR = 2.3387326777819246\n",
      "Iteration #5000: SSR = 2.237965289062857\n",
      "Iteration #5500: SSR = 2.1655020049097575\n",
      "Iteration #6000: SSR = 2.1133926110175105\n",
      "Iteration #6500: SSR = 2.075919993490542\n",
      "Iteration #7000: SSR = 2.0489728931802604\n",
      "Iteration #7500: SSR = 2.029594844272687\n",
      "Iteration #8000: SSR = 2.015659809514191\n",
      "Iteration #8500: SSR = 2.0056389248190225\n",
      "Iteration #9000: SSR = 1.9984327619486453\n",
      "Iteration #9500: SSR = 1.9932507061827076\n",
      "Final estimates: w0 = 29.93069181464587; w1 = -0.013509936729851188\n"
     ]
    }
   ],
   "source": [
    "for i in np.arange(10000):\n",
    "    w1 = update_w1(x, y, w0) \n",
    "    w0 = update_w0(x, y, w1) \n",
    "    if i % 500 == 0:\n",
    "        print(\"Iteration #{}: SSR = {}\".format(i, compute_SSR(x, y, w0, w1)))\n",
    "print(\"Final estimates: w0 = {}; w1 = {}\".format(w0, w1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try plotting the result again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVNX5x/HPAwvSVNDdGBRwFbBXWBANlkWNEeyFWIgNoyIWNBEl/rBgiRAjxBINsaFEDLagRoMSNho7iwqCiiKKoKhYooLSn98f5+4wLFtml5m5M7vf9+t1XzNz5u6dZy7sffaUe465OyIiIgBN4g5ARERyh5KCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEhCQdwB1FVhYaEXFxfHHYaISF6ZPn36l+5eVNt+eZcUiouLKS8vjzsMEZG8YmbzU9lPzUciIpKgpCAiIglKCiIikqCkICIiCUoKIiKSoKSQDaNGQVnZumVlZaFcRCSHKClkQ48e0L//2sRQVhZe9+gRb1wiIpXk3X0Keam0FCZODIlg0CC4/fbwurQ07shERNahmkK2lJaGhHDNNeFRCUFEcpCSQraUlYUawvDh4bFyH4OISA5QUsiGij6EiRNhxIi1TUlKDCKSY5QUsmHatHX7ECr6GKZNizcuEZFKzN3jjqFOSkpKXBPiiYjUjZlNd/eS2vZTTUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREEjKaFMzsIzN7y8zeNLP1bkO24GYzm2tmM82sWybjERGRmmVjPYVSd/+ymvcOBbpG217A7dGjiIjEIO7moyOB+zx4BWhrZu1jjklEpNHKdFJw4Bkzm25mZ1Xx/lbAgqTXC6OydZjZWWZWbmblixcvzlCoIiKS6aTwM3fvRmgmGmxm+1V636r4mfWmbXX3se5e4u4lRUVFmYhTRETIcFJw90+jxy+Ax4CelXZZCHRMet0B+DSTMYmISPUylhTMrLWZbVzxHPg5MKvSbo8Dp0SjkHoB37r7okzFJCIiNcvk6KMtgMfMrOJzHnD3f5nZOQDufgfwFNAXmAv8AJyewXhERKQWGUsK7j4P2L2K8juSnjswOFMxiIhI3cQ9JFVERHKIkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIQsaTgpk1NbM3zOzJKt47zcwWm9mb0XZmpuMREZHqFWThMy4E3gE2qeb9v7v7eVmIQ0REapHRmoKZdQD6AXdm8nNERCQ9Mt18NAYYCqypYZ9jzWymmT1sZh0zHE/qRo2CsrJ1y8rKQrmISAOVsaRgZocBX7j79Bp2ewIodvfdgCnAuGqOdZaZlZtZ+eLFizMQbRV69ID+/dcmhrKy8LpHj3X3U/IQkQYkkzWFnwFHmNlHwINAHzMbn7yDu3/l7sujl38Fuld1IHcf6+4l7l5SVFSUwZCTlJbCxIkhEVxxRXicODGUJ0s1eYiI5IGMJQV3H+buHdy9GDgBmOruA5L3MbP2SS+PIHRI547SUhg0CK65JjxWTggV+6SSPERE8kDW71MwsxFmdkT08gIzm21mM4ALgNOyHU+Nysrg9tth+PDwWLmZqEIqyUNEJA+Yu8cdQ52UlJR4eXl55j+oohmo4q/+yq+r2nfQoJA8VFMQkRxjZtPdvaS2/XRHc3WmTVv34l7RTDRt2rr7JSeLESPWNiVVV6sQEclhjaamsGwZrF4NrVunOaBRo0KncnLNoKwsJI+hQ9P8YSIi9aOaQiV33w0dO8KwYfDJJ2k88NCh6zcVlZYqIYhIXmo0SaFnT+jTJ/xhX1wMJ520fkuQiEhj12iSQkkJPPwwzJ0L558PTz4ZEkXv3vDII6FpSUSksWs0SaHCNtvATTfBwoUwZgx8+ikcdxx06RLKv/027ghFROKTUlIws3ZmtrOZbWtmDSKRbLIJXHghvP8+PPpo6G/4zW/C45AhMG9e3BGKiGRftRd4M9vUzH5nZm8BrwB/ASYC883sITNrEAPxmzaFo4+G55+H8nI44gi47Tbo2hWOOQb++1/IswFaIiL1VtNf/Q8DC4B93X17d+8dzT/UEbgBONLMBmYlyizp3h3Gj4ePPoLLLoPnnoP99gsjTsePhxUr4o5QRCSzGs19CvXxww9w//2h7+Hdd2HLLWHwYDj7bNh886yEICKSFmm7T8GCAWZ2RfS6k5n1TEeQua5Vq5AAZs+Gp56CXXaByy8P/Q7nnBMShYhIQ5JKp/Gfgb2BE6PX3wO3ZSyiHNSkCRx6KEyeDLNmwcknw733wo47Qt++8Mwz6ncQkYYhlaSwl7sPBpYBuPs3QPOMRpXDdt4Z/vpXWLAgTHX0+utwyCGw665w553w449xRygiUn+pJIWVZtYUcAAzK6Lm5TUbhaKiMKP2/Pmh1lBQAL/+NXTqFJZV+OyzuCMUEam7VJLCzcBjwE/M7DrgBeD6jEaVRzbaCE49Fd54I8yDt88+cO21ITmceiq8+WbcEYqIpK7WpODufwOGAr8HFgFHuftDmQ4s35jBAQfApEkwZ07ooH7kEdhzzzA/3uOPw5pGX78SkVyX6t3JnwP/BV4CWppZt8yFlP+6doVbbglTafzhD+Hu6COPhO23h1tvhSVL4o5QRKRqqQxJvQaYSWhG+mO03ZjhuBqEtm3ht7+FDz6Av/8dCgvDZHwdOsAll8DHH8cdoYjIulKpKfQHOrv7Ae5eGm19Mh1YQ1JQEBZje/nlsB1yCIweDdtuC7/8JbzySrTjqFHrr9hWVhbKRUSyIJWkMAtom+lAGotevUKtYd48uPjicO/D3nuH8rIlPfDkpTwrlvrs0SPeoEWk0ah1mgszKwEmEZLD8opydz8is6FVLZvTXGTDkiVhSOuf/hTWeuhfVMY9P/SnybmDaHHP7euuEy0iUk/pXI5zHDCSMAneH5M2SYM2beC888KIpccfh8W7lHLj0kG0+MM1PL3NIOZ2VEIQkewpSGGfL9395oxH0sg1aQKHHw6Htylj1Ru3M6nTcPaZdju/7FpKm8NLueiiMOTVLO5IRaQhS6WmMN3Mfm9me5tZt4ot45E1RlEfQsGjEzlyxgiaPjyRJ1r1p8lzZfTpA926wbhxsHx57YcSEamPVJLCnkAvwl3MGpKaSdOmrdOHsNmxpbR+ciIPDZ3GnXfCypVw2mmw9dZh3qXFi+MNV0QaHq2nkEfcYcqUMJz16afDFBsDBoTlQ3fZJe7oRCSXpdrRXG1SMLMB7j7ezC6u6n13v2kDY6yXxpwUkr3zThixdN99YWbWgw+Giy4K90A0aRCraItIOqVj9FHr6HHjKrY2GxyhbJAdd4Q77ghTeF93XVjnoW/fMLX3HXeEVeNEROoqlfsUfubuL9ZWli2qKVRtxQp46KHQtDR9Omy2WZiUb/Bg2GqruKMTkbil8z6FW1Iskxg1bx5WhJs2DZ5/HvbfH264AYqLQ7+D8qiIpKLa+xTMbG9gH6CoUr/CJkDTVD8gWqCnHPjE3Q+r9N5GwH1Ad+Ar4Jfu/lHK0ct6zGDffcM2b16YrfWuu+Bvf4PevUO/w5FHQtOU/wVFpDGpqabQnNB3UMC6/QnfAcfV4TMuBN6p5r2BwDfu3gUYTbhzWtJk221Dc9LChXDTTeHx2GPD1N5jxsB338UdoYjkmlT6FLZ29/n1OrhZB8I0GdcBF1dRU5gMXOXuL5tZAfAZUOQ1BKU+hfpbtSosAjR6NLz4Imy8MZx5ZpjOe5tt4o5ORDIpbX0K9U0IkTGEVduqW3NsK2BB9DmrgG+BzSvvZGZnmVm5mZUv1h1b9VZQEGoKL7wAr70WptW45Rbo0mVteZ7dtiIiaZaxEe1mdhjwhbtPr2m3KsrWuyy5+1h3L3H3kqKiorTF2Jj16BH6GT78EIYODTNs7Lsv9OwJDzwQ7p4WkcYnk7c5/Qw4wsw+Ah4E+pjZ+Er7LAQ6AkTNR5sCX2cwJqmkQwf4/e/D/Q633w7ffx9GMRUXh/Kv9a8h0qikshzndmb2bzObFb3ezcz+r7afc/dh7t7B3YuBE4Cp7j6g0m6PA6dGz4+L9lEDRgxat4ZzzoG334Z//hN22gl+97uQNAYNClN7i0jDl0pN4a/AMGAlgLvPJFzk68XMRphZxQI9dwGbm9lc4GLgsvoeV9KjSZNwZ/Szz8LMmXDiiXDPPbDDDtCvX5h7SWlbpOFKJSm0cvfXKpWtqsuHuPt/KkYeufsV7v549HyZux/v7l3cvae7z6vLcRulLK7jvOuu4R6Hjz+Gq68ON8AdfDDstlsoX7Ys7R8pIjFLJSl8aWadiTqAzew4YFFGo5Lq9egR1m3O4jrOP/kJXHFFSA733BNqE2eeCZ06wZVXwuefZ+yjRSTLUrlPYVtgLOHu5m+AD4EBcd15rPsUWJsIBg0KvcNZXsfZPYQwZgw8+SQ0awYnnRSm8N5996yFISJ1kM77FOa5+0FAEbCDu/fWVBQxKy0NCeGaa8JjFhMChKk0+vQJa0rPmQO//nXIS3vsAQceCE88AWuquzNFRHJaKqOPrjeztu6+1N2/N7N2ZnZtNoKTapSVhRrC8OHhsXIfQxZ17Qq33hqm0Bg5Et57D444InRM33YbLFkSW2giUg+p9Ckc6u7/q3jh7t8AfTMXktSooulo4sSwJufEiev2McSkXbtwE9y8eTBhQnh93nnQsSNcemm4D0JEcl8qSaFpNJspAGbWEtiohv0lkyqt40xpaXg9bVrmPrMOI56aNYMTToBXX4WXXgqjlW68McytVFEuIrkrlaQwHvi3mQ00szOAZwmT3Ekchg5dvw+htDSUZ0o9RzztvXfIV/PmhU7op5+GXr1gn33CgkCr6jSwWUSyodbRRwBmdihwIGGuomfcfXKmA6uORh/FJA0jnr7/Pgxp/dOfQqLo1CnM0HrmmdC2bYbiFhEgvSuv4e5Pu/tv3f03cSYEiVEaRjxtvDFccEHojP7HP0KT0iWXhH6HCy6ADz7IQNwiUiepjD7qZWbTzGyJma0ws9VmpuVZGps0jnhq2jSs/vaf/4T1pI8+Gu64I4xkOuooeO45TaUhEpdUagq3AicC7wMtgTPRGs2NSwZHPHXrBvfdB/Pnw+WXhzUdDjgAuncP5StWbHj4IpK6VJuP5gJN3X21u98DZPduKYlXFkY8tW8fWqYWLICxY2H5cjj1VNh6a7j2Wvjyy7R9lIjUIJVpLp4HDgLuJCyXuQg4zd1jmdBAHc2Ngzs880xYOnTyZGjRAgYMCKOYdt457uhE8k86O5p/Fe13HrCUsCjOsRsWnkjNzOCQQ+Bf/4LZs+GUU2D8eNhll7Xl6ncQSb+U1miOprj+zt2vdveLgS2yEJsIEBb8+ctfQtPStdfCW2/BoYeGGsPYsfDjj3FHKNJwVJsUzKypmZ1oZr81s12issPM7CVC57NIVhUWhs7ojz6C+++Hli3h7LPDkNbLL4dPP407QpH8V1NN4S7CSKPNgZvN7B7gRmCUu++ZjeBEqtK8eehfKC8Pw1f33TesJ11cDL/6Fbz+etwRiuSvmpJCCXCwuw8jTIB3PHCAu/8jK5GJ1MIM9tsPHnsM3n8/3FP3j3+E4az77x+er16dwoGyuJqdSK6rKSmscPc1EJbNBN5z98+yE5ZI3XTuHKbPWLgQ/vjHcN/D0UfDdtuF8u+/r+GHY1jNTiRXVTsk1cx+AOZWvAQ6R68NcHffLSsRVqIhqZKKVatCTWH06DBb6yabhDmWzj8/NDOtJ+bV7EQyLR1DUncEDo+2w5JeHxY9iqyVY00wBQVw3HHw4othuu6+fUONoXNnOP74UL7O30Mxr2YnkitqSgofR8NRq9wAzMyyFKfkuhxugunZMyz88+GHYQK+KVOgd2/Ya69QvnIlObWanUis3L3KDfgPcD7QqVJ5c6APYU2F06r7+Uxt3bt3d8lRU6e6Fxa6Dx8eHqdOjTuiKi1Z4n7bbe5du7qD+/GFU31p60L/blIUb8X3yNH4ReoDKPcUrrE11RR+AawGJpjZp2b2tpnNI0yMdyIw2t3vzVSykjxUWxNMjjQxtW4N554L774LTzwBB7WdRr+lE/npiaUMHgzvbZWF1exEclUqmQNoBrQH2qayfyY31RTSaOTI9f8anjo1lNdHbTWFyn+B59Bf5DNmuJ9+unvz5qH20K+f+5Qp7mvWxB2ZSHqQYk0h1gt8fTYlhTRK50U61WPleBPTZ5+5X3mle1FR+O3YbTf3u+92//HHuCMT2TCpJoWUps6WBqpiCuz+/eGKK9aumVCfkTepTq+d46N8ttgCrroKPv4Y7rorjFA644wwhfdVV8Hnn8cdoUiGpZI5cmlTTSEDhg8PfxYPH575z0qlppDuZq0NsGZNaEbq1y+coubNQzPTzJlZD0Vkg5Cu5iPClNntUjlYNjYlhTTLZnNOXZuYcqzv4d133QcNcm/VKvzmHHig+5NPuq9eHWtYIilJZ1K4lnAn80TCiCRL6cDQAngNmAHMBq6uYp/TgMXAm9F2Zm3HVVJIo2xffOtSA8jhvoevvnK/4Qb3rbYKv0HbbReGuC5ZEndkItVLW1IIx8KAQ4AHowRxPdA5hZ9pEz1vBrwK9Kq0z2nAranEULEpKaRRDjXTVCmbzVr1sGKF+wMPuPfoEcJs18790kvdFyyIOzKR9aWaFFJdo9kJS3F+BqwC2gEPm1m1A8yjOJZEL5tFm9bKyiVDh67f0VtaGsrjlgd3GDdrBieeGKbReOEF6NMH/vCHMLfSiSfCa6/FHaFI3dWaFMzsAjObDowCXgR2dfdBQHdqWZYzWqjnTeAL4Fl3f7WK3Y41s5lm9rCZdaz7V5AGp2KKjIkTYcSItSOkcjAxQJjC+2c/g4cfhrlz4cIL4amnwjQaFeWrVsUdpUhqUqkpFALHuPsh7v6Qu68E8DCt9mE1/aC7r3b3PYAOQM+KFdySPAEUe5hxdQph6oz1mNlZZlZuZuWLFy9OIWTJa6kOb81B22wTpu5esADGjIFFi8IEfF26wE03wbffxh2hSM2qnTo77R9kdiWw1N1vrOb9psDX7r5pTcfR1NmST1avhscfDwni+eehTZtw38MFF4QZW0WyJR1TZ29oAEVm1jZ63hI4CHi30j7tk14eAbyTqXhEsqLS/E5Nm8LRbct4rt8oysvhqKPgz3+Grl3DIkDPP19pCm+RmGXyjub2QJmZzQSmEfoUnjSzEWZ2RLTPBWY228xmABcQRiOJ5K8aphDv3h3uvz+sCjdsWEgI++8PJSWhfMWKeEMXgSw2H6WLmo8k56W4itsPP4RkMGZMmLG1fXsYPBjOPhsKC2OIWxq02JuPRBqcVKf+TnF+p1atQgKYPRuefhp23RX+7/+gY8dQ/o4aUyUGSgoiqUp1dbk63mPRpAn84hcweTLMmgUDBsC4cbDTTnDooaE8zyr0ks9SucMtlzbd0SyxytKaEV984T5ihPsWW4S7pXfayX3sWPcffkjT95BGB02dLZIBtTUNpekei6KiUNGYPz/UGpo3h7POgk6dQvmiRWn6PiKVqKNZpC5S7ERON3d47jkYPTosIVpQACecABddBHvumfGPlwZAHc0i6Rbj9BtmcMABMGkSvPcenHMOPPooTOg2iiG7lzFpUrhRLhFnlte9loZDSUEkVTky/UaXLnDzzbBwIXQ/uwfDZ/VnzFFlbL89PHp+GX58FZ3fIilS85FInls9pYyVx/TngU0Gcfgnt3NG64lsf04p558flhEVATUfiTQaTQ8qpcWQQZzxyTUsP30QrQ8rZcyYMLdS//7w8stxRyj5RElBJN8l3RfR4YnbefDsMubNg4svhmeegX32gV694MEHYeXKuIOVXKekIJLPqun87vRBGaNGhX6HW2+Fr74KC/9su23og/7mm7gDl1ylpCCSz2rp/G7TJsynNGdOmMK7a1e49NIwlcZ558H778cYu+QkdTSLNDJvvhkm4ZswITQnHXYYDBkS8olZ3NFJpqijWUSqtMcecO+94W7p4cPhlVfgwAPDTXD33gvLl8cdocRJSUEk21KdbTXDfvpTuPpq+PhjuPPOsI706aeHYawjRsAXX2Q1HMkRSgoi2ZbqbKtZ0qIFDBwIb70VRit16wZXXhnmWRo4MMzcKo2HkoJItlV0BvfvD1dcsXb0UBbmUKqJGRx8MDz1FLz9dqg1TJgQ1nmoKF+zJtYQJQuUFETikOJCPHHZccdw68OCBXD99SFJ9OsX1ni44w5YujTuCCVTlBRE4lDHhXjisvnmYT3pDz+E8ePDENdBg8KQ1mHD4JNPshBEjvTBNBZKCiLZFuNsq/XVvDmcfHK4/eG//w0Vm1GjoLg4lGd0lHgqfTBKHGmjpCCSbanOtprOC12ajmUGvXvDI4/A3LnhBrgnngjX54ryxBTe6ZJKH0yOdd7ntVSWZ8ulTctxSqORpqU9036sSr791v2mm9yLi8PSocXF4fW3327wodc1fHj4gOHDq36/tqVSGzlSXI4z9ot8XTclBWlU0nmhS9exRo6scm3q1b8f6Y884t67d7iybLyx+5Ah7vPm1T/k5OOnFHttiSMV1Xw/Hzmy/sfMAUoKIg1FOi506TxWCrWOadPcTzrJvaDAvUkT96OPdn/+efc1ayodK5ULcKq1nHQlvQzWquKkpCDSEGSzplCXv5BTjGvhQvdhw9w32yxcbbp3dx8/3n358krHqekCnM7Ekep3bIBNUUoKIvku230Kdf28mmodlS6+S5e6Txoy1UcVjnRw33JL9+uvd//yS0/PBbiuF/tUvmM6a2g5QElBJN+ls2073X8h17ZfNRff1VOm+j//6X7QQeHq07Kl+9lnuy8elMULcCrfUTWF/NmUFEQyLNVRPhvYxj9zpvvAge4/bzbVv6DQ/9Z5uC/ftNDX/DsLF+CavqP6FOK/0NdlU1IQyaBU/kKuSw0mhQSzevNCH3faVN9iC/cDmOpfNS30J38z1X/4YcO/TnWfmba+ldrk0EgmJQURqZt0/4VcxwSzbJn7Pfe4D9x2ql/CyMSPLlpUv4+vMaZs1QJyqNYRe1IAWgCvATOA2cDVVeyzEfB3YC7wKlBc23GVFEQyJJ1/1W7AxXDNmrDb4Ye7m7k3b+5+6qnub75Z9zDWk8v9NBmWC0nBgDbR82bRRb9XpX3OBe6Inp8A/L224yopiOSBNF1833vPffBg91atwtWqtNR90iT31avTGGt95dlIptiTwjofAq2A14G9KpVPBvaOnhcAXxKtG13dpqQg0vh8/XXIJx06hKtWly7ut9zi/v33MQeWRyOZciIpAE2BN4ElwMgq3p8FdEh6/QFQWMV+ZwHlQHmnTp0ydc5EJMetWOH+4IPue+0Vrl5t27pfcon7/PkxBpUnI5lSTQoZnSXV3Ve7+x5AB6Cnme1SaRer6seqOM5Ydy9x95KioqJMhCoieaBZM/jlL+GVV+Cll8KKcH/8I2y77dryrKptXYxUZ8TNIRYSSBY+yOxKYKm735hUNhm4yt1fNrMC4DOgyGsIqqSkxMszOnm7iOST+fPhllvgr3+F776DXr3goovgmGOgoCCDH5y8LkZp6fqvc4yZTXf3ktr2y1hNwcyKzKxt9LwlcBDwbqXdHgdOjZ4fB0ytKSGIiFS29dZw442wcCHcfDMsXhxqDZ07h/L//S9DH5yHtYBUZKymYGa7AeMI/QpNgInuPsLMRhDath43sxbA/cCewNfACe4+r6bjqqYgIjVZvRqefBJGj4bnnoPWreH00+HCC6FLl7iji0+qNYWsNR+li5KCiKTqjTdgzBiYMAFWrYLDDw9NS/vvH1aRa0xibz4SEYnbnnvCuHGh3+Hyy+HFF0MrT7duoXz58rgjzD1KCiLS4LVvD9dcAwsWwNixsGIFnHYaFBeH8sWL444wdygpiEij0bIl/PrXMGsWTJ4Me+wBV1wBnTqF8tmz444wfkoKItLomMHPfw5PPx0SwSmnwPjxsMsua8vXrIk7yngoKYhIo7bTTvCXv4SmpeuuC7WIvn1h551D+Q8/xB1hdikpiIgAhYXwu9/BRx/B/fdDq1ZwzjnQsWMo/+STuCPMDiUFEZEkzZvDgAFQXh7uc9hvP7jhhtApPWAATJ8ed4SZpaQgIlIFs5AQHnsM3n8fBg+GSZOgpGRt+erVcUeZfkoKIiK16Nw53AS3cGGYgO/jj8PcStttB3/6U5hzqaFQUhARSdGmm8LFF8PcufDQQ/DTn8KQIaHf4eKLQ39EvlNSEBGpo4ICOO64cIf0q69Cv35hMr7OnUP5Cy9Ans0glKCkICKyAXr2hAceCLWESy6BqVNh333Xlq9cGXeEdaOkICKSBh06hFFKCxbAbbeFfoaTT4ZttgnlX38dd4SpUVIQEUmj1q3h3HPhnXfCFN477ADDhoV+h3PPhTlz4o6wZkoKIiIZ0KRJ6GuYMgVmzAgL/9x1V0gSFeW52O+gpCAikmG77QZ33x2Gsl55ZVic7eCDYffdQ/myZXFHuJaSgohIlmyxBVx1VUgOd98dygYODLO0XnklfP55rOEBSgoiIlnXokVYInTGjNCM1LMnjBgRksPpp8PMmfHFpqQgIhITMzjwwNAhPWcOnHkmTJwYmpUqyrM9hbeSgohIDthuuzCUdcGCMIR1zpywpvQOO4TypUuzE4eSgohIDtlsM7j0Uvjww3DzW9u2cN554T6ICRMy//lKCiIiOahZMzjxxDCNxosvwkEHhRvhMq0g8x8hIiL1ZQb77BO2bFBNQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREEpQUREQkwTwXV3mogZktBubHHUc1CoEv4w6inhR7PPI5dsjv+Btb7Fu7e1FtO+VdUshlZlbu7iVxx1Efij0e+Rw75Hf8ir1qaj4SEZEEJQUREUlQUkh5mwt2AAAIZUlEQVSvsXEHsAEUezzyOXbI7/gVexXUpyAiIgmqKYiISIKSQg3M7G4z+8LMZiWV7W5mL5vZW2b2hJltkvTeMDOba2ZzzOyQpPJfRGVzzeyyXIzfzA42s+lR+XQz65P0M92j8rlmdrOZWS7FnvR+JzNbYma/TSrL+rmvx/+b3aL3Zkfvt4jKc/q8m1kzMxsXlb9jZsOSfiaO897RzMqiWGab2YVR+WZm9qyZvR89tovKLTqvc81sppl1SzrWqdH+75vZqTkY+8lRzDPN7CUz2z3pWBt27t1dWzUbsB/QDZiVVDYN2D96fgZwTfR8J2AGsBGwDfAB0DTaPgC2BZpH++yUg/HvCWwZPd8F+CTpZ14D9gYMeBo4NJdiT3r/EeAh4LfR61jOfR3PewEwE9g9er050DQfzjtwEvBg9LwV8BFQHON5bw90i55vDLwX/V6OAi6Lyi8DRkbP+0bn1YBewKtR+WbAvOixXfS8XY7Fvk9FTMChSbFv8LlXTaEG7v488HWl4u2B56PnzwLHRs+PJPyCLHf3D4G5QM9om+vu89x9BfBgtG/G1SV+d3/D3T+NymcDLcxsIzNrD2zi7i97+F93H3BULsUOYGZHEX55ZyftH8u5r2PsPwdmuvuM6Ge/cvfVeXLeHWhtZgVAS2AF8B3xnfdF7v569Px74B1gq+izx0W7jWPteTwSuM+DV4C20Xk/BHjW3b92928I3/kXuRS7u78UxQbwCtAher7B515Joe5mAUdEz48HOkbPtwIWJO23MCqrrjwu1cWf7FjgDXdfToh1YdJ7ccZfZexm1hq4FLi60v65dO6rO+/bAW5mk83sdTMbGpXn/HkHHgaWAouAj4Eb3f1rcuC8m1kxofb7KrCFuy+CcPEFfhLtlpO/synGnmwgocYDaYhdSaHuzgAGm9l0QjVvRVReVXuv11Ael+riB8DMdgZGAmdXFFVxjLjiry72q4HR7r6k0v75EHsB0Bs4OXo82swOJD9i7wmsBrYkNJn+xsy2JebYzawNoSlxiLt/V9OuVZTF+jtbh9gr9i8lJIVLK4qq2K1OsRfUZWcBd3+XUOXHzLYD+kVvLWTdv7o7ABXNMdWVZ10N8WNmHYDHgFPc/YOoeCFrq6YQY/w1xL4XcJyZjQLaAmvMbBkwnRw597X8v3nO3b+M3nuK0KY/ntw/7ycB/3L3lcAXZvYiUEL4SzWW825mzQgX1b+5+6NR8edm1t7dF0XNQ19E5dX9zi4EDqhU/p9Mxg11jh0z2w24k9DX9FVUXNN1KDWZ7DxpCBuh4yy50+0n0WMTQjvvGdHrnVm3o3keodOnIHq+DWs7fnbOwfjbRrEdW8UxphE64io6PPvmUuyVfuYq1nY0x3bu63De2wGvEzpqC4ApQL98OO+Ev07vieJrDbwN7BbXeY/iuA8YU6n8D6zbWTsqet6PdTuaX4vKNwM+jP5t2kXPN8ux2DsR+i33qbT/Bp/7jP8Hy+cNmEBoL11JyMADgQsJIwPeA24gugEw2v9yQs//HJJGihBGObwXvXd5LsYP/B+hffjNpK3iYlBCaFf+ALg1+TvnQuyVfu4qoqQQ17mvx/+bAYQO8lkVv/T5cN6BNoTRXrMJCeGSmM97b0JTycyk/8N9CSO6/g28Hz1uFu1vwG1RjG8BJUnHOoNw0Z0LnJ6Dsd8JfJO0b3m6zr3uaBYRkQR1NIuISIKSgoiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKUiDZWYtzew5M2ttZu+a2a5J7w01szuyGMtOZjbDzN6I5rbJ9Oe1iL5700x/ljQsSgrSkJ0BPOruS4EhwJ+jOfS3IsztNKzGn05BHS66xwAPu/ue7v5RGo5XI3dfRpjZ9Lh0HE8aDyUFyTtm1iNaXKRFVAuYbWa7VLHrycAkAHf/F+FO3VOA0cBVHk09bGaXmdlr0TGvSPqcJywsODTbzM6MygrM7H9mdq2ZvUaYFC45tm5m9mp0rEfMbFMzOwI4DzjHzKZU2n+945nZ1WY2zcxmmdkdZmFxHTPbzsymRjWO1ytqHNXFD/wjOgciqcvG7efatKV7A64FbiRMUzCsivebA59VKtuSMHVDWVJZX+DPhCkPmgD/IppPhrVTCrQiTOPQjjC3jAPHVBPX20Dv6Pn1hOmkK+IdUsX+6x2PdadhmEA0ZQphgr/Do+ctorhqir8A+CLufytt+bVpllTJVyMIE8YtAy6o4v1C4H/JBe7+qZlNBZ5MKv45YeWqN6LXbQhrHLwEXBT9lQ9htsnOhHlmVhBmk12HmW0OtHD3F6KiccD9KXyXysc70MwuIVz4C4HpZvYKUOjuT0TfZVn0mdXG7+6rzMzNrKW7/5hCHCJKCpK3NiNcAJsRLp5LK73/Y1Re2Zpoq2DAte5+V/JOZnYQYWnKXu7+o5m9kHS8H929qknD6ruGcuJ4ZtaKMPldN3f/xMyuTfrc6j5zvfiTNAeW1zMuaYTUpyD5aiwwHPgbYVGgdXjoL2hqZlUlhmSTgYHR6m2YWQczKwQ2Bb6OEsLOQI/aAvKwJsKPZrZPVPQr4LlUv1CkJSFpfWlmG7N2udRvorLDozhbRAmkuvgxsy0Ia22vqeJzRKqkmoLkHTM7BVjl7g9Eo3VeMrM+7j610q7PEKYknrLeQSLu/pSZ7QC8EvXnfk9YPOafwFlmNgN4l7A0Yip+BdxuZi2Jpl2uw1fD3b8ys3GEKbPnV/rck4G/mNl1hCanY2uI/0ugNPoeIinT1NnSYJnZnsDF7v6ruGOJg5lNAn7j7nPjjkXyh5qPpMFy9zeAssZ4A5eZbUS4L0IJQepENQUREUlQTUFERBKUFEREJEFJQUREEpQUREQkQUlBREQS/h/2fSvHQ0IAcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_fit(x_test, predict(x_test, w0, w1), x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does more than 10 iterations considerably improve fit in this case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Linear algebra solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In lectures, we saw that it's possible to solve for the optimal weights $\\mathbf{w}^\\star$ analytically. The solution is\n",
    "$$\\mathbf{w}^* = \\left[\\mathbf{X}^\\top \\mathbf{X}\\right]^{-1} \\mathbf{X}^\\top \\mathbf{y}$$\n",
    "where\n",
    "$$\\mathbf{X} = \\begin{pmatrix} \n",
    "        1 & x_1 \\\\ 1 & x_2 \\\\ \\vdots & \\vdots \\\\ 1 & x_n \n",
    "    \\end{pmatrix} \n",
    "  \\quad \\text{and} \\quad \n",
    "  \\mathbf{y} = \\begin{pmatrix} \n",
    "          y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_n\n",
    "      \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "We construct $\\mathbf{X}$ in the code block below, remembering to include the $x_0 = 1$ column for the bias (intercept)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.000e+00 1.896e+03]\n",
      " [1.000e+00 1.900e+03]\n",
      " [1.000e+00 1.904e+03]\n",
      " [1.000e+00 1.908e+03]\n",
      " [1.000e+00 1.912e+03]\n",
      " [1.000e+00 1.920e+03]\n",
      " [1.000e+00 1.924e+03]\n",
      " [1.000e+00 1.928e+03]\n",
      " [1.000e+00 1.932e+03]\n",
      " [1.000e+00 1.936e+03]\n",
      " [1.000e+00 1.948e+03]\n",
      " [1.000e+00 1.952e+03]\n",
      " [1.000e+00 1.956e+03]\n",
      " [1.000e+00 1.960e+03]\n",
      " [1.000e+00 1.964e+03]\n",
      " [1.000e+00 1.968e+03]\n",
      " [1.000e+00 1.972e+03]\n",
      " [1.000e+00 1.976e+03]\n",
      " [1.000e+00 1.980e+03]\n",
      " [1.000e+00 1.984e+03]\n",
      " [1.000e+00 1.988e+03]\n",
      " [1.000e+00 1.992e+03]\n",
      " [1.000e+00 1.996e+03]\n",
      " [1.000e+00 2.000e+03]\n",
      " [1.000e+00 2.004e+03]\n",
      " [1.000e+00 2.008e+03]\n",
      " [1.000e+00 2.012e+03]]\n"
     ]
    }
   ],
   "source": [
    "X = np.hstack((np.ones_like(x), x))\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we can express $\\mathbf{w}^\\star$ explicitly in terms of the matrix inverse $(\\mathbf{X}^\\top \\mathbf{X})^{-1}$, this isn't an efficient way to compute $\\mathbf{w}$ numerically. It is better instead to solve the following system of linear equations:\n",
    "$$\\mathbf{X}^\\top\\mathbf{X} \\mathbf{w}^\\star = \\mathbf{X}^\\top\\mathbf{y}$$\n",
    "\n",
    "This can be done in numpy using the command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.solve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "?np.dot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.88952457e+01]\n",
      " [-1.29806477e-02]]\n"
     ]
    }
   ],
   "source": [
    "w = np.linalg.solve(np.transpose(X)@X, np.transpose(X)@y)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting this solution, as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYFNXVx/HvEWRHRSEJOiCiEF8XVARF1OiIS1SCC0JcUFQMigtxF42jIrhAXAhEMe4YF0SjcY24MHGJogyLCIoKKIKi4r6LwHn/uDVNM8z09AzdXd0zv8/z1NPd1TVVZ0qpM3XurXvN3REREQFYL+4AREQkfygpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKYiISIKSgoiIJCgpiIhIgpKCiIgkNIw7gJpq3bq1d+jQIe4wREQKyvTp0z9z9zbVbVdwSaFDhw6UlZXFHYaISEExs0XpbKfykYiIJCgpiIhIgpKCiIgkKCmIiEiCkoKIiCQoKeTC6NFQWrrmutLSsF5EJI8oKeRC9+7Qv//qxFBaGj537x5vXCIiFRTccwoFqbgYJk0KiWDIEBg/PnwuLo47MhGRNehOIVeKi0NCGDEivCohiEgeUlLIldLScIdQUhJeK7YxiIjkASWFXChvQ5g0CS6/fHUpSYlBRPKMkkIuTJu2ZhtCeRvDtGnxxiUiUoG5e9wx1Ei3bt1cA+KJiNSMmU13927Vbac7BRERSVBSEBGRBCUFERFJUFIQEZEEJQUREUlQUhARkQQlBRERSVBSEBGRBCUFERFJyGpSMLP3zewNM5tlZms9hmzBWDObb2azzaxrNuMREZHUcjGfQrG7f1bFdwcCnaJlV2B89CoiIjGIu3x0CHCXB1OBjcysbcwxiYjUW9lOCg48bWbTzWxwJd9vBixO+rwkWiciIjHIdvlod3f/yMx+BTxjZvPc/YWk762Sn1lr2NYooQwGaN++fXYiFRGR7N4puPtH0eunwMPALhU2WQK0S/pcBHxUyX5udvdu7t6tTZs22QpXRKTey1pSMLPmZtay/D2wPzCnwmaPAsdFvZB6AF+7+9JsxSQiIqlls3z0a+BhMys/zr3u/pSZnQLg7jcBTwIHAfOBH4ATshiPiIhUI2tJwd0XAjtUsv6mpPcOnJatGEREpGbi7pIqIiJ5RElBREQSlBRERCRBSUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBRERCRBSUFERBKUFEREJEFJQUREErKeFMysgZnNNLPHK/nueDNbZmazouWkbMcjIiJVa5iDY/wZeAvYoIrv73f303MQh4iIVCOrdwpmVgQcDNyazeOIiEhmZLt8NAY4H1iVYpu+ZjbbzB40s3ZZjid9o0dDaema60pLw/rabCciUgCylhTMrDfwqbtPT7HZY0AHd+8CPAtMqGJfg82szMzKli1bloVoK9G9O/Tvv/qCX1oaPnfvXrvtREQKgLl7dnZsdhVwLLACaEJoU3jI3QdUsX0D4At33zDVfrt16+ZlZWWZDrdy5Rf4IUNg/HiYNAmKi2u/nYhITMxsurt3q267rN0puPuF7l7k7h2AI4EpFROCmbVN+tiH0CCdP4qLw4V+xIjwWtWFPt3tRETyXM6fUzCzy82sT/RxqJnNNbPXgaHA8bmOJ6XS0vCXf0lJeK3YdlDT7URE8lzWykfZkrPyUXlJqLwUVPFzTbcTEYlR7OWjgjdt2poX9uLi8HnatNptJyJSAHSnICJSD+hOQUREakxJQUREEupNUpg3D265BX74Ie5IRETyV1pJwcxamdm2ZtbRzAoykUyaBIMHQ1ERnH8+LFoUd0QiIvmnygu8mW1oZheZ2RvAVOAfwCRgkZk9YGYF1d+ypARefBH22w+uuw46doTDD4f//hcKrK1dRCRrUv3V/yCwGNjT3X/r7nu4ezd3bwdcDRxiZoNyEmUGmMEee8D998P778OwYfDCC6EH6Q47qLQkIgL1vEvqjz/CxIkwdizMmgWtWsGf/gSnngqbb56RQ4iI5IWMdUm1YICZXRJ9bm9mu2QiyLg1bQonnAAzZoS7hn33hWuvDaWlvn3h+edVWhKR+iWdRuMbgd2Ao6LP3wI3ZC2iGJjBnnuGxuj33oMLLggJYe+9Yccd4dZbVVoSkfohnaSwq7ufBvwE4O5fAo2yGlWM2rWDK6+ExYvhtttCwvjTn8L6YcPggw/ijlBEJHvSSQq/RHMdOICZtSH1TGp1QtOmcOKJMHNmKC3tsw9ccw1ssQUccURYp9KSiNQ16SSFscDDwK/M7ArgJeDKrEaVR8pLSw88AAsXhmccSkthr71Caem220KDtYhIXVBtUnD3ewjzLF8FLAUOdfcHsh1YPmrfHq66CpYsCe0MACedFB6Iu/BClZZEpPCl+3TyJ8CLwMtAUzPrmr2Q8l/TpjBoUOjG+vzz4VmH0aNDr6V+/VRaEpHC1bC6DcxsBGFGtAVE7QrR6z7ZC6swmMHvfheWDz6AG28MD8E9+GB4IG7oUDjqqJBEREQKQTp3Cv2BLd19b3cvjpZ6nxAqat8err469Fq65ZZwpzBoUOi1dNFFYb2ISL5LJynMATbKdiB1RbNmoZ1h1qwwrtJee8GoUaHXUr9+YfylKktLo0evPb9zaWlYLyKSA+kkhauAmWY22cweLV+yHVihMwsJ4V//Cr2WzjkHnnsulJq6doU77oCffqrwQ927h/mdyxND+XzP3bvnPH4RqZ+qHfvIzOYSRkh9g6TnE9z9+eyGVrlCno7zhx/g3nvDWEtvvAGtW4cH44YMCWUmYHUiGDIExo9fc/5nEZFayuR0nJ+5+1h3L3X358uXDMRY75SXll5/PVz799xzdWmpf3946SXwvYtDQhgxIrwqIYhIDqWTFKab2VVmtpuZdS1fsh5ZHWYWxlV66KHVpaVnnw1J4uTOpfx0/XhWXFgS7hQqtjGIiGRROuWjyq5KHlcPpEIuH6Xyww8wpaSUnn/rT9+Vk5jTupir9i/lhKf60+BBlZBEZN2kWz6q9jkFd9fVKAeaNYPev56GPz2JS6yYsWPh5InFTGQSJ5w5jQ43FtOzZ7jLEBHJlirvFMxsgLvfbWZnV/a9u1+X1ciqUFfvFCrz/vurH4j76qvQa2noUPjjH6FJk7ijE5FCkomG5ubRa8tKlhbrHKFUq0OH8IjCkiXwj3/Azz/D8ceHnkoXXwwffhh3hCJS16TTprC7u/+vunW5Up/uFCpyD+3O48bBI49AgwZhhrihQ2G33VRaEpGqZbJL6rg010mWmYV5HR5+GBYsgDPPhMmTYffdw/NtEyZU8kCciEgNVJkUoi6o5wBtzOzspOUyoEHOIpRKbbEF/PWvobR0001hTofjjw9jMJWUqLQkIrWT6k6hEaHtoCFrtid8AxyR7gHMrIGZzTSzxyv5rrGZ3W9m883sVTPrUJPgBZo3h5NPhjlzwrMOPXvCFVeE9oijjoJXXtEw3iKSvnTaFDZ390W1PkDovdQN2MDde1f47lSgi7ufYmZHAoe5+x9T7a8+tymka+HC0Gvp1lvh66+hWzc444zQa6lx47ijE5E4ZKxNYR0TQhFwMHBrFZscAkyI3j8I9DJTc+m66tgxzCe9ZEl4KPr772HgwFBauuQS+OijuCMUkXyV7sxrtTWGMJXnqiq+3wxYDODuK4CvgU2yHFO90aIFnHIKzJ0bSks9esDIkbD55nD00TB1qkpLIrKmrCUFM+sNfOru01NtVsm6tS5TZjbYzMrMrGzZsmUZi7G+MINevUI31vnzQxfWJ58M3Vh32QX++c/wDISISLVJwcw6m9lzZjYn+tzFzC5OY9+7A33M7H1gIrCPmd1dYZslQLtovw2BDYEvKu7I3W92927u3q1NmzZpHFqq0rEjXHttKC3deGMoLR13XCgtXXopLF0ad4QiEqd07hRuAS4EfgFw99nAkdX9kLtf6O5F7t4h2n6Kuw+osNmjwMDo/RHRNipo5ECLFmFk7rlz4ZlnYNddw2jd7durtCRSn6WTFJq5+2sV1q2o7QHN7HIz6xN9vA3YxMzmA2cDw2q7X6kdM9h3X3j0UXj33dBL6YknQmlp113h7rtVWhKpT9KaZMfMtiSq9ZvZEUCNigzu/t/y7qjufom7Pxq9/8nd+7n7Vu6+i7svrGH89U8W53Hecku47rrw4NsNN8C338Kxx4aG6csuU2lJpD5IJymcRpiOc2sz+xA4ExiS1aikajmYx7lFCzj1VHjzTXj66bDryy8PyeGYY+DVVzN2KBHJM9U+vJbY0Kw5sJ67f5vdkFLTw2vEMo/z/Pnh7uH22+Gbb0KvpaFDoV8/aNQoq4cWkQzI2MNrZnalmW3k7t+7+7dm1srMRmYmTKmV4tzP47zVVnD99aHX0g03hCelBwxYXVr6+OOshyAiOZBO+ehAd/+q/IO7fwkclL2QpFqlpeEOoST38zi3bLm6tDR5Muy8MwwfHnotDRgAr1XskiAiBSWdpNDAzBIj5phZU0Aj6MSlvHQ0aVIo9E+atGYbQ46stx7svz88/ji8805IFI89Fnos9egB994Ly5fnNCQRyYB0ksLdwHNmNsjMTgSeYfV4RZJr06at2YZQXBw+T5uWvWNW0+OpUycYMyaUlv7+9zB16DHHhNLS8OEqLYkUkrQams3sQKAXYViKp919crYDq4oammOQfHdSXLz25wpWrQoPxI0dG4bTWH/9MELr0KEZ7SQlIjWQbkNz2r2P8oWSQkxq2ePp3XfD3cMdd4TnHnr0CMmhb1/1WhLJpUz2PuphZtPM7DszW25mK83sm8yEKQWjlj2eOnWCv/0tPBA3bhx88UUYRmPzzUOTyCefZDluEamRdNoU/g4cBbwLNAVOQnM01z/r2OOpZUs4/XR46y34z39gp53CAHzt24cB+bLZJCIi6Utr6Gx3nw80cPeV7n4HkP2O8ZI/Mtjjab314Pe/D20Nb78d5nv497/Dw3C77Qb33adeSyJxSicp/GBmjYBZZjbazM4Cmmc5LsknWerx1LlzKC0tWRIapT//PJSWOnQIVSqVlkRyL605moFPgEbAWYQ5D26M7h5yTg3NddeqVeGBuLFj4amnQkP0kUeGkVu7Vds8JiKpZHSO5mg002/cfbi7nw38OiNRiiRZbz048MDQ5jBvHpx8Mjz0UOjG2rMnTJwIv/wSd5QidVuVScHMGpjZUWZ2rpltF63rbWYvExqfRbLmt78NdwwffhhKTMuWwVFHhdLSyJHw6adxRyhSN6W6U7iN0NNoE2Csmd0BXAOMdvedchGcyAYbhOca3n47NE536RI6QLVrBwMHwvRUM4CLSI1V2aYQzcncxd1XmVkT4DNgK3ePddACtSnI22+HB+LuvBO++y6UloYOhcMPD09Pi8jaMtGmsNzdV0GYIQ14J+6EIAKhtDRuXOi1NGZMKCUdeWQoLV1xRSg11UgWZ7MTKTSpksLWZjY7Wt5I+vyGmc3OVYAiVdlwQ/jzn8OdwxNPwPbbw8UXh9LS8cfDjBlp7igHs9mJFIpU5aPNU/2guy/KSkTVUPlIUpk3b3Vp6fvvYffdQ2npsMOqKS3FMJudSC5lonz0QdQdtdIlOohlLGKRDNh665AUPvwwzBT38cdhhNYttoArr0xRWophNjuRfJQqKZSa2Rlm1j55pZk1MrN9zGwCMDC74UnByLO6/IYbwplnhgmAHn8ctt0W/vKXUFo64QSYObPCD8Q4m51IPkmVFH4PrATuM7OPzOxNM1tIGBjvKOB6d78zBzFKIcjTuvx668HBB4cnpd98EwYNggcegK5dYc89Q5VoxTP5MZudSD5Id5Kd9YHWwI/J8zXHQW0Keay6uvzo0SFJJK8rLQ1jKJ1/fs7C/Oqr0OYwbhwsXAgjNxhNh37d2f+qYtq0iS8ukWxKt00Bdy+oZeedd3bJYyUl7hBeK5oyxb116/Ba2eccW7HC/bHH3PfbL4TcuLH7CSe4z5gRSzgiWQWUeRrX2LSGzpY6KtPtANXV5ctHV+3fHy65JOWUnrnQoAH07g1PPw1z58KJJ4ZwunaF3/0ulJlWrIglNJH4pJM58mnRnUIGZfIv95rsK9XdRMy+/NL9uuvcO3YMIRYVuV95pfuyZXFHJrJuyNSdgpmdbmatcpCfJNcy+Zd7unMu5Hkvn402grPOCr2WHn00dHG96CIoKgqN1LNmxR2hSJZVlzWAkcB8YBKhR5Klk22ytehOIQty9Zd7uncTo0atvW7KlLA+BnPnup9yinuzZuE07bmn+wMPuP/ySyzhiNQKmbpTcPeLgU6EUVOPB941syvNbMtsJSrJoVz+5Z7u3USedW/dZptwapYsgWuvDa/9+kHHjnD11fDZZ7GEJZId6WSOkGTYARgDzAPGAzMJw2hXtX0T4DXgdWAuMLySbY4HlgGzouWk6uLQnUIG5VlvoDWUx1JSkj8xRVascH/kEfdevcKdQ5Mm7oMGuc+aFXdkIlUjg20KQ81sOjAa+B+wvbsPAXYG+qb40Z+Bfdx9B2BH4Pdm1qOS7e539x2j5dbq4pEMytLcyxmRx8NONGgAffrAs8/CnDlh8L377oMdd4S99oJ//Uu9lqRwpdMltTVwuLsf4O4PuPsvAB6G1e5d1Q9Fyem76OP60VL9k3KSO+efv/bFtrg4Px7YyvMG6XLbbru6tHTNNfDBB3DEEaG0NGoUfP553BGK1Ew6bQqXeBUjorr7W6l+NprScxbwKfCMu79ayWZ9oyG5HzSzdmlFLXVbeRtCAQ070aoVnHMOzJ8P//43dOoEw4aFXkt/+hPM1mDzUiCy+vCau6909x2BImCX8rmekzwGdHD3LsCzwITK9mNmg82szMzKltV4BhUpOPlc1qpGgwZwyCHw3HPwxhthytB77oEddoC994aHHlJpSfJbWmMfZeRAZpcC37v7NVV83wD4wt03TLUfjX0kea2S8Z2+eaSUmTdPY+Dc81m0CNq3h9NOg5NOgo03jjFWqVcyMZ/CugbQxsw2it43BfYl9FxK3qZt0sc+QMpylEjeq6Q77QYn9Wevc7uzYEEoLW21FVxwgUpLkp+yWT5qS5iTYTYwjdCm8LiZXW5mfaJthprZXDN7HRhK6KIqUrhSPCWeXFqaPRuOPXZ1aam4GB5+WKUliV/OykeZovKRFIRLLgndaUtKQmN5Fb74Am67DW64ARYtgs03D6WlQYNUWpLMir18JFLnpDuqbA260268MZx3HixYEO4UOnYMPYKLimDw4NBYLZJLSgoi6Upn+I1adqdt0AAOPRSmTAmlpQED4O67oUsX2Gef0BaxcmUWfzeRiJKCSLrSGVU2A91pt98ebr4ZFi8OD8AtWACHHQZbbhkekPvyywz/XiJJ1KYgUlNpthdkyooV8NhjMHYs/Pe/0LRpaKQ+4wzYruKTPyJVUJuCSDbEMPxGw4bhTqG0FF5/HY45Bu66K9xR9Oql0pJklpKCSLryYPiNLl3gllvCWEujRkGvGaP522GlbLVVUmlpXaZUlXpPSUEkXXk0/MYmm4ReShc80J2nNujPIRuUct55cHTbUr49qD/vtYln7gkpfGpTECl00R3Mx4cPofld4zli5SSe/qWYXr1Cu0Pv3qF3k9RvalMQqS+iuSd+c/MIWp43hHuXFnP11WGe6UMPDcNqXHutei1JepQURApdhcbvTWaXcsEFsHAhPPhgGIDv3HPDA3FDhsCbb8YdsOQzJQWRQpai8bthQ+jbF55/HmbOhCOPhDvvDBMD7bsvPPqoei3J2pQURApZmo3fO+4YxlhavBiuugrefjsMztepE1x3HXz1VQyxS15SQ7NIPbRiRXi+YexYePFFaNYsTAh0+umwzTZxRyfZoIZmEalSw4ZhLukXXlhdWrr99lBa2m+/8AS1Skv1k5KCSK6lO9pqjpSXlpYsgSuvhHnzoE8f6NwZrr9epaX6RklBJNfSGW01Bq1bw4UXhl5LkybBppvC2WeHXkunngpvaV7EekFJQSTX0hltNUbrrw/9+oW2hhkzQni33x7aGvbfHx5/HFatijtKyRYlBZE4RA+cMWJEeM2ThFDRTjuFhLB4MVxxRXjG4Q9/UGmpLlNSEIlDDKOtros2beCii+C998JNTdu2q0tLp50W2iGyJs/aYOo6JQWRXMuD0VZrK7m0NH16eH/rrfB//wcHHABPPJGF0lI6bTBKHBmjpCCSa+mOtprJC10WLppdu8Idd4TS0siRMHduGHyvc2cYMwa+/rrWu15TOm0wedp4X5DcvaCWnXfe2UXqhSlT3Fu3Dq+VfY5rX1VYvtz9/vvdd9/dHdybN3c/7TT3t97K0AFKSsKOS0oq/778dyopyfjvVhcAZZ7GNTb2i3xNFyUFqVcyeaHL1L5GjVr7Z6dMCesjZWXuAwe6N2oUrjIHHOD+xBPuK1dmOfbqEkc9pqQgUldk8kKXiX3V4K7jk0/cR4xw33TTcNittnIfM8b9q6+iDdJIMGkfL4dJrxApKYjUBbm8U6jJxbCGcS1f7j5xonvPnuGq06KF++mnuy+6M40LfiYTRyb3VWCUFEQKXa7bFGp6vFR3HSkuvtOmuR933OrS0vndp/hPG7T2VRevQ+JLN6Hl+q4jjygpiBS6TJYxanrRrO5iWN12aVx8y0tLbdu6DyckmKkHlPjXX9f816uReto+oaQgIrWTbi+fDP21vXzyFP+xZWu/bbMS/5TWfmCTKX7GGe5vv53B36mietiTSUlBRGounYthTe5gaphg3rpxin/duLXv13CKg/uBB7o/+eQ69FpKdcx1uMtJWx41WispiEjNZLqBdR0SzDcXj/LLLw+lJXDv1Ml97Fhf99JSOr9jJi/kedRoHXtSAJoArwGvA3OB4ZVs0xi4H5gPvAp0qG6/SgoiWZKHF8Off3a/9173Hj3C1aplS/ehQ93feafmIbl7PH+550kpKh+SggEtovfrRxf9HhW2ORW4KXp/JHB/dftVUhApAFm4+L76qvuAAe7rr++J0tJ//pPh0lJNZLKMlgOxJ4U1DgLNgBnArhXWTwZ2i943BD4jmje6qkVJQaR+W7rUffhw99/8JlzBOnd2HzfO/ZtvchxIgXVvzYukADQAZgHfAaMq+X4OUJT0eQHQOtU+lRRExD2Ulu65x33XXT0zpaXayGWj9TpKNylkdZRUd1/p7jsCRcAuZrZdhU2ssh+ruMLMBptZmZmVLVu2LBuhikiBadQIjj4apk6FV1+FQw4JU1N07gwHHwyTJ+dghrjqJktKd0TcPGIhgeTgQGaXAt+7+zVJ6yYDl7n7K2bWEPgYaOMpgurWrZuXlZVlP2ARKTgffwz/+EdIDp98Ar/9LZxxBhx3HLRsmYUDlg/RPWRIOGgeTatakZlNd/du1W2XtTsFM2tjZhtF75sC+wIV52d6FBgYvT8CmJIqIYiIpPKb38Cll8IHH8Ddd8OGG8Lpp4cZ4s48E+bPz+DBCniypFSyWT5qC5Sa2WxgGvCMuz9uZpebWZ9om9uATcxsPnA2MCyL8YhIPdGoERxzTCgrTZ0a5pW+8cZQWurdG55+Gtb5z88CLA2lI2flo0xR+UhEamPp0lBauummUFraeutwF5G10lKeib18JCKST9q2hcsug0WLQmmpZcvVpaWzzspwaamAKSmISL3SuHEoLb32Wigt9e4Nf/97KC394Q/wzDMZKC0VMCUFEam3dt0V7rkn3D2UlIREsf/+sO22oQ3iu+/ijjD3lBREpN7bdFMYPjz0WvrnP6F5czjttFBaOvtsWLAg7ghzR0lBRCTSuDEMGBDuGF55BQ46CMaNg06doE8fePbZul9aUlIQEanADHr0gHvvDaWliy8O7Q/77RdKS+PH193SkpKCiEgKm24ank1bvBgmTIBmzeDUU0Np6ZxzYOHCuCPMLCUFEZE0NG4cnmmYNg1efhkOPBDGjoWttgrjLtWV0pKSgohIDZjBbrvBfffB++/DX/4S2h/22w+22y48HPf993FHWXtKCiIitbTZZmGA1A8+CKWlJk3C2HhFRXDuufDee3FHWHNKCiIi66hJk1BaKiuD//0PDjgA/vY32HLLUFp67rnCKS0pKYiIZIgZ9OwJEyeuWVrad1/Yfvsw9lK+l5aUFEREsiC5tHTnnWHk1lNOyf/SkpKCiEgWNWkCAwfC9Onw0kthGI0xY0Jp6dBDYcqU/CotKSmIiOSAGey+O9x/fygtXXRRaH/o1SuUlm6+OT9KS0oKIiI5VlQEI0eGB+LuuCOUlk4+Gdq1g/POC0kjLkoKIiIxadIEjj8+lJZefDE863D99aG0dNhhYWbPXJeWlBRERGJmBnvssbq0NGxYSBL77ANdusAtt8APP+QmFiUFEZE8UlQEV1wRSku33w4NG8LgwWH9ffdl//hKCiIieahpUzjhBJgxA154ITzrsMUW2T9uw+wfQkREassM9twzLLmgOwUREUlQUhARkQQlBRERSVBSEBGRBCUFERFJUFIQEZEEJQUREUlQUhARkQTzfBrIOw1mtgxYFHccVWgNfBZ3ELWk2ONRyLFDYcdf32Lf3N3bVLdRwSWFfGZmZe7eLe44akOxx6OQY4fCjl+xV07lIxERSVBSEBGRBCWFzLo57gDWgWKPRyHHDoUdv2KvhNoUREQkQXcKIiKSoKSQgpndbmafmtmcpHU7mNkrZvaGmT1mZhskfXehmc03s7fN7ICk9b+P1s03s2H5GL+Z7Wdm06P1081sn6Sf2TlaP9/MxpqZ5VPsSd+3N7PvzOzcpHU5P/e1+P+mS/Td3Oj7JtH6vD7vZra+mU2I1r9lZhcm/Uwc572dmZVGscw1sz9H6zc2s2fM7N3otVW03qLzOt/MZptZ16R9DYy2f9fMBuZh7MdEMc82s5fNbIekfa3buXd3LVUswO+ArsCcpHXTgL2i9ycCI6L32wCvA42BLYAFQINoWQB0BBpF22yTh/HvBGwavd8O+DDpZ14DdgMM+A8gIojhAAAHWUlEQVRwYD7FnvT9v4AHgHOjz7Gc+xqe94bAbGCH6PMmQINCOO/A0cDE6H0z4H2gQ4znvS3QNXrfEngn+nc5GhgWrR8GjIreHxSdVwN6AK9G6zcGFkavraL3rfIs9p7lMQEHJsW+zudedwopuPsLwBcVVv8WeCF6/wzQN3p/COEfyM/u/h4wH9glWua7+0J3Xw5MjLbNuprE7+4z3f2jaP1coImZNTaztsAG7v6Kh//r7gIOzafYAczsUMI/3rlJ28dy7msY+/7AbHd/PfrZz919ZYGcdweam1lDoCmwHPiG+M77UnefEb3/FngL2Cw69oRoswmsPo+HAHd5MBXYKDrvBwDPuPsX7v4l4Xf+fT7F7u4vR7EBTAWKovfrfO6VFGpuDtAnet8PaBe93wxYnLTdkmhdVevjUlX8yfoCM939Z0KsS5K+izP+SmM3s+bABcDwCtvn07mv6rx3BtzMJpvZDDM7P1qf9+cdeBD4HlgKfABc4+5fkAfn3cw6EO5+XwV+7e5LIVx8gV9Fm+Xlv9k0Y082iHDHAxmIXUmh5k4ETjOz6YTbvOXR+srqvZ5ifVyqih8AM9sWGAWcXL6qkn3EFX9VsQ8Hrnf37ypsXwixNwT2AI6JXg8zs14URuy7ACuBTQkl03PMrCMxx25mLQilxDPd/ZtUm1ayLtZ/szWIvXz7YkJSuKB8VSWb1Sj2hjXZWMDd5xFu+TGzzsDB0VdLWPOv7iKgvBxT1fqcSxE/ZlYEPAwc5+4LotVLWH1rCjHGnyL2XYEjzGw0sBGwysx+AqaTJ+e+mv9vnnf3z6LvniTU9O8m/8/70cBT7v4L8KmZ/Q/oRvhLNZbzbmbrEy6q97j7Q9HqT8ysrbsvjcpDn0brq/o3uwTYu8L6/2Yzbqhx7JhZF+BWQlvT59HqVNeh9GSz8aQuLISGs+RGt19Fr+sR6rwnRp+3Zc2G5oWERp+G0fstWN3ws20exr9RFFvfSvYxjdAQV97geVA+xV7hZy5jdUNzbOe+Bue9FTCD0FDbEHgWOLgQzjvhr9M7oviaA28CXeI671EcdwFjKqz/K2s21o6O3h/Mmg3Nr0XrNwbei/7btIreb5xnsbcntFv2rLD9Op/7rP8PVsgLcB+hXvoLIQMPAv5M6BnwDnA10QOA0fZ/IbT8v01STxFCL4d3ou/+ko/xAxcT6sOzkpbyi0E3Ql15AfD35N85H2Kv8HOXESWFuM59Lf6/GUBoIJ9T/o++EM470ILQ22suISGcF/N534NQKpmd9P/wQYQeXc8B70avG0fbG3BDFOMbQLekfZ1IuOjOB07Iw9hvBb5M2rYsU+deTzSLiEiCGppFRCRBSUFERBKUFEREJEFJQUREEpQUREQkQUlBREQSlBSkzjKzpmb2vJk1N7N5ZrZ90nfnm9lNOYxlGzN73cxmRmPbZPt4TaLfvUG2jyV1i5KC1GUnAg+5+/fAmcCN0Rj6mxHGdrow5U+noQYX3cOBB919J3d/PwP7S8ndfyKMbHpEJvYn9YeSghQcM+seTS7SJLoLmGtm21Wy6THAIwDu/hThSd3jgOuByzwaetjMhpnZa9E+L0k6zmMWJhyaa2YnResamtlXZjbSzF4jDAqXHFtXM3s12te/zGxDM+sDnA6cYmbPVth+rf2Z2XAzm2Zmc8zsJrMwuY6ZdTazKdEdx4zyO46q4gf+HZ0DkfTl4vFzLVoyvQAjgWsIwxRcWMn3jYCPK6zblDB0Q2nSuoOAGwlDHqwHPEU0ngyrhxRoRhjGoRVhbBkHDq8irjeBPaL3VxKGky6P98xKtl9rf6w5DMN9REOmEAb4+0P0vkkUV6r4GwKfxv3fSkthLRolVQrV5YQB434ChlbyfWvgq+QV7v6RmU0BHk9avT9h5qqZ0ecWhDkOXgbOiv7KhzDa5JaEcWaWE0aTXYOZbQI0cfeXolUTgH+m8btU3F8vMzuPcOFvDUw3s6lAa3d/LPpdfoqOWWX87r7CzNzMmrr7j2nEIaKkIAVrY8IFcH3CxfP7Ct//GK2vaFW0lDNgpLvflryRme1LmJqyh7v/aGYvJe3vR3evbNCw2s6hnNifmTUjDH7X1d0/NLORScet6phrxZ+kEfBzLeOSekhtClKobgZKgHsIkwKtwUN7QQMzqywxJJsMDIpmb8PMisysNbAh8EWUELYFulcXkIc5EX40s57RqmOB59P9hSJNCUnrMzNryerpUr+M1v0hirNJlECqih8z+zVhru1VlRxHpFK6U5CCY2bHASvc/d6ot87LZraPu0+psOnThCGJn11rJxF3f9LMtgamRu253xImj3kCGGxmrwPzCFMjpuNYYLyZNSUadrkGvxru/rmZTSAMmb2ownGPAf5hZlcQSk59U8T/GVAc/R4iadPQ2VJnmdlOwNnufmzcscTBzB4BznH3+XHHIoVD5SOps9x9JlBaHx/gMrPGhOcilBCkRnSnICIiCbpTEBGRBCUFERFJUFIQEZEEJQUREUlQUhARkYT/B2dz8deCv0W/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w0, w1 = w\n",
    "plot_fit(x_test, predict(x_test, w0, w1), x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should verify that the sum of squared residuals $SSR(w_0, w_1)$, match or beats the earlier iterative result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.9799837689231596\n"
     ]
    }
   ],
   "source": [
    "print(compute_SSR(x, y, w0, w1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** The error we computed above is the *training* error. It doesn't assess the model's generalization ability, it only assesses how well it's performing on the given training data. In later worksheets we'll assess the generalization ability of models using held-out evaluation data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Solving using scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have a good understanding of what's going on under the hood, you can use the functionality in `sklearn` to solve linear regression problems you encounter in the future. Using the `LinearRegression` module, fitting a linear regression model becomes a one-liner as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "lr = LinearRegression().fit(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LinearRegression` module provides access to the bias weight $w_0$ under the `intercept_` property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'lr' 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-70-e33b6092e6ca>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mlr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintercept_\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'lr' is not defined"
     ]
    }
   ],
   "source": [
    "lr.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the non-bias weights under the `coef_` property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.01298065]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should check that these results match the solution you obtained previously. Note that sklearn also uses a numerical linear algebra solver under the hood."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}