{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Data Fitting with Least Squares"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as npla\n",
        "import scipy.linalg as spla\n",
        "import matplotlib.pyplot as pt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Suppose we are modeling a relationship between $x$ and $y$, and the \"true\" relationship is $y = a+bx$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f1323d14ac8>]"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH+NJREFUeJzt3Xd4FWXexvHvQyDU0DskhA4BQgtFsKCgIqJS9F1dK66i\n23SL0lGxYl1dXQt2d3Fdl4QigiIqiopIEdIIEEJJKEkgJoT05DzvH2T3YlEEciaZU+7PdXmRE4Zn\nfmPCzTDMucdYaxERkcBRy+0BRETEWQp2EZEAo2AXEQkwCnYRkQCjYBcRCTAKdhGRAKNgFxEJMAp2\nEZEAo2AXEQkwtd3YacuWLW1kZKQbuxYR8VubNm06bK1tdbrtXAn2yMhINm7c6MauRUT8ljFm75ls\np0sxIiIBRsEuIhJgFOwiIgFGwS4iEmAU7CIiAeaMg90Y84YxJssYk3jC55obYz4xxuys/LFZ9Ywp\nIiJn6mzO2N8Cxp70uRnAp9ba7sCnla9FRMRFZxzs1tovgZyTPn0V8Hblx28DExyaS0QkoPxQUMq8\nD5I4WlxW7fvy9hp7G2vtQYDKH1ufakNjzFRjzEZjzMbs7Gwvdysi4h+stXwYf5CL//IFf1+3l+/S\nTj4/dl6NvfPUWrsAWAAQExOjJ2iLSMDLPFrM3CWJrErOJLpjE/7+q2H0bte42vfrbbBnGmPaWWsP\nGmPaAVlODCUi4s+stby/MZ2HP9xGabmHWeN6cevIztQOqZkbEb0N9mXAzcD8yh+Xej2RiIgf23ek\nkJmL4/k69QhDOzfn8cnRdG7ZsEZnOONgN8b8ExgFtDTGZAD3czzQ3zfG/ArYB1xTHUOKiPi6Co/l\nrW/28NTH2wmpZXh4Ql9+OTSCWrVMjc9yxsFurb3uFD812qFZRET80o7MfKYtimdLei4X9mzFIxP7\n0b5pfdfmcaW2V0QkEJSWe3j5i1288FkqDeuG8OwvBnDVgPYYU/Nn6SdSsIuIVMHW9Fymx8aTciif\nK/q35/4romjZqK7bYwEKdhGRs1JUWsGzq3fw6to0WoXV5dWbYrg4qo3bY/0PBbuIyBlat+sIM+Pi\n2XOkkOuGhjNzXG8a16vj9lg/omAXETmN/OIyHluZwrvr9xHRvAHv3jaMEd1auj3WKSnYRUR+xmcp\nmcyKSyQrv5jbzu3Mny/pSf3QELfH+lkKdhGRn3DkWAkPLk9m6ZYD9GwTxss3DmZAeFO3xzojCnYR\nkRNYa1m29QDzPkgmv7iMu0d357cXdiO0tv88l0jBLiJS6VBeMXOWJLB6Wxb9w5vyxORoerYNc3us\ns6ZgF5GgZ63lvQ3pPPrhNso8HuZc3pspIzsT4kIdgBMU7CIS1PYcLmBmXALr0o5wTpcWzJ/cj04t\nara0y2kKdhEJShUeyxtf7ebpT7ZTp1YtHpvUj2uHhLteB+AEBbuIBJ3th/KZtmgrWzPyGNO7DQ9P\n6EvbJvXcHssxCnYRCRql5R7+9nkqL65JpXG9Ovz1uoFcEd0uIM7ST6RgF5GgsCU9l2mLtrIj8xgT\nB3Zg7vgomjcMdXusaqFgF5GAVlhazjOrdvDG17tp07geb9wSw0W9fKu0y2kKdhEJWN+kHmZGXAL7\ncgq5flgEMy7rRZgPlnY5TcEuIgEnr6iMx1Zs470N6US2aMB7U4czvEsLt8eqMQp2EQkonyRnMmdJ\nAtn5Jdxxfhf+eHEP6tXx7dIupynYRSQgHD5WwgPLklgef5BebcN49aYYojv6R2mX0xTsIuLXrLUs\n3XKAeR8kUVBSwZ8v7sGdo7pSJ8R/SrucpmAXEb91ILeIOUsS+Swli4ERx0u7urfxv9IupynYRcTv\neDyWd7/bx/yVKVR4LHPHR3HLiEi/Le1ymoJdRPzK7sMFzIiNZ/3uHM7t1pLHJvUjvHkDt8fyKY4E\nuzHmj8BtgAUSgCnW2mIn1hYRASiv8PD6V7t55pMdhNauxROTo7kmpmPA1QE4wetgN8Z0AO4Coqy1\nRcaY94Frgbe8XVtEBGDbwaNMWxRPwv48Lolqw0MT+tKmceCUdjnNqUsxtYH6xpgyoAFwwKF1RSSI\nlZRX8MJnqby0ZhdNG9Thb78cxLh+bXWWfhpeB7u1dr8x5ilgH1AErLLWrjp5O2PMVGAqQEREhLe7\nFZEAt2nvD0yPjSc16xiTBnVg7uVRNAvQ0i6neX2jpzGmGXAV0BloDzQ0xtxw8nbW2gXW2hhrbUyr\nVq283a2IBKiCknLmfZDE1S9/Q2FJOW9NGcIz/zdAoX4WnLgUMwbYba3NBjDGxAEjgH84sLaIBJGv\ndh5mRlw8GT8UcdM5nZg2theN6urmvbPlxP+xfcBwY0wDjl+KGQ1sdGBdEQkSeYVlPLIimfc3ZtCl\nZUPev+MchnZu7vZYfsuJa+zrjTGLgM1AOfA9sMDbdUUkOHyUeIi5SxPJKSjlN6O6ctfo7kFX2uU0\nR/6OY629H7jfibVEJDhk5RfzwLIkViQcIqpdY968ZQh9OzRxe6yAoItXIlKjrLXEbd7Pg8uTKSqr\n4N5LezL1/C5BXdrlNAW7iNSY/blFzIpL4Isd2Qzu1IzHJ0fTrXUjt8cKOAp2Eal2Ho9l4fq9zF+Z\nggXmXdmHG4d3opZKu6qFgl1EqtWu7GPMiI1nw54fOK97Sx6dqNKu6qZgF5FqUVbh4dW1aTy7eif1\n64Tw1DX9mTyog+oAaoCCXUQcl7g/j+mx8SQdOMrYPm15cEIfWoeptKumKNhFxDHFZRU8/9lOXv4i\njWYNQnnp+kFc1q+d22MFHQW7iDhi454cpsXGk5ZdwORBHZk7vjdNG6jfxQ0KdhHxSkFJOU98lMI7\n3+6lfZP6vHPrUM7voaI/NynYRaTKvtiRzay4BA7kFXHzOZHce2lPGqq0y3X6CojIWcstLOWh5duI\n3ZxB11YN+fcd5xATqdIuX6FgF5GzsjLhIHOXJpFbWMpvL+zK7y9SaZevUbCLyBnJOlrMfUuT+Cjp\nEH07NObtW4fQp71Ku3yRgl1Efpa1ln9vyuDh5ckUl3uYPrYXt5/Xmdoq7fJZCnYROaX0nEJmLU5g\n7c7DDI1szvzJ/ejSSqVdvk7BLiI/UuGxvLNuD09+vB0DPHRVH64fptIuf6FgF5H/kZqVz7RF8Wze\nl8sFPVrx6KR+dGha3+2x5Cwo2EUEOF7a9coXu/jrp6k0qBvCM//Xn4kDVdrljxTsIkJCRh73LtpK\nyqF8Lo9uxwNX9KFVWF23x5IqUrCLBLHisgqeXb2TV9em0aJhKK/cOJhL+7R1eyzxkoJdJEh9tzuH\n6bHx7D5cwC9iwpl1eW+a1K/j9ljiAAW7SJDJLy7j8Y9S+Me3+whvXp+Ftw1jZLeWbo8lDlKwiwSR\nz7dnMTsugYNHi/nVuZ358yU9aBCqGAg0jnxFjTFNgdeAvoAFbrXWrnNibRHxXk5BKQ8tT2bx9/vp\n3roRsb8ewaCIZm6PJdXEqT+qnwM+stZebYwJBfSkWhEfYK1lefxBHliWRF5RGXeN7s5vL+xK3doq\n7QpkXge7MaYxcD5wC4C1thQo9XZdEfFO5tFiZi9OZPW2TKI7NmHh7cPo1bax22NJDXDijL0LkA28\naYzpD2wC7rbWFjiwtoicJWst/9qQziMrtlFa7mH2uN5MGRmp0q4g4sRXujYwCHjJWjsQKABmnLyR\nMWaqMWajMWZjdna2A7sVkZPtO1LI9a+tZ0ZcAlHtGvPxH87n9vO7KNSDjBNn7BlAhrV2feXrRfxE\nsFtrFwALAGJiYqwD+xWRShUey1vf7OGpj7cTUsvw6MR+XDskXKVdQcrrYLfWHjLGpBtjelprtwOj\ngWTvRxORM7Ej83hp15b0XC7q1ZpHJvalXROVdgUzp+6K+T2wsPKOmDRgikPrisgplJZ7eGnNLl74\nfCeN6tbmuWsHcGX/9irtEmeC3Vq7BYhxYi0ROb2t6blMj40n5VA+V/RvzwNXRNGikUq75Di95UzE\njxSVVvCX1Tt4bW0arcPq8dpNMYyJauP2WOJjFOwifmLdriPMjItnz5FCrhsazsxxvWlcT6Vd8mMK\ndhEfd7S4jPkrU3h3/T46tWjAu7cPY0RXlXbJqSnYRXzYp9symb04kaz8Ym4/rzN/urgn9UNVByA/\nT8Eu4oOOHCth3gfJLNt6gJ5twnj5xsEMCG/q9ljiJxTsIj7EWsuyrQeY90Ey+cVl/GFMd34zqhuh\ntfXOUTlzCnYRH3Ewr4g5ixP5NCWLAeFNeeLqaHq0CXN7LPFDCnYRl3k8ln9u2MdjK1Io93iYc3lv\npozsTIjqAKSKFOwiLtpzuIAZcfF8m5bDiK4tmD8pmogWepyBeEfBLuKCCo/l9a/SeHrVDkJDavHY\npOOlXaoDECco2EVq2PZD+UxbtJWtGXmM6d2ahyf0o22Tem6PJQFEwS5SQ0rKK3jx8128uCaVsHp1\neP66gYyPbqezdHGcgl2kBny/7wemx8azI/MYEwd2YO74KJo3DHV7LAlQCnaRalRYWs7Tq3bwxte7\nadu4Hm/cEsNFvVTaJdVLwS5STb5JPcyMuAT25RRyw/AIpo/tRZhKu6QGKNhFHJZXVMZjK7bx3oZ0\nIls04L2pwxnepYXbY0kQUbCLOOiT5EzmLEkgO7+EOy7owh/H9KBeHZV2Sc1SsIs44PCxEh5YlsTy\n+IP0ahvGqzfFEN1RpV3iDgW7iBestSzZsp95HyRTWFLBPZf04I4LulInRKVd4h4Fu0gVHcgtYvbi\nBD7fns3AiKY8MTma7irtEh+gYBc5Sx6PZeF3+5i/YhseC/dfEcVN50SqtEt8hoJd5CzsPlzA9Nh4\nvtudw7ndWvLYpH6EN1dpl/gWBbvIGSiv8PD6V7t55pMd1K1diyeujuaawR1VByA+ScEuchrJB44y\nPTaehP15XBLVhocn9KV1Y5V2ie9yLNiNMSHARmC/tXa8U+uKuKWkvIIXPkvlpTW7aNqgDi9eP4jL\n+rbVWbr4PCfP2O8GtgGNHVxTxBWb9h4v7UrNOsakytKuZirtEj/hSLAbYzoClwOPAH9yYk0RNxSU\nlPPUqu289c0e2jepz1tThjCqZ2u3xxI5K06dsT8LTAN0E6/4rbU7s5kZl0DGD0XcOLwT0y/rRaO6\n+mco8T9ef9caY8YDWdbaTcaYUT+z3VRgKkBERIS3uxVxTF5hGY+sSOb9jRl0admQ9+84h6Gdm7s9\nlkiVOXE6MhK40hgzDqgHNDbG/MNae8OJG1lrFwALAGJiYqwD+xXx2keJB5m7NImcglJ+Paord4/u\nrtIu8XteB7u1diYwE6DyjP2ek0NdxNdk5Rdz/9IkViYeIqpdY968ZQh9OzRxeywRR+gCogQVay2x\nm/fz0PJkisoquPfSnkw9v4tKuySgOBrs1to1wBon1xRxSsYPhcxanMiXO7KJ6dSM+ZOj6da6kdtj\niThOZ+wS8Dweyzvr9vDEx9sBmHdlH24c3olaKu2SAKVgl4CWmnWMGbHxbNz7A+d1P17a1bGZSrsk\nsCnYJSCVVXhY8GUaz326k3q1a/Hk1dFcrdIuCRIKdgk4ifvzmB4bT9KBo4zr15YHruxD6zCVdknw\nULBLwCguq+Cvn+7klS/TaN4wlJdvGMTYvu3cHkukxinYJSBs2JPD9Nh40rILuGZwR+ZcHkWTBnXc\nHkvEFQp28WvHSsp54qMU3lm3lw5N6/P3Xw3lvO6t3B5LxFUKdvFbX+zIZlZcAgfyipgyMpJ7LulJ\nQ5V2iSjYxf/kFpby4PJk4jbvp1vrRiy6cwSDOzVzeywRn6FgF79hrWVl4iHuW5pIbmEZv7+oG7+7\nqBt1a6u0S+RECnbxC1lHi5m7NJGPkzLp26Ex79w6jKj2eliXyE9RsItPs9by700ZPLw8mZJyDzMu\n68Vt53amtkq7RE5JwS4+Kz2nkJlxCXyVepihkc2ZP7kfXVqptEvkdBTs4nMq/lPa9dF2ahl4aEJf\nrh8aodIukTOkYBefkpqVz7RF8Wzel8uonq14ZGI/OjSt7/ZYIn5FwS4+oazCwytf7OKvn6bSsG4I\nz/5iAFcNaK/SLpEqULCL6xIy8rh30VZSDuVzeXQ75l3Zh5aN6ro9lojfUrCLa4rLKvjL6h28tnY3\nLRqG8sqNg7m0T1u3xxLxewp2ccX6tCPMiEtg9+ECrh0SzsxxvWlSX6VdIk5QsEuNyi8uY/7KFBau\n30d48/osvG0YI7u1dHsskYCiYJca83lKFrMWJ3DoaDG3juzMPZf2oEGovgVFnKbfVVLtcgpKefCD\nJJZsOUD31o2I/fUIBkWotEukuijYpdpYa1kef5AHliWRV1TGXaO789sLu6q0S6SaKdilWmQeLWbO\nkkQ+Sc4kumMTFt4+jF5tVdolUhO8DnZjTDjwDtAW8AALrLXPebuu+CdrLf/akM4jK7ZRWu5h9rje\nTBkZqdIukRrkxBl7OfBna+1mY0wYsMkY84m1NtmBtcWP7D1SwMy4BL7ZdYRhnZvz+ORoIls2dHss\nkaDjdbBbaw8CBys/zjfGbAM6AAr2IFHhsbz59W6eWrWdOrVq8cjEvlw3RKVdIm5x9Bq7MSYSGAis\nd3Jd8V3bD+UzPTaeLem5jO7Vmocn9qVdE5V2ibjJsWA3xjQCYoE/WGuP/sTPTwWmAkRERDi1W3FJ\nabmHl9bs4oXPdxJWrw7PXTuAK/urtEvEFzgS7MaYOhwP9YXW2rif2sZauwBYABATE2Od2K+4Y2t6\nLtMWxbM9M5+rBrTnvvFRtFBpl4jPcOKuGAO8Dmyz1j7j/Ujiq4pK/1PalUbrsHq8dlMMY6LauD2W\niJzEiTP2kcCNQIIxZkvl52ZZa1c4sLb4iHW7jjAjLp69Rwq5bmgEM8f1onE9lXaJ+CIn7or5CtCF\n1QB1tLiMx1ak8M/v9tGpRQPevX0YI7qqtEvEl+mdp3JKq5Mzmb0kgez8Em4/rzN/urgn9UNVByDi\n6xTs8iNHjpUw74Nklm09QM82YSy4MYb+4U3dHktEzpCCXf7LWsuyrQd4YFkSx0rK+eOYHvx6VFdC\na6sOQMSfKNgFgIN5RcxZnMinKVn0D2/Kk1dH06NNmNtjiUgVKNiDnMdj+eeGfTy2IoVyj4c5l/dm\nysjOhKgOQMRvKdiD2J7DBcyIi+fbtBxGdG3B/EnRRLRo4PZYIuIlBXsQKq/w8MbXu3l61Q5CQ2ox\nf1I/fjEkXHUAIgFCwR5kUg4dZfqieLZm5DGmdxsentCXtk3quT2WiDhIwR4kSsorePHzXby4JpXG\n9erw/HUDGR/dTmfpIgFIwR4Evt/3A9Nj49mReYwJA9pz3xV9aN4w1O2xRKSaKNgDWGFpOU+v2sEb\nX++mbeN6vHnLEC7s1drtsUSkminYA9TXqYeZERdPek4RNw7vxLSxPQlTaZdIUFCwB5i8ojIe/XAb\n/9qYTueWDfnX1OEM69LC7bFEpAYp2APIx0mHmLskkcPHSrjjgi78cUwP6tVRaZdIsFGwB4Ds/BIe\nWJbEhwkH6dU2jNdvHkK/jk3cHktEXKJg92PWWhZ/v58HlydTWFLBPZf04I4LulInRKVdIsFMwe6n\n9ucWMXtxAmu2ZzMooilPXB1Nt9Yq7RIRBbvf8XgsC9fvZf7KFDwW7hsfxc0jIlXaJSL/pWD3I2nZ\nx5gRm8B3e3I4r3tLHp3Yj/DmKu0Skf+lYPcD5RUeFqxN49nVO6lXuxZPXB3NNYM7qg5ARH6Sgt3H\nJR3IY3psPIn7jzK2T1senNCH1mEq7RKRU1Ow+6jisgqe/2wnL3+RRrMGobx0/SAu69fO7bFExA8o\n2H3Qpr05TFsUz67sAiYP6sjc8b1p2kClXSJyZhTsPqSgpJwnP97O2+v20L5Jfd6+dSgX9Gjl9lgi\n4mccCXZjzFjgOSAEeM1aO9+JdYPJ2p3ZzIxLYH9uETcN78S9Y3vRqK7+3BWRs+d1chhjQoC/ARcD\nGcAGY8wya22yt2sHg7zCMh76MJlFmzLo0qoh799xDkMim7s9loj4MSdOCYcCqdbaNABjzHvAVYCC\n/TQ+SjzI3KVJ5BSU8ptRXblrdHeVdomI15wI9g5A+gmvM4BhDqwbsLLyi7l/aRIrEw/Rp31j3rxl\nCH07qLRLRJzhRLD/1Ltk7I82MmYqMBUgIiLCgd36H2stsZv389DyZIrKKpg2tie3n9dFpV0i4ign\ngj0DCD/hdUfgwMkbWWsXAAsAYmJifhT8gS49p5BZixNYu/MwMZ2a8fjV0XRt1cjtsUQkADkR7BuA\n7saYzsB+4Frglw6sGxA8Hsvfv93L4x+lYIAHr+rDDcM6UUulXSJSTbwOdmttuTHmd8DHHL/d8Q1r\nbZLXkwWA1KxjzIiNZ+PeHzi/RysendiXjs1U2iUi1cuRG6WttSuAFU6sFQjKKjws+DKN51bvpH5o\nCE9f059JgzqotEtEaoTeAeOwxP15TFsUT/LBo4zr15Z5V/alVVhdt8cSkSCiYHdIcVkFz326kwVf\nptG8YSgv3zCYsX3buj2WiAQhBbsDNuzJYfqieNIOF3DN4I7MuTyKJg3quD2WiAQpBbsXjpWU88RH\nKbyzbi8dm9Xn778aynndVdolIu5SsFfRmu1ZzF6cyIG8IqaMjOSeS3rSUKVdIuIDlERn6YeCUh5a\nnkzc9/vp1roRi+4cweBOzdweS0TkvxTsZ8hay4qEQ9y/LJHcwjJ+d2E3fj+6G3Vrq7RLRHyLgv0M\nZB0tZs6SRFYlZ9KvQxPeuXUYUe0buz2WiMhPUrD/DGst/96YwUMfJlNa7mHGZb247dzO1FZpl4j4\nMAX7KaTnFDIzLoGvUg8ztHNz5k/qRxeVdomIH1Cwn6TCY3n7mz08+fF2QmoZHp7Ql18OjVBpl4j4\nDQX7CXZm5jM9Np7N+3IZ1bMVj07sR/um9d0eS0TkrCjYgdJyD698sYvnP0ulYd0Qnv3FAK4a0F6l\nXSLil4I+2BMy8rh30VZSDuVzRf/23H9FFC0bqbRLRPxX0AZ7cVkFf1m9g1e/TKNVWF1evSmGi6Pa\nuD2WiIjXgjLY16cdYUZcArsPF3Dd0HBmjutN43oq7RKRwBBUwZ5fXMb8lSksXL+PiOYNePe2YYzo\n1tLtsUREHBU0wf5ZSiazFyeSebSY287tzJ8v6Un9UNUBiEjgCfhgzyko5cEPkliy5QDdWzfixV+P\nYGCESrtEJHAFbLBba/kg/iAPLEviaFEZd4/uzm8u7KrSLhEJeAEZ7Ifyjpd2rd6WSf+OTXj89mH0\naqvSLhEJDgEV7NZa3tuQzqMfbqPM42H2uN7cem5nQlQHICJBJGCCfe+RAmbEJrAu7QjDuzRn/qRo\nIls2dHssEZEa5/fBXuGxvPn1bp5atZ06tWrx6MR+XDskXKVdIhK0vAp2Y8yTwBVAKbALmGKtzXVi\nsDOx/VA+02Lj2Zqey+herXl4Yl/aNVFpl4gEN2+fGPEJ0NdaGw3sAGZ6P9LplZZ7eHb1DsY/v5b0\nnEL+et1AXrs5RqEuIoKXZ+zW2lUnvPwWuNq7cU5vS3ou0xfFsz0znysrS7taqLRLROS/nLzGfivw\nLwfX+5HnP93JX1bvoHVYPV6/OYbRvVXaJSJystMGuzFmNdD2J35qtrV2aeU2s4FyYOHPrDMVmAoQ\nERFRpWEjWjTg2qERzLisl0q7REROwVhrvVvAmJuBO4HR1trCM/k1MTExduPGjV7tV0Qk2BhjNllr\nY063nbd3xYwFpgMXnGmoi4hI9fL2rpgXgDDgE2PMFmPMyw7MJCIiXvD2rphuTg0iIiLO8PaMXURE\nfIyCXUQkwCjYRUQCjIJdRCTAKNhFRAKM129QqtJOjckG9lbxl7cEDjs4jpt0LL4nUI4DdCy+yptj\n6WStbXW6jVwJdm8YYzaeyTuv/IGOxfcEynGAjsVX1cSx6FKMiEiAUbCLiAQYfwz2BW4P4CAdi+8J\nlOMAHYuvqvZj8btr7CIi8vP88YxdRER+hl8GuzHmIWNMfGWj5CpjTHu3Z6oqY8yTxpiUyuNZbIxp\n6vZMVWGMucYYk2SM8Rhj/PLuBWPMWGPMdmNMqjFmhtvzVJUx5g1jTJYxJtHtWbxhjAk3xnxujNlW\n+b11t9szVZUxpp4x5jtjzNbKY5lXrfvzx0sxxpjG1tqjlR/fBURZa+90eawqMcZcAnxmrS03xjwO\nYK2d7vJYZ80Y0xvwAK8A91hr/epJKsaYEI4/kP1iIAPYAFxnrU12dbAqMMacDxwD3rHW9nV7nqoy\nxrQD2llrNxtjwoBNwAQ//ZoYoKG19pgxpg7wFXC3tfbb6tifX56x/yfUKzUE/O9Pp0rW2lXW2vLK\nl98CHd2cp6qstdustdvdnsMLQ4FUa22atbYUeA+4yuWZqsRa+yWQ4/Yc3rLWHrTWbq78OB/YBnRw\nd6qqsccdq3xZp/K/asstvwx2AGPMI8aYdOB64D6353HIrcBKt4cIUh2A9BNeZ+CnIRKIjDGRwEBg\nvbuTVJ0xJsQYswXIAj6x1lbbsfhssBtjVhtjEn/iv6sArLWzrbXhHH+A9u/cnfbnne5YKrc57QPB\n3XYmx+HHzE98zm//JhhIjDGNgFjgDyf9bd2vWGsrrLUDOP638qHGmGq7TObVE5Sqk7V2zBlu+i7w\nIXB/NY7jldMdS+UDwcdz/IHgPhsmZ/E18UcZQPgJrzsCB1yaRSpVXo+OBRZaa+PcnscJ1tpcY8wa\nYCxQLf/A7bNn7D/HGNP9hJdXAiluzeKtEx4IfqUeCO6qDUB3Y0xnY0wocC2wzOWZglrlPzi+Dmyz\n1j7j9jzeMMa0+s8db8aY+sAYqjG3/PWumFigJ8fvwtgL3Gmt3e/uVFVjjEkF6gJHKj/1rT/e4WOM\nmQg8D7QCcoEt1tpL3Z3q7BhjxgHPAiHAG9baR1weqUqMMf8ERnG8RTATuN9a+7qrQ1WBMeZcYC2Q\nwPHf6wCzrLUr3Juqaowx0cDbHP/eqgW8b619sNr254/BLiIip+aXl2JEROTUFOwiIgFGwS4iEmAU\n7CIiAUbBLiISYBTsIiIBRsEuIhJgFOwiIgHm/wE2ODkym2Hv2QAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f132c146a20>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "a = 4\n",
        "b = 2\n",
        "\n",
        "def f(x):\n",
        "    return a + b*x\n",
        "\n",
        "plot_grid = np.linspace(-3, 3, 100)\n",
        "\n",
        "pt.plot(plot_grid, f(plot_grid))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "But suppose we don't know $a$ and $b$, but instead all we have is a few noisy measurements (i.e. here with random numbers added):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f132c121a90>]"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXNxAgQAh7JmGvECIjDMFZHIioDG211lnF\nTu2wDBkiLqrWuurAbX9YS0mQISiCWxEBCzeDACGMhBECSAjZyf3+/khskaJA7knOHe/n4+Ej3MPh\nnM8x4c3JybnvY6y1iIhI8AhzewAREXGWgl1EJMgo2EVEgoyCXUQkyCjYRUSCjIJdRCTIKNhFRIKM\ngl1EJMgo2EVEgkx9N3baunVr26VLFzd2LSISsDZs2HDQWtvmVOu5EuxdunRh/fr1buxaRCRgGWN2\nnc56uhQjIhJkFOwiIkFGwS4iEmQU7CIiQUbBLiISZE472I0xrxhjDhhjUo9b1tIY874xZlv1xxa1\nM6aIiJyuMzljfw0YfcKyqcBqa21PYHX1axERcdFpB7u19hPg8AmLrwJer/7168A4h+YSEQkq3xSW\ncd/SNI6WlNf6vny9xt7OWrsPoPpj2+9b0RgzyRiz3hizPi8vz8fdiogEBmst73j2cfFfP+bva3bx\nVdaJ58fOq7N3nlpr5wHzABITE/UEbREJerlHS5j5dior03NJiI7i7z8fRt8OzWp9v74Ge64xpoO1\ndp8xpgNwwImhREQCmbWWBeuzeeCdzZRVeLlnTB9uHdmV+vXq5kZEX4N9CXATMLf642KfJxIRCWC7\nDxUxbZGHzzMPMbRrS/48MYGurZvU6QynHezGmH8AFwCtjTE5wL1UBfoCY8zPgd3ANbUxpIiIv6v0\nWl77YiePvbeFemGGB8bF89OhsYSFmTqf5bSD3Vp73ff81iiHZhERCUhbcwuYvNDDxuwjXNi7DQ+O\n70/H5hGuzeNKba+ISDAoq/Dy/MfbeeaDTJo0rMcTPxnAVQM6Ykzdn6UfT8EuIlIDm7KPMCXJQ8b+\nAq44qyP3XhFH66YN3R4LULCLiJyR4rJKnli1lRc/zaJNZENevDGRi+PauT3WdyjYRURO05rth5iW\n7GHnoSKuGxrDtDF9adYo3O2x/oeCXUTkFApKynl4RQZvrt1NbMvGvHnbMEb0aO32WN9LwS4i8gM+\nyMjlnuRUDhSUcNs5XfnjJb2JaFDP7bF+kIJdROQkDh0rZc6ydBZv3EvvdpE8f8NgBsQ0d3us06Jg\nFxE5jrWWJZv2ct/SdApKyrlrVE9+fWEPGtQPnOcSKdhFRKrtzy9hxtsprNp8gLNimvPIxAR6t490\ne6wzpmAXkZBnreWtddk89M5myr1eZlzel1tGdqWeC3UATlCwi0hI23mwkGnJKazJOsTZ3Voxd2J/\nOreq29IupynYRSQkVXotr3y2g7+8v4XwsDAentCfa4fEuF4H4AQFu4iEnC37C5i8cBObcvK5qG87\nHhgXT/uoRm6P5RgFu4iEjLIKL3/7MJNnP8qkWaNwnrpuIFckdAiKs/TjKdhFJCRszD7C5IWb2Jp7\njPEDOzFzbBwtmzRwe6xaoWAXkaBWVFbB4yu38srnO2jXrBGv3JzIj/r4V2mX0xTsIhK0vsg8yNTk\nFHYfLuL6YbFMvawPkX5Y2uU0BbuIBJ384nIeXr6Zt9Zl06VVY96aNJzh3Vq5PVadUbCLSFB5Pz2X\nGW+nkFdQyh3ndeP3F/eiUbh/l3Y5TcEuIkHh4LFSZi9JY5lnH33aR/LijYkkRAdGaZfTFOwiEtCs\ntSzeuJf7lqZRWFrJHy/uxS8u6E54vcAp7XKagl1EAtbeI8XMeDuVDzIOMDC2qrSrZ7vAK+1ymoJd\nRAKO12t586vdzF2RQaXXMnNsHDeP6BKwpV1OU7CLSEDZcbCQqUke1u44zDk9WvPwhP7EtGzs9lh+\nxZFgN8b8HrgNsEAKcIu1tsSJbYuIAFRUenn5sx08/v5WGtQP45GJCVyTGB10dQBO8DnYjTGdgDuB\nOGttsTFmAXAt8Jqv2xYRAdi87yiTF3pI2ZPPJXHtuH9cPO2aBU9pl9OcuhRTH4gwxpQDjYG9Dm1X\nREJYaUUlz3yQyXMfbad543D+9tNBjOnfXmfpp+BzsFtr9xhjHgN2A8XASmvtyhPXM8ZMAiYBxMbG\n+rpbEQlyG3Z9w5QkD5kHjjFhUCdmXh5HiyAt7XKazzd6GmNaAFcBXYGOQBNjzM9OXM9aO89am2it\nTWzTpo2vuxWRIFVYWsF9S9O4+vkvKCqt4LVbhvD4jwco1M+AE5diLgJ2WGvzAIwxycAI4P8c2LaI\nhJDPth1karKHnG+KufHszkwe3YemDXXz3ply4v/YbmC4MaYxVZdiRgHrHdiuiISI/KJyHlyezoL1\nOXRr3YQFd5zN0K4t3R4rYDlxjX2tMWYh8DVQAfwbmOfrdkUkNLybup+Zi1M5XFjGry7ozp2jeoZc\naZfTHPkex1p7L3CvE9sSkdBwoKCE2UvSWJ6yn7gOzXj15iHEd4pye6ygoItXIlKnrLUkf72HOcvS\nKS6v5E+X9mbSed1CurTLaQp2Eakze44Uc09yCh9vzWNw5xb8eWICPdo2dXusoKNgF5Fa5/Va5q/d\nxdwVGVjgviv7ccPwzoSptKtWKNhFpFZtzzvG1CQP63Z+w7k9W/PQeJV21TYFu4jUivJKLy9+msUT\nq7YREV6Px645i4mDOqkOoA4o2EXEcal78pmS5CFt71FG92vPnHH9aBup0q66omAXEceUlFfy9Afb\neP7jLFo0bsBz1w/isv4d3B4r5CjYRcQR63ceZnKSh6y8QiYOimbm2L40b6x+Fzco2EXEJ4WlFTzy\nbgZvfLmLjlERvHHrUM7rpaI/NynYRaTGPt6axz3JKezNL+ams7vwp0t700SlXa7TZ0BEztiRojLu\nX7aZpK9z6N6mCf+642wSu6i0y18o2EXkjKxI2cfMxWkcKSrj1xd257c/UmmXv1Gwi8hpOXC0hFmL\n03g3bT/xnZrx+q1D6NdRpV3+SMEuIj/IWsu/NuTwwLJ0Siq8TBndh9vP7Up9lXb5LQW7iHyv7MNF\n3LMohU+3HWRol5bMndifbm1U2uXvFOwi/s6zAFbPgfwciIqGUbMg4ce1ustKr+WNNTt59L0tGOD+\nq/px/TCVdgUKBbuIP/MsgKV3Qnlx1ev87KrXUGvhnnmggMkLPXy9+wjn92rDQxP606l5RK3sS2qH\ngl3En62e899Q/1Z5cdVyh4O9vNLLCx9v56nVmTRuWI/Hf3wW4weqtCsQKdhF/Fl+zpktr6GUnHz+\ntHATGfsLuDyhA7Ov6EebyIaO7kPqjoJdxJ9FRVddfjnZcgeUlFfyxKptvPhpFq2aNOCFGwZzab/2\njmxb3KNgF/Fno2Z99xo7QHhE1XIffbXjMFOSPOw4WMhPEmO45/K+REWE+7xdcZ+CXcSffXsd3cG7\nYgpKyvnzuxn835e7iWkZwfzbhjGyR2uHBhZ/oGAX8XcJP3bsB6UfbjnA9OQU9h0t4efndOWPl/Si\ncQPFQLBx5DNqjGkOvATEAxa41Vq7xolti4jvDheWcf+ydBb9ew892zYl6ZcjGBTbwu2xpJY49U/1\nk8C71tqrjTENAD2pVsQPWGtZ5tnH7CVp5BeXc+eonvz6wu40rK/SrmDmc7AbY5oB5wE3A1hry4Ay\nX7crIr7JPVrC9EWprNqcS0J0FPNvH0af9s3cHkvqgBNn7N2APOBVY8xZwAbgLmttoQPbFpEzZK3l\nn+uyeXD5ZsoqvEwf05dbRnZRaVcIceIzXR8YBDxnrR0IFAJTT1zJGDPJGLPeGLM+Ly/Pgd2KyIl2\nHyri+pfWMjU5hbgOzXjvd+dx+3ndaifUPQvgr/Ewu3nVR88C5/chNeLEGXsOkGOtXVv9eiEnCXZr\n7TxgHkBiYqJ1YL8iUq3Sa3nti5089t4W6oUZHhrfn2uHxNReaZcLHTZy+nwOdmvtfmNMtjGmt7V2\nCzAKSPd9NBE5HVtzq0q7NmYf4Ud92vLg+Hg6RNVyaVcddtjImXPqrpjfAvOr74jJAm5xaLsi8j3K\nKrw899F2nvlwG00b1ufJawdw5Vkd66a0q446bKRmHAl2a+1GINGJbYnIqW3KPsKUJA8Z+wu44qyO\nzL4ijlZN67C0q5Y7bMQ3esuZSAApLqvkr6u28tKnWbSNbMRLNyZyUVy7uh+kFjtsxHcKdpEAsWb7\nIaYle9h5qIjrhsYwbUxfmjVyqbSrFjpsxDkKdhE/d7SknLkrMnhz7W46t2rMm7cPY0R3PyjtcrDD\nRpylYBfxY6s35zJ9USoHCkq4/dyu/OHi3kQ0UB2A/DAFu4gfOnSslPuWprNk0156t4vk+RsGMyCm\nudtjSYBQsIv4EWstSzbt5b6l6RSUlPO7i3ryqwt60KC+6gDk9CnYRfzEvvxiZixKZXXGAQbENOeR\nqxPo1S7S7bEkACnYRVzm9Vr+sW43Dy/PoMLrZcblfbllZFfq1VYdgAQ9BbuIi3YeLGRqsocvsw4z\nonsr5k5IILaVHmcgvlGwi7ig0mt5+bMs/rJyKw3qhfHwhKrSrjqpA5Cgp2AXqWNb9hcweeEmNuXk\nc1Hftjwwrj/toxq5PZYEEQW7SB0prajk2Q+38+xHmUQ2Cufp6wYyNqGDztLFcQp2kTrw793fMCXJ\nw9bcY4wf2ImZY+No2aSB22NJkFKwi9SiorIK/rJyK698voP2zRrxys2J/KiPC6VdElIU7CK15IvM\ng0xNTmH34SJ+NjyWKaP7EOlWaZeEFAW7iMPyi8t5ePlm3lqXTZdWjXlr0nCGd2vl9lgSQhTsIg56\nPz2XGW+nkFdQyh3nd+P3F/WiUbhKu6RuKdhFHHDwWCmzl6SxzLOPPu0jefHGRBKiVdol7lCwi/jA\nWsvbG/dw39J0ikorufuSXtxxfnfC66m0S9yjYBepob1Hipm+KIUPt+QxMLY5j0xMoKdKu8QPKNhF\nzpDXa5n/1W7mLt+M18K9V8Rx49ldVNolfkPBLnIGdhwsZEqSh692HOacHq15eEJ/YlqqtEv8i4Jd\n5DRUVHp5+bMdPP7+VhrWD+ORqxO4ZnC06gDELynYRU4hfe9RpiR5SNmTzyVx7XhgXDxtmwVQaZdn\nAayeA/k5EBUNo2bpIdRBzrFgN8bUA9YDe6y1Y53arohbSisqeeaDTJ77aDvNG4fz7PWDuCy+fWCd\npXsWwNI7oby46nV+dtVrULgHMSfP2O8CNgPNHNymiCs27Koq7co8cIwJ1aVdLQKxtGv1nP+G+rfK\ni6uWK9iDliPBboyJBi4HHgT+4MQ2RdxQWFrBYyu38NoXO+kYFcFrtwzhgt5t3R6r5vJzzmy5BAWn\nztifACYDuolXAtan2/KYlpxCzjfF3DC8M1Mu60PThgH+Y6io6KrLLydbLkHL57fHGWPGAgestRtO\nsd4kY8x6Y8z6vLw8X3cr4pj8onImL9zEDS9/RYN6YSy442zuHxcf+KEOVT8oDY/47rLwiKrlErSc\n+ModCVxpjBkDNAKaGWP+z1r7s+NXstbOA+YBJCYmWgf2K+Kzd1P3MXNxGocLy/jlBd25a1TP4Crt\n+vY6uu6KCSnGWucy1hhzAXD3qe6KSUxMtOvXr3dsvyJn6kBBCfcuTmNF6n7iOjTjkasTiO8U5fZY\nIj/IGLPBWpt4qvWC4HtNkdNnrSXp6z3cvyyd4vJK/nRpbyad102lXRJUHA12a+1HwEdOblPEKTnf\nFHHPolQ+2ZpHYucWzJ2YQI+2Td0eS8RxOmOXoOf1Wt5Ys5NH3tsCwH1X9uOG4Z0JU2mXBCkFuwS1\nzAPHmJrkYf2ubzi3Z1VpV3QLlXZJcFOwS1Aqr/Qy75Msnly9jUb1w3j06gSuVmmXhAgFuwSd1D35\nTEnykLb3KGP6t2f2lf1oGxlApV0iPlKwS9AoKa/kqdXbeOGTLFo2acDzPxvE6PgObo8lUucU7BIU\n1u08zJQkD1l5hVwzOJoZl8cR1Tjc7bFEXKFgl4B2rLSCR97N4I01u+jUPIK//3wo5/Zs4/ZYIq5S\nsEvA+nhrHvckp7A3v5hbRnbh7kt60yQY+l1EfKS/BRJwjhSVMWdZOslf76FH26Ys/MUIBndu4fZY\nIn5DwS4Bw1rLitT9zFqcypGicn77ox785kc9aFg/iEq7RBygYJeAcOBoCTMXp/JeWi7xnZrxxq3D\niOuoh3WJnIyCXfyatZZ/bcjhgWXplFZ4mXpZH247pyv1Vdol8r0U7OK3sg8XMS05hc8yDzK0S0vm\nTuxPtzYq7RI5FQW7+J3Kb0u73t1CmIH7x8Vz/dBYlXaJnCYFu/iVzAMFTF7o4evdR7igdxseHN+f\nTs0jTv0HReQ/FOziF8orvbzw8XaeWp1Jk4b1eOInA7hqQEeVdonUgIJdXJeSk8+fFm4iY38Blyd0\n4L4r+9G6aUO3xxIJWAp2cU1JeSV/XbWVlz7dQasmDXjhhsFc2q+922OJBDwFu7hibdYhpiansONg\nIdcOiWHamL5ERai0S8QJCnapUwUl5cxdkcH8tbuJaRnB/NuGMbJHa7fHEgkqCnapMx9mHOCeRSns\nP1rCrSO7cvelvWjcQF+CIk7T3yqpdYcLy5izNI23N+6lZ9umJP1yBINiVdolUlsU7FJrrLUs8+xj\n9pI08ovLuXNUT359YXeVdonUMgW71IrcoyXMeDuV99NzSYiOYv7tw+jTXqVdInXB52A3xsQAbwDt\nAS8wz1r7pK/blcBkreWf67J5cPlmyiq8TB/Tl1tGdlFpl0gdcuKMvQL4o7X2a2NMJLDBGPO+tTbd\ngW1LANl1qJBpySl8sf0Qw7q25M8TE+jSuonbY4mEHJ+D3Vq7D9hX/esCY8xmoBOgYA8RlV7Lq5/v\n4LGVWwgPC+PB8fFcN0SlXSJucfQauzGmCzAQWOvkdsV/bdlfwJQkDxuzjzCqT1seGB9PhyiVdom4\nybFgN8Y0BZKA31lrj57k9ycBkwBiY2Od2q24pKzCy3MfbeeZD7cR2SicJ68dwJVnqbRLxB84EuzG\nmHCqQn2+tTb5ZOtYa+cB8wASExOtE/sVd2zKPsLkhR625BZw1YCOzBobRyuVdon4DSfuijHAy8Bm\na+3jvo8k/qq47NvSrizaRjbipRsTuSiundtjicgJnDhjHwncAKQYYzZWL7vHWrvcgW2Ln1iz/RBT\nkz3sOlTEdUNjmTamD80aqbRLxB85cVfMZ4AurAapoyXlPLw8g398tZvOrRrz5u3DGNFdpV0i/kzv\nPJXvtSo9l+lvp5BXUMrt53blDxf3JqKB6gBE/J2CXf7HoWOl3Lc0nSWb9tK7XSTzbkjkrJjmbo8l\nIqdJwS7/Ya1lyaa9zF6SxrHSCn5/US9+eUF3GtRXHYBIIFGwCwD78ouZsSiV1RkHOCumOY9enUCv\ndpFujyUiNaBgD3Fer+Uf63bz8PIMKrxeZlzel1tGdqWe6gBEApaCPYTtPFjI1GQPX2YdZkT3Vsyd\nkEBsq8ZujyUiPlKwh6CKSi+vfL6Dv6zcSoN6Ycyd0J+fDIlRHYBIkFCwh5iM/UeZstDDppx8Lurb\njgfGxdM+qpHbY4mIgxTsIaK0opJnP9zOsx9l0qxROE9fN5CxCR10li4ShBTsIeDfu79hSpKHrbnH\nGDegI7Ou6EfLJg3cHktEaomCPYgVlVXwl5VbeeXzHbRv1ohXbx7ChX3auj2WiNQyBXuQ+jzzIFOT\nPWQfLuaG4Z2ZPLo3kSrtEgkJCvYgk19czkPvbOaf67Pp2roJ/5w0nGHdWrk9lojUIQV7EHkvbT8z\n307l4LFS7ji/G7+/qBeNwlXaJRJqFOxBIK+glNlL0ngnZR992kfy8k1D6B8d5fZYIuISBXsAs9ay\n6N97mLMsnaLSSu6+pBd3nN+d8Hoq7RIJZQr2ALXnSDHTF6Xw0ZY8BsU255GrE+jRVqVdIqJgDzhe\nr2X+2l3MXZGB18KssXHcNKKLSrtE5D8U7AEkK+8YU5NS+GrnYc7t2ZqHxvcnpqVKu0TkuxTsAaCi\n0su8T7N4YtU2GtUP45GrE7hmcLTqAETkpBTsfi5tbz5Tkjyk7jnK6H7tmTOuH20jVdolIt9Pwe6n\nSsorefqDbTz/cRYtGjfguesHcVn/Dm6PJSIBQMHuhzbsOszkhR625xUycVA0M8f2pXljlXaJyOlR\nsPuRwtIKHn1vC6+v2UnHqAhev3Uo5/dq4/ZYIhJgHAl2Y8xo4EmgHvCStXauE9sNJZ9uy2Nacgp7\njhRz4/DO/Gl0H5o21L+7InLmfE4OY0w94G/AxUAOsM4Ys8Ram+7rtkNBflE597+TzsINOXRr04QF\nd5zNkC4t3R5LRAKYE6eEQ4FMa20WgDHmLeAqQMF+Cu+m7mPm4jQOF5bxqwu6c+eonirtEhGfORHs\nnYDs417nAMMc2G5w8CyA1XMgPweiomHULA50vZJ7F6exInU//To249WbhxDfSaVdIuIMJ4L9ZO+S\nsf+zkjGTgEkAsbGxDuw2AHgWwNI7oby46nV+NhWLf8vjlSmsrhjB5NG9uf3cbirtEhFHOZEoOUDM\nca+jgb0nrmStnWetTbTWJrZpEyJ3eqye899Qr1a/soQ/hL3FirvO5VcX9FCoi4jjnEiVdUBPY0xX\nY0wD4FpgiQPbDXz5OSdd3MabR/c2Tet4GBEJFT4Hu7W2AvgN8B6wGVhgrU3zdbvBoLxpx5MuN1HR\ndTyJiIQSR64DWGuXW2t7WWu7W2sfdGKbgay80svfPsxkypHxFHPCO0bDI2DULHcGE5GQoAu8Dkvd\nk89Vz3zOo+9toaTvBMrHPAlRMYCp+njFU5DwY7fHFJEgprc2OqSkvJInV29j3idZtGzSgOd/NpjR\n8e2BwTD0p26PJyIhRMHugHU7DzNloYesg4VcMziaGZfHEdU43O2xRCREKdh9cKy0gkfezeCNNbuI\nbhHB338+lHN7hsitnCLitxTsNfTRlgNMX5TK3vxibhnZhbsv6U0TlXaJiB9QEp2hbwrLuH9ZOsn/\n3kOPtk1Z+IsRDO7cwu2xRET+Q8F+mqy1LE/Zz71LUjlSVM5vLuzBb0f1oGF9lXaJiH9RsJ+GA0dL\nmPF2KivTc+nfKYo3bh1GXMdmbo8lInJSCvYfYK3lX+tzuP+ddMoqvEy9rA+3ndOV+up3ERE/pmD/\nHtmHi5iWnMJnmQcZ2rUlcyf0p5v6XUQkACjYT1Dptbz+xU4efW8L9cIMD4yL56dDYwkLO1k7sYiI\n/1GwH2dbbgFTkjxE5yzj04iFtKrMw6yJhsazVAMgIgFDwQ6UVXh54ePtPP1BJhPDv+D+iFeoX1lS\n9Zv52VUPywCFu4gEhJD/KWBKTj5XPvMZf3l/K5fGt+eBZsn/DfVvlRdXPTRDRCQAhOwZe0l5JX9d\ntZUXP8miTWRDXrwxkYvj2sHsPSf/A9/z0AwREX8TksG+NusQU5NT2HGwkOuGxjBtTF+aNaou7YqK\nrrr8ciI9HENEAkRIBXtBSTlzV2Qwf+1uYls25s3bhjGiR+vvrjRq1ncfQA16OIaIBJSQCfYPMnKZ\nviiV3KMl3HZOV/54SW8iGpykDuDbH5CunlN1+SUquirU9YNTEQkQQR/shwvLmLM0jbc37qVn26Y8\n+8sRDIw9RWlXwo8V5CISsII22K21LPXsY/aSNI4Wl3PXqJ786sLuKu0SkaAXlMG+P7+qtGvV5lzO\nio7iz7cPo097lXaJSGgIqmC31vLWumweemcz5V4v08f05dZzulJPdQAiEkKCJth3HSpkalIKa7IO\nMbxbS+ZOSKBL6yZujyUiUucCPtgrvZZXP9/BYyu3EB4WxkPj+3PtkBiVdolIyPIp2I0xjwJXAGXA\nduAWa+0RJwY7HVv2FzA5ycOm7COM6tOWB8bH0yEqoq52LyLil3ztinkfiLfWJgBbgWm+j3RqZRVe\nnli1lbFPf0r24SKeum4gL92UqFAXEcHHM3Zr7crjXn4JXO3bOKe2MfsIUxZ62JJbwJVndeTeK+Jo\n1bRhbe9WRCRgOHmN/Vbgnw5u7388vXobf121lbaRjXj5pkRG9W1Xm7sTEQlIpwx2Y8wqoP1Jfmu6\ntXZx9TrTgQpg/g9sZxIwCSA2NrZGw8a2asy1Q2OZelmf/5Z2iYjIdxhrrW8bMOYm4BfAKGtt0en8\nmcTERLt+/Xqf9isiEmqMMRustYmnWs/Xu2JGA1OA80831EVEpHb5elfMM0Ak8L4xZqMx5nkHZhIR\nER/4eldMD6cGERERZ4T8M09FRIKNgl1EJMgo2EVEgoyCXUQkyCjYRUSCjM9vUKrRTo3JA3bV8I+3\nBg46OI6bdCz+J1iOA3Qs/sqXY+lsrW1zqpVcCXZfGGPWn847rwKBjsX/BMtxgI7FX9XFsehSjIhI\nkFGwi4gEmUAM9nluD+AgHYv/CZbjAB2Lv6r1Ywm4a+wiIvLDAvGMXUREfkBABrsx5n5jjKe6UXKl\nMaaj2zPVlDHmUWNMRvXxLDLGNHd7ppowxlxjjEkzxniNMQF594IxZrQxZosxJtMYM9XteWrKGPOK\nMeaAMSbV7Vl8YYyJMcZ8aIzZXP21dZfbM9WUMaaRMeYrY8ym6mO5r1b3F4iXYowxzay1R6t/fScQ\nZ639hctj1Ygx5hLgA2tthTHmzwDW2ikuj3XGjDF9AS/wAnC3tTagnqRijKlH1QPZLwZygHXAddba\ndFcHqwFjzHnAMeANa2282/PUlDGmA9DBWvu1MSYS2ACMC9DPiQGaWGuPGWPCgc+Au6y1X9bG/gLy\njP3bUK/WBAi8f52qWWtXWmsrql9+CUS7OU9NWWs3W2u3uD2HD4YCmdbaLGttGfAWcJXLM9WItfYT\n4LDbc/jKWrvPWvt19a8LgM1AJ3enqhlb5Vj1y/Dq/2ottwIy2AGMMQ8aY7KB64FZbs/jkFuBFW4P\nEaI6AdnHvc4hQEMkGBljugADgbXuTlJzxph6xpiNwAHgfWttrR2L3wa7MWaVMSb1JP9dBWCtnW6t\njaHqAdq/cXfaH3aqY6le55QPBHfb6RxHADMnWRaw3wkGE2NMUyAJ+N0J360HFGttpbV2AFXflQ81\nxtTaZTJhP4dpAAABS0lEQVSfnqBUm6y1F53mqm8C7wD31uI4PjnVsVQ/EHwsVQ8E99swOYPPSSDK\nAWKOex0N7HVpFqlWfT06CZhvrU12ex4nWGuPGGM+AkYDtfIDbr89Y/8hxpiex728EshwaxZfHfdA\n8Cv1QHBXrQN6GmO6GmMaANcCS1yeKaRV/8DxZWCztfZxt+fxhTGmzbd3vBljIoCLqMXcCtS7YpKA\n3lTdhbEL+IW1do+7U9WMMSYTaAgcql70ZSDe4WOMGQ88DbQBjgAbrbWXujvVmTHGjAGeAOoBr1hr\nH3R5pBoxxvwDuICqFsFc4F5r7cuuDlUDxphzgE+BFKr+rgPcY61d7t5UNWOMSQBep+prKwxYYK2d\nU2v7C8RgFxGR7xeQl2JEROT7KdhFRIKMgl1EJMgo2EVEgoyCXUQkyCjYRUSCjIJdRCTIKNhFRILM\n/wPOntV2zWfLkAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f1323d6d668>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "npts = 5\n",
        "\n",
        "np.random.seed(22)\n",
        "points = np.linspace(-2, 2, npts) + np.random.randn(npts)\n",
        "values = f(points) + 0.3*np.random.randn(npts)*f(points)\n",
        "\n",
        "pt.plot(plot_grid, f(plot_grid))\n",
        "pt.plot(points, values, \"o\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What's the system of equations for $a$ and $b$?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---------------\n",
        "Now build the system matrix $M$--a Vandermonde matrix:\n",
        "\n",
        "(We'll call it $M$ because $V$ conventionally used for the result of the SVD)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.        , -2.09194992],\n",
              "       [ 1.        , -2.46335065],\n",
              "       [ 1.        ,  1.08179168],\n",
              "       [ 1.        ,  0.76067483],\n",
              "       [ 1.        ,  1.50887086]])"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M = np.array([\n",
        "    1+0*points,\n",
        "    points,\n",
        "    ]).T\n",
        "M"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What's the right-hand side vector?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-------------\n",
        "Now solve the least-squares system:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U, sigma, VT = npla.svd(M)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "(5, 5)\n",
            "(2,)\n",
            "(2, 2)\n"
          ]
        }
      ],
      "source": [
        "print(U.shape)\n",
        "print(sigma.shape)\n",
        "print(VT.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Determine $x$, from $U\\Sigma V^T \\mathbf x = \\mathbf b$, where $\\mathbf b$ is `values`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = VT.T@((U.T@values)[:2] / sigma)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Recover the computed $a$, $b$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a_c, b_c = x"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f1323ae1e10>"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdcleX/x/HXxRJREQcuEME9yBXuvTWVMss0zVQUV+XI\nbZnmyNLMNLe2yb6WlThCXLgXKuUkFypOHCDIhuv3B+avciHnHA7n8Hk+HjzgnHPf1/W5Ed7eXOe+\nr0tprRFCCGE9bMxdgBBCCOOSYBdCCCsjwS6EEFZGgl0IIayMBLsQQlgZCXYhhLAyEuxCCGFlJNiF\nEMLKSLALIYSVsTNHp0WLFtWenp7m6FoIISzWoUOHbmqtXZ+2nVmC3dPTk9DQUHN0LYQQFkspdSEz\n28lQjBBCWBkJdiGEsDIS7EIIYWXMMsb+KCkpKURGRpKYmGjuUszC0dERd3d37O3tzV2KEMLC5Zhg\nj4yMpECBAnh6eqKUMnc52Uprza1bt4iMjMTLy8vc5QghLFymh2KUUl8qpW4opY7947nCSqlNSqnT\n9z8XymohiYmJFClSJNeFOoBSiiJFiuTav1aEEMb1LGPsXwPt//PcOGCL1roCsOX+4yzLjaH+t9x8\n7EII48p0sGutdwC3//P0i8A397/+BnjJSHUJIYRVuRV/i+FBw4lJjDF5X4ZeFVNca30V4P7nYo/b\nUCnlr5QKVUqFRkVFGdit8UVHR7Nw4UJzlyGEsDJaa346/hNVF1ZlwcEF7Liww+R9ZtvljlrrpVpr\nH621j6vrU++IzXaPC/a0tDQzVCOEsAZXYq/w8qqX6fZzNzwKenDI/xCdK3U2eb+GXhVzXSlVUmt9\nVSlVErhhjKLMYdy4cZw9e5aaNWtib2+Po6MjhQoV4tSpUwQHB9OpUyeOHct433j27NnExcUxefJk\nzp49y9ChQ4mKisLJyYlly5ZRuXJlMx+NEMKctNZ8eeRL3g1+l6S0JGa1mcXw+sOxs8meCxEN7SUQ\neBOYef/zGoMrAoYHDSfsWpgxmnqgZomazG0/97Gvz5w5k2PHjhEWFkZISAgdO3bk2LFjeHl5ERER\n8dj9/P39Wbx4MRUqVGD//v0MGTKErVu3GrV2IYTlOHfnHP5r/dlyfgtNyzRleeflVChSIVtryHSw\nK6VWAs2BokqpSOADMgJ9lVLKD7gIvGqKIs2hbt26T72mPC4ujj179vDqq/9/2ElJSaYuTQiRA6Wl\npzH/wHwmbp2IrbJlUcdF+D/vj43K/hv8Mx3sWusej3mplZFqeeBJZ9bZJV++fA++trOzIz09/cHj\nv683T09Px8XFhbAw4/51IYSwLMdvHMcv0I/9l/fzQoUXWNxxMaULljZbPTJXzH0FChQgNjb2ka8V\nL16cGzducOvWLZKSkli3bh0Azs7OeHl58dNPPwEZ42p//PFHttUshDCv5LRkpm6fSu2ltTlz+wzf\nd/medT3WmTXUIQdNKWBuRYoUoVGjRnh7e5M3b16KFy/+4DV7e3smTZpE3bp1KVWq1L/eHA0ICGDw\n4MFMmzaNlJQUunfvTo0aNcxxCEKIbHTw8kH8Av04euMo3b2783n7zymW77FXfGcrpbXO9k59fHz0\nfxfaOHnyJFWqVMn2WnIS+R4IkfPFp8QzOWQyn+79lBL5S7Co4yJ8K/lmS99KqUNaa5+nbSdn7EII\nkUkhESEMWDuAM7fPMKD2AGa1mUVBx4LmLushEuxCCPEUd5PuMmbTGJYcWkLZQmXZ0nsLLb1amrus\nx5JgF0KIJ1j/13oGrhvI1birjKw/kqktp+Jk72Tusp5Igl0IIR4h6l4UwzcO54ejP+BdzJtfXvuF\num51zV1WpkiwCyHEP2it+fHYj7wT9A4xiTF80OwDJjSZgIOtg7lLyzQJdiGEuO/y3csMXj+YtX+t\npa5bXVb4rsC7mLe5y3pmcoPSY0yePJnZs2c/9vXffvuNEydOZGNFQghT0Vqz7NAyqi6syuZzm/m0\n7afs6bfHIkMdLDjYAwLA0xNsbDI+BwRkb/8S7EJYhzO3z9Dq21b4r/Pn+ZLPc3TwUUY2GImtja25\nS8syiwz2gADw94cLF0DrjM/+/oaH+/Tp06lYsSKNGzcmPDwcgGXLllGnTh1q1KhB165diY+PZ8+e\nPQQGBjJ69Ghq1qzJ2bNnH7mdECLnSktP49M9n1J9UXUOXT3E0k5L2dJ7C+UKlzN3aYbTWmf7x/PP\nP6//68SJEw899zhlymidEen//ihTJtNNPCQ0NFR7e3vre/fu6ZiYGF2uXDk9a9YsffPmzQfbTJw4\nUc+bN09rrfWbb76pf/rppwevPW67Z/Es3wMhRNYdvX5U11laRzMZ7bvSV0fGRJq7pEwBQnUmMtYi\n3zy9ePHZns+MnTt30qVLF5ycMq5P9fXNuEX42LFjvPfee0RHRxMXF0e7du0euX9mtxNCmE9yWjIz\nds5gxs4ZuDi6sLLrSl6r9prVLSZvkcHu4ZEx/PKo5w3xqH/cPn368Ntvv1GjRg2+/vprQkJCHrlv\nZrcTQpjHgcsH6LemH8ejjtOrei8+a/cZRZ2KmrxfrTV/rfuLsq3KYu9kb/L+wELH2KdPB6f/3Pjl\n5JTxfFY1bdqUX3/9lYSEBGJjY1m7di0AsbGxlCxZkpSUFAL+MYj/32l+H7edEMK87iXf492N79Jg\nRQNikmJY12Md33X5LltC/fqf1/mu9Xf86Psjh5cfNnl/f7PIM/aePTM+T5yYMfzi4ZER6n8/nxW1\na9fmtddeo0aNGhQrVow6deoAMHXqVOrVq4erqyv16tV7EObdu3dnwIABzJs3j59//vmx2wkhzGfr\n+a0MWDuAc3fOMej5QXzc5mOc8zibvN97UffY9v42Di87jKOLIx2+6IDPwKdOymg0Mm1vDiLfAyGM\nIzoxmtHBo1l+ZDnlC5dneeflNPNslm39f9PyGy7suECdoXVo/kFz8hbOa5R2ZdpeIUSuFBgeyOD1\ng7kWd43RDUczpfkU8tobJ1if5PSG07jXdydv4by0/bQtdo52uFZxNXm/j2KRY+xCCPFfN+7doPvP\n3Xnxxxcp6lSU/f3380mbT0we6lEnowjoEMAPHX9g//z9AJSsVdJsoQ5yxi6EsHBaa344+gPDgoYR\nmxzL1BZTGdtoLPa2pr0CJeF2AiFTQji44CAO+R1oO6ctdYfmjNkfJdiFEBbrUswlBq8fzPrT66nv\nXp8Vviuo6lo1W/oO7B9I+JpwavvXpsWHLcjnmi9b+s0MCXYhhMVJ1+ksPbSUMZvGkKbT+KzdZ7xd\n922Tz+9ydtNZXKu44uzuTKsZrWg+uTnFqxd/+o7ZTIJdCGFRTt86zYC1A9h+YTuty7ZmaaeleBXy\nMmmft07fIvjdYP5a+xf1htej/WftKVrZ9NfBZ5VR3jxVSo1QSh1XSh1TSq1USjkao93sNm/ePKpU\nqUJPQy6If4KIiAi8vS1zGlAhzC01PZVZu2dRfXF1wq6FscJ3BcG9gk0a6okxiQSPDmZhtYVEbIug\n1cxWtJ7Z2mT9GYvBZ+xKKTfgHaCq1jpBKbUK6A58bWjb2W3hwoVs3rwZd3d3c5cihPiHP6//Sb81\n/Th09RAvVX6JBS8soFSBUibvd+PIjYR9FUbNPjVpNaMV+UvkN3mfxmCsoRg7IK9SKgVwAq4Yqd1s\nM2jQIM6dO0eHDh3o06cPO3fu5Ny5czg5ObF06VKqV6/O5MmTyZ8/P6NGjQLA29ubdevWAdChQwca\nN27Mnj17cHNzY82aNeTNm5dDhw7Rr18/ANq2bWu24xPCEiWlJjFtxzRm7p5J4byFWfXKKl6p+opJ\nJ+2K2B5BgZIFKFKxCM0mNaPOkDqUet70/4kYk8HBrrW+rJSaDVwEEoBgrXXwf7dTSvkD/gAemZit\n6+vmXz/0XMVOFWk4qmGWXu8T0ueJ/S1evJigoCC2bdvGlClTqFWrFr/99htbt26ld+/ehIWFPXH/\n06dPs3LlSpYtW0a3bt1YvXo1vXr1om/fvsyfP59mzZoxevToJ7YhhPh/ey/txS/Qj5M3T9K7Rm/m\ntJ1DEaciJusvOiKaTaM3ceLnE9TsU5MXv3oRlzIuuJRxMVmfpmLwGLtSqhDwIuAFlALyKaV6/Xc7\nrfVSrbWP1trH1dV8F+5nxq5du3jjjTcAaNmyJbdu3SImJuaJ+3h5eVGzZk0Ann/+eSIiIoiJiSE6\nOppmzTJuZf67TSHE48UlxzE8aDiNvmxEXHIcv/f8nW9e+sZkoZ4Um8SWCVv4ovIXnN5wmuYfNueF\nhS+YpK/sYoyhmNbAea11FIBS6hegIfC9IY0+7Qzb0Nef5FHz5yilsLOzIz09/cFziYmJD77OkyfP\ng69tbW1JSEhAa2118zwLYUqbz21mwNoBRERHMLTOUD5q9REF8hQwaZ8hk0PYN2cf1XtVp9VHrXB2\nN/0kYaZmjKtiLgL1lVJOKiPFWgEnjdCu2TRt2vTB1LshISEULVoUZ2dnPD09OXw4Y+rNw4cPc/78\n+Se24+LiQsGCBdm1axeATOcrxGPcSbiD3xo/2nzXBgdbB3b02cEXL3xhslC/tOcS18KuAdB4bGP8\n9vrR5bsuVhHqYJwx9v1KqZ+Bw0AqcARYami75jR58mT69u1L9erVcXJy4ptvvgGga9eufPvtt1Sr\nVo169epRsWLFp7b11Vdf0a9fP5RS8uapEI/w68lfGbJhCFH3ohjfeDyTmk3C0c40V0zHXIph89jN\nHFt5jEovVqL7b93JVywf+YrlnLtGjUGm7c1B5HsgcpNrcdd4+/e3+fnEz9QsUZMVviuoXbK2SfpK\niU9h9ye72f3JbtDQYFQDGo9tjEN+B5P0Zyoyba8QIkfSWvPdn98xPGg48SnxzGg5g1ENR5l00q69\nc/ayfcp2qnWrRutPWlvklS7PQoJdCJFtLsZcZOC6gQSdCaJh6Yas8F1B5aKVTdLXldArpCWnUbph\naeoNq4dnc088Ghu4MLKFyFHBnpuvIjHHkJgQ2SVdp7M4dDFjN49Fa838DvMZUmcINsr4S0LEXoll\ny4Qt/PHNH3g29+TNbW+Sp0CeXBPqkIOC3dHRkVu3blGkSJFcF+5aa27duoWjo0VOsSPEE4XfDKf/\n2v7suriLNmXbsLTzUjxdPI3eT2piKnvn7GXnjJ2kp6TTaGwjmkxoYvR+LEGOCXZ3d3ciIyOJiooy\ndylm4ejoKHPUCKuSkpbCp3s/ZXLIZJzsnfj6xa/pXaO3yU7cDq84zNaJW6ncpTJtZrWhcLnCJunH\nEuSYYLe3t8fLy7RTbwohsseRq0fwC/TjyLUjvFzlZRa8sIAS+UsYvZ9rYdeIvxlP2dZlqd2/NsW8\ni+HZzNPo/VgaWfNUCGE0iamJTNwykTrL6nAl9go/v/ozq7utNnqo37txj7X+a1lSewmbxmwCwC6P\nnYT6fTnmjF0IYdl2X9yNX6Af4bfCebPGm8xpN4fCeY07HJKWnMb+efvZMXUHKfEp1B9en2aTmhm1\nD2sgwS6EMEhcchzjN49nwcEFeBT0YGOvjbQtZ5q7rI//dJxNozdR4YUKtJ3TlqKVcu4qRuYkwS6E\nyLKNZzbiv86fSzGXeKvuW8xoNYP8DsZdjOLG8RvcOXuHSr6V8O7ujbO7swy5PIUEuxDimd1OuM3I\njSP55o9vqFy0Mjv77qSRRyOj9hF/K56QD0IIXRyKi6cLFTpWwMbWRkI9EyTYhRDPZPWJ1QzdMJRb\nCbeY0HgC7zd736iTdqWlpBG6KJSQySEkxSThM9iH5lOaY2Mr13pklgS7ECJTrsZe5a3f3+KXk79Q\nu2RtgnoFUbNETaP3c27zOYKGBVG2dVnafdaOYt7FjN6HtZNgF0I8kdaar8O+ZmTwSBJSEpjZaibv\nNnwXOxvjxcfN8JtcC7uG92velG9fnj7b++DRxCPX3YVuLBLsQojHioiOwH+tP5vObaKJRxOW+y6n\nYpGnr0OQWYnRiWz/cDsH5h/AqagTlV+sjJ2jHWWaljFaH7mRDFoJkcMFBICnJ9jYZHzOjoW40tLT\nmLd/Ht4LvdkbuZcFLywgpE+I0UI9PTWd0MWhzK8wn31z91GjTw0Ghg3EzlHONY1BvotC5GABAeDv\nD/HxGY8vXMh4DNCzp2n6PBl1Er9AP/ZG7qV9+fYs6bQEj4LGnRnxSugV1g9eT5mmZWg3tx0la5U0\navu5XY5ZQUkI8TBPz4ww/68yZSAiwrh9paSl8MnuT/hwx4fkd8jP3HZz6VW9l9HGue+cu0PE9ghq\n9a0FZKw76t7AXcbRn4GsoCSEFbh48dmez6pDVw7RL7Aff17/k27VujGv/TyK5y9ulLaTYpPYOX0n\n+z7bh72TPVVeroJjQUdKNyxtlPbFwyTYhcjBPDwefcbuYaSRkYSUBKZsn8LsPbMplq8Yv772Ky9V\nfskobet0Tdg3YWwZv4V71+9R480atJrRCseCsu6AqUmwC5GDTZ/+7zF2ACenjOcNtfPCTvwC/Th9\n+zR+tfyY3XY2Lo7GWwv05qmbrO2/Frd6bvRY2wO3Om5Ga1s8mQS7EDnY32+QTpyYMfzi4ZER6oa8\ncXo36S7jNo9jUegivFy82PzGZlqVbWWUeqMvRHN6/WnqDKmDa1VX/Pb6UapOKRlHz2YS7ELkcD17\nGu8KmN9P/87AdQOJvBvJiPojmNpiKvkc8hncbvK9ZHZ/vJs9s/agbBSVu1SmQMkCuNWVs3RzMEqw\nK6VcgOWAN6CBflrrvcZoWwhhuJvxNxmxcQTf//k9VV2rssdvD/Xd6xvcrk7XHP3hKJvHbSb2cize\n3b1p/XFrCpQsYISqRVYZ64z9cyBIa/2KUsoBcDJSu0IIA2itWXV8FW///jZ3Eu8wqekkJjSZQB67\nPEZp/+7luwT2D6SYdzFe+d8reDQy7vXuImsMDnallDPQFOgDoLVOBpINbVcIYZgrsVcYvH4wgeGB\n+JTyYYvvFp4r/pzB7d69fJdjK4/RcFRDCpYuiN8eP0rULIGykXH0nMIYZ+xlgSjgK6VUDeAQMExr\nfc8IbQshnpHWmhVHVjAqeBRJaUnMbjObYfWHGTxpV0pCCntm72H3zN2kp6VT6cVKFKlQhJK15a7R\nnMYYc8XYAbWBRVrrWsA9YNx/N1JK+SulQpVSoVFRUUboVgjxX+funKP1d60ZsHYANUvU5OjgowbP\nxKi15viq4yyovICQSSGU71CeoSeHEnSgSLbPYSMyxxhn7JFApNZ6//3HP/OIYNdaLwWWQsaUAkbo\nVwhxX1p6GvMPzGfi1onYKluWdFpC/9r9sVGGn7sl3E5g7YC1uHi58NI3L+HZ3NMsc9iIzDPKXDFK\nqZ1Af611uFJqMpBPaz36cdvLXDFCGM/xG8fxC/Rj/+X9dKzQkcWdFuPu7G5Qm3HX4ji8/DBNJjRB\n2Siu/3kd12quD1Yxys45bMT/y+65Yt4GAu5fEXMO6GukdoUQj5GclszMXTOZtmMaznmcCXg5gB7e\nPQy6GSg1KZV9c/exc/pOUhNTqfBCBUrWLknx6v+eNya75rARWWOUYNdahwFP/V9ECGEcBy8fxC/Q\nj6M3jtLduzvz2s/DNZ9rltvTWhO+Jpzgd4O5c+4OFTtXpO2nbSlSocgjtzf1HDbCMLLQhhAWJD4l\nntHBo6m/oj63E24T2D2QlV1XGhTqACnxKawfvB47Rzt6BfeiR2CPx4Y6ZExr4PSfu1WMNYeNMJxM\nKSCEhQiJCGHA2gGcuX2GAbUHMKvNLAo6Fsxye/E349k/fz9N32uKQz4Hem/tTZEKRbCxe/r5ninm\nsBHGI8EuRA4XkxjD2M1jWXJoCeUKlWNr76208GqR5fbSUtI4uOAgIZNDSI5LxqulF57NPHGt8mxn\n/cacw0YYlwS7EDnYur/WMWjdIK7GXeXdBu/yYYsPcbLP+owdpzecZuPIjdwKv0W5tuVoO6ctxaoV\nM2LFIieQYBciB4q6F8WwoGGsPLYS72Le/PLaL9R1q2tQm2kpafz+zu8oG0WPtT2o0LGCTKdrpSTY\nhchBtNb8eOxH3gl6h5jEGCY3m8z4JuNxsHXIUnsJdxLY++leGo9vjEM+B3r+3hOXMi7YOtgauXKR\nk0iwC5FDRN6NZPD6waz7ax313OqxwncF1YpVy1Jb6anphC4JJWRSCInRibjVc6NS50pPvNJFWA8J\ndiHMLF2ns+zQMkZvGk1qeipz2s7hnXrvYGuTtbPqc5vPsXHERm4cu4FnC0/az23/0A1GwrpJsAth\nRmdun2HA2gGERITQ0qslyzovo2yhslluT6drNo/dTEp8Ct1Wd6Nyl8oyjp4LSbALYQZp6Wl8tu8z\n3t/2Pg62DizttJT+tftnKYQTYxLZ/cluGoxsgFMRJ7qt7kb+Evmxc5Rf79xK/uWFyGbHbhyj35p+\nHLxykM4VO7Oo4yLcnJ99bdD0tHTCvgpj68St3Iu6R7FqxXju9edw8XQxQdXCkkiwC5FNklKT+GjX\nR8zYOYOCjgX5seuPdKvWLUtn6Rd2XCBoWBDXwq5RulFpXt/wOqWeL2WCqoUlkmAXIhvsj9yPX6Af\nx6OO06t6Lz5r9xlFnYpmub3tU7aTcDuBrj92pVq3ajKOLv5Fgl0IE7qXfI/3t73P3H1zcXN2Y12P\ndXSs2PGZ20mOS2bXx7t43v95CpYuyEvfvkTeQnmxd7I3QdXC0kmwC2EiW89vZcDaAZy7c47BPoOZ\n2Xomznmcn6kNna7547s/2DJ+C3FX43B2c8ZnkA/Obs/WjshdJNiFMLLoxGhGB49m+ZHllC9cnpA3\nQ2jm2eyZ27m05xJBw4O4cvAKbnXdeO2X13Cvb9jKSCJ3kGAXwogCwwMZvH4w1+KuMabhGCY3n0xe\n+7xZamvfZ/uIvRxLl++68Nzrz6FsZBxdZI4EuxBGcOPeDd75/R3+d/x/VC9enTXd1+BT6tkWFUuJ\nT2HP7D1UfbUqrlVceWHBC9g72eOQP2vzxIjcS4JdCANorQk4GsCwoGHEJccxrcU0xjQag71t5t/U\n1Fpz7MdjbB67mbuX7mLnaIdrFVfyFctnwsqFNZNgFyKLLsVcYtD6QWw4vYH67vVZ4buCqq5Vn6mN\nK6FXCBoexKXdlyhRqwQvB7xMmSZlTFSxyC0k2IV4Ruk6nSWhSxizeQzpOp3P23/O0DpDszRp1+Hl\nh7l9+jadl3emZp+a2NjKMsTCcEprne2d+vj46NDQ0GzvVwhDnb51mv5r+7Pjwg5al23N0k5L8Srk\nlen9UxNT2fvZXrxaeOFe352EOwnY2NqQxzmPCasW1kIpdUhr/dQ3b+SMXYhMSE1P5bO9nzEpZBKO\ndo586fslfWr2yfQdn1prTv16iuBRwUSfj6bxhMa413cnb6GsXTEjxJPI331CPMUf1/6g/vL6jNk8\nhvbl23NiyAn61uqb6VC/FnaNb1t+y6quq3DI58Abm9+g1fRWJq76/wUEgKcn2NhkfA4IyLauhZkY\n7YxdKWULhAKXtdadjNWuEOaSlJrEtB3TmLl7JoXzFuanV3+ia5Wuzzwvy/FVx7l+9DodF3Wkdv/a\n2Nhl3/lUQAD4+0N8fMbjCxcyHgP07JltZYhsZrQxdqXUSMAHcH5asMsYu8jp9l7ai1+gHydvnuSN\n6m/wWbvPKOKUuWXl0pLT2D9/P8WrF6dcm3IkxyWTlpJmlmEXT8+MMP+vMmUgIiK7qxGGyuwYu1FO\nHZRS7kBHYLkx2hPCXOKS4xgeNJxGXzbiXso9fu/5O992+TZToa61JnxtOAu9F7Jp1CbCA8MBcMjv\nYLax9IsXn+15YR2MNRQzFxgDFDBSe0Jku01nN+G/zp+I6AiG+AxhZuuZFMiTuR/pG8dvsHHERs5t\nOkfRykV5fcPrVOhQwcQVP52Hx6PP2D08sr8WkX0MPmNXSnUCbmitDz1lO3+lVKhSKjQqKsrQboUw\nmjsJd/Bb40fb79viYOvAjj47WNBxQaZDXWvNmd/PcOXgFdp/3p5Bfw7KEaEOMH06ODn9+zknp4zn\nhfUyeIxdKfUR8AaQCjgCzsAvWutej9tHxthFTvHLyV8YumEoUfeiGN1wNB80/wBHO8en7peWkkbo\n4lAKlCxA1VeqkpacRtLdJJyKOj113+wWEAATJ2YMv3h4ZIS6vHFqmTI7xm7UG5SUUs2BUfLmqcjp\nrsVd460Nb7H65GpqlqjJCt8V1C5ZO1P7ntl4ho0jNnLz5E1q9K7BS9+8ZOJqhcggNygJ8Qhaa779\n41tGbBxBfEo801tOZ3TD0ZmatOvWX7cIfjeYv9b9RaFyhXjtt9eo5FspG6oW4tkYNdi11iFAiDHb\nFMJYLkRfYOC6gWw8u5FGpRux3Hc5lYtWzvT+F3df5MKOC7SZ1Ya6b9fFLo+cF4mcSX4yhdVL1+ks\nOLCA8VvGAzC/w3yG1BmCjXrytQPpaekcXn4YWwdbavWtRc03a1KxU0Xyucp0uiJnk2AXVu3UzVP0\nD+zP7ku7aVuuLUs7LaWMy9OnxT2/9TxBw4O4cfQGlV+qTK2+tVA2SkJdWAQJdmGVUtJSmL1nNlO2\nTyGvfV6+evEr3qzx5lOnA7hz7g7Bo4I59espCpYpyKs/vUqVrlWyqWohjEOCXVidI1eP4Bfox5Fr\nR3il6ivM7zCfEvlLZGrf60evczb4LC2mtaDByAbY5838SkhC5BQS7MJqJKYm8uH2D/lk9ye45nNl\ndbfVvFzl5Sfuo9M1Yd+EkXQ3ifrD6lPJtxLDzg+TIRdh0STYhVXYdXEX/QP7E34rnL41+/Jp208p\nlLfQE/e5uOsiQcODuHroKl4tvaj3Tj2UknF0Yfkk2IVFi02KZfyW8Sw4uIAyBcsQ3CuYNuXaPHGf\nmIsxbBqzieP/O04BtwJ0+b4Lz73+3DNPxytETiXBLizWxjMb8V/nz6WYSwyrN4xpLaeR3yH/U/eL\njogmPDCcppOa0mhMIxzyOWRDtUJkHwl2YXFuJ9xmxMYRfPvHt1QpWoXd/XbToHSDx26vteboD0e5\nc+4Ozd5vRpmmZRhxaQRORXLevC5CGIMEu7AYWmtWn1zN0A1DuZ1wm/eavMd7Td8jj93jF4K+fOAy\nQcOCiNyowUEUAAAUnUlEQVQXiVs9NxqPa4ytva2EurBqEuzCIlyNvcrQDUP59dSv1C5Zm+BewdQo\nUeOx28dejWXz2M38+d2f5C+RH98vfan5Zk2UjYyjC+snwS5yNK01X4d9zcjgkSSmJvJx648Z2WAk\ndjZP/tGNvxnPydUnaTSuEU0mNCFPgcef1QthbSTYRY51/s55/Nf5s/ncZpp4NGG573IqFqn4yG21\n1pz46QRXQq/Q5pM2FH+uOCMiR5htSTohzEmCXeQ4aelpLDiYMWmXjbJh4QsLGegz8LGTdl09fJWg\n4UFc3HmR4jWK03xyc+yd7CXURa4lwS5ylJNRJ/EL9GNv5F46lO/A4k6L8Sj46AU670XdY8v4LRz5\n8ghORZzouLgjtfvXxsbWKGu0C2GxJNhFjpCSlsInuz/hwx0fUsChAN93+Z7Xn3v9iTcNJcclc3zV\nceqPqE+z95vh6PL0Je2EyA0k2IXZHbpyiH6B/fjz+p90q9aN+R3mUyxfsYe201oTviacMxvP0GlR\nJwp5FWLExRES6EL8hwS7MJuElAQmh0zm072fUixfMX597Vdeqvzo9UOvH73OxhEbOb/lPK5VXUm4\nnUDewnkl1IV4BAl2YRY7Luygf2B/Tt8+Tf9a/ZnVdhYuji4PbZdwO4Gt723l0JJDOLo40mF+B3wG\n+WBjJ+PoQjyOBLvIVneT7jJ201gWH1qMl4sXm9/YTKuyrR67fXpaOsdXHafO0Do0n9ycvIXlShch\nnkaCXWSbDac3MHDdQC7fvczwesOZ1nIa+RweniL39IbTHP3hKF2+7UI+13wMOzeMPM5yg5EQmSXB\nLkzuZvxNhgcNJ+BoAFVdq/KT30/Ud6//0HZRJ6MIHhnMmaAzFK5QmNirsTi7OUuoC/GMJNiFyWit\nWXV8FW///jZ3Eu8wqekkJjSZ8NCkXUl3k9g2aRsHvjiAQ34H2n7alrpv1cXWwdZMlQth2STYhUlc\nib3CkPVDWBO+Bp9SPmzx3cJzxZ979MYKTvx8gtr9a9NiagtZwUgIAxkc7Eqp0sC3QAkgHViqtf7c\n0HaFZdJas+LICkYFjyIpLYnZbWYzrP6whybtOrf5HKGLQun6Y1fyFMjDW6fewiG/LHghhDEY44w9\nFXhXa31YKVUAOKSU2qS1PmGEtoUFOXv7LP7r/Nl6fivNyjRjue9yyhcu/69tbp+5TfCoYMLXhOPi\n6UJ0RDRFKhSRUBfCiAwOdq31VeDq/a9jlVInATdAgj2XSEtP4/P9n/Pe1vewt7VnccfFDHh+wL8m\n7UqJTyFkcgj75u7DLo8dLWe0pMGIBtg5ymigEMZm1N8qpZQnUAvYb8x2Rc517MYx+gf2Z//l/XSq\n2IlFHRfh7uz+0HY2djaErwmneq/qtJzekgIlC5ihWiFyB6MFu1IqP7AaGK61vvuI1/0BfwAPj0fP\n1icsR3JaMjN3zWTajmkUdCzIDy//QHfv7v+atOvCjgvs/mQ3r/zvFRzyOeB/2F8WjhYiGxgl2JVS\n9mSEeoDW+pdHbaO1XgosBfDx8dHG6FeYx8HLB+kX2I9jN47x+nOvM7fdXFzzuT54PToimk2jN3Hi\n5xM4l3bmztk7FK9eXEJdiGxijKtiFLACOKm1nmN4SSKnik+J54NtHzBn3xxK5i9JYPdAOlfq/OD1\n1KRUdkzdwZ7Ze1A2iuZTmtNwVEPsnezNWLUQuY8xztgbAW8AR5VSYfefm6C13mCEtkUOERIRQv/A\n/py9cxb/2v580uYTCjoW/Nc2tva2nAk6Q9WuVWn9cWuc3Z3NVK0QuZsxrorZBcjS71YqJjGGMZvG\nsPTwUsoVKsfW3ltp4dXiweuX9l5i+5TtvBzwMk5FnOi7sy/2eeUMXQhzkmvNxGOtDV/LoPWDuBZ3\njXcbvMuHLT7Eyd4JgJhLMWwZt4WjPxwlf8n83D59G6ciThLqQuQAEuziIVH3ohgWNIyVx1biXcyb\n3177jTpudQBIT01n54yd7Jq5C52uaTKxCY3HNZYbjITIQSTYxQNaa1YeW8k7v7/D3aS7TGk+hXGN\nx+Fg+/+hrWwVF3depGKnirT5pA0ung8vjiGEMC8JdgFA5N1IBq8fzLq/1lHXrS5f+n5JtWLVALgS\neoUt47fgu8KXgh4F6bG2h9wxKkQOJr+duVy6TmfZoWWM3jSa1PRU5rSdwzv13sHWxpbYq7FsnbiV\nsK/CyFcsH7fP3KagR0EJdSFyOPkNzcXO3D7DgLUDCIkIoaVXS5Z1XkbZQmXRWrNr5i52Tt9JalIq\nDUc3pOl7TWXBCyEshAR7LpSansrcfXN5f9v7ONg6sKzzMvxq+T2YDkApxdXDV/Fq5UXb2W0pXL6w\nmSsWQjwLCfZc5uj1o/gF+nHwykF8K/my8IWFuDm7ce2Pa2watYl2c9tRrFoxunzXBbs88uMhhCWS\n39xcIik1iY92fcSMnTNwcXThx64/0q1aN+JvxrNu0DoOLzuMYyFHos9HU6xaMQl1ISyY/PbmAvsj\n9+MX6MfxqOP0fK4nc9vPpahTUfZ9vo+QSSGkxKdQ9+26NPugGXkL5TV3uUIIA0mwW7F7yfd4f9v7\nzN03FzdnN9a/vp4XKrzw4PWbp25SulFp2s1pR9HKRc1YqRDCmCTYrdSWc1sYsHYA56PPM8RnCB+1\n/oiks0l83/57mn3QjNINStNhXgds7W3NXaoQwsgk2K1MdGI0o4JHseLICioUrsD2PtupU6AO20Zt\nI3RRKHkK5OFuZMY6KBLqQlgnCXYr8tup3xiyfgjX711nTMMxTG4+mRNfnWDeuHkkxSRR2782Lae2\nxKmok7lLFUKYkAS7Fbged523f3+bn078RPXi1VnbYy3Pl3oegOgL0ZSsVZJ2n7WjePXiZq5UCJEd\nJNgtmNaa7//8nuEbhxOXHMe0FtPwc/Vj28BtFBhUgIodK9JiSguUrfrXWqRCCOtmY+4CRNZcjLlI\nxx860vu33lQqUokD3Q/gs8aHZdWXEbE9gvioeABs7Gwk1IXIZeSM3cKk63QWhy5m7OaxpOt05rab\nS9PwpgQ3DCb+Zjy1+tWi5fSW5C+e39ylCiHMRILdgvx16y/6B/Zn58WdtCnbhiWdluBVyIu9e/ZS\ntHJR2s9tT8naJc1dphDCzCTYLUBqeiqz98xmcshk8trnZWntpTh/5cxd7sIbUG9YPeqPqC9DLkII\nQMbYc7ywa2HUW16P8VvG07l0Z76+8TXXX7nO2eCzJMcmA2BjK+PoQoj/J2fsOVRiaiJTt0/l490f\nU9SpKCucV3Bn3B3CroVR/Y3qtPqoFc5uzuYuUwiRA0mw50B7Lu3BL9CPUzdP8Wb1N5nTfg4XfrzA\nYc/DvPbba7jXczd3iUKIHEyCPQeJS45jwpYJfHHgCyqnV2b+ofnUKVSHwnkLU+jNQtR8sybKRoZc\nhBBPZpRgV0q1Bz4HbIHlWuuZxmg3N9l0dhP+6/y5cuMKo8+PpsCvBYghhtRmqQAS6EKITDM42JVS\ntsACoA0QCRxUSgVqrU8Y2nZucCfhDiODR/J12Ne0uNKCgYEDSbqWROXXKtP649a4lHExd4lCCAtj\njDP2usAZrfU5AKXUj8CLgAT7U/xy8heGbhhKVFwU4xuP5/XY19lzaA/tf2qPR2MPc5cnhLBQxgh2\nN+DSPx5HAvWM0K5VCAiAiRPh4kXw8IDp06HVi9d4a8NbBO8P5tU9r9K6YWt6tOqB1ppqvtVk2EUI\nYRBjBPujUkg/tJFS/oA/gIdH7jgbDQgAf3+Iz5i2hQsXoF//VBx/fZ/no2/z7u53sU23pWT7jLtF\nlVKP/m4KIcQzMEawRwKl//HYHbjy34201kuBpQA+Pj4PBb81mjjx/0P9b+6JF/H9tTIu6e5UebkK\nbWa1oVDZQuYpUAhhlYwR7AeBCkopL+Ay0B143QjtWryLF//5SAOKVGxJTM9H760v49XCy0yVCSGs\nmcFTCmitU4G3gI3ASWCV1vq4oe1ag5JuKeQjDl8CaU8QABcpw+8eAyXUhRAmY5Tr2LXWG4ANxmjL\nGqSkpfBJyCeULXSdJpElsSOVvdQHNE5OiukzZCBdCGE6MgmYkR25eoT249pz49UbtD5ahAI13Pit\n1BC2qDaUKaNYuhR69jR3lUIIayZTChhJYmoiU7ZNYdbeWVTSlWjr2pauK7tSvl15xpu7OCFEriLB\nbgRbw7ay4u0VJNxNoPeHvfm07ae4zHKRqXSFEGYhwW6A6Lhopr87HdtvbSmfXJ5SPUvh39lfbjAS\nQpiVBHsW/fzrz+wetBuXGy6k1Eqh75d98azpae6yhBBC3jx9VjfjbtL7194M2D4A7KDWl7WYemiq\nhLoQIseQM/ZMir8dz4rhKzgcepiV3Vcy1ncs7816D0d7R3OXJoQQ/yLB/hTpqelsm7+NkA9CsI21\nRTVRHOhzgFoetcxdmhBCPJIE+xNcP3qdFS+vIOVMCpe9LlNxXkWW9l6KnY1824QQOZck1COkp6Vz\n4e4FhuwYgsc9D24NvcX0ydOpVLSSuUsTQoinkmD/h6S7SWyftp0Dmw7wUdePsLW1xXeNLwN9BmKj\n5H1mIYRlkGAn4ww97OswgscFk3QziTCPGGw/CyfujhsfL1E4T5dpAIQQliPXB/vtM7dZ1W0V149c\n55LHJba1yUPkzs9JTsz41ly4kLFYBki4CyEsQ64dX0hPSwfgtD7NyTsn+bnrz9z79B7Jp+Y9CPW/\nxcdnLJohhBCWINedsSfHJbNr5i7C14YT8VEEs0NnU+KtEizqtAjfSr7YdHv0fv9eNEMIIXKuXBPs\nOl3zZ8CfbBm3hdgrsZx//jz/2/E//Br5MavNLAo6FgQyFpy+cOHh/XPJMq1CCCuQK4L9buRdVnVd\nxeUDl0mukMy3ft/iUN2B9Z3X09Kr5b+2nT793wtQAzg5ZTwvhBCWwKqDPT01HRs7G/IVy8dd7hLS\nI4QdFXcwvMFwpracipO900P7/P0G6cSJGcMvHh4ZoS5vnAohLIVVBntKQgp7Zu/h6PdH6bqrK6N3\njibghQCqulZlj+8e6rnXe+L+PXtKkAshLJdVBbvWmuOrjrN5zGZiLsaQv01+6n1Rj6t2V/mg2QeM\nbzyePHZ5zF2mEEKYlNUEe/zNeP7X5X9c3HWRIs8V4a/3/+IH2x+oU6oO633X81zx58xdohBCZAuL\nD/a0lDRs7W3JWzgveQrmoeB7BXkv73sk62Rmt5jN8PrDsbWxNXeZQgiRbSw22FMTU9k3dx8HFx5k\n4OGBXLW5yvJXl7MtYhvN3ZqzrPMyyhcub+4yhRAi21lcsGutOfXbKYLfDSb6fDQVO1dk4e6FTDo2\nCXtbe5Z0WkL/2v1l0i4hRK5lULArpWYBnYFk4CzQV2sdbYzCHiUpNokfX/yRiG0RuFZzpfGqxkyI\nmcCBsAN0qtiJRR0X4e7sbqruhRDCIhh6WrsJ8NZaVwf+AsYbXtLjOeR3IH/x/LSd35aoz6PocKoD\n5+6cY2XXlQR2D5RQF0IIDDxj11oH/+PhPuAVw8p5MqUUpWeXxi/Qj2O7jtHDuweft/8c13yupuxW\nCCEsijEHovsBvxuxvYdM2zGNBisacCfhDmt7rOWHrj9IqAshxH889YxdKbUZKPGIlyZqrdfc32Yi\nkAoEPKEdf8AfwCOLM2qVK1SOAbUH8HHrjx9M2iWEEOLflNbasAaUehMYBLTSWsc/bXsAHx8fHRoa\nalC/QgiR2yilDmmtfZ62naFXxbQHxgLNMhvqQgghTMvQMfYvgALAJqVUmFJqsRFqEkIIYQBDr4qR\nWzuFECKHkdszhRDCykiwCyGElZFgF0IIKyPBLoQQVkaCXQghrIzBNyhlqVOlooALWdy9KHDTiOWY\nkxxLzmMtxwFyLDmVIcdSRmv91HlUzBLshlBKhWbmzitLIMeS81jLcYAcS06VHcciQzFCCGFlJNiF\nEMLKWGKwLzV3AUYkx5LzWMtxgBxLTmXyY7G4MXYhhBBPZoln7EIIIZ7AIoNdKTVVKfXn/Rklg5VS\npcxdU1YppWYppU7dP55flVIu5q4pK5RSryqljiul0pVSFnn1glKqvVIqXCl1Rik1ztz1ZJVS6kul\n1A2l1DFz12IIpVRppdQ2pdTJ+z9bw8xdU1YppRyVUgeUUn/cP5YpJu3PEodilFLOWuu7979+B6iq\ntR5k5rKyRCnVFtiqtU5VSn0MoLUea+aynplSqgqQDiwBRmmtLWolFaWULRkLsrcBIoGDQA+t9Qmz\nFpYFSqmmQBzwrdba29z1ZJVSqiRQUmt9WClVADgEvGSh/yYKyKe1jlNK2QO7gGFa632m6M8iz9j/\nDvX78gGW97/TfVrrYK116v2H+wB3c9aTVVrrk1rrcHPXYYC6wBmt9TmtdTLwI/CimWvKEq31DuC2\nueswlNb6qtb68P2vY4GTgJt5q8oanSHu/kP7+x8myy2LDHYApdR0pdQloCcwydz1GInJFwQXj+UG\nXPrH40gsNESskVLKE6gF7DdvJVmnlLJVSoUBN4BNWmuTHUuODXal1Gal1LFHfLwIoLWeqLUuTcYC\n2m+Zt9one9qx3N/mqQuCm1tmjsOCqUc8Z7F/CVoTpVR+YDUw/D9/rVsUrXWa1romGX+V11VKmWyY\nzKAVlExJa906k5v+AKwHPjBhOQZ52rHcXxC8ExkLgufYMHmGfxNLFAmU/sdjd+CKmWoR990fj14N\nBGitfzF3PcagtY5WSoUA7QGTvMGdY8/Yn0QpVeEfD32BU+aqxVD/WBDcVxYEN6uDQAWllJdSygHo\nDgSauaZc7f4bjiuAk1rrOeauxxBKKde/r3hTSuUFWmPC3LLUq2JWA5XIuArjAjBIa33ZvFVljVLq\nDJAHuHX/qX2WeIWPUqoLMB9wBaKBMK11O/NW9WyUUi8AcwFb4Eut9XQzl5QlSqmVQHMyZhG8Dnyg\ntV5h1qKyQCnVGNgJHCXjdx1ggtZ6g/mqyhqlVHXgGzJ+tmyAVVrrD03WnyUGuxBCiMezyKEYIYQQ\njyfBLoQQVkaCXQghrIwEuxBCWBkJdiGEsDIS7EIIYWUk2IUQwspIsAshhJX5P72SnKqBAbSpAAAA\nAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f1323b11160>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "def f_c(x):\n",
        "    return a_c + b_c * x\n",
        "\n",
        "pt.plot(plot_grid, f(plot_grid), label=\"true\", color=\"green\")\n",
        "pt.plot(points, values, \"o\", label=\"data\", color=\"blue\")\n",
        "pt.plot(plot_grid, f_c(plot_grid), \"--\", label=\"found\",color=\"purple\",)\n",
        "\n",
        "if 0:\n",
        "    # show residual components\n",
        "    pt.vlines(points,\n",
        "                np.minimum(f_c(points), values),\n",
        "                np.maximum(f_c(points), values),\n",
        "                color=\"red\", lw=2)\n",
        "\n",
        "pt.legend(loc=\"best\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* If we enable 'show residual components above', what will appear?\n",
        "* Is it possible for the residual to involve the 'true model'?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "\n",
        "# - Little red bars connecting the data to the computed fit line.\n",
        "# - No. The residual is only computed against the known data, not the unknown true model."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What should happen if we change the number of data points?\n",
        "* What happens if there are lots of outliers?\n",
        "* What should happen if we don't add noise?\n",
        "* What about a bigger model?\n",
        "* What about different functions in the model?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "# - The fit gets better\n",
        "# - Least squares \"punishes\" outliers with the square of their distance. \n",
        "#   So outliers have a habit of severely wrecking a least squares fit.\n",
        "# - We should recover the exact data.\n",
        "# - No problem--just add coefficients.\n",
        "# - No problem--just use a different Vandermonde matrix."
      ]
    }
  ],
  "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.5.2+"
    },
    "widgets": {
      "state": {},
      "version": "1.1.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}