worksheet02a_solutions.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": 1,
   "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": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "arange([start,] stop[, step,], dtype=None)\n",
       "\n",
       "Return evenly spaced values within a given interval.\n",
       "\n",
       "Values are generated within the half-open interval ``[start, stop)``\n",
       "(in other words, the interval including `start` but excluding `stop`).\n",
       "For integer arguments the function is equivalent to the Python built-in\n",
       "`range <http://docs.python.org/lib/built-in-funcs.html>`_ function,\n",
       "but returns an ndarray rather than a list.\n",
       "\n",
       "When using a non-integer step, such as 0.1, the results will often not\n",
       "be consistent.  It is better to use ``linspace`` for these cases.\n",
       "\n",
       "Parameters\n",
       "----------\n",
       "start : number, optional\n",
       "    Start of interval.  The interval includes this value.  The default\n",
       "    start value is 0.\n",
       "stop : number\n",
       "    End of interval.  The interval does not include this value, except\n",
       "    in some cases where `step` is not an integer and floating point\n",
       "    round-off affects the length of `out`.\n",
       "step : number, optional\n",
       "    Spacing between values.  For any output `out`, this is the distance\n",
       "    between two adjacent values, ``out[i+1] - out[i]``.  The default\n",
       "    step size is 1.  If `step` is specified as a position argument,\n",
       "    `start` must also be given.\n",
       "dtype : dtype\n",
       "    The type of the output array.  If `dtype` is not given, infer the data\n",
       "    type from the other input arguments.\n",
       "\n",
       "Returns\n",
       "-------\n",
       "arange : ndarray\n",
       "    Array of evenly spaced values.\n",
       "\n",
       "    For floating point arguments, the length of the result is\n",
       "    ``ceil((stop - start)/step)``.  Because of floating point overflow,\n",
       "    this rule may result in the last element of `out` being greater\n",
       "    than `stop`.\n",
       "\n",
       "See Also\n",
       "--------\n",
       "linspace : Evenly spaced numbers with careful handling of endpoints.\n",
       "ogrid: Arrays of evenly spaced numbers in N-dimensions.\n",
       "mgrid: Grid-shaped arrays of evenly spaced numbers in N-dimensions.\n",
       "\n",
       "Examples\n",
       "--------\n",
       ">>> np.arange(3)\n",
       "array([0, 1, 2])\n",
       ">>> np.arange(3.0)\n",
       "array([ 0.,  1.,  2.])\n",
       ">>> np.arange(3,7)\n",
       "array([3, 4, 5, 6])\n",
       ">>> np.arange(3,7,2)\n",
       "array([3, 5])\n",
       "\u001b[0;31mType:\u001b[0m      builtin_function_or_method\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": 3,
   "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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "786.0197711453593\n"
     ]
    }
   ],
   "source": [
    "def update_w0(x, y, w1):\n",
    "    return (y - w1*x).mean()\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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.3998724072997095\n"
     ]
    }
   ],
   "source": [
    "def update_w1(x, y, w0):\n",
    "    return ((y - w0)*x).sum()/(x**2).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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_test = np.arange(1890, 2020)[:, None]\n",
    "\n",
    "def predict(x_test, w0, w1): \n",
    "    return w0 + w1 * x_test"
   ]
  },
  {
   "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": 10,
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5104.580163103217\n"
     ]
    }
   ],
   "source": [
    "def compute_SSR(x, y, w0, w1): \n",
    "    return ((y - w0 - w1*x)**2).sum()\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": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration #0: SSR = 5102.897935763439\n",
      "Iteration #500: SSR = 3670.1237324136036\n",
      "Iteration #1000: SSR = 2639.7951157414004\n",
      "Iteration #1500: SSR = 1898.870931140823\n",
      "Iteration #2000: SSR = 1366.0616521953518\n",
      "Iteration #2500: SSR = 982.9108067637796\n",
      "Iteration #3000: SSR = 707.3815017083582\n",
      "Iteration #3500: SSR = 509.24440022600925\n",
      "Iteration #4000: SSR = 366.76114975678183\n",
      "Iteration #4500: SSR = 264.29938710660775\n",
      "Iteration #5000: SSR = 190.6176541954955\n",
      "Iteration #5500: SSR = 137.6320562088593\n",
      "Iteration #6000: SSR = 99.52934766658781\n",
      "Iteration #6500: SSR = 72.12913980966421\n",
      "Iteration #7000: SSR = 52.42525472313143\n",
      "Iteration #7500: SSR = 38.2559064637811\n",
      "Iteration #8000: SSR = 28.066523547010338\n",
      "Iteration #8500: SSR = 20.73919117263098\n",
      "Iteration #9000: SSR = 15.470000666868614\n",
      "Iteration #9500: SSR = 11.680849695214238\n",
      "Final estimates: w0 = 56.89456037769841; w1 = -0.027293057729518334\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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm8lfP6//HX1W4SSWo7osjwQ5OmHZVjKI4hdIwZkiknhUNykEMdxUFxyHB0RMhwkAiZh7aZaqdpVw7xjZPpRGRIpq7fH59771a7Paxda+17rbXfz8fjfqy17vVZ97rWzb6vPvfnvq+PuTsiIiIAdeIOQEREMoeSgoiIlFJSEBGRUkoKIiJSSklBRERKKSmIiEgpJQURESmlpCAiIqWUFEREpFTduAOorubNm3vr1q3jDkNEJKvMnj37S3fPr6pd1iWF1q1bU1RUFHcYIiJZxcw+SqadTh+JiEgpJQURESmlpCAiIqWUFEREpJSSgoiIlFJSEBGRUkoKIiJSqtYkheXL4fzzYeXKuCMREclctSYpvPQS3HQTtG0L06bFHY2ISGaqNUnh+OPh7behWTPo2xdOPDH0HkREZK1akxQAunWDoiIYNQqmTIE2beDf/wb3NH/x2LFQWLjuusLCsF5EJIPUqqQAUL8+jBwJc+bAzjtD//5w+OHw3/+m8Uu7dYN+/dYmhsLC8LpbtzR+qYhI9dW6pFCiXTt44w24/nqYPj28vu02WLMmDV/WqxdMnhwSwciR4XHy5LBeRCSD1NqkAJCXF65IKi4O/2gfPBh694b330/Dl/XqBUOGwBVXhEclBBHJQLU6KZTYcUd48UW44w6YOxd23x2uuw5+/TWFX1JYCOPHw4gR4bHsGIOISAZQUoiYwcCBsGgRHHQQXHgh9OgB8+enYOMlYwiTJ8Po0WtPJSkxiEiGUVIoY5ttYOpUeOgh+Ogj6No1DAP89NNGbHTWrHXHEErGGGbNSknMIiKpYp726zFTq6CgwGtq5rWvvoKhQ+G++8JNbxMnQvfuNfLVIiIpZWaz3b2gqnbqKVSiWTO491546in49lvo2TMMTP/wQ9yRiYikh5JCEvr0gYULw9VJ48ZBhw6hbIaISK5RUkjS5pvDrbfCK69A3bpwwAFwxhnwzTdxRyYikjpKCtW0zz4wbx5cfDHcfXcYa3jssbijEhFJjbQmBTNbamYLzGyuma03OmzBTWa2xMzmm1mXdMaTKptsAtdcAzNmwFZbwZFHhitMv/gi7shERDZOTfQUerl7pwpGvQ8B/l+0DALG10A8KdO1a7iq9Mor4fHHQ6/h3ntroMCeiEiaxH366I/APR68DWxhZi1ijqla6tWDSy8Nd0LvuiucfDIceih8/HHckYmIVF+6k4IDz5vZbDMbVM772wKJ9UmXReuyTps28NprcOONYTC6XbswMJ2WAnsiImmS7qSwl7t3IZwmOtvM9inzvpXzmfVOvpjZIDMrMrOi5Rk8M05eHpx7biiw16MHnH027LcfvPde3JGJiCQnrUnB3T+NHv8HTAX2KNNkGdAq4XVL4NNytjPB3QvcvSA/Pz9d4abMDjvAc8/BXXfBggWhwN6YMSkusCcikgZpSwpmtqmZNS55DhwIFJdp9gRwcnQVUndgpbt/lq6YapIZnHpqKLDXpw8MHw577hnGHkREMlU6ewq/A143s3nATOApd3/WzAab2eCozdPAh8AS4HbgrDTGE4sWLeDRR8P0n598AgUFYWB69eq4IxMRWZ8K4tWgFStg2DCYNAl22y0U2OvZM+6oRKQ2UEG8DLTlluEu6GefhVWr4Pe/DwPT338fd2QiIoGSQgwOOihcoXT22XDLLdC+PTz/fNxRiYgoKcSmcWO4+WZ49VVo2DAkitNOC6eYRETioqQQs9//PlyRdMkloURG27bwyCNxRyUitZWSQgZo2BCuuirUUWrRAo45Jiyffx53ZCJS2ygpZJDOnWHmzJAgnnwy9BruvlsF9kSk5igpZJh69cKppLlzQ1I47TQ4+GBYujTuyESkNlBSyFC77RYGoW+5Bd58M1yhdPPNKrAnIumlpJDB6tQJl60WF6+9p2GffeDdd+OOTERylZJCFth+e3jmmXAn9KJF0LFjGHf45Ze4IxORXKOkkCXMwgQ+ixdD376hftIee8A778QdmYjkEiWFLPO738HDD4cie59/HhLDJZfAjz/GHZmI5AIlhSx15JHhVNIpp8A110CnTvD663FHJSLZTkkhizVtGiqtvvAC/Pwz7L03nHMOfPdd3JGJSLZSUsgBBxwQZng777wwL3S7dmFgWkSkupQUcsRmm8G4cfDGG+F5nz5hYPqrr+KOTESyiZJCjunRA+bMgcsugwceCHdFP/ywSmWISHKUFHJQgwZwxRVQVAStWkG/fnDUUfBZTsx+LSLppKSQwzp2hLffhrFjw2xvbdrAnXeq1yAiFUt7UjCzPDObY2ZPlvPeqWa23MzmRssZ6Y6ntqlbFy68EObNC0li4EA48ED48MO4IxORTFQTPYXzgMWVvP+Qu3eKljtqIJ5aaZddoLAQxo+HGTOgQ4cwMP3bb3FHJiKZJK1JwcxaAocC2XmwHzs2HEkTFRaG9VmoTh0YPBgWLoR994Xzzw+F9hYtijsyEckU6e4pjAMuAior+Hy0mc03sylm1irN8VRPt25hlLYkMRQWhtfduq1tk4WJo1UreOopuO8+eP/9MLnPFVeEG+BEpHZLW1Iws8OA/7n77EqaTQNau/vuwIvApAq2NcjMisysaPny5WmItgK9esHkySERjBwZHidPDutLJJM4MpAZ9O8feglHHRV+Xrdu4YolEanF3D0tC3A1sAxYCnwOrALuq6R9HrCyqu127drVa9yIEe4QHsszfbp78+bh/ebNw+ss8/jj7tts416njvuFF7qvWhV3RCKSSkCRJ3HsTltPwd0vcfeW7t4aOB6Y7u4nJbYxsxYJL/tS+YB0PEpGZ0eMCI9lTxVB6DkMGRLOwQwZsm5PIkv07RvGGgYOhGuvhd13h1deiTsqEalpNX6fgpmNNrO+0ctzzWyhmc0DzgVOrel4KlVyKmjyZBg9eu2ppPLGEKpKHFlgiy1gwgR46aUw7ed++4Uc9+23cUcmIjUmme5EJi01evpozJj1TwVNnx7WJ75OPGVU9nWW+v5792HDwumkli3dn3wy7ohEZGOQ5Okj8yy7vbWgoMCLMmk0dOzYMEKbeMqosBBmzYKLLoovrhSZMSOcUlq4MAxMjxsHzZvHHZWIVJeZzXb3girbKSlIVX7+OcwJfdVV0KQJ3HwzHHdcuIJJRLJDsklBtY+kSvXrw+WXw+zZsMMOcMIJcMQR8MkncUcmIqmmpCBJ69AB3noLrrsuzPbWti3cfrsK7InkEiUFqZa8PLjgApg/H7p0gUGDYP/94YMP4o5MRFJBSUE2yM47h0tXb7stnFbq0AGuv14F9kSynZKCbLA6dUJPYeHC0Fu44ALo2ROKi+OOTEQ2lJKCbLSWLeGJJ8L0nx9+GE4rjRqlAnsi2UhJQVLCDI4/HhYvhmOPDVcrde0KM2fGHZmIVIeSgqRU8+Zw//0wbRp8/TX06BFOK61aFXdkIpIMJQVJi8MOC2MNf/pTGIDu0CFrS0KJ1CpKCpkkCyfsqUyTJvCvf4WfUKcO9O4dBqZXrow7MhGpiJJCJsnSCXuqst9+MG8eXHghTJwYbnqbNi3uqESkPEoKmSSZmd6yVKNGocMzYwY0axbmbzjhBKjJifREpGpKCpkmBybsqUxBQZjyc/RoeOQRaNMmDEyrVIZIZlBSyDQ5MmFPZerXDz9vzpxwZ/RJJ8Hhh8N//xt3ZCKSVFIws6Zm1s7MdjQzJZJ0SXamtxzRrh288QbccEP4ie3ahYHpNWvijkyk9qrwAG9mTczsr2a2AHgbuA2YDHxkZg+bWW6d18gEs2atO4ZQMsYwa1a8caVRXh4MHQoLFsAee4QzZr17w/vvxx2ZSO1U4SQ7ZvYCcA8wzd2/KfNeV2AAsMDdJ6Y9ygSaZCd3ucNdd8GwYfDTT6GzdP75ULdu3JGJZL+NnmTH3f/g7veWTQjRe7PdfWgyCcHM8sxsjpk9Wc57DczsITNbYmYzzKx1VduT3GUGp58OixbBQQeF2Ux79AhlukWkZlQ5PmDBSWY2Mnq9nZntUY3vOA9YXMF7A4Gv3X1n4AZgTDW2Kzlqm21g6tRw5uzjj0MNpZEjQ+9BRNIrmUHjW4EewAnR6++AfyazcTNrCRwK3FFBkz8Ck6LnU4D9zTTzr4Rew7HHhl7DCSeEK3Q7dw4zv4lI+iSTFPZ097OB1QDu/jVQP8ntjwMuAiq6nmRb4L/Rdn8FVgLNkty21ALNmsE998DTT8P338Nee4WB6R9+iDsykdyUTFL4xczyAAcws3wqPsiXMrPDgP+5++zKmpWzbr2RbzMbZGZFZla0XLfA1kqHHBIK7J11Ftx4I7RvDy++GHdUIrknmaRwEzAV2MrM/g68DlyVxOf2Avqa2VLgQaC3md1Xps0yoBWAmdUFmgArym7I3Se4e4G7F+Tn5yfx1ZKLGjeGW26BV1+FevXgD3+AgQPhm/UuhRCRDVVlUnD3+wmngK4GPgOOcPeHk/jcJe7e0t1bA8cD0939pDLNngBOiZ4fE7VRwYPK5Fgl1Q2x996hwN7w4TBpUiiw99hjcUclkhuSvTv5C+A14E1gEzPrsqFfaGajzaxv9HIi0MzMlgDDgOEbut1aI0crqVbXJpvA1VeHAntbbQVHHhl2wxdfxB2ZSHar8Oa10gZmVwCnAh+w9ny/u3vv9IZWPt28xtpEMGRIqI+UI5VUN9Qvv8C114Z5oTfdFMaNgwEDwhVMIhJs9M1rCfoBO7n7fu7eK1piSQgSyfFKqtVVrx789a8wd26ounrKKdCnT7jHQUSqJ5mkUAxske5ApBpqQSXVDdGmDbz2Gtx0U3hs1w7++U8V2BOpjmSSwtXAHDN7zsyeKFnSHZhUII5Kqlk0uF2nDvz5z1BcHEpknHMO7Lsv/Oc/cUcmkh2SSQqTCOUnrgH+kbBIHOKopJqFg9utW8Nzz4UCe8XF0LEjXHNNGH8QkYolM9D8irvvW0PxVEkDzTHJ4sHtzz+Hs8+GRx8NpTImTgyPIrVJKgeaZ5vZ1WbWw8y6lCwpiFGySRYPbm+9dZj6c8oU+PTT0MG59FJYvTruyEQyTzJJoTPQnXAXc8mpo+vSGZRkoBwY3D766FBgb8AAuOoq6NQpzPwmImslc0dzr3IWXZJam+TQNKFbbhnGGZ57LvQU9t4bzj03FNsTkcqn4zwpehxW3lJzIUrscnCa0AMPDAPQ55wT6im1bw/PPx93VCLxq6ynsGn02LicZbM0xyWZ5KKL1h9D6NUrrM9im2229p6Ghg3DbG+nnQYr1ivJKFJ7VDj7rbvfFj190d3XOfNqZnulNSqRGrTXG2OZf2M3Rr3aizFj4JlnYPKQQvbZZFbWJz6R6kpmoPnmJNeJZKdu3ah/Uj/+fkAhRUXQt3EhbS7vx9+e7MZnn8UdnEjNqrCnYGY9gJ5AfpkxhM2BvHQHJlJjSsZI+vWj05Ah3PbNeB7802TG3NOLm9rCDTeEekoqsCe1QWU9hfqEsYO6rDue8C1h7gORtbKoFEa5Eu7DsCFDOGFCL+bNCwPQp50WxhuWLo07SJH0qzApuPsr7j4K6O7uoxKW6939/RqMUbJBFpbCWEc592Hsuiu88kooqvfWWyFB3HyzCuxJjnP3rFq6du3qkqGmT3dv3tx9xIjwOH36+m3GjFl//fTpYX1cSuIuiavsa3dfutT94IPdwb1nT/dFi2KKVWQDAUWexDE22ZnXJBel+pRPMqUwMrFHkcR9GNtvD08/DffcA+++G+6G/vvfVWBPclAymSOTFvUUUiiJfyFv0PYq6ylUp12G+vxz9379Qq+hY0f32bPjjkikaiTZU6i6AewCvAQUR693By5LZuPpWJQUUixVB+jqJpgRI8L/fiNGrP9eJp5iKsfUqe5bb+2el+d+8cXuq1bFHZFIxVKZFF4B9gDmJKwrTuJzDYGZwDxgITCqnDanAsuBudFyRlXbVVJIg8oO0MmqzoG8qkSU6h5MGq1Y4T5wYNh9u+zi/uqrcUckUr5UJoVZ0WNiUpibxOcM2Cx6Xg+YQbiSqWxSuCWZQEsWJYUUq+lTOcke8LPsFNMLL7i3bh3+os46y33lyrgjEllXskkhmYHmL81sJyAc6c2OAaq8zzOKo6T2ZL1oqXxGH6lZcVQ/Tba4XpbN33DAAaHA3tCh4YrW9u1DuQyRrFNV1gB2BF4EVgGfAK8DrZPJOIQ7n+cC3wNjynn/VEKCmQ9MAVpVtU31FFIok8/dZ1lPIdGbb7q3bRt6DQMGuH/5ZdwRiaTw9JGvPYBvCjROtn2Zz24BFALty6xvBjSIng8Gplfw+UFAEVC03XbbpW2nSYbIojGF9USJdvXqkM/q1nU/osl0n3viGF+zJu7gpDZLNilUefrIzK4ysy3c/Qd3/87MmprZldXsjXwDvAwcXGb9V+7+U/TydqBrBZ+f4O4F7l6Qn59fna+WbJTN8zdE92E0eLOQ0aPh3fGF3PlDP4b+uxtHHRWmAxXJZMmMKRwSHdQBcPevgT5VfcjM8s1si+j5JsABwLtl2rRIeNkXWJxM0JLjsnn+hoTieowcyU6X9GPzZybTZ2wvnn0W2raFiRPBNbomGSqZpJBnZg1KXkQH+AaVtC/RAig0s/nALOAFd3/SzEabWd+ozblmttDM5gHnEsYYRDJTsneAlxkkzzugFxdeCPPnQ8eOcMYZ8Ic/wIcf1lzoIkmr6vwScBFhcHkgcHr0/KJkzk2lY9FAs8QmBZfT/vab+/jx7o0buzdq5H7DDe6//lqDv0FqLVI50AwcAlwH/AM4KJnPpGtRUpBYpejGu48/dj/00PAXuOee7sXFNRS/1FrJJoWkCuK5+zPu/hd3v8Ddn0t9f0UkS1R1/0SSg+StWsG0aXD//bBkCXTuHDb588819DtEKpDM1UfdzWyWmX1vZj+b2W9m9m1NBCeSccqZd2Ed1RgkN4MTT4TFi+Hoo2HkSCgoqOZFVtk+uZFknGR6CrcAJwDvA5sAZ6A5mqU2StMd4Pn58MAD8Pjj8NVX0L17yCGrViXx4UwsRS5ZLdnTR0uAPHf/zd3vAjK75oBIOqT5/om+fWHRIhg4EK69Nlyp9PLLVXyozCWwpUkrw8uCSOZKJimsMrP6wFwzG2tm5xPubhapXWrg/okmTWDCBHjppTDtZ69eMHgwrFxZyYeyrE6UZLZkksKAqN05wA9AK+DodAYlUtv17g0LFsAFF8Dtt0O7dvDUUxU0rmqcQ6QaqkwK7v6Ru69292/dfZS7DwN+VwOxidRqjRrBddfBW29B06Zw2GHQvz8sX57QKI5Kt5LTKkwKZpZnZieY2V/MrH207jAze5Mw+CwiG6oaVw3tsQfMng2XXw4PPxxKZTz4YFQqI5vrRElGMq+gCIuZ3U04VTQT2BP4COgBDHf3x2oqwLIKCgq8qKgorq8XSY3Ef+H36rX+6woUF4eB6Jkz4fDDw9mibbetwbgla5nZbHcvqLJdJUmhGNjd3deYWUPgS2Bnd/88taFWj5KC5IySRDBkSDi6J3nV0G+/wY03wmWXQb164UqlM86AOkldSyi1VbJJobL/jX529zUA7r4aeC/uhCCSUzbwqqG8PBg2LAxEd+0KZ54J++8f7ozOSbpBr0ZVlhR2M7P50bIg4fWCqPKpiGyMjbxqaKedwqWrt98O77wDu+8O//hH6EmsJ5sPrMneoJfNvzGTVFQUCdi+siWZwkrpWFQQT3JCMoXzqjFd6rJl7ocfHgrsdevmvmDBBnxfJktmetZs/41pxsZWSSUab6j0w0m0SfWipCA5IZkDfjUPcmvWuD/4oHt+vnu9eu5/+5v7Tz+V2X4q5r2Oa27vESPCIWvEiIrbpOI3ZvLc5RshFUnhZeDPwHZl1tcHegOTgFOT+ZJULkoKUqtswEFu+XL3/v3DX3e7du5vv53wZjIH1mRjSsW/yJM9AFdnP1T2G9OQjLNFKpJCQ+As4A3gU2AR8CHh0tTbgU7JfEGqFyUFqXU28ED+5JPuLVu6m7mff777qqeqOLBW51/IVR2kq3uwr+wAXJ2DdIrmu0hZryqDbHRSWKcR1CNMr7lFMu3TuSgpSK2ykQenlSvdhwxx34/p/mWd5j77uhQdfN0rT1apPJCnMsEk833J/L4slNKkkEmLkoLUGik8jfHBmWO8/zbTHdzPOMP96699407TVGfgd2NP+SSrOj2dqr4v28dfyhF7UohOP80E5gELgVHltGkAPAQsAWYAravarpKC1BopPqCsWuV+0UXudeq4b7ON++OPV9Aw2QNmMsmqpgaHqyNVp5iq810ZMD6RCUnBgM2i5/Wig373Mm3OAv4VPT8eeKiq7SopiGycWbPcO3QIf/3HHef+xRcJbyZzgE7l4HBNHzST+b5kf186BsnTKGVJgVAyu2kyG6tkG42Ad4A9y6x/DugRPa9LKKVR6WWuSgoiG++nn9yvuMK9fn33Zs3c77vPfc1LMfwLuaZPr6Ty+1LdY0qzVCaFK6PTO5OBg6tzbwKQB8wFvgfGlPN+MdAy4fUHQPPKtqmkIJI6Cxe6d+8ejgQTdx3jXzyYogNmBp1LT6tUj62kUaqvPjLgIODBKEFcBeyUzGejz28BFALty6xfWE5SaFbO5wcBRUDRdtttl9YdJ1Lb/Pqr+7hx7o0auTdu7H7rre6//RZ3VFkkVVdhpVmySSHZOZod+DxafgWaAlPMLKmiIu7+DeFmuIPLvLWMUJ4bM6sLNAFWlPP5Ce5e4O4F+fn5yXyliCQpLw/OOy+U5d5zTzjrrFCb7/33444sC1RVvyob57uoKmsA5wKzCef/jwXqRevrAB9U8rl8ovsagE2A14DDyrQ5m3UHmidXFY9OH4mkz5o17hMnujdp4t6wYTjb88svcUeVoTKoF5AMUthTaA4c5e4HufvD7v5LlEzWAIdV8rkWQGFUUXUW8IK7P2lmo82sb9RmItDMzJYAw4DhScQjImliBqefDosWwcEHw8UXQ/fuMG9e3JFloGzsBSShwkl2MpUm2RGpGe4wZQqccw6sWAHDh4eJfRo0iDsy2RCpmGRHRGoxMzj22NBrOPFEuPJK6NwZ3nor7sgknZQURKRSzZrBpEnwzDPwww+w114wdCh8/33ckUk6KCmISFIOPjhcoXTWWWGO6A4d4IUX4o5KUk1JQUSS1rgx3HILvPoq1K8PBx4IAwfC11/HHZmkipKCiFTb3nuHK5KGDw+nltq2halT445KUkFJQUQ2SMOGcPXVMHMmbL01HHUU9OsHX3wRd2SyMZQURGSjdOkSEsNVV8ETT0CbNnDPPeGSVsk+SgoistHq1YNLLoG5c0NSOOUUOOQQ+OijuCOT6lJSEJGU2W03eO01uPlmeP11aN8e/vlPWLMm7sgkWUoKIpJSdeqEu6CLi6Fnz/B8333hP/+JOzJJhpKCiKRF69bw7LNw992wcCF07AjXXAO//BJ3ZFIZJQURSRuzML6waBEcfngYd9hzT5gzJ+7IpCJKCiKSdltvDQ8/DI88Ap9+Ct26wV//CqtXxx2ZlKWkICI15qijYPFiOPnkcI9Dp07wxhtxRyWJlBREpEY1bQp33gnPPRd6CnvvDX/+M3z3XdyRCSgpiEhMDjwwXKH05z+Hy1bbtw+JQuKlpCAisdlss1Bx9fXXoVGjUIn11FPDpD4SDyUFEYldz57hiqRLL4X77w93RU+ZEndUtZOSgohkhIYNw+xus2ZBy5Zh1rejj4bPPos7stolbUnBzFqZWaGZLTazhWZ2Xjlt9jOzlWY2N1pGpiseEckOnTrBjBnhRrenngplue+6SwX2ako6ewq/Ahe4exugO3C2mbUtp91r7t4pWkanMR4RyRJ168LFF8P8+WGGt9NPh4MOgqVL444s96UtKbj7Z+7+TvT8O2AxsG26vk9Ecs8uu8DLL4erk956K1yhdNNN8NtvcUeWu2pkTMHMWgOdgRnlvN3DzOaZ2TNm1q4m4hGR7FGnTpgXeuFC2GcfOO+8cG/D4sVxR5ab0p4UzGwz4BFgqLt/W+btd4Dt3b0jcDPwWAXbGGRmRWZWtHz58vQGLCIZabvtwhjDvfeGiqudOsHf/64Ce6mW1qRgZvUICeF+d3+07Pvu/q27fx89fxqoZ2bNy2k3wd0L3L0gPz8/nSGLSAYzg5NOCr2EI46Ayy6DggKYPTvuyHJHOq8+MmAisNjdr6+gzdZRO8xsjyier9IVk4jkhq22gocegqlTYfnyUHl1+HD48ce4I8t+6ewp7AUMAHonXHLax8wGm9ngqM0xQLGZzQNuAo5314VnIpKcI44IZblPPRXGjAlzNrz6atxRZTfLtmNwQUGBFxUVxR2GiGSYl16CP/0J/u//wsD01VfD5pvHHVXmMLPZ7l5QVTvd0SwiOWH//WHBAjj/fBg/Ply++vTTcUeVfZQURCRnbLopXH89vPkmNG4Mhx4KAwbAl1/GHVn2UFIQkZzTvTu88w6MHAkPPhhKZUyerFIZyVBSEJGc1KABjBoVLlfdfns47jg48sgwHahUTElBRHLa7ruHEhnXXhsm8WnbFiZOVK+hIkoKIpLz6taFv/wlDER36gRnnAEHHAAffhh3ZJlHSUFEao2dd4bp0+G228K8De3bww03qMBeIiUFEalV6tSBQYPCTW+9e8OwYbDXXqHgnigpiEgt1bIlTJsG//43fPABdO4Mo0fDzz/HHVm8lBREpNYygxNOCL2GY46Bv/0tFNibNSvuyOKjpCAitV5+fugxPPEErFgR7nO48EJYtSruyGqekoKISOTww8PYwp/+BNddFy5nffnluKOqWUoKIiIJmjSBf/0rXKUE0KsXnHkmrFwZb1w1RUlBRKQcvXrB/Pnh/oY77oB27eDJJ+OOKv2UFEREKtCoUbgT+q23oGnTcHrpxBPDxD65SklBRKQKe+wRaiiNGgVTpoRSGQ88kJulMpQURESSUL9+qLo6Zw7stFPoMfTtC8uWxR2+TlbgAAALvElEQVRZaikpiIhUQ7t28MYbYd6Gl14KrydMgDVr4o4sNZQURESqKS8vzPBWXBxudjvzzDDz25IlcUe28dKWFMyslZkVmtliM1toZueV08bM7CYzW2Jm882sS7riERFJtR13hBdfhNtvD5P6dOgQ7m/49de4I9tw6ewp/Apc4O5tgO7A2WbWtkybQ4D/Fy2DgPFpjEdEJOXMQinuRYvgwAPDndA9e4Yy3dkobUnB3T9z93ei598Bi4FtyzT7I3CPB28DW5hZi3TFJCKSLttuC489Fqb/XLoUunQJtZR++inuyKqnRsYUzKw10BmYUeatbYH/JrxexvqJQ0QkK5iFaT8XLYLjjw9VV7t2hRllj3wZLO1Jwcw2Ax4Bhrr7t2XfLucj6135a2aDzKzIzIqW5/JdIyKSE5o3h3vvhaeeCuUxevQI8zb88EPckVUtrUnBzOoREsL97v5oOU2WAa0SXrcE1ptW290nuHuBuxfk5+enJ1gRkRTr0ycU2Bs8OMzw1qFDuIw1k6Xz6iMDJgKL3f36Cpo9AZwcXYXUHVjp7p+lKyYRkZq2+eZw663wyithrugDDghVWL/5Ju7IypfOnsJewACgt5nNjZY+ZjbYzAZHbZ4GPgSWALcDZ6UxHhGR2OyzD8ybBxddBHfeGUplPP543FGtzzzLincUFBR4UVFR3GGIiGywoiIYODBUYT3uOLjpJthqq/R+p5nNdveCqtrpjmYRkRpWUBASw5VXwtSp0KYN3HdfZhTYU1IQEYlBvXpw6aUwdy7suisMGACHHgoffxxvXEoKIiIxatMGXnsNbrwxDEa3awfjx8dXYE9JQUQkZnl5cO65ocBe9+5w1lmw337w3ns1H4uSgohIhthhB3j++XB10oIF0LEjjB1bswX2lBRERDKIGZx2WiiVccghcPHFsOee4XLWmqCkICKSgVq0gEcfDdN/fvJJuGJp3Lj0f6+SgohIBjv66NBr6N8/TAOabnXT/xUiIrIxttwS7r67Zr5LPQURESmlpCAiIqWUFEREpJSSgoiIlFJSEBGRUkoKIiJSSklBRERKKSmIiEiprJt5zcyWAx/FHUcFmgNfxh3EBlLs8cnm+BV7PDYk9u3dPb+qRlmXFDKZmRUlM91dJlLs8cnm+BV7PNIZu04fiYhIKSUFEREppaSQWhPiDmAjKPb4ZHP8ij0eaYtdYwoiIlJKPQURESmlpFAFM7vTzP5nZsUJ6zqa2VtmtsDMppnZ5gnvXWJmS8zsP2Z2UML6g6N1S8xseKbFbmZ/MLPZ0frZZtY74TNdo/VLzOwmM7NMij3h/e3M7Hsz+0vCuoze79F7u0fvLYzebxitz+j9bmb1zGxStH6xmV2S8Jk49nsrMyuMYlloZudF67c0sxfM7P3osWm03qL9usTM5ptZl4RtnRK1f9/MTsnA2PtHMc83szfNrGPCtjZu37u7lkoWYB+gC1CcsG4WsG/0/HTgiuh5W2Ae0ADYAfgAyIuWD4AdgfpRm7YZFntnYJvoeXvgk4TPzAR6AAY8AxySSbEnvP8I8DDwl+h1Nuz3usB8oGP0uhmQlw37HTgReDB63ghYCrSOcb+3ALpEzxsD70V/k2OB4dH64cCY6HmfaL8a0B2YEa3fEvgwemwaPW+aYbH3LIkJOCQh9o3e9+opVMHdXwVWlFm9K/Bq9PwF4Ojo+R8JfyQ/ufv/AUuAPaJlibt/6O4/Aw9GbTMmdnef4+6fRusXAg3NrIGZtQA2d/e3PPxfdw9wRCbFDmBmRxD+eBcmtM/4/Q4cCMx393nRZ79y99+yZL87sKmZ1QU2AX4GviW+/f6Zu78TPf8OWAxsG333pKjZJNbuxz8C93jwNrBFtN8PAl5w9xXu/jXhNx+cSbG7+5tRbABvAy2j5xu975UUNkwx0Dd6fizQKnq+LfDfhHbLonUVrY9DRbEnOhqY4+4/EeJclvBexsVuZpsCFwOjyrTPhv2+C+Bm9pyZvWNmF0XrM36/A1OAH4DPgI+B69x9BRmw382sNaH3OwP4nbt/BuHgC2wVNcvIv9ckY080kNDjgRTErqSwYU4Hzjaz2YSu3s/R+vLO+Xol6+NQUewAmFk7YAxwZsmqcraRabGPAm5w9+/LtM+G2OsCvwf6R49Hmtn+ZEfsewC/AdsQTpdeYGY7EnPsZrYZ4VTiUHf/trKm5ayL9e+1GrGXtO9FSAoXl6wqp1m1Yq9bncYSuPu7hG4/ZrYLcGj01jLW/Zd3S6DklExF62tUJbFjZi2BqcDJ7v5BtHoZa7umkJmx7wkcY2ZjgS2ANWa2GphN5u/3ZcAr7v5l9N7ThHP695H5+/1E4Fl3/wX4n5m9ARQQ/qUay343s3qEg+r97v5otPoLM2vh7p9Fp4f+F62v6O91GbBfmfUvpzNuqHbsmNnuwB2EsaavotWVHYOSk87Bk1xZCINniQNvW0WPdQjnek+PXrdj3YHmDwkDP3Wj5zuwdvCnXYbFvkUU19HlbGMWYSCuZMCzTybFXuYzl7N2oDkb9ntT4B3CQG1d4EXg0GzY74R/nd4VxbcpsAjYPa79HsVxDzCuzPprWXewdmz0/FDWHWieGa3fEvi/6L9N0+j5lhkW+3aEMcueZdpv9L5P+/9g2b4ADxDOmf5CyMIDgfMIVwe8B1xDdBNg1P5Swuj/f0i4WoRwpcN70XuXZlrswGWE88NzE5aSg0EB4bzyB8Atib83E2Iv87nLiZJCNuz3qP1JhAHy4pI/+mzY78BmhKu9FhISwoUx7/ffE06VzE/4f7gP4Yqul4D3o8cto/YG/DOKcQFQkLCt0wkH3SXAaRkY+x3A1wlti1K173VHs4iIlNJAs4iIlFJSEBGRUkoKIiJSSklBRERKKSmIiEgpJQURESmlpCA5y8w2MbNXzGxTM3vXzDokvHeRmf2rBmNpa2bzzGxOVNsm3d/XMPrteen+LsktSgqSy04HHnX3H4ChwK1RDf1tCbWdLqn000moxkH3KGCKu3d296Up2F6l3H01obLpManYntQeSgqSdcysWzS5SMOoF7DQzNqX07Q/8DiAuz9LuFP3ZOAG4HKPSg+b2XAzmxltc2TC90yzMOHQQjM7I1pX18y+MbMrzWwmoShcYmxdzGxGtK1HzKyJmfUFzgEGm9mLZdqvtz0zG2Vms8ys2Mz+ZRYm1zGzXcxsetTjeKekx1FR/MBj0T4QSV5N3H6uRUuqF+BK4DpCmYJLynm/PvB5mXXbEEo3FCas6wPcSih5UAd4lqieDGtLCjQilHFoSqgt48BRFcS1CPh99PwqQjnpkniHltN+ve2xbhmGB4jKpRAK/B0ePW8YxVVZ/HWB/8X930pLdi2qkirZajShYNxq4Nxy3m8OfJO4wt0/NbPpwJMJqw8kzFw1J3q9GWGOgzeB86N/5UOoNrkToc7Mz4Rqsusws2ZAQ3d/PVo1Cbg3id9Sdnv7m9mFhAN/c2C2mb0NNHf3adFvWR19Z4Xxu/uvZuZmtom7/5hEHCJKCpK1tiQcAOsRDp4/lHn/x2h9WWuipYQBV7r7xMRGZnYAYWrK7u7+o5m9nrC9H929vKJhGzqHcun2zKwRofhdF3f/xMyuTPjeir5zvfgT1Ad+2sC4pBbSmIJkqwnACOB+wqRA6/AwXpBnZuUlhkTPAQOj2dsws5Zm1hxoAqyIEkI7oFtVAXmYE+FHM+sZrRoAvJLsD4psQkhaX5pZY9ZOl/p1tO7wKM6GUQKpKH7M7HeEubbXlPM9IuVST0GyjpmdDPzq7v+OrtZ508x6u/v0Mk2fJ5QkfnG9jUTc/Wkz2w14OxrP/Y4wecxTwCAzmwe8S5gaMRkDgPFmtglR2eVq/DTc/Sszm0Qomf1Rme/tD9xmZn8nnHI6upL4vwR6Rb9DJGkqnS05y8w6A8PcfUDcscTBzB4HLnD3JXHHItlDp48kZ7n7HKCwNt7AZWYNCPdFKCFItainICIipdRTEBGRUkoKIiJSSklBRERKKSmIiEgpJQURESn1/wHL1MSHDOa7wgAAAABJRU5ErkJggg==\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": 14,
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "Solve a linear matrix equation, or system of linear scalar equations.\n",
       "\n",
       "Computes the \"exact\" solution, `x`, of the well-determined, i.e., full\n",
       "rank, linear matrix equation `ax = b`.\n",
       "\n",
       "Parameters\n",
       "----------\n",
       "a : (..., M, M) array_like\n",
       "    Coefficient matrix.\n",
       "b : {(..., M,), (..., M, K)}, array_like\n",
       "    Ordinate or \"dependent variable\" values.\n",
       "\n",
       "Returns\n",
       "-------\n",
       "x : {(..., M,), (..., M, K)} ndarray\n",
       "    Solution to the system a x = b.  Returned shape is identical to `b`.\n",
       "\n",
       "Raises\n",
       "------\n",
       "LinAlgError\n",
       "    If `a` is singular or not square.\n",
       "\n",
       "Notes\n",
       "-----\n",
       "\n",
       ".. versionadded:: 1.8.0\n",
       "\n",
       "Broadcasting rules apply, see the `numpy.linalg` documentation for\n",
       "details.\n",
       "\n",
       "The solutions are computed using LAPACK routine _gesv\n",
       "\n",
       "`a` must be square and of full-rank, i.e., all rows (or, equivalently,\n",
       "columns) must be linearly independent; if either is not true, use\n",
       "`lstsq` for the least-squares best \"solution\" of the\n",
       "system/equation.\n",
       "\n",
       "References\n",
       "----------\n",
       ".. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,\n",
       "       FL, Academic Press, Inc., 1980, pg. 22.\n",
       "\n",
       "Examples\n",
       "--------\n",
       "Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``:\n",
       "\n",
       ">>> a = np.array([[3,1], [1,2]])\n",
       ">>> b = np.array([9,8])\n",
       ">>> x = np.linalg.solve(a, b)\n",
       ">>> x\n",
       "array([ 2.,  3.])\n",
       "\n",
       "Check that the solution is correct:\n",
       "\n",
       ">>> np.allclose(np.dot(a, x), b)\n",
       "True\n",
       "\u001b[0;31mFile:\u001b[0m      /usr/lib64/python3.6/site-packages/numpy/linalg/linalg.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.linalg.solve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.88952457e+01]\n",
      " [-1.29806477e-02]]\n"
     ]
    }
   ],
   "source": [
    "w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting this solution, as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.979983768923164\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": 19,
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([28.89524574])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the non-bias weights under the `coef_` property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.01298065]])"
      ]
     },
     "execution_count": 21,
     "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}