13-OrdinaryLeastSquares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares\n",
    "\n",
    "If you recall, in an early notebook, we introduced the notion of fitting a line, or a linear model, to some data. We can then use the definition of that line to predict new data points. However, what we didn't dig into at that point was how we learn and evaluate our linear model, and how to learn the best model.\n",
    "\n",
    "To do so, we need a way to measure how good our model is, or an error measurement, that we can use to evaluate our model. Together with a procedure to update models, we can try and learn models that minimize error - that is to say, models that best fit the data. \n",
    "\n",
    "Ordinary least squares is on such approach for learning and evaluating models. OLS seeks to minimize the sum squared errors. Squared errors are calculated as the square of the difference between the model prediction of a data point, and the data point itself. One way to think about this is as an error function - OLS defines how we will calculate the error of a model, given the data. The model with the lowest error, defined in terms of OLS, is the best model. When we talk about fitting a model with OLS, we mean finding the solution that has the lowest OLS error - the lowest value for the sum of squared errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "Ordinary least squares (OLS) means minimizing the error of the sum of squares between the predictions made by the model, and the observed data. \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "Find more information on OLS on\n",
    "<a href=\"https://en.wikipedia.org/wiki/Ordinary_least_squares\" class=\"alert-link\">wikipedia</a>\n",
    "check out this cool \n",
    "<a href=\"http://setosa.io/ev/ordinary-least-squares-regression/\" class=\"alert-link\">interactive tool</a>\n",
    "and/or check out this \n",
    "<a href=\"https://www.datarobot.com/blog/ordinary-least-squares-in-python/\" class=\"alert-link\">tutorial</a>\n",
    "about doing OLS in Python.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will create a minimal dataset, and explore fitting a simple linear model fit to it, using OLS.\n",
    "\n",
    "In this case, we will be using `numpy` for measuring least squares. Note that for real datasets, this is unlikely to be how you apply models, since it will usually be more practical to use `scikit-learn` or `statsmodels` to manage all the components of model fitting full datasets, but the underlying math is all the same. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create some data\n",
    "# x is an evenly space array of integers\n",
    "x = np.arange(0, 6)\n",
    "\n",
    "# y is some data with underlying relationship y = (theta) * x\n",
    "#  For this example, the true relation of the data is y = 2x\n",
    "true_rel = 2\n",
    "y = true_rel * x\n",
    "\n",
    "# Add some noise to the y dimension\n",
    "noise = np.random.normal(0, 1, len(x))\n",
    "y = y + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAM50lEQVR4nO3dX2isd53H8c9ncyJOj0pcOitNTtm4ILnxwsggux7worWbrhabi72oUFERzs2uW3eXSHMl3glZRC9EONSqi6VlabNniwvGoi1ScKtzTo6mbcwqbtWTdPeMSFYrA03j14szOSRp/s7zzDzzTd4vCEmeM5nn+1D6Zvg9zzzjiBAAIJ8/qXoAAEB3CDgAJEXAASApAg4ASRFwAEjqTD93duutt8b4+Hg/dwkA6V2+fPnXEVHfvb2vAR8fH1ez2eznLgEgPdu/2Gs7SygAkBQBB4CkCDgAJEXAASApAg4ASfX1KhQAOE0uLa5qbmFFa+ttjY7UNDM1oenJsdKen4ADQA9cWlzV7PyS2hubkqTV9bZm55ckqbSIs4QCAD0wt7ByM95b2hubmltYKW0fBBwAemBtvX2s7d0g4ADQA6MjtWNt7wYBB4AemJmaUG14aMe22vCQZqYmStsHJzEBoAe2TlRyFQoAJDQ9OVZqsHdjCQUAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEjq0IDbftj2ddvPb9v2p7afsv3Tzve39nZMAMBuR3kF/jVJd+/a9qCk70TEOyR9p/M7AKCPDg14RHxP0m92bb5X0tc7P39d0nTJcwEADtHtGvjbIuJlSep8/7P9Hmj7gu2m7War1epydwCA3Xp+EjMiLkZEIyIa9Xq917sDgFOj24D/n+3bJKnz/Xp5IwEAjqLbgD8p6aOdnz8q6T/KGQcAcFRHuYzwUUnflzRh+5rtT0j6nKS7bP9U0l2d3wEAfXTmsAdExIf3+ac7S54FAHAMvBMTAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABI6tD7gQNAWS4trmpuYUVr622NjtQ0MzWh6cmxqsdKi4AD6ItLi6uanV9Se2NTkrS63tbs/JIkEfEusYQCoC/mFlZuxntLe2NTcwsrFU2UHwEH0Bdr6+1jbcfhCDiAvhgdqR1rOw5HwAH0xczUhGrDQzu21YaHNDM1UdFE+XESE0BfbJ2o5CqU8hBwAH0zPTlGsEvEEgoAJEXAASApAg4ASRFwAEiqUMBt/6PtF2w/b/tR228sazAAwMG6DrjtMUn/IKkREe+UNCTpvrIGAwAcrOgSyhlJNdtnJN0iaa34SACAo+g64BGxKulfJP1S0suS/j8ivr37cbYv2G7abrZare4nBQDsUGQJ5a2S7pX0dkmjks7avn/34yLiYkQ0IqJRr9e7nxQAsEORJZT3S/qfiGhFxIakeUnvLWcsAMBhigT8l5L+0vYtti3pTknL5YwFADhMkTXw5yQ9LumKpKXOc10saS4AwCEK3cwqIj4j6TMlzQIAOAbeiQkASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUoUCbnvE9uO2f2J72fZflTUYAOBgZwr+/RclfSsi/tb2GyTdUsJMAIAj6Drgtt8i6X2SPiZJEfGqpFfLGQsAcJgiSyh/Iakl6au2F20/ZPvs7gfZvmC7abvZarUK7A4AsF2RgJ+R9G5JX46ISUm/l/Tg7gdFxMWIaEREo16vF9gdAGC7IgG/JulaRDzX+f1x3Qg6AKAPug54RPyvpF/ZnuhsulPSi6VMBQA4VNGrUD4p6ZHOFSg/l/Tx4iMBAI6iUMAj4qqkRkmzAACOgXdiAkBSBBwAkiLgAJAUAQeApAg4ACRFwAEgKQIOAEkRcABIioADQFJF30oPoEuXFlc1t7CitfW2Rkdqmpma0PTkWNVjIRECDlTg0uKqZueX1N7YlCStrrc1O78kSUQcR8YSClCBuYWVm/He0t7Y1NzCSkUTISMCDlRgbb19rO3AXgg4UIHRkdqxtgN7IeBABWamJlQbHtqxrTY8pJmpiX3+Ang9TmJiIJy2KzK2ju00HTPKR8BRudN6Rcb05NiJPj70HksoqBxXZADdIeCoHFdkAN0h4KgcV2QA3SHgqBxXZADd4SQmKscVGUB3CDgGAldkAMfHEgoAJEXAASApAg4ASRFwAEiqcMBtD9letP3NMgYCABxNGa/AH5C0XMLzAACOoVDAbZ+T9EFJD5UzDgDgqIq+Av+CpE9L+sN+D7B9wXbTdrPVahXcHQBgS9cBt32PpOsRcfmgx0XExYhoRESjXq93uzsAwC5FXoGfl/Qh2y9JekzSHba/UcpUAIBDdR3wiJiNiHMRMS7pPknfjYj7S5sMAHAgrgMHgKRKuZlVRDwj6ZkyngsAcDS8AgeApAg4ACRFwAEgKQIOAEkRcABIioADQFIEHACSIuAAkBQBB4CkCDgAJEXAASApAg4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKQIOAAkRcABICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgqa4Dbvt220/bXrb9gu0HyhwMAHCwMwX+9jVJ/xwRV2y/WdJl209FxIslzQYAOEDXr8Aj4uWIuNL5+XeSliWNlTUYAOBgpayB2x6XNCnpuT3+7YLtpu1mq9UqY3cAAJUQcNtvkvSEpE9FxG93/3tEXIyIRkQ06vV60d0BADoKBdz2sG7E+5GImC9nJADAURS5CsWSviJpOSI+X95IAICjKPIK/Lykj0i6w/bVztcHSpoLAHCIri8jjIhnJbnEWQAAx8A7MQEgKQIOAEkRcABIioADQFIEHACSKnIzK/TIpcVVzS2saG29rdGRmmamJjQ9yW1mAOxEwAfMpcVVzc4vqb2xKUlaXW9rdn5Jkog4gB1YQhkwcwsrN+O9pb2xqbmFlYomAjCoCPiAWVtvH2s7gNOLgA+Y0ZHasbYDOL0I+ICZmZpQbXhox7ba8JBmpiYqmgjAoOIk5oDZOlHJVSgADkPAB9D05BjBBnAollAAICkCDgBJEXAASIqAA0BSBBwAkiLgAJAUAQeApAg4ACRFwAEgqYF/JyYfbgAAexvogPPhBgCwv4FeQuHDDQBgfwMdcD7cAAD2N9AB58MNAGB/hQJu+27bK7Z/ZvvBsobawocbAMD+uj6JaXtI0pck3SXpmqQf2n4yIl4sazg+3AAA9lfkKpT3SPpZRPxckmw/JuleSaUFXOLDDQBgP0WWUMYk/Wrb79c623awfcF203az1WoV2B0AYLsiAfce2+J1GyIuRkQjIhr1er3A7gAA2xUJ+DVJt2/7/ZyktWLjAACOqkjAfyjpHbbfbvsNku6T9GQ5YwEADtP1ScyIeM3230takDQk6eGIeKG0yQAAB3LE65ate7czuyXpF13++a2Sfl3iOBlwzKcDx3zyFT3eP4+I151E7GvAi7DdjIhG1XP0E8d8OnDMJ1+vjneg30oPANgfAQeApDIF/GLVA1SAYz4dOOaTryfHm2YNHACwU6ZX4ACAbQg4ACSVIuC9vu/4oLH9sO3rtp+vepZ+sH277adtL9t+wfYDVc/Ua7bfaPsHtn/UOebPVj1Tv9gesr1o+5tVz9IPtl+yvWT7qu1mqc896GvgnfuO/7e23Xdc0ofLvO/4oLH9PkmvSPrXiHhn1fP0mu3bJN0WEVdsv1nSZUnTJ/y/sSWdjYhXbA9LelbSAxHxXxWP1nO2/0lSQ9JbIuKequfpNdsvSWpEROlvXMrwCvzmfccj4lVJW/cdP7Ei4nuSflP1HP0SES9HxJXOz7+TtKw9bk18ksQNr3R+He58DfarqRLYPifpg5IeqnqWkyBDwI9033GcDLbHJU1Keq7aSXqvs5RwVdJ1SU9FxIk/ZklfkPRpSX+oepA+Cknftn3Z9oUynzhDwI9033HkZ/tNkp6Q9KmI+G3V8/RaRGxGxLt041bM77F9opfLbN8j6XpEXK56lj47HxHvlvQ3kv6us0RaigwB577jp0BnHfgJSY9ExHzV8/RTRKxLekbS3RWP0mvnJX2osyb8mKQ7bH+j2pF6LyLWOt+vS/p33VgWLkWGgHPf8ROuc0LvK5KWI+LzVc/TD7brtkc6P9ckvV/ST6qdqrciYjYizkXEuG78f/zdiLi/4rF6yvbZzol52T4r6a8llXZ12cAHPCJek7R13/FlSf920u87bvtRSd+XNGH7mu1PVD1Tj52X9BHdeEV2tfP1gaqH6rHbJD1t+8e68SLlqYg4FZfVnTJvk/Ss7R9J+oGk/4yIb5X15AN/GSEAYG8D/wocALA3Ag4ASRFwAEiKgANAUgQcAJIi4ACQFAEHgKT+CP5tfQFxuz9uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data\n",
    "plt.plot(x, y, '.', ms=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observing the data above, we can see that there is some relation between the x and y dimension. \n",
    "\n",
    "We would like to measure what that relation is. That's where OLS comes in. \n",
    "\n",
    "OLS is a procedure to find the model (in this case, line) that minimizes the squared distances between each observed data point and the model prediction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape that data to play nice with numpy\n",
    "x = np.reshape(x, [len(x), 1])\n",
    "y = np.reshape(y, [len(y), 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit an OLS Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy provides us with a function to calculuate the OLS solution. In this case, we are fitting the model:\n",
    "\n",
    "$$y = \\theta x $$\n",
    "\n",
    "Note that we are not fitting an intercept here (no 'b' value, if you think of 'y = ax + b'). \n",
    "\n",
    "In this simple model, we are therefore implicitly assuming an intercept value of zero. \n",
    "\n",
    "You can fit intercepts (and linear models with more parameters) with OLS, you just need to add them in. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the (Ordinary) Least Squares best fit line using numpy\n",
    "#   This gives us a fit value (theta), and residuals (how much error we have in this fit)\n",
    "theta, residuals, _, _ = np.linalg.lstsq(x, y, rcond=None)\n",
    "\n",
    "# Pull out theta value from array\n",
    "theta = theta[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated theta value is: 2.1194\n"
     ]
    }
   ],
   "source": [
    "# Check what the OLS derived solution for theta is:\n",
    "print('Estimated theta value is: {:1.4f}'.format(theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The true relationship between y & x is: \t2.0000\n",
      "OLS calculated relationship between y & x is: \t2.1194\n"
     ]
    }
   ],
   "source": [
    "# Check how good our OLS solution is\n",
    "print('The true relationship between y & x is: \\t{:1.4f}'.format(true_rel))\n",
    "print('OLS calculated relationship between y & x is: \\t{:1.4f}'.format(theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks good! The absolute error between the true value, and our estimate is quite small!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The residuals for the model fit are: \t3.2984\n"
     ]
    }
   ],
   "source": [
    "# Check what the residuals are. Residuals are the error of the model fit\n",
    "print('The residuals for the model fit are: \\t{:1.4f}'.format(residuals[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVzVVf7H8dcBVMBdSE0RF1xwYVNERc21XDNbrCxrrBzHaWwZ22zKrGaabLOZ0snMX1mjZo1pi1vupuKGguYuqClpKqgoArKd3x8HiOWy3nu5XPg8Hw8eci/fe++5qG++nO/5fI7SWiOEEML5uDh6AEIIIcpHAlwIIZyUBLgQQjgpCXAhhHBSEuBCCOGk3Cryxby9vXWrVq0q8iWFEMLp7dmzJ15rfVPB+ys0wFu1akVkZGRFvqQQQjg9pdQvlu6XKRQhhHBSEuBCCOGkJMCFEMJJVegcuCXp6enExcWRmprq6KEIO3J3d8fHx4caNWo4eihC2NWczbEE+tQn3M+7yGMiYuPZH5fIpH5+Vr2WwwM8Li6OunXr0qpVK5RSjh6OsAOtNQkJCcTFxdG6dWtHD0cIuwr0qc/kRVHMeiDEYohHxMbnft1aDp9CSU1NxcvLS8K7ClNK4eXlJb9liWoh3M+bWQ+EMHlRFBGx8fm+lje8iztDLy2HBzgg4V0NyN+xqE4shbitwxsqSYCXZM7m2EI/yQqKiI1nzubYChqREEIUL2+Iz1xz1ObhDU4S4DlzSkWFeM5PtkCf+uV6fldXV4KDg+ncuTNBQUHMnDmTrKysYh9z6tQpFi1aVK7XE0JUD+F+3ozr4csHG2IY18PXpuENThLg9p5T8vDwIDo6moMHD7J27VpWrlzJa6+9VuxjJMCFECWJiI1nwc7TPDmwLQt2ni5xJqHMtNYV9tGtWzdd0KFDhwrdV5RtMRd1yOtr9LaYixZvl1ft2rXz3Y6NjdWNGjXSWVlZ+uTJk7pPnz46JCREh4SE6G3btmmtte7Ro4euV6+eDgoK0jNnzizyOPG7svxdC+HsbJlXQKS2kKlOFeBa//5NeO/HIzYJb60LB7jWWjdo0ED/9ttv+vr16zolJUVrrfWxY8d0znvYuHGjHjFiRO7xRR0nficBLqqLosK6vCFeVIA7fB14WeWdU3pyYFubzynl0Nl7haanpzN58mSio6NxdXXl2LFjFo8v7XFCiKqtuGndvNPBtrig6XQBXnBOqaefl81D/MSJE7i6utK4cWNee+01mjRpwr59+8jKysLd3d3iY95///1SHSeEqNr2xyUWG845Ib4/LrF6BXjBn2w9/bxsvjTn4sWLTJo0icmTJ6OUIjExER8fH1xcXPj888/JzMwEoG7duly7di33cUUdJ4SoXkpTHh/u5111CnlKw9KvJcWtTimLlJSU3GWEgwcP5rbbbmP69OkAPP7443z++ef07NmTY8eOUbt2bQACAwNxc3MjKCiI999/v8jjhBDCXlTOXG9FCA0N1QU3dDh8+DAdO3Ys9nElLRW0R4WTsL3S/F0LIQpTSu3RWocWvN8pzsDLMqckhBDVhVPMgVfknJIQQthcVha42P582SnOwIUQwilpDdGLYHYYXLdxFSYS4EIIYR/nD8Fnw+HbP4NHQ7hx1eYv4RRTKEII4TRuJMHmt2DHf6BWXRj1IQSPs8sUigS4EELYgtZwZDmsmgpX4yDkIRj8GtT2sttLyhQKZlu3O+64g3bt2uHn58dTTz1FWloaAJs2bWLkyJGFHrN8+XJCQkIICgqiU6dOfPzxx4WOmT9/Pkop1q9fn3vfsmXLUEqxZMmSUo+vqDGU5phNmzZRv359goODCQ4OZvDgwQDMmTOHL774InecZ8+eLfV4hBAFXDoJi+6Fr8aBRwN4dA3cMcuu4Q1yBo7Wmrvuuos///nPfPfdd2RmZjJx4kReeukl3nnnHYuPSU9PZ+LEiezatQsfHx9u3LjBqVOnLB4bEBDAl19+yaBBgwBYvHgxQUFB9no7FvXt25fly5fnu2/SpEm5n8+fP58uXbrQrFmzCh2XEE4v4wZs+wC2vAsubjDknxD2J3CtmGit9mfgGzZswN3dnUceeQQwmzu8//77fPrppyQnJ1t8zLVr18jIyMDLy/x0rVWrFh06dLB4bN++fdm1axfp6ekkJSURExNDcHBw7tfXr19PSEgIAQEBPProo9y4cQOA1atX4+/vT58+fVi6dGnu8devX+fRRx+le/fuhISE8N1335Xrfb/66qu8++67LFmyhMjISB588EGCg4NJSUkp1/MJUe3EboSPwmHjP6D9UJi8G3r9pcLCGyrhGfgPPxS+z88POnWCjAxYtarw1zt0gPbtITUV1q7N/7Xbby/+9Q4ePEi3bt3y3VevXj18fX2JiYmx+JhGjRoxatQoWrZsyaBBgxg5ciRjx47FxcJFCqUUgwcP5scffyQxMZFRo0Zx8uRJwGzoPH78eNavX0/79u15+OGH+eijj5g0aRJ//OMf2bBhA23btuW+++7Lfb433niDgQMH8umnn3LlyhXCwsJyp0WKsmXLltwfGmPGjOGll17K/do999zDrFmzePfddwkNLVToJYQo6Npv8OPf4MA30KgNjPsG2hb/f9BeSjwDV0p9qpS6oJQ6kOe+RkqptUqp49l/NrTvMO1Ha21xw92i7s8xb9481q9fT1hYGO+++y6PPvpokcfef//9LF68mMWLFzN27Njc+48ePUrr1q1p3749AH/4wx/46aefOHLkCK1bt6Zdu3YopRg3blzuY9asWcOMGTMIDg6mf//+pKamcvr06WLfY9++fYmOjiY6OjpfeAshyiAzA3bMgQ9D4fBy6P8i/Hm7w8IbSncGPh+YBXyR576pwHqt9Qyl1NTs2y/YYkDFnTG7uRX/dXf3ks+4C+rcuTPffPNNvvuuXr3KmTNn8PPzIyEhocjHBgQEEBAQwEMPPUTr1q2ZP3++xePCwsI4cOAAHh4euWENv/cct6SoHx5aa7755ptCUzbnz58v8rmEEFaKi4Tlf4Xf9oPfIBj+DniVXCFubyWegWutfwIuFbj7DuDz7M8/B0bbeFwVZtCgQSQnJ+euyMjMzOSZZ55h/PjxeHp6WnxMUlISmzZtyr0dHR1Ny5Yti32dN998k3/+85/57vP39+fUqVO5UzX//e9/6devH/7+/pw8eZLY2FgAvvzyy9zHDBkyhA8//DA3/KOiosr2hi0o2BpXCJEt+RL88BTMG2wqKcd8bqZMyhDeGRmmkt4eynsRs4nW+hxA9p+NizpQKTVRKRWplIq8ePFiOV/OfpRSLFu2jP/973+0a9eO9u3b4+7uni9s169fj4+PT+5HVFQUb7/9Nh06dCA4OJjp06cXefadY9iwYQwYMCDffe7u7nz22WeMGTOGgIAAXFxcmDRpEu7u7sydO5cRI0bQp0+ffD8cpk2bRnp6OoGBgXTp0oVp06ZZ/T0YP348kyZNkouYQuTQGqIWwqxQ2Ptfc3Fy8i7oPBqKmVotKDkZvv4aDh60zzBL1U5WKdUKWK617pJ9+4rWukGer1/WWpc4D17edrKiapC/a+EUzh+CFVPg9HZo0QNGzISmXcr0FKmpZkoXYOtWaNsWmjYt/5CKaidb3lUo55VSN2utzymlbgYulH9oQghRCdxIgs0zYPt/wL0+jJoFwQ+WqQQ+JQV274bYWLj3XqhdG/r0sd+Qyxvg3wN/AGZk/1m+xchCCOFoWsPhH2D1VLj6K3R92JTAezYq9VNkZcGBA7BnD2RmQpcuUKOGHcecrcQAV0p9CfQHvJVSccB0THB/rZR6DDgNjLHnIIUQwi4unYRVz8PxNdAkAMbMhxZhZXqKjAxYuhSuXAFfX+jZExo0KPlxtlBigGutxxbxpUE2HosQQlSMjBuw7d+w5b3sEvg3IWximaook5PB09Msb27bFm66CVq0sOOYLah0lZhCCGFXsRthxTNwKRY632n6l9QrfR+gGzfMVMmhQ3DHHSa4u3a143iLIQEuhKgerp4zJfAHl2aXwC+FtqWfSMjKgsOHITIS0tKgY0eoW9eO4y2Fat/MKiEhIbfVatOmTWnevHnu7ZyWsrawbt263Lau/v7+TJ06tVSPGT26+BqpvXv3snr16tzby5YtK7KLYnlkZGTQoIgJvdmzZ7Nw4UKbvZYQdpGZATs+glnd4cgK6P+37BL40oe31vD997BtG3h5wd13m9UlOUsFHaXan4F7eXkRHR0NmA59derU4dlnn813jNYarbXFZlVlMWDAAL799luSk5MJCgrizjvvpEePHlY95969ezlw4ABDhw4F4M4777Tq+criL3/5S4W9lhDlcmY3rPgr/PZzuUrgk5KgTh1Tu9OhAwQHQ6tW9htuWVX7M/CixMTE0KVLFyZNmkTXrl05c+ZMvjPRxYsXM2HCBMD0IbnrrrsIDQ0lLCyMHTt2FPvcnp6eBAUF8euvvwKmNH/8+PGEhYUREhLCDxZaMu7YsYNevXoREhJC7969OX78OCkpKbz++ussXLiQ4OBglixZwrx583j66acBOHnyJAMGDCAwMJBbb72VuLg4AMaNG8dTTz1FeHg4bdq0YdmyZQD8+uuv9OnTh+DgYLp06UJERETu60+dOpWgoCB69erFhQtm2f/LL7/Mv/71LwD69OnD008/Ta9evQgICKBgwZYQFSqnBP7/boXrCXDvF2UqgU9Lg507YfFiyG4eSseOlSu8obKdga+aan5S2lLTABg2o1wPPXToEJ999hlz5swhIyOjyOOefPJJnn/+eXr27MmpU6cYOXIkBw4cKPL4S5cuceLECfpkr/B//fXXGTp0KPPnz+fy5cv06NGDW2+9Nd9jOnbsyNatW3F1dWX16tW8/PLLfPXVV7zyyiscOHAgN0jnzZuX+5jHH3+cCRMm8OCDDzJ37lyefvrp3J2ALly4wLZt2/j555+59957ufPOO1mwYAG33347L7zwApmZmbll9YmJifTr148ZM2YwZcoUPv30U4tTQDdu3GD79u1s2LCBCRMm5P5mI0SFycqCfYtg7SuQcsWUwPefavamLAWt4dgx2LXLFOW0bw9Nmth5zFaoXAFeyfj5+dG9e/cSj1u3bh1Hjx7NvX358mVSUlLw8PDId9zGjRsJDAzkyJEjTJs2jcaNTQuZNWvWsGrVKmbMMD9oLLWIvXLlCg8//HBug6vS2LlzZ+5OPA8//HC+vimjR49GKUVgYGDubwLdu3fnT3/6E6mpqYwePZqgoCAyMjLw8PBg2LBhAHTr1o0tW7ZYfL2cVrkDBw7kwoULJCUlUadOnVKPVwirnD9oVpdYUQK/di2cOmVCe+hQs8KkMqtcAV7OM2V7qV27du7nLi4u+dq/pqam5n6utWbXrl3UrFmz2OfLmQM/cuQIffv2ZfTo0QQEBKC15ttvv8XPL/+vd3lD/KWXXmLIkCE8/vjjxMTE5M55l1etWrXyjR9M8G7atIkVK1bw4IMP8uKLL3Lffffle1+urq5F/jZSsAVucf3UhbAZK0vgk5LAwwNcXc0Zd5s2Zl23M5A58FJycXGhYcOGHD9+nKysrNx5Y4DBgwcze/bs3NslTR34+/vz/PPP8/bbbwOmRewHH3yQ+3VLLWITExNp3rw5QL7Oh8W1gu3Zsydff/01AAsWLOCWW24pdly//PILTZs2ZeLEiYwfP77MrWq/+uorwGyk3KRJk3w/AIWwOa3h0HcwOwwiPoSQB+GJPdD1oVKFd0aGWc+dt1tgq1bOE94gAV4mb731FkOHDmXQoEH4+Pjk3j979my2bdtGYGAgnTp14pNPPinxuR5//HHWr1/P6dOnmT59OsnJyQQEBNC5c2deffXVQse/8MILPPfcc/Tu3Tvf/QMHDmTfvn2EhIQU2ul+1qxZzJ07l8DAQL766ivef//9Yse0fv16goKCcvfafOKJJ0p8H3nVq1eP8PBwnnjiiVJ9D4Qot0snYOEY+Pph8GgEj62FUR+Wun9JTIwJ7j17oGVLc9btjErVTtZWpJ1s1dWnTx9mzZqVb8PmguTvWlitYAn8gJfKXAK/daupovT2hvBw69q8VhRbt5MVQoiKFbsBVjxbrhL45GQzx12rlpnn9vY267qd/TKNBLiwia1btzp6CKKqsqIEPjPTtHndu9cEd+/e0Lix+agKKkWAl7QDvHB+FTlVJyqXOZtjCfSpT7ifd5HHRMTGsz8ukUn98qzEysyA3Z/AhjcgM82UwPd+CmqUrn791CnYsQOuXjXz3F3KtqLQKTj8Iqa7uzsJCQnyH7wK01qTkJCAu6MbRwiHCPSpz+RFUUTExlv8ekRsPJMXRRHoU//3O8/shk/6m00WWoTB49uh/wulDu+9e2HNGjNtMnw4DBkC9euX/Dhn4/AzcB8fH+Li4qiMGx4L23F3d8+3ckdUH+F+3sx6IITJi6KY9UBIvjPxnPDOvT/5Eqx7FfZ+DnWbmRL4jqNKNVmdmmqWBtapY5YC1qwJnTqVaUc0p+PwVShCiOqhYFjnu926Uf4S+J5/LnUJfFaWWVWyZ8/vFZRVjaxCEUI4VN4z8XE9fFmw87QJ7zrn4bMH4cyOMpfAx8VBRITZzqx5c7CyuafTkQAXQlSYcD9vxvXw5YMNMTzT72bCY943vbrd68MdsyHogVLPeRw9Cps3Q716Zo67ZUs7D74SkgAXQlSYiNh4Fuz4hQ+DThO2czJwCbr+AQa/WqoqyrQ0s6a7QQNo3drc7tTJXKysjiTAhRAVIiI2nhkLV7G+6f9oeHQz1xt2ZHziM0zsfD/hJYS31uaMe9cus5HwPfeYi5QBARU0+EpKAlwIYXfbj/3Kvi9f41u1DJf4mjDkTWqHTWTiqSsWV6fkde6cmedOSDBl7+HhFTz4SkwCXAhhVwe3fEuz9VPpxTnodBcMeSO3BL64JYYAZ87AqlVmaeCgQeBX+t3QqgUJcCGEfVw9Bz++SOeDy0ip2xLusFwCnxPi++MSCffzJiMDLl82myk0b27OuP39wU3SqhD5lgghbCszA3bNhY3/zC2B9yihBD7cz5twP29iYsxelFlZMHasCe2qWAJvKxLgQgjbObMLlk+B8z9D28FmF/hGJTfbvnDBzHNfuPB7m1c54y6ZVd8ipdRfgQmABn4GHtFapxb/KCFElZN8CdZNh71flLkE/tIl+PZbs61Zv36ma6D0tiudcge4Uqo58CTQSWudopT6GrgfmG+jsQkhKrusLIheaErgUxOh1+RSlcBnZsL589CsGTRqBLfcYnbFKWFbWVGAtb+kuAEeSql0wBM4a/2QhBBO4bcDZhf4MzugRU8YOROadC7xYSdPmjav16/DAw+Ydd3+/hUw3iqo3AGutf5VKfUucBpIAdZordfYbGRCiMrpxjXYNKPMJfCXLpl57rNnzVn3sGEmvEX5WTOF0hC4A2gNXAH+p5Qap7VeUOC4icBEAF9fXyuGKoRwqJxd4Fe/CNfOlqkEPiUFli6FGjXMrjgdO1btNq8VxZoplMHASa31RQCl1FIgHMgX4FrrucBcMO1krXg9IYSjJMTCquchZh00CTAXKVt0L/YhWVmmW6Cvr7lAOXCgWdddq1YFjbkasCbATwM9lVKemCmUQYA0+xaiKklPhW3/gi0zwbUmDJ0B3f9Y4i7wZ87A9u2mzevdd4OXl7lIKWzLmjnwnUqpJcBeIAOIIvtMWwhRBcSsh5XPwqUT0Pmu7F3gby72IVeumAuUp0+bNq9Dh5rwFvZh1SoUrfV0YLqNxiKEqAyuns3eBX4ZNPKDh5aB38ASH5aRAd9/b6ZOevY0FZQyz21fUuskhDAyM2DXx9kl8Okw4CUIf7LYEnitzbLA1q1N5eTAgeaM28OjAsddjUmACyHg9E5YMQXOH4C2t8Lwt0ssgT971sxzJySYqRJfX5B9qyuWBLgQ1VneEvh6zeHe/0LH24utZb92zcxznzxp2rwOHmzCW1Q8CXAhqqOCJfDhT0C/qVCrTokPXb3ahHhoKAQGStMpR5JvvRDVzW8HzHTJmZ3g2wtGvFdsCbzWcOKE2TTYzc00nKpd23wIx5IAF6K6yFsC79EA7vgPBI0tdqlI3javffuaCsrGjStwzKJYEuBCVHUFS+C7jYdB04stgb9+3WwgfPy46VfSvz+0a1dhIxalJAEuRFWWEAsrn4PY9dC0dCXwAJs3m82EQ0IgONj0MBGVjwS4EFVRoRL4t6D7hGJL4E+cMLu+e3qaHXFcXaFu8W29hYNJgAtR1eQtge9yN9z2RrEl8AkJZp773Dno2tWsLmnQoALHK8pNAlyIquLqWTPPfejbUpXAp6RAZCQcPgzu7uYipWys4FwkwIVwdnlL4LMyYMDL0PtJcCu+b+vu3XDsGAQEQLdusp2ZM5IAF8KZFSqBfwcatS768NOmerJRo98LcWS6xHlJgAvhjJIvmSrKqP+WqgT+yhXTt+TMGbPre//+5mKlbGnm3CTAhXAmWVkQvQDWTocbV023wH4vFFkCf+MG7NkDBw+apYC9ekHnkvcdFk5CAlwIZ1GoBH4mNOlU7EMOHjQfHTuaKRP3ojvDCickAS5EZXfjGmx8E3bOMSXwoz8yJfBFTJecPWv+bNbMzHG3bCm74lRVEuBCVFZamyWBq1+Ea79ll8C/UmQJ/NWrps3rqVOmL3ezZqb5lIR31SUBLkRllK8EPhDuWwA+oRYPTU+HqCjYv9/0pQoLM0sDRdUnAS5EZZKeClvfNx9utWDY2xD6WLEl8CdPQnS0WV0SFiYrS6oTCXAhKouYdbDiWbh8ErrcA0PegLpNLR56/rzpGNimjekS2KgReHtX8HiFw0mAC+Foib/Cjy+alq9ebeHh76BNf4uHXr8OO3dCTAw0bGg2E1ZKwru6kgAXwlEy02Hnx7DpTVMCP/Bls67bQgl8RoaZ446ONtc2u3aFoKBit64U1YAEuBCOcHoHLJ8CFw5Cu9vMXHcxJfAXL5rGU23aQI8e0uZVGBLgQlSk6wmw7hWIWgD1fOC+heA/wuKpdHy82cqsUye4+Wa45x4z1y1EDglwISpCVpbpW7JuuinM6f0U3PK8xRL4lBTTKfDIEbOipH17s55bwlsUZFWAK6UaAPOALoAGHtVab7fFwISoMn772UyXxO0C3/DsXeALl8BnZsKBA7B3r/k8MNDMdbvJaZYogrX/NP4NrNZa36OUqgnIClQhcty4Znp075wDHo1g9BwIur/IK4/Xr5szbx8f03Sqfv0KHq9wOuUOcKVUPeAWYDyA1joNSLPNsIRwYlrDwWXw499MCXzoIzBwmsUS+MuXTSFO165Qrx6MGSPBLUrPmjPwNsBF4DOlVBCwB3hKa33dJiMTwhklxJr9KGM3FFsCf+OGWVVy6JBp8+rvb+a7JbxFWVgT4G5AV+AJrfVOpdS/ganAtLwHKaUmAhMBfH19rXg5ISqxQiXw70D3x8DFNd9hWVlmD8rISEhLkzavwjrWBHgcEKe13pl9ewkmwPPRWs8F5gKEhoZqK15PiMrp+Dpz1n35JASMgdv+UWQJfHq62WDBywvCw2VlibBOuQNca/2bUuqMUqqD1vooMAg4ZLuhCVHJ5SuBb1dkCfzVq2ZThZ49oVYtuOsusy+lENaydhXKE8DC7BUoJ4BHrB+SEJVcKUvg09JMm9effwZXV7Oe28tLwlvYjlUBrrWOBiw3KRaiKspXAj8Ehr8NDVvlO0RrOHYMdu0yRTnS5lXYi5QICFEaZSiB1xr27TPLAocOhZtucsB4RbUgAS5EcSyVwPd7AWrWzndYUpKZLunZ0ywLHDlSzriF/UmAC1GUvCXwLXubEvjGHfMdkpFhWrzu22dut2kDzZtLeIuKIQEuREGpV80FyhJK4GNizOYK16+Dn59p8yoXKEVFkgAXIoelEvhBr4BHQ4uHHz0KHh4waBA0tbzsWwi7kgAXAvKXwN8cZC5S+nTLd0hysqmg7NrVnGkPGmTWdcuuOMJRJMBF9ZaekqcE3t1iCXxmplnLHRVlPm/WDNq2lfJ34XgS4KL6KlQC/wbUbZLvkF9+ge3bTTVly5ZmlYk1DafmbI4l0Kc+4X5F70IcERvP/rhEJvXzK/8LiWrBxdEDEKLCJf4KXz0EC+8GFzd4+Hu4e16h8AYT4K6uMHw4DBlifbfAQJ/6TF4URURsvMWvR8TGM3lRFIE+0pZQlEzOwEX1kZluVpZsfBN0punRHf5EvhL41FQzz92+PTRubDZWcHUFFxud6oT7eTPrgRAmL4pi1gMh+c7Ec8K74P1CFEUCXFQPeUvg2w+FYW/lK4HPyjK9uffsMT1MGjQwAV6jhu2HYinEJbxFeUiAi6rtegKsfQWiF0D9FnD/IugwPN/Skbg4iIiAK1d+386soeWVgzaTN8TH9fBlwc7TEt6izCTARdWUlQVRX8C6V7NL4J+Gfs8XKoEHiI83hw8ZYi5UVpRwP2/G9fDlgw0xPDmwrYS3KDO5iCkcZs7m2CIv5uWIiI1nzubYsj3xuf3w6W3ww1PQuDNM2ga3vpYb3mlpsGMHnDhhDg8IMHtRVmR4g3lvC3ae5smBbVmw83SJ3wshCpIAFw5j8xUZqVdh1VSY2w8unYQ7P4bxy6GxP2AKLY8cgcWLYf9+SEgwD3N1NR8VKe+c95TbOuROp0iIi7KQABcOk3ceuGBwlemintZw4BuY1d2sMun2CDwRma9/yfnzsHQp/PSTuUB5113Qvbu93lnxLL234r4XQhRFAlw4lKXgKlN4x8fAf++EJY+afSj/uB5GzizUvyQpyewEP2gQjBoF3g6abi7uvUmIi7JSWlfcPsOhoaE6MjKywl5POI+cYCv1ioz0FNgyE7b9y5TAD3oFQh/NLYHPafNaq5aZ4865z83Bl+2lElOUh1Jqj9a60O5nsgpFVAplWpFxfG12CfwpCLg3exd4U0WptWnzumuXafPaMU/7bkeHN1CqUA7385YVKaJUKsE/aSEKr8jo6edVOMQS42D1i3D4e/BuD3/4AVrfkvvlhATYsgUuXDDbmA0eDE0KV8cLUWVIgAuHKzgv3NPPK/88cWY67PgINs0AnWWmS3o9AW418z1PZqaZ6+7fH9q1kzavooWldD8AABA3SURBVOqTOXDhUEVd1Mu5//PBmQREvQYXDhUqgc/MNMsBb9wwXQJz7qvoJYFC2JvMgYtKp9gVGU1hVavFNPlxCTdqN6PW/YvMLvDZTp40xTjXrkHr1mbuWykJb1G9SIALh9kfl1g4vLOyYO/nsO5VmqQlEdf5T6z2epgJ/l0A05f7p5/g7Flo1AhGjDCbCAtRHUmAC4cptCLj3H5Y/lf4NRJa9oER7+HT2J8JeQ5xcYHEROjTB/z9bdfmVQhnJAEuHC/1Kmz8J+z62OwCf+fHEHgfKEVWFhw8aCopBw82e1GOHSvBLQRIgAtH0hoOLoXVf4Ok86YQZ9C03CrKM2fMdmY5bV7T001/bglvIQyrA1wp5QpEAr9qrUdaPyRRLcTHwMpn4MQmuDkYxi6C5mYX+ORk2LzZBHj9+jB0KPj6Ona4QlRGtjgDfwo4DNSzwXOJqq5gCfzwd/OVwIM5y752zSwN7NJFzriFKIpVAa6U8gFGAG8AU2wyIlF1FVECrzUcOWxK4EeMMAE+ZowU4ghREmvPwP8FPA/ULeoApdREYCKAr/weXD0lxsHqqXD4h0Il8GfPmnnuhARo2tRsKuzpKeEtRGmUO8CVUiOBC1rrPUqp/kUdp7WeC8wFU4lZ3tcTTqiYEvi0NDPPffKkWVkyeDC0aePoAQvhXKw5A+8NjFJKDQfcgXpKqQVa63G2GZpwar9shxVTLJbAg5kmSU6G0FAIDKwcnQKFcDbl/m+jtX4ReBEg+wz8WQlvwfX47F3gF/6+C7z/CLSG48dM75Lbbzd9ukeNkqkSIawh5z3CNvKUwJOWBH3+Crc8BzVrc/68mee+cAEaNzbNp2rVkvAWwlo2CXCt9SZgky2eSzihc/tg+ZR8JfA09icrCzZtMKtLPD1hwABo21aCWwhbkTNwUX6pV2HjG7BrLnh65ZbAaxQKs347KwtCQiA42Mx7CyFsRwJclF3OLvA/vmRK4Ls/BgNfBo+GnDgBu3fDsGFQr55ZXSKEsA8JcFE2RZTAx8dDxFr47TfT5jUtzdEDFaLqkwAXpZOeAlveg23/BjePfCXwW7fCoUPg7g59+5o2rzLPLYT9SYCLkh1bY0rgr/xi2rze+nd0nSa5Ie3mBgEB0K0b1KxZ/FMJIWxHAlwULTEOVr0AR5bnK4E/fRq2r4RbboGbb/59P0ohRMWSABeFZabDjv/AprfylcBfSapJxEqIi4MGDRw9SCGEBLjI75cIs6b74mFoPyy7BL4lkZEQFWWWAvbqBZ07S5tXIRxNAlwY1+NhzTTYtyi7BP5LstoPRylQgIcHdOxoepe4uzt6sEIIkAAXWVmwdz6sey1fCfyvF2uzfalpNNW+vTnjFkJULhLg1dm5fdm7wO+BVn1h+LtcdfdnxyY4dQrq1pWzbSEqMwnw6ig1ETa8Abs/yS6BnwuB93LgoGLHDjO3HRZmlga6upb8dEIIx5AAr05yS+D/BkkXoPtj6IHT0LUa4KLMGXfbtia8PT0dPVghREkkwKuL+OOw4hk4uTm7BH4x5926ErEaWrQwFydbtjQfQgjnIAFe1VkogU/yf5Rdka65bV4bNnT0IIUQ5SEBXpUd+xFWPvd7Cfxt/+D4ucZsWWJmU7p2haAgafMqhLOSAK+Krpwxu8AfWQ7eHeAPy8lo0Rc3N9Mp0NcXevQwc95CCOclAV6V5JbAzzCn2IOmE99+MhE7a1L3lNkRx8tLenQLUVVIgFcVp7aZi5QXD0OH4aQOmMGuYy058p1Zy92+vaMHKISwNQlwZ5d00ewCv28R1PeF+7/kTO3hrFsFmZmmkrJrV2nzKkRVJAHurPKVwF+HPlNI6/ksNevUxivZLA3s3h3q13f0QIUQ9iIB7ozORsOKKbkl8Il932PrkQ5kboBRo8zSQJnnFqLqkwB3JgVK4NNv/4RdKWM4tEFRo4YpxtFatjMTorqQAHcGhUrgJxAf9DIr1jcgLQ06dTLbmUnjKSGqFwnwyi5vCXyzENLuWUzNVl1pkGHmuYODzdpuIUT1IwFeWaWnwE/vmhL4Gp6kDHyPn5If4fIOV8a0MBsJDxzo6EEKIRyp3AGulGoBfAE0BbKAuVrrf9tqYNVanhL4zC73E93470Qda4yrK4SEyBy3EMKw5gw8A3hGa71XKVUX2KOUWqu1PmSjsVU/BUrgr49ZztK9fUk5Ah06mGWB0uZVCJGj3AGutT4HnMv+/JpS6jDQHJAAL6vMdNg+Gza/BVqT3u9VavT9C56uNWl92YT3TTc5epBCiMrGJnPgSqlWQAiw08LXJgITAXx9fW3xclXLqW1mTffFI2S0HcHOhjM4/osv9/cCdzfo08fRAxRCVFZWB7hSqg7wDfC01vpqwa9rrecCcwFCQ0O1ta9XZSRdhLXTYN+X6Pq+xHRfzE+/DYPfzMoSN7m8LIQogVUxoZSqgQnvhVrrpbYZUhWXlQl75sP61yAtmYxeU1hy4Tmu/uqJn59p81qnjqMHKYRwBtasQlHA/wGHtdYzbTekKixPCXymb19cb38Pt5s60DYSfHygaVNHD1AI4UysOQPvDTwE/KyUis6+729a65XWD6uKSU2EDf+A3fPQnt4c7vgJEVfHcLeboiGmBF4IIcrKmlUoWwFZkVwcreHnJfDj39DXLxLfegJrMl4m5VoDAgKhdm1HD1AI4czkUpm9XDwGK5+Bkz+hm4WwrvHXnLwRQqtW0LMn1Kvn6AEKIZydBLitpSXDlndh2wfoGp6oEe+huj1Cs8OudKxv5rqFEMIWXBw9gCrl6Gr4Tw/Y8h4XGt/NwoaRnGk6AVxc6dy5+PCeszmWiNj4Yp8+IjaeOZtjbTxoIYSzkgC3hStnYPGD8OV9pGZ6sKrpCr5z/ZjWAY1LXUEZ6FOfyYuiigzxiNh4Ji+KItBHttgRQhgS4NbISIOt78PsMIjdwKHmr7Kgzla0bx/uuQd69y59j+5wP29mPRBiMcRzwnvWAyGE+3nb4Y0IIZyRzIGX16mt2bvAH0F3GIEaNgP3S77c6gotW5bvKfOGeE5YS3gLIYoiAV5WeUrgUz18+clrMT5Bw+jUANo0sP7p84b4uB6+LNh5WsJbCGGRBHhpZWXCns/Q61+HtGR+bvAMu92fpV0nT1q3tu1Lhft5M66HLx9siOHJgW0lvIUQFkmAl8bZKFg+Bc7u5Ur9vqypOxOPFu25Ixy87ZCtEbHxLNh5micHtmXBztP09POSEBdCFCIBXpyUK7DxDfTueeDpjbprHqne99A9VdGmjX1esuCcd08/L5kDF0JYJKtQLNEa9n+NntUdvWseh2pPYG/4bggcw83NKi68ofjVKUKI6k0CvKCLx9Cf3w5L/0hChg/LvDdwvus7dAiywRXKYhS32kRCXAhhiUyh5MhTAp/p4sn2+jOJbzme3r1dadLE/i+/Py6x2GmSnBDfH5coUylCCACU1hW3SU5oaKiOjIyssNcrtaOryVr5HC6JpyFoLIlhr3P+emPatZMd4IUQjqeU2qO1LtR4unqfgV85jV45FXVsBYk1/DkRsIJud/ahPiAF60KIyq56BnhGGuyYTdamt83y7rqvcbXz4/QIr+nokQkhRKlVvwDPUwL/i/tIDrd5k+D+vjRr5uiBCSFE2VSfAE+6QOaqabgeXAwNfEm75ytSXYcytAO4yFocIYQTqvoBnpVJVuRnZK19HZWezNEmz9LhsWeoWdOTjo4emxBCWKFSnnvabHODs1Hc+M9gXFY+w3kVxNaACJqMnQY1PW04WiGEcIxKGeBWb26QcgVWPIueO4DMS3FsbzaPzAe/p9/d7Wlg33ocIYSoMJVyCsVSX+wcxfbH1pr0vf/Ddd1LuKTGo7tP5HTzl+gRWF/muYUQVU6ljTVL5ePFhbe+cJTrH91OjR/+yOUsH7Ie24DLiLfxD5bwFkJUTZXyDDxHqTY3SEvm2sp3qL3vQ9zwZL/vTJrfPh6Xm1wdN3AhhKgAlTrAoYTNDY6uIvOH56mbdJoT9R7AZcjrBHYu5S7CQgjh5Cp9gFva3KB7/WQyfngBj19W4nqTP2d6rsS3Z2/cKv27EUII27Eq8pRSQ4F/A67APK31DJuMKluhzQ1a1mXHx38n1HMxbsqFjAGv4dbnL7RwrWHLlxVCCKdQ7gBXSrkCs4FbgThgt1Lqe631IVsMrGB4X9q7hS4/PkO451F+vNYXt1H/ZFC/QFu8lBBCOCVr1meEATFa6xNa6zRgMXCHLQZVMLzTlj5Fo+9H4pKZwtl+X1Hnz/N57qfzsrmBEKJas2YKpTlwJs/tOKBHwYOUUhOBiQC+vr6leuKCmxvUbNyaKwHPUnvYMzTz9KQZyOYGQohqr9wbOiilxgBDtNYTsm8/BIRprZ8o6jGVdkMHIYSoxIra0MGaKZQ4oEWe2z7AWSueTwghRBlYE+C7gXZKqdZKqZrA/cD3thmWEEKIkpR7DlxrnaGUmgz8iFlG+KnW+qDNRiaEEKJYVq0D11qvBFbaaCxCCCHKQNo8CSGEk5IAF0IIJyUBLoQQTkoCXAghnFS5C3nK9WJKXQR+KefDvYHqVjsv77l6kPdc9Vn7fltqrQv1yq7QALeGUirSUiVSVSbvuXqQ91z12ev9yhSKEEI4KQlwIYRwUs4U4HMdPQAHkPdcPch7rvrs8n6dZg5cCCFEfs50Bi6EECIPCXAhhHBSThHgSqmhSqmjSqkYpdRUR4/H3pRSnyqlLiilDjh6LBVBKdVCKbVRKXVYKXVQKfWUo8dkb0opd6XULqXUvuz3/Jqjx1RRlFKuSqkopdRyR4+lIiilTimlflZKRSulbLqjTaWfA8/ePPkYeTZPBsbaavPkykgpdQuQBHyhte7i6PHYm1LqZuBmrfVepVRdYA8wuor/HSugttY6SSlVA9gKPKW13uHgodmdUmoKEArU01qPdPR47E0pdQoI1VrbvHDJGc7A7bZ5cmWltf4JuOTocVQUrfU5rfXe7M+vAYcxe65WWdpIyr5ZI/ujcp9N2YBSygcYAcxz9FiqAmcIcEubJ1fp/9zVmVKqFRAC7HTsSOwveyohGrgArNVaV/n3DPwLeB7IcvRAKpAG1iil9mRv8m4zzhDgysJ9Vf5MpTpSStUBvgGe1lpfdfR47E1rnam1DsbsJxumlKrS02VKqZHABa31HkePpYL11lp3BYYBf8meIrUJZwhw2Ty5GsieB/4GWKi1Xuro8VQkrfUVYBMw1MFDsbfewKjsOeHFwECl1ALHDsn+tNZns/+8ACzDTAvbhDMEuGyeXMVlX9D7P+Cw1nqmo8dTEZRSNymlGmR/7gEMBo44dlT2pbV+UWvto7Vuhfl/vEFrPc7Bw7IrpVTt7AvzKKVqA7cBNltdVukDXGudAeRsnnwY+Lqqb56slPoS2A50UErFKaUec/SY7Kw38BDmjCw6+2O4owdlZzcDG5VS+zEnKWu11tViWV010wTYqpTaB+wCVmitV9vqySv9MkIhhBCWVfozcCGEEJZJgAshhJOSABdCCCclAS6EEE5KAlwIIZyUBLgQQjgpCXAhhHBS/w8yOGXBfvqJQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the raw data, with the true underlying relationship, and the OLS fit\n",
    "fig, ax = plt.subplots(1)\n",
    "ax.plot(x, y, 'x', markersize=10, label='Data')\n",
    "ax.plot(x, 2*x, '--b', alpha=0.4, label='OLS Model Fit')\n",
    "ax.plot(x, theta*x, label='True Relationship')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predict New Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The prediction for a new x of 2.5 is 5.299\n"
     ]
    }
   ],
   "source": [
    "# With our model, we can predict the value of a new 'x' datapoint\n",
    "new_x = 2.5\n",
    "pred_y = theta * new_x\n",
    "print('The prediction for a new x of {} is {:1.3f}'.format(new_x, pred_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3iUVdrH8e9JAiShkwgIIZRQQkmDECCAVKWKWFBRdFFZllUsiw0XWdRdV2y4q7Ai8iq6gOgiWGgCoQiEDgHpJJQQQSABQklC2nn/OElMmdSZyWSS+3NduczMPDNzBuSXJ+c5932U1hohhBDOx8XRAxBCCFE2EuBCCOGkJMCFEMJJSYALIYSTkgAXQggn5Vaeb+bt7a1btGhRnm8phBBOb/fu3fFa61vy31+uAd6iRQt27dpVnm8phBBOTyl12tL9MoUihBBOSgJcCCGclAS4EEI4qXKdA7ckLS2NuLg4UlJSHD0UYUfu7u74+PhQrVo1Rw9FCLuavTGGQJ+6hPt5F3pMZEw8++MSmdDHz6r3cniAx8XFUbt2bVq0aIFSytHDEXagtSYhIYG4uDhatmzp6OEIYVeBPnWZuHAvMx8KsRjikTHxOY9by+FTKCkpKXh5eUl4V2JKKby8vOS3LFElhPt5M/OhECYu3EtkTHyex3KHd1Fn6CXl8AAHJLyrAPk7FlWJpRC3dXhDBQnw4szeGFPgJ1l+kTHxzN4YU04jEkKIouUO8Rmrj9o8vMFJAjx7TqmwEM/+yRboU7dMr+/q6kpwcDAdO3YkKCiIGTNmkJmZWeRzTp06xcKFC8v0fkKIqiHcz5sx3Xz5cF00Y7r52jS8wUkC3N5zSh4eHkRFRXHw4EHWrFnDihUreP3114t8jgS4EKI4kTHxzN8eyzP9WzN/e2yxMwmlprUut68uXbro/A4dOlTgvsJsib6oQ95YrbdEX7R4u6xq1qyZ53ZMTIxu0KCBzszM1CdPntS9evXSISEhOiQkRG/ZskVrrXW3bt10nTp1dFBQkJ4xY0ahx4nflebvWghnZ8u8AnZpC5nqVAGu9e9/CO//dMQm4a11wQDXWut69erp3377Td+4cUMnJydrrbU+duyYzv4M69ev18OGDcs5vrDjxO8kwEVVUVhYlzXECwtwh68DL63cc0rP9G9t8zmlbDprr9C0tDQmTpxIVFQUrq6uHDt2zOLxJT1OCFG5FTWtm3s62BYXNJ0uwPPPKXX387J5iJ84cQJXV1caNmzI66+/TqNGjdi3bx+ZmZm4u7tbfM4HH3xQouOEEJXb/rjEIsM5O8T3xyVWrQDP/5Otu5+XzZfmXLx4kQkTJjBx4kSUUiQmJuLj44OLiwtffPEFGRkZANSuXZtr167lPK+w44QQVUtJyuPD/bwrTyFPSVj6taSo1SmlkZycnLOMcODAgdxxxx1MmzYNgCeffJIvvviC7t27c+zYMWrWrAlAYGAgbm5uBAUF8cEHHxR6nBBC2IvKnustD6GhoTr/hg6HDx+mffv2RT6vuKWC9qhwErZXkr9rIURBSqndWuvQ/Pc7xRl4aeaUhBCiqnCKOfDynFMSQgiby8wEF9ufLzvFGbgQQjglrSFqIcwKgxs2rsJEAlwIIezj/CH4fCh892fwqA83r9r8LZxiCkUIIZzGzeuw8W3Y9h+oURtGfATBY+wyhSIBLoQQtqA1HFkGKyfD1TgIeQQGvg41vez2ljKFgtnW7a677qJNmzb4+fnx7LPPkpqaCsCGDRsYPnx4gecsW7aMkJAQgoKC6NChA5988kmBY+bNm4dSioiIiJz7li5dilKKxYsXl3h8hY2hJMds2LCBunXrEhwcTHBwMAMHDgRg9uzZfPnllznjPHv2bInHI4TI59JJWHg/fD0GPOrB46vhrpl2DW+QM3C01txzzz38+c9/5vvvvycjI4Px48czZcoU3n33XYvPSUtLY/z48ezYsQMfHx9u3rzJqVOnLB4bEBDAV199xYABAwBYtGgRQUFB9vo4FvXu3Ztly5bluW/ChAk538+bN49OnTrRpEmTch2XEE4v/SZs+RA2vQcubjDonxD2J3Atn2it8mfg69atw93dncceewwwmzt88MEHfPbZZyQlJVl8zrVr10hPT8fLy/x0rVGjBu3atbN4bO/evdmxYwdpaWlcv36d6OhogoODcx6PiIggJCSEgIAAHn/8cW7evAnAqlWr8Pf3p1evXixZsiTn+Bs3bvD444/TtWtXQkJC+P7778v0uV977TXee+89Fi9ezK5du3j44YcJDg4mOTm5TK8nRJUTsx4+Dof1/4C2g2HiTujxVLmFN1TAM/Affyx4n58fdOgA6emwcmXBx9u1g7ZtISUF1qzJ+9iddxb9fgcPHqRLly557qtTpw6+vr5ER0dbfE6DBg0YMWIEzZs3Z8CAAQwfPpzRo0fjYuEihVKKgQMH8tNPP5GYmMiIESM4efIkYDZ0Hjt2LBEREbRt25ZHH32Ujz/+mAkTJvDHP/6RdevW0bp1ax544IGc13vzzTfp378/n332GVeuXCEsLCxnWqQwmzZtyvmhMWrUKKZMmZLz2H333cfMmTN57733CA0tUOglhMjv2m/w01/hwLfQoBWM+RZaF/1v0F6KPQNXSn2mlLqglDqQ674GSqk1SqnjWf+tb99h2o/W2uKGu4Xdn23u3LlEREQQFhbGe++9x+OPP17osQ8++CCLFi1i0aJFjB49Ouf+o0eP0rJlS9q2bQvAH/7wB37++WeOHDlCy5YtadOmDUopxowZk/Oc1atXM336dIKDg+nbty8pKSnExsYW+Rl79+5NVFQUUVFRecJbCFEKGemwbTZ8FAqHl0HfV+DPWx0W3lCyM/B5wEzgy1z3TQYitNbTlVKTs26/bIsBFXXG7OZW9OPu7sWfcefXsWNHvv322zz3Xb16lTNnzuDn50dCQkKhzw0ICCAgIIBHHnmEli1bMm/ePIvHhYWFceDAATw8PHLCGn7vOW5JYT88tNZ8++23BaZszp8/X+hrCSGsFLcLlv0FftsPfgNg6LvgVXyFuL0Vewautf4ZuJTv7ruAL7K+/wIYaeNxlZsBAwaQlJSUsyIjIyOD559/nrFjx+Lp6WnxOdevX2fDhg05t6OiomjevHmR7/PWW2/xz3/+M899/v7+nDp1Kmeq5r///S99+vTB39+fkydPEhMTA8BXX32V85xBgwbx0Ucf5YT/3r17S/eBLcjfGlcIkSXpEvz4LMwdaCopR31hpkxKEd7p6aaS3h7KehGzkdb6HEDWfxsWdqBSarxSapdSatfFixfL+Hb2o5Ri6dKl/O9//6NNmza0bdsWd3f3PGEbERGBj49PztfevXt55513aNeuHcHBwUybNq3Qs+9sQ4YMoV+/fnnuc3d35/PPP2fUqFEEBATg4uLChAkTcHd3Z86cOQwbNoxevXrl+eEwdepU0tLSCAwMpFOnTkydOtXqP4OxY8cyYcIEuYgpRDatYe8CmBkKe/5rLk5O3AEdR0IRU6v5JSXBN9/AwYP2GWaJ2skqpVoAy7TWnbJuX9Fa18v1+GWtdbHz4GVtJysqB/m7Fk7h/CFYPglit0KzbjBsBjTuVKqXSEkxU7oAmzdD69bQuHHZh1RYO9myrkI5r5S6VWt9Til1K3Ch7EMTQogK4OZ12Dgdtv4H3OvCiJkQ/HCpSuCTk2HnToiJgfvvh5o1oVcv+w25rAH+A/AHYHrWf8u2GFkIIRxNazj8I6yaDFd/hc6PmhJ4zwYlfonMTDhwAHbvhowM6NQJqlWz45izFBvgSqmvgL6At1IqDpiGCe5vlFJPALHAKHsOUggh7OLSSVj5EhxfDY0CYNQ8aBZWqpdIT4clS+DKFfD1he7doV694p9nC8UGuNZ6dCEPDbDxWIQQonyk34Qt/4ZN72eVwL8FYeNLVUWZlASenmZ5c+vWcMst0KyZHcdsQYWrxBRCCLuKWQ/Ln4dLMdDxbtO/pE7J+wDdvGmmSg4dgrvuMsHdubMdx1sECXAhRNVw9ZwpgT+4JKsEfgm0LvlEQmYmHD4Mu3ZBaiq0bw+1a9txvCVQ5ZtZJSQk5LRabdy4MU2bNs25nd1S1hbWrl2b09bV39+fyZMnl+g5I0cWXSO1Z88eVq1alXN76dKlhXZRLIv09HTqFTKhN2vWLBYsWGCz9xLCLjLSYdvHMLMrHFkOff+aVQJf8vDWGn74AbZsAS8vuPdes7oke6mgo1T5M3AvLy+ioqIA06GvVq1avPDCC3mO0VqjtbbYrKo0+vXrx3fffUdSUhJBQUHcfffddOvWzarX3LNnDwcOHGDw4MEA3H333Va9Xmk89dRT5fZeQpTJmZ2w/C/w2y9lKoG/fh1q1TK1O+3aQXAwtGhhv+GWVpU/Ay9MdHQ0nTp1YsKECXTu3JkzZ87kORNdtGgR48aNA0wfknvuuYfQ0FDCwsLYtm1bka/t6elJUFAQv/76K2BK88eOHUtYWBghISH8aKEl47Zt2+jRowchISH07NmT48ePk5yczBtvvMGCBQsIDg5m8eLFzJ07l+eeew6AkydP0q9fPwIDA7n99tuJi4sDYMyYMTz77LOEh4fTqlUrli5dCsCvv/5Kr169CA4OplOnTkRGRua8/+TJkwkKCqJHjx5cuGCW/b/66qv861//AqBXr14899xz9OjRg4CAAPIXbAlRrrJL4P/vdriRAPd/WaoS+NRU2L4dFi2CrOahtG9fscIbKtoZ+MrJ5ielLTUOgCHTy/TUQ4cO8fnnnzN79mzS09MLPe6ZZ57hpZdeonv37pw6dYrhw4dz4MCBQo+/dOkSJ06coFfWCv833niDwYMHM2/ePC5fvky3bt24/fbb8zynffv2bN68GVdXV1atWsWrr77K119/zd/+9jcOHDiQE6Rz587Nec6TTz7JuHHjePjhh5kzZw7PPfdczk5AFy5cYMuWLfzyyy/cf//93H333cyfP58777yTl19+mYyMjJyy+sTERPr06cP06dOZNGkSn332mcUpoJs3b7J161bWrVvHuHHjcn6zEaLcZGbCvoWw5m+QfMWUwPedbPamLAGt4dgx2LHDFOW0bQuNGtl5zFaoWAFewfj5+dG1a9dij1u7di1Hjx7NuX358mWSk5Px8PDIc9z69esJDAzkyJEjTJ06lYYNTQuZ1atXs3LlSqZPNz9oLLWIvXLlCo8++mhOg6uS2L59e85OPI8++mievikjR45EKUVgYGDObwJdu3blT3/6EykpKYwcOZKgoCDS09Px8PBgyJAhAHTp0oVNmzZZfL/sVrn9+/fnwoULXL9+nVq1apV4vEJY5fxBs7rEihL4NWvg1CkT2oMHmxUmFVnFCvAyninbS82aNXO+d3FxydP+NSUlJed7rTU7duygevXqRb5e9hz4kSNH6N27NyNHjiQgIACtNd999x1+fnl/vcsd4lOmTGHQoEE8+eSTREdH58x5l1WNGjXyjB9M8G7YsIHly5fz8MMP88orr/DAAw/k+Vyurq6F/jaSvwVuUf3UhbAZK0vgr18HDw9wdTVn3K1amXXdzkDmwEvIxcWF+vXrc/z4cTIzM3PmjQEGDhzIrFmzcm4XN3Xg7+/PSy+9xDvvvAOYFrEffvhhzuOWWsQmJibStGlTgDydD4tqBdu9e3e++eYbAObPn89tt91W5LhOnz5N48aNGT9+PGPHji11q9qvv/4aMBspN2rUKM8PQCFsTms49D3MCoPIjyDkYXh6N3R+pEThnZ5u1nPn7hbYooXzhDdIgJfK22+/zeDBgxkwYAA+Pj4598+aNYstW7YQGBhIhw4d+PTTT4t9rSeffJKIiAhiY2OZNm0aSUlJBAQE0LFjR1577bUCx7/88su8+OKL9OzZM8/9/fv3Z9++fYSEhBTY6X7mzJnMmTOHwMBAvv76az744IMixxQREUFQUFDOXptPP/10sZ8jtzp16hAeHs7TTz9doj8DIcrs0glYMAq+eRQ8GsATa2DERyXuXxIdbYJ7925o3tycdTujErWTtRVpJ1t59erVi5kzZ+bZsDk/+bsWVstfAt9vSqlL4DdvNlWU3t4QHm5dm9fyYut2skIIUb5i1sHyF8pUAp+UZOa4a9Qw89ze3mZdt7NfppEAFzaxefNmRw9BVFZWlMBnZJg2r3v2mODu2RMaNjRflUGFCPDidoAXzq88p+pExTJ7YwyBPnUJ9/Mu9JjImHj2xyUyoU+ulVgZ6bDzU1j3JmSkmhL4ns9CtZLVr586Bdu2wdWrZp67U+lWFDoFh1/EdHd3JyEhQf6BV2JaaxISEnB3dOMI4RCBPnWZuHAvkTHxFh+PjIln4sK9BPrU/f3OMzvh075mk4VmYfDkVuj7conDe88eWL3aTJsMHQqDBkHdusU/z9k4/Azcx8eHuLg4KuKGx8J23N3d86zcEVVHuJ83Mx8KYeLCvcx8KCTPmXh2eOfcn3QJ1r4Ge76A2k1MCXz7ESWarE5JMUsDa9UySwGrV4cOHUq1I5rTcfgqFCFE1ZA/rPPcbtkgbwl89z+XuAQ+M9OsKtm9+/cKyspGVqEIIRwq95n4mG6+zN8ea8K71nn4/GE4s63UJfBxcRAZabYza9oUrGzu6XQkwIUQ5Sbcz5sx3Xz5cF00z/e5lfDoD0yvbve6cNcsCHqoxHMeR4/Cxo1Qp46Z427e3M6Dr4AkwIUQ5SYyJp75207zUVAsYdsnApeg8x9g4GslqqJMTTVruuvVg5Ytze0OHczFyqpIAlwIUS4iY+KZvmAlEY3/R/2jG7lRvz1jE59nfMcHCS8mvLU2Z9w7dpiNhO+7z1ykDAgop8FXUBLgQgi723rsV/Z99TrfqaW4xFeHQW9RM2w8409dsbg6Jbdz58w8d0KCKXsPDy/nwVdgEuBCCLs6uOk7mkRMpgfnoMM9MOjNnBL4opYYApw5AytXmqWBAwaAX8l3Q6sSJMCFEPZx9Rz89AodDy4luXZzuMtyCXx2iO+PSyTcz5v0dLh82Wym0LSpOeP29wc3SasC5I9ECGFbGemwYw6s/2dOCbxHMSXw4X7ehPt5Ex1t9qLMzITRo01oV8YSeFuRABdC2M6ZHbBsEpz/BVoPNLvANyi+2faFC2ae+8KF39u8yhl38az6I1JK/QUYB2jgF+AxrXVK0c8SQlQ6SZdg7TTY82WpS+AvXYLvvjPbmvXpY7oGSm+7kilzgCulmgLPAB201slKqW+AB4F5NhqbEKKiy8yEqAWmBD4lEXpMLFEJfEYGnD8PTZpAgwZw221mV5xitpUV+Vj7S4ob4KGUSgM8gbPWD0kI4RR+O2B2gT+zDZp1h+EzoFHHYp928qRp83rjBjz0kFnX7e9fDuOthMoc4FrrX5VS7wGxQDKwWmu92mYjE0JUTDevwYbppS6Bv3TJzHOfPWvOuocMMeEtys6aKZT6wF1AS+AK8D+l1Bit9fx8x40HxgP4+vpaMVQhhENl7wK/6hW4drZUJfDJybBkCVSrZnbFad++crd5LS/WTKEMBE5qrS8CKKWWAOFAngDXWs8B5oBpJ2vF+wkhHCUhBla+BNFroVGAuUjZrGuRT8nMNN0CfX3NBcr+/c267ho1ymnMVYA1AR4LdFdKeWKmUAYA0uxbiMokLQW2/As2zQDX6jB4OnT9Y7G7wJ85A1u3mjav994LXl7mIqWwLWvmwLcrpRYDe4B0YC9ZZ9pCiEogOgJWvACXTkDHe7J2gb+1yKdcuWIuUMbGmjavgweb8Bb2YdUqFK31NGCajcYihKgIrp7N2gV+KTTwg0eWgl//Yp+Wng4//GCmTrp3NxWUMs9tX1LrJIQwMtJhxydZJfBp0G8KhD9TZAm81mZZYMuWpnKyf39zxu3hUY7jrsIkwIUQELsdlk+C8weg9e0w9J1iS+DPnjXz3AkJZqrE1xdk3+ryJQEuRFWWuwS+TlO4/7/Q/s4ia9mvXTPz3CdPmjavAwea8BblTwJciKoofwl8+NPQZzLUqFXsU1etMiEeGgqBgdJ0ypHkj16Iqua3A2a65Mx28O0Bw94vsgReazhxwmwa7OZmGk7VrGm+hGNJgAtRVeQugfeoB3f9B4JGF7lUJHeb1969TQVlw4blOGZRJAlwISq7/CXwXcbCgGlFlsDfuGE2ED5+3PQr6dsX2rQptxGLEpIAF6IyS4iBFS9CTAQ0LlkJPMDGjWYz4ZAQCA42PUxExSMBLkRlVKAE/m3oOq7IEvgTJ8yu756eZkccV1eoXXRbb+FgEuBCVDa5S+A73Qt3vFlkCXxCgpnnPncOOnc2q0vq1SvH8YoykwAXorK4etbMcx/6rkQl8MnJsGsXHD4M7u7mIqVsrOBcJMCFcHa5S+Az06Hfq9DzGXArum/rzp1w7BgEBECXLrKdmTOSABfCmRUogX8XGrQs/PBYUz3ZoMHvhTgyXeK8JMCFcEZJl0wV5d7/lqgE/soV07fkzBmz63vfvuZipWxp5twkwIVwJpmZEDUf1kyDm1dNt8A+LxdaAn/zJuzeDQcPmqWAPXpAx+L3HRZOQgJcCGdRoAR+BjTqUORTDh40X+3bmykT98I7wwonJAEuREV38xqsfwu2zzYl8CM/NiXwhUyXnD1r/tukiZnjbt5cdsWprCTAhaiotDZLAle9Atd+yyqB/1uhJfBXr5o2r6dOmb7cTZqY5lMS3pWXBLgQFVGeEvhAeGA++IRaPDQtDfbuhf37TV+qsDCzNFBUfhLgQlQkaSmw+QPz5VYDhrwDoU8UWQJ/8iRERZnVJWFhsrKkKpEAF6KiiF4Ly1+Ayyeh030w6E2o3djioefPm46BrVqZLoENGoC3dzmPVzicBLgQjpb4K/z0imn56tUaHv0eWvW1eOiNG7B9O0RHQ/36ZjNhpSS8qyoJcCEcJSMNtn8CG94yJfD9XzXrui2UwKenmznuqChzbbNzZwgKKnLrSlEFSIAL4Qix22DZJLhwENrcYea6iyiBv3jRNJ5q1Qq6dZM2r8KQABeiPN1IgLV/g73zoY4PPLAA/IdZPJWOjzdbmXXoALfeCvfdZ+a6hcgmAS5EecjMNH1L1k4zhTk9n4XbXrJYAp+cbDoFHjliVpS0bWvWc0t4i/ysCnClVD1gLtAJ0MDjWuutthiYEJXGb7+Y6ZK4HeAbnrULfMES+IwMOHAA9uwx3wcGmrluNznNEoWw9n+NfwOrtNb3KaWqA7ICVYhsN6+ZHt3bZ4NHAxg5G4IeLPTK440b5szbx8c0napbt5zHK5xOmQNcKVUHuA0YC6C1TgVSbTMsIZyY1nBwKfz0V1MCH/oY9J9qsQT+8mVTiNO5M9SpA6NGSXCLkrPmDLwVcBH4XCkVBOwGntVa37DJyIRwRgkxZj/KmHVFlsDfvGlWlRw6ZNq8+vub+W4Jb1Ea1gS4G9AZeFprvV0p9W9gMjA190FKqfHAeABfX18r3k6ICqxACfy70PUJcHHNc1hmptmDctcuSE2VNq/COtYEeBwQp7XennV7MSbA89BazwHmAISGhmor3k+Iiun4WnPWffkkBIyCO/5RaAl8WprZYMHLC8LDZWWJsE6ZA1xr/ZtS6oxSqp3W+igwADhku6EJUcHlKYFvU2gJ/NWrZlOF7t2hRg245x6zL6UQ1rJ2FcrTwIKsFSgngMesH5IQFVwJS+BTU02b119+AVdXs57by0vCW9iOVQGutY4CLDcpFqIyylMCPwiGvgP1W+Q5RGs4dgx27DBFOdLmVdiLlAgIURKlKIHXGvbtM8sCBw+GW25xwHhFlSABLkRRLJXA93kZqtfMc9j162a6pHt3syxw+HA54xb2JwEuRGFyl8A372lK4Bu2z3NIerpp8bpvn7ndqhU0bSrhLcqHBLgQ+aVcNRcoiymBj442myvcuAF+fqbNq1ygFOVJAlxUbQsWwJQpEBsLzZrBk3eD26rfS+AH/A086lt86tGj4OEBAwZAY8vLvoWwKwlwUXUtWADjx0NSkrkdGwtT/w1j2sEbEeDTJc/hSUmmgrJzZ3OmPWCAWdctu+IIR5EAF1XXlCm/h3e2NCAiGT77PbwzMsxa7r17zfdNmkDr1lL+LhxPAlxUXbGxlu8/cybn29OnYetWU03ZvLlZZWJNw6nZG2MI9KlLuF/huxBHxsSzPy6RCX38yv5GokpwcfQAhCh3ib/C149AnUIez9V07fRpU0U5dCgMGmR9t8BAn7pMXLiXyJh4i49HxsQzceFeAn2kLaEongS4qDoy0iDyI5jZFY6vhqcfLLDeT3t6cuTRN7lwwdzu0QPuvddssmAL4X7ezHwoxGKIZ4f3zIdCijxDFyKbBLioGmK3wSd9YPWr0LI3PLUd/r4Q5syB5s3RSpF6a3M2jZnDpmYP5wR4tWrgYuN/JZZCXMJblIXSuvw6vIaGhupdu3aV2/sJwY0EWPM3iJoPdZvBkLeh3dA8S0fi4iAyEq5c+X07s/qWVw7aVHZoj+nmy/ztsRLeolBKqd1a6wJ9p+QipqicMjNh75ew9rWsEvjnoM9LBUrgAeLjzeGDBpkLleUl3M+bMd18+XBdNM/0by3hLUpNplCEw8zeGFPoxbxskTHxzN4YU7oXPrcfPrsDfnwWGnaECVvg9tdzwjs1FbZtgxMnzOEBAWYvyvIMbzCfbf72WJ7p35r522OL/bMQIj8JcOEwNl+RkXIVVk6GOX3g0km4+xMYuwwa+gOmS+CRI7BoEezfDwkJ5mmuruarPOWe8550R7tCL2wKURQJcOEwNluRoTUc+NasLtk+G7o8Bk/vytO/5Px5WLIEfv4Z6tUzu+J07WqvT1Y0S5+tqD8LIQojAS4cyuoVGfHR8N+7YfHjZh/KP0bA8BkF+pdcv252gh8wAEaMAG8HTTcX9dkkxEVpySoUUSGUekVGWjJsmgFb/gVu7qbpVOjjObvAZ7d5rVHDzHFn3+fm4Mv2UokpykJWoYgKrVQrMo6vydoF/hQE3J+1C3wjwMymREeb7cxu3ID2udp3Ozq8gRKFcrift6xIESVSAf6XFqLgiozufl4FQywxDla9Aod/AO+28IcfoeVtOQ8nJMCmTXDhgtnGbOBAaNSonD+IEOVIAlw4XP554e5+XnnniTPSYNvHsPiuiWkAABBWSURBVGE66EwzXdLjaXCrnud1MjLMXHffvtCmjbR5FZWfzIELhyrsol72/V8MzCBg7+tw4RC0HWwqKbN2gc/IMMsBb940XQKz7yvvJYFC2JvMgYsKp8gVGY1hZYtFNPppMTdrNqHGgwvNLvBZTp40xTjXrkHLlmbuWykJb1G1SIALh9kfl1gwvDMzYc8XsPY1GqVeJ67jn1jl9Sjj/DsBpi/3zz/D2bPQoAEMG2Y2ERaiKpIAFw5TYEXGuf2w7C/w6y5o3guGvY9PQ3/G5TrExQUSE6FXL/D3t32nQCGciQS4cLyUq7D+n7DjE7ML/N2fQOADoBSZmXDwoKmkHDjQ7EU5erQEtxAgAS4cSWs4uARW/RWunzeFOAOm5lRRnjljtjPLbvOalmaf/txCOCurA1wp5QrsAn7VWg+3fkiiSoiPhhXPw4kNcGswjF4ITc1GwklJsHGjCfC6dWHw4Dy7nAkhstjiDPxZ4DCF7zAoxO/yl8APfS9PCTyYs+xr18zSwE6d5IxbiMJYFeBKKR9gGPAmMMkmIxKVVyEl8FrDkcOmBH7YMBPgo0ZJIY4QxbH2DPxfwEtA7cIOUEqNB8YD+MrvwVVTYhysmgyHfyxQAn/2rJnnTkiAxo0hJcXsMyzhLUTxyhzgSqnhwAWt9W6lVN/CjtNazwHmgKnELOv7CSdURAl8aqqZ5z550qwsGTgQWrVy9ICFcC7WnIH3BEYopYYC7kAdpdR8rfUY2wxNOLXTW2H5JIsl8GCmSZKSIDQUAgMrRqdAIZxNmf/ZaK1fAV4ByDoDf0HCW3AjPmsX+AVmF/isEnit4fgx07vkzjtNn+4RI2SqRAhryHmPsI1cJfCkXodef4HbXoTqNTl/3sxzX7gADRua5lM1akh4C2EtmwS41noDsMEWryWc0Ll9sGxSnhJ4GvqTmQkb1pnVJZ6e0K8ftG4twS2ErcgZuCi7lKuw/k3YMQc8vXJK4DUKhVm/nZkJISEQHGzmvYUQtiMBLkovexf4n6aYEviuT0D/V8GjPidOwM6dMGQI1KljVpcIIexDAlyUTiEl8PHxELkGfvvNtHlNTXX0QIWo/CTARcmkJcOm92HLv8HNI08J/ObNcOgQuLtD796mzavMcwthfxLgonjHVpsS+CunTZvX2/+OrtUoJ6Td3CAgALp0gerVi34pIYTtSICLwiXGwcqX4ciyPCXwsbGwdQXcdhvceuvv+1EKIcqXBLgoKCMNtv0HNrydpwT+yvXqRK6AuDioV8/RgxRCSICLvE5HmjXdFw9D2yFZJfDN2bUL9u41SwF79ICOHaXNqxCOJgEujBvxsHoq7FuYVQL/FZlth6IUKMDDA9q3N71L3N0dPVghBEiAi8xM2DMP1r6epwT+14s12brENJpq29accQshKhYJ8Krs3L6sXeB3Q4veMPQ9rrr7s20DnDoFtWvL2bYQFZkEeFWUkgjr3oSdn2aVwM+BwPs5cFCxbZuZ2w4LM0sDXV2LfzkhhGNIgFclOSXwf4XrF6DrE+j+U9E16uGizBl369YmvD09HT1YIURxJMCrivjjsPx5OLkxqwR+EefdOhO5Cpo1Mxcnmzc3X0II5yABXtlZKIG/7v84O3a55rR5rV/f0YMUQpSFBHhlduwnWPHi7yXwd/yD4+casmmxmU3p3BmCgqTNqxDOSgK8MrpyxuwCf2QZeLeDPywjvVlv3NxMp0BfX+jWzcx5CyGclwR4ZZJTAj/dnGIPmEZ824lEbq9O7VNmRxwvL+nRLURlIQFeWZzaYi5SXjwM7YaS0m86O44158j3Zi1327aOHqAQwtYkwJ3d9YtmF/h9C6GuLzz4FWdqDmXtSsjIMJWUnTtLm1chKiMJcGeVpwT+BvSaRGr3F6heqyZeSWZpYNeuULeuowcqhLAXCXBndDYKlk/KKYFP7P0+m4+0I2MdjBhhlgbKPLcQlZ8EuDPJVwKfduen7EgexaF1imrVTDGO1rKdmRBVhQS4MyhQAj+O+KBXWR5Rj9RU6NDBbGcmjaeEqFokwCu63CXwTUJIvW8R1Vt0pl66mecODjZru4UQVY8EeEWVlgw/v2dK4Kt5ktz/fX5OeozL21wZ1cxsJNy/v6MHKYRwpDIHuFKqGfAl0BjIBOZorf9tq4FVablK4DM6PUhUw7+z91hDXF0hJETmuIUQhjVn4OnA81rrPUqp2sBupdQarfUhG42t6slXAn9j1DKW7OlN8hFo184sC5Q2r0KIbGUOcK31OeBc1vfXlFKHgaaABHhpZaTB1lmw8W3QmrQ+r1Gt91N4ulan5WUT3rfc4uhBCiEqGpvMgSulWgAhwHYLj40HxgP4+vra4u0ql1NbzJrui0dIbz2M7fWnc/y0Lw/2AHc36NXL0QMUQlRUVge4UqoW8C3wnNb6av7HtdZzgDkAoaGh2tr3qzSuX4Q1U2HfV+i6vkR3XcTPvw2B38zKEje5vCyEKIZVMaGUqoYJ7wVa6yW2GVIll5kBu+dBxOuQmkR6j0ksvvAiV3/1xM/PtHmtVcvRgxRCOANrVqEo4P+Aw1rrGbYbUiWWqwQ+w7c3rne+j9st7Wi9C3x8oHFjRw9QCOFMrDkD7wk8AvyilIrKuu+vWusV1g+rkklJhHX/gJ1z0Z7eHG7/KZFXR3Gvm6I+pgReCCFKy5pVKJsBWZFcFK3hl8Xw01/RNy4S33Icq9NfJflaPQICoWZNRw9QCOHM5FKZvVw8Biueh5M/o5uEsLbhN5y8GUKLFtC9O9Sp4+gBCiGcnQS4raUmwab3YMuH6GqeqGHvo7o8RpPDrrSva+a6hRDCFlwcPYBK5egq+E832PQ+Fxrey4L6uzjTeBy4uNKxY9HhPXtjDJEx8UW+fGRMPLM3xth40EIIZyUBbgtXzsCih+GrB0jJ8GBl4+V87/oJLQMalriCMtCnLhMX7i00xCNj4pm4cC+BPrLFjhDCkAC3RnoqbP4AZoVBzDoONX2N+bU2o317cd990LNnyXt0h/t5M/OhEIshnh3eMx8KIdzP2w4fRAjhjGQOvKxObc7aBf4Iut0w1JDpuF/y5XZXaN68bC+ZO8Szw1rCWwhRGAnw0spVAp/i4cvPXovwCRpCh3rQqp71L587xMd082X+9lgJbyGERRLgJZWZAbs/R0e8AalJ/FLveXa6v0CbDp60bGnbtwr382ZMN18+XBfNM/1bS3gLISySAC+Js3th2SQ4u4crdXuzuvYMPJq15a5w8LZDtkbGxDN/eyzP9G/N/O2xdPfzkhAXQhQgAV6U5Cuw/k30zrng6Y26Zy4p3vfRNUXRqpV93jL/nHd3Py+ZAxdCWCSrUCzRGvZ/g57ZFb1jLodqjmNP+E4IHMWtTcovvKHo1SlCiKpNAjy/i8fQX9wJS/5IQroPS73Xcb7zu7QLssEVyiIUtdpEQlwIYYlMoWTLVQKf4eLJ1roziG8+lp49XWnUyP5vvz8uschpkuwQ3x+XKFMpQggAlNblt0lOaGio3rVrV7m9X4kdXUXmihdxSYyFoNEkhr3B+RsNadNGdoAXQjieUmq31rpA4+mqfQZ+JRa9YjLq2HISq/lzImA5Xe7uRV1ACtaFEBVd1Qzw9FTYNovMDe+Y5d21X+dqxyfpFl7d0SMTQogSq3oBnqsE/rT7cA63eovgvr40aeLogQkhROlUnQC/foGMlVNxPbgI6vmSet/XpLgOZnA7cJG1OEIIJ1T5Azwzg8xdn5O55g1UWhJHG71Auyeep3p1T9o7emxCCGGFCnnuabPNDc7u5eZ/BuKy4nnOqyA2B0TSaPRUqO5pw9EKIYRjVMgAt3pzg+QrsPwF9Jx+ZFyKY2uTuWQ8/AN97m1LPfvW4wghRLmpkFMolvpiZyuyP7bWpO35H65rp+CSEo/uOp7YplPoFlhX5rmFEJVOhY01S+XjRYW3vnCUGx/fSbUf/8jlTB8yn1iHy7B38A+W8BZCVE4V8gw8W4k2N0hN4tqKd6m57yPc8GS/7wya3jkWl1tcHTdwIYQoBxU6wKGYzQ2OriTjx5eofT2WE3UewmXQGwR2LOEuwkII4eQqfIBb2tyga90k0n98GY/TK3C9xZ8z3Vfg270nbhX+0wghhO1YFXlKqcHAvwFXYK7WerpNRpWlwOYGzWuz7ZO/E+q5CDflQnq/13Hr9RTNXKvZ8m2FEMIplDnAlVKuwCzgdiAO2KmU+kFrfcgWA8sf3pf2bKLTT88T7nmUn671xm3EPxnQJ9AWbyWEEE7JmvUZYUC01vqE1joVWATcZYtB5Q/v1CXP0uCH4bhkJHO2z9fU+vM8Xvz5vGxuIISo0qyZQmkKnMl1Ow7olv8gpdR4YDyAr69viV44/+YG1Ru25ErAC9Qc8jxNPD1pArK5gRCiyivzhg5KqVHAIK31uKzbjwBhWuunC3tOhd3QQQghKrDCNnSwZgolDmiW67YPcNaK1xNCCFEK1gT4TqCNUqqlUqo68CDwg22GJYQQojhlngPXWqcrpSYCP2GWEX6mtT5os5EJIYQoklXrwLXWK4AVNhqLEEKIUpA2T0II4aQkwIUQwklJgAshhJOSABdCCCdV5kKeMr2ZUheB02V8ujdQ1Wrn5TNXDfKZKz9rP29zrXWBXtnlGuDWUErtslSJVJnJZ64a5DNXfvb6vDKFIoQQTkoCXAghnJQzBfgcRw/AAeQzVw3ymSs/u3xep5kDF0IIkZcznYELIYTIRQJcCCGclFMEuFJqsFLqqFIqWik12dHjsTel1GdKqQtKqQOOHkt5UEo1U0qtV0odVkodVEo96+gx2ZtSyl0ptUMptS/rM7/u6DGVF6WUq1Jqr1JqmaPHUh6UUqeUUr8opaKUUjbd0abCz4FnbZ58jFybJwOjbbV5ckWklLoNuA58qbXu5Ojx2JtS6lbgVq31HqVUbWA3MLKS/x0roKbW+rpSqhqwGXhWa73NwUOzO6XUJCAUqKO1Hu7o8dibUuoUEKq1tnnhkjOcgdtt8+SKSmv9M3DJ0eMoL1rrc1rrPVnfXwMOY/ZcrbS0cT3rZrWsr4p9NmUDSikfYBgw19FjqQycIcAtbZ5cqf9xV2VKqRZACLDdsSOxv6yphCjgArBGa13pPzPwL+AlINPRAylHGlitlNqdtcm7zThDgCsL91X6M5WqSClVC/gWeE5rfdXR47E3rXWG1joYs59smFKqUk+XKaWGAxe01rsdPZZy1lNr3RkYAjyVNUVqE84Q4LJ5chWQNQ/8LbBAa73E0eMpT1rrK8AGYLCDh2JvPYERWXPCi4D+Sqn5jh2S/Wmtz2b99wKwFDMtbBPOEOCyeXIll3VB7/+Aw1rrGY4eT3lQSt2ilKqX9b0HMBA44thR2ZfW+hWttY/WugXm3/E6rfUYBw/LrpRSNbMuzKOUqgncAdhsdVmFD3CtdTqQvXnyYeCbyr55slLqK2Ar0E4pFaeUesLRY7KznsAjmDOyqKyvoY4elJ3dCqxXSu3HnKSs0VpXiWV1VUwjYLNSah+wA1iutV5lqxev8MsIhRBCWFbhz8CFEEJYJgEuhBBOSgJcCCGclAS4EEI4KQlwIYRwUhLgQgjhpCTAhRDCSf0/n4hzn1UUotYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax.plot(new_x, pred_y, 'or')\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recalculate Model Error\n",
    "\n",
    "We can use the same procedure as we just used to predict new point to predict the models predictions of the data point we know, and compare these \n",
    "\n",
    "Note that this is the same as what the fitting procedure does, and, as well see, this should lead to calculating the same error as the model returned to us. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate model predictions for our observed data points\n",
    "preds = theta * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residuals, as returned by the OLS fit, are the just the sum of squares between the model fit and the observed data points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-calculate the residuals 'by hand'\n",
    "error = np.sum(np.subtract(preds, y) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error, returned by the model \t\t: 3.2984\n",
      "Error, as recalculated residuals \t: 3.2984\n"
     ]
    }
   ],
   "source": [
    "# Check that our residuals calculation matches the scipy implementation\n",
    "print('Error, returned by the model \\t\\t: {:1.4f}'.format(residuals[0]))\n",
    "print('Error, as recalculated residuals \\t: {:1.4f}'.format(error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "Note: In practice, you won't use numpy for OLS. Other modules, like statsmodels, have implementations of OLS more explicitly for linear modelling. </div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "See the 'LinearModels' notebook and/or \n",
    "<a href=\"http://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.html\" class=\"alert-link\">OLS in StatsModels</a>.\n",
    "</div>"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}