OLS_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Randomly generate 100 values for each $\\mathbf{y}$ and $\\mathbf{x}$.\n",
    "\n",
    "We will also generate a residual for each value to simulate some error in the linear relationship between $\\mathbf{x}$ and $\\mathbf{y}$. We will use predetermined values of $\\alpha = 2$ and $\\beta = 0.3$ for generating these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "# define alpha, beta\n",
    "alpha = 2\n",
    "beta = 0.3\n",
    "\n",
    "# generate x, y\n",
    "x = 2 * np.random.randn(100) + 3   # Array of 100 values with mean = 2, stddev = 3\n",
    "res = 0.5 * np.random.randn(100)       # Generate 100 residual terms\n",
    "y = alpha + beta * x + res                  # Actual values of Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's do OLS!\n",
    "\n",
    "The aim of OLS is to estimate values of $\\alpha$ and $\\beta$ to estimate response variables $\\hat y$ with the linear model:\n",
    "\n",
    "$$\\hat{y}_i = \\alpha + \\beta x_i$$\n",
    "\n",
    "For OLS, we need to estimate $\\alpha$ and $\\beta$ that minimize the squared errors between the observed reponse variables and the predicted response variables for each predictor $X$:\n",
    "\n",
    "$$SSE = \\sum_{i=1}^{N}\\left(y_i-\\hat y_i \\right)^2$$\n",
    "\n",
    "You can estimate the parameters as follows:\n",
    "\n",
    "$$\\beta = \\frac{\\sum_{i=1}^{N} \\left(x_i-\\bar{\\mathbf{x}} \\right) \\left(y_i-\\bar{\\mathbf{y}} \\right)}{\\sum _{i=1}^{N} \\left( x_i - \\bar{\\mathbf{x}} \\right)^2} = \\frac{Cov(\\mathbf{x},\\mathbf{y})}{Var(\\mathbf{x})}$$\n",
    "\n",
    "$$ \\alpha = \\bar{\\mathbf{y}} - \\beta \\bar{\\mathbf{x}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "betaHat = 0.32867, alphaHat = 1.95155\n"
     ]
    }
   ],
   "source": [
    "xbar = np.mean(x)\n",
    "ybar = np.mean(y)\n",
    "\n",
    "# estimate beta and alpha from data\n",
    "betaHat = np.sum((x-xbar)*(y-ybar))/np.sum((x-xbar)**2)\n",
    "alphaHat = ybar - betaHat*xbar\n",
    "\n",
    "print('betaHat = %.5f, alphaHat = %.5f' % (betaHat, alphaHat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How well do these match up with the actual values that we used to generate the data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the estimated parameters, let's predict values of $\\hat{\\mathbf{y}}$ given the values $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = alphaHat + betaHat*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot our values of $\\hat{\\mathbf{y}}$ against our training values of $\\mathbf{y}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debgcZZn38e+PcIAEkSAcBzghBAVBwIFgBDTODCAjW4C4jiuO45hBcVyG3dcR8B0FZBw3lLyIjCIqA4IRIxgZQcUFNAugGFAENDkJEpawJUiW+/2j6pA+vZ0+3V1Vvfw+13UuTldVVz3dOTx31f1sigjMzMw2K7oAZmbWGRwQzMwMcEAwM7OUA4KZmQEOCGZmlnJAMDMzwAHBciJpmqSQtHn6+npJ78jhumdLujzr63SirD67pK9I+o92n9eK54Bgz5J0v6S1kp6U9GdJ/y3pOVlcKyKOioivNlimw7MoQ5YkbZ1+j9eN4z3/KOmnWZbLrB4HBCt3bEQ8BzgAeBnwkfIDlPDfTn2vB/4CvFrSTkUXxqwR/p/aqoqIYeB6YF8AST+S9HFJPwPWAC+QtK2kL0taKWlY0n9ImpAeP0HSf0p6SNK9wDGl50/P988lr98taamkJyT9VtIBkr4GTAW+m95tn5Yee7Ckn0taLel2SYeUnGc3ST9Oz3MDsEOtz5heb1bJ683T8h4gaStJl0t6OL3OryT91Ti+wncAc4E7gLeWXXcXSddIWpWe/0JJL06Pf3n6WVfX+J5GPUVI+qykZZIel7RI0t80Urh6nz19fZWkByQ9JuknkvapcZ6Kp5o0Nbh7+vuW6d/Bn9KnzrmSJqb7dpA0P/1+H5F0s280iuUv36qStAtwNLCkZPPbgTnANsAfga8C64HdgenAq4GRyuvdwKx0+wySO+Za13oDcDZwAvBc4Djg4Yh4O/An0qeWiPikpCHge8B/AM8DTgGuljSYnu4bwCKSQPB/SSrmWr4JvLnk9RHAQxGxOH3ftsAuwPbAicDaOucq/TxTgUOAr6c/J5TsmwDMJ/n+pgFDwBURsTS9xi/Szzq5kWsBvwL2J/kuvgFcJWmrBt5X77NDcjOwB/B8YHH6OZpxPvCitIy7k3zej6b7TgaWA4PAXwEfBjyXToEcEKzcvPTu9KfAj4FPlOz7SkTcGRHrSSqgo4APRsRTEfEg8GngTemxbwQ+ExHLIuIR4Nw61/xn4JMR8atI3BMRf6xx7NuA6yLiuojYGBE3AAuBo9OK+GXAv0fEXyLiJ8B361z3G8Bxkialr9+SbgNYRxIIdo+IDRGxKCIer3OuUicAd0TEb0kq3n0kTU/3HQjsDJyafm9PR0TT7QYRcXlEPBwR6yPiU8CWwJ4NvLXeZyciLo2IJyLiLyTBej9J246nbJJEcmPwoYh4JCKeIPl7GvkbWQfsBOwaEesi4ubw5GqFckCwcrMjYnJE7BoR742I0rviZSW/7woMACvTR/7VwP8juaOEpNIrPb5WBQ/JXfgfGizfrsAbRq6ZXveVJBXLzsCjEfFUI9eNiHuApcCxacV4HJsqxa8BC4ArJK2Q9ElJAw2W8QTSO+qIWEESWEeeVHYB/pgG1ZZJOjlN/zyWfhfbUidNNqLeZ0/TfedJ+oOkx4H707eNed4yg8AkYFHJv9X30+0AFwD3AD+QdK+kM8Z5fmuzzYsugHWV0ru3ZSSNpjvUqNxWklR+I6bWOe8y4IUNXHPk2K9FxLvLD5S0K7CdpK1LgsLUKucoNZI62Qz4bVpREhHrgHOAcyRNA64D7ga+XOdcSHoFSarlTEknp5u3IXlKOCUt/1RJm1f53qqV8ymSSnXEjiXX+hvgdOBVwJ0RsVHSo4DqlbFE1c9O8rRwPHA4STDYFqh13lHlk7Rjyb6HSNJs+6RtUqOkTwwnAyenbRQ3SfpVRPywwfJbm/kJwZoSESuBHwCfkvRcSZtJeqGkv0sPuRJ4v6QpkrYD6t39XQKcIumlSuyeVu4AfwZeUHLs5SR3tUekd7JbSTpE0pQ0zbSQpBLfQtIrgWPH+ChXkLR9vIeSlImkQyW9JM35P06S3tjQwFfzDuAGYG+SvPn+JA3zk0hSbL8kCZbnKemaupWkmSWfdYqkLUrOdxvwWkmT0obad5Xs24akDWcVsLmkj5K0wTSq6mdPz/sX4OG03J+ofOuzbicJdvunbRdnj+yIiI3Al4BPS3o+gKQhSUekv89K/61F8h1voLHv2DLigGCtOAHYAvgtyR3kt0hSN5BUBAtIKozFwDW1ThIRVwEfJ6mUngDmkbRRQNL28JE05XBKRCwjuXv9MElFuAw4lU1/y28BDgIeAc4CLqv3AdLA9gvgFcD/lOzaMf08j5OkVn5MEoxIe8rMLT9XWiG+Efh8RDxQ8nMfSQrqHRGxgSRI7U7SYL4c+If0FDcCdwIPSHoo3fZp4BmSYPFVRjfuLiBp/P0dSWrsaUan6eqq89kvS883TPJve0udc/wO+Bjwv8DvSdqeSp1Okha6JU0//S+b2jj2SF8/mZbjixHxo0bLb+0nt+GYmRn4CcHMzFIOCGZmBjggmJlZygHBzMyALhyHsMMOO8S0adOKLoaZWVdZtGjRQxExWO+YrgsI06ZNY+HChUUXw8ysq0iqN1sA4JSRmZmlMg0IShY3+bWk2yRV3Nano1I/J+keSXconXrXzMzyl0fK6NCIeKjGvqNIRivuQTK69KL0v2ZmlrOiU0bHA5elUx7fAkyWV5cyMytE1gEhSKa2XSRpTpX9Q4yee2V5um0USXMkLZS0cNWqVRkV1cysv2UdEGZGxAEkqaGTJP1t2f5q0+lWTK4UERdHxIyImDE4WLfXlJmZNSnTNoR0cRAi4kFJ3yZZLeonJYcsZ/Sc+VOAFVmWycxsvOYtGeaCBXezYvVadp48kVOP2JPZ0yuSGV0vsyeEdK73bUZ+J5l3/Tdlh10LnJD2NjoYeCydktfMrCPMWzLMmdf8muHVawlgePVazrzm18xbUrHmT9fLMmX0V8BPJd1OsijI9yLi+5JOlHRiesx1wL0k86V/CXhvhuUxMxu3Cxbczdp1o9ftWbtuAxcsuLugEmUns5RRRNwL7Fdl+9yS3wM4KasymJm1asXqtePa3s2K7nZqZtbRdp48cVzbu5kDgplZHacesScTByaM2jZxYAKnHrFnjXd0r66b3M7MLE8jvYn6oZeRA4KZ2RhmTx/qyQBQzgHBzKxD5T3+wQHBzKwDjYx/GOnyOjL+AcgsKLhR2cysAxUx/sEBwcysAxUx/sEBwcysAxUx/sEBwcysAxUx/sGNymZmHaiI8Q8OCGZmHSrv8Q9OGZmZGeCAYGZmKQcEMzMDcggIkiZIWiJpfpV9h0h6TNJt6c9Hsy6PmZlVl0ej8geApcBza+y/OSJm5VAOMzOrI9MnBElTgGOAS7K8jpmZtS7rlNFngNOAjXWOebmk2yVdL2mfagdImiNpoaSFq1atyqSgZmb9LrOAIGkW8GBELKpz2GJg14jYD/g8MK/aQRFxcUTMiIgZg4ODGZTWzMyyfEKYCRwn6X7gCuAwSZeXHhARj0fEk+nv1wEDknbIsExmZlZDZo3KEXEmcCYkvYmAUyLibaXHSNoR+HNEhKQDSQLUw1mVycysWY0uVtP2RW1+/3t40YuS31euhB13bP5cY8h96gpJJwJExFzg9cB7JK0H1gJviojIu0xmZvU0ulhNWxe1uekmOOyw0du2377JT9CYXAamRcSPRrqWRsTcNBgQERdGxD4RsV9EHBwRP8+jPGZm49HoYjVtWdTmF78AaXQwmDsXImBgYNxlHw9PbmdmNoZGF6tpaVGbxYvhpS8dve3ww+GGGxoqYzs4IJhZrvJeOL4ddp48keEqlXr5YjWNHjfKnXfCvvuO3nbggXDrrU2VtRWey8jMcjOSYx9evZZgU4593pLhootWV6OL1YxrUZvf/z5JDZUGg732SlJDBQQD8BOCmeWoXo49q6eEdjyRNLpYTUPH/fGPMG3a6AvstBOsWDHuz9ZuDghmlpu8F45vZ6+fRherqXncihUwVLZ94kRYs2Zc5ciSA4KZ5aapHHsTRp4Kql0r6yeSCqtWwfOfX7m9A3vYuw3BzHKTx8Lxpe0UtWT1RDLK6tVJG0F5MNi4sSODAfgJwcxylMfC8dXaKcq1+4lklEceqT6AbMMG2Kyz78EdEMwsV1kvHN/I3f+aZ9Yzb8lwe8vx+OOw7baV29evhwkTKrd3oM4OV2Zm49TI3f+ja9a1r7vrE08kqaHyYPD000lqqEuCATggmFmPqdZOoSrHjXtKiXJPP50EgueWLQb51FNJINhyy+bPXRCnjMysp5S2UwyvXssEiQ01GnGbalxetw622KJy+0MPZT75XNb8hGBmPWf29KFnnxRqBQMYZ+Pyhg3JE0F5MBgeTp4IujwYgJ8QzDpGN87x08nG6m3UcHfXiOq9g/7wB3jBC1ooYedxQDDrAG2dR9+A+umgoUYDrqq0Ptx6azL5XA9yQDDrAEXM8dNJsng6qjUqemjyRH52xmFV3lGiWiC4/no48siWytTpMm9DkDRB0hJJ86vsk6TPSbpH0h2SDsi6PGadqB1z/MxbMszM825ktzO+x8zzbuz4GURHZDUD6qF7DVb0LhozTSRVBoMrrkjSRj0eDCCfRuUPAEtr7DsK2CP9mQNclEN5zDpOrcbNRhs9u3VaaWjTKmNl5i0Z5upFw5Q2Jwt43UtrDIqrFghGVin7h39ouhzdJtOAIGkKcAxwSY1Djgcui8QtwGRJO2VZJrNO1OocP1lUqnnJYgbUat9HADfdtWr0gdUCwSc+kQSCf/mXpq/frbJuQ/gMcBqwTY39Q8CyktfL020rSw+SNIfkCYKpU6e2v5RmHWDLzTd7thLbbtIAZx27T8N59Dymlc6qF1S9GVCbveaY30e1NoJTToELLhhX2XtNZk8IkmYBD0bEonqHVdlW0Wk4Ii6OiBkRMWNwcLBtZTTrBCPpntVr1z277dE16zj72jsbTvm0mnJqtIxZpKRqPR0dutdg09es9bnvO39WZTB4xzuSJ4IODQZ5tg1lmTKaCRwn6X7gCuAwSZeXHbMc2KXk9RSg+GWDzHJUq7/86rWNz7dTa7qGQ/dqzw1Ulimp2dOHOPe1L2Fo8kRE0gvo3Ne+hJvuWtX0Ncu/j/vPn8X9588afdCRRyaB4CtfafkzZCXvtqHMUkYRcSZwJoCkQ4BTIuJtZYddC7xP0hXAQcBjEbESsz5SL63TaNfT2dOHWPjHR/j6LX969hE7gKsXDTNj1+e1nNrJOiVVbQbUD/3PbU1fc+Rcsw+YUrnzuc+Fxx4bfyELkHd35NynrpB0oqQT05fXAfcC9wBfAt6bd3nMijZWWqfRSvemu1ZV5FvbdRefdUqq7deUqgeDiK4JBpD/kqO5BISI+FFEzEp/nxsRc9PfIyJOiogXRsRLImJhHuUx6yTV0j2lGq10s6w88ljprC3XrNZrCJJA0KGrlNWTdyD25HZmbdRMA+BIDn27SQMV+8ZT6WZZedTK82c5inpc1+yxQDAi70DsqSusL+QxcVwr8xGN5NBbKeehew1y+S1/qrq9WUVPuDfm6mrVggB0dRAolceSo6UcEKzn5TVxXDsaAFtZXrJi0NUY28fSCRPu1QxIPR4ISmW95Ggpp4ys5+U1ijfvBsCsr1/06OdqXS5nHzClJ1NDncJPCNbz8qqo6424zUO7r190gCsNSBVjCEY4CLSVA4L1vKwq6vJ0xqF7DXL1ouFRd9VZ98QpdeoRe45K8bR6/aID3IrVax0IcuaUkfW8LHpqVEtnXL1omNe9dCjXnjil2t0TqIiups+Skmkmykw7fT4zz/1h9tfvU35CsJ6XRU+NWvn1m+5aNfbiKxlqZwNk3j1cgJqNxdNOT5ZTyfOJqx85IFhfaLSibLSbZdH59bzk1sNl661hzZqKzfMWL+eCBXcjrzOdCwcEs9R4ulkWnV/vGfvtB3fcUbk9bSOYjdeUzpPbEMxS4+lmWWh+vRe85jVJeqg8GLj7aKH8hGCWGk8aqJD8ei9417vg0ksrtzsIdAQHBLPUeNNAeY4g7WQNtbuccw6cfXblmzdurD3q2HLnlJFZymmg8RtzAZcvfCGp8MuDwYYNyVOBg0FH8ROCWcppoPGr1e6y+JMXMfuKj1e+4emnYcstxzxv0ZPq9avMAoKkrYCfAFum1/lWRJxVdswhwHeA+9JN10TEx7Iqk9lY8kwDtaPSK7riLG9fecX9t/GN//lI5YGPPALbbdfQOTthUr1+leUTwl+AwyLiSUkDwE8lXR8Rt5Qdd/PI4jlm/aIdlV4nVJwj7S57//lervvK+ysPuP9+2HXXcZ0z72UjbZPM2hDS1dCeTF8OpD/uSmBGe2YSLXo2UoCz9p3I/efPqggGN31zQdJGMM5gAP0z6K8TZdqGIGkCsAjYHfhCRNxa5bCXS7odWAGcEhF3VjnPHGAOwNSpUzMssVk+2lHpFVpxPvIIbL89ry7bfPI/foK/ef8JLd3Jd9Ogv6JTdu2WaS+jiNgQEfsDU4ADJe1bdshiYNeI2A/4PDCvxnkujogZETFjcLD51Z/MOkU7lrssYuF71q5NegZtv/3o7V/8IkTwqf8+s+UKsVt6e43Zw6oL5dLtNCJWAz8Cjizb/vhIWikirgMGJO2QR5nMitRqpTdvyTBrnllfdd+aZ9a3v1IaGS8wadLo7aeckqSG3vOetl2qiPWbm9EJKbt2y7KX0SCwLiJWS5oIHA6cX3bMjsCfIyIkHUgSoB7OqkxmnaKVLq7ljcnlHl2zrr2Ny9XGCsyeDd/+duvnrqEbBv31YltHlm0IOwFfTdsRNgOujIj5kk4EiIi5wOuB90haD6wF3hThMezWH5qt9KrdmZZrS6+caoFgaAiWL2/+nD2km9o6GpVZQIiIO4DpVbbPLfn9QuDCrMpg1m0aaaRs9A606TvVPlrAvhXtXqGuE3ikslmHaHRcQa0703LjvlN1IBiXXhzZ7oBglpOx7v4bHZBV7c60nNLjGuJA0LRuaOsYDwcEy12v9d1uRCN3/402UpbemQ6vXosYPeJTwFsPnjr2d+pAYGUcECxXnTDdQhEaufsfTyNl6Z3puANsDoGgH4N+L3BAsFz16zw1jdz9N9tI2XDaYpyBoNlKvR+Dfq8EQK+HYLnqxb7bjWhkVHFmA7Kk6sGgznKVrYzC7cUBW/X00ohlPyFYrnqx73YjGr37b2sjZQupoVae5Pot6PfSU6+fECxX3TJPTbvlOh1DE08E5Vqp1AuZY6lAvRQA/YRguerFvtuNyryLYhsbi1t5kuvFAVv19NJTrwOC5a7X+m63S9MNkxn0GmqlUu+3oN9LAdABwbpKO3tzdFLPkKZ65mTYfbTVSr2fgn4vBUB121xyM2bMiIULFxZdDCtAtVk+Jw5MaCoXX2vG0O0mDXDWsftk/j9zeTBa88x6Hl2zruK4ockT+dkZh43e6AFl1gRJiyJiRr1j3KhsXaOd3RlrzRg6MnV0ll0Gq3VTrBYMoKxhsg2NxWb1OGVkXaOdvTnqvSfrLoONTF89YufJE/vuiaCTUnn9xgHBukY7e3OMNWNoll0GGz33/efPqr6jgUDQrZVqP45y7iSZpYwkbSXpl5Jul3SnpHOqHCNJn5N0j6Q7JB2QVXms+7VzDEO1c5XKsstgrXNPnjjA0OSJ3H/+rKrBYOa5P2Te4rEXp+nmkbP9Nsq502T5hPAX4LCIeFLSAPBTSddHxC0lxxwF7JH+HARclP7XrEI7e3OMvOfsa+9k9drR+fusuwzW6qZ429lHVD1+2unzk18avFvu5pGzvTTIqxtluWJaAE+mLwfSn/Jn3eOBy9Jjb5E0WdJOEbEyq3JZd2tnd8aRc+WdXikPbPfVSA3NPPeHFWmtRir2bq5Ue2mQVzfKtA0hXU95EbA78IWIuLXskCFgWcnr5em2UQFB0hxgDsDUqVMzK68Vo+h8d1F95n925quq70jbCFac8b2qu8eq2Lu5Uu2lQV7dKNNupxGxISL2B6YAB0rat+yQat0nKlrMIuLiiJgRETMGBwezKKoVpJvz3U2TmH3AlIrN8xYvH9Vg3OycQFnMFzVvyTAzz7uR3c74HjPPuzGzf59c53yyCrn0MoqI1ZJ+BBwJ/KZk13Jgl5LXU4AVeZTJOkM357vHbe+9YenSis3TTvsuSAw1sFRmo+sjQPtGzrbS86eZp79+GuXcaTILCJIGgXVpMJgIHA6cX3bYtcD7JF1B0pj8mNsP+ku9fHfRqaS2ectb4JvfrNj8wlO/w4bNNt3J11sqc7zfQTsr1WaDtruQdp8snxB2Ar6atiNsBlwZEfMlnQgQEXOB64CjgXuANcA7MyyPdaBa+e7JkwZyrUwyCT7nnANnn12x+fCz53NPlTg41lKZzWjH52q2kbqvnv56xJgBQdL7gK9HxKPjOXFE3AFMr7J9bsnvAZw0nvNab6mWFhmYIFavWVfRmJRVZdL2O9mvfQ1OOKFi80Hv/Sqb7zKFQ/caZHjRcOYNp+36XM02Undzb6d+1cgTwo7AryQtBi4FFkS3zYhnHas8LTJ50gBPPr2+smdBqrQyadddfSspkdLrn7/9w7xyzhsrjjtuzhe5Y7u0d9zqtVy9aJjXvXSIm+5alVk6bN6SYU6+8nY2lP2v2kxQbbYtY/KkgapzNE2eNNDwtS1fYwaEiPiIpH8HXk2S0rlQ0pXAlyPiD1kX0HpfaVpk5nk31pzoDTbdlbbzrr6ZO9nS67/g4eXceP6JlQctWMDMxZtXHUtw012rKmcxbZORspUHgxHjvUNvti2j1m2jbyc7V0NtCBERkh4AHgDWA9sB35J0Q0SclmUBrb/Uq6xK70rbmZ9uJiVywYK7mfjYIyz9/Fsrd15yCbzrXQCsuLG5sQStGGvyvGbGIzTTlvHY2uqBvdZ2K14jbQjvB94BPARcApwaEeskbQb8HnBAsLapVTlPkEb1R29nfnrcKZG1a6sOKvvMzDfz2Ve+lfvedcyz24oYJNZoUG238hTathMHKqYFge4YINevGhmYtgPw2og4IiKuioh1ABGxEagxHaNZc2oNqvrUG/cbdYfazoXcGx4MtWFDMhX1pEmjNl/74r9l2unz+cwr31px/SwGiY2l1ndQHlSraXYAWrUBhk89s56BzUaPPfWo487WSBvCR+vsqxxlY9aCRvPV7Z7iYMyUSJU1CW7feU+Of/un6l6/iOUVa303jQSDZttlqqWp1m0Itps0wKQtNu/+sSR9wktoWtfKZeBatcVpttqKeT+/Z9RMqeNZejOPcjdzjZnn3Vg1vVV1Gc8yu53xvao9wwTcd94xVfZY3hpZQtML5FjXynSKgzqrlFVbj/npdRsbOm1eo3eb+W5aaZfp5gn1bBOvqWxWqoF1i1tZxKWTF4BppV2miLYSaz8HBDMY1wL2rdxJd/Lo3VYqdc9S2hucMrL+1sQC9q2kR8Z6b6vtC628v9UGcM9S2v3cqGz9qYlAMKJaG0IjvXjGei/Q9HlbLZf1vkYalZ0ysv4yjtRQLa2kR+q9t9X2hU5un7Du4JSR9YcWngiqaSU9Uuu9rbYvdHL7hHUHBwTrbW0OBJDdOILxtE1UK4O7flqrnDKy3tSG1FA1Wa4B3Wgvn1plOHSvQXf9tJZkFhAk7SLpJklLJd0p6QNVjjlE0mOSbkt/ak6TYdaQNgaCavP6ZJmnb7RtolYZbrprlbt+WkuyTBmtB06OiMWStgEWpdNl/7bsuJsjwpPkWWvanBqqNaK41rTS7crTN9I24bYCy0pmASEiVgIr09+fkLQUGALKA4JZ8zIIBBcsuLtqLn7tug1MkKouPJNnnr5T1qG23pNLG4KkaSTrK99aZffLJd0u6XpJ+9R4/xxJCyUtXLVqVYYlta6RQRtBaW6+lg0RFXn6gQniqb+sH/eU0c2q1dYQgbudWksyDwiSngNcDXwwIh4v270Y2DUi9gM+D8yrdo6IuDgiZkTEjMHBwWwLbJ0to8ZiGHulMdiUlx/J0283aQACVq9d1/ZG5lpqtTXUWolsePXa3IKVdbdMRypLGgDmAwsi4r8aOP5+YEZEPFTrGI9U7lMZdB8tV2sK5xHVRv22MmV0Pe2cvrqURy73r0JHKksS8GVgaa1gIGnH9DgkHZiW5+GsymRdKMMngnL12gFq9djJooG32a6t1VJJ5ZxCsnqyTBnNBN4OHFbSrfRoSSdKOjE95vXAbyTdDnwOeFN02+RKlo0cA8GIU4/Yk4EJldcc2Ew179DbuZTniGa7tpankmpxbySrJbOAEBE/jQhFxF9HxP7pz3URMTci5qbHXBgR+0TEfhFxcET8PKvyWGcb6fNfRCAYMXv6EFtvUdnxbt3GqFkZZ7EOQCtPHbOnD/GzMw7jvvOOYSiDYGW9zSOVrXDzlgxz3AG78LMzX1W5M4dAUKpWw2ytyjiLdQDa9dThRWtsvDyXkRVr332ZfeedFZunnT4/aZjNuTjNzAfU7nUATj1iz6rTWI+3Im91fQPrPw4Ils9i9eVe8xqYV9nLeNrp85/9vYhcd7sq41a0syL3ojU2Hg4IfS6vRd+f9ZGPwMc/XrF52mnfrWg7aCRF0u5g1il31a7IrQgOCH2uXo+WtlZIc+fCe95TuX39eubd8QATm7grH08wG0/gcGVs/coBoc9lPlHad74Ds2dXbl+7FrbaCmj+rnys7pkj55s8aYAnn17Puo1J47Tn+DGrzgGhz2W2qMovfwkHHVS5/dFHYfLkis3N3JXXClrlM5M+uqay51AmT0FmXc7dTvtc27sm/v73SVtAeTBYtizpPlolGDSrVtCaII05JxGM/RRUbT0Es17mgNDn2taP/oEHkkDwoheN3r50aRIIpkxpW5lH1Apm1aanrqbeU1CWK6OZdSqnjKy1RtSnnoLnPKdy+y9+AQcf3FrBxlCr7aHWegalxnoKyq2x3ayDOCBYc9atgy22qNx+7bVw7LG5FaNWMCsfSzAwQWy9xeY8tnZdQ43WXpXM+pEDgo1PBGxWJdN45ZXwhjfkX54q2jGWILPGdrMO5oBgjas26dwXvgDvfSDu0rUAAAz7SURBVG/+ZRlDq2MJOmHEslneHBBsbNUCwSc/Caeemn9ZctIpI5bN8uSAYLVVCwTnnAMf/eioTe2aPqKQOZXq8Ihl6zeZBQRJuwCXATsCG4GLI+KzZccI+CxwNLAG+MeIWJxVmaxB1QLBe9+bpIfKtGsupNznVCq7dicFIrOiZDkOYT1wckS8GDgYOEnS3mXHHAXskf7MAS7KsDw2lmqL07zmNUlDcpVgAM2v7tXoec6+tnJq7HbyeAOzTbJcMW3lyN1+RDwBLAXKb7uOBy6LxC3AZEk7ZVWmXtPqSNq6q5S97GVJILjmmrrnaFf3zFrHr167LtPKuV0BzawX5NKGIGkaMB24tWzXELCs5PXydNvKsvfPIXmCYOrUqVkVs6u0mmKZt2SY2QdMoXzauae3H2Srhx5suBzt6p5Z6zxAw4PBmkn9eLyB2SaZT10h6TnA1cAHI+Lx8t1V3lIx70BEXBwRMyJixuDgYBbF7Dot3dlKzD6gciqJaafP51WnXDGucrRrLqR6xzdSOTeb+mnXcpVmvSDTgCBpgCQYfD0iquUelgO7lLyeAqzIsky9otbddN3Ks8YC9tNOn//sSmXjvTNu11xIs6cPsd2kgar7Gqmcmw2QXnfYbJMsexkJ+DKwNCL+q8Zh1wLvk3QFcBDwWESsrHGspeYtGUZUeZSiRuVZrdcQo5errPv+MdTqnjneFM5Zx+7T9GCwZlM/Hm9gtkmWbQgzgbcDv5Z0W7rtw8BUgIiYC1xH0uX0HpJup+/MsDw944IFd1cNBqIs9VIjEBDBvCXDTa1S1qhm2jhaqZxbacvweAOzRGYBISJ+SvU2gtJjAjgpqzL0qlp3vUFaqdYJBCNaqXwbufNvdrbQZitnTzVh1jqPVO5Cte6Gv/f1U+D8WZVvqLE+QDOVb6N3/nn33nHqx6x1DghdqPxu+NKrzuawexdWHtjgQjHj0eidfxGzhTr1Y9Yar5jWhUZ69py55BruP39WZTCIyCQYQP11jEsHx7n3jln38RNCN7r+emYffXTl9oyCQKl6A8iqpY+cwjHrHoocKpF2mjFjRixcWCU90mWamlDtxz+GQw4ZvW377eGhhzIrZ7nyNoRqhiZP5GdnHJZbmcxsbJIWRcSMesc4ZVSAcY+qveWWpOdQaTAYmXQux2AAowei1eJpH8y6kwNCARoeVbtkSRIIXv7yTdsOP7yhSeeyUvpkM6FG91ZP+2DWndyGUIAxu2T+9rewzz6jdx54INxaPjdgvsrTRRuqpBvdcGzWvRwQClCrYfZl6x+pHFS2555w1105lay+ak82ABMkNka44disyzkgFKB8HMHOjz/Izy/6p9EH7bgjrOysaZ1qPdlsjOC+847JuTRm1m4OCAUYuYO+9Kqfc+25bxy9c+JEWLOmgFKNrYjBZmaWHzcqF2HNGmbPmV0ZDDZu7NhgAB5sZtbr/ISQp6efhle/Gm6+efT2jRtrT0jXQTzYzKy39U1AaGogWLs88wzMmgU33LBp27/+K3z2s10RCEp5viCz3tUXAaHV9Yebtm4dvO518N3vbtr27nfD3LmwWX9k6woNxGY2LpnVSpIulfSgpN/U2H+IpMck3Zb+fDSrsrS0/nAzNmyAN74RtthiUzB4+9th/Xq4+OK+CgbNrHNsZsXIsmb6CnDkGMfcHBH7pz8fy6oguc3Nv3EjnHACbL45XHVVsu0Nb0ieFC67DCZMqP/+HpN7IDazlmQWECLiJ8AjWZ1/PGp1i2xbd8mNG2HOnKTC/9rXkm3HHpu0HVx5ZRIg+lDei+SYWWuKzl28XNLtkq6XtE+tgyTNkbRQ0sJVq1aN+yKZdZeMgPe/PwkEX/pSsu3v/z7pTXTttTAw0Nr5u1zmgdjM2qrIgLAY2DUi9gM+D8yrdWBEXBwRMyJixuDg4LgvVDpDp0imZz73tS9pvnEzAk47LWkL+Pznk20zZyZjCH7wA9hyy+bO2wbzlgwz87wb2e2M741asKYIHrdg1l0Ky2VExOMlv18n6YuSdoiITOZzblt3ybPOgo+VNHe89KXJOgVbb936uVtUWG+qGjxuway7FBYQJO0I/DkiQtKBJE8rDxdVnjGdey58+MObXu+zD/ziF7DNNsWVqUyj6x3nyeMWzLpHZgFB0jeBQ4AdJC0HzgIGACJiLvB64D2S1gNrgTdFJy7f9ulPw7/926bXL3whLFwIkycXV6Ya3IhrZq3ILCBExJvH2H8hcGFW12/ZF78IJ5206fXQENx+e7JkZYfy5HNm1oqiexl1nksvTaaTGAkG228Pf/4zLF/e0cEA3IhrZq3pzw7y1Vx+eTKaeMTWW8Pvfgc771xcmcbJjbhm1goHhKuuSqaZGDFhAtx3H+yyS3FlakF5UBgZFeygYGZj6d+A8J3vwOzZo7fdey/stlsx5WmTTut6ambdo//aEL7//aSNoDQY/O53yWCzLg8G4PmDzKx5/fWE8LznwaOPbnp9552w997FlScD7npqZs3qnyeEiE3tArfdlrzusWAAnj/IzJrXPwFBSsYRRMB++xVdmsy466mZNau/UkZ9wF1PzaxZDgg9yPMHmVkz+idlZGZmdTkgmJkZ4IBgZmYpBwQzMwMcEMzMLOWAYGZmQIYBQdKlkh6U9Jsa+yXpc5LukXSHpAOyKouZmY0tyyeErwBH1tl/FLBH+jMHuCjDspiZ2RgyCwgR8RPgkTqHHA9cFolbgMmSdsqqPGZmVl+RbQhDwLKS18vTbRUkzZG0UNLCVatW5VI4M7N+U2RAUJVtUe3AiLg4ImZExIzBwcGMi2Vm1p+KDAjLgdJ1KqcAKwoqi5lZ3ytycrtrgfdJugI4CHgsIlYWWJ6eNm/JsGdANbO6MgsIkr4JHALsIGk5cBYwABARc4HrgKOBe4A1wDuzKku/8zrLZtaIzAJCRLx5jP0BnJTV9W2TeussOyCY2QiPVO4DXmfZzBrhgNAHvM6ymTXCAaEPeJ1lM2uEl9DsA15n2cwa4YDQJ7zOspmNxSkjMzMDHBDMzCzlgGBmZoADgpmZpRwQzMwMcEAwM7OUkimFuoekVcAf23S6HYCH2nSubuPP3p/82fvTDsDWEVF3QZmuCwjtJGlhRMwouhxF8Gf3Z+83/uxjf3anjMzMDHBAMDOzVL8HhIuLLkCB/Nn7kz97f2ros/d1G4KZmW3S708IZmaWckAwMzPAAQFJF0i6S9Idkr4taXLRZcqSpCMl3S3pHklnFF2ePEnaRdJNkpZKulPSB4ouU54kTZC0RNL8osuSN0mTJX0r/X99qaSXF12mPEj6UPq3/htJ35S0Vb3j+z4gADcA+0bEXwO/A84suDyZkTQB+AJwFLA38GZJexdbqlytB06OiBcDBwMn9dnn/wCwtOhCFOSzwPcjYi9gP/rge5A0BLwfmBER+wITgDfVe0/fB4SI+EFErE9f3gJMKbI8GTsQuCci7o2IZ4ArgOMLLlNuImJlRCxOf3+CpFLoi1WDJE0BjgEuKboseZP0XOBvgS8DRMQzEbG62FLlZnNgoqTNgUnAinoH931AKPNPwPVFFyJDQ8CyktfL6ZMKsZykacB04NZiS5KbzwCnARuLLkgBXgCsAv47TZldImnroguVtYgYBv4T+BOwEngsIn5Q7z19ERAk/W+aQyv/Ob7kmP9DklL4enElzZyqbOu7fseSngNcDXwwIh4vujxZkzQLeDAiFhVdloJsDhwAXBQR04GngJ5vP5O0HUkGYDdgZ2BrSW+r956+WFM5Ig6vt1/SO4BZwKuitwdmLAd2KXk9hTEeIXuNpAGSYPD1iLim6PLkZCZwnKSjga2A50q6PCLqVg49ZDmwPCJGnga/RR8EBOBw4L6IWAUg6RrgFcDltd7QF08I9Ug6EjgdOC4i1hRdnoz9CthD0m6StiBpYLq24DLlRpJI8shLI+K/ii5PXiLizIiYEhHTSP7Nb+yjYEBEPAAsk7RnuulVwG8LLFJe/gQcLGlS+rf/KsZoTO+LJ4QxXAhsCdyQfGfcEhEnFlukbETEeknvAxaQ9Di4NCLuLLhYeZoJvB34taTb0m0fjojrCiyT5eNfga+nN0L3Au8suDyZi4hbJX0LWEySDl/CGFNYeOoKMzMDnDIyM7OUA4KZmQEOCGZmlnJAMDMzwAHBzMxSDghmLUhnUL1P0vPS19ulr3ctumxm4+WAYNaCiFgGXAScl246D7g4Iv5YXKnMmuNxCGYtSqfDWARcCrwbmJ7OJmvWVTxS2axFEbFO0qnA94FXOxhYt3LKyKw9jiKZYnjfogti1iwHBLMWSdof+HuSVdg+JGmngotk1hQHBLMWpLNIXkSytsKfgAtIFiUx6zoOCGateTfwp4i4IX39RWAvSX9XYJnMmuJeRmZmBvgJwczMUg4IZmYGOCCYmVnKAcHMzAAHBDMzSzkgmJkZ4IBgZmap/w+MSUi2J9bVugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)   # scatter plot showing actual data\n",
    "plt.plot(x, yhat, 'r')     # regression line\n",
    "plt.title('Predicted vs. Actual values')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('y')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the Pearson correlation between $\\hat{\\mathbf{y}}$ and $\\mathbf{y}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearson's r = 0.790\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import pearsonr\n",
    "\n",
    "corr, _ = pearsonr(yhat, y)\n",
    "print(\"Pearson's r = %.3f\" % corr)"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}