13-OrdinaryLeastSquares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Least Squares"
   ]
  },
  {
   "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": [
    "#### Ordinary Least Squares\n",
    "\n",
    "Ordinary least squares can be thought of as an error function - a formalization of how we will calculate the error of a model, given the data. \n",
    "\n",
    "Typically, we want to minimize this error term. So the OLS solution to, for example, a linear model, is the model which has the lowest value for the sum of squared errors, calculated as the square of the difference between the model prediction of a data point, and the data point itself.\n",
    "\n",
    "Here, we'll create a minimal dataset, and explore fitting a simple linear model fit to it, using OLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "%matplotlib inline\n",
    "\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "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, 0.5, len(x))\n",
    "y = y + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADRtJREFUeJzt3V9sZHd5xvHnsWctyFKEcaawzWZt\nLKGoERcNHtFpV+IiISiFiHDBRaINggrLN/wJUAmFK8RdLyoEFyska0kLwkqEklRFFNFEsBGKxAQ8\nydImLC2raaYsSVljjCBVJe8wby88W62WzXo958wcz+vvR7Ls8Z717x2t9qujn8/McUQIADD5pqoe\nAABQDoIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASCJ2jgXu/HGG2NhYWGcSwLAxGu3\n27+KiPpux4016AsLC1pfXx/nkgAw8Wx3r+c4tlwAIAmCDgBJ7Bp02w/ZvmD7+cu+90bbT9r+2eDz\n7GjHBADs5nrO0P9B0l1XfO9BSd+NiLdK+u7gMQCgQrsGPSK+L+nXV3z7HklfHXz9VUnvL3kuAMAe\nDbuH/qaIeFmSBp//uLyRAADDGPkvRW2v2F63vb6xsTHq5QBgX2l3t3Ty9Dm1u1sjX2vY69B/aftI\nRLxs+4ikC692YESsSlqVpEajwf3uABwY7e6WTpxqabvX10xtSmvLTS3Nj+4akmHP0L8p6UODrz8k\n6Z/KGQcA8mh1NrXd66sf0sVeX63O5kjXu57LFh+W9ANJt9g+b/sjkv5W0p22fybpzsFjAMBlmotz\nmqlNadrSodqUmotzI11v1y2XiLjvVf7ojpJnAYBUluZntbbcVKuzqebi3Ei3W6Qxv5cLABw0S/Oz\nIw/5Jbz0HwCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAk\nQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCS\nIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiiUNBtf8r2C7aft/2w7deUNRgAYG+GDrrtmyR9QlIj\nIt4maVrSvWUNBgDYm6JbLjVJr7Vdk3SDpJeKjwQgs3Z3SydPn1O7u1X1KOnUhv2LEfEL238n6b8k\n/a+kJyLiiSuPs70iaUWSjh07NuxyABJod7d04lRL272+ZmpTWltuaml+tuqx0iiy5TIr6R5Jb5H0\nJ5IO277/yuMiYjUiGhHRqNfrw08KYOK1Opva7vXVD+lir69WZ7PqkVIpsuXyLkn/GREbEXFR0uOS\n/rKcsQBk1Fyc00xtStOWDtWm1Fycq3qkVIbectHOVkvT9g3a2XK5Q9J6KVMBSGlpflZry021Optq\nLs6x3VKyInvoz9h+VNKzknqSnpO0WtZgAHJamp8l5CNS5AxdEfE5SZ8raRYAQAG8UhQAkiDoAJAE\nQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiC\noANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDlSo3d3SydPn1O5uVT0KEqhVPQBwULW7WzpxqqXt\nXl8ztSmtLTe1ND9b9ViYYJyhAxVpdTa13eurH9LFXl+tzmbVI2HCEXSgIs3FOc3UpjRt6VBtSs3F\nuapHwoRjywWoyNL8rNaWm2p1NtVcnGO7BYURdKBCS/OzhBylYcsFAJIg6ACQBEEHgCQIOgAkQdAB\nIIlCQbf9BtuP2v6p7bO2/6KswQAAe1P0ssUvSfpORHzA9oykG0qYCQAwhKGDbvv1kt4p6cOSFBHb\nkrbLGQsAsFdFtlwWJW1I+nvbz9k+ZfvwlQfZXrG9bnt9Y2OjwHIAgGspEvSapLdL+nJE3CbpfyQ9\neOVBEbEaEY2IaNTr9QLLAQCupUjQz0s6HxHPDB4/qp3AAwAqMHTQI+K/Jf3c9i2Db90h6SelTAUA\n2LOiV7l8XNLa4AqXjqS/Lj4SAGAYhYIeEWckNUqaBQBQAK8UBYAkCDoAJEHQASAJgg4ASRB0AEiC\noANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB\n0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIo\nHHTb07afs/2tMgYCAAynjDP0BySdLeHnAAAKKBR020clvVfSqXLGAQAMq+gZ+hclfUZSv4RZAAAF\nDB1023dLuhAR7V2OW7G9bnt9Y2Nj2OUAALsocoZ+XNL7bL8o6RFJt9v++pUHRcRqRDQiolGv1wss\nBwC4lqGDHhGfjYijEbEg6V5J34uI+0ubDACwJ1yHDgBJ1Mr4IRHxlKSnyvhZAIDhcIYOAEkQdABI\ngqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEHHvtHubunk\n6XNqd7eqHgWYSKW8fS5QVLu7pROnWtru9TVTm9LaclNL87NVjwVMFM7QsS+0Opva7vXVD+lir69W\nZ7PqkYCJQ9CxLzQX5zRTm9K0pUO1KTUX56oeCZg4bLlgX1ian9XaclOtzqaai3NstwBDIOjYN5bm\nZwk5UABbLgCQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEg\nCYK+T3H3HgB7NfS7Ldq+WdLXJL1ZUl/SakR8qazBDjLu3gNgGEXO0HuS/iYi/lRSU9JHbd9azlgH\nG3fvATCMoYMeES9HxLODr38n6aykm8oa7CDj7j0AhlHKDS5sL0i6TdIzZfy8g4679wAYRuGg236d\npMckfTIifnuVP1+RtCJJx44dK7rcgcHdewDsVaGrXGwf0k7M1yLi8asdExGrEdGIiEa9Xi+yHADg\nGoYOum1L+oqksxHxhfJGAgAMo8gZ+nFJH5R0u+0zg4/3lDQXAGCPht5Dj4inJbnEWQAABfBKUQBI\ngqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAk\nQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCS\nIOgAkARBB4AkCDoAJDERQW93t3Ty9Dm1u1tVjwIA+1at6gF20+5u6cSplrZ7fc3UprS23NTS/GzV\nYwHAvlPoDN32Xbb/3fY52w+WNdTlWp1Nbff66od0sddXq7M5imUAYOINHXTb05JOSvorSbdKus/2\nrWUNdklzcU4ztSlNWzpUm1Jzca7sJQAghSJbLu+QdC4iOpJk+xFJ90j6SRmDXbI0P6u15aZanU01\nF+fYbgGAV1Ek6DdJ+vllj89L+vNi41zd0vwsIQeAXRTZQ/dVvhd/cJC9Ynvd9vrGxkaB5QAA11Ik\n6Ocl3XzZ46OSXrryoIhYjYhGRDTq9XqB5QAA11Ik6D+S9Fbbb7E9I+leSd8sZywAwF4NvYceET3b\nH5P0L5KmJT0UES+UNhkAYE8KvbAoIr4t6dslzQIAKGAiXvoPANidI/7gwpTRLWZvSOoO+ddvlPSr\nEseZBDzng4HnnF/R5zsfEbteVTLWoBdhez0iGlXPMU4854OB55zfuJ4vWy4AkARBB4AkJinoq1UP\nUAGe88HAc85vLM93YvbQAQDXNkln6ACAa5iIoI/jRhr7ie2HbF+w/XzVs4yD7Zttn7Z91vYLth+o\neqZRs/0a2z+0/ePBc/581TONi+1p28/Z/lbVs4yD7Rdt/5vtM7bXR7rWft9yGdxI4z8k3amdNwT7\nkaT7IqLU913fT2y/U9Irkr4WEW+rep5Rs31E0pGIeNb2H0lqS3p/8n9jSzocEa/YPiTpaUkPRESr\n4tFGzvanJTUkvT4i7q56nlGz/aKkRkSM/Lr7SThD//8baUTEtqRLN9JIKyK+L+nXVc8xLhHxckQ8\nO/j6d5LOauf99tOKHa8MHh4afOzvs6sS2D4q6b2STlU9S0aTEPSr3Ugj9X/2g8z2gqTbJD1T7SSj\nN9h6OCPpgqQnIyL9c5b0RUmfkdSvepAxCklP2G7bXhnlQpMQ9Ou6kQYmn+3XSXpM0icj4rdVzzNq\nEfH7iPgz7dxL4B22U2+v2b5b0oWIaFc9y5gdj4i3a+f+yx8dbKmOxCQE/bpupIHJNthHfkzSWkQ8\nXvU84xQRv5H0lKS7Kh5l1I5Let9gT/kRSbfb/nq1I41eRLw0+HxB0j9qZxt5JCYh6NxII7nBLwi/\nIulsRHyh6nnGwXbd9hsGX79W0rsk/bTaqUYrIj4bEUcjYkE7/4+/FxH3VzzWSNk+PPhFv2wflvRu\nSSO7em3fBz0iepIu3UjjrKRvZL+Rhu2HJf1A0i22z9v+SNUzjdhxSR/UzhnbmcHHe6oeasSOSDpt\n+1+1c9LyZEQciMv4Dpg3SXra9o8l/VDSP0fEd0a12L6/bBEAcH32/Rk6AOD6EHQASIKgA0ASBB0A\nkiDoAJAEQQeAJAg6ACRB0AEgif8D1ngsRS+G5UsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1138e4eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data\n",
    "f = plt.figure()\n",
    "plt.plot(x, y, '.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observing the data above, we can see that there is some relation between the x and y dimension, but we want a way to measure what that relation is. OLS is a procedure to do so: 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": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reshape 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": [
    "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'). In this simple model, we are therefore implicitly assuming an intercept value of zero. You can fit intercepts (and linear models with more parameters) with OLS, you just need to add them in. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "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)\n",
    "\n",
    "# Pull out theta value from array\n",
    "theta = theta[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.98695402961\n"
     ]
    }
   ],
   "source": [
    "# Check what the OLS derived solution for theta is:\n",
    "print(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The true relationship between y & x is: \t 2\n",
      "OLS calculated relationship between y & x is: \t 1.98695402961\n"
     ]
    }
   ],
   "source": [
    "# Check how good our OLS solution is\n",
    "print('The true relationship between y & x is: \\t', true_rel)\n",
    "print('OLS calculated relationship between y & x is: \\t', theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3701226131131277"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check what the residuals are\n",
    "residuals[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XlcVXX+x/HX9wIKoqCxmIkrluW+\nkAtWk0tWZmrLjEu2zTROi2Jlzdi0TM5My1SjLVj+bC/TanLKtMUll1TMRCGXzIXc0FRA2UGB+/39\nATqmgAgXLtz7fj4ePeDe+733fO7x4dvT93zP+RhrLSIiUvc53F2AiIi4hgJdRMRDKNBFRDyEAl1E\nxEMo0EVEPIQCXUTEQyjQRUQ8hAJdRMRDKNBFRDyEb01uLDQ01LZu3bomNykiUuetX78+1VobdrZx\nNRrorVu3Jj4+viY3KSJS5xlj9lRknKZcREQ8hAJdRMRDKNBFRDxEjc6hl6agoIDk5GTy8/PdXUqt\n5e/vT0REBH5+fu4uRURqMbcHenJyMo0aNaJ169YYY9xdTq1jrSUtLY3k5GTatGnj7nJEpBY765SL\nMeYtY8xhY8zmU547zxiz2Bizo+Rnk8oWkJ+fT0hIiMK8DMYYQkJC9H8wInJWFZlDfwe45rTnJgPf\nWGsvBL4peVxpCvPyaf+I1B0zViQRl5Ra7pi4pFRmrEhy+bbPGujW2m+BI6c9PRx4t+T3d4ERLq7r\nV9y5g0REzkWXiGDGz044mVnZ2VBQ8L/X45JSGT87gS4RwS7fdmVXuTS11v4CUPIzvKyBxphxxph4\nY0x8SkpKpTZ2+g46XVV3kI+PD926daNjx4507dqVqVOn4nQ6y33P7t27mT17dqW2JyKeKzoylNgx\n3bnvgwRmfXWU/3zsZMfX3wD/y6rYMd2Jjgx1+barfdmitXamtTbKWhsVFnbWK1dLdWIHlRbqrthB\nAQEBJCYmsmXLFhYvXsyXX37JlClTyn2PAl1EytIhJJRRTXvz9jtxDD44iA7rb2Tjd4urNcyh8oF+\nyBjTDKDk52HXlVS60kK9Ov61Cw8PZ+bMmcTGxmKtZffu3Vx++eX06NGDHj16EBcXB8DkyZNZuXIl\n3bp1Y9q0aWWOExHvkpUFn35ynMF5M1jU4Q4Cj2/jq3ZPcsciW61hDpVftvg5cDvwbMnPeS6rqByn\nhvrY3i2ZtXZvteygtm3b4nQ6OXz4MOHh4SxevBh/f3927NjB6NGjiY+P59lnn+WFF15gwYIFAOTm\n5pY6TkS8Q34++PtDo4z1jM4cj3/mj9DxRj4M+CPPrjpKzIBW1RrmUIFAN8bMAa4EQo0xycDfKA7y\nj40xfwD2Ar+tziJPFR0ZytjeLXl56U5iBrSrth1krQWKL3waP348iYmJ+Pj4sH379lLHV3SciHiW\noiLYsAG2/pDD78Kewj/xNfwbng+j5hDn14uZsxOIGdCOWWv30icyxL1H6Nba0WW8NNDFtVRIXFIq\ns9burdYd9PPPP+Pj40N4eDhTpkyhadOm/PDDDzidTvz9/Ut9z7Rp0yo0TkQ8x6FDsGIFBB5axm9z\nJuKfvAeifg+DniRuf8GvpoT7RIbU2jl0tzh1zvzBwe3LPFFaFSkpKdx9992MHz8eYwwZGRk0a9YM\nh8PB+++/T1FREQCNGjUiKyvr5PvKGicinmnNGvjq06P0TL6P69JGENDAF+74AoZOOyPMofzFHa5S\nZwK9tBOgrtpBeXl5J5ctDho0iMGDB/O3v/0NgHvvvZd3332XPn36sH37dgIDAwHo0qULvr6+dO3a\nlWnTppU5TkQ8kLWEpXzGmCO9aJs1By57AO5ZDa0vK3exRnWHujkxV1wToqKi7OknCrdu3coll1xS\n7vvOtpqlutd21gYV2U8iUn2OHSs+Km9/wS802/AQ/LQAzu8Cw2OhWdeT42asSKJLRHC5WRSXlMrG\n5Azu/k1khbZtjFlvrY062zi335yrIjYmZ5Qb1if+1duYnOGxgS4i7rNrF6xaaWl15D3CVz8O9hgM\nehL6TgCfX8doRUI6OjK0WrKqTgS6O3eQiHiv3FxYvRrSdiRxdc5EwnNWQqvLYNjLEFKxo+uaVCcC\nXUTEHfbsKiR483QGZj2N8asHQ1+EHreDo3aeflSgi4icIisLMjIgwncTF68aj8lIhPZD4Lp/Q9AF\n7i6vXAp0ERHAWtiyBdavzadnznM0T38JE9AEbn4bOt4AdeA21gp0EfF66enFFwiZvXHclBNDw/wd\n0HUMXP0UNDjP3eVVWO2cCKphycnJDB8+nAsvvJDIyEgmTpzI8ePHAVi+fDlDhw494z0LFiyge/fu\ndO3alQ4dOvB///d/Z4x55513CAsLo1u3bnTr1o3bbrsNgCeeeIIlS5YA8OKLL5Kbm1uN305EypOV\nBZ//J5P2OycxLO1aGvofg7H/hRteq1NhDgp0rLXceOONjBgxgh07drB9+3ays7N59NFHy3xPQUEB\n48aNY/78+fzwww8kJCRw5ZVXljp25MiRJCYmkpiYyHvvvQfA3//+dwYNGgQo0EXc5URXx0a/LGT0\nkT60z3oTet8D96yBdm65s0mVef2Uy9KlS/H39+fOO+8EiptdTJs2jTZt2pR5T/SsrCwKCwsJCQkB\noH79+rRv377C27zjjjsYOnQoBw4c4MCBA/Tv35/Q0FCWLVtW9S8kIuUqLCy+mdbOH1K5KfAv1N/+\nCX5hl8CYd6HFpe4ur0pqXaDPn3/mc5GR0KFD8R/EV1+d+Xr79nDRRcX/4i5e/OvXrr++/O1t2bKF\nnj17/uq5oKAgWrZsyc6dO0t9z3nnncewYcNo1aoVAwcOZOjQoYwePRpHKUuZPvroI1atWgXAxIkT\nT/7DARATE8PUqVNZtmwZoaFaQy9S3Q4ehBXLLWEHP+a3OZPxLcqCKx+Byx4E33ruLq/Kal2g1zRr\nbalNmMt6/oQ33niDTZs2sWTJEl544QUWL17MO++8c8a4kSNHEhsb68qSRaQS4uJgV+I++uc8wAU5\niyHiUhj2CoR7zi01al2gl3dE7etb/uv+/mc/Ij9dx44dmTt37q+ey8zMZN++fURGRpKWllbmezt3\n7kznzp259dZbadOmTamBLiK1gNNJ8/1v0DttCg6HE655FnqNA4ePuytzKa8/KTpw4EByc3NPnrAs\nKipi0qRJ3HHHHTRo0KDU92RnZ7N8+fKTjxMTE2nVqlWltn/6bXhFxDXy82HZMji4aRu8fQ2ttjyM\nT+temHu/gz73eFyYQy08Qq9pxhg+/fRT7r33Xv7xj3/gdDoZMmQITz/99Mkx33zzDREREScfz5kz\nh+eee44//elPBAQEEBgYWOmj83HjxnHttdfSrFkznRQVcZGff4a4lce5OO1FwrOfh/qBMGIGdB1V\nJy4Qqqw6cftc0X4SqYjcXFi1CnK2r2dA9gSC87dAxxvh2n9Bw3B3l1dpHnX7XBGRitiblEOzjU/T\nKftVaNQURsyBi4e4u6wao0AXkTotM7P4ZlotCpbT/psYTPYe6HknXDUF/IPdXV6NUqCLSJ1kLWze\nDD+sPUrfzMcgcxbmvMjivp6tL3N3eW6hQBeROufo0eKbaQXumcfN2Q9TvzAV+t0PV04GvwB3l+c2\nCnQRqVMyM+HLjw9yWcYkWuWW9PUc9jFc0M3dpbmdAl1E6oS8PAjwtwTteI9RqY/jONnXczz4+Lm7\nvFrB6wM9LS2NgQOL76x28OBBfHx8CAsLA+D777+nXj3X3N9hyZIl3HTTTbRp0waApk2bsnDhQqZP\nn07jxo255ZZbeOuttxgyZAjnn3++S7Yp4gkKC2H9etiTkMQNvhPxS16JTy3u6+lOXh/oISEhJCYm\nAvDkk0/SsGFDHnrooV+NsdZirS315lvnon///nz22We/eu6+++47+ftbb71Fjx49FOgiJX75Bb5d\nXkirX17lpuyncfj5wdBp0OOOWtvX0520R8qwc+dOOnXqxN13302PHj3Yt28fjRs3Pvn6hx9+yF13\n3QXAoUOHuPHGG4mKiqJXr1589913Fd7OY489xosvvshHH31EYmIiI0eOpFu3bicbbIh4q9WrYfV/\nNzFozyD6ZD6Oz4X9MfethajfK8zLULuO0L+aDAc3ufYzz+8M1z5bqbf++OOPvP3228yYMYPCwsIy\nx8XExPDnP/+ZPn36sHv3boYOHcrmzZvPGLds2TK6dSs+cTNq1CgmT5588rWRI0fyyiuvEBsbe3KM\niNcqyKftrufom/oSpkHd6uvpTrUr0GuZyMhILr307De8X7JkCdu2bTv5+OjRo+Tl5REQ8OvlU6VN\nuYhIsfz84lvcdgpaQ/iaCTRL2wFdR8PVT9e5VnDuUrsCvZJH0tUlMDDw5O8Oh4NT73uTf6J/FcVz\n7K48gSribZKS4PuVmXRNmUJ4zhsQ3BLGzoV2g9xdWp1SpYkoY8wDxpgtxpjNxpg5xhh/VxVW2zgc\nDpo0acKOHTtwOp18+umnJ18bNGgQ06dPP/n4xEnWc6Vb6YonmrEiibik1FJfy8mBhQvh29mfcfXu\nS7kkp6Sv571rFOaVUOlAN8Y0B2KAKGttJ8AHGOWqwmqjf/3rX1xzzTUMHDjwV7fTnT59OqtXr6ZL\nly506NCB119/vVKff+edd3LXXXfppKh4lC4RwYyfnVBqqO/fnkrTtbdzZ9Ht+Ac3wvxhUfH/qddv\n6IZK675K3z63JNC/A7oCmcBnwMvW2kVlvUe3z6087Sepy+KSUhk/O4HYMd3pFBZK+lFLy4z/UPDF\nnyE/i4Nd76XFsMfAt767S62Vqv32udba/caYF4C9QB6wqLwwFxHvFR0Zyiuju/OHqTu4IyKVsb6P\nQ84ifuJCCofNpnvPaHeX6BGqMuXSBBgOtAEuAAKNMWNLGTfOGBNvjIlPSUmpfKUiUmcdOQKHN57H\nM03juS/vSsLyvuV5cydZt3yhMHehqpwUHQTsstamWGsLgP8CZ/zJWGtnWmujrLVRJy6pL2VMFcrw\nfNo/UpdlZsI3H22n37ZrGeF8jNQmXRiQ9yw+fe8h+sKm7i7Po1Ql0PcCfYwxDYwxBhgIbD3XD/H3\n9yctLU2hVQZrLWlpafj7e+wCIvFQublA4XGCEp7nppR+hJmf2BH9HDdmPcSN/aOZtXZvmatfpHKq\nMoe+1hjzCbABKAQSgJnn+jkREREkJyej6Ziy+fv7/2pVjUhtVlgI8fFwaMN6riuagG/aFhwdb2Dd\nJZP506f7iL2lO9GRofSJDDl5ojQ6MtTdZXsEtzeJFhHPceAArF6ey0X7n6JLzqvQsClm6FTi/HqV\nGt6nrn5RqJetoqtcdIcbEXGJVatgw9zlXLurL11zYjE9b8eMX1tmmEPx6pfYMd3LXKcu56Z2Xfov\nInVT3lHab3uMsLRZ2CZtYfj/+npuTE4q9wj8RKhvTM7QUXoVacpFRColL6/4Zlpd631O6NqHICcV\noid4fV/P6lDtFxaJiPfauRMSvj1IVMpDhObNL+7rOUZ9Pd1NgS4iFZadDatWWgJ+ep/h2Y/hh/p6\n1iYKdBGpsENbf6ZL4kQuOPYttlU/zPUvQ2g7d5clJRToIlKujAzIOFpIywOv0XbpU+Ao7utp1Nez\n1lGgi0ipnE7YtAmS4jbzm/TxcCwB034IXPdvCLrA3eVJKRToInKGI0dg5bJ8Wux6nhE5L2IC1Nez\nLlCgi8ivZGbC6jlruDI9huCC7errWYco0EUEKL6ZVgOfLIJWTuH6lNexQS1hpPp61iUKdBEvV1AA\n69ZB1oZFDMp/AJ/s/dD7HsyAx9QKro5RoIt4seRk+H5ZKp0PTCY67z/Y0Ith5CJo0cvdpUklKNBF\nvNSqlZbj6//DkKzJ1LeZ8JvJmMsfVF/POkyBLuKN0vfRedODBKcvwjaPwgyPhXA1Ia/rFOgiXiIv\nD1avctLT+SZN1j9JsHXCNc9ieo0Dh4+7yxMXUKCLeIHt22HLiu30TZlAk+PfQdv+cP2L0KS1u0sT\nF9J1uyI1aMaKpLM2cohLSmXGiiSXbC87G77+ooCM+c8z7EA/ws1PMOI1uPVThbkHUqCL1KAuEcHl\nduc50ZKtS0SwS7aXunEDvRKv5NKsf+LocB2OCeug2xhd7emhFOgiNai8lmuu6q+Zng57k3Jh4aO0\nXjyQxvXTYNQczG/fgYbhVfwGUpsp0EVqWGmh7oowdzohMRHiPlhBkw/7wppY6HE7jvFr4eIhrvwK\nUkvppKiIG5wa6mN7t2TW2r1VCvPUVFiz9CgX7n6cIbnv42zcFkb8r6+neAcFuoibREeGMrZ3S15e\nupOYAe0qHeaZmZA453MGpj9EgDMV+t2PQ309vZKmXETcJC4plVlr9xIzoB2z1u496+qX0+XkAFkH\nCfpqLIPSbsU/tCnmj0vhqikKcy+lI3QRNzh9zrxPZEiF59ALCmDd95bC+FlclvsojqLivp4O9fX0\nejpCF6lhpZ0ALW/1y6mSk+HLD3bRcsVwrjg6Hpp2hnvi4LIHFOaiQBepSeWtZjlbqH+7vJB9H7/C\ndXv6cgEJMHQajjvmQ0hkTZUvtZymXERq0MbkjHKnVU6E+sbkjF+PObiZnonjCcxMwF40BDNUfT3l\nTMZaW2Mbi4qKsvHx8TW2PZG6LDcX4lYe49K85wneOA0CmsC1z6mvpxcyxqy31kadbZyO0EVqoW3b\nIGnZd0SnTiC4UH09pWIU6CK1SFYWxC3PImLrFK7NfQMb1AKGqa+nVEyVAt0Y0xh4A+gEWOD31to1\nrihMxBulr1tEv8QHCCzaD73vxqG+nnIOqnqE/hLwtbX2ZmNMPaCBC2oS8Srp6ZB1KI0WmyfTYtPH\nOEMuxoxQX085d5UOdGNMEHAFcAeAtfY4cNw1ZYl4PqcTfki0pK/6hL4Zf8HaTMxvJuNQX0+ppKoc\nobcFUoC3jTFdgfXARGttzqmDjDHjgHEALVu2rMLmRDxHaip8vziZjrsfpPuxhTibRWFuUF9PqZqq\nXFjkC/QAXrPWdgdygMmnD7LWzrTWRllro8LCwqqwORHPkJnhZNv7r3PV9t60KFoJ1zyL44+LFOZS\nZVU5Qk8Gkq21a0sef0IpgS4ixbKzoWHedoLmx9AvYw3ONgNwDHsRmrRyd2niISod6Nbag8aYfcaY\n9tbabcBA4EfXlSbiGY4fh3XfFVA//kV6Zj+HqdcARryGo+toXSAkLlXVVS4TgA9KVrj8DNxZ9ZJE\nPMfevfDjkg1cemACIYWbcXa4ATPkObWCk2pRpUC31iYCZ70cVcQbrVyaS9D6p7g651WcgU3h5jk4\n1ApOqpGuFBWpDj+voFd8DPVzd+PscSc+g6eAf7C7qxIPp0AXcZHcXPhueTq90h+j4fb3qX9eW/jd\nFzjU11NqiAJdxAV++gkOfPM5fdIeIsAW9/VEfT2lhinQRaogMxO+X3qItj89xID8zykK74IZ8TFc\n0M3dpYkXUqCLVJa15KyaxeUbH8WPfOzAJ/GJVl9PcR8Fusg5OnoUcpJ3EbFhIs12raAooh+OES9D\naDt3lyZeToEuUkFOJyQmFHH829eIyvwntp4fZug0fHrcAQ615xX3U6CLVEBKCiQu3EzXPRMIL9hA\nUbshmGHq6ym1iwJd5Cwyjxxj33vPMzBrGta/MQx/Gx/19ZRaSIEuUobsbGh4ZC1Bn0+gR9Y2ijqP\nxuda9fWU2kuBLnKa48chfnUWwfF/p0PO65jgFjB2Lj7q6ym1nAJd5BR790LS14u59ND9BBbtx/a6\nGzNQfT2lblCgi5RYvSiN8ITJ9M/7mMImF2NuXIRRX0+pQxToItbC5rn0iv8zPgWZ2Csm43uF+npK\n3aNAF6+VkwPx3yTTK+VBAvYtxK95FAx7BZp2cHdpIpWiQBevYy38tNVJ+pI36Xt0Cr6OIrj6Gej9\nJ3D4uLs8kUpToItXycyE9Yt2cPH2CVxyfA2FrfrjGPEiNGnt7tJEqkyBLt6jqIDj37zEFZufg3oB\nMOI1fNXXUzyIAl083pEjkP/zBi5YF0PooU0Uth+B7/XPq6+neBwFunisoiL4IT4X35XP0Ck7Ftuo\nKWbUbHwvvs7dpYlUCwW6eKTDh+HHL7+l+74Ygot2UdjtDnyvngIBjd1dmki1UaCLx8k8nM6Rdx7n\nytz3KAhqCzcswLfN5e4uS6TaKdDFY2RlQaPk+QR98RCN8lIo6jMRv4GPqK+neA0FutR5x4/D+hWH\naLrhYRrlzYPzO2PGfISP+nqKl1GgS522e5flwFcf0CPlUeqZPIr6P4HPZTHq6yleSYEuddaaL3fR\nIvF+oo8v53izaBw3vQyhF7q7LBG3UaBL3eMsgu9eo9f6p8D44BwylXpRd6qvp3g9BbrUGdnZkLho\nC5cemED91PX4XHQNXDcVgpu7uzSRWkGBLrWetfDT5mMcW/IC0RlTcdZvDDe9CZ1u0mX7IqdQoEut\nlpEBm75cS8ekCTQp3EZBh1H4DX1GfT1FSlHlQDfG+ADxwH5r7dCqlyRS4lgWzi/+Tr/tr1PYIAJG\nzsXvQvX1FCmLK47QJwJbgSAXfJYIaWlw/MclNIu/nyYZyRT2/BN+gx9XX0+Rs6hSoBtjIoDrgKeA\nB11SkXitoiLYtDaNwFWPcGHuR9jQ9pg/LMJXfT1FKqSqR+gvAn8GGrmgFvFihw5afp4/l26//IX6\nNoPCfn/Bt/8k9fUUOQeVDnRjzFDgsLV2vTHmynLGjQPGAbRs2bKymxMPlrU/mfz3J9E3/2uOhfXE\ncXMsDvX1FDlnVTlC7wcMM8YMAfyBIGPMLGvt2FMHWWtnAjMBoqKibBW2Jx4mM8NJ0Pa3aLT4SQKL\niii66hnq91VfT5HKqnSgW2sfAR4BKDlCf+j0MBcpzbFj8MPSHbRIiCHoeBy0vRLH9S+pr6dIFWkd\nutSo3T8XcHTBS/Q48hz4BeAc9iqO7mN0gZCIC7gk0K21y4Hlrvgs8Vzr5iXQZtMEWhdu4li7EdQf\n/hw0auruskQ8ho7Qpfodz4XlzxCVEEth/XCcN35A/Q66Bk3E1RToUm2ys2Hzgm/puT8Gv6xdmB63\n43fV39XXU6SaKNDF5ayFnxLTcSx5nD4573G8YVu4fT60ucLdpYl4NAW6uFR6Omz/fAEdd08iwJnC\n8UsnUm+w+nqK1AQFurhO1iF85j5Mr33zyG/cGcfvPqTeBd3dXZWI11CgS5WlplicCR8QvuFRGhXk\nUfCbJ/C/Qn09RWqaAl0qragINq/cRcia+4k4thzbMhoz7GX81NdTxC0U6FIpBw8UcXDeDDoe/ifG\n4UPB1VPx662+niLupECXc5adtAXHhxPoVrCevBbXEHCz+nqK1AY6nJIKy0g7BkufouEHVxDi2EPh\niDcJ+P2HlQ7zGSuSiEtKLXdMXFIqM1YkVerzRbyNAl3O6tgx2DBvLc7XLodvn4NON+MT8z2+3W6u\n0j1YukQEM352QpmhHpeUyvjZCXSJCK70NkS8iQJdyrVrWxa7pj9M94SraeCXS9HouXDj/0FgSJU/\nOzoylNgx3UsN9RNhHjumO9GRoVXelog3UKBLmRI/WULox31pn/k6+V3GUf/+Nfi0d22T5tJCXWEu\nUjk6KSpnykmDhY/QbfNH5Ddqj715IQGtelfb5k4N9bG9WzJr7V6FuUglKNDlpKxMy855c+my/y/4\nHE+HK/6M/xUP1Uhfz+jIUMb2bsnLS3cSM6CdwlykEhTogrWwfV0yAUsn0T3/a/JDeuBz5+fQtGON\n1RCXlMqstXuJGdCOWWv30icyRKEuco4U6F4u/aiTPXPf4pL9T+JjCjl25dP4X3F3jfb1PH3OvE9k\niObQRSpBJ0W9WeoO6s2+jq7Jkyhs2hOf8Wuof+V9bg1zKH/1i4iUTYHuhVIPFXBk3gvwWj8aZG2h\n4LrpNPjTZ3Bemxqto7zVLAp1kXOnKRcvUlgIW5cmcMG6CZxXsAl7yXDMkOfxc1Nfz43JGeVOq5wI\n9Y3JGZp6EakABbqX+GVvLumfPUPHI7EU1Aun4KZZ+HW+3q013f2byLOOiY4MVZiLVJAC3Qtkb/6W\nBp/G0KxoFzntbyNwxD/U11PEAynQPVjGoXSC1z5Oww3vUdCoDUXD5xPYTn09RTyVAt0D5efDjvkL\naLt1EtYexkTH4HflI1CvgbtLE5FqpED3MLs3H4KvHqZzzjxygzrh/O0cfFr0cHdZIlIDFOiewlq2\nzP6AyJ2P4kceuX0fp8GgierrKeJFFOie4MguWHA/HX9eTk5IX+qNfIUG4errKeJtdGFRLVGZ7j2Z\n6UX8+GYszlejITkervs3gfd9iUNhLuKVFOi1xLl077EWtq/awrFXr6LDvkfJa3o53LcWLr1LTZpF\nvJj+9tcSFe3ec0lQI7bNeIp2S64g2Lmb/KFvEnjXRxAc4abKRaS2UKDXImft3uOXRIP3L+fiQ8+R\n3eYm6j2wDv+oqvX1FBHPUemTosaYFsB7wPmAE5hprX3JVYV5q9K69zw98CI6fv8v2PY69YOaU/C7\nTwjqcJW7SxWRWqYqq1wKgUnW2g3GmEbAemPMYmvtjy6qzWud6N7z0pIkHgk/Qr9vHqBhUTL20j9i\nBj2BX/1G7i5RRGqhSge6tfYX4JeS37OMMVuB5oACvYriklL5dPF2PvD/iMuy53O0XiSFtyzEr231\n9fUUkbrPJevQjTGtge7A2lJeGweMA2jZsqUrNufR4namMO+dV/kq8F0akMG2thO4bXc002wk0e4u\nTkRqtSqfFDXGNATmAvdbazNPf91aO9NaG2WtjQoLC6vq5jzairhNFM4ayb98p+IX0grnH1fQ/rZ/\nMu2WXmr0ICJnVaVAN8b4URzmH1hr/+uakrxPfp6TZdNj6fHVYPo5tsDgp6h/3xJ8m3cC1L1HRCqm\nKqtcDPAmsNVaO9V1JXmXPet34L84hv75cRwI7kfgndMh5MxWcOreIyJnU5U59H7ArcAmY0xiyXN/\ntdZ+WfWyvEBRATvfe4nWe57D6eNPzqDpXNDvlnLXlKt7j4iUpyqrXFYBuqKlMg4kwLwJtDu0iYzm\nwwga9QL13NTXU0Q8h+62WIMy03JJ+fgZ2h6OxQSGwchZBF/i3r6eIuI5FOg1wOmEXcu+JWxNDJGF\nu8hsdxtBN6mvp4i4lgK9mh0bfYfFAAAIMElEQVT9JZ30/zxB5JF3ya3fmvybPifokt+4uywR8UAK\n9Oq0dQGNPp9EcN5h0jtMoPGIv6qvp4hUGwV6NUjdfQj/pQ/TcO88fJt2omDUHBq3Ul9PEaleCnQX\nKiyw7J43m4gtf8WPXJz9H8dx2UT81NdTRGqAAt1FDm7djZ1/P+1yl5HRuC+O371MvQsucndZIuJF\nFOhV5Swib8UMQr79JxgHGZe9QPCAP6gVnIjUOAV6FWTs+JHg5eMJ2L+e3IjB1LthKsEhLdxdloh4\nKQV6JeRlHeOXj/9Nq31TcfoH4bjxDRp0Vis4EXEvBfo52hf3PY2WjadtwTaORPyOxqOehYYh7i5L\nRESBXmHHstn7zt9p8ctM8vyakz38E87rrr6eIlJ7KNArYscSWHA/LTKSORL5R8777RMYf/X1FJHa\nRYFejsyDR8j4+BFaHPkQQi/C/P5rQlr2cXdZIiKlUqCXwllk2ff1XMLX/4XmznSOdH6Y84Y9BH7+\n7i5NRKRMCvTTHN2zn/xPHqRV1tdkNOwBN8/jvNad3F2WiMhZKdBPcDph/dsEL/objQoLSYt6ipAh\n94DDx92ViYhUiAIdSN22k8BvYgg4vBpHmysovOYlQpq2dXdZIiLnxKsDvSC/gAOfvEzznf/C6fDH\neX0sjh5jqacLhESkDvLaQD+UmIDfVxNodWwTqeHDCB71PI7zznd3WSIileZ9gV6Qx7GFzxAWH8sx\n31CODn6f0Ohh7q5KRKTKvCrQMzeuJGh5DPWP/EzWhbfSYPg/CGjYxN1liYi4hFcEeu6RdI589AQR\nh96lKKg1Prd9TqO26uspIp7F4wN9/zcLaBI3iQuKDnO43QRCf/tXqK++niLieTw30LMPc/Dth2me\n9hkZAZ0oGjGH8Pbq6ykinqtWttWZsSKJuKTUcsfEJaUyY0XSGc9bp4WEDyD2UsKPfsmhTo8RNGk5\njRTmIuLhamWgd4kIZvzshDJDPS4plfGzE+gSEfyr5zP37CZl2g0w714IuxjHPatpevPDGF81aRYR\nz1crAz06MpTYMd1LDfUTYR47pjvRkaEAOAuLSP5kOgHv9KVJ9jpSe70Ad34FYWrSLCLeo1YGOpQe\n6qWF+dFtP3L034OJ2PxXMhpfRtHd3xE65I9q0iwiXqdWnxQ9NdTH9m7JrLV7/xfmhcdg5b9pvHIq\nxwgi5Yo3COuvvp4i4r2qFOjGmGuAlwAf4A1r7bMuqeoU0ZGhjO3dkpeX7iRmQDuiI0NJS/yehssm\nUD/jJ0zn3+Ez8BnCGoe6etMiInVKpQPdGOMDTAeuApKBdcaYz621P7qqOCieZpm1di8xA9rxSdx2\nrvrpNTodeYc8v+Y4R/8HR/vB6JSniEjVjtB7ATuttT8DGGM+BIYDLgv0U+fMLzycyO/NZILS9rMl\n5Fba//5pHA3V11NE5ISqnDlsDuw75XFyyXMucWqY90z8B2ELb8Lh14A1fT/itoxhrDt0zFWbEhHx\nCFUJ9NLOPtozBhkzzhgTb4yJT0lJqdAHn76apf75bcno8hCBk1bS75qry1zSKCLizaoy5ZIMtDjl\ncQRw4PRB1tqZwEyAqKioMwK/NBuTM361NJF+Ezn1EqITq182Jmf8b4yIiJcz1lYoY898ozG+wHZg\nILAfWAeMsdZuKes9UVFRNj4+vlLbExHxVsaY9dbaqLONq/QRurW20BgzHlhI8bLFt8oLcxERqV5V\nWodurf0S+NJFtYiISBXo+ngREQ+hQBcR8RAKdBERD6FAFxHxEAp0EREPUel16JXamDEpwJ5Kvj0U\n8LZLQ/WdvYO+s+er6vdtZa0NO9ugGg30qjDGxFdkYb0n0Xf2DvrOnq+mvq+mXEREPIQCXUTEQ9Sl\nQJ/p7gLcQN/ZO+g7e74a+b51Zg5dRETKV5eO0EVEpBx1ItCNMdcYY7YZY3YaYya7u57qZox5yxhz\n2Biz2d211ARjTAtjzDJjzFZjzBZjzER311TdjDH+xpjvjTE/lHznKe6uqaYYY3yMMQnGmAXurqUm\nGGN2G2M2GWMSjTHVev/wWj/lUtKMejunNKMGRru6GXVtYoy5AsgG3rPWdnJ3PdXNGNMMaGat3WCM\naQSsB0Z4+J+xAQKttdnGGD9gFTDRWvudm0urdsaYB4EoIMhaO9Td9VQ3Y8xuIMpaW+3r7uvCEfrJ\nZtTW2uPAiWbUHsta+y1wxN111BRr7S/W2g0lv2cBW3Fhf9rayBbLLnnoV/Jf7T66cgFjTARwHfCG\nu2vxRHUh0Ku1GbXULsaY1kB3YK17K6l+JVMPicBhYLG11uO/M/Ai8GfA6e5CapAFFhlj1htjxlXn\nhupCoFeoGbXUfcaYhsBc4H5rbaa766lu1toia203ivvx9jLGePT0mjFmKHDYWrve3bXUsH7W2h7A\ntcB9JVOq1aIuBHqFmlFL3VYyjzwX+MBa+19311OTrLXpwHLgGjeXUt36AcNK5pQ/BAYYY2a5t6Tq\nZ609UPLzMPApxdPI1aIuBPo64EJjTBtjTD1gFPC5m2sSFyo5QfgmsNVaO9Xd9dQEY0yYMaZxye8B\nwCDgJ/dWVb2stY9YayOsta0p/nu81Fo71s1lVStjTGDJiX6MMYHAYKDaVq/V+kC31hYCJ5pRbwU+\n9vRm1MaYOcAaoL0xJtkY8wd311TN+gG3UnzElljy3xB3F1XNmgHLjDEbKT5oWWyt9YplfF6mKbDK\nGPMD8D3whbX26+raWK1ftigiIhVT64/QRUSkYhToIiIeQoEuIuIhFOgiIh5CgS4i4iEU6CIiHkKB\nLiLiIRToIiIe4v8BgwXYOv/rlLgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118ba4588>"
      ]
     },
     "metadata": {},
     "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 Fit')\n",
    "ax.plot(x, theta*x, label='True Fit')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The prediction for a new x of 2.5 is 4.967\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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4VFX+x/H3mRBICCRgCqKhBkXp\nJVKCulJERQQsuxSx7bqsBYKKuriW1d21rLqgblB+2JWirqwi2ADpBJFAIkWkRFpASgLpCSSZ8/sj\ngaUkISSTzGTm83oen2Rmzsz9zuXhw/Xcc+/XWGsREZHaz+HuAkRExDUU6CIiXkKBLiLiJRToIiJe\nQoEuIuIlFOgiIl5CgS4i4iUU6CIiXkKBLiLiJerU5MbCwsJsy5Yta3KTIiK13tq1a1OtteFnG1ej\ngd6yZUsSEhJqcpMiIrWeMWZXRcZpykVExEso0EVEvIQCXUTES9ToHHppCgoKSElJIT8/392leKyA\ngAAiIyPx9/d3dyki4sHcHugpKSk0bNiQli1bYoxxdzkex1pLWloaKSkptGrVyt3liIgHO+uUizHm\nHWPMQWPMxpOeO88Ys8AYs63kZ+PKFpCfn09oaKjCvAzGGEJDQ/V/MCJyVhWZQ38PuPa05yYC31lr\nLwK+K3lcaQrz8mn/iNQeU5cmE5+cWu6Y+ORUpi5Ndvm2zxro1tplwOHTnh4KvF/y+/vAMBfXdQp3\n7iARkXPRKTKEsTMTT2RWdjYUFPzv9fjkVMbOTKRTZIjLt13ZVS5NrLW/ApT8jChroDFmjDEmwRiT\ncOjQoUpt7PQddLqq7iA/Pz+6dOlC+/bt6dy5M5MmTcLpdJb7np07dzJz5sxKbU9EvFdMVBhxo7py\n/4xEpn99hP984mTbN98B/8uquFFdiYkKc/m2q33ZorV2mrU22lobHR5+1itXS3V8B5UW6q7YQYGB\ngSQlJbFp0yYWLFjAV199xTPPPFPuexToIlKWdqFhjGjSk3ffi2fg/gG0W3sT679fUK1hDpUP9APG\nmKYAJT8Puq6k0pUW6tXxr11ERATTpk0jLi4Oay07d+7kiiuuoFu3bnTr1o34+HgAJk6cyPLly+nS\npQuTJ08uc5yI+JasLPjs02MMzJvK/HZ3EnRsC1+3eZo759tqDXOo/LLFL4A7gBdKfs5xWUXlODnU\nR/dszvTVu6tlB7Vu3Rqn08nBgweJiIhgwYIFBAQEsG3bNkaOHElCQgIvvPACL7/8MvPmzQMgNze3\n1HEi4hvy8yEgABpmrGVk5lgCMn+C9jfxUeAfeWHFEWL7tajWMIcKBLoxZhZwFRBmjEkB/kpxkH9i\njPkDsBv4bXUWebKYqDBG92zOa4u2E9uvTbXtIGstUHzh09ixY0lKSsLPz4+tW7eWOr6i40TEuxQV\nwbp1sPnHHH4X/iwBSW8Q0OB8GDGLeP8eTJuZSGy/NkxfvZteUaHuPUK31o4s46X+Lq6lQuKTU5m+\nene17qBffvkFPz8/IiIieOaZZ2jSpAk//vgjTqeTgICAUt8zefLkCo0TEe9x4AAsXQpBBxbz25zx\nBKTsgujfw4Cnid9bcMqUcK+oUI+dQ3eLk+fMHxrYtswTpVVx6NAh7rnnHsaOHYsxhoyMDJo2bYrD\n4eDDDz+kqKgIgIYNG5KVlXXifWWNExHvtGoVfP3ZEbqn3M/1acMIrF8H7vwSBk8+I8yh/MUdrlJr\nAr20E6Cu2kF5eXknli0OGDCAgQMH8te//hWA++67j/fff59evXqxdetWgoKCAOjUqRN16tShc+fO\nTJ48ucxxIuKFrCX80OeMOtyD1lmz4PIH4d6V0PLychdrVHeom+NzxTUhOjrann6icPPmzVx66aXl\nvu9sq1mqe22nJ6jIfhKR6nP0aPFRedsLfqXpuofh53lwficYGgdNO58YN3VpMp0iQ8rNovjkVNan\nZHDPb6IqtG1jzFprbfTZxrn95lwVsT4lo9ywPv6v3vqUDK8NdBFxnx07YMVyS4vDHxCx8kmwR2HA\n09B7HPidGqMVCemYqLBqyapaEeju3EEi4rtyc2HlSkjblsw1OeOJyFkOLS6HIa9BaMWOrmtSrQh0\nERF32LWjkJCNU+if9RzGvy4MfgW63QEOzzz9qEAXETlJVhZkZEBknQ1csmIsJiMJ2g6C6/8FwRe4\nu7xyKdBFRABrYdMmWLs6n+45L3Jh+quYwMZwy7vQ/kaoBbexVqCLiM9LTy++QMjsjufmnFga5G+D\nzqPgmmeh/nnuLq/CPHMiqIalpKQwdOhQLrroIqKiohg/fjzHjh0DYMmSJQwePPiM98ybN4+uXbvS\nuXNn2rVrx//93/+dMea9994jPDycLl260KVLF26//XYAnnrqKRYuXAjAK6+8Qm5ubjV+OxEpT1YW\nfPGfTNpun8CQtOtoEHAURv8XbnyjVoU5KNCx1nLTTTcxbNgwtm3bxtatW8nOzubxxx8v8z0FBQWM\nGTOGuXPn8uOPP5KYmMhVV11V6tjhw4eTlJREUlISH3zwAQB/+9vfGDBgAKBAF3GX410dG/76LSMP\n96Jt1tvQ8164dxW0ccudTarM56dcFi1aREBAAHfddRdQ3Oxi8uTJtGrVqsx7omdlZVFYWEhoaCgA\n9erVo23bthXe5p133sngwYPZt28f+/bto2/fvoSFhbF48eKqfyERKVdhYfHNtLb/mMrNQX+m3tZP\n8Q+/FEa9D80uc3d5VeJxgT537pnPRUVBu3bFfxBff33m623bwsUXF/+Lu2DBqa/dcEP529u0aRPd\nu3c/5bng4GCaN2/O9u3bS33Peeedx5AhQ2jRogX9+/dn8ODBjBw5EkcpS5k+/vhjVqxYAcD48eNP\n/MMBEBsby6RJk1i8eDFhYVpDL1Ld9u+HpUss4fs/4bc5E6lTlAVXPQaXPwR16rq7vCrzuECvadba\nUpswl/X8cW+99RYbNmxg4cKFvPzyyyxYsID33nvvjHHDhw8nLi7OlSWLSCXEx8OOpD30zXmQC3IW\nQORlMOTfEOE9t9TwuEAv74i6Tp3yXw8IOPsR+enat2/P7NmzT3kuMzOTPXv2EBUVRVpaWpnv7dix\nIx07duS2226jVatWpQa6iHgAp5ML975Fz7RncDiccO0L0GMMOPzcXZlL+fxJ0f79+5Obm3vihGVR\nURETJkzgzjvvpH79+qW+Jzs7myVLlpx4nJSURIsWLSq1/dNvwysirpGfD4sXw/4NW+Dda2mx6RH8\nWvbA3Pc99LrX68IcPPAIvaYZY/jss8+47777+Pvf/47T6WTQoEE899xzJ8Z89913REZGnng8a9Ys\nXnzxRf70pz8RGBhIUFBQpY/Ox4wZw3XXXUfTpk11UlTERX75BeKXH+OStFeIyH4J6gXBsKnQeUSt\nuECosmrF7XNF+0mkInJzYcUKyNm6ln7Z4wjJ3wTtb4Lr/gkNItxdXqV51e1zRUQqYndyDk3XP0eH\n7NehYRMYNgsuGeTusmqMAl1EarXMzOKbaTUrWELb72Ix2bug+11w9TMQEOLu8mqUAl1EaiVrYeNG\n+HH1EXpnPgGZ0zHnRRX39Wx5ubvLcwsFuojUOkeOFN9MK2jXHG7JfoR6hanQ5wG4aiL4B7q7PLdR\noItIrZKZCV99sp/LMybQIrekr+eQT+CCLu4uze0U6CJSK+TlQWCAJXjbB4xIfRLHib6eY8HP393l\neQSfD/S0tDT69y++s9r+/fvx8/MjPDwcgB9++IG6dV1zf4eFCxdy880306pVKwCaNGnCt99+y5Qp\nU2jUqBG33nor77zzDoMGDeL88893yTZFvEFhIaxdC7sSk7mxznj8U5bj58F9Pd3J5wM9NDSUpKQk\nAJ5++mkaNGjAww8/fMoYay3W2lJvvnUu+vbty+eff37Kc/fff/+J39955x26deumQBcp8euvsGxJ\nIS1+fZ2bs5/D4e8PgydDtzs9tq+nO2mPlGH79u106NCBe+65h27durFnzx4aNWp04vWPPvqIu+++\nG4ADBw5w0003ER0dTY8ePfj+++8rvJ0nnniCV155hY8//pikpCSGDx9Oly5dTjTYEPFVK1fCyv9u\nYMCuAfTKfBK/i/pi7l8N0b9XmJfBs47Qv54I+ze49jPP7wjXvVCpt/7000+8++67TJ06lcLCwjLH\nxcbG8uijj9KrVy927tzJ4MGD2bhx4xnjFi9eTJcuxSduRowYwcSJE0+8Nnz4cP79738TFxd3YoyI\nzyrIp/WOF+md+iqmfu3q6+lOnhXoHiYqKorLLjv7De8XLlzIli1bTjw+cuQIeXl5BAaeunyqtCkX\nESmWn198i9sOwauIWDWOpmnboPNIuOa5WtcKzl08K9AreSRdXYKCgk787nA4OPm+N/nH+1dRPMfu\nyhOoIr4mORl+WJ5J50PPEJHzFoQ0h9Gzoc0Ad5dWq1RpIsoY86AxZpMxZqMxZpYxJsBVhXkah8NB\n48aN2bZtG06nk88+++zEawMGDGDKlCknHh8/yXqudCtd8UZTlyYTn5xa6ms5OfDtt7Bs5udcs/My\nLs0p6et53yqFeSVUOtCNMRcCsUC0tbYD4AeMcFVhnuif//wn1157Lf379z/ldrpTpkxh5cqVdOrU\niXbt2vHmm29W6vPvuusu7r77bp0UFa/SKTKEsTMTSw31vVtTabL6Du4quoOAkIaYP8wv/j/1eg3c\nUGntV+nb55YE+vdAZyAT+Bx4zVo7v6z36Pa5laf9JLVZfHIqY2cmEjeqKx3Cw0g/Ymme8R8KvnwU\n8rPY3/k+mg15AurUc3epHqnab59rrd1rjHkZ2A3kAfPLC3MR8V0xUWH8e2RX/jBpG3dGpjK6zpOQ\nM5+fuYjCITPp2j3G3SV6hapMuTQGhgKtgAuAIGPM6FLGjTHGJBhjEg4dOlT5SkWk1jp8GA6uP4/n\nmyRwf95VhOct4yVzF1m3fqkwd6GqnBQdAOyw1h6y1hYA/wXO+JOx1k6z1kZba6OPX1JfypgqlOH9\ntH+kNsvMhO8+3kqfLdcxzPkEqY070S/vBfx630vMRU3cXZ5XqUqg7wZ6GWPqG2MM0B/YfK4fEhAQ\nQFpamkKrDNZa0tLSCAjw2gVE4qVyc4HCYwQnvsTNh/oQbn5mW8yL3JT1MDf1jWH66t1lrn6RyqnK\nHPpqY8ynwDqgEEgEpp3r50RGRpKSkoKmY8oWEBBwyqoaEU9WWAgJCXBg3VquLxpHnbRNONrfyJpL\nJ/Knz/YQd2tXYqLC6BUVeuJEaUxUmLvL9gpubxItIt5j3z5YuSSXi/c+S6ec16FBE8zgScT79yg1\nvE9e/aJQL1tFV7noDjci4hIrVsC62Uu4bkdvOufEYbrfgRm7uswwh+LVL3Gjupa5Tl3OjWdd+i8i\ntVPeEdpueYLwtOnYxq1h6P/6eq5PSS73CPx4qK9PydBRehVpykVEKiUvr/hmWp3rfkHY6ochJxVi\nxvl8X8/qUO0XFomI79q+HRKX7Sf60MOE5c0t7us5Sn093U2BLiIVlp0NK5ZbAn/+kKHZT+CP+np6\nEgW6iFTYgc2/0ClpPBccXYZt0Qdzw2sQ1sbdZUkJBbqIlCsjAzKOFNJ83xu0XvQsOIr7ehr19fQ4\nCnQRKZXTCRs2QHL8Rn6TPhaOJmLaDoLr/wXBF7i7PCmFAl1EznD4MCxfnE+zHS8xLOcVTKD6etYG\nCnQROUVmJqyctYqr0mMJKdiqvp61iAJdRIDim2nV98siePkz3HDoTWxwcxiuvp61iQJdxMcVFMCa\nNZC1bj4D8h/EL3sv9LwX0+8JtYKrZRToIj4sJQV+WJxKx30Ticn7DzbsEhg+H5r1cHdpUgkKdBEf\ntWK55dja/zAoayL1bCb8ZiLmiofU17MWU6CL+KL0PXTc8BAh6fOxF0ZjhsZBhJqQ13YKdBEfkZcH\nK1c46e58m8ZrnybEOuHaFzA9xoDDz93liQso0EV8wNatsGnpVnofGkfjY99D675wwyvQuKW7SxMX\n0nW7IjVo6tLkszZyiE9OZerSZJdsLzsbvvmygIy5LzFkXx8izM8w7A247TOFuRdSoIvUoE6RIeV2\n5znekq1TZIhLtpe6fh09kq7isqx/4Gh3PY5xa6DLKF3t6aUU6CI1qLyWa67qr5meDruTc+Hbx2m5\noD+N6qXBiFmY374HDSKq+A3EkynQRWpYaaHuijB3OiEpCeJnLKXxR71hVRx0uwPH2NVwySBXfgXx\nUDopKuIGJ4f66J7Nmb56d5XCPDUVVi06wkU7n2RQ7oc4G7WGYf/r6ym+QYEu4iYxUWGM7tmc1xZt\nJ7Zfm0qHeWYmJM36gv7pDxPoTIU+D+BQX0+fpCkXETeJT05l+urdxPZrw/TVu8+6+uV0OTlA1n6C\nvx7NgLTbCAhrgvnjIrj6GYW5j9IRuogbnD5n3isqtMJz6AUFsOYHS2HCdC7PfRxHUXFfT4f6evo8\nHaGL1LDSToCWt/rlZCkp8NWMHTRfOpQrj4yFJh3h3ni4/EGFuSjQRWpSeatZzhbqy5YUsueTf3P9\nrt5cQCIMnozjzrkQGlVT5YuH05SLSA1an5JR7rTK8VBfn5Jx6pj9G+meNJagzETsxYMwg9XXU85k\nrLU1trHo6GibkJBQY9sTqc1ycyF++VEuy3uJkPWTIbAxXPei+nr6IGPMWmtt9NnG6QhdxANt2QLJ\ni78nJnUcIYXq6ykVo0AX8SBZWRC/JIvIzc9wXe5b2OBmMER9PaViqhToxphGwFtAB8ACv7fWrnJF\nYSK+KH3NfPokPUhQ0V7oeQ8O9fWUc1DVI/RXgW+stbcYY+oC9V1Qk4hPSU+HrANpNNs4kWYbPsEZ\neglmmPp6yrmrdKAbY4KBK4E7Aay1x4BjrilLxPs5nfBjkiV9xaf0zvgz1mZifjMRh/p6SiVV5Qi9\nNXAIeNcY0xlYC4y31uacPMgYMwYYA9C8efMqbE7Ee6Smwg8LUmi/8yG6Hv0WZ9NozI3q6ylVU5UL\ni+oA3YA3rLVdgRxg4umDrLXTrLXR1tro8PDwKmxOxDtkZjjZ8uGbXL21J82KlsO1L+D443yFuVRZ\nVY7QU4AUa+3qksefUkqgi0ix7GxokLeV4Lmx9MlYhbNVPxxDXoHGLdxdmniJSge6tXa/MWaPMaat\ntXYL0B/4yXWliXiHY8dgzfcF1Et4he7ZL2Lq1odhb+DoPFIXCIlLVXWVyzhgRskKl1+Au6pekoj3\n2L0bflq4jsv2jSO0cCPOdjdiBr2oVnBSLaoU6NbaJOCsl6OK+KLli3IJXvss1+S8jjOoCdwyC4da\nwUk10pWiItXhl6X0SIilXu5OnN3uwm/gMxAQ4u6qxMsp0EVcJDcXvl+STo/0J2iw9UPqndcafvcl\nDvX1lBqiQBdxgZ9/hn3ffUGvtIcJtMV9PVFfT6lhCnSRKsjMhB8WHaD1zw/TL/8LiiI6YYZ9Ahd0\ncXdp4oMU6CKVZS05K6ZzxfrH8Scf2/9p/GLU11PcR4Euco6OHIGclB1ErhtP0x1LKYrsg2PYaxDW\nxt2liY9ToItUkNMJSYlFHFv2BtGZ/8DW9ccMnoxftzvBofa84n4KdJEKOHQIkr7dSOdd44goWEdR\nm0GYIerrKZ5FgS5yFpmHj7Lng5fonzUZG9AIhr6Ln/p6igdSoIuUITsbGhxeTfAX4+iWtYWijiPx\nu059PcVzKdBFTnPsGCSszCIk4W+0y3kTE9IMRs/GT309xcMp0EVOsns3JH+zgMsOPEBQ0V5sj3sw\n/dXXU2oHnZoX3zZjBrRsCQ4H+U2a4f9oX/ruu4WA4AaYP8zHMegFhbnUGgp08V0zZsCYMbBrF1hL\nwMEUzv/vEuyxa6hz/zI1aZZaR4EuPsv52OPFd9Q6iSkAM/MHNWmWWkmBLj7HWtj8kxOzZ1fpA3bv\nrtmCRFxEgS4+JTMTlszeRqPPB2FCylhH3rx5zRYl4iIKdPEdRQUc++5lrtzYhwg2w6P3Qv36p46p\nXx+efdY99YlUkZYtitc7fBjyf1nHBWtiCTuwgcK2w6hzw0vFfT1bxMDjjxdPszRvXhzmt97q7pJF\nKkWBLl6rqAh+TMilzvLn6ZAdh23YBDNiJnUuuf5/g269VQEuXkOBLl7p4EH46atldN0TS0jRDgq7\n3Emda56BwEbuLk2k2ijQxetkHkzn8HtPclXuBxQEt4Yb51Gn1RXuLkuk2inQxWtkZUHDlLkEf/kw\nDfMOUdRrPP79H1NfT/EZCnSp9Y4dg7VLD9Bk3SM0zJsD53fEjPoYP/X1FB+jQJdabecOy76vZ9Dt\n0OPUNXkU9X0Kv8tj1ddTfJICXWqtVV/toFnSA8QcW8KxpjE4bn4Nwi5yd1kibqNAl9rHWQTfv0GP\ntc+C8cM5aBJ1o+9SX0/xeQp0qTWysyFp/iYu2zeOeqlr8bv4Wrh+EoRc6O7SRDyCAl08nrXw88aj\nHF34MjEZk3DWawQ3vw0dblZfT5GTKNDFo2VkwIavVtM+eRyNC7dQ0G4E/oOfV19PkVJUOdCNMX5A\nArDXWju46iWJlDiahfPLv9Fn65sU1o+E4bPxv0h9PUXK4ooj9PHAZiDYBZ8lQloaHPtpIU0THqBx\nRgqF3f+E/8An1QpO5CyqFOjGmEjgeuBZ4CGXVCQ+q6gINqxOI2jFY1yU+zE2rC3mD/Opo1ZwIhVS\n1SP0V4BHgYYuqEV82IH9ll/mzqbLr3+mns2gsM+fqdN3glrBiZyDSge6MWYwcNBau9YYc1U548YA\nYwCaqxOMlCJrbwr5H06gd/43HA3vjuOWOBxN2rm7LJFapypH6H2AIcaYQUAAEGyMmW6tHX3yIGvt\nNGAaQHR0tK3C9sTLZGY4Cd76Dg0XPE1QURFFVz9Pvd5/Aoefu0sTqZUqHejW2seAxwBKjtAfPj3M\nRUpz9Cj8uGgbzRJjCT4WD62vwnHDq9C4pbtLE6nVtA5datTOXwo4Mu9Vuh1+EfwDcQ55HUfXUbpA\nSMQFXBLo1tolwBJXfJZ4rzVzEmm1YRwtCzdwtM0w6g19ERo2cXdZIl5DR+hS/Y7lwpLniU6Mo7Be\nBM6bZlCvna5BE3E1BbpUm+xs2DhvGd33xuKftQPT7Q78r/6b+nqKVBMFurictfBzUjqOhU/SK+cD\njjVoDXfMhVZXurs0Ea+mQBeXSk+HrV/Mo/3OCQQ6D3HssvHUHai+niI1QYEurpN1AL/Zj9Bjzxzy\nG3XE8buPqHtBV3dXJeIzFOhSZamHLM7EGUSse5yGBXkU/OYpAq5UX0+RmqZAl0orKoKNy3cQuuoB\nIo8uwTaPwQx5DX/19RRxCwW6VMr+fUXsnzOV9gf/gXH4UXDNJPx7qq+niDsp0OWcZSdvwvHROLoU\nrCWv2bUE3qK+niKeQIdTUmEZaUdh0bM0mHEloY5dFA57m8Dff1TpMJ+6NJn45NRyx8QnpzJ1aXKl\nPl/E1yjQ5ayOHoV1c1bjfOMKWPYidLgFv9gfqNPllirdg6VTZAhjZyaWGerxyamMnZlIp8iQSm9D\nxJco0KVcO7ZksWPKI3RNvIb6/rkUjZwNN/0fBIVW+bNjosKIG9W11FA/HuZxo7oSExVW5W2J+AIF\nupQp6dOFhH3Sm7aZb5LfaQz1HliFX1vXNmkuLdQV5iKVo5OicqacNPj2Mbps/Jj8hm2xt3xLYIue\n1ba5k0N9dM/mTF+9W2EuUgkKdDkhK9Oyfc5sOu39M37H0uHKRwm48uEa6esZExXG6J7NeW3RdmL7\ntVGYi1SCAl2wFrauSSFw0QS65n9Dfmg3/O76Apq0r7Ea4pNTmb56N7H92jB99W56RYUq1EXOkQLd\nx6UfcbJr9jtcuvdp/EwhR696joAr76nRvp6nz5n3igrVHLpIJeikqC9L3UbdmdfTOWUChU264zd2\nFfWuut+tYQ7lr34RkbIp0H1Q6oECDs95Gd7oQ/2sTRRcP4X6f/oczmtVo3WUt5pFoS5y7jTl4kMK\nC2HzokQuWDOO8wo2YC8dihn0Ev5u6uu5PiWj3GmV46G+PiVDUy8iFaBA9xG/7s4l/fPnaX84joK6\nERTcPB3/jje4taZ7fhN11jExUWEKc5EKUqD7gOyNy6j/WSxNi3aQ0/Z2gob9XX09RbyQAt2LZRxI\nJ2T1kzRY9wEFDVtRNHQuQW3U11PEWynQvVB+PmybO4/Wmydg7UFMTCz+Vz0Gdeu7uzQRqUYKdC+z\nc+MB+PoROubMITe4A87fzsKvWTd3lyUiNUCB7i2sZdPMGURtfxx/8sjt/ST1B4xXX08RH6JA9waH\nd8C8B2j/yxJyQntTd/i/qR+hvp4ivkYXFnmIynTvyUwv4qe343C+HgMpCXD9vwi6/yscCnMRn6RA\n9xDn0r3HWti6YhNHX7+adnseJ6/JFXD/arjsbjVpFvFh+tvvISravefS4IZsmfosbRZeSYhzJ/mD\n3ybo7o8hJNJNlYuIp1Cge5Czdu/xT6b+h1dwyYEXyW51M3UfXENAdNX6eoqI96j0SVFjTDPgA+B8\nwAlMs9a+6qrCfFVp3Xue638x7X/4J2x5k3rBF1Lwu08Jbne1u0sVEQ9TlVUuhcAEa+06Y0xDYK0x\nZoG19icX1eazjnfveXVhMo9FHKbPdw/SoCgFe9kfMQOewr9eQ3eXKCIeqNKBbq39Ffi15PcsY8xm\n4EJAgV5F8cmpfLZgKzMCPuby7LkcqRtF4a3f4t+6+vp6ikjt55J16MaYlkBXYHUpr40BxgA0b97c\nFZvzavHbDzHnvdf5Ouh96pPBltbjuH1nDJNtFDHuLk5EPFqVT4oaYxoAs4EHrLWZp79urZ1mrY22\n1kaHh4dXdXNebWn8BgqnD+efdSbhH9oC5x+X0vb2fzD51h5q9CAiZ1WlQDfG+FMc5jOstf91TUm+\nJz/PyeIpcXT7eiB9HJtg4LPUu38hdS7sAKh7j4hUTFVWuRjgbWCztXaS60ryLbvWbiNgQSx98+PZ\nF9KHoLumQOiZreDUvUdEzqYqc+h9gNuADcaYpJLn/mKt/arqZfmAogK2f/AqLXe9iNMvgJwBU7ig\nz63lrilX9x4RKU9VVrmsAHRFS2XsS4Q542hzYAMZFw4heMTL1HVTX08R8R6622INykzL5dAnz9P6\nYBwmKByGTyfkUvf29RQR76FeaB7HAAAIT0lEQVRArwFOJ+xYvIzwVbFEFe4gs83tBN+svp4i4loK\n9Gp25Nd00v/zFFGH3ye3Xkvyb/6C4Et/4+6yRMQLKdCr0+Z5NPxiAiF5B0lvN45Gw/6ivp4iUm0U\n6NUgdecBAhY9QoPdc6jTpAMFI2bRqIX6eopI9VKgu1BhgWXnnJlEbvoL/uTi7PskjsvH46++niJS\nAxToLrJ/807s3Adok7uYjEa9cfzuNepecLG7yxIRH6JArypnEXlLpxK67B9gHGRc/jIh/f6gVnAi\nUuMU6FWQse0nQpaMJXDvWnIjB1L3xkmEhDZzd1ki4qMU6JWQl3WUXz/5Fy32TMIZEIzjpreo31Gt\n4ETEvRTo52hP/A80XDyW1gVbOBz5OxqNeAEahLq7LBERBXqFHc1m93t/o9mv08jzv5DsoZ9yXlf1\n9RQRz6FAr4htC2HeAzTLSOFw1B8577dPYQLU11NEPIsCvRyZ+w+T8cljNDv8EYRdjPn9N4Q27+Xu\nskRESqVAL4WzyLLnm9lErP0zFzrTOdzxEc4b8jD4B7i7NBGRMinQT3Nk117yP32IFlnfkNGgG9wy\nh/NadnB3WSIiZ6VAP87phLXvEjL/rzQsLCQt+llCB90LDj93VyYiUiEKdCB1y3aCvosl8OBKHK2u\npPDaVwlt0trdZYmInBOfDvSC/AL2ffoaF27/J05HAM4b4nB0G01dXSAkIrWQzwb6gaRE/L8eR4uj\nG0iNGELIiJdwnHe+u8sSEak03wv0gjyOfvs84QlxHK0TxpGBHxIWM8TdVYmIVJlPBXrm+uUEL4ml\n3uFfyLroNuoP/TuBDRq7uywREZfwiUDPPZzO4Y+fIvLA+xQFt8Tv9i9o2Fp9PUXEu3h9oO/9bh6N\n4ydwQdFBDrYZR9hv/wL11NdTRLyP9wZ69kH2v/sIF6Z9TkZgB4qGzSKirfp6ioj38si2OlOXJhOf\nnFrumPjkVKYuTT7jeeu0kDgD4i4j4shXHOjwBMETltBQYS4iXs4jA71TZAhjZyaWGerxyamMnZlI\np8iQU57P3LWTQ5NvhDn3QfglOO5dSZNbHsHUUZNmEfF+HhnoMVFhxI3qWmqoHw/zuFFdiYkKA8BZ\nWETKp1MIfK83jbPXkNrjZbjrawhXk2YR8R0eGehQeqiXFuZHtvzEkX8NJHLjX8hodDlF93xP2KA/\nqkmziPgcjz4penKoj+7ZnOmrd/8vzAuPwvJ/0Wj5JI4SzKEr3yK8r/p6iojvqlKgG2OuBV4F/IC3\nrLUvuKSqk8REhTG6Z3NeW7Sd2H5tiIkKIy3pBxosHke9jJ8xHX+HX//nCW8U5upNi4jUKpUOdGOM\nHzAFuBpIAdYYY76w1v7kquKgeJpl+urdxPZrw6fxW7n65zfocPg98vwvxDnyPzjaDkSnPEVEqnaE\n3gPYbq39BcAY8xEwFHBZoJ88Z37RwSR+byYSnLaXTaG30fb3z+FooL6eIiLHVeXM4YXAnpMep5Q8\n5xInh3n3pL8T/u3NOPzrs6r3x9yeMYQ1B466alMiIl6hKoFe2tlHe8YgY8YYYxKMMQmHDh2q0Aef\nvpql3vmtyej0MEETltPn2mvKXNIoIuLLqjLlkgI0O+lxJLDv9EHW2mnANIDo6OgzAr8061MyTlma\nSJ/xnHwJ0fHVL+tTMv43RkTExxlrK5SxZ77RmDrAVqA/sBdYA4yy1m4q6z3R0dE2ISGhUtsTEfFV\nxpi11tros42r9BG6tbbQGDMW+JbiZYvvlBfmIiJSvaq0Dt1a+xXwlYtqERGRKtD18SIiXkKBLiLi\nJRToIiJeQoEuIuIlFOgiIl6i0uvQK7UxYw4Buyr59jDA1y4N1Xf2DfrO3q+q37eFtTb8bINqNNCr\nwhiTUJGF9d5E39k36Dt7v5r6vppyERHxEgp0EREvUZsCfZq7C3ADfWffoO/s/Wrk+9aaOXQRESlf\nbTpCFxGRctSKQDfGXGuM2WKM2W6MmejueqqbMeYdY8xBY8xGd9dSE4wxzYwxi40xm40xm4wx491d\nU3UzxgQYY34wxvxY8p2fcXdNNcUY42eMSTTGzHN3LTXBGLPTGLPBGJNkjKnW+4d7/JRLSTPqrZzU\njBoY6epm1J7EGHMlkA18YK3t4O56qpsxpinQ1Fq7zhjTEFgLDPPyP2MDBFlrs40x/sAKYLy19ns3\nl1btjDEPAdFAsLV2sLvrqW7GmJ1AtLW22tfd14Yj9BPNqK21x4Djzai9lrV2GXDY3XXUFGvtr9ba\ndSW/ZwGbcWF/Wk9ki2WXPPQv+c+zj65cwBgTCVwPvOXuWrxRbQj0am1GLZ7FGNMS6Aqsdm8l1a9k\n6iEJOAgssNZ6/XcGXgEeBZzuLqQGWWC+MWatMWZMdW6oNgR6hZpRS+1njGkAzAYesNZmurue6mat\nLbLWdqG4H28PY4xXT68ZYwYDB621a91dSw3rY63tBlwH3F8ypVotakOgV6gZtdRuJfPIs4EZ1tr/\nuruemmStTQeWANe6uZTq1gcYUjKn/BHQzxgz3b0lVT9r7b6SnweBzyieRq4WtSHQ1wAXGWNaGWPq\nAiOAL9xck7hQyQnCt4HN1tpJ7q6nJhhjwo0xjUp+DwQGAD+7t6rqZa19zFobaa1tSfHf40XW2tFu\nLqtaGWOCSk70Y4wJAgYC1bZ6zeMD3VpbCBxvRr0Z+MTbm1EbY2YBq4C2xpgUY8wf3F1TNesD3Ebx\nEVtSyX+D3F1UNWsKLDbGrKf4oGWBtdYnlvH5mCbACmPMj8APwJfW2m+qa2Mev2xRREQqxuOP0EVE\npGIU6CIiXkKBLiLiJRToIiJeQoEuIuIlFOgiIl5CgS4i4iUU6CIiXuL/AdUR2Ci5j4GyAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118ba4588>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ax.plot(new_x, pred_y, 'or')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We can also see what the model would predict for all the points we did observe\n",
    "preds = theta * x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Residuals are the just the sum of squares between the model fit and the observed data points\n",
    "# Re-calculate the residuals 'by hand'\n",
    "error = np.sum(np.subtract(preds, y) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error from : 1.37012261311\n",
      "Error from : 1.37012261311\n"
     ]
    }
   ],
   "source": [
    "# Check that our residuals calculation matches the scipy implementation\n",
    "print('Error from :', residuals[0])\n",
    "print('Error from :', 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 [default]",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}