01_LinearRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reproduce results\n",
    "set.seed(1234567)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data-Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Random data in which y is a noisy function of x\n",
    "X <- runif(100, -5, 5)\n",
    "y <- X + rnorm(100) + 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Linear Regression in R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A linear regression assumes that the true function (we are trying to replicate) can be written as:\n",
    "\n",
    "\\begin{equation*} \n",
    "   y = X\\beta + \\epsilon\n",
    "\\end{equation*}\n",
    "\n",
    "In other words the true function is a linear combination of its parameters.\n",
    "\n",
    "The vector of residuals given an estimate of $\\beta$ is thus:\n",
    "\n",
    "\\begin{equation*} \n",
    "   e = y - X\\beta^{OLS}\n",
    "\\end{equation*}\n",
    "\n",
    "Where OLS indicates that this is the Ordinary Least Squares estimate of $\\beta$\n",
    "\n",
    "The OLS estimator by definition minimises the sum of squared residuals:\n",
    "\n",
    "\\begin{equation*} \n",
    "   e'e = (y - X\\beta^{OLS})'(y - X\\beta^{OLS})\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*} \n",
    "   e'e = y'y - y'X\\beta^{OLS} - \\beta^{OLS'}X'y + \\beta^{OLS'}X'X\\beta^{OLS'}\n",
    "\\end{equation*}\n",
    "\n",
    "Since the tranpose of a scalar is a scalar: $y'X\\beta^{OLS} = (y'X\\beta^{OLS})' = \\beta^{OLS'}X'y$ we get the Residual Sum of Squares (RSS):\n",
    "\n",
    "\\begin{equation*} \n",
    "   RSS = e'e = y'y - 2\\beta^{OLS'}X'y + \\beta^{OLS'}X'X\\beta^{OLS'}\n",
    "\\end{equation*}\n",
    "  \n",
    "We take the derivative of this w.r.t to beta-hat and set it equal to 0:\n",
    "\n",
    "\\begin{equation*} \n",
    "   0 = -2X'y + 2X'X\\beta^{OLS}\n",
    "\\end{equation*}\n",
    "\n",
    "The chart below shows an example of such a function; for the example we assume that wages are a linear function of height:\n",
    "\n",
    "```R\n",
    "X <- runif(100, -5, 5)  # Height\n",
    "y <- X + rnorm(100) + 3  # Wages\n",
    "```\n",
    "\n",
    "The line-of-best-fit here is a line that minimises the squared sum of residuals, it has a slope of 0.95 and an intercept of 2.95."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Call:\n",
      "lm(formula = y ~ X)\n",
      "\n",
      "Coefficients:\n",
      "(Intercept)            X  \n",
      "     2.9455       0.9519  \n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD///8kdWheAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAgAElEQVR4nO3diZqiPLdA4Yhzaavc/822TArIkISdCdb7nL+/6uoq\nofqwmiSgqhzAYir0DgBrQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQ\nAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQ\nAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEjzVNv812p9rrZX\n6lF+8Hp/0av86KHU3n5fZ/frNLFL38/q/xg/nz1Z72PSCGmew5AuSl3KD27vL7p3PyWq3od/\nmZrYJYGQmg1szkZ/bCMOQ/qn1KH84PT+ourf8oNS/+z3dWa/ml1xFtL839BKbfTHNmJycJge\nSFk9ovtE+h7jZUaPYGRJSNOPaPyNK7PRH9tI9+B4KnUsP6hOHcUfXt9TndOz87W34/uj/fnZ\n+lzxn/v7m06P1qO9z0S3vDwz7asz0b0+M7UfIH++v+xw/Tz465yprP6j1+X9kOp4+z5k/UDN\nrOtZzrnKb/2eVItf/x0+e/3zgw5tqflcd2c6j6R31l6lTf7QhnqHxrWaxJyrI744isujJ2tH\nc2hGgv/ydkjn6pOtkupu3n/w/uicl2Xd+w/wr/rwUD/QM/v+UfNxPUIsXKoHej9e2eit3N3f\nkM6tvf75QQe29JlmdXem80iEhAn9Q+NQHOu3egj2nT3tv1/7bu3wKvM45u2Qaq2FrVf1h9n7\nu/flI1Rjve4DZN9vzVu/LbZfndFe7126Ng9ZnYI+s673Hz2HQvrdl9+Q2lvqf27opyIkTFCq\ne+y8j9Ts2vp3Ont/cM+qM0n1Jfvy4G1+9z3ksnt1gLcevBwgPouTyLn4pkd1buk8wK36xlvW\nrvRVJaSqL3y118z3ZYvfZvadfWh+pGKvz90Jmer9pJ0t1d/b25n+I200I0LS0Q+pOMA+/5ar\n6ihrxmjdA6kfUpHaq/sl5Wr3pfijf8VZ5dI6tTTfeqxXxm/Vdx4/6xPH6vxwunf391LsUjUC\n+1fs1yUfCmlgX/ohdbaUfz7X3pn+IxESRv2EVPyb3/xT/vlc51/+92nrdj58zgkTS8/lQKw5\nhxzqgVj3AbLvJrrDsqyMRvVbehZNv08Rt895bnQfpkPqbKn+2t7O9B+JkDDq9+AoDt9L7w87\n0dz23eOx/ZW9h3sfmfV47lB+VAXafoDfTbQe/Nwc6q1Vg/37QYr/FTOvrDPkmg2p89FPVr8/\nLyHVNvpjG/k5OF7FgCp7NX/4an2VaiYSan+6PnRCOpcr30WW13IBsFhx6z5A7ySQdb//datW\n+A6dh7zXZ6N6LdAupGzgTzkjDdvoj23k5+Co1rs/c6RyUnNvTyT29URCJ6R7+WDFwsXj81H3\nAX7nSL1JUX7vrmBU86P7Z540vg/TIXW31P3cbfgRCQmj+gfH+/jMHtn3CC1WG4pVu2veS0fr\njPRZXqtXlttfUz1As1CmPqe7Yp3sVp6E9p/1gPbyW/NA5Ziv/YDFf16j+9IPqbOl3qqdGv6p\nmg1sDiHNU215NaspTiT77h+25iKHckB1z7RCKu9MKD86qdZdE60HGLmOVJT83o3D83MRtnGu\nT5jFA57bG83qT+iF1NnSz+cGf6pMdfdkMwhpXi+keoX68llXbt8jUB1P/5q2PncR5eMh3VS9\ngv65GaH/ANXo73MzQf3b6oBtFhsO7Yf8Vz/QXTV3wNbfelLts8tsSJ0tdT/XurNhaAObQ0jz\nuiF97rU7fm4ZuO1Vdm4vOeSP9wGVnR7V186E9P6iesn7+6Sk7gOUvzvcP9/5Or9HdMd6+lLO\njw7XzkO274X9/AzNTrcvd82E1NlS87nuzvQf6di9WWIzCGmhfhYOub0x3FBUOxMBQlrIQ0j1\nOO1xaM6FIUW1MxEhpIU8hHT4Diz7697+RbUzESGkhTyE9HmqRAzrYVHtTEQIaSEfc6TXpbgE\nnPVvTg0jqp2JByEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCA\nkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCA\nkAABhAQIICRAACEBAggJEOAhJAUkxuIolw8nwCYASYQECCAkQAAhAQIICRBASIAAQgIEEBIg\ngJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkYNzf35/eFxISMKasSC8lQgLG\n/LV+nUFIwIi/3n+nEBIwgpAAAYQESGCOBAhg1Q4QwXUkwCdCAgQQEiCAkAABhAQIICRAACEB\nAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEB\nAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIWAPtNwRzhZCQPoO3qHSFkJA+gzdN\ndoWQkLy/3n+1vkf4BEZISJ55SPJjQUJC8ixCMvz6eYSE9Jl2YTMWnEFISJ/pSE0jpME51MTE\nipCwBmZrB7MhDZY5mSshYYPmxoKDfz75TYSEDZoZCw6esaZPY4SETZocCxISIICQAAnMkQAB\nzRyqMwBk1Q4w9VdmlHfL4ToStkLyZlSTGyYICSsiejOq0Y1EhIQVEb0ZlZCwUbI3o0Yb0vOk\nskueX/cqOzvaBDZN+K7uSOdIr0y9XS/Fr+rgZBPYNumQDGZcO/OHtz3Kz+p9Hjpn6vTKX+XH\n8pvAxkk/YU93DXC383hGyspvVOpV/idzsQlsXKCXE9p5Hdop9f21+U/nj1ssN4HNC/ECd8Ww\nLsAZqfj1xRkJK7Erp0ceQ2rmSOdX/bH8JgDf6lUGVu2ABZrVOq4jAdZ2n1Vv7mwAbLUuHhES\nYKl9EZaQACu7zs0MhATY6N0TREiAhf69dYQEGNv93KNKSICpgVu9CQkw83s6ygkJMDT8zCNC\nAkyMPIOPkAB9g8O6AiEB2safUE5IWCM3z++beGEGQsL6uHnG+eiwrkBIWB/5Ny3P514niJCw\nOg7etHz29bYICavjIKTJYV2BkLA68iHNv/wjIWF9pOdIGi+jSkhYH9lVu9lhXYGQsEaC15H0\nXtWbkIApmq+OT0jAOK1hXYGQgFH6b9ZCSECtP7HSPh3lhATUfpb6jN47jJCAUv/ik9l78BES\nUOjdDmEyrCsQElDohmT8lrCEBBQ6IZm/tTIhAaXvHMl0WFcgJKD0WbWzyIiQgI/qOpJVR4QE\ntNkM6wqEBHxZZkRIQIt1R4QENGyHdQVCAioLMiIkoLLkdJQTElBalhEhAYWlHRESsHRYVyAk\nbN7yjAgJkOiIkLBxAsO6AiFh02QyIiRsm1RHhIQNExrWFQgJmyWXESFhuyQ7IiRslOCwrkBI\n2CTZjAgJ2yTdESFhg4SHdQVCwubIZ0RI2BwHp6OckJAm+/eIdZIRISFFC9613FFHhIQE9d/K\nSJubYV2BkJCc3lsZ6XOWESEhQbYhOeyIkJAeu5DcDesKhIT02MyRnGZESEiRxaqd444ICUky\nvI7kdlhXICSsn/OMCAkb4KEjQkLKdEZ47od1BUJCmv7+/rTWHLxkREhIU6uhyZD8nI5yQkKa\nyo7mryf5yoiQkKTmVDRzSvLXESEhRVoheRvWFQgJCWrGdH/5REfe9qZASEhRdSqaWrXz2xEh\nIUl1Q6PXkbwO6wqEhDRNXov1nREhYY129q+NYouQsDa7nf1ro1gjJKzMzv61URYgJKzLzv61\nUZYgJKxJsVpHSD43gTUqV+sIyecmsEL1qjdzJI+bwOp8LsIueEVja4SElWhfhOU6krdNYF28\n3xPUQ0iIlNFZJXBGhIRImc1zgndESIiTycpb6GFdgZAQI5NrQRIZLV6dICTEyCAkgY4E1st9\nhvQ6Z+9fL3ulDjdHm8BKaIckMqwTuILrMaRnplT+ev9SODjZBFZD89gWmR1J3FPkMaSTOr7e\nv5ye76ZO6uxiE1gNvdGWzCpDYiEp9ap/eY/yVOZiE1iR+fm/1GpdciG9f8lU6ze9P26x3AS2\nRG7RO6050kk98vxS/FKckSYnSYSEWYIXj9JatXuo7PzIj9m7pPte3V1sAmv3Ge8JX4RN6jrS\nPfuO3S5uNoFV+545IriXocvvBdnbaV9UdLw8nW0CK9bMZWK4J6iHOxuQjGb0FV9GhISENNMj\n30/a00FISEYZ0G7n+dUY9BAS0vFXDuti7IiQ4Neidea/v533VzXRREjwaOGVz12AV8fXREjw\naNm9ODGu1jUICf4suzs05o4ICR4tCSnCi7BthAR/FoQUd0aEBK+s50ixd0RI8Mly1S7yYV2B\nkOCVzfp1/BkREkLSqyqFjggJwWi+vkkSHRESgtFaeUgjI0JCMDpr4YmcjnJCQjAaISWTESEh\nmPmQEuqIkBDMzBwpnWFdgZAQyvSqXVIZERJCmriOlFhHhIQYpTWsKxAS4pNcRoSECCXYESEh\nNukN6wqEhI8oXllkWUbBfgRCQk3gvU0ELOoo4I9ASKgJvNvWYguHdQF/BEJCReL9H5daODsK\n+SMQEirhQ1q8ykBIATaBnuAhLV+sI6QAm0Df/ATD6ZKYxKI3cyT/m0Df3JKX0yUxmYtHdrso\n8u8DIeFj+ohy+c/9b0aWR/fPt80+jtC/D4QEPS4nID8dSZ39NB5H6N8HQoKe+oBzMLYbGNZJ\nnf005n29/9oiJOhp/m0XT2lgdiR1dGs8DiHBs7/W/wkaWmUgJCGEFKG/6nyUi5Y0vFrnMSTm\nSPCuGdXJhTS26O1vjsSqHbwTX7gbvXjkcdWO60jwTvaENHkRVmpNw9cTlAgJ+kRvbkjyibCj\nCAkm5P6BX1dHhIQg0nxhhgmEhADWlhEhIYDVnY5yQoJ/K8yIkODdKjsiJPi1xmFdgZDg00oz\nIiR41evIw20H3NmA1ekN6zy8Lqq/l14lJPjSH9Z5eM0ffy8rREjwZLgjp4e5xxe6IyR48bta\nR0jyexFgE+hxOyt3+MIMEwjJ/SbQ4XhWPvGEcuZIDhGSb06PuJGLsKzaOUdInjkdA41fhDUb\nTloNPrmOBI8chiR1T1Acbyg4ipCQuwxJ7J4gf9MdK4SEgqvDVLijeEsiJBTcDJxmh3X6MxhC\nskFI/jmYlc9nlGvXS0g2CGkNZod1RuNJ5kgWCCl986t1ZicZVu0sEFLyNFYZTEdrvi4JWSEk\nuKCzWhf7tMcIIUGe5kXYyKc9RggJ4nQvHkU+7TFCSJBmcBE26mmPEUKCrLW+3tYMQoKobWZE\nSJDlsKO4h4GEBDntYZ3wcR/7wgQhQUwno1z2uI99qZyQIKSzyiB93Ed/8ZaQIKMzO2of9yJj\nPEKyQkjJ6b0a8ee/QmM8QrJCSInpXzxqhdT5vT3mSDYIKS2/i97NcS92JmHVzgYhJWXgzVqa\n474TkvVsqfxGriNFuQlIGXmzluq4b4VkfU6J/WRUIiQsM/NmLd/fWs9yYp8elQgJi8y9Wcvn\ndGI9W4p+wa5ESFhA581a6rkNIQl8S4SbgACTN2shpMlvue7z/LlX+39SO/S7CcTK7M1amCNN\nfMtdvX+XqTfRkggpAaZv1sKq3cS3HNQtf6h9flMHk0eY2yohxc/izVqWXUeK28KQiiQe6qzR\nRucRCCl1G31C+TiBkI7qrhWS6hLeK/hERn2Lh3aPu8pyraHdv4yQ0jI6oqKjH8sXG5S6FGeb\n+/w3vo7q8CwfYWir2pXBj9E5PsO6AYuXv7NihpTvb1rfelOq+ELmSCkYW3UmoyGeL8g+D+r4\nIqQUjF0HpaNB3u9suKjsTkgJGA6JYd2IxSHdj+XK3VP72x/7+TkQIYU3GBIZjVka0qFaGlCZ\nfkn5iZBSMDBHoqNRC0O6qsOryOKqTmK7lBNSFH5W7RjWTVgYUqbqpQPZFWtCikL3OhIZTRG4\ns4GQNkH/TY+ivy/OhYUh7eszUnHjqiBCio3usC6JO7VdkJkj3TN1FdulnJCioz2sS+K5Qy4s\nXbU71vf0mDyLwnATCG6nO2BL49msLohcR1JHvTuELDeBsHY77QEbIbn9lgg3AV07gwEbIbn9\nlgg3sWU/47SJgdtuZ5QHcyS7b2k98+FwFtspQnLpZ5w2NXArVxlMQmLVzupbOk8iykLuFXT9\nnDQmziK7zh/prTdsMaPlQ7tTVjyj756pf/lRiZ2TCMmdnyrGM/lcPNrsgE3fwpDO6lH+96EO\n+UvuoiwhuaMf0vfi0WYHbPoEbhH6fCB3mxAhuaMdUve9LMlo2uKbVpszUkZIidCbI3Grt5nF\nQ7tmjnQ2fJFI6b2CrtlVu+LsQ0aGRJ7YV90ipOTutyMkpyavI5Uf7hjIGRK6Rag4LRUvyyWE\nkML5K4d1hGSIOxvQ8de7CAs9hISOv95FWOiRCunfcemezG4CPjSrdYRkZmlIZyevMkxIoey4\ni8HO4uXvhsZrf9ttAv6Ut3pzF4ONxRdkb/lBPZ8H3rEvFVNPmZj9CowRuEXo8j4bPWSfa05I\nrsw+ZQJ2BEK6FxdimSOlYXwGxD1BiywM6fge2j3VPv9HSEmYeMqE3x1ZneVvNFbfJsRLFstx\nN0nReMqEnx1ZnaXL35fidycl95y+302s0/gx6nLZbCSk4WEd63cGuLMhiKlj1OmFnOGnTBh8\nLYYtCCk7Xv+9RHfmZxOrNXGMun1Fq6GCJzuiJD0LQqpe8OR4uRu8NZLhJtZq6hi1P371ZjT9\nrxpdrSMkEwtCev27nvbVbQ3H8+0ReK+S4iIkyxnN+GodIZlYPEd6XE8H7rWb0zsNTB6jllMT\nu2+bWvVmjmRAZrHhfiCkCb8ni6lj1O7UYnX+mL4Iy6qdAYGQ/l32nJEm/WYzfYzqTHb6X2MT\n0uxFWK4jaVsY0vNavq/L/iJ6z+rKQho8xpcdo2PnOKOQuJlB0JLFhvs5K9btTnfxRfANhCTx\nkAPzLv1tcG+dqIXL38er/Np3Tkg2j2g4oyEjWYtCOju6HruykOSXv4bTNBgtcjqSxhnJA/Hl\nr7+6GduHJCNxi+ZIJ+ZImmSXv8pHe//P7DG/+0BH8li1S9BfXpVkEtL3rMiwzgWuI6WnGtUZ\nnuQ+8zQycoI7G9Jjswr4+R46coN77dLzicLgnFR/pdawjvsZLHD3d4LqsV1usBRYfZ3Oi+Nz\nh50Vno+UoNaxrn3Al8sMWu9L3voV2niGbJK+oy/tU9Jf8WYtGickw8dFhddsSJTxAb/TfAKt\n6eOiREiJMj3gdVfrCMkOIaXKaC5jcBGWOZIVQkqVyeqaycUjVu2sEFK6dK/3mN4TxHUkC4S0\netzL4AMhrR0deUFI68at3p4Q0qqRkS+EtGZ05A0hrRfDOo8IabXcZ8Q6+RchrZXzjrhy20ZI\n6+RhWMe9RG2EtEoeZkfc3dpBSGvkY5WBkDoIaX38rNYRUgchrY6vRW/mSG2EtDbeLh6xatdG\nSOvSG9a5vdLDdaQvQlqVfkY55wxPCGlF+qsMnVkMpw+nCCkSAsd5f3bUXlfj5OQYIUVB4jj/\nWWXohNT5jCFOZvMIKQrLl5IHLh61Qlpy0YeTmQ5CisHyi5uDi97fPBeFZP2dW0JIMVgc0vDF\no++5ZMEGuINBCyHFYOHBWg7rBicyn0/an1YISQshRWHR8KnKKJ+cyNhPdAhJCyFFYcmEvhzW\nzZdovfTGHEkHIUXC9jjftTpyc7izaqeDkNJWrzK4HX9xHWkeISWtWa1jIhMaISWsdRGWiUxg\nhJSu9sUjJjKBEVKqfm71JqOQCClRvIpqXAgpTXQUGUJKEa/qHR1CShAZxYeQ0kNHEfIZ0uuk\n1OFeP8jkoxDSOIZ1UfIY0itThWP1IFsIaXJF2nK5mozi5DGks7q+a7pmh/JB1h/S5DVS2wuo\ndBQpjyFl1Tc+s/1zGyG1fjX7w9aXdVtjWBctjyE17bwOh6GQVJvlJmIyeR9p+cm/udFd/7RF\nRvHyGNJevZqPDhs4I82FpDG4a5+23tHRUcQ8hnRVp/qjpzoQUn1S0nuE99ftdtxMFzGfy9/n\nTz33mdHbGkKamSP91X84ta73/e9f9cIMgnsHWV4vyD6OzUfP0wZC0lq10wrpr54eUVK0uLPB\noanFhGalYX5s9/ddrSOkaBGSZyavNNectnY8kzx6hORVa7indUn2r1xlaM5fdBQvQvKqcxrS\nu0no8yqqPAU2ZoTkk8UQrXnZutmrtwiKkHwyDsny5R+JzjtC8sm0CLuXfxycfdGWW4Tk1cxS\nXf8e1c536Yf0+8W8WpdrhOSVyVMrbF/+cSg7Xj/SNULybGKI1T3arV/+cSAkrkM5R0jR6B7t\n9i//SEghEFI02kf7omfw/Y7jCMk5QopG62hf9syjgYEgcyTXCCke33tUFz/S7+J3rj/JggVC\nikd9tLt5YQauI7lFSDEp71ENvROwQUhxsT0dccIJjJCiYp1RzhQoLEKKie2wjkW54AgplN/B\nmPUqA5eJwiOkMAYGY/arDIQUHiGF8TsYW7BaR0jhEVIQP4f+sotHzJGCI6Qg+iEtvHjEql1w\nhBREL6TlF2G5jhQYIYXRHozxZi0rQEhhtAZjZLQGhBRKMxijo1UgpLBCDeuYUwkjpKCCZZSz\nyieLkEIKNazjupM4QvJF8N66pbgTQh4h+SF6b91ShCSPkPz4GUyFvHhESPIIyYvfe+sC7UiF\nOZI4QvJC+N66pVi1E0dIXvReRTX8RViuIwkjJD8699aF3BG4QUh+RHZvHSckaYRkzfBgrL88\ngmEdUyQHCMmS5cEYQUYs2rlASJbsDsahjrwPs7iM5AAh2bE6GIeGdQGGWYTkACHZsTkYB4d1\nAYZZhOQAIdmxOBgnOgpSEh1JIiRLpgfjyGpdmJBYtRNHSJYMD8ax1TqJkCxWK7iOJI2Q5o0c\ndQYH48TFo8XDLE4vUSCkOQIH6tTFo8UPz4QnCoQ0Z/mBOnMRdtkwiyW4OBDSjMUHquN7gggp\nDoQ0Y+mB6vqeIEKKAyHNWHigur+3jjlSFAhpzpID1cet3qzaRYGQ5iw4UD3d6s1FoQgQ0jzb\nAzWKp0zAD0JyxXRYx3klaYTkiHFGOTOdlBGSG6bDOtbeEkdILhiv1nE1KHWE5ID5KgMhpY6Q\n5Fms1hFS6ghJmt1FWOZIiSMkYZYXj1i1SxwhiVpwTxDXkZJGSJK4l2GzCEkQHW0XIY0zHGzF\n8KreCIWQxphO/8lo0whpjOGCNB1tGyGNMLtEyrBu6whphFFIZLR5hDTCJKSRVyPmwtCGENIY\n7TnS8LCOWxW2hZDG6JYwMqzj5rltIaRxWmOzyY4oaTMIaZHR1TpC2hhCWmLuTSYIaTMIyd7k\nxSPmSNtCSNZGR3XfZQpW7TaDkGyNrTK0KiKj7SAkOzOrDCS0NYRkZXR2xCLDRhGSDVbr0ENI\n5qZW6whpowjJ2PSt3syRtomQTM2+tXLOsvcGEZIZjWfwsey9RYRkhGfwYRghmaAjjCCktulR\nGS/MgFGE9DWzTkBGGOc1pH+Xoyocz/9cbWKJ6ZVrOsIEjyG99urr4GQTi0xeS2VYh0keQzqr\n7PYoP3reM3V2sYlFpkIiI0zzGFKmHp+PHypzsYlFxkPidIQ5HkNSauw39WdaLDexzNgciYww\nizPS18iqHR1hnt850v1ZfhTnHCkfvI7EsA46fC5/H1pjt/3LySbEaWXE3XXwex3pXF5Hyo6X\nKK8jDdDpiPu9wZ0Nk/SGdTwDCdsNSWc0pjc74jmxyLcaktZoTHOVgZCQbzak1q8jtFfrCAn5\nRkPSOPYNFr2ZI4GQRphcPGLVDoQ0yPQiLNeRsMmQeOYRpG00pInRGPcEwcI2Q5oYjYXMiCFi\nurYa0piAHbFokTJCatMa1rk6b7CMnjJCauln1Eum/K2z8wYXdpNGSF+9jnrJdN6Fj5DQlXJI\nsmOsn2FdL5nub+UPd0JKWrohCY+xfmZHvQO7+U/3t5KYI6Us4ZBavy73u8rgPyRW7RKWbEii\nI6Gh1brhkJyeN7iOlC5CyscuHg3PkThvYAghjV6EHV6147yBIcmGJDbGmrgIO3QdCRiScEgy\nYyxuUYWEdEOSOUPQEUSkHNJyPGUCQjYdEhlByoZD4nQEOdsN6ZMRi3FYbrMhNR1xgRUSNhrS\nd1jHraKQsM2QvrMjnrwAEcmHZDPDaa0yEBJEJB6SzQyns1pHSBCRekitXzX1n1Bu/gjAj7RD\nsjif/L7ASc6qHRbbWEiDz+AjIyy2rZC4lwGOpB2S4QyHjuBK6iEZzHA031qZgR4sJB6SwYGv\nmVHO0gMsJB+SJs1bvVkMh52NhMQ7lMOtbYTEO5TDsS2EpP8MPkKCpQ2EZPQO5a1fAX3rD8no\n4hGrdrCz9pCMX5iB60iwsfKQuJcBfqw7JDqCJ2sOidfbgjcrDomM4M96Q6IjeLTWkBjWwauV\nhrRjFRterTOkHddV4dcaQ9rtuNMHnq0wpB33nsK71YVUrDIQEnxbW0jlYh0hwbeVhVQvejNH\ngmerCmnHmx4hkDWF1HlxfDKCTysKiXsZEM5qQuKeIIS0lpD6b9bCyA5erSSkbkesNcC3VYTU\nH9ax+g3f1hDSz3uH9f4LOLeCkH5WGQgJ3iUf0sBqHSHBu9RDGlz0Zo4E39IOaeTiEat28C3p\nkMavwXIdCX6lHJLme4fRFNxLNyTdt4TNGeXBvWRDMnoPvj9OTHAr1ZBM38uSExOcSjMk7Vu9\nu5eUCAmuJBmS/jMm6nSacxElwZEUQ7J/L0tCgiPphWT2DL7u5IiQ4EhyIRk/EbZcrmOOBLdS\nC8nyCeWs2sGttEJa8MIMXEeCS0mFxOubIFYphURHiFY6IfF6W4hYMiGREWLmM6TXSanDvX6Q\nyUf5+UNOR4ibx5BemSocqwcxComMEDmPIZ3V9V3TNTuUD2ISEh0hdh5DyqpvfGb7p1FIDOsQ\nP48hNYEikQUAAAWySURBVO28DgeTkMgICfAY0l69mo8O+iENdcRdCoiNx5Cu6lR/9FQHzZCG\nhnVm980RHXzwufx9/tRzV3ohLX75R25WhR9eL8g+js1Hz9PPo6i26lMTHWmWxNMn4EfEdzaM\nvYxq779TeEIfPIk3pLHVOkJChEKEND0/qjcx8XLErV9nEBI8iTSkqYuwJgsIzJHgR5whzVyE\n1V/SZtUOfsQZ0tgXWFwU4joSfEgpJE4viFZSIbV+BaIS7/L3D5bgEC9CAgQQEiAgoZCYIyFe\nSYXEqh1ilVJIXBRCtNIKCYgUIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiA\nAEICBBASIICQAAGEBAggJEAAIQECIg0JSIzFUS4fjg/hdzv4HrADEezBV0z7YiD8bgffA3Yg\ngj34imlfDITf7eB7wA5EsAdfMe2LgfC7HXwP2IEI9uArpn0xEH63g+8BOxDBHnzFtC8Gwu92\n8D1gByLYg6+Y9sVA+N0OvgfsQAR78BXTvhgIv9vB94AdiGAPvmLaFwPhdzv4HrADEezBV0z7\nYiD8bgffA3Yggj34imlfDITf7eB7wA5EsAdfMe2LgfC7HXwP2IEI9uArpn0BkkVIgABCAgQQ\nEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAhIO6V/Yfb/uVXZ+\nBdr4OQu48ULQn74R+AjoiGhXDL2yoPt+Lt+1IAtzLB3Kje+DbLsU9KdvBD4CuiLaFUNHmzff\nEPNQp/dRdFWnEBv/p7JH/sjUvxAbLwT96T/CHgE9Ee2KmZvVu9iIOVYbD7MPZ3XPi7+BS4iN\nF4L+9I3AR0BPRLti5KkOMfw1htmHo3rmxWnhGGLjLUH/PxDJEdCIaFeMHNQzgr/GlzqE2KyK\n4YQQ7KevxXEEfES0KyYu6hb8OMqLWcI9xGYjCSnQT1+J5Aj4iGhXDJSjmvB/jc8szOAqjpBC\n/fSlSI6Ar4h2xcC+WHgN/tf4ygINbaIIKdhPX4rjCGiJaFc01G85fSrHFEH+Gttven0IdSEn\niyGkYD99IeARMCKiXdFQH8ZL3sZdZg/envvD0/vmK9Wq3TPkql3An74Q8AgYEc2OmIjhr/Ee\ncMnqUv57fFfnYHsQ8qcvxHAEdEWzI+ZCX8UIt/HgdzYE/em/4smIkGydgv6LuC83He5gDvvT\nf4TefltEu2Iq6F9j2KHFq7z7O8imS5EMrEJvvy2iXQHSRUiAAEICBBASIICQAAGEBAggJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAA\nIQECCAkQQEiAAEICBBASIICQAAGEBAggJEAAIQECCAkQQEiR+7wrXe/t6cZ/e3e9SxhASJEz\nDmnP/0tD4G89cmMhjX1ZVG+suiH8rUeOkNLA33rkeiFd9yq7fn97ztS5/PD9v7PKLvUbjgfa\n1y3j7zxy3ZCOZSeH5reH4nenKqTyj66EFAh/55FTX3l+V4dX/jqoexXSXWWP/JFVIb3/5Kr2\nDO0C4W89cp2Qjur1/tRLHatejkVR75zKkP7leTPIg3/8rUeuM7RrNdVq5vshIQXD33rkCCkN\n/K1HrhdS5/OEFA/+1iPXCamaFPV+eyekCPC3HrlOSLdimS6/NosNnVW75muUegbc3c0ipMh1\nryOVF45U9uz8thvS/v3nwfZ2uwgpcr93NqjT8/Pbc6YO/7oh/dsTUgCElL7yTgeERUgJU+qW\n56+jOofeERBSyi7VDImRXAQIKWXXg1J7zkcxICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQI\nICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQI\nICRAACEBAggJEPAf+kYLYQQHHwgAAAAASUVORK5CYII=",
      "text/plain": [
       "Plot with title \"Explain Wages with Height\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fit a model (regress weight on height)\n",
    "fit <- lm(y ~ X)\n",
    "print(fit)\n",
    "# Coefficients:\n",
    "#(Intercept)            X  \n",
    "#     2.9455       0.9519 \n",
    "\n",
    "# beta-hat\n",
    "fit_params <- fit$coefficients\n",
    "\n",
    "# Plot\n",
    "plot(x=X, y=y, cex = 1, col = \"grey\",\n",
    "     main = \"Explain Wages with Height\", xlab = \"Height\", ylab = \"Wages\")\n",
    "\n",
    "# Draw the regression line (intercept, slope)\n",
    "abline(a=fit_params[[1]], b=fit_params[[2]], col=\"blue\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Linear Regression from Scratch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that $X'X$ is a positive definite matrix (our variables are not a perfect linear combination of each other & we have more observations than variables) we can find a closed-form solution for $\\beta^{OLS}$:\n",
    "\n",
    "\\begin{equation*} \n",
    "   \\beta^{OLS} = (X'X)^{-1}X'y\n",
    "\\end{equation*}\n",
    "\n",
    "In R we can run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           [,1]\n",
      "intcpt 2.945535\n",
      "       0.951942\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAOVBMVEUAAAAA/wBNTU1oaGh8\nfHyMjIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD//wD///+mzl89AAAACXBI\nWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3diWKqPLuGYUStrXVV4fwPdglJmBECLxngvvb+\nWjtB2p9nZSQkOYDVEt8FAPaAIAECCBIggCABAggSIIAgAQIIEiCAIAECCBIggCABAggSIIAg\nAQIIEiCAIAECCBIggCABAggSIIAgAQIIEiCAIAECCBIggCABAggSIIAgAQIIEiCAIAECCBIg\ngCABAggSIIAgAQIIEiCAIAECCBIggCABAggSIIAgAQIIEiCAIAECCBIggCABAggSIIAgTUua\npr931ue0c5L8lS9e7296la/+kuS8vKyT5fr6UKT6s/N/jd5nvxaXMWoEadqGQfpOku/yxe/7\nmx7tT4nSZfiXJh+KJBAkc4LDOeivbWXDIP1Lkkv54uv9Terf8kuS/Fte1olymaJsFqTpv9BO\nHfTXtmJzcdheSKlu0VUhfbfxUqsjWFkTpM9HtP7BnTnor22lfXE8k+RavlBVR/HFn3dX5+vZ\n+t7f6/vV+fZsfK5493j/0Ndf42jvmug3L2ums6qJHrpmah4gf76/7fJTHfx1S5NUf+n1/T5k\ncv2tD6kPZHpdz7LPVf5oXakWb/9dqlL3ftGhM5nPtQvTOtK8WnuXDvlLW+pcGj+qE3NTV3xx\nFZdXT9oMzcW0BP/lzSDd1CcbSdK5eX/h/eqWl8l6dA/wT7286AM90/pL5rVuIRa+1YHexysz\n+lsWtx+kW6PUvV904ExVN6tdmNaRCBI+6F4al+Ja/9VNsLr3dK6/9521y6uMxzVvBklrDGy9\n1BfT90+fyyOotl77AGn9o3njw+L8qkZ7vYv0Yw6pqqCq1/X+0nMoSP2y9IPUPFP3c0O/FUHC\nB0nSvnbeV2r60/h3On2/eKSqJlHfci4vXvNRfcmlD3WBNw5eNhCfRSVyK37oT9UtrQP8qh/8\nTZspfakIJeobX80x83OZxToz51YZzK9UlPrW7pAlnd+0dSb9s53CdI900BgRpDm6QSousOrf\n8kRdZaaN1r6QukEqovZqf0s52v1dfOlfUat8N6oW86NXPTL+q37yWo1PXFX98PVol/e7KJJq\ngf0ryvWdDwVpoCzdILXOlFefaxameySChFG9IBX/5pt/yqvPtf7lf1dbv7dLVSd8GHouG2Km\nDrnohlj7AGl9inazLC1Dk3Sz9Cwy/a4ifqt6brQMn4PUOpP+3k5hukciSBjVvziKy/e788VW\naH7P7eux+Z2dw72vTN2eu5SvVECbB+ifonHwm7nUG6MG5/dBiv+KnlfaanJNBqn1qher/u9L\nkLSD/tpWehfHq2hQpS/zxVfjuxLTkUjOXz9/c4J0K0e+i1j+lAOAxYhb+wCdSiBt//zrV43w\nXVqHfOjaSI8FLgtSOvBVaqRhB/21rfQuDjXeXfWRyk7No9mROOuOxJwgPcqDFQMXf9Wr9gH6\nfaROpyh/tEcwVP/oUfWTxsvwOUjtM7U/9zt8RIKEUd2L4319pn9pfYUWow3FqN1P3onOrBqp\nGl7TI8vN71EHMANlSVXdFeNkv2UldK7GA5rDb+ZAZZuvecDi3Wu0LN0gtc7UGbVLhn8rc4LD\nIUjTkqZc9WqKiuTc/mKjL3IpG1SPdFaQypUJ5auvpLFqonGAkXmkIsnvYlye1SSscdMVZnHA\nW/Okqf7EvCC1ztT73OBvlSbtkhwGQZrWCZIeof6uxpWbawTU9fTPZKtaRZSPB+k30SPo1WKE\n7gFU669aTKA/VBesGWy4NA/5Tx/okZgVsPpHv5Jm7TIZpNaZ2p9rrGwYOsHhEKRp7SBVa+2u\n1ZKB33OS3ppDDvnf+4JKv/7U904E6f1Nesi7vimpfYDyo8uj+snX7d2iu+ruS9k/uvy0Dtlc\nC1v9DqbQzemuiSC1zmQ+1y5M90jX9mKJwyBIK3VjsaFtF4ZbCqowASBIKzkIkm6n/V1MXehT\nUIUJCEFayUGQLnXDsjvu7V5QhQkIQVrJQZCqWyVCGA8LqjABIUgruegjvb6LKeC0uzjVj6AK\nEw6CBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQI\nIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACC\nBAggSIAAggQIIEiAAAdBSoDILLjK5YPj4RSAJIIECCBIgACCBAggSIAAggQIIEiAAIIECCBI\ngACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEjAuPv9Pu8bCRIwpkzRvCgRJGDM\nvfF2AkECRtw77z8hSMAIggQIIEiABPpIgABG7QARzCMBLhEkQABBAgQQJEAAQQIEECRAAEEC\nBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAA\nQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBwh7MfiDYVggS4mfxiMqtECTE\nz+KhyVshSIjevfN+1s8IV2AECdGzD5J8W5AgIXoLgmT5/dMIEuJnm4slbcEJBAnxs22pzQjS\nYB/qQ8eKIGEP7MYOJoM0mMyPcSVIOKCptuDg1z/+EEHCAU20BQdrrM/VGEHCIX1sCxIkQABB\nAiTQRwIEmD5UqwHIqB1g617GKG8nh3kkHIXkYlSbBRMECTsiuhjVaiERQcKOiC5GJUg4KNnF\nqMEG6fmVpN95/nNO0ttGp8ChCa/qDrSP9EqTt5/v4m1y2eQUODbpIFn0uE72h196ld+Sdz10\nS5OvV/4qX8ufAgcnfcPe3DHA08lhjZSWP5gkr/JdusUpcHCethM6OW3aJUn91rxrfblh4Slw\neD42uCuadR5qpOLtixoJO3Equ0cOg2T6SLeXfi1/CsA1PcrAqB2wghmtYx4JWOxUjXqzsgFY\nqjF5RJCAhZqTsAQJWOTUWsxAkIAlOmuCCBIwW1a96q6tI0jAXJlJ0qm3RpUgAXOZIA0s9SZI\nwEyZ+q9fHeUECZitrI6y4TuPCBIwT9Z420OQgFnKBL2bdcNJIkjALEWATuZFD0HCHsnf31fl\naDhJBAn7s8Ud51k9WkeQcAzyDy3Ps5b+1wkSdkdwV646MxP7bREk7I5ckKraZ3AStokgYXfk\ngzS9/SNBwv5I9ZGKFBX/zdhGlSBhf6RG7VSQJpt1BYKEPVo/j5SZkYahMbo+ggQMyMw9E6ds\nbHldC0ECBmR5Oc5QNOtmVUkECegztx6p1zOSRJAArdGxao0yzKmSCBJQag71lRXS5zVBHQQJ\nKDUmn+xzRJCAUnM5hNmYYdbAt0KQgEIjSNnHO4+GESSgoGsiFaRPdx4NI0hA6V6Oc9+te0cK\nQQJK9/s7SMWo3ZyldT0ECcjVxNG9yNKyHBEkIK+X1s1b6j2AIAFqaV3x/mQx4t1CkHB4WXXP\nxMlq7qiJIOHoMnMnrO0kbBNBwtGpCqmK0LIkESQcnKqNstOJIAHLZeomvgWTsE0ECcfWaNYt\nHbErECQcWqZ39V6bJIKEQ8v0M/gIErBctnCNag9BwjGZB1merG89GkSQcEjZkvvJPyBIOCSz\n/aPU8QgSjkjtW7d0qfcAgoQjUrdMCB6QIOGA1t0yMYQg4XjKOyaW3zIxhCDheOpdvcUQJBzO\nkn3rphAkHE69MQNBApaSm4VtIEg4GMHJowaChBgtf0bsJjEiSIjRiqeWb5QjgoQINR5lZGeb\nZl2BICE698772TaLEUFChJYGacMcESTExypI1QD3ds26AkFCfCz6SJlJ0qYxIkiIkcWonQnS\nxjkiSIjS3HkkdQPfxs26AkHCnuk7Ybc/EUHCjmWNt9siSIjYRAtPPfPoJLc0dRxBQpzu9/uH\nMYes2jl1xTOPbBAkxKiRoaEgqRskzLPDcgdVEkFCjMocjc8n6SDl8jfwjSFIiJCpikaqpKy7\npff2VRJBQoQIUjCnQMxMm+6e93KUVSMNLka9KwQJMVJV0cCoXVbt6p05mT8yCBJipDPUm0cy\n7bhy+0eCRJAwZXguVldIJ/PaGYKEHTGzRyplLqskgoQdkX14mA2ChP3wlyOChH1QGcqL7tFd\nB2jpzneLECTsQTmwoEfr7nqUgSARJFgqR+v02rq7HmUgSAQJdnqdo8xxjggS9mBglGHx5uDL\nECTEr/dI2OV77C9FkBCz+kZY2UfCWiNICNScWiWr9ttSH/lDkBCkeXtAVqPe5kNvCBKCNGtX\nYlUdbfFIWGsECSGasU9+NjzqvfB8a0cnCBJCNBmkLDOLuwWqoRVPADRcBul1S99vv89Jcvnd\n6BTYiakg6ZUM5faPUmeLJUjPNEny1/tN4bLJKbAbE9d2pmskke0fFz8BsMFhkL6S6+v95uv5\nztRXctviFNiNz60tVSFlQqPekQUpSV76zbuVl6RbnAI78qn/r7d/1BFam6TogvR+kyaNDzpf\nblh4ChxDZu4+EhpuiKuP9JX85fl38aaokT52kggSRuipI9k7YeMatftL0ttffk3fSXqck8cW\np8DeZXe9AaReWyc1BRvVPNIjrdtu39ucArt2L25+LW+AVWsZ3O4B+ZHbCdnfr3ORouv3c7NT\nYL/K6sj1pt4zsbIBscjumcdtgiYQJESiyk3m+qa9OQgSIqGmYIs1QQQpoFMgJplepVqOMtzD\nadDVCBKcWjTOnJmpo3umhu7CQ5Dg0MKZT1UZFWuCMrU7vnzJ1iJIcGjZWhw9B5tVzzwiSOGc\nAh4sXB1a9ZACHPY2CBLcWRYknSO9mCG8DJUIEtxZEaTqlgmCFNYp4MOSPlJ3pXeYSSJIcGjJ\nqJ1e6h3iCEMDQYJTFvNI1V17/rd/nEaQ4M/E3eR6qDuA7R+nEST48rGdl5l9vYPvHSkECb58\nHHnIstau3sEjSPDk41i4qpDqXb2DR5DgyecgtfbbigBBgiefglQtqYsmSQQJvnzoI+k7+AjS\nSgTpCMZH7Rq7eseSJIIEf8bmkRq7ehOkNQjSoUWwtK6HICE4WTXoHUmKcoKE8ERYHxEkBCee\nSdgmgoTK6p3kJayLkbdfgSBBE3i2iYBVOfL4KxAkaAJP21ptZbPO469AkKBIPP9xrZW9I5+/\nAkGC4j9Iq0cZCJKHU6DDe5DWD9YRJA+nQNd0B2PTITGJQW/6SO5Pga6pIa8thsSqCVeZyaNl\nRRT594EgofL5itrgn/tq6UI/Rguv7t6PTR5H6N8HgoR5tuiAmCD1ciRV+804jtC/DwQJ8+gL\nTrJtl6ldvQeadVK134x+X+f9UgQJ85h/2+WipB+tPNA7krq6ZxyHIMGxe+P/ReidVAcWeBMk\nIQQpQHdVH+VSSSoTdDoN3SrhMEj0keCcadUJBSnXGzOMVkku+kiM2sE50YE7k6N8sEpyN2rH\nPBKck6go6geGmdG6wdtgpcY0XN2gRJAwn0BFkVVjDDHeUD6OIMHG6n/gzcNadpYjggSH9JNa\n8lg3ZviAIMGZTM/BxvOwlvkIEpyp+kf7yxFBgjOZmYTdQ5+oiyDBFf00y1NMG6jORpDgiNoR\nP9NPmdgbggRHVHWke0meyyKPIMEN9SxLnSEdJAfLDljZgL0wlVB3DtbBvqjutl4lSNhYNQtr\nltaZLzjY88fdtkIECRszQTpVHysOdqFzuNEdQcK2VIoGltYRJPlSeDgFOrbqlbero+YJO+83\nQJC2PwVaNuuVfxjupo+0OYLk2ha7P1ZvT8OLghi12xxBcmyDNlDjlol8bAbWrjm5qPHJPBIc\n2jBI6s6j9UsZwnig4CiChHyLIOnBut6g92LuujuLECQUxC9THaRT8+M1HA7ALUKQUBBtOGX1\nLqoTGzPM78EQpCUIkntivfIsM7fCTt3BZ5NegrQEQYpXZrZmmL6Dz6o9SR9pAYIUr7IRl+l9\ngj4Gya6SYdRuAYIULVUhZb2l3gNsW2uupoQWIUgQpR8eZgIkGaSgESRIytReQfN2UQ2822OF\nIEFStd/WnHmjwLs9VggSBNX7bc2bgA2622OFIEGQeVjLzCppRwgS5NRrguZWSbtBkCBn04e1\nhN0MJEgQ04qR8HUf+sAEQYKYxsYM4td96EPlBAlCWs8Ok77ug5+8JUiQ0donqHndi7TxCNIi\nBCk67f226uteqI1HkBYhSHGoBua6j4RtBKn18XL0kZYgSFGoJov62z+a616sJmHUbgmCFAUT\npG51dK9bdK0gLe4tlT/IPFKQp8Bq5iETne6RiZB6YT67ok4JvTIqESQspR4Z1m3WdToz9YeL\nezmhd49KBAkLDe/q3e0UVdXJ4t5S8AN2JYKEZdT+JqdukvqXve7bECSBHwnwFFjJ7OrdSdLo\nZU+QPv7IzznPn+fk/E+qQP1TIERVjoaTNHDV00f68COP5P1RmryJJokgBS+rR+s6QRobY2PU\n7sOPXJLf/C8557/JxeYIU2clSMHSqfl059HohM+6eaSwrQxSEYm/5DYjG60jEKRYVXOwA8+y\nPDSBIF2Tx6wgJW3CpYID9a7evksSnNVNu79Hkuazmnb/UoIUl26LqnoIHznqWT/YkCTfRW3z\nmP7B1zW5PMsjDJ11dsrgRr+Pr4NEs27A6uHvtOgh5effWT/6myTFN9JHikFv1Fk37KYe1nJM\njidkn5fk+iJIMejMg2Z6W+/ph7Uck/OVDd9J+iBIEWgHKdM745/Mo8rRtjpIj2s5cvec/eN/\n5+k+EEHyrxkk86yWXC+tI0k9a4N0UUMDSTo/SfkXQYpB3Ucyj5jIzbPDCFLPyiD9JJdXEYuf\n5EusSDlBCkI9aqcewrflLqrxWxmkNNFDB7Ij1gQpCHoeST87zIzWEaIhAisbCNLOqSBVo3WT\njyoPfl3cFlYG6axrpGLhqiCCFIiyQ6TezHwIXwwrtbcg00d6pMmPWJFyghSKzDzoaPZoXRT3\nDm1h7ajdVa/psbmLwvIU8KWsfdSo9/0+a7QujrtZtyAyj5Rc560QWngKeKKnjyxG6wjStj8S\n4CkwyYzW6VxkM9JBkLb9kQBPcWS9gbXhkTY9yFCvEpqdpOPlaP3wd+VyEysUQdpSb2BtZKRN\nDS+0HsK34OBHIRekpLjBz1+pMFev0hi4XyLXo96n8r1Vg415pEU/8pUWd/Q90uRffk3E6iSC\ntJ1eKnqfqG+EPemPDttgm29lkG7JX/n+L7nkL7lJWYK0nZlB0ktUc1UlHbXBNp/AEqHqhdwy\nIYK0nckgqeqoM+h90AbbfKsXrZoaKSVIkZjqI/VWBLFKdYbVTTvTR7pZbhIpXSrMNTFql5W1\nD+GxJHJjn1oilMittyNImxqdRyoro+IlG5zYEloiVFRLxbZcQgiSF+WSoHs5WncnSXZY2QAt\nq/bbUh/BBkGCkpmGnRr1pkqyIxWkf9e1JZk8BTbV2UaVINlZG6TbJrsMEyTHTLOOUe+lVg9/\nGzP2/l52CmxP71pXVkcZqxiWWD0h+5tfkufzwhP7YjG0RqG5q3fGKoYlBJYIfb9roz/Ze80J\n0lYGV83Vo3Xs/riUQJAexUQsfaQ4DK7jLoN0YhvVNVYG6fpu2j2Tc/6PIEVh8M4iRhkErH/Q\nmF4mxJbFcrbrpAwFqfWwlu4TyuktzbR2+Pu7+Ogrkbunr3+KfRq/Rre8+Wc4SCPVEXchWWBl\ngxefrtFNb0ftH7x+JGy3Scd9sRZWBCm9/vx7iRamd4rd+nCNbruj1dCDYU/VK4cF2ZsVQVIb\nnly/HxaPRrI8xV59ukaXX7/zejTmERP6w/FRBoJkY0WQXv9+vs5qWcP19vvnuVRR2SJIdj2a\nzCRp/AnlBMnG6j7S38/XhbV2UzqVxcdrdGHXxObH6upnPEf0kazIDDY8LgTpg35l8ekaXTZY\nZlN/1HcefcoRo3Y2BIL07/tMjfRRPzafr9E5nZ3u91gFqd7V2+4kGLUySM+f8rku52/RNas7\nC9LgNb7uGh2r4+YEKcvqKglS1gw2PG5pMW739RAfBD9AkCQOOdDvmpEj9ZSJ/MRSIEkrh7+v\nP/Jj3zlBWnLEuT0aPc6Q6d29IWNVkG4bzcfuLEjyw1/D0ZzVWlT38J1OBEkWNZID4sNfd52Z\nBYcsn2ap75ggSXJW9ZG+6CPNJDv8VR7t/Z/dMcvnG5U9pBNtO3GM2kXonqsk2QTpXj4pzPaZ\nsJiLeaT4qFadZSWnYqObddwJK46VDfFZMgp4L2PEnbCbYa1dfKogWdRJKkfV9o8fqyPWMyzA\n6u8I6bZdbjEUqBYFzamOWGG3CPcjRahxrc+84NUTycvRhsmDN95iNu6QjVLd+pp738T9fjq9\nczQZJO5CWoY9GyI184JXwSnyc3qHL5teqUqQliFIkZp3wevgmFsmsuk5WIK0DEGK1ay+TGb/\naGX6SIsQpFjNGV2rlgJZTB4xarcIQYrX9HyPCpLtJCzzSAsQpB3LGm+xLYK0Y2WETgTJBYK0\nS9WG+KcTVZITBGmPMjPMUC31xsYI0h5l6lktY7vjQx5B2iE9WDdrqTdkEKQdsp6EXYhx8hpB\n2h9Ho97M3DYRpP1Ro3XbB6nxFgRpR+qq6JRvXiWxurWFIO1HVk0fuRj1JkgtBGk3MjN7dHIy\n6k2QWgjSXqhnteR6v63yM9uekD5SE0HaiUwHyd1+W4zaNRGkfejXR6VtZ3qYR6oRpH3IzEOP\nmp+kznCHIO2AelSLilEzSa1eDNXHpghSIFZc56pjVLTq2kFqjqtROW2MIAVh1XWuKiRTH9VJ\nagWp9Rnr0pHAKQQpCGuu8/IZE5m6g28sSGsmfajM5iBIIVg1uTn+zKM6nquCtLxoB0KQQrDw\nOtcrgorukVlmNzJqtyJIrGCYhSCFYNnFalpyH249qno3y6sVgjQLQQrCoutcB2njh7UQpFkI\nUhAWXOdZpt+UaxnuZaPu0wEWD73RR5qDIAXC+pGw+obyk86R6h1tcbkzajcHQYqSnjrqNeu2\nudyZR5pGkKJk1qie9DAdHRnfCFKM9OK6xjaqdGQ8I0gxqteomiTRkfGMIEWotdTbjHjTkfGK\nIMWnuybId3mQE6QYmcmjzt1H8IkgRafe1Xv6GeVwhSBFh3ZdiAhSLExkiFGQXAbp9ZUkl4c+\nyMejEKSeqhl3On3+RnjhMEivNClc1UGOEKSPI9J2w9VV9UOMwuQwSLfk552mn/RSHmT/Qfo4\nR2o7gZrpKokcBcphkFL1g8/0/DxGkBpv7b7Y+DadtXJpXUazLlwOg2Sy87pchoKUNC08RUg+\nriMtP3mfat3V1ZZ+mKVg8SDLYZDOycu8uhygRpoK0ozGXVVtvSujd+gYpAuYwyD9JF/61TO5\nECRdKc06Qvb+vtMpu5OkYLkc/r5V6XlMtN72EKSJPtJdf/HTuF6u2nTvquiu1gQRpGA5nZD9\nu5pXz68DBGnWqN3nIJX9orvOUU6VFC5WNmzo02CCGWn43LYrgnQ3GzNQJQWMIDlms9NcWQG9\nqyEWBYWPIDnVaO7NmZLNyrE6chQBguSU3ROLytgUd8KW38ktsCEjSC7ZbfajclS8v8+YvYVX\nBMmluUEyW+KrNUGZ9WZbhM45guTSzERkZlcTM1g378eqswz1vsjWtgiSUxNDdfpqN0EyS1Rt\ng9T/Znbr2hpBcmrWJK1eodobrLPKUfvb2T9yawTJsQ9NLHO1t6sj9UP5/AplIEjsaLw5ghSM\n+naJPM+7t0xYdHEIkg8EKRhmaV15I+xpxdRrvx1HkDZHkIJxVxszFIPdevvHpQfqNwTpI22N\nIAWjWNud1c/gW7MaqNcQZNRuawQpFJnemOF+6sweyWAeaVsEKRBls65colp9xmdxYIkgBUI/\nzHLpUm8qHM8IUhgyFaTTsr2C6AJ5R5D8q1Yy6AeV2x+BQTnvCJIvVWMsM9NHJ/3eOklME/lH\nkPxoNMZUkLLMNOvsqySC5B9B8qNujOlVQWtuKCdI/hEkLxqX/thab/vDkSOPCJIXemt8s9R7\n7fYmjNp5R5C8MJ0jHaRqieriSVjmkTwjSH6onfEzmQoJ/hEkP8rFQOoBE6uXqCIABMmTLDMP\nD9Mf+y0OViJInqiFDMXskf7Y7enpUwkjSH5kWWfI22nviFE+cQTJj8wsBvLSpGPeSRxBcqXV\nmCoroJPqJnkoSuc91iNIbnSeq6xH67J61yCnhem8x3oEyY2qMZXphXVqV28/VRJBkkeQnKgv\n3ay7sM7HJCx9JHEEyYkqSL2VDH4GGxi1k0aQnAgsSMwjiSNIbqj7Je7VM/iwMwTJjXux++O9\nfgaf9+JQIckiSItZXYyZ3khVP4Nvs0LNQhdJHkFayO5iHHsGnx8M2skjSAvZXYxjz+Arj+G6\nbmAaaQMEaRnb55P3po+qA7lvZhGkDRCkZSyDlOnqqN+k89DMIkgbIEjLWF2MVQ+pHyQvFzV9\nJHkEaSGbi7Fcoqpmj7pJ8hMkRu3EEaSFLC7GTI8yDFVJEkFaMFrBPJI0gjRt5KqbcTHq5yp/\nWhK0uplF9RIEgjRl+YWafRr0Xn94c4DGW3hDkKYsvlCrbetO9WcGT7CqOmEILgwEacLyC1Xf\ntbfxUm+CFAaCNGHxhaonYU/T37kKQQoDQZqwPEjm4WEbo48UBII0ZeGFWjbjTg62ZGDULggE\naYr9hapmXsul3k72NmFSKAAEaZrlhar3CcrNJKz3u/jgAEESp4N0Mk+ZmBkk6pWoESRp/e1N\n5iSJnk7kCJI0FaSJSdgext4iR5CE6ZEGy0lYZoNiR5BkLdwniCDFjiDJKmePcoJ0OARJlNlv\na2GSyFG0CJKoz7dMfMCoXeQIkgSTmjUrvZlHihpBEpCZJG291BvBIkgCTJDI0XERpHFzG1tm\nN2JydGAEacz87n97YwYcEkEaM3tAOmu8xVERpBHzp0jLOVgeHnZwBGnEjCDVVVEAzzyCXwRp\nxHSQ9FideehRf1tvJoYOhCCNmewjVds/Dj48jKUKx0KQxkwlQQ96j915xOK5YyFI4z63zVSQ\nxhYFsZz7YHVLeosAAAurSURBVAjSQnqkYWz2iCAdDEFaSN3BNzpUR5AOhiDZM885On26E5Y+\n0rEQJGuZ3mNreNDbdK0YtTsWgmRNP65lcNC7GSDmkY6EINnS222N3QlLk+6YCJKF6nHK43fC\nMshwUARpPtM5YrQOPQRpPvMsy0+jdQTpoAjSbKpNNz5ap9BHOiaCNJu5o7z+aAjD3sdEkObS\nm3pP71vHsPcREaQZZozW4eAI0rTGaB03wmIYQWoabpVVa4K4oxxjCFJtZJygNVpHkjDIaZD+\nfV+TwvX2b6tTrDEycq17SAt3x8cxOAzS65zULpucYpWBuVRzwwSjDJjgMEi3JP39K189H2ly\n2+IUq/SD1BxlID74xGGQ0uSvev2XpFucYpVekMraZ/IOPiB3GqQkGftAf6Zh4SnW6faRyrV1\n1SgDScIH1Ei1zqhdkaHGaB1Jwgdu+0iPZ/kqzD5S3plH6qwJIkj4wOXw96XRdju/NjmFoKy3\nJmh8xTer6w7P7TzSrZxHSq/fQc4jtZnRusmKiPXeYGXDqEyvCZoRpMZbHNVRgzTZGqtH6yaS\nxD2xyI8apE+tMb0X8ew1QQQJ+WGD1HjbkdmuCSJIyA8apPFr3+xvMm+0rnkUcnRsBKnJPDvM\n6hHljNqBILWp2qfeRnXu8YjR4R0ySKN3Huk7YVnDAFsHDdLIvbCZro5IEiwdM0jDrbGyXXeq\nbuZzjyZivI4apAEqR43bYh1j0CJmBKlSzMGafYI+BWmreoNh9JgRpHotQ2/uqBOZ8sPN6g0m\ndqNGkKq1DN01QZ3ItJ7CR5DQFnOQhNpYA2sZ1OEbb3sfyl/uBClq8QZJqo01tpahc2Gbd+0P\nJdFHilnEQWq8XUMFqb+WwX2QGLWLWLRBkmoJjT7KcjhIm9YbzCPF6+hB+vCQieE+EvUGhhAk\n1T0aCtLgqB31BoZEG6S1bazMvD2N7/84NI8EDIk4SKvaWGOTR8Ai8QZpXQ0xNnkELBJzkFYw\nk0eWd/ABIw4dJGIEKccM0ujkEbDMIYPUnjxiMA7rHTZIZvKICVZIOGKQWpNHLBWFhGMGqZ48\n4uYFiIg+SPY9nNbkEUGCiMiDZNHDqcboWpNHBAkiYg9S4+1nmUlSZ/KIPhIkxB0ki/rEBKk7\nCcuoHSQcJUgf1gQxj4T1jhUk1gRhI3EHaVYPp97wkTVB2ErsQZru4ZhnhxXNuhlJoqGHJSIP\n0vSFX95p9OGG8u7RcoYesED0QZqiK6QPN5Q3MRiOZfYepPKRRzyhHFvbfZBynlAOB3YeJP3s\nsLnDdQQJC+07SNWzw+ySRI5ga7dBUrtt6d3x5weJUTssstcglUN15hl8FkliHgmL7DRIalgh\ny6rROhY1YFN7DVLxf/1nWQIb2WeQdHXEGlW4stMgFf+xbx3c2WWQzP6PJAmu7DFImV5bR78I\nzuwySO8cFaPYVElwZodBKmaP1LwqVRJc2WWQGPSGa/sLEjmCB7sLUjHKwCJuuLa3IJVzsAQJ\nru0sSHotA3dDwLFdBalaE8TdEHBsT0FqbY5PjOBS/EGqnzKxQUGAeaIPUrU5PjmCR3sJUvdh\nLbTs4FTUQcrqp0y0Ps9YA1yLOUhmV+/eHXyMfsO1yINk7jxqYT4WzkUcpKyxb13rCwQJzkUf\npFP/KRMECc7FG6SsGmzIh5NEjuBOnEHSY3XZ6TQcJEbt4FiUQcrMRqqjdx4xjwS3Yg7SSb2e\nQKbgQIxB0k8oP9UfjaOVByciC1Jm3lS7ek8FybylYsKW4gpSVq0KmrkxQzUSTsWETcUXpHIO\ntjcJO6I9pUSQsJWogqTqormjDAUdHVMXkSRsJMYgzRpl0NpVEUHCRmIKkn6a5an54ZR254gg\nYSOxBcl++8dyuI4+ErYVUZDMnbBLtk9l1A7biitIp/nDDF3MI2FL8QSpcUM5W3ojNDEFaeZi\nBsC9aILEI8oRsmiCxLZ1CJnLIL2+kuTy0Af5eJTeF9n+EWFzGKRXmhSu6iBWQSJGCJzDIN2S\nn3eaftJLeRCbIJEjhM5hkFL1g8/0/LQKEs06hM9hkEx2XpeLTZCIESLgMEjn5GVeXeYHaShH\nrFJAaBwG6Sf50q+eyWVmkIaadXbr5ggdXHA5/H2r0vNI5gVpsFlns5Kbxapww+mE7N/VvHp+\n9Y6SNKlPfcjRzCRx+wTcCHhlw8honU2QuKEPjoQbpLHROoKEAPkI0uf+kT7F+Ki3RXONIMGR\nQIP0aRLWZgCBPhLcCDNIE5Ow84e0GbWDG2EGaewbFkwKMY8EF2IKEtULghVVkBpvgaCEO/zd\nwxAcwkWQAAEECRAQUZDoIyFcUQWJUTuEKqYgMSmEYMUVJCBQBAkQQJAAAQQJEECQAAEECRBA\nkAABBAkQQJAAAQQJEECQAAEECRBAkAABBAkQQJAAAQQJEECQAAEECRAQaJCAyCy4yuWD44L/\nYnsvAQUIoAS1kMpiwX+xvZeAAgRQglpIZbHgv9jeS0ABAihBLaSyWPBfbO8loAABlKAWUlks\n+C+29xJQgABKUAupLBb8F9t7CShAACWohVQWC/6L7b0EFCCAEtRCKosF/8X2XgIKEEAJaiGV\nxYL/YnsvAQUIoAS1kMpiwX+xvZeAAgRQglpIZbHgv9jeS0ABAihBLaSyWPBfbO8loAABlKAW\nUlmAaBEkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQ\nJEBAxEH657fsP+ckvb08nfyWejx5wetvb3i+AloCKoqlV+q17LfyqQWpn2vpUp787OXcJa+/\nveH5CmgLqCiWrkseviHmL/l6X0U/yZePk/9L0r/8L03++Th5wetvX/F7BXQEVBQ7v4ueYiPm\nqk7upwy35JEXf4FvHycveP3tDc9XQEdARbHyTC4h/Bn9lOGaPPOiWrj6OHmD1/8FArkCjICK\nYuWSPAP4M76Si4/TJiFUCN5+ey2MK6ASUFFsfCe/3q+jvOglPHycNpAgefrtlUCugEpARbFQ\ntmr8/xmfqZ/GVRhB8vXblwK5AmoBFcXCuRh49f5nfKWemjZBBMnbb18K4wpoCKgoM+hHTn+V\nbQovf8bmQ68vviZy0hCC5O23L3i8AkYEVJQZ9GW85jHuMiV4e54vT+enV9So3dPnqJ3H377g\n8QoYEUxBbITwZ3x4HLL6Lv89fiQ3byXw+dsXQrgC2oIpiD3fsxj+Tu59ZYPX374WTowI0lJf\nXv9FPJen9ncx+/3tK77P3xRQUWx5/TP6bVq8ytXfXk5dCqRh5fv8TQEVBYgXQQIEECRAAEEC\nBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAA\nQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEECBBAkQABBAgQQJEAAQQIEECRAAEEKXPVU\nus7j6cY/fGxdJAwgSIGzDtKZ/0l94K8euLEgjX1bUA9WPRD+6oEjSHHgrx64TpB+zkn6U394\nS5Nb+fL93y1Jv/UDxz2V9cj4mweuHaRrmZOL+fBSfPSlglR+6YcgecLfPHBJLc8fyeWVvy7J\nQwXpkaR/+V+qgvT+yk9ypmnnCX/1wLWCdE1e70+9kqvKy7VI1DtOZZD+5blp5ME9/uqBazXt\nGplqZKZ+SZC84a8eOIIUB/7qgesEqfV5ghQO/uqBawVJdYo6Hz4IUgD4qweuFaTfYpgu/zGD\nDa1RO/M9SfL0WNzDIkiBa88jlRNHSfpsfdgO0vn9dW+lPS6CFLj+yobk61l9eEuTy792kP6d\nCZIHBCl+5UoH+EWQIpYkv3n+uiY33wUBQYrZt+oh0ZILAEGK2c8lSc7URyEgSIAAggQIIEiA\nAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAAIIECCBIgACCBAgg\nSIAAggQIIEiAAIIECCBIgACCBAggSIAAggQIIEiAgP8P760Uus+6sgAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"Explain Wages with Height\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Matrix of predictors (we only have one in this example)\n",
    "X_mat <- as.matrix(X)\n",
    "# Add column of 1s for intercept coefficient\n",
    "intcpt <- rep(1, length(y))\n",
    "\n",
    "# Combine predictors with intercept\n",
    "X_mat <- cbind(intcpt, X_mat)\n",
    "\n",
    "# OLS (closed-form solution)\n",
    "beta_hat <- solve(t(X_mat) %*% X_mat) %*% t(X_mat) %*% y\n",
    "print(beta_hat)\n",
    "# 2.945535, 0.951942\n",
    "\n",
    "# Plot\n",
    "plot(x=X, y=y, cex = 1, col = \"grey\",\n",
    "     main = \"Explain Wages with Height\", xlab = \"Height\", ylab = \"Wages\")\n",
    "\n",
    "# Draw the previous regression line\n",
    "abline(a=fit_params[[1]], b=fit_params[[2]], col=\"blue\")\n",
    "# Current regression line\n",
    "abline(a=beta_hat[[1]], b=beta_hat[[2]], col=\"green\")\n",
    "\n",
    "# To get y-hat:\n",
    "y_hat <- X_mat %*% beta_hat\n",
    "points(x=X, y=y_hat, pch = 2, col='yellow')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we mentioned the assumption that the matrix $X'X$ is not-singular; below are two examples when this will fail and we cannot calculate the inverse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>100</li>\n",
       "\t<li>2</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 100\n",
       "\\item 2\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 100\n",
       "2. 2\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 100   2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>2</li>\n",
       "\t<li>2</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 2\n",
       "\\item 2\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 2\n",
       "2. 2\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 2 2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "NULL"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "NULL"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dim(X_mat)  # 100 by 2\n",
    "inv <- solve(t(X_mat) %*% X_mat)\n",
    "dim(inv)  # Possible to invert\n",
    "\n",
    "# 1. We have a column that is a (perfect) linear combo of another\n",
    "X_mat_fail <- cbind((10*X_mat[,1])+8, X_mat)\n",
    "inv <- try(solve(t(X_mat_fail) %*% X_mat_fail), silent = TRUE)\n",
    "dim(inv)  # Couldn't invert\n",
    "\n",
    "# 2. We have more variables than observations\n",
    "X_mat_fail <- matrix(runif(100, -5, 5), nrow=10, ncol=20)\n",
    "inv <- try(solve(t(X_mat_fail) %*% X_mat_fail), silent = TRUE)\n",
    "dim(inv)  # Couldn't invert"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Linear Regression with Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we can also use an interative method known as **Gradient Descent**, this is a generic method for continuous optimisation. With GD we randomly initialise $\\beta^{GD}$ and then calculate the residual (error) and move in the opposite direction to the gradient by a small amount proportional to a parameter we call the **learning-rate**. GD is a bit like rolling a ball down a hill - it will gradually converge to a stationary-point. If the function is convex with a small enough step-size (learning-rate) and high-enough number of iterations we are guaranteed to find a global minimiser. **Stochastic Gradient Descent** is usually used for neural-networks to avoid getting stuck in a local minimum due to a non-convex cost function (along with other methods).\n",
    "\n",
    "The general-formula for GD:\n",
    "\n",
    "1. Find a cost-function\n",
    "2. Randomly initialise your $\\beta$ vector\n",
    "3. Get the derivative of the cost-function given $\\beta$\n",
    "4. Move the $\\beta$ vector in the opposite direction to the gradient\n",
    "\n",
    "In the case of this linear-regression:\n",
    "\n",
    "Our cost-function is the Mean Squared Error (MSE) which is:\n",
    "\n",
    "$MSE = \\frac{RSS}{N}$\n",
    "\n",
    "and:\n",
    "\n",
    "$RSS: y'y - 2\\beta^{OLS'}X'y + \\beta^{OLS'}X'X\\beta^{OLS'}$\n",
    "\n",
    "The derivative of the RSS is: $-2X'y + 2X'X\\beta^{OLS}$\n",
    "\n",
    "This can be simplified to: $2X'(X\\beta^{OLS} - y)$\n",
    "\n",
    "So, we can write our 'delta' as:\n",
    "\n",
    "\\begin{equation*} \n",
    "    \\frac{dLoss}{d\\beta} = \\frac{2}{N}\\sum_ix_i(x_i\\beta^{OLS} - y)\n",
    "\\end{equation*}\n",
    "\n",
    "This means our equation for $\\beta^{OLS}$ becomes:\n",
    "\n",
    "\\begin{equation*} \n",
    "    \\beta^{OLS} = \\beta^{OLS} - \\frac{lr}{N}\\sum_ix_i(x_i\\beta^{OLS} - y)\n",
    "\\end{equation*}\n",
    "\n",
    "```R\n",
    "for (j in 1:epochs)\n",
    "{\n",
    "    residual <- (X_mat %*% beta_hat) - y\n",
    "    delta <- (t(X_mat) %*% residual) * (1/nrow(X_mat))\n",
    "    beta_hat <- beta_hat - (lr*delta)\n",
    "}\n",
    "```\n",
    "\n",
    "With learning-rate set to 0.1 and epochs set to 200 we converge to the same result: 2.95, 0.95. We can track how the line-of-best has been gradually fitted with this method by plotting it at each iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gradient_descent <- function(X, y, lr, epochs)\n",
    "{\n",
    "  X_mat <- cbind(1, X)\n",
    "  # Initialise beta_hat matrix\n",
    "  beta_hat <- matrix(0.1, nrow=ncol(X_mat))\n",
    "  for (j in 1:epochs)\n",
    "  {\n",
    "    residual <- (X_mat %*% beta_hat) - y\n",
    "    delta <- (t(X_mat) %*% residual) * (1/nrow(X_mat))\n",
    "    beta_hat <- beta_hat - (lr*delta)\n",
    "    # Draw the regression line each epoch\n",
    "    abline(a=beta_hat[[1]], b=beta_hat[[2]], col=\"grey\")\n",
    "  }\n",
    "  # Return \n",
    "  beta_hat\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      [,1]\n",
      "  2.945535\n",
      "X 0.951942\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAOVBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD//wD///+lBszRAAAACXBI\nWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2diWKjOrOEcZxNySQOev+HHQPaNyRokID67j02\nNhhEfmq61GpExwEAq+lqNwCAMwAhAUAAhAQAARASAARASAAQACEBQACEBAABEBIABEBIABAA\nIQFAAIQEAAEQEgAEQEgAEAAhAUAAhAQAARASAARASAAQACEBQACEBAABEBIABEBIABAAIQFA\nAIQEAAEQEgAEQEgAEAAhAUAAhAQAARASAARASAAQACEBQACEBAABEBIABEBIABAAIQFAAIQE\nAAEQEgAEQEgAEAAhAUAAhAQAARASAARASAAQACHN05nMb5v1neCl637Ghcdzo8e49NN1L8vb\nOtuut0ST9Lf5p+F9+7a4jYcGQppnQyF9dN3HuPD13Ojb/ooU0YZ/ty7RJAIhyQNcjouedhEb\nCulf193HhbfnRtO/5feu+7e8rTPtkk3ZTEjzf6GTctHTLqLk4ii9kG7C0SmRPj3erWgPRawR\nUnqPxT88GRc97SLsi+O3617HhSl0DCs/n12dt19r26/X59LL+6/x3fD2/fzR24+xt2ck+uJj\nZHqZItG3iEzmDvjvc7P7p9r54/3W3cSqx8dzl93rl96l2JHsdf2Ofa7xpzqoDq//7qrV3omG\njiS/sxtj7Skvap+SS550Ic6l8Tl1Yt6nK364iser52aK5i6d4D9uCul9+tJQktDNc8Vz6Z2P\nyvp2d/BvWryLHf3e9Cq5LBziwMe0o+f+Ro1+jc31hfRutNo70cCRVDfLboy1JwgJJHAvjftw\nrX8JC6Z7Ty9626fW7o9RHq/cFJLASGw9ppW3569fxj1MXs/ewU3/lBsfh+NPEe3xbNKn3OUU\nglSv67nqNyQkvy2+kMwjud+FzgpCAgm6zr52nlfq7dP4d/r2XPi+TZFk2uRlvHjlJ33J3b6n\nC9zY+WgQf4cg8j786GeKLdYOvqYfft1MlT4mCXXThg8zZ/4yalFr5sVqgzylodXvdoesc87U\nOpL4rdMYd08XlRGElIMrpOECU/+Wd9NVJj2afSG5Qhqk9rA3GbPdH8Oqf0NU+TBCi/zpq8iM\nf02/fFX5idcpPrx92+39GJo0ObB/Q7s+eEhIgba4QrKOxNV3ZmPcPUFIIIonpOHffPlPufrO\n+pf/Gba+3u8qJiRSz6MRkzHkLoyYvYObPoRty26jaDpXS7+Dpp8h4kvFuWgb0kKyjiS2dRrj\n7glCAlH8i2O4fD+clZZovl7s69Hc0tnd88oUfu4+Lk0CNXfgH8LY+bu81I2swctzJ8N/Q8/r\nZlmuWSFZS56s/POFkAQXPe0ivIvjMRiq20OufBhbdbIj0b28ff7kCOl9zHwPsvwcE4BDxs3e\ngRMEbvbvH19Thu9u7fJbRCORC1wmpFtgLSJSmIuedhHexTHlu1UfaezUfJsdiRfRkcgR0ve4\nsyFx8aOW7B34fSSnU8S/7QzG1D/6Vv2keBvSQrKPZH/3Fd4jhASiuBfH8/q8/dz0FTpkG4as\n3Sd3pJMVkVR6TWSWzW2mHchEWafC3ZAn+xqD0IvKB5jpN7mj0fOZOxzeHtG2uEKyjuRk7brw\nWckDXA4IaZ7OhE+9miGQvNgrjb7IfTRU37csIY2VCePSW2dUTRg7iIwjDUp+NuP+qwZhJe8i\nYA47fDcPehNf5AnJOpL3XfCsbp3dkssAIc3jCElkqD9UXtmsEZiup39SW6qKiMeF9NWJDLoq\nRnB3MLk/VUwgPk4XrEw23M1d/hM7+u5kBaz46VtnRpdZIVlHsr8zKhtCB7gcENI8tpBUrd2r\nKhn4eulu72bKgf88L6jb28+07YyQnhuJlLe+Kcnewfjp/q1++Xh/OrpX0X0Z+0f3T2uXZi2s\nOgfZaHO4a0ZI1pHkd3Zj3D292sUSlwFCWokriw3ZtjC8kKYa0wAQ0kp2EJLwaT93GQtr0lRj\nGgJCWskOQrprY+nmvfenqcY0BIS0kh2EpG6VaCEf1lRjGgJCWskefaTHxzAEfHOLU+vQVGPa\nAUICgAAICQACICQACICQACAAQgKAAAgJAAIgJAAIgJAAIABCAoAACAkAAiAkAAiAkAAgAEIC\ngAAICQACICQACICQACAAQgKAAAgJAAIgJAAIgJAAIABCAoAACAkAAiAkAAiAkAAgAEICgAAI\nCQACICQACICQACAAQgKAAAgJAAIgJAAIgJAAIABCAoAACAkAAnYQUgfAwVhwldMLp8IhAKAE\nQgKAAAgJAAIgJAAIgJAAIABCAoAACAkAAiAkAAiAkAAgAEICgAAICQACICQACICQACAAQgKA\nAAgJAAIgJAAIgJAAIABCAoAACAmAOIyxvA0hJABijCrKkxKEBEAMZrzOACEBEIE57ykgJAAi\nQEgAEAAhAUAB+kgAEICsHQAkYBwJgD2BkAAgAEICgAAICQACICQACICQACAAQgKAAAgJAAIg\nJAAIgJAAIABCAoAACAkAAiAkAAiAkAAgAEICgAAICQACICQACICQACAAQgKAAAgJAAIgJAAI\ngJAAIABCAoAACAkAAiAkAAiAkAAgAEICgAAICQACICQACICQACAAQgKAAAgJAAIgJHAGsh8I\nthUQEjg+BY+o3AoICRyfgocmbwWEBA4Pc96zfkMcwCAkcHjKhUTvBSEkcHgWCKlw+3kgJHB8\nSnWxxAvOACGB41Pq1DKEFOxDJTpWEBI4A2W5g1khBZWZlCuEBC7InBcMrk/+CEICF2TGCwYj\nVjqMQUjgkiS9IIQEAAEQEgAUoI8EAAGyD2UZQGTtACiFjTLitnIwjgSuAmUxaknBBIQETgRp\nMWpRIRGEBE4EaTEqhAQuCm0xarNC+n3rbh+cf750t/eNDgEuDXFVd6N9pMete/L5Mbx2900O\nAa4NtZAKelx/5btfepW/d8849H7r3h78MS7THwJcHOob9nJzgH9/C4659Cq/jT/susf4dtvi\nEODiVJpO6G9R0n3pVd51+lW+WasNFh4CXJ4aE9wNOqoQkYbXByISOAl/o452FJLsI70/xDL9\nIQDYm0lGewoJWTtwPqSOdhQSxpHA2fhTOtpTSE0dAoDVaBlBSAAsxdQRhATAIv4sHUFIACzB\nlhGEBEA+vVpydQQhAZBLL5X05+kIQgIgFykkX0YMN/YBkEk//RcMR4hIAGQyhqM+EI44IhIA\nuUy2rg/pCBEJgExGHT1tnaukQURINgCQySCk0db1TjgadQRrB87IohCRROnIVtJo62DtwCnZ\n4o7zXmfrtJCkrUNEAmeE/qHlvLewbd1wJAgJnA7CWbl6VRQUzHprWwdrB04HnZCG6DMuRAZh\nRVSCtQNnhF5I4UFYaeswjgTOCVUfaVDR8N+srUNEAmeEKms3CWnW1iHZAM7K+nGkXmYaer8o\nKGDrYO0ACNDLeyb+er8oyLF1SDYAEKHnY55hsHVOSPJtHawdAGHkrUeDcsyQxGO2DtYOAIHR\nsbKyDGZICtk6WDsAFGaqbwxIqZogy9bB2gGgMAaf4jriMVsHIQEwYJZDyIkZwlkGz9aZPy8A\nQgInxBBSH7nzaEw1BGwd+kgASEQkmoTk33mUsnXoIwGgYOMQLIv0juZsHYQEwMgz/gxTBKVL\nvYO2jnP0kQAY6aeJtnJKvT1bhz4SABOytG6u1DucrUMfCYCRXpR6/82UekdtHYQELk+v7pn4\nc8eOsmydWiwDQgKnopd3wiYGYWWWIWjr0EcCQAYkNau3NwgbGzxCHwkAzRSN+r8/X0h+qXfU\n1kFI4OL000188UFYc/AoaOuGVwgJXBvD1jlzemfbOuOtBAgJnIdezOodUlKurUOyAVyeXjyD\nzxFSqCYoauuQbABXpw/VqPJ0qXfQ1kFI4JqIB1mqWb29DpJX6h21dcMrhAQuSR+9nzx6B1/U\n1hlvJUBI4PjI6R+ZR7mtQ7IBXJVp3rpEqXd8Ygbf1iHZAK7KdMtEKBzFJ2ZI2joICVyQ+C0T\nS2zd8Aohgesx3jHh3zKx1NYNQEjgeuhZvT1bNz8xQzhbh6wduByxeevMmqDAHXxRW4esHbgk\nemIGoybIKPUutXXI2oErkpwdny+0dRASuBiBwaPCiRlcWze8QkjgiDB5GRcTkpFzB1/A1nHj\n1bN1AxASOB5M176VkppGlacHj5LZOggJHA87JhQwb+tm7uAL2rqFjYGQQF0Wj91k2DoetnXM\nEpDXAKS/wQFZKiRPR3zZxAy+jmDtwAEpElIvF7a1dRASOB4FfaReKinzYS1LbN3wCiGB41GQ\ntZNCCto6TmPrBiAkcERYZvJ7uoEvZusYXbYOQgJnRtwJG7F1nMbWWV/lAyGBoyCmCoplGYhs\nHdLf4ODMOLzpmUd/vffMI1pbh/Q3OC5CEZGcw/T8PfHsMO9O2A1sHYQEjoihoZCQhhsk9LPD\nmB2SyG3d8AohgSMiPJladhBC4v4NfCFbx7SOxhdDR3ZUiti6AQgJHBB5qUdCUu9O6e1Mo+oM\nwvK4rePmWypbByGBA7JYSHZNEJGtG4CQwAGR8UKJQNOrTENg2jq+ZGKGWVvnLOYCIYHq6Kve\nydr1albv3k3WJSdmWGPrMI4EjorQkDeONPk4Mf2jl/Uut3XM15HVBPUdhASOSXgsVgQkMf2j\nNW0dt7MMebZOHUwd1lvAOBI4HXL0SCS9jVwdycQMEVs3vEJI4EQkHh5m2rpATdAaWzcAIYHz\nENNRxh18a2zdAIQETsGkIT52j7yKIMY3tXUDEBI4A2OKQd7B5yS+/Zogx9YlSr2zbN0AhATO\nwJitC2QZsku919i6AQgJnACvc9Tvaevs77OBkEBrpLJ129s6jsoGcArij4Rdcwdftq1DZQM4\nOvpGWO+RsPm2btrVOlsHIYFGCdcA2fRqvq2pJshT0i62bniFkECThCq7fVTWW2Tr6ti6AQgJ\nNIl9YUeYwlHokbC72roBCAm0iHdJ+/TxrDfjpbZOFhDZB820dQMQEmiRWSH1vbgRNjz9Y5mt\ns3xkua0b2FNIj/fb8/XjpevuXxsdApyEOSGJSoZx+seAjgptnYhb9hHLBmF3FNLvrev44/ky\ncN/kEOA0zPSRehGRQtM/+raOaR2NL76O5BahcDRr6wZ2FNJb9/p4vrz9PjX11r1vcQhwGtJZ\nuykg9X7Wm5fOt2UHo4BmMmzdwI5C6rqHeHm6vO62xSHAiUiNI4npH/tQ1psbd/DlZesCvaMi\nWzewq5CeL7fO+OCsNlh4CHAN+ildJyXUGzIqmW9L6cjWU7GtG9jV2v1w/jG8DBEp2UmCkEAE\nMXQUqlHldpZhfr4t9b3jI0ttXeSrOZZe5T/d7f2Hv96eSvp+6b63OAQ4Oz0TE0CK2rrEHXxF\ng7Bsja3jO1d/f9+0d/vY5hDg1DyFM0pJ3sHXs8QdfJm2jlvvi2zdsH7fAdmvt5dBRa8fv5sd\nApyXMRyFJ/Xm2YNHrq2z3oPhKMvWobIBHIWnjqI38M2XesdtXbBTVGLrhlcICRyEaS5ivrQm\nyJQUsa0bgJDAQZiGYAM1QcW2jptvBLZuAEICB6AXVap9qCZo0R18lLZuAEICu5KqV4jSy6Gj\nSUPuLeVL5tuitHUDEBLYkXQFXZQpGA01QX4pQwO2bgBCAjtiJ51zEWOwvZqN2L4TttjWMV9H\namGJrRuAkMB+eBdwHqqHFJ4df97WRQaPAjry2+psGl4PIYE9WSYkoSNRzOCXem8xMcOcrfOU\nBSGB/VghJHXLhNE/qmjrvG8gJLAjS/pIbqW3LvVuxdYNQEhgR5Zk7fTDWtzBo1Zs3QCEBHal\nYBypl6+x6R9bsXUDEBKoR0pVfS8n9vamfzQGj5qwdQMQEqhF0uf1cl7v+Ye1VLd1AxASqEUy\n89D31qzesZqgJmzdAIQEKpHMhU8BSc/qHasJasLWDUBIoBJpIVnzbblKaszWDUBIoBIpIU23\n8LmTmzi2TicbKtu6AQgJ1CLRR1J38HmzBBm2zghHdW1dxuoQEBIgIZ61M2b1Tsy35WQZ6tm6\nYTUiEqhHbBzJmNU7Md8W97MMdWzdsBpCAs1RPt/WtrZurhhjWA0hgeboVdI79w6+bW1durXT\nLiEk0BoL7+AjsXWLswwQEmiM2UFYx9Z5Ylpj67zWzGYZxAKEBBSxvv+uBGqCCqZRlWrjho7k\nnrewdWoJQgKCZA3pbgTDUa6tM05hD1tn/gBCAgLbD9XBt3WsxNYZShM73NTWmR8gJDDh/Qte\ngVipd9LWaTEx9XWWrVs/eGQAIYGJ+kKKZhnStk63WvbxpEJC4cjQkXnsVbZuAEICE9WFFM8y\nlE6j6mQbNrd1AxASEMz3kTbN6iXu4MuxdVw5Ok9H9LYu8HeAkIBgLmu3RVZPPPFoZvAobesM\nP8dEAONB9cRtnXtSpeGIQ0jAIB1xNsjqqWeHhWydnWWI2zqzdbqbpFrt58Htkwj8+1AcjjiE\nBHLZog8lhRS2dXIMKc/W6a/MxjohzN7K/KX3TYzIWggJ5CEvVUIh9dOs3ptPzMB4wta53y6w\ndQMQEshDX8tUUhKPVk7UBMk4lGnruP1uqCSerfNiWLLJ8XOHkEAmzPh/EsRMqv4EJ25N0DJb\nF1BJwNY5mywLRxxCAtno/gaNksaA9PfX9+6TLGUxg5n7XmbrPDsasHXmz9PBNrkWQgLZeCmx\nVYw6Ys4Tj7St49YskAttnchUcHet+UvGbU3FSK+FkEA2pIk7qSPW997s+Fwlv9VVXmzr5AoW\nqBiyTyLr34e5niGEBPKhCEhiCNbI1vXB2fGZMYa0MFs3fhUaRwrpYsbWzZ0WhATyIShu6FWO\nYX52/HW2zvoqbOuS3+SvHYCQQAmrk9/yYS01J2Yo1VHOOUNIYD/Ek1p4+cQMi2ydPvC2tm4A\nQgK70Ysx2AUTM9SzdZkhGEICu6H6R6smZuD6jcTWEYQjDiGB/ejlIGz4YS1RW+eJidbWpZqc\n3yOEkMBeiKdZ/kUe1hK1dZzbizvaupzTmoCQwE5Mj7LsxVMmvCzDgW3dAIQEdmIKR5OGAg9r\nObCtG4CQwD5Mz7IUGgo8rCVu6wLhKN/WRSobIt/krvSBkMDWiFRdaHJ8WVKeOd9Wsa1LVaQm\nI07xuDOEBDZGjcLK0jonzWBkGchtndq+1NYVneEAhAQ2RgpJDh65tk4FDnJbF1AcD2zpryvX\nEYQENmZSUaxE1bxXImbrmCWgkmxdtHuUtnUlpyeBkMAEW/LvcAZ2OPKzdSoORW2dbB83PuVl\n66T03DMjtnUDEBIYILhBIozINIRktGBihtJBWGEZ7SbRZhkEEBIYsP+5J6FXr39eURA3boPV\n4YjU1slF5mpjg3DEISQwEu1MLMe4ZcJN1elB2HCWIWLrlEUrGoR1Q8wW4YhDSGBkQyFNdx65\nMwV5pd6zto6pgJVp64LnQ59lEEBIgG8hJJGsCyS9ebDUO27rrDeWb+tCqtnG1g1ASGCAvI8k\nhKSyddboUaLU27d11luRrfMatZGtG4CQwABp1q7Xs6i6o0fcsXVaODPZOrd71JKtG4CQwATZ\nOFLfy1th/Tv4HFuncnfztk6GTCdAEdm69ecOIQFaejk1Q+gOPsZMW8dKbJ2Rl6C3ddlnFwVC\nArSMjzzqxTxByTv4pJ1kXiSy5KKikg5ctLaOJBRDSICUKSD1kVJvMQjLdbaOz9o6uWc1JERt\n64rOLwaEBEgRDw+TAkrdwWdkvmtm63LPLA2EBCjpp7mCAnfw2bbOGDwK9ZEcWxeMQu3YugEI\nCVCi5tuK1ATZg7BKWJwnbV0oHDVk6wYgJECInm8rXBPkDcKGUt++rVP7bzQccQgJkCIf1tJH\naoJcWzf8Zltbl2gsoYwgJECJrgnq4zVBOls3/MQT0162jjIccQgJUDL7sBbL1nGugo9eTNg6\nxtyv1ti63HPKBEICZMSzdWZN0ApbpzN8zlrvQ/CLzHWLgJAAGfMPa1lu6+SKFm3dAIQEiPCf\nHWbbOrZmYgbf6TVk6xbuE0ICAXwZKVsn9ROydVy/xm2d3iKQCfc+BL8wVm2gI8YQkQAJYVvH\nmKxeSN3Bl5Gtkxv6a90P1k4CbCIjWDuwhl4uhB8JWz4xQ9DWyc8ssNb7YPwkxAbhaNolhAQW\nIx9lGcsy6GydTrpFSr1Ttk6sYe7G7qbRL7JWLUXsEkICi5FCCumIBWb1Fr2lUE1QytaJX6qM\nQ9rW7ZtlUEeDkMBS5EMmAraOmdM/qqyDjCvDjwtsnfqhs9b7EPzCWEWuI2OPEBJYyhiO/CeU\nG5M/SvlwHZ7822FnbZ1RIc7X6Kj4BGewhAkhgYVEZ/WWMYlLPyfTdpwvsHXyh7b2uPch+EXW\nqoXYO4SQwDKm+U3+/PlNuIpBk60zQhPn1iL337xsnR23Gg1HHEICS5GzertK0sZOJaytbB2z\nQ8tcts5cxZbqiDwc+TtcK6TPF85/X7qXf8sbNXcI0CJKR7aSEqXefh9p3tYpLWqH520a+8JY\nteD8UoQOtVJI393z0617QqokCKl5ep2tM5/AJ5Uks3Xc0JLoMnGeb+u8H3J30+gXWasWEdzf\nSiHduy/+073wr+5esoe5o0JIzSKGYIN3HgVtHTe6R1JGlp/LsHVyn9zbIvqFXkOso8j+Vgpp\nkMRP956hDWsPENJRUWOwiVJvq5iBW+m2VDgK2rrAYlVbFz0SgZBeu+8sIXU2uYcA7aBn9Q6o\naMrHFU7MMGfr/MWa4Sixu9XW7ue7u/Esa/fvBiEdC8dP6Yfwld7B5w/C5ts67l6/JTqaOcEy\nkqpcn2zouo8h2nzP//Dx2t1/xz2EjpqtMrAPTHdvBEJIc3fwndTWpfe2Ov19G3pI/OUr66df\nXTdsiD7SEbAz01wZu6yHteTYukA4atfWze1t5wHZ33v3+oCQjoBzlfdiWu/gw1pSts4Tkzt4\ntMbW7RWO5kW5e2XDR3f7hpAOgC2kfnp+2BSO/EdMcHPwyNaR3skGti7eekod5exstZC+X8fM\n3W/2z39e5vtAEFJ9TCHJZ7VwUVrnTP5o2DreiK3LOsVMsna2Vkj3KTXQ3fKVxN8gpCOgY4l8\nxASXzw7zJvZWfi6eZdjR1u0ejvhqIX1298cgi8/urXxHeYcAldBZu+khfLkPa3FsXXLwaBtb\nl3uGGWRrcqWQbp1IHdBmrCGkJhDjSOLZYTJb58yOz2ZsHdevh7N1BaGNoLIBQjo5k5BUti74\nXFht6xgzC+pobF3B4BGlrSvZ1UohvYiINBSuEgIhNcI4ADu9uL4uYutE5oHW1rnNai0ccao+\n0vet+yzfUd4hQD2mLAOPZut8W7dgYgZCHRGGo9Jdrc3avYqanpK7KAoPAWoxRB/ei9q6ULbO\ntXWceZFIymgXW1d2egnKFUkyjtS95lUILTwEqIQYPirI1kkF2QmHfWwdZTgq/wnmbAAxZLaO\nlWfrcidmKNZRm+GIQ0iXhLmXivfFiEgyME9IYVsnFeT0kTxbFwpH2bYufk7RNYUsDGyr09+K\n+/uiBsweAhCjB1pjX0xM8z/mZuvkfzLnwKO2Th/XW6hv6xbviE5I3XCDHxEQ0oY43f7AF0JC\nfMwy9G6t94yts8eRDmbrlu9orbV7uw139H3fun/8tSOLSRDSdvhXtPuFvhF2ytY5nq5sYoYd\nbF39cMRXC+m9+xnff7o7f9ANykJI25EpJJ1l6J0sQ9TWeWLaydbFVhSyTo8EJUJqga5MCELa\njlkhTeEoON1Wbm3djraOKhyt3c/qolUZkW4Q0kGY6yN5FUHG4FF7ti56mmWs3s9qayf7SO+F\nk0RStwrkMpO168ekm/+QCR6bb6vU1pHWBKVONB+CsEZyY99UItTR1dtBSJsSHUcag9Gw6E9w\nIhINKhzZCe8SW8f9pbq2jmQ3RCVCQ1gapuUiAkKqwlgSxKZsXbisLn++rcPYOiI1orIBCHo1\n31a4jkFnGZqYRpUqHJHsBUICkl4au8DDWoyYZNu68ZeerbPfNrJ1BacWh67QlUpI/17XtmT2\nEGBT3GlUe1dCsfm2fFtnvW1i62gEQCej9UJ632SWYQhpZ6Sti90yMTd4tI2t2zYcUcqIIP0t\nyZj7e9khwPaIWet6P8vg2bp4OMq1dSxbR7HmUiiAVkYEA7Jf/N79/t7xxL6j4OW+uT2rd/Bu\niVRNUKmtM5vCIx+cDZ0TSJ5fHsQyIikR+nhGox/ae80hpK3wRmMHdLbOUdJCW7dtTdDsOc5D\nHY44iZC+h4FY9JGOgRFDNKOQ/rInZpi1dfbBeFBQC20dhQQ2kNFqIb0+rd1v98L/QUiHwAsV\nA6ksw4ytS2YZtrB16dPLYRMZUTxoTJQJYcpiOkLdGKI9O+8D1sNavELvtK0zd7a9rSMJR+t3\nEWRt+vtj+PTW0d3T5x/inMTFEuzGUB3VeR+I3TIxiYjN1gTtZ+tmTi6Dzf6FQmVDHVJiCXZj\nyA7s7Vw/EtbPMpzM1m0no1VCur1+/nuQNsY7xGlJiCXYjaE7sKfg0HxbzMzWqYrvgK1jloDa\ntnVbymiVkKYJT14/vgsejVR4iLOSEstyIeV1rcRWvfgYtHVTIZAUUGmpd5u2blMZrRLS49/n\n28tU1vD6/vVTuVWHYgshlXWteqkk/wnlbrZOC2dzW3fUcMQJ+kg/n2931NrN4QSLpFgW9pFK\nfjaGn3HJ1xFzBmGVQHewdclTW8HmMqJKNnzfIaQEfrBIXfXLsnYlgUzfeeTpSOUYJvGfwtbt\nICMSIf37eEFESuLLJi2WnM6Ou02RkJxZvT1bx85k6/aQ0Woh/X6Oz3V5+SCtWT2ZkILXeF5m\nILrLSIzLuWx6dc9EwNbZWYY9s3WR1h4iHPF1yYbv99uQt3v7Jk+CX0BIFLsMXJkZOpqeMsH/\ngnMR6yxDE7ZurQz2ktHq9PfrJ33um0NIS/aY27USeYbem9WbxyZmqGjrZs8myX4yWimk943G\nY08mpKV5uLkdurvMcovTPXx/f76QjJjEjbd6tm51OFr38zIQkXZgWR4uuUMrNJQwPs1S3DGR\nnNV7M1uXW94FbFgAACAASURBVBN0oHDEV/aR3tBHymRdaiG4N1kAl00/TFk39pDknbC9PcGJ\nMXgUzDIQ2TrvZCInWXRy3q/3lRGydodEJadLfjOWAEWeCevaOjscHc/W7S4jjCMdEXFJl10t\nUzGDsHX2nbC+rZNHobB1ySzDNrZufxmhsuGIhFMNM78ZZbTNfFtpWxdqufocPoXDhSOOWrsj\noq/g/Etm0lFg+seArRPBzgtHu9m67LMK/biKjFD9fUj0dZ192UxFQTkTM8julxGDtrF1W4Sj\nWjLC/UiHxNBQ5oUzZuyG19CM3ratY3vZusi55Z3QBj9eB+6QPSTMd1ZJhgA0DML6Y7Butk7d\nPXE8W1cvHHHM2XBYMoU03XXUi4kZ7JAUzNbJ7hEzFRS2daFwVM/WVZURhHRY8oQkboSVt0w4\nISmVrTMVFLF1XlsqhqPKMoKQjovT6Q/Tlz9amUWydfbbBrZuVTha/lMiIKSjkpO166f/yqZR\ntV8KsgwrbV3yTJJUD0ccQjoy8+NIk5CiNUERW8eiWYbtbN0KLbQgIwjp1IhMQzRR59k6L8Hg\nZBm2s3Ulp2X/sgkZQUinZhTSX/iuIxGTHFvHuS+mlm1dKzKCkE5KL18ynsFnZ+s8Mam0RoNZ\nhmZkBCGdk16mGVSpd2TwqMDWyX0H1JOpo/OGIw4hnZN+elZLfHZ8Pfq6o60LNnWxGpqSEYR0\nSkSyLqPUu76tKzw1//iNACGdkNQgrLZ1bLWtk4HNXet9iF34y8PRwt9tB4R0PpJZb3fwaLmt\nM4ac7LXuBx678E8TjjiEdEambF1wHlV/8Gi5rWNajJwnwxGtjlqUEYR0Knr1Kp9QHsvW+VmG\n4mwdC9W1cv8Dsa1rU0YQ0pno1fBRJOvt2jqmdcQ5L7F1xjeM8x1tXaMygpBORC9Hj/78rHfY\n1nGKbB3jyXAUvvRPFo44hHQepme1cDHfliUkclundzKno4itW3SGDcsIQjoNvRDSfvNtiWTF\nXrauaRlBSGfBj0c6HEVt3fjLFbbOGUfaNMvQtowgpLPQy4cehXREb+u496X34UrhiENIp2B6\nVMsko2CptxQUma3z1rofeEQxixTRvowgpGZgy6+WqWM0uDpXSHaWQXZoguGoOFvHeTockWUZ\njiAjCKkRVJXBEqaAJONRuNQ7aOs4txdjts4Q+SpbtywcLfhNBSCkJnAu4CLGZ0z00x18nreb\ngtHUTTJ7SvJoGbaOqYi21tYtOLljhCMOIbWB56RKiDzzSNk6EZi0hAptHTN+5Dc51G4yW3cY\nGUFIbbBQSKIiaOge9WII1hqEVcsRW2eHo0i2Tstsb1t3IBlBSG2wTEjyPtjUM/hkOFGdMC/f\nnZuts2S0h607kowgpEZY1EcSQpq9g4/LwSQVl0qydQGRF9u6s4cjDiE1woKsXd9zPTs+m5+Y\nQQzM8hJb52bKnbXeB7JwdDQZQUjNUDiOJEu9CSdmiNk6Gcy8te4HHpZMuSiOJyMI6aCIoaO8\nWb1lWPFs3bQvneUWn+VB9DfGlV1u64pP7oAygpAOiqxR/fNKghxbZ2gmZuust+BCxNZtk2U4\nYjjiENIxEcV1GdOoRrN1ebbOWnBtnd0mmizDQWUEIR0TXaMa6BtZtm7YepWty87Whdp5GRlB\nSIfEKvWOZ+tq27pSWRxYRhDSEXFrgpSMltk6I6cgj0Bj6wpP68gygpCOiHq0sh+Oltg6uVta\nW3epcMQhpAOiZ/U2n1HOjZikbJ1TbpqydaFwtMbWlZ3T0WUEIR2QcKk3E+HIsnU839ap3Vew\ndceXEYR0HHr5HuwecXm/xAFt3QlktK+QHm9dd/8WO0nuBULy6KWS/rx5gpSWKOfb2s3WnSEc\n8V2F9Lh1A6/TTq4gJJa6SJIrPaZp67gqUQ2MvwZs3XQg/bKHrSsTxklktKuQ3rvPp5o+b/dx\nJ+cXUrKiu7TcuxchKaKjhmxd9inxE8loVyHdph/+3l5+ryEk47VspbGZuNTG0rq+LVtHEI5K\nNm6cHYUktfO430NC6kwWHqIlfOPkrmRz7k6HLfEwy1Jb54kpbuvC4Shl62KnlcuJwhHfVUgv\n3UMu3S8QkeaElGHuVNh6BqNh2CgwjSonmW8r2E/a1tadS0a7CumzexNLv90dQrLixMwengoa\nn8E3U+p9HFt3Nhntm/5+V+r5nnFvZxDSTB9JXt+JK0p6uiEUTTVByWlUD2TrTiejnQdkf17l\n0u/bBYSU8m66oDSxAzF6xISOrJKg8MQMjq2zw9GWtu7a4YijsmFTUskEmWlIe7tBSExPzNAv\nmJhhH1uXOAt30zPKCELaHeb4rOS2Y7KOrZuYYQdbV6CNk8oIQtoZw+7lDMmON+7lPoNvna2L\nhaMcHSXPYemmBwNC2hXbbM398zwWBfnPDvNsne5wObbODkdJW2e0cTNbd9pwxCGkffGMVZJJ\nR9Gnk2fXBBXaulCD4+3OF8eZZQQh7UuukHr5+ifuhA0PwvK4rVOfW7F155YRhLQvmUISs6iG\nZyOet3XiK8fWST1R2rqCcJS74VGBkHZlJlUnrnYpJH824gJb52frrGErEluXOFPnvHK3PCwQ\n0q5kDdKKCtXIREHZ2TrjUMrdLbZ1a7IMF5ARhLQ7iVSdvNrtcORkGuazdfKzjj/MeV9g60In\nMnuuRdsdHAipGYYLTvSNOI/dMsHmbR1TilThyAlIe9q6a8gIQmoIpp/BN5Z6RwaPcgdhVXbB\n2MLN3y22dQhHLhBSM7BpYoaeqekfg4NHSVs3fmd1xFQwCujIOHahrcs8o8vICEJqiJ6NEzPI\nrLd/y0SGrbPF5PSOnEzHYluXqY8ryQhCaodeTMzA/vzRo6Ct04Ounq1z3qaIxOxvltu6rNO5\nlowgpGYYbd1YoiqV01vhyLd1PBSJpp05OqLN1mWdzsVkBCE1g3iYZbrUO2HrjNFWt4YhEIU2\ntnVXC0ccQmqFfhLSH1s+MUMkHIUWtrV1F5QRhNQCqpJB3H8UztYx19aNvzV1ZCbldDgitHVZ\nCrmkjCCkesgLX06iOj6Dry+cmMHQ0fSfpSZrodDWLQ5HORudEAipDoYZm4TU99LWlU7MoNZz\nXbqgDuMvRHSEcLQOCKkOOnCIqqA1EzM4o0aLbF2oeelvAid1WRlBSJUwEmrRWm/H1hlpbh2O\nrKx3ZVt3ZRlBSJWYrrmei1LvxY9Wlo7OyNrVsnWXlhGEVAnZORJCUiWqgUHY2Zog8VrX1l07\nHHEIqRZsGjqKBSRuxKSYrbNrgrazdRkaubyMIKRaiGKgWIkqE+EoZeu4ft00W5dxLrObnB8I\nqRJ9Lx8e5pg6lWXItnVKVBvYupxwNLvFFYCQKtH3cvTIFRJnTClpM1vHrD142wd+FwHhaAJC\nqkPfOylvq3eUZesitXWhcOTqSNYUrbZ1kJEEQqrD6Ot4H55GNc/WiT25tk4dImHrmPt94GOG\nSiAjDYS0F9a//2MH6W/KN3hphq1tndJRIunAA8Iq3uBKQEj7oIdMRcJ7mtW793tHW9s68c5Y\nSGJOg2fPCEggpH1QgaMXhXXTrN5WSNrH1nE70jmrAr8Lnw9kZAMh7YK+1Hu3sK7vncz3DoOw\nTsquOBxBRh4Q0i4oIXmVDL0ZjnawddP3zFxbHo5m1l8RCGkXZoVEZ+vC4cjSETdTH8W2DuEo\nBIS0D2PvaNJR6hl8+fNtRW0d95dsadg6KLV1kFEYCGkf2DD7I5t9Bp9l6/g2tm56Zf6q8BfO\naaRP87pASItxCwOS9GIi1Zln8O1i69zuF49s6wMZRYGQFmL31+coeQYfty/zDWydfEU4ogNC\nWohzMc+QegafkWqwbZ0+Aqmtc+KatzoCZJQEQlpGNHUcpPeHj6yaICazDPO2Tophqa1T754z\nhYzWACEto1BIvQhHXpZB1AR52bqorXMOmx+OnE3ZbPoueyXgENJSioQkjF0sy2AMHlm2zhNT\nlq0LNNLeVPpHW2WpE0E4mgdCWkhJH2ksUZ1Gj3onHjGpJBUitI70EehsnV4TX+0AGeUAIS2k\nIGvXiyxDqtRbhaMFtk73duZtnfhkd5BSpwEZ5QEhzRMZMMoYRxLPVY7W1uXVBFkC8kTDHO1x\nc2t709An7qrMOce5UwQTENIcZQNGJn086a3v4JO7D9i6zMEj+5fW1s6mgU+hL7JWARsIaY6S\nzpCFmrZO6si888iwdbKkwbN15oFns3VWfEzaOuf04meGcFQAhDRDUXrOohfzMqwq9U7bulDi\nwV1eausgoyIgpBkWC0kMwv4xD2XrzMGjpbbO3sRZXB6OIKNCIKQZlgtJPjwsJCOd7i6/g88f\nWbV8XYmtizYeMioFQppjYR9p0NEoI/L5thxBqX2ZX7qb+p8QjmiBkOYoz9qNSe9pniBfSJS2\nTi0FvuSu+LPDEWS0BAhpnowBIxMxTxCXg7DurN6Ets7c0vrS+5CtI8hoGRASOUJIf/IpE72t\no5itE0NKRbaO24IgsHWQ0UIgJGr86U2MWb2jtk6tKbV1+sAEtg7haDEQEjWTkEKDsClbp+5H\n4pzc1uWGI8hoBRASMSLTkHoGXyBbxwwJLbJ1djhK9JWi4QgyWgWEREt0niBn8CiSrWPqOi+0\ndWYbkrYuGo7mzw0kgJBoGUePWPoOvkC2zolMa2zdTJYh3GyEo7VASKTI+bb8O/iMLENoENbs\nI8nvl9k6qz15tg4yWg+EREpsniDf1unukBGTqtg6yIgCCImCXr5n3cGnfRznWkxyHGlvWwcZ\nkQAhEdBLJS0o9bbExPU7ua1DONoUCIkAKSRfR3ml3qaONrN1oYZDRmRASHFY5nUmZyMO6ci1\ndeN+uQo7+iVs60LhaImtC54JZEQIhBTD6v4nsSdmILV1ujXeQpGtCzUbMqIEQorhWKw4Yqqg\nhmxdTpYB4YgWCCmCHxRijGOwoTthk7bOExOhrXPOJHAKkBE1EFKEDCH16tWvZciamKGWrYOM\n6IGQIswLSeTq5EOPeldGjDVg64LhKHpKYDEQUozZPpKa/jH08DCu5FTX1gXOCzraAggpxlzW\nTiS9I9M/Tm8JWxcMR14YDIejXFsXaD1ktBEQUpz0ONIkpOT0j1wt+7aO61c3HIU6RYtsXeCM\nEicE1gAhLURkGiLFDE1k69wmQ0YbAiEtZLqDLz79Y+yWCc5r2TrIaEsgpHJE72jMMiSnf5QV\n3SFbFwlHm9k6hKNtgZCKmabbCiW9rYkZZM4uYuvEztxwRGPrPNVARlsDIRUjHtcSSnqzQDjK\nsXVy11S2zmkxZLQ9EFIpYrqt0J2w2tbJftJaWxcLR0W2DjLaAQipAFnJELkTVtk6PpOtE7vL\nsHXGwZfaOoSjXYCQ8pGdo2S2jpdm6+Tel9q6dDiCjHYCQspHPssyla1jcvBoN1tnNTGj0g5s\nAoSUzeTpEtk6KSDG3XC0m63jqY9gQyCkbOQd5eFsnbR1+vPutg7hqCIQUi5iUu90tm5KeksJ\n7Wzr7E+Q0a5ASBlkZOvipd5VbB1ktDcQ0jxGtm5+dvxglmFvWwcZ7Q6EZMKC/5KrmqC52fG9\npHcdW4dwVAEISaOCiY2VrUvOjl/L1rH4dmAndhXSv4/XbuD1/d9Wh1iD47kkooeUMTt+LVsX\nXQN2Y0chPV46zX2TQ6yCOe9c3zCROTu+/LUXjra0dYlSVrAfOwrpvbt9/YxLv9+37n2LQ6zC\nF5KZZXCLgmRZXdTW6X0R2zqEoxbZUUi37kct/3S3LQ6xCk9IQ/RJ1gTVsXXBRrs/ATuzo5C6\nLvZBfGOw8BDrcPtIY22dyjJ489a1Zesgo7ogImmcrN2gISNbF8oyZM+3tbmtg4wqs28f6ft3\nXGqzj8S5PY7k1ASZ89YpESVsHdev29g6hKOW2DP9fTe828tjk0MQ0ns1Qb1p61Spt3k33/A7\nUluXGjyKbASqsO840vs4jnR7/WhyHMlGZuvmaoKMeMQ5ta0zGxTxeJBRE6CyIUIvaoL6VE0Q\n17dMMNfWqXC0ia0Lfw2qcVUhsbl/yHW2LlETxIzQpDMVgXC0la1DOGqFawopUlU3IuYizqoJ\nUmNJ+svw4NE2tg4yaoeLCsl4dehLaoLMbJ0xlpSRZSCwdZBRS1xSSI7tMpDzm8xk6zjXAUiG\npJ1sndP1Ao0AIZn0zsQMqWydEo7ZO9re1sWaDuoCIZn00+DRX0hHrq0bd6AydvvYOoSjZrmk\nkGJ9pF7cCetN/7jNHXw5ti4kMMioQS4qpHDWbhRSsNQ7Zutig7DRwaNSW+d/gIya5JpCUrkB\ni9HXjeEoa2IGxr1ItNLWMbnvwEaeNkFTXFVIASYdTSoyp3/cy9bJQ/HQRghHjQMhKYYxWDlP\nUHJiBjWCxLkvpuW2jmmV+htZAQ+0B4SkaxkSD2sRAnIHYQltnRXqnI04ZNQ8EJKqZZi5g0+q\niXM1CEuarWNWJPMtHmTUNEcWUihhsIBALYORZRAJAJlxYH44yrJ1wXAUXPZUg3B0BI4rpEgK\nu5hoLYM9MYMUlApJZbaO+0tuMsHQoKsjyKh5Diwk43UNk5ACOjJsHVddGKaWyQdhnd6XswK0\nzWGFFOqTLyHyKEs9eKRqUsUbV1m8VbYuXLLAOPfD0cozBHtwdSFFHjJhJBdUPR3nOtvg9ZGK\nbV3oXMzdyWWEo2MAIU3do9CjLI0byrn6qEtUKW1deAVkdBwOK6S1faRevv758z9qW6fiEOdy\nIIlbkYjM1vkfGGR0IA4spFXprMTgUWhiBm50VxbZulg4Stk6yOhAHFdIqqeyiNTgkVMT5Opo\nOrRaVO9JW2c22zqF2AqEo2NxZCGtQA4e5dzBZ2elDQVtZ+vYmn8jQA0uLaToIKwzeLS/rYOM\njsY1hRQfPHJsHXN0NL5sbOsgoyNySSHFB49cWzdsGLJ1KhxtYuvWnBuow2WFFBg8KrB1Ykfb\nZOvAAbmikGKDR03YujUnBupxTSFFB49q27p1JwbqcXghlXcpkoNHNW0dwtGBObiQCsoberkQ\nHDyS4SjX1slwRGfrIKNDc3QhGa9peqmk+OBRia2z3hbbOmNr6OjYHFtIjs1KIYUUDkeVsnV6\nCTI6OFcRUqomKGTr2J62DuHo+FxLSPm2ju9n6xh0dAKOLaSsPpIIRsNrWEZLbF0qChXZOild\ncGyOLqT5rJ18dthg6+I1QcrWTd/71amOrbMa4Hxlts9tq7OC4XaJk3BwIc33L4bBokhNkFnq\nLWydLhQa961ftrF1jCMcnYTDC2kOEZByaoLUMqO3dQhHJ+fsQhofeRSpCZK2jqlsnQxNoT7S\nOlvHvQ8M4ehMnF5IPFwTFJ6YQSrInv5xC1untQpOwcmFJJ4d5qXrQraO62ydGXs2sXXzORJw\nLM4tJPXsMC/LYFs5lf1jRhpwC1uHcHRWTiukabYtMTu+LSTu2DozHBk64pyT2zq1FuHoZJxV\nSGOqTj6Dz1aSZ+u0cFSmwQ5Hrq0LhaNcW2fvGJyFkwppHD16vqpsnfPwMMvWDT9wB2G5fvVt\nHfeWsm2dFefAeTirkIb/Cz3LUqadc2qC6G0dwtFpOaeQRDiiuIOP0Na5kQ2ciJMKafgvNm+d\njEhJW6fCEZmtY44iwak4pZDk/I9u0rvA1ok90dk6d0/gXJxRSL2oretjDz3a39YhHJ2dUwpJ\n3FDe+zWqNLYuGI4Sts5JW4ATckIh6RpVM+lNaeu4v5SqrWOBLcDJOKWQIvPW1bB1jLHAFuB0\nnE9IcR3VsXWBb8H5OJ2QYoNHlWxdaAtwQs4mpOA8QbyWrQttAU7JyYREM40qka3zfwZOy6mE\ntJetC4cdbukI4ehanElI2bbOqMPmnJfbOuOYEVvHoKOLcXwh6adM+DISto4VTKNKYetMGUFH\n1+DwQlKT4ydsnbJy+9g6Fg5e4MycRUjRLIM935afZTB15FbELbJ1jAV+D87OoYXU66dMzGYZ\nmBOOQrbOeltm6xjC0TU5spDkrN5Zd/DtYusYs5fBZTi4kOSdR2EdST/H9rF1evfOBuD8HFhI\nvTFvnZ+tE52V2XBEZ+sYj4gOXIDDC+kv8JQJ19aFaoJCti4VheZsnR2OoKOrcVwhTZkGIaHy\nO/hCtk7uutzW6SFedwNwDY4pJJGr6//+PCFNlo7JLPQuts4pjYCOLsghhdTLiVSjd/Dp7lHI\n1iVrgsptnaVKhKNrcmQh/bmmLmjrjPSD30cisHUIR4AfU0jiCeVyYgYvW2eWeuuCO05g6/xw\n5OgQMroqBxNSL1/UrN4sNKu3tnVM6ElV7qyydWZTmF4JWweOJaReVQVlTswgl2VgcsLRSlvn\nhDPYugtzPCGNY7CBQVjm2Tqur3cmTR/nng9bauuYq6Pi0wTn4VBCmmJRJMvAPVvHVbaO2NaZ\ng6/uO7gmRxSSn2WID8IypoMH474AFtm6QK8IMro4RxKSeJqlenjYjK0zvJ38yNwLf6Gtc3tF\nCEeX52hCCmcZZm2d1UfiK20dwhHwOJCQ5J2wgVLv0MQMjJuCml4obJ2/MXQEDiakPz/NoCKO\nDkGcq2Wje8TcQLLI1jFvY9g6wI8kJOOGcvcR5TIiqYChR3mIs3XM+x4yAgNHEpJfzMCdmqAl\n820V2DpfRwhHYOIwQgo+olyFo7St414AWWbrvEXICAgOI6SsWb1lKDIVlLJ1oXAUz9a53yMc\nAcWeQnq8dd39W+wkuRdvZWieIMfWueEox9ZxbymerfOWISOg2VFIj1s38DrtpEhIyekfQ4NH\nxLYutAwdAYMdhfTefT7V9Hm7jzspEdLih7XQ2DoGWwfm2FFIt+mHv7eX3yIhLbB1nphW2DpL\nRrB1IMyOQpLaedzvJUJaYus4txdX2DoW+B7hCLjsKKSX7iGX7vlCIrF1rgCybR0L2rq5MwXX\nY0chfXZvYum3u2cKKf6wFpZp6+R2nOeGI2s1sgwgiz3T3+9KPd9dnpCW2jr1nS4YCuqI+0t2\n2EKWAeSx64Dsz6tc+n3z9tKZTF+terSyileqfGgi09bpXVrfQ0YgSMOVDRFbZ94yMZet0xuU\n2jpmf48sA0jTrpDmHq2sKr5jto5rASUrU0PdIFdH3pYAWNQQUrp/JA4x92hllrJ17uBRma0z\nO1bG9whHIE6jQko9WpnpAVkesXViN8yWVZ6ts8sh9MaQEUjQppBmsnXi87Bl3NapTAPz8tbh\ncGTqBbYOlNGmkIK2jnMloCJb5+Wt58MRbB0o5BhC4ipjN5+tk7uw3vJsnW/iYOtAHocQkp1l\nyMjW+YEkw9YFKhcQjkAmbaa/Q7bOGFnNtnV6f2rJW7B2YKkG4Qhk076QVJZB1ta54YjE1iEc\ngVU0LyS/JmgDW+cnyAPpOwDiNC4ky9bpSfBpbV1gvBa2DpTRtpDMwSOuB4VobV3oe9g6UEbT\nQgrYOhWTOPfFxH0FzNq6oL4QjkAhDQvJzdYZ6iGzdSykL39UFoAZ2hWSbes2ydaZ4uHOImQE\nSmhWSL6tm1bYClpl6zzx6EWEI1BGo0JK2Tq5xUpb56fo9CJkBAppU0jdDrbOWHYXoSNQSptC\n2t7WmQdzFmHrQDmNCill6zwxLcjWGcuwdYCANoXUpWwd5/Ziqa1DOAL0tCmkDW0di1g8hCOw\nhjaF1G1l6xjCEdiENoWUtHWhcJRp66wpiH0dEbQcXJQ2hdSlbB3Xr2W2jsHWga1oU0hb2Do2\nE46gI7CCNoXUkds6efeFtylsHaCgUSGlbJ0bjkytiCXufsXskOPaOoQjsJI2hZS0dXob86OV\nUnAWEI7A1rQppI7S1ukbAd1NxRGgI7CaNoVEaOs86+bZOtqmg2vSppCGZ40x8d8zOnVqUTyH\njFlvxoL7FWP2amvZWwfAYhZc5fTC2YP6za7eAjSggRZoWmpLAfWbXb0FaEADLdC01JYC6je7\negvQgAZaoGmpLQXUb3b1FqABDbRA01JbCqjf7OotQAMaaIGmpbYUUL/Z1VuABjTQAk1LbSmg\nfrOrtwANaKAFmpbaUkD9ZldvARrQQAs0LbWlgPrNrt4CNKCBFmhaaksB9ZtdvQVoQAMt0LTU\nlgLqN7t6C9CABlqgaaktBdRvdvUWoAENtEDTUlsAOCwQEgAEQEgAEAAhAUAAhAQAARASAARA\nSAAQACEBQACEBAABEBIABEBIABAAIQFAAIQEAAEQEgAEQEgAEAAhAUDAgYX0r27bP1+62/uj\n0sHfbxUPPlD17CWVrwCLhppSyONWte3v41MLbnWupft48Jcqxx6pevaSyleATUNNKeR1ycM3\nyPjp3p5X0Wf3VuPg/7rbD/+5df9qHHyg6tkr6l4BDg01pYyvRU+xIeN1OnidNrx333z4C3zU\nOPhA1bOXVL4CHBpqShG/3b2FP2OdNrx2v3wIC681Dm5Q9X+BRq4ASUNNKeLe/TbwZ3x09xqH\n7VoICNXOXtDGFaBoqCklfHRf1a8jPvQSvmscthEhVTr7iUauAEVDTSlgdDX1/4y/tzrmqg0h\n1Tr7kUauAE1DTSngZUi8Vv8zPm6VrE0TQqp29iNtXAEGDTUlA/HI6bfRU1T5M5oPvb7XGsi5\ntSCkamc/UPEKiNBQUzIQl/Gax7jTtODJ78v9d/fDT0xZu9+aWbuKZz9Q8QqI0ExDSmjhz/hd\nMWX1Mf57/N29V2tBzbMfaOEKsGmmIeXUHsWod/DqlQ1Vz17TjowgpKW8Vf0X8WU8dL2Lue7Z\nK2of36ShppRS9c9Y11o8xurvKoceacRY1T6+SUNNAeC4QEgAEAAhAUAAhAQAARASAARASAAQ\nACEBQACEBAABEBIABEBIABAAIQFAAIQEAAEQEgAEQEgAEAAhAUAAhAQAARASAARASAAQACEB\nQACEBAABEBIABEBIABAAIQFAAIQEAAEQEgAEQEgAEAAhAUAAhAQAARASAARASAAQACEBQACE\nBAABEBIABEBIjaOeSuc8ni7+8XvrJoEAEFLjFAvpBf+T1gB/9caJCSm2WVMPVr0Q+Ks3DoR0\nDPBXUgX+eQAAAYxJREFUbxxHSJ8v3e1Tf3y/de/j4vO/9+72IR44XqmtVwZ/88axhfQ66uQu\nP96HT2+TkMZVnxBSJfA3b5xOw/l3d3/wx737noT03d1++M9tEtJzzWf3AmtXCfzVG8cS0mv3\neH716F4nvbwOinrKaRTSP86lyQP7g79641jWztCUoRm9CCFVA3/1xoGQjgH+6o3jCMn6HkJq\nB/zVG8cS0tQpcj5+Q0gNgL9641hC+hrSdPxTJhusrJ3cput+Kzb3skBIjWOPI40DR93t1/po\nC+nlub5aa68LhNQ4fmVD9/arPr7fuvs/W0j/XiCkCkBIx2esdAB1gZAOTNd9cf547d5rNwRA\nSEfmY+ohwck1AIR0ZD7vXfeCeNQCEBIABEBIABAAIQFAAIQEAAEQEgAEQEgAEAAhAUAAhAQA\nARASAARASAAQACEBQACEBAABEBIABEBIABAAIQFAAIQEAAEQEgAEQEgAEAAhAUAAhAQAARAS\nAARASAAQACEBQACEBAABEBIABEBIABAAIQFAwH//uNeWu3Om+QAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"Explain Wages with Height\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "plot(x=X, y=y, cex = 1, col = \"grey\",\n",
    "     main = \"Explain Wages with Height\", xlab = \"Height\", ylab = \"Wages\")\n",
    "\n",
    "beta_hat <- gradient_descent(X, y, 0.1, 200)\n",
    "print(beta_hat)\n",
    "# 2.945535, 0.951942\n",
    "\n",
    "# Draw the regression line\n",
    "abline(a=beta_hat[[1]], b=beta_hat[[2]], col=\"red\")\n",
    "\n",
    "# To get y-hat:\n",
    "y_hat <- X_mat %*% beta_hat\n",
    "points(x=X, y=y_hat, pch = 2, col='yellow')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}