Applications_IP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications of Inner Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear regression\n",
    "\n",
    "A common example of the use of least squares is Linear Regression.  In this scenario, we suppose that two variables are related through a linear relationship, or we wish to approximate their relationship as linear.  Although data may show that the relationship is not exactly linear, or the data may be subject to experimental error, we might still hope to find the linear model that best fits the data.  This is commonly known as the \"line of best fit\", or the \"least squares line\".  \n",
    "\n",
    "For our example, suppose that we are studying the price of beer and nachos at various entertainment venues, and we collect data from a representative sample of venues.\n",
    "\n",
    "|Beer|Nachos|\n",
    "|----|------|\n",
    "|\\$4.5|\\$12 |\n",
    "|\\$4.75|\\$13|\n",
    "|\\$5.5|\\$14.8|\n",
    "|\\$5.9|\\$15.5|\n",
    "|\\$6.25|\\$15.9|\n",
    "|\\$7 |\\$17.25|\n",
    "\n",
    "We plot the data and notice that there is strong correlation between the two variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAApy0lEQVR4nO3de3RcZ3nv8e/jWxzHuI4xmDhO4pALNBfiRFkhnLQcLgVCShvaA22yAs0pFzct5JRCS8NJ1ykFukoPPRza0lWalkvSZoXSAiUN4VZaSrMOhdiqFcXKYFlEEZJVK4NlK0rGI4/nOX/sLb8TRbJGszV6t/b8Pmt5jWZmz55n/2asV+++vK+5OyIiIs1aEbsAERFZXtRwiIjIgqjhEBGRBVHDISIiC6KGQ0REFkQNh4iILEjbGg4zO8vM/sXMHjGzvWb26+njm8zsG2bWn96ePsfrrzWz75vZfjO7rV11iojIwli7ruMwszOAM9y928yeBewGXg/8d+CQu384bRBOd/ffnvHalcA+4FXAMPAgcKO797WlWBERaVrbehzuPuru3enPTwCPAGcC1wN3povdSdKYzHQVsN/df+DuU8Bn09eJiEhkq5biTcxsO3A58F1gi7uPQtK4mNlzZ3nJmcAPG+4PAy+eY907gZ0A69at67rggguY7kWtWLGCWq3G6tWrmZqaYu3atTz11FOsW7fuGbdHjx5lzZo1HDt2jFWrVlGv16fXT71eZ9WqVUxNTXHKKadw9OhRTj311Geso1KpsHbtWqampli1ahXHjx9nxYqkba7X66xcuZJarcaaNWvmXUe1WmXNmjXUajVWrFihbdI2aZu0TW3Zpr1795bd/Tmz/X6dS9sbDjNbD3weeJe7T5hZUy+b5bFZ96m5+x3AHQCXXnqp79mzp8VKi2VgYIDzzjsvdhnRKYdAWQTKIjCzxxb6mraeVWVmq0kajbvd/QvpwwfT4x/Tx0HGZnnpMHBWw/1twIH53m/VqiXpQC0LmzZtil1CLiiHQFkEyiKbdp5VZcAngUfc/aMNT90L3Jz+fDPwpVle/iBwgZmda2ZrgBvS153UdNdR4KmnnopdQi4oh0BZBMoim3b2OK4B3gy8wsz2pP+uAz4MvMrM+knOmvowgJltNbP7Ady9BrwT+BrJQfXPufveNtZaONP7TTudcgiURaAssmnbvh13f4DZj1UAvHKW5Q8A1zXcvx+4fyHv2eTxk46wevXq2CXkgnIIlEWgLLIpVLOrXVXB5ORk7BJyQTkEyiJQFtkUquHQwfFg8+bNsUvIBeUQKItAWWRTqIZjamoqdgm5MTw8HLuEXFAOgbIIlEU2hWo4TjnllNgl5Mb5558fu4RcUA6BsgiURTaFajiOHj0au4Tc2LtXJ6GBcmikLAJlkU3bBjmM4corr/Rdu3bFLkNEZNkws93ufuVCXlOoHocu6gl2794du4RcUA6BsgiURTbqcYiIdDD1ONTjOEF/USWUQ6AsAmWRjXocIiIdrON7HJVKJXYJudHb2xu7hFxQDoGyCJRFNoVqONauXRu7hNy48MILY5eQC8ohUBaBssimUA2HrhwPhoaGYpeQC8ohUBaBssimUA2HxqoKtmzZEruEXFAOgbIIlEU2hWo4jh8/HruE3Dh8+HDsEnJBOQTKIlAW2RSq4dDkLIGO9ySUQ6AsAmWRjX7TiojIghSq4dBEToEGfEwoh0BZBMoim7YdTTazTwGvA8bc/ZL0sb8FXpAushE47O47ZnntIPAEcByoNXtxysqVKzPXXRQbN26MXUIuKIdAWQTKIpt29jg+A1zb+IC7/6K770gbi88DXzjJ61+eLtv0FY21Wq2VOgvp4MGDsUvIBeUQKItAWWTTth6Hu3/bzLbP9pyZGfALwCsW8z3XrFmzmKtb1s4+++zYJeSCcgiURaAssol1jOMngYPu3j/H8w583cx2m9nOZleq/ZbBvn37YpeQC8ohiJ7F3XfD9u2wYkVye/fd0UqJnsUyF6vhuBG45yTPX+PuVwCvBd5hZi+da0Ez22lmu8xs1+HDhymXy4yOjjIyMsL4+DgDAwNUKhX6+vqo1+t0d3cDYXTM7u5u6vU6fX19VCoVBgYGGB8fZ2RkhNHRUcrlMoODg0xOTlIqlajVavT09DxtHdO3vb29VKtV+vv7mZiYYGhoiLGxMcbGxhgaGmJiYoL+/n6q1eqJsXJmrqOnp4darUapVGJycpLBwcGWtuniiy8u3Da18jldeumlhdumVj+ntWvXRtumQx//OPW3vx0eewzc4bHHqL/tbRz/67+O8t3bvHlzbj+npf7utaKto+Omu6rumz44nj62ChgButx93hnjzez9wKS7/9F8y1500UXeahBFs3v3brq6umKXEZ1yCKJmsX170mjMdM45MDi41NXoe9FguYyO+1NAaa5Gw8xOM7NnTf8MvBp4uJkVr1u3btGKXO70nyKhHIKoWcw1NlSkMaP0vcimbQ2Hmd0DfAd4gZkNm9lb06duYMZuKjPbamb3p3e3AA+YWQ/wPeDL7v7VZt5TEzkFmqgmoRyCqFnMdTA60kFqfS+y0UROItJ+d98NO3dC4x9369bBHXfATTfFq0uWza6qttFETsH0wblOpxyCqFncdFPSSJxzDpgltxEbDX0vsilUj6Orq8vVBU3UajUNM49yaKQsAmURdHyPo1qtxi4hN/bv3x+7hFxQDoGyCJRFNoVqOHTleLBt27bYJeSCcgiURaAssilUw6GxqoJyuRy7hFxQDoGyCJRFNoVqODSRU7B+/frYJeSCcgiURaAssinUb9oiHejP6tixY7FLyAXlECiLQFlkU6iGQwJNapVQDoGyCJRFNoVqOLSrKtDwKwnlECiLQFlkU6jftDo4Hhw6dCh2CbmgHIKTZpGjIc+Xgr4X2RTqCpjVq1fHLiE3tm7dGruEXFAOwZxZzBwO5LHHkvtQ2OFA9L3IplA9jqmpqdgl5Majjz4au4RcUA7BnFncfvvTx5CC5P7tt7e/qEj0vcimUEOOaJDDoF6v65gPyqHRnFmsWJFMrjSTGRT0ILK+F0HHDzmiYdWDPXv2xC4hFzo2h1mOWcyZRc6GPF8KHfu9WCSFajh0pkRwxRVXxC4hFzoyh+ljFg3TtLJzJ1c88sjsy//+7ydDnDdaty55vKA68nuxiArVcKjHEWiU4ERH5jDHMYvqb/7m7MvnbMjzpdCR34tFpGMcIkXTgccspHU6xqEexwnd3d2xS8iFjsxhjmMT1S1blriQ/OrI78UiKlTDoWMcwY4dO2KXkAsdmcMcxyxWf+QjcerJoY78XiyitjUcZvYpMxszs4cbHnu/mY2Y2Z7033VzvPZaM/u+me03s9uafc+jR48uRumFUCqVYpeQCx2ZwxzHLEo6IHxCR34vFlHbjnGY2UuBSeAud78kfez9wKS7/9FJXrcS2Ae8ChgGHgRudPe++d7ziiuucHVBE5VKhVNPPTV2GdEph0BZBMoiyNUxDnf/NtDKgDBXAfvd/QfuPgV8Fri+mRdqqOTgwIEDsUvIhULksEjjSBUii0WiLLKJcYzjnWb2ULor6/RZnj8T+GHD/eH0sVmZ2U4z22Vmu370ox9RLpcZHR1lZGSE8fFxBgYGqFQq9PX1Ua/XTxwUmz4dr7u7m3q9Tl9fH5VKhYGBAcbHxxkZGWF0dJRyuczg4CCTk5OUSiVqtRo9PT1PW8f0bW9vL9Vqlf7+fiYmJhgaGmJsbIyxsTGGhoaYmJigv7+farVKb2/vrOvo6emhVqtRKpWYnJxkcHCwpW06/fTTC7dNrXxOmzZtWtbb9PjHPoa//e3PuCbjBx/60IK36cknn8zFNuXhu+fuhdumVj+nVrT1dFwz2w7c17CragtQBhz4IHCGu79lxmveCLzG3d+W3n8zcJW73zrf+1122WU+/YF1upGREc48c872tmMs+xy2b08ai5nOOQcGBxe0qmWfxSJSFkGudlXNxt0Puvtxd68Df0myW2qmYeCshvvbAPUrF0jj8CSWfQ5DQwt7/CSWfRaLSFlks6TpmdkZDXd/Dnh4lsUeBC4ws3PNbA1wA3Bvk+vPXmRBaIj5xLLPYRHHkVr2WSwiZZFNO0/HvQf4DvACMxs2s7cC/9vMes3sIeDlwG+ky241s/sB3L0GvBP4GvAI8Dl339vMe2o6yGBycjJ2Cbmw7HNYxHGkln0Wi0hZZNO2iZzc/cZZHv7kHMseAK5ruH8/cP9C33PVqkLNS5XJ5s2bY5eQC8s+h+nxom6/Pdk9dfbZSaPRwjhSyz6LRaQssinUjj5N5BQMDw/HLiEXCpHDTTclB8Lr9eS2xcEHC5HFIlEW2RSq4TjllFNil5Ab559/fuwSckE5BMoiUBbZFKrh0JAjwd69TR0WKjzlECiLQFlko2HVRUQ6WO6v42g3DaseaKKahHIIlEWgLLJRj0NEpIOpx6Eexwn6iyqhHAJlESiLbNTjEBHpYB3f46hUKrFLyI3pkTU7nXIIlEWgLLIpVMOxdu3a2CXkxoUXXhi7hFxQDoGyCJRFNoVqOHTleDDUwuipRaQcAmURKItsCtVwaKyqYMuWLbFLyAXlECiLQFlkU6iG4/jx47FLyI3Dhw/HLiEX2p7DIk3ruhT0nQiURTaF+hNdk7MEOt6TaGsOd98NO3fC9Gng6bSuQMsDEbaTvhOBsshGv2lFWnX77aHRmPbUU8njIgVWqIZDEzkFGvAx0dYcFnFa16Wg70SgLLIpVMOxcuXK2CXkxsaNG2OXkAttzWERp3VdCvpOBMoim0I1HLVaLXYJuXHw4MHYJeRCW3NYxGldl4K+E4GyyKadc45/yszGzOzhhsc+YmYlM3vIzL5oZhvneO1gOjf5HjNregyRNWvWLELlxXB2Tv/qXWptzeGmm+COO+Ccc8Asub3jjlweGAd9Jxopi2za2eP4DHDtjMe+AVzi7i8C9gHvO8nrX+7uOxYyhor2Wwb79u2LXUIutD2HRZrWdSnoOxEoi2za1nC4+7eBQzMe+7q7T+9P+ndg22K+56mnnrqYq1vWLr300tgl5IJyCJRFoCyyiXmM4y3AV+Z4zoGvm9luM9t5spWY2U4z22Vmu4aHhymXy4yOjjIyMsL4+DgDAwNUKhX6+vqo1+t0d3cDYVjl7u5u6vU6fX19VCoVBgYGGB8fZ2RkhNHRUcrlMoODg0xOTlIqlajVavT09DxtHdO3vb29VKtV+vv7mZiYYGhoiLGxMcbGxhgaGmJiYoL+/n6q1eqJQdZmrqOnp4darUapVGJycpLBwcGWtmnXrl2F26ZWPqfdu3cXbpta/Zy+9a1vFW6bWv2cHnjggcJtU6ufUyvaOqy6mW0H7nP3S2Y8fjtwJfDzPksBZrbV3Q+Y2XNJdm/dmvZgTkrDqouILEzbhlU3sxVmdrmZ/bSZvcLMWh7oxcxuBl4H3DRbowHg7gfS2zHgi8BVzaxbEzkFmqgmkSmHZTScSDP0nQiURTYn7XGY2XnAbwM/BfQDjwNrgQuBp4C/AO5091mvvJvZ4zCza4GPAv/V3R+f4zWnASvc/Yn0528AH3D3r863MepxyKKZOZwIJKfa5visKZFWtKPH8SHgb4Dz3P017v4md39DelbUzwI/Brx5jmLuAb4DvMDMhs3srcDHgWcB30hPtf1EuuxWM7s/fekW4AEz6wG+B3y5mUYDNJFTo+l9rJ2u5RwKOJyIvhOBssimUFPHdnV1ubqgiVqtpmHmyZDDihUw2/8Ns+TU22VI34lAWQRLMnWsmV1jZteamS30te1WrVZjl5Ab+/fvj11CLrScwzIbTqQZ+k4EyiKbeRsOM7vLzC5Of76FZHfTrcAn21zbgunK8WDbtkW9RGbZajmHZTacSDP0nQiURTYnbTjM7ByS02afSH/+FZJG4x3A1WZ2tpltaH+ZzdFYVUG5XI5dQi60nMMyG06kGfpOBMoim/l28r2M5AD4tcApwEbg+cB5wMr0+T3AQ22qb0E0kVOwfv362CXkQqYcbrppWTcUM+k7ESiLbE7acLj7nWb2EuCNJI3GJ9z9rvQ02be6+11LUGPTinSgP6tjx47FLiEXlEOgLAJlkU0zpxX8GvAaYMrdv5k+9mzgt9pWlWSmSa0SyiFQFoGyyGbehiO9uO8rMx4bAnI3zZl2VQXrZh7Y7VDKIVAWgbLIZr6D4/9oZj9jZqtnee75ZvYBM3tL+8pbGB0cDw4dOjT/Qh1AOQTKIlAW2czX43g78G7gY2Z2iDDkyHZgAPi4u3+prRUuwOrVz2jfOtbWrVtjl5ALyiFQFoGyyOakPQ53/093f6+7n0dygPyDJA3JJe7+qjw1GgBTU1OxS8iNRx99NHYJuaAcAmURKItsCjXkiAY5DOr1uo75oBwaKYtAWQRLMuRInmlY9WDPnj2xS8gF5RAoi0BZZKMeh4hIB1uqQQ5PN7MXLfR1S0E9jkCjBCeUQ6AsAmWRTVM9DjP7Fsn8G6tIhhh5HPhXd393O4tbKPU4REQWpp09jh9z9wng54FPu3sXyayAuaIeRzA9OX2nUw6BsgiURTbNNhyrzOwM4BeA+9pYTya6GjTYsWNH7BJyQTkEyiJQFtk023B8APgaMODuD5rZ80nmIM+Vo0ePxi4hN0qlUuwSckE5BMoiUBbZNNVwuPvfufuL3P1X0/s/cPf/drLXmNmnzGzMzB5ueGyTmX3DzPrT29PneO21ZvZ9M9tvZrc1uzGayCk499xzY5eQC8ohUBaBssimqYbDzLaZ2RfThuCgmX3ezOabQuszJPN4NLoN+Ka7XwB8M70/871WAn8GvBa4CLjRzC5qpk4NlRwcOHAgdgm5oBwCZREoi2ya3VX1aeBeYCtwJvCP6WNzcvdvAzNHErseuDP9+U7g9bO89Cpgf9qrmQI+m75uXpp8Pti0aVPsEnJBOQTKIlAW2TTbcDzH3T/t7rX032eA57TwflvcfRQgvX3uLMucCfyw4f5w+tiszGynme0ys10HDx6kXC4zOjrKyMgI4+PjDAwMUKlU6Ovro16vnzibYvo87u7ubur1On19fVQqFQYGBhgfH2dkZITR0VHK5TKDg4NMTk5SKpWo1Wr09PQ8bR3Tt729vVSrVfr7+5mYmGBoaIixsTHGxsYYGhpiYmKC/v5+qtUqvb29s66jp6eHWq1GqVRicnKSwcHBlrbpySefLNw2tfI5PfXUU4XbplY/p4GBgcJtU6uf04EDBwq3Ta1+Tq1o9jqOfyLZ9XRP+tCNwC+7+yvned124D53vyS9f9jdNzY8P+7up894zRuB17j729L7bwaucvdb56vzsssu8+kPrNONjo5yxhlnxC4jOuUQKItAWQTtvI7jLSSn4v4nMAq8IX1soQ6mp/WS3o7NsswwcFbD/W1AUzskzayFkopJQ8wnlEOgLAJlkU2zZ1UNufvPuvtz3P257v56d3+shfe7F7g5/flmYLZh2R8ELjCzc81sDXBD+rp5aTrIYHJyMnYJuaAcAmURKItsmjqabGbPIZnUaXvja9x9zl6Hmd0DvAzYbGbDwO8CHwY+Z2ZvJZl69o3psluBv3L369y9ZmbvJLluZCXwKXff29TG6OD4CZs3b45dQi4oh0BZBMoim2Z/034J+Dfgn4DjzbzA3W+c46lnHBdx9wPAdQ337wfub7K2EzSRUzA8PMwLX/jC2GVEpxwCZREoi2yabTjWuftvt7WSRXDKKafELiE3zj///Ngl5IJyCJRFoCyyafbg+H1mdt38i8WlIUeCvXub2rtXeMohUBaBssjmpKfjmtkTgAMGnAZUgWPpfXf3DUtRZLM0rLqIyMIs+um47v4sd9+Q3q5w91Mb7ueq0QANq95IE9UklEOgLAJlkU2zFwD+HPDP7n4kvb8ReJm7/0Nbq1sg9ThERBamnRcA/u50owHg7odJTq/NFfU4Av1FlVAOgbIIlEU2zfY4HnL3F814rNfdL21bZS1Qj0NEZGHa2ePYZWYfNbPzzOz5ZvZ/gdw12ZVKJXYJuTE9QFqnUw6BsgiURTbNNhy3AlPA3wJ/BxwF3tGuolq1du3a2CXkxoUXXhi7hFxQDoGyCJRFNs2OVfWku9/m7le6e5e7v8/dn2x3cQulK8eDoaGh2CXkgnIIlEWgLLJZyFhV7wUuBk78We/ur2hTXS3RWFXBli1bYpeQC8ohUBaBssim2V1VdwMl4Fzg94BBklFsc+X48aaG0eoIhw8fjl1CLiiHQFkEyiKbZhuOZ7v7J4Fj7v6v6ai4V7exrpasWNHs5hSfjvcklEOgLAJlkU2z+3aOpbejZvbTJBMrbWtPSSIikmfNNhwfMrMfA94D/CmwAfiNtlXVIk3kFGjAx4RyCJRFoCyyaarhcPf70h+PAC9vXznZrFy5MnYJubFx48bYJeSCcgiURaAssjlpw2Fm/+skT7u7f3CR68mkVqvFLiE3Dh48yIYNuRuHcskph0BZBMoim/l6HLNdq3Ea8Fbg2UCuGo41a9bELiE3zj777Ngl5IJyCJRFoCyymW9Y9f8z/Q+4AzgV+GXgs8DzW3lDM3uBme1p+DdhZu+asczLzOxIwzIn6/mcoP2Wwb59+2KXkAvKIVAWgbLIZt5BDs1sE/Bu4CbgTuCP3X18Ud7cbCUwArzY3R9rePxlwG+6++sWsj4NcigisjCLPsihmX2E5EK/J4BL3f39i9VopF4JDDQ2GlloWPVAw0YnlEOgLAJlkc18V8y9B9gK/A5wIN2tNGFmT5jZxCK8/w3APXM89xIz6zGzr5jZxXOtwMx2mtkuM9t15MgRyuUyo6OjjIyMMD4+zsDAAJVKhb6+Pur1Ot3d3UD44nR3d1Ov1+nr66NSqTAwMMD4+DgjIyOMjo5SLpcZHBxkcnKSUqlErVajp6fnaeuYvu3t7aVardLf38/ExARDQ0OMjY0xNjbG0NAQExMT9Pf3U61WT4zOOXMdPT091Go1SqUSk5OTDA4OtrRNl19+eeG2qZXPqaurq3Db1OrntGHDhsJtU6uf09atWwu3Ta1+Tq1oaj6OdjCzNSQXEl7s7gdnPLcBqLv7pJldR7J77IL51nnRRRd5q0EUze7du+nq6opdRnTKIVAWgbIIWtlVFbPhuB54h7u/uollB4Er3b18suV0jENEZGHaOZFTO9zIHLupzOx5Zmbpz1eR1Pmj+VaoiZyC6a5yp1MOgbIIlEU2UXocZrYO+CHw/Om5zM3sFgB3/4SZvRP4VaAGVIB3u/v/m2+9XV1droNeiVqtpmHmUQ6NlEWgLIJl0+Nw96fc/dnTjUb62Cfc/RPpzx9394vd/TJ3v7qZRgOgWq22q+RlZ//+/bFLyAXlECiLQFlkU6hxyHXleLBtmwYvBuXQSFkEyiKbQjUcGqsqKJdPeh5Bx1AOgbIIlEU2hWo4NJFTsH79+tgl5IJyCJRFoCyyKdRv2linFufRsWPH5l+oAyiHQFkEyiKbQjUcEmhSq4RyCJRFoCyyKVTDoV1Vwbp162KXkAvKIVAWgbLIplC/aXVwPDh06FDsEnJBOQTKIlAW2RSq4Vi9enXsEnJj69atsUvIBeUQKItAWWRTqIZjamoqdgm58eijj8YuIReUQ6AsAmWRTbRBDttBgxwG9Xpdx3xQDo2URaAsgmUz5Ei7aCKnYM+ePbFLyAXlECiLQFlkox6HiEgHU49DPY4TNEpwQjkEyiJQFtmoxyEi0sHU41CP44TpOYY7nXIIlEWgLLJRj6OgdNZIQjkEyiJQFkHH9ziOHj0au4TcKJVKsUvIBeUQKItAWWRTqIZDEzkF5557buwSckE5BMoiUBbZRGk4zGzQzHrNbI+ZPWPfkiX+xMz2m9lDZnZFM+vVUMnBgQMHYpeQC8ohUBaBssgm5mztL3f3uabhei1wQfrvxcCfp7cnpcnng02bNsUuIReUQ6AsAmWRTV53VV0P3OWJfwc2mtkZ871IY+wHOsMsoRwCZREoi2xiNRwOfN3MdpvZzlmePxP4YcP94fSxZzCznWa2y8x2Pf7445TLZUZHRxkZGWF8fJyBgQEqlQp9fX3U6/UTp+FNXwDU3d1NvV6nr6+PSqXCwMAA4+PjjIyMMDo6SrlcZnBwkMnJSUqlErVajZ6enqetY/q2t7eXarVKf38/ExMTDA0NMTY2xtjYGENDQ0xMTNDf30+1WqW3t3fWdfT09FCr1SiVSkxOTjI4ONjSNplZ4baplc9pxYoVhdumVj+ngwcPFm6bWv2cjhw5UrhtavVzakWU03HNbKu7HzCz5wLfAG519283PP9l4A/c/YH0/jeB97r7SS/33LFjh2sMmkS5XGbz5s2xy4hOOQTKIlAWwbI5HdfdD6S3Y8AXgatmLDIMnNVwfxsw79Es7aoKJicnY5eQC8ohUBaBsshmyRsOMzvNzJ41/TPwauDhGYvdC/xSenbV1cARdx+db906OB7or6mEcgiURaAssonR49gCPGBmPcD3gC+7+1fN7BYzuyVd5n7gB8B+4C+BX2tmxZrIKRgeHo5dQi4oh0BZBMoim0INOdLV1eUa9TJRq9XUA0M5NFIWgbIIls0xjnbRkCPB3r17Y5eQC8ohUBaBssimUD0ODXIoIrIwHd/j0EU9gXbZJZRDoCwCZZGNehwiIh1MPQ71OE7QX1QJ5RAoi0BZZKMeh4hIB+v4HkelUoldQm5Mj3PT6ZRDoCwCZZFNoRqOtWvXxi4hNy688MLYJeSCcgiURaAssilUw6Erx4OhoaHYJeSCcgiURaAssilUw6ErQYMtW7bELiEXlEOgLAJlkU2hGo7jx4/HLiE3Dh8+HLuEXFAOgbIIlEU2hWo4Vqwo1OZkouM9CeUQKItAWWSj37QiIrIghWo4NJFToAEfE8ohUBaBssimUA3HypUrY5eQGxs3boxdQi4oh0BZBMoim0I1HLVaLXYJuXHw4MHYJeSCcgiURaAssilUw7FmzZrYJeTG2WefHbuEXFAOgbIIlEU2hWo4tN8y2LdvX+wSckE5BMoiUBbZLPkgh2Z2FnAX8DygDtzh7n88Y5mXAV8CHk0f+oK7f2C+dWuQQxGRhVkugxzWgPe4+48DVwPvMLOLZlnu39x9R/pv3kYDNKx6Iw0bnVAOgbIIlEU2S95wuPuou3enPz8BPAKcuRjrXrdu3WKsphC6urpil5ALyiFQFoGyyCbqMQ4z2w5cDnx3lqdfYmY9ZvYVM7u4mfWpxxHoL6qEcgiURaAssonWcJjZeuDzwLvcfWLG093AOe5+GfCnwD+cZD07zWyXme06cuQI5XKZ0dFRRkZGGB8fZ2BggEqlQl9fH/V6ne7ubiB8cbq7u6nX6/T19VGpVBgYGGB8fJyRkRFGR0cpl8sMDg4yOTlJqVSiVqvR09PztHVM3/b29lKtVunv72diYoKhoSHGxsYYGxtjaGiIiYkJ+vv7qVarJ+YDmLmOnp4earUapVKJyclJBgcHW9qmyy+/vHDb1Mrn1NXVVbhtavVz2rBhQ+G2qdXPaevWrYXbplY/p1ZEmQHQzFYD9wFfc/ePNrH8IHClu5dPttzFF1/se/fuXZwil7menh4uu+yy2GVEpxwCZREoi2BZHBw3MwM+CTwyV6NhZs9Ll8PMriKp80fzrVsDlwUXX9zU3r3CUw6BsgiURTYxdlVdA7wZeIWZ7Un/XWdmt5jZLekybwAeNrMe4E+AG7yJrlG1Wm1f1cvM/v37Y5eQC8ohUBaBsshmyWc+cvcHAJtnmY8DH1/ounXleLBt27bYJeSCcgiURaAssinUleMaqyool096OKhjKIdAWQTKIptCNRyayClYv3597BJyQTkEyiJQFtkU6jdtjDPE8urYsWOxS8gF5RAoi0BZZFOohkMCTWqVUA6BsgiURTaFaji0qyrQ8CsJ5RAoi0BZZFOo37Q6OB4cOnQodgm5oBwCZREoi2wK1XCsXr06dgm5sXXr1tgl5IJyCJRFoCyyKVTDMTU1FbuE3Hj00UfnX6gDKIdAWQTKIpsoY1W1iyZyCur1uo75oBwaKYtAWQTLYqyqdtKw6sGePXtil5ALyiFQFoGyyEY9DhGRDqYeh3ocJ2iimoRyCJRFoCyyUY9DRKSDqcehHscJ0zN+dTrlECiLQFlkox5HQemskYRyCJRFoCyCju9xHD16NHYJuVEqlWKXkAvKIVAWgbLIplANhyZyCs4999zYJeSCcgiURaAssilUw6GhkoMDBw7ELiEXlEOgLAJlkU2hGo5Vq5Z8Jtzc2rRpU+wSckE5BMoiUBbZRGk4zOxaM/u+me03s9tmed7M7E/S5x8ysyuaWa/G2A90hllCOQTKIlAW2Sx5w2FmK4E/A14LXATcaGYXzVjstcAF6b+dwJ8vaZEFoDNGEsohUBaBssgmRnpXAfvd/QfuPgV8Frh+xjLXA3d54t+BjWZ2xnwrNrPFr3aZ0hDzCeUQKItAWWQT46DAmcAPG+4PAy9uYpkzgdGZKzOznSS9EoCqmT28eKUua5uBcuwickA5BMoiUBbBCxb6ghgNx2zdgplXITazTPKg+x3AHQBmtmuhF7IUlbJIKIdAWQTKIjCzBV81HWNX1TBwVsP9bcDMc+OaWUZERCKI0XA8CFxgZuea2RrgBuDeGcvcC/xSenbV1cARd3/GbioREVl6S76ryt1rZvZO4GvASuBT7r7XzG5Jn/8EcD9wHbAfeAr45SZXf0cbSl6ulEVCOQTKIlAWwYKzKNQghyIi0n46mVlERBZEDYeIiCxIIRqO+YYw6RRmdpaZ/YuZPWJme83s12PXFJuZrTSz/zCz+2LXEpOZbTSzvzezUvr9eEnsmmIxs99I/388bGb3mNna2DUtFTP7lJmNNV7vZmabzOwbZtaf3p4+33qWfcPR5BAmnaIGvMfdfxy4GnhHB2cx7deBR2IXkQN/DHzV3V8IXEaHZmJmZwL/A7jS3S8hOUHnhrhVLanPANfOeOw24JvufgHwzfT+SS37hoPmhjDpCO4+6u7d6c9PkPxyODNuVfGY2Tbgp4G/il1LTGa2AXgp8EkAd59y98NRi4prFXCqma0C1tFB14i5+7eBQzMevh64M/35TuD1862nCA3HXMOTdDQz2w5cDnw3cikxfQx4L9DpwyY/H3gc+HS62+6vzOy02EXF4O4jwB8BQyRDGB1x96/HrSq6LdPXyaW3z53vBUVoOJoenqRTmNl64PPAu9x9InY9MZjZ64Axd98du5YcWAVcAfy5u18OPEkTuyOKKN1/fz1wLrAVOM3M3hS3quWnCA2HhidpYGarSRqNu939C7Hriega4GfNbJBk9+UrzOxv4pYUzTAw7O7Tvc+/J2lIOtFPAY+6++Pufgz4AvBfItcU28Hp0cfT27H5XlCEhqOZIUw6giXjyn8SeMTdPxq7npjc/X3uvs3dt5N8J/7Z3TvyL0t3/0/gh2Y2PQrqK4G+iCXFNARcbWbr0v8vr6RDTxRocC9wc/rzzcCX5nvBsp9rda4hTCKXFcs1wJuBXjPbkz72P939/nglSU7cCtyd/nH1A5ofxqdQ3P27Zvb3QDfJWYj/QQcNP2Jm9wAvAzab2TDwu8CHgc+Z2VtJGtY3zrseDTkiIiILUYRdVSIisoTUcIiIyIKo4RARkQVRwyEiIguihkNERBZEDYfIHMzsuJntMbMeM+s2s7ZdKGZmHzOzl8547P0z7l9qZp9pVw0izVLDITK3irvvcPfLgPcBf7AYK01HdG68vwm4Oh2ADjP7CTN7ELjFzL5nZq8AcPdeYJuZnb0YdYi0Sg2HSHM2AOPTd8zst8zsQTN7yMx+r+HxN6W/7PeY2V9MNxJmNmlmHzCz7wIz58J4A/DVhvsfJWmoPgG8Ctjf8Nw/0lnDgEsOqeEQmdupaQNQIhma/YMAZvZq4AKSIf13AF1m9lIz+3HgF4Fr3H0HcBy4KV3XacDD7v5id39gxvtcAzQOxjgFPA/A3Y+4+1DDc7uAn1y8TRRZuGU/5IhIG1XSBoB0xry7zOwS4NXpv/9Il1tP0pC8COgCHkyGQeJUwoBxx0kGn5zNGSTDnk/bCfwhcE06vtT73H0wfW6MZFRXkWjUcIg0wd2/Y2abgeeQDOX/B+7+F43LmNmtwJ3u/r5ZVnHU3Y/PsfoKcGL6UnfvA37GzH4fGCEZuPKV6dNr0+VFotGuKpEmmNkLSQbR/BHJgJpvSec9wczONLPnkky7+Yb05+m5nM9pYvWPAOc3vNcl6Y/HSAbje1bDshcCDyMSkXocInM7tWGUYQNuTnsNX0+PZ3wn3SU1CbzJ3fvM7HfS51eQ/OJ/B/DYPO/zZeBXCFPc/pqZ7SCZZ+b1wLsaln15urxINBodVyQHzOwB4HWNc4Gb2fvd/f0N908B/hX4CXevLXmRIintqhLJh/cAM6/P+NaM+2cDt6nRkNjU4xARkQVRj0NERBZEDYeIiCyIGg4REVkQNRwiIrIgajhERGRB/j+RcJ4MzwwL+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import laguide as lag\n",
    "\n",
    "x = np.array([4.5, 4.75, 5.5, 5.9, 6.25, 7])\n",
    "y = np.array([12, 13, 14.8, 15.5, 15.9, 17.25])\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "ax.plot(x,y,'ro');\n",
    "ax.set_xlim(0,10);\n",
    "ax.set_ylim(0,20);\n",
    "ax.set_xlabel('Beer (\\$)');\n",
    "ax.set_ylabel('Nachos (\\$)');\n",
    "\n",
    "ax.grid('True',ls=':')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we propose a linear equation $n  = x_1b + x_2$ that relates the nacho price, $n$, to the beer price, $b$.  If there were only two data points, we would solve a system of two equations to determine the parameters $x_1$ and $x_2$.  The larger number of data points leads to an inconsistent overdetermined system.\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "4.5x_1 + x_2 & = & 12 \\\\\n",
    "4.75x_1 + x_2 & = & 13 \\\\\n",
    "5.5x_1 + x_2 & = & 14.8 \\\\\n",
    "5.9x_1 + x_2 & = & 15.5 \\\\\n",
    "6.25x_1 + x_2 & = & 15.9 \\\\\n",
    "7x_1 + x_2 & = & 17.25 \\\\\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "We find the least squares solution by solving the normal equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.04147727]\n",
      " [3.20732008]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[4.5, 1],[4.75, 1],[5.5, 1],[5.9, 1],[6.25, 1],[7, 1]])\n",
    "B = np.array([[12],[13],[14.8],[15.5],[15.9],[17.25]])\n",
    "\n",
    "A_LS = A.transpose()@A\n",
    "B_LS = A.transpose()@B\n",
    "X_hat = lag.SolveSystem(A_LS,B_LS)\n",
    "print(X_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that this least squares solution produces a line that fits the data reasonably well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAys0lEQVR4nO3dfXxcZZn/8c+VtmlaSi21UEhLH8AiltY+hEWU1UVdXUAU2dUVt4uoaGUVVFRcWPytPIiwrj9Ff65iVxRwK+jiAxUQwQcEFIE2NqRNh6ah6ZA0JISmDUMnk07m+v1xJr2HkDSTOTNzzpy53q9XX5mHM2eu+c5p7tzn4b5FVTHGGGPyVRN0AcYYYyqLNRzGGGMmxBoOY4wxE2INhzHGmAmxhsMYY8yEWMNhjDFmQkrWcIjIsSLyexHZJiJbReRT2cdni8gDItKa/XnEGK8/Q0SeEpEdInJ5qeo0xhgzMVKq6zhE5BjgGFVtFJHDgU3Au4EPAntU9YZsg3CEqv7riNdOArYDbwM6gCeA96tqS0mKNcYYk7eS9ThUtUtVG7O3XwC2AfOAc4Bbs4vditeYjHQKsENVn1bVQeCO7OuMMcYEbHI53kREFgGrgMeAuaraBV7jIiJHjfKSecAzOfc7gNeNse61wFqA6dOnNyxZsoThXlRNTQ3pdJopU6YwODhIXV0d+/fvZ/r06S/7OTAwQG1tLQcOHGDy5MlkMpnh9ZPJZJg8eTKDg4NMnTqVgYEBpk2b9rJ1JJNJ6urqGBwcZPLkyQwNDVFT47XNmUyGSZMmkU6nqa2tHXcdqVSK2tpa0uk0NTU19pnsM9lnKvFnSqUGmTbN+0yJxHREBjjyyLqK/kz5fE9bt27tVdUjR/v9OpaSNxwiMgP4KfBpVe0XkbxeNspjo+5TU9V1wDqA5cuX6+bNmwusNFra2to4/vjjgy4jcJaDY1k4I7Nobob3vhfWrYPVqwMsLAAismuirynpWVUiMgWv0Vivqj/LPtydPf4xfBykZ5SXdgDH5tyfD+we7/0mTy5LB6oizJ49O+gSQsFycCwLZ2QW8+bBkUfC4GBABVWYUp5VJcDNwDZV/VrOUxuAC7K3LwDuGuXlTwBLRGSxiNQC52Vfd0jDXUcD+/fvD7qEULAcHMvC2b9/P7//PXzkI6AKs2fDww/DqacGXVllKGWP4zTgfOAtIrI5++8s4AbgbSLSinfW1A0AIlIvIvcCqGoauBj4Nd5B9Z+o6tYS1ho5w/tNq53l4FgWTk1NDbEYPPQQdHcHXU3lKdm+HVV9hNGPVQC8dZTldwNn5dy/F7h3Iu+Z5/GTqjBlypSgSwgFy8GxLOA3vwERWLFiCh/7GHzwgzBtWtBVVZ5I/Qliu6qcRCIRdAmhYDk41Z7F0BBceilcf72XRU2NNRqFitTRZDs47syZMyfoEkLBcnCqNYsHH4Q3vAFqa2HDBjj6aBgaqs4siiVSPY5BOyXioI6OjqBLCAXLwanGLJqa4M1vhm9/27u/eLHXy6jGLIqpZEOOBKGhoUE3bdoUdBmhkE6nrQeG5ZCrmrLYvRvq673bP/4xnHMO1NW556spi/GIyCZVPXkir4lUj2NgYCDoEkJj61Y7CQ0sh1zVksXXvw4nngjxuHf/fe97aaMB1ZNFqUSqyZ1mR7oOWrFiRdAlhILl4EQ9i3QaJk+Gc8+FvXth7tyxl416FqUWqR6HXeDk2C47j+XgRDULVfjwh2HtWu/+okVw9dUwderYr4lqFuUSqR7H9OnTgy4hNBoaGoIuIRQsByeqWYjAggVeA6Lq3R9PVLMoF+txRJT9ReWxHJwoZbFnD3zgAzD8ka66yutl5HsNcJSyCEKkGg7rcTj2F5XHcnCilIWIN1zIX/5S2OujlEUQItVwJJPJoEsIjebm5qBLCAXLwan0LJ57Dr78ZW931BFHQCzmDVJYiErPImiRajjqRp5zV8VOOOGEoEsIBcvBqfQsfvELb5fU8JQ7fv67V3oWQYtUw2FXjjvx4ZPYq5zl4FRiFj098MQT3u0LL4StW2HVKv/rrcQswiRSZ1XZlaDO3EOdxF5FLAenErN43/u8C/meesq7RmPJkuKstxKzCJNI9TiGhoaCLiE09u7dG3QJoWA5OJWSRXc3pFLe7RtvhF/+0ms0iqlSsgirSDUcNlGNY8d7PJaDUwlZdHfD0qXeQXCAFSu8+8VWCVmEme3bMcYE7sABmDLFGybkssvg3e8OuiJzKJH6E90mcnJswEeP5eCENYsHHoBXvQp27fLuX365N0hhKYU1i0pRsh6HiHwfOBvoUdVl2cd+DLw6u8gsYK+qrhzlte3AC8AQkM53yN9Jkyb5rjsqZs2aFXQJoWA5OGHNYskSeM1rvBn6yiWsWVSKUvY4bgHOyH1AVd+nqiuzjcVPgZ8d4vVvzi6b9zjx6XS6kDojqbu7O+gSQsFycMKUxQ9/CJ/6lHd70SK47z447rjyvX+YsqhEJWs4VPUhYM9oz4mIAP8I3F7M96ytrS3m6iraggULgi4hFCwHJ0xZPPWUN1xIUHuMwpRFJQrqGMcbgW5VbR3jeQXuF5FNIrI235Xafktn+/btQZcQCpaDE2QWqnDzR/9M4zHvgJoavvjDV/HgR9f7uvrbD9su/Amq4Xg/h+5tnKaqq4EzgU+IyJvGWlBE1orIRhHZuHfvXnp7e+nq6qKzs5O+vj7a2tpIJpO0tLSQyWRobGwE3OiYjY2NZDIZWlpaSCaTtLW10dfXR2dnJ11dXfT29tLe3k4ikSAWi5FOp2lqanrJOoZ/Njc3k0qlaG1tpb+/n3g8Tk9PDz09PcTjcfr7+2ltbSWVSh0cK2fkOpqamkin08RiMRKJBO3t7QV9ppNOOilyn6mQ72n58uWR+0yFfk91dXWBfaZdX/0u//69Bax79p2gypR4G6z9CEM//GEg296cOXNC+z2Ve9srREnnHBeRRcDdwwfHs49NBjqBBlUdd8Z4EbkKSKjqV8dbdunSpVpoEFGzadMmGwEUyyFXubNQhQ0b4F3vAlm8iPZdsIA4NeT8zlm4ENrby1bTMNsunEqZc/xvgdhYjYaIHCYihw/fBt4ObMlnxTasumP/KTyWg1PuLO66y7se4667gHicRex6aaMBbmLwMrPtwp+SNRwicjvwKPBqEekQkQuzT53HiN1UIlIvIvdm784FHhGRJuBx4B5VvS+f97SJnBybqMZjOTjlyCKTcR2Ic86Bn//c+8lYB6MDOkht24U/Jd1VVW4nn3yybty4MegyjKla//Iv3u6pbdtg5sycJ9av9yYFz/3jbvp0WLcO1qwpe53GqZRdVSVjEzk5wwfnqp3l4JQqi0zGGzIEvImVrr4aDj98xEJr1niNxMKF3vR9CxcG2mjYduFPpHocDQ0Nal1QTzqdtmHmsRxylSKLgQE44wx44xvh2muLuuqSsu3CqfoeR2p4LGbDjh07gi4hFCwHpxRZ1NXBa18Lxx9f9FWXlG0X/kSq4bArx5358+cHXUIoWA5OsbLYuRPOOssNSvjNb8IHP1iUVZeNbRf+RKrhsLGqnN7e3qBLCAXLwSlWFiKwZYs3bEilsu3Cn0g1HDaRkzNjxoygSwgFy8Hxk8VTT8F//Id3e9EiaGuDt7+9OHUFwbYLfyL1mzZKB/r9OjB8mkuVsxwcP1n86Efwla/As89696dMKVJRAbHtwp9INRzGsUmtPJaDM9Estm2DJ5/0bv/bv8HWrXD00SUoLAC2XfgTqfPRbFeVY8OveCwHZyJZDA3B2Wd7F3b//vcwdWp0Gg2w7cKvSP2mtYPjzp49o06FUnUsB+eQWaxfD4sW0SavQhcuYtId6/nRj+COO8pXXznZduFPpBqOKZW+47WI6uvrgy4hFCwHZ8wsssOBbN41i6Vs5eb438Latbxux3rmzi1vjeVi24U/kWo4BgcHgy4hNHbu3Bl0CaFgOThjZTHwb9fA/v2soIl/5xrexQZvTKkrryxzheVj24U/kRpyxAY5dDKZjB3zwXLINVoWN98MX/rIThpZzRHsfekLRLyBqCLItgun6occsWHVnc2bNwddQihUbQ7ZYxbU1Hg/168fNYuVK+EN05tQ5OXriPC83FW7XRRJpBoOO1PCWb16ddAlhEJV5jA8hPmuXd40fLt2wdq1rN62DVW46ipvBFuAhgZYv+5FZk8fMc7b9Olw3XVlL71cqnK7KKJINRzW43BslGBPVeZw5ZUvnfcCYP9+Up/7HCLeREvt7V6bAoRuyPNyqMrtoojsGIcxUVNTk9MqwCBTuJ4ruIDbWKQ7SafBRhQ3w+wYh/U4DmpsbAy6hFCoyhxGHJt4lqP5Kp/jJ4dfAFijAVW6XRRRpBoOO8bhrFy5MugSQqEqc7juOlLTZnEn/wDAAp5hW91qPvftVwVcWHhU5XZRRCVrOETk+yLSIyJbch67SkQ6RWRz9t9ZY7z2DBF5SkR2iMjl+b7nwMBAMUqPhFgsFnQJoVCVOaxZw3f+4QHey500sQIWLmT+964iZgeED6rK7aKISnaMQ0TeBCSA21R1Wfaxq4CEqn71EK+bBGwH3gZ0AE8A71fVlvHec/Xq1WpdUE8ymWTatGlBlxG4asphYMAbvXbRIkil4I9/hLe8xT1fTVmMx7JwQnWMQ1UfAgoZEOYUYIeqPq2qg8AdwDn5vNCGSnZ2794ddAmhEIkcRrkmYzTvfKf3b2jIG5Qwt9GAiGRRJJaFP0Ec47hYRJ7M7so6YpTn5wHP5NzvyD42KhFZKyIbRWTj888/T29vL11dXXR2dtLX10dbWxvJZJKWlhYymczBg2LDp+M1NjaSyWRoaWkhmUzS1tZGX18fnZ2ddHV10dvbS3t7O4lEglgsRjqdpqmp6SXrGP7Z3NxMKpWitbWV/v5+4vE4PT099PT0EI/H6e/vp7W1lVQqRXNz86jraGpqIp1OE4vFSCQStLe3F/SZjjjiiMh9pkK+p9mzZ1f0Z3ruxhvRj370ZddkPP2lLwHwpz81ksl4n+nSSwf5zGc6efHF0T/Tiy++GIrPFIZtT1Uj95kK/Z4KUdLTcUVkEXB3zq6quUAvoMC1wDGq+uERr3kv8Heq+pHs/fOBU1T1kvHeb8WKFTr8hVW7zs5O5s0bs72tGhWfw6JFbnLvXAsX8uyf2/mbv4FLL4WLLhp/VRWfRRFZFk6odlWNRlW7VXVIVTPAf+PtlhqpAzg25/58wPqVE2Tj8HgqPod4fMzH586Fv/5rOOGE/FZV8VkUkWXhT1nTE5Fjcu6eC2wZZbEngCUislhEaoHzgA15rt9/kRFhQ8x7Kj6HEddkPM5f8Tc8SN+8ZYh4gxSOPJYxlorPoogsC39KeTru7cCjwKtFpENELgS+IiLNIvIk8Gbg0uyy9SJyL4CqpoGLgV8D24CfqOrWfN7TpoN0EolE0CWEQsXncN113rhRWZMYokvq6fjE9RNeVcVnUUSWhT+RGnLETsd1EokEM2bMCLqMwEUhhz984QGavv1HPrn3GliwgKFrv8yk8/9pwuuJQhbFYlk4oT/GUWo2kZPT0dERdAmhEIUcbul8G98+6ipSyQy0txfUaEA0sigWy8KfSPU4Ghoa1Ea99KTTaSbboEQVm8PvfgfHHeedVNXf740v5XdEnUrNohQsC6fqexw25IizdWteh4UirxJz2LsXzj0Xrr3Wuz9zpv9GAyozi1KxLPyJVI/DhlU3lWzLFli2zLv96KOwYkVxGgxjDqXqexw2rLpju+w8lZLDhg2wfDncf793//WvL36jUSlZlINl4Y/1OIwJ0IsvwmGHweAgfOMbcPHFYGPvmXKyHof1OA6yv6g8Yc7hyivh1FO9kWxra+Gyy0rbaIQ5i3KzLPyJ1GkFNpGT09DQEHQJoRDmHE47zRvJtlzCnEW5WRb+RKrHkUwmgy4hNIZH1qx2YcphYAA+9CG46Sbv/llnwQ03eEOgl0OYsgiaZeFPpBqOurq6oEsIjRPyHfku4sKUw9Sp0N0Nzz8fzPuHKYugWRb+RKrhsCvHnfhYo6pWmaBz2LMHLrnE+ykCd9/tHdsIQtBZhIll4U+kGg67EtSZO3du0CWEQtA57NrljWD7hz9494MczTvoLMLEsvAnUg3HUDmPNIbc3r17gy4hFEqewyjTuj7/PPzv/3pPr1rlNR7nnlvaMvJh24RjWfgTqYbDJmdx7HiPp6Q5rF8Pa9e+bFrXa/4pxgUXQE+Pt9iRR5auhImwbcKxLPyx37TGFOrKKyF77VAvr6STeti/n6u3/SOPPgpHHRVwfcaUSKQOCthETo4N+OgpaQ7ZA6xpJnEqf+YEtnMv72BWxxZmrSjd2xbKtgnHsvAnUg3HpEmTgi4hNGbNmhV0CaFQyhxemP8aDn+mhckM8VU+x/G0eU+MmO41LGybcCwLfyK1qyqdTgddQmh0d3cHXUIolCqHJ5+E4/Zu4pe1/wDAu7mL5WzxRia87rqSvKdftk04loU/pZxz/Psi0iMiW3Ie+08RiYnIkyLycxGZNcZr27Nzk28WkbxHLaytrS1C5dGwIKR/9ZZbsXMYHhP0xBPhHX9fx/HXXgALF3oXaSxcCOvWwZo1RX3PYrFtwrEs/Cllj+MW4IwRjz0ALFPV1wLbgSsO8fo3q+rKiYzaaPstne3btwddQigUM4df/AJOP90NSnjLLbD08++E9nbIeNO6hrXRANsmclkW/pSs4VDVh4A9Ix67X1WH9yf9GZhfzPecZuNRH7R8+fKgSwiFYuYwdSocOBDckCF+2TbhWBb+BHmM48PAr8Z4ToH7RWSTiKw91EpEZK2IbBSRjR0dHfT29tLV1UVnZyd9fX20tbWRTCZpaWkhk8nQ2NgIuGGVGxsbyWQytLS0kEwmaWtro6+vj87OTrq6uujt7aW9vZ1EIkEsFiOdTtPU1PSSdQz/bG5uJpVK0draSn9/P/F4nJ6eHnp6eojH4/T399Pa2koqlTo4yNrIdTQ1NZFOp4nFYiQSCdrb2wv6TBs3bozcZyrke9q0aVPBn2lgIMWXv/wMP/iB99ozz4RvfnMT9fXBfqZCv6cHH3wwtN9Tube9Rx55JHKfqdDvqRAlnchJRBYBd6vqshGPXwmcDPy9jlKAiNSr6m4ROQpv99Yl2R7MIdlETqaYVOGMM7zDF7/6lffTmKgp2UROIlIjIqtE5B0i8hYRKXigFxG5ADgbWDNaowGgqruzP3uAnwOn5LNum8jJsYlqPBPNQRX+53+ygxL+aD13tLyWe349GVm8yLtSvILZNuFYFv4csschIscD/wr8LdAKPAfUAScA+4HvAreq6qhX3o3scYjIGcDXgL9R1efGeM1hQI2qvpC9/QBwjareN96HsR6H8au1FZYuhavO3cyV95x28MpwwDvVNsRnTRlTiFL0OL4E/A9wvKr+nar+s6q+J3tW1LuAVwDnj1HM7cCjwKtFpENELgS+BRwOPJA91fam7LL1InJv9qVzgUdEpAl4HLgnn0YDbCKnXMP7WKtdPjmowvDfG0uWwMMPwxWPnfvSRgO8+0GNiV4Etk04loU/JT3GUW4NDQ1qXVBPOp22YebJL4evfx0+9zloaoJlw0fjamrcRRu5RLxTbyuQbROOZeGU7BjHiDc5TUTOEAnfocJUKhV0CaGxY8eOoEsIhbFyUIV9+7zbw9O5Ll2as8BYF4hV8IVjtk04loU/4zYcInKbiJyUvX0R3u6mS4CbS1zbhNmV4878+UW9RKZijZXDe97jzZGRycCsWfDRj46YZOm667xjGrlCPJxIPmybcCwLfw7ZcIjIQrzTZl/I3v4YXqPxCeBUEVkgIjNLX2Z+bKwqp7e3N+gSQiE3h9w9T+ecA+973yFOsV2zxjsQXiHDieTDtgnHsvBnvJ18p+MdAD8DmArMAo4DjgcmZZ/fDDxZovomxCZycmbMmBF0CaEwnMNzz3kNxaWXwjvfCR/4QB4vXrOmohuKkWybcCwLfw7ZcKjqrSLyeuC9eI3GTap6W/Y02QtV9bYy1Ji3KB3o9+vAgQNBlxAKwzm84hXecCGJRMAFBci2Ccey8Cef0wo+DvwdMKiqv80+9krgspJVZXyzSa2grQ2+8IWZ3HKLN87UQw9V99Xftk04loU/4+7bUdWMqv4qp9FAVeOq+lhpS5s421XlTB95YLcKbd8O9947jSezO1KrudEA2yZyWRb+jHdw/Jci8k4RmTLKc8eJyDUi8uHSlTcxdnDc2bNnz/gLRVBrqzf8OcCZZ8KDD+7ir/4q0JJCo1q3idFYFv6Mt6vqo8BngBtFZA9uyJFFQBvwLVW9q6QVTsCUKS9r36pWfX190CUE4rLLYNMmr9GYOhVOPPHooEsKjWrdJkZjWfhzyB6Hqj6rqp9X1ePxDpBfi9eQLFPVt4Wp0QAYHBwMuoTQ2LlzZ9AllM1TT3mDEgL813/BY495jQZUVw7jsSwcy8KfSA05YoMcOplMpiqO+fT1eZdY/NM/eVd/j1QtOeTDsnAsC6csQ46EmQ2r7mzevDnoEkpquIdxxBHwve/BF784+nJRz2EiLAvHsvDHehym4tx3nzdkyO9/jx34Nsancg1yeISIvHairysH63E4URwlePhvnFNP9S7oPvbY8V8TxRwKZVk4loU/efU4RORBvPk3JuMNMfIc8AdV/Uwpi5so63FE1403wq9/DffcM2IwQmOML6XscbxCVfuBvwd+oKoNeLMChor1OJzhyemjYsYMOPxwmOhcXVHLwQ/LwrEs/Mm3x9EMvB24FbhSVZ8QkSezMwGGhvU4nEo/a+TAAbjhBli50huUULWwK78rPYdisiwcy8IpZY/jGuDXQFu20TgObw7yUBkYGAi6hNCIxWJBl+CLKtx5J/w2O9BNocOFVHoOxWRZOJaFPyU7q0pEvg+cDfSo6rLsY7OBH+Nded4O/KOq9o3y2jOAb+AN3f49Vb0hn/dcvXq1WhfUk0wmmTZtWtBlTMjgIHznO/Cxj0FdHfT3w0yfs71UYg6lYlk4loVTsh6HiMwXkZ+LSI+IdIvIT0VkvCm0bsGbxyPX5cBvVXUJ8Nvs/ZHvNQn4L+BMYCnwfhFZOnK50dhQyc7u3buDLmHC/vQn+PSn4a7seAR+Gw2ozBxKxbJwLAt/8t1V9QNgA1APzAN+mX1sTKr6EDByJLFz8I6TkP357lFeegqwQ1WfVtVB4I7s68Zlk887s2fPDrqEvAwOekOEAJx+OjQ2ehMuFUul5FAOloVjWfiTb8NxpKr+QFXT2X+3AEcW8H5zVbULIPvzqFGWmQc8k3O/I/vYqERkrYhsFJGN3d3d9Pb20tXVRWdnJ319fbS1tZFMJmlpaSGTyRw8m2L4PO7GxkYymQwtLS0kk0na2tro6+ujs7OTrq4uent7aW9vJ5FIEIvFSKfTNDU1vWQdwz+bm5tJpVK0trbS399PPB6np6eHnp4e4vE4/f39tLa2kkqlaG5uHnUdTU1NpNNpYrEYiUSC9vb2gj7Tiy++WBGfae3aBG99q/L44ztJJpNMnVrc72n//v2h/p7Kue21tbVF7jMV+j3t3r07cp+p0O+pEPmeVfUbvF1Pt2cfej/wIVV96zivWwTcnXOMY6+qzsp5vk9VjxjxmvcCf6eqH8nePx84RVUvGa/OFStW6PAXVu26uro45phjgi5jVKmU92/mTIjH4ckn4eyzS/NeYc6h3CwLx7JwSnlW1YeBfwSeBbqA92Qfm6huETkGIPuzZ5RlOoDca4LnA3ntkJRqn6knR1iHmD9wwLvy++Mf9+4vWFC6RgPCm0MQLAvHsvAnr4YjO+Pfu1T1SFU9SlXfraq7Cni/DcAF2dsXAKMNy/4EsEREFotILXBe9nXjsukgnUTIJtce/mqmTIEPfcgbMqQcwpZDkCwLx7LwJ6+jySJyJN6kTotyX6OqY/Y6ROR24HRgjoh0AF8EbgB+IiIXAnG8OT4QkXq8027PUtW0iFyMd93IJOD7qro1rw9jB8cPmjNnTtAlHPTUU96ghDffDKecAp/8ZPneO0w5BM2ycCwLf/L9TXsX8DDwG2Aonxeo6vvHeOplx0VUdTdwVs79e4F786ztIJvIyeno6ODEE08MugwAjj4aXvEKCOL6zDDlEDTLwrEs/Mn34PhmVV1Z+nL8aWhoUBv10pNOpwPtgT36qNfDWLfOG5Sw0CFD/Ao6hzCxLBzLwinlwfG7ReSs8RcLlg054mzdmtfevRK+PzzwAHR0ePeDOm8h6BzCxLJwLAt/DtnjEJEXAAUEOAxIAQey91VVi3Btb/HYIIfBevhhb3fU297m9TD274fDDgu6KmPMoRS9x6Gqh6vqzOzPGlWdlnM/VI0G2LDqucq9yy6TgU99Cq66yu2WCkOjYbsuHcvCsSz8yfcYx7nA71R1X/b+LOB0Vf1FSaubIOtxlN8f/wgNDd6ghDt3wpFHenNnGGMqQymPcXxxuNEAUNW9eKfXhor1OJxy/EUVi8Eb3+jNzgeweHH4Gg37y9KxLBzLwp98exwvm7RJRJpVdXnJKiuA9TjK49lnvVNsAe64w5toKQy7pYwxE1fKHsdGEfmaiBwvIseJyNeB0DXZyYnOKxphwwOkFdtNN8GSJfD00979884Ld6NRqhwqkWXhWBb+5Hsi8yXA/8GbhEmA+4FPlKqoQtXV1QVdQmiccMIJRV3f0BBMmgTveAfs2uV6HGFX7BwqmWXhWBb+5DtW1YuqermqnqyqDap6haq+WOriJsquHHfi8XhR1qMKF18MH/ygd//YY+H662H69KKsvuSKlUMUWBaOZeHPRMaq+jxwEnDwz3pVfUuJ6iqIXQnqzJ07tyjrEfF6F9One6fc1uS7czMkipVDFFgWjmXhT76/BtYDMWAxcDXefOFPlKimgg0N5TWMVlXYu3dvwa/t74e1a93MfF/4AnzlK5XXaIC/HKLGsnAsC3/y/VXwSlW9GTigqn/Ijop7agnrKkhNJf5mKxE/x3tUveFCHn+8iAUFxI57OZaFY1n4k+++nQPZn10i8g68iZXml6YkE4S9e+G734XLLvNGst26tXKOYxhjyivfP9G/JCKvAD4LfA74HnBpyaoqkE3k5Ex0wMcNG+DKK10vIyqNhg186VgWjmXhT149DlW9O3tzH/Dm0pXjz6RJk4IuITRmzZo17jJ79sCOHd7kSuef7/2M2hQF+eRQLSwLx7Lw55ANh4j8+yGeVlW9tsj1+JJOp4MuITS6u7uZOfPQ41CuWQNbtkBbG9TWRq/RgPxyqBaWhWNZ+DNej2O0azUOAy4EXgmEquGora0NuoTQWLBgwaiPP/+8d6V3XR3853/CgQNeoxFVY+VQjSwLx7LwZ7xh1f/v8D9gHTAN+BBwB3BcIW8oIq8Wkc05//pF5NMjljldRPblLHOons9Btt/S2b59+8see/55WLYMrr7au79sGaxaVebCymy0HKqVZeFYFv6Me4xDRGYDnwHWALcCq1W1r9A3VNWngJXZdU8COoGfj7Low6p69kTWPW3atELLipzly934k+k0TJ4Mr3ylN2fGmWcGWFiZ5eZQ7SwLx7Lw55A9DhH5T7wL/V4AlqvqVX4ajVG8FWhT1V3FWJkNq+4MDxv90EMvHZTw8sthxYoACyszGz7bsSwcy8Kf8U7H/SxQD3wB2J3drdQvIi+ISH8R3v884PYxnnu9iDSJyK9E5KSxViAia0Vko4hs3LdvH729vXR1ddHZ2UlfXx9tbW0kk0laWlrIZDI0NjYCbsNpbGwkk8nQ0tJCMpmkra2Nvr4+Ojs76erqore3l/b2dhKJBLFYjHQ6TVNT00vWMfyzubmZVCpFa2sr/f39xONxenp66OnpIR6P09/fT2trK6lU6uDonCPX0dTURDqdJhaLkUgkaG9vL+gzrVq1ipaWFo45ZoBjjtlPb+++iv9MhXxPDQ0NkftMhX5PM2fOjNxnKvR7qq+vj9xnKvR7KkRe83GUgojU4l1IeJKqdo94biaQUdWEiJwFfENVl4y3zqVLl2qhQUTJ//4v3HlnDz/+8VFBlxK4TZs20dDQEHQZoWBZOJaFU8r5OErhTKBxZKMBoKr9qprI3r4XmCIic8Zb4fSoXLXm07Zt0N5+FC+Gbvzi8rNfDo5l4VgW/gTZcLyfMXZTicjRIiLZ26fg1fn8eCu0iZw8V1wB3/52U6gnWCqX4V0GxrLIZVn4E8iuKhGZDjwDHDc8l7mIXASgqjeJyMXAvwBpIAl8RlX/NN56Gxoa1A56edLptA0zj+WQy7JwLAunYnZVqep+VX3lcKORfewmVb0pe/tbqnqSqq5Q1VPzaTQAUqlUqUquODt27Ai6hFCwHBzLwrEs/InUOOR25bgzf74NXgyWQy7LwrEs/IlUw2FjVTm9vb1BlxAKloNjWTiWhT+RajhsIidnxowZQZcQCpaDY1k4loU/kfpNG9Q1KWF04MCB8ReqApaDY1k4loU/kWo4jGOTWnksB8eycCwLfyLVcNiuKscuhvRYDo5l4VgW/kTqN60dHHf27NkTdAmhYDk4loVjWfgTqYZjypQpQZcQGvX19UGXEAqWg2NZOJaFP5FqOAYHB4MuITR27twZdAmhYDk4loVjWfgT2Oi4pXDyySfrxo0bgy4jFDKZjB3zwXLIZVk4loVTMUOOlIpN5ORs3rw56BJCwXJwLAvHsvDHehzGGFPFrMdhPY6DbJRgj+XgWBaOZeGP9TiMMaaKWY/DehwHDc8xXO0sB8eycCwLf6zHEVF21ojHcnAsC8eycKq+xzEwMBB0CaERi8WCLiEULAfHsnAsC38i1XDYRE7O4sWLgy4hFCwHx7JwLAt/Amk4RKRdRJpFZLOIvGzfkni+KSI7RORJEVmdz3ptqGRn9+7dQZcQCpaDY1k4loU/Qc7W/mZVHWsarjOBJdl/rwO+k/15SDb5vDN79uygSwgFy8GxLBzLwp+w7qo6B7hNPX8GZonIMeO9yMbYd+wMM4/l4FgWjmXhT1ANhwL3i8gmEVk7yvPzgGdy7ndkH3sZEVkrIhtFZONzzz1Hb28vXV1ddHZ20tfXR1tbG8lkkpaWFjKZzMHT8IYvAGpsbCSTydDS0kIymaStrY2+vj46Ozvp6uqit7eX9vZ2EokEsViMdDpNU1PTS9Yx/LO5uZlUKkVrayv9/f3E43F6enro6ekhHo/T399Pa2srqVSK5ubmUdfR1NREOp0mFouRSCRob28v6DOJSOQ+UyHfU01NTeQ+U6HfU3d3d+Q+U6Hf0759+yL3mQr9ngoRyOm4IlKvqrtF5CjgAeASVX0o5/l7gOtV9ZHs/d8Cn1fVQ17uuXLlSrUxaDy9vb3MmTMn6DICZzk4loVjWTgVczququ7O/uwBfg6cMmKRDuDYnPvzgXGPZtmuKieRSARdQihYDo5l4VgW/pS94RCRw0Tk8OHbwNuBLSMW2wB8IHt21anAPlXtGm/ddnDcsb+mPJaDY1k4loU/QfQ45gKPiEgT8Dhwj6reJyIXichF2WXuBZ4GdgD/DXw8nxXbRE5OR0dH0CWEguXgWBaOZeFPpIYcaWhoUBv10pNOp60HhuWQy7JwLAunYo5xlIoNOeJs3bo16BJCwXJwLAvHsvAnUj0OG+TQGGMmpup7HHZRj2O77DyWg2NZOJaFP9bjMMaYKmY9DutxHGR/UXksB8eycCwLf6zHYYwxVazqexzJZDLoEkJjeJybamc5OJaFY1n4E6mGo66uLugSQuOEE04IuoRQsBwcy8KxLPyJVMNhV4478Xg86BJCwXJwLAvHsvAnUg2HXQnqzJ07N+gSQsFycCwLx7LwJ1INx9DQUNAlhMbevXuDLiEULAfHsnAsC38i1XDU1ETq4/hix3s8loNjWTiWhT/2m9YYY8yERKrhsImcHBvw0WM5OJaFY1n4E6mGY9KkSUGXEBqzZs0KuoRQsBwcy8KxLPyJVMORTqeDLiE0uru7gy4hFCwHx7JwLAt/ItVw1NbWBl1CaCxYsCDoEkLBcnAsC8ey8CdSDYftt3S2b98edAmhYDk4loVjWfhT9kEOReRY4DbgaCADrFPVb4xY5nTgLmBn9qGfqeo1463bBjk0xpiJqZRBDtPAZ1X1NcCpwCdEZOkoyz2sqiuz/8ZtNMCGVc9lw0Z7LAfHsnAsC3/K3nCoapeqNmZvvwBsA+YVY93Tp08vxmoioaGhIegSQsFycCwLx7LwJ9BjHCKyCFgFPDbK068XkSYR+ZWInJTP+qzH4dhfVB7LwbEsHMvCn8AaDhGZAfwU+LSq9o94uhFYqKorgP8H/OIQ61krIhtFZOO+ffvo7e2lq6uLzs5O+vr6aGtrI5lM0tLSQiaTobGxEXAbTmNjI5lMhpaWFpLJJG1tbfT19dHZ2UlXVxe9vb20t7eTSCSIxWKk02mamppeso7hn83NzaRSKVpbW+nv7ycej9PT00NPTw/xeJz+/n5aW1tJpVIH5wMYuY6mpibS6TSxWIxEIkF7e3tBn2nVqlWR+0yFfE8NDQ2R+0yFfk8zZ86M3Gcq9Huqr6+P3Gcq9HsqRCAzAIrIFOBu4Neq+rU8lm8HTlbV3kMtd9JJJ+nWrVuLU2SFa2pqYsWKFUGXETjLwbEsHMvCqYiD4yIiwM3AtrEaDRE5OrscInIKXp3Pj7duG7jMOemkvPbuRZ7l4FgWjmXhTxC7qk4DzgfeIiKbs//OEpGLROSi7DLvAbaISBPwTeA8zaNrlEqlSld1hdmxY0fQJYSC5eBYFo5l4U/ZZz5S1UcAGWeZbwHfmui67cpxZ/78+UGXEAqWg2NZOJaFP5G6ctzGqnJ6ew95OKhqWA6OZeFYFv5EquGwiZycGTNmBF1CKFgOjmXhWBb+ROo3bRBniIXVgQMHgi4hFCwHx7JwLAt/ItVwGMcmtfJYDo5l4VgW/kSq4bBdVY4Nv+KxHBzLwrEs/InUb1o7OO7s2bMn6BJCwXJwLAvHsvAnUg3HlClTgi4hNOrr64MuIRQsB8eycCwLfyLVcAwODgZdQmjs3Llz/IWqgOXgWBaOZeFPIGNVlYpN5ORkMhk75oPlkMuycCwLpyLGqiolG1bd2bx5c9AlhILl4FgWjmXhj/U4jDGmilmPw3ocB9lENR7LwbEsHMvCH+txGGNMFbMeh/U4Dhqe8avaWQ6OZeFYFv5YjyOi7KwRj+XgWBaOZeFUfY9jYGAg6BJCIxaLBV1CKFgOjmXhWBb+RKrhsImcnMWLFwddQihYDo5l4VgW/kSq4bChkp3du3cHXUIoWA6OZeFYFv5EquGYPLnsM+GG1uzZs4MuIRQsB8eycCwLfwJpOETkDBF5SkR2iMjlozwvIvLN7PNPisjqfNZrY+w7doaZx3JwLAvHsvCn7A2HiEwC/gs4E1gKvF9Elo5Y7ExgSfbfWuA7ZS0yAuyMEY/l4FgWjmXhTxDpnQLsUNWnVXUQuAM4Z8Qy5wC3qefPwCwROWa8FYtI8autUDbEvMdycCwLx7LwJ4iDAvOAZ3LudwCvy2OZeUDXyJWJyFq8XglASkS2FK/UijYH6A26iBCwHBzLwrEsnFdP9AVBNByjdQtGXoWYzzLeg6rrgHUAIrJxoheyRJVl4bEcHMvCsSwcEZnwVdNB7KrqAI7NuT8fGHluXD7LGGOMCUAQDccTwBIRWSwitcB5wIYRy2wAPpA9u+pUYJ+qvmw3lTHGmPIr+64qVU2LyMXAr4FJwPdVdauIXJR9/ibgXuAsYAewH/hQnqtfV4KSK5Vl4bEcHMvCsSycCWcRqUEOjTHGlJ6dzGyMMWZCrOEwxhgzIZFoOMYbwqRaiMixIvJ7EdkmIltF5FNB1xQ0EZkkIn8RkbuDriVIIjJLRO4UkVh2+3h90DUFRUQuzf7/2CIit4tIXdA1lYuIfF9EenKvdxOR2SLygIi0Zn8eMd56Kr7hyHMIk2qRBj6rqq8BTgU+UcVZDPsUsC3oIkLgG8B9qnoisIIqzURE5gGfBE5W1WV4J+icF2xVZXULcMaIxy4HfquqS4DfZu8fUsU3HOQ3hElVUNUuVW3M3n4B75fDvGCrCo6IzAfeAXwv6FqCJCIzgTcBNwOo6qCq7g20qGBNBqaJyGRgOlV0jZiqPgTsGfHwOcCt2du3Au8ebz1RaDjGGp6kqonIImAV8FjApQTpRuDzQLUPm3wc8Bzwg+xuu++JyGFBFxUEVe0EvgrE8YYw2qeq9wdbVeDmDl8nl/151HgviELDkffwJNVCRGYAPwU+rar9QdcTBBE5G+hR1U1B1xICk4HVwHdUdRXwInnsjoii7P77c4DFQD1wmIj8c7BVVZ4oNBw2PEkOEZmC12isV9WfBV1PgE4D3iUi7Xi7L98iIv8TbEmB6QA6VHW493knXkNSjf4W2Kmqz6nqAeBnwBsCrilo3cOjj2d/9oz3gig0HPkMYVIVxBtX/mZgm6p+Leh6gqSqV6jqfFVdhLdN/E5Vq/IvS1V9FnhGRIZHQX0r0BJgSUGKA6eKyPTs/5e3UqUnCuTYAFyQvX0BcNd4L6j4uVbHGsIk4LKCchpwPtAsIpuzj/2bqt4bXEkmJC4B1mf/uHqa/IfxiRRVfUxE7gQa8c5C/AtVNPyIiNwOnA7MEZEO4IvADcBPRORCvIb1veOux4YcMcYYMxFR2FVljDGmjKzhMMYYMyHWcBhjjJkQaziMMcZMiDUcxhhjJsQaDmPGICJDIrJZRJpEpFFESnahmIjcKCJvGvHYVSPuLxeRW0pVgzH5sobDmLElVXWlqq4ArgCuL8ZKsyM6596fDZyaHYAOEflrEXkCuEhEHheRtwCoajMwX0QWFKMOYwplDYcx+ZkJ9A3fEZHLROQJEXlSRK7Oefyfs7/sN4vId4cbCRFJiMg1IvIYMHIujPcA9+Xc/xpeQ3UT8DZgR85zv6S6hgE3IWQNhzFjm5ZtAGJ4Q7NfCyAibweW4A3pvxJoEJE3ichrgPcBp6nqSmAIWJNd12HAFlV9nao+MuJ9TgNyB2McBI4GUNV9qhrPeW4j8MbifURjJq7ihxwxpoSS2QaA7Ix5t4nIMuDt2X9/yS43A68heS3QADzhDYPENNyAcUN4g0+O5hi8Yc+HrQX+AzgtO77UFarann2uB29UV2MCYw2HMXlQ1UdFZA5wJN5Q/ter6ndzlxGRS4BbVfWKUVYxoKpDY6w+CRycvlRVW4B3ish1QCfewJVvzT5dl13emMDYripj8iAiJ+INovk83oCaH87Oe4KIzBORo/Cm3XxP9vbwXM4L81j9NuBVOe+1LHvzAN5gfIfnLHsCsAVjAmQ9DmPGNi1nlGEBLsj2Gu7PHs94NLtLKgH8s6q2iMgXss/X4P3i/wSwa5z3uQf4GG6K24+LyEq8eWbeDXw6Z9k3Z5c3JjA2Oq4xISAijwBn584FLiJXqepVOfenAn8A/lpV02Uv0pgs21VlTDh8Fhh5fcaDI+4vAC63RsMEzXocxhhjJsR6HMYYYybEGg5jjDETYg2HMcaYCbGGwxhjzIRYw2GMMWZC/j8MkjoaIPnNKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_model = np.linspace(2,10,40)\n",
    "y_model = x_model*X_hat[0,0] + X_hat[1,0]\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "ax.plot(x,y,'ro');\n",
    "ax.plot(x_model,y_model,'b',ls=':');\n",
    "\n",
    "ax.set_xlim(0,10);\n",
    "ax.set_ylim(0,20);\n",
    "ax.set_xlabel('Beer (\\$)');\n",
    "ax.set_ylabel('Nachos (\\$)');\n",
    "\n",
    "ax.grid('True',ls=':')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises \n",
    "\n",
    "**Exercise 1:** The sales of a company (in million dollars) has been shown in the table below:\n",
    "\n",
    "\n",
    "|Year|Sales|\n",
    "|----|------|\n",
    "|2004|13 |\n",
    "|2006|29|\n",
    "|2008|35|\n",
    "|2010|45|\n",
    "|2012|48|\n",
    "\n",
    "\n",
    "$(a)$ Find the least squares regression line. Plot the regression line and the data points to show that the line fits the data reasonably well.\n",
    "\n",
    "$(b)$ Use the regression line model to estimate the sales of the company in 2020.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** The concept of least squares can also be used to find polynomials of best fit that are of degree greater than one. Given the points $(-1,0),(0,1),(1,3),$ and $(2,9)$, find the parabola $y = ax^2 + bx + c$ that best fits through the points in the sense of least squares. Then plot the points on top of a graph of the parabola to observe the fit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** Lets revisit a problem we saw earlier to see how the technique of least squares compares to our previous method of interpolation. Here are the average radii of the orbits of the planets in our solar system, and their average orbital velocity around the Sun.\n",
    "\n",
    "|Planet                           | Distance from Sun (million km)  | Orbital Velocity (km/s)         |\n",
    "| ------------------------------- | ------------------------------- | ------------------------------- |\n",
    "|Mercury                          | 57.9                            | 47.4                            |\n",
    "|Venus                            | 108.2                           | 35.0                            | \n",
    "|Earth                            | 149.6                           | 29.8                            |\n",
    "|Mars                             | 228.0                           | 24.1                            |\n",
    "|Jupiter                          | 778.5                           | 13.1                            |\n",
    "|Saturn                           | 1432.0                          | 9.7                             |\n",
    "|Uranus                           | 2867.0                          | 6.8                             |\n",
    "|Neptune                          | 4515.0                          | 5.4                             |\n",
    "\n",
    "$(a)$ Use the method of least squares to find the cubic polynomial of best fit for these points. Plot the data points together with the graph of the polynomial to observe the fit.\n",
    "\n",
    "$(b)$ The dwarf planet Ceres has an average distance from the sun of 413.5 million km. What does the polynomial suggest the value would be? How closely does this agree with what the points seem to suggest?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Information Retrieval\n",
    "\n",
    "In the earlier section on the applications of [vector spaces](Applications_VS.ipynb) we discussed how documents could be represented as vectors, and how matrix algebra might be used to identify which documents most closely match a set of keywords defined by a search.  Here we describe an alternative that makes use of the fact that dot products can be used to compute the angles between vectors, which can be interpreted as a measure of how closely the vectors align.\n",
    "\n",
    "Let's again consider the specific example of a database which contains a list of webpages with content related to information retrieval, and that the set of searchable keywords is small.\n",
    "\n",
    "**\\{algorithm, engine, information, google, computations, matrix, optimization, retrieval, search, theory \\}**  \n",
    "\n",
    "Again, suppose we look at just a few titles in such a database.\n",
    "\n",
    "- Search Engine Algorithms\n",
    "- How Google Search Optimization Works\n",
    "- Information Retrieval Theory\n",
    "- Matrix Models of Information Retrieval\n",
    "- Theory Behind Search Engines\n",
    "- Computations in Information Retrieval\n",
    "\n",
    "Previously, each webpage title was represented by a $1\\times 10$ row vector that contained entries of ones for keywords that were included in the title, and zeros for those that were not.  (*For example, Search Engine Algorithms was represented as $[1, 1, 0, 0, 0, 0, 0, 0, 1, 0]$.*)  We now make one minor adjustment to this representation.  Each of these row vectors will be scaled by its magnitude to make it a unit vector.  Thus *Search Engine Algorithms* will be represented by the following vector.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5774 0.5774 0.     0.     0.     0.     0.     0.     0.5774 0.    ]]\n"
     ]
    }
   ],
   "source": [
    "T1_raw = np.array([[1, 1, 0, 0, 0, 0, 0, 0, 1, 0]])\n",
    "T1_scaled = T1_raw/lag.Magnitude(T1_raw.transpose())\n",
    "print(np.round(T1_scaled,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will represent the database as a matrix with one such row for each webpage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5774 0.5774 0.     0.     0.     0.     0.     0.     0.5774 0.    ]\n",
      " [0.     0.     0.     0.5774 0.     0.     0.5774 0.     0.5774 0.    ]\n",
      " [0.     0.     0.5774 0.     0.     0.     0.     0.5774 0.     0.5774]\n",
      " [0.     0.     0.5774 0.     0.     0.5774 0.     0.5774 0.     0.    ]\n",
      " [0.     0.5774 0.     0.     0.     0.     0.     0.     0.5774 0.5774]\n",
      " [0.     0.     0.5774 0.     0.5774 0.     0.     0.5774 0.     0.    ]]\n"
     ]
    }
   ],
   "source": [
    "T1 = np.array([[1, 1, 0, 0, 0, 0, 0, 0, 1, 0]])\n",
    "T2 = np.array([[0, 0, 0, 1, 0, 0, 1, 0, 1, 0]])\n",
    "T3 = np.array([[0, 0, 1, 0, 0, 0, 0, 1, 0, 1]])\n",
    "T4 = np.array([[0, 0, 1, 0, 0, 1, 0, 1, 0, 0]])\n",
    "T5 = np.array([[0, 1, 0, 0, 0, 0, 0, 0, 1, 1]])\n",
    "T6 = np.array([[0, 0, 1, 0, 1, 0, 0, 1, 0, 0]])\n",
    "\n",
    "D = np.vstack((T1,T2,T3,T4,T5,T6))\n",
    "D_scaled = lag.ScaleMatrixRows(D)\n",
    "print(np.round(D_scaled,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform a search of the database, we construct a $10 \\times 1$ vector built in a similar way from the keywords in our search.  Suppose that we want to search for entries matching the words \"information\", \"retrieval\", and \"matrix\".  The unit vector representing our keywords is $[0, 0, 0.5774, 0, 0, 0.5774, 0, 0.5774, 0, 0]^T$.  \n",
    "\n",
    "At the [start of this chapter](Inner_Products.ipynb) we gave the following formula relating the dot product of two vectors with the cosine of the angle between them.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\cos{\\theta} = \\frac{U\\cdot V}{||U||||V||}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "When $U$ and $V$ are unit vectors the formula simplifies to $\\cos{\\theta} = {U\\cdot V}$.  This is the scenario we have constructed by scaling the search vector and database matrix rows to have unit magnitude.  The dot product of our search vector with any row of the database matrix tells us directly the cosine of the angle between the vectors.  In this application, the value of the cosine will always be between zero and one since the entries are all positive values.  The closer the value is to one, the *smaller* the angle between the vectors and the better the directions of the two vectors agree.  \n",
    "\n",
    "The result of the multiplying the search vector by the database matrix is a vector that contains the dot products of the search vector with each of the rows.  The location of the largest entry corresponds to the row in the matrix that best agrees with the search vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.    ]\n",
      " [0.    ]\n",
      " [0.6667]\n",
      " [1.    ]\n",
      " [0.    ]\n",
      " [0.6667]]\n"
     ]
    }
   ],
   "source": [
    "S_raw = np.array([[0,0,1,0,0,1,0,1,0,0]])\n",
    "S_scaled = S_raw.reshape((10,1))/lag.Magnitude(S_raw.transpose())\n",
    "\n",
    "result = D_scaled@S_scaled\n",
    "print(np.round(result,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fourth entry in the vector is one, which means that the fourth webpage in the database contains all of the keywords in the search vector and contains none of the *searchable keywords* that are not in the search vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise 1:** Write a function that performs the same task as the $\\texttt{laguide}$ function $\\texttt{ScaleMatrixRows}$.  That is, it should accept an $m \\times n$ matrix as an argument, and return a matrix where each entry has been divided by the magnitude of the row in which it is located.  Test your function by comparing results with $\\texttt{ScaleMatrixRows}$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** The search method described in this section often yields the same results as the one described in [Chapter 3](Applications_VS.ipynb).  Construct an example where the two methods yield *different* results.  \n",
    "(*Hint: If we ignore the actual words involved, we see the question is one of matrix multiplication.  Let $A$ be a $4\\times 6$ database matrix consisting of ones and zeros and let $A_u$ be the matrix returned by $\\texttt{ScaleMatrixRows(A)}$.  Let $X$ be a $6\\times 1$ search vector consisting of ones and zeros and let $X_u = X/||X||$.  To build the required example, arrange the ones and zeros in $A$ and $X$ such that the largest entry in $AX$ appears in a different location than the largest entry in $A_uX_u$.*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3:** Suppose that we want to search for entries matching the words \"information\", \"retrieval\", and \"computations\". Use the information discussed in this section to determine the webpage in the database which contains all or most of the keywords in the search vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4:** The following example of a database with webpages containing the content related to Global Warming. The set of searchable keywords is: **\\{causes, effects, global, warming, fossil-fuels, climate , potentials, change, actions, prevent \\}**\n",
    "\n",
    "Some titles in this database are as follows:\n",
    "\n",
    "- Causes of Global Warming\n",
    "- Global Warming vs Climate change\n",
    "- Fossil Fuels - A Major cause of Global Warming\n",
    "- Actions to prevent Global Warming\n",
    "- Potentials of Global Warming\n",
    "- Effects of Global Warming\n",
    "\n",
    "Suppose the keywords we want to search are \"potentials\", \"global\" and \"warming\". Create the search vector and database matrix. Analyse the resulting vector to find the webpage containing the most relevant information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Code solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "- Berry, Michael W. and Murray Browne. *Understanding Serach Engines: Mathematical Modeling and Text Retrieval*. 2nd ed., SIAM, 2005 \n",
    "\n",
    "- Leon, Steven J. *Linear Algebra with Applications*. 9th ed., Pearson., 2015"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}