L06-1-IntroToRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"instructions\" style=\"border-radius: 5px; background-color:#f5f5f5;\" >\n",
    "<h1>Instructions</h1>\n",
    "<p>Look for the <b>3 Your Turn</b> sections to complete the code and/or answer questions.<p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 6 - Introduction to Regression\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The method of regression is one of the oldest and most widely used analytics methods. The goal of regression is to produce a model that represents the ‘best fit’ to some observed data. Typically the model is a function describing some type of curve (lines, parabolas, etc.) that is determined by a set of parameters (e.g., slope and intercept). “Best fit” means that there is an optimal set of parameters according to an evaluation criteria we choose.\n",
    "\n",
    "A regression model attempts to predict the value of one variable, known as the **dependent variable**, **response variable** or **label**, using the values of other variables, known as **independent variables**, **explanatory variables** or **features**. Single regression has one label used to predict one feature. Multiple regression uses two of more feature variables. \n",
    "\n",
    "Virtually all machine learning models, including some of the latest deep learning methods, are a form of regression. These methods often suffer from the same problems, including overfitting and mathematically unstable fitting methods. \n",
    "\n",
    "Linear regression is the foundational form of regression. To solve linear regression, normally the **method of least squares** is used.  Here, we are looking to find parameters of the function that minimizes the square of the error between the predictions, which are the function output, and the observed target values.\n",
    "\n",
    "## History\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/a/a3/Tobias_Mayer.jpg\" title=\"Tobias Mayer\" style=\"float:left; margin-right: 7px; height: 150px;\"><img src=\"https://upload.wikimedia.org/wikipedia/commons/e/e3/Pierre-Simon_Laplace.jpg\" title=\"Pierre-Simon Laplace\" style=\"float: right; margin-left: 10px; height: 150px;\">\n",
    "Regression is based on the method of least squares or the method of minimum mean square error. The ideas around least squares or averaging errors have occured over nearly three centuries. The fist known publication of a 'Method of Averages' was by the German astronomer Tobias Mayer in 1750. \n",
    "\n",
    "Lapace used a similar method which he published in 1788.\n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/0/03/Legendre.jpg\" style=\"height: 150px; float: left; margin-right: 10px;\" title=\"Caricature of Legendre - published method of least squares\"><img src=\"https://upload.wikimedia.org/wikipedia/commons/9/9b/Carl_Friedrich_Gauss.jpg\" style=\"height: 150px; float: right; margin-left: 10px;\" title=\"Carl Friedrich Gauss - early adoptor of least squares\">The first publication of the method or least squares was by the French mathematician Adrien-Marie Legendre in 1805. \n",
    "\n",
    "It is very likely that the German physicist and mathematician Gauss developed the method of least squares as early as 1795, but did not publish the method until 1809, aside from a reference in a letter in 1799. Gauss never disputed Legendre's priority in publication. Legendre did not return the favor, and opposed any notion that Gauss had used the method earlier. \n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/0/0b/Francis_Galton.jpg\" style=\"height: 150px; float: left; margin-right: 10px;\" title=\"Francis Galton - inventor of regression\">The first use of the term **regression** was by Francis Galton, a cousin of Charles Darwin, in 1886. Galton was interested in determining which traits of plants and animals, including humans, could be said to be inherited. \n",
    "\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/b/be/Karl_Pearson%3B_Sir_Francis_Galton.jpg\" style=\"height: 350px; float: right; margin-left: 10px;\" title=\"Karl Pearson(left) and Francis Galton(right)\">\n",
    "While Galton invented a modern form regression, it fell to Karl Pearson to put regression and multiple regression on a firm mathematical footing. Pearson's 1898 publication proposed a method of regression as we understand it today. \n",
    "\n",
    "Many others have expanded the theory of regression in the 120 years since Pearson's paper. Notably, Joseph Berkson published the logistic regression method in 1944, one of the first classification algorithms. In recent times the interest in machine learning has lead to a rapid increase in the numbers and types of regression models.\n",
    "\n",
    "In this and subsequent sections, we will introduce and expand the concepts of linear regression and several methods associated with it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Linear Regression\n",
    "\n",
    "We will focus here on **linear models** which are foundational\n",
    "- Derived with linear algebra\n",
    "- Basis of many machine learning models\n",
    "- Understanding linear models is basis for understanding behavior of many statistical and ML models\n",
    "- Basis of time series models\n",
    "\n",
    "### Linear model of a straight line\n",
    "\n",
    "Let's have a look at the simplest case of a regression model for a straght line. If we have one feature and one label, there are some number of values pairs, ${x_i,y_i}$, we can define a line that best fits that data.  \n",
    "\n",
    "<img src=\"https://library.startlearninglabs.uw.edu/DATASCI410/img/ymxb.jpg\" alt=\"y_equals_mx_plus_b\" style=\"width: 450px;\"/>\n",
    "<center>**Single regression model**</center>\n",
    "\n",
    "$$where\\\\\n",
    "slope = m = \\frac{rise}{run} = \\frac{\\delta y}{\\delta x}\\\\\n",
    "and\\\\\n",
    "y = b\\ at\\ x = 0$$\n",
    "\n",
    "\n",
    "If we have a number of values pairs, ${x_i,y_i}$, we can write the equation for the line with the errors as:\n",
    "\n",
    "$$y_i = mx_i + b + \\epsilon_i \\\\\n",
    "where \\\\\n",
    "\\epsilon_i = error$$\n",
    "\n",
    "We can visualize these errors as shown in the figure below.\n",
    "\n",
    "<img src=\"https://library.startlearninglabs.uw.edu/DATASCI410/img/LSRegression.jpg\" alt=\"LS Regression\" style=\"width: 450px;\"/>\n",
    "<center>**Example of Least Squares Regression**</center>\n",
    "\n",
    "We want to solve for $m$ and $b$ by minimizing the error, $\\epsilon_i$. We call this **least squares regression**.\n",
    "\n",
    "$$min \\Sigma_i \\epsilon^2 = min \\Sigma_i{ (y_i - (mx_i + b))^2}$$\n",
    "\n",
    "There are lots of computationally efficient algorithms for finding minimums of equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A first regression model\n",
    "\n",
    "Let's give regression a try. The code in the cell below computes data pairs along a straight line. Normally distributed noise is added to the data values. Run this code and examine the head of the data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.680006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.204082</td>\n",
       "      <td>-2.282142</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.408163</td>\n",
       "      <td>-0.522395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.612245</td>\n",
       "      <td>1.340155</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.816327</td>\n",
       "      <td>1.379675</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x         y\n",
       "0  0.000000  0.680006\n",
       "1  0.204082 -2.282142\n",
       "2  0.408163 -0.522395\n",
       "3  0.612245  1.340155\n",
       "4  0.816327  1.379675"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Paramters of generated data\n",
    "n_points = 50\n",
    "x_start, x_end = 0, 10\n",
    "y_start, y_end = 0, 10\n",
    "y_sd = 1\n",
    "\n",
    "# Generate data columns\n",
    "x_data = np.linspace(x_start, x_end, n_points)\n",
    "y_error = np.random.normal(loc=0, scale=y_sd, size=n_points)\n",
    "y_data = np.linspace(y_start, y_end, n_points) + y_error\n",
    "\n",
    "# Put data in dataframe\n",
    "sim_data = pd.DataFrame({'x':x_data, 'y':y_data})\n",
    "# Or alternatively:\n",
    "# sim_data = pd.DataFrame(data=np.column_stack((x_data, y_data)), columns=['x','y'])\n",
    "\n",
    "sim_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, you can visualize these data by exectuting the code in the cell below. Notice that the points nearly fall on a straight line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'x vs y')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ4klEQVR4nO3dfZBdd13H8c83KS69XRwaSi6Qdu+Cdqid+MQuz6Nm03WmPAxF1HHxgsUx7vhQiEpGq+vIH7LKaEatD6hLKFbvwuoEkIodoUuTYaxaSQpCa8AykA2RkpSJWLbRQJuvf+y9dbN77+7ec885v3PO7/2aubN7z557z/d3b3K+5/d4zN0FAIjPttABAADCIAEAQKRIAAAQKRIAAESKBAAAkSIBAECkSAAAECkSAABEigQAAJEiASB6ZvZtZnbOzF7Qfv4cM/uqme3psu+tZnZ4zbbbzOwP27+/ycy+YGZfN7Mvmlmzy3s8y8zOm9kzVm0bM7NHzOwpaZcP6MVYCgKQzOynJf2SpDFJH5T0GXc/0GW/hqQTkp7l7o+a2XZJpyX9kKTPSHpY0gvd/XNm9mxJO9z9wS7vc5ekv3P3P20//31Jl7n7m7MpIbAeNQBAkru/S9JDku6T9GxJMz32W5J0v6TXtjftlXTe3f+l/fyipN1mdrm7P9zt5N92h6Q3SFI7ibxe0l+lUBRgy0gAwP97l6Tdkv7I3S9ssN97tXLClqQfbz+Xuz8m6cck/Yykh83s783suh7v8SFJ15vZ8yT9oKT/dvd/TaEMwJbRBARIMrNhSf8m6YikV0j6Tnc/12PfZ0o6JelaSQ9Ieqm7n1izz+WS3i7pRe7+fT3e58+10nx0naQT7v72lIoDbAk1AGDFbZKOu/s+SX8v6c967ejuj0g6Kuk9kr7YOfmbWd3MXmNmV0i6IGlZ0hMbHPMvJb1J0msktVIoA9AXEgCiZ2Y3SbpRK0030kpn8Au6jeBZ5b2SJts/O7ZJequkL0s6J+kHJP1crzdw93u10mdwv7ufTBo/kBRNQEBAZnaPpPe6+6HQsSA+JAAgEDN7oaS7JV3j7l8PHQ/iQxMQEICZ3SFpUdIvcPJHKNQAACBS1AAAIFKXhQ6gH1dddZWPjo4meu1jjz2mK664It2ACo4yx4Eyx2GQMh8/fvyr7v7MtdtLlQBGR0d17NixRK89evSo9uzZk25ABUeZ40CZ4zBImc1sqdv2zJuAzOx2MztrZg+s2rbDzO42s4faP6/MOg4AwKXy6AP4C61MslntVkkfc/drJX2s/RwAkKPME4C7f1wrsyJXu0krqyGq/fO1WccBALhULsNAzWxU0ofdfXf7+dfc/emr/v5f7t61GcjMpiVNS1K9Xh9bWFhIFMPy8rKGh4cTvbasKHMcKHMcBinzxMTEcXcfX7u98J3A7j4naU6SxsfHPWknCJ1GcaDMcaDM6Qg1D+BM+25Jav88GygOAAhufn5eo6Oj2rZtm0ZHRzU/P5/LcUMlgDsl3dz+/Wat3BwDAKIzPz+v6elpLS0tyd21tLSk6enpXJJAHsNA3yfpnyU938xOm9lPSXqHpB80s4e0cjekd2QdBwAU0czMjM6fP3/JtvPnz2tmputdSVOVeR+Au7++x59uyPrYAFB0p06d6mt7mlgLCAACGhkZ2XB7p39g7969qfcPkAAAFN7aTtLFxcXQIaVmdnZWtVrtkm21Wk2zs7OZ9w+QAAAUWreT4MGDBzc9CYYaWdOvZrOpubk5NRoNmZkajYbm5ubUbDYz7x8gAQAotG4nwQsXLmx4Egw5siaJZrOpkydP6uLFizp58qSazZXbUWfdP0ACAFBoSU6CIUfWpGmz/oFBkQAAFFqSk2DIkTVp2qh/IA0kAACF1u0kODQ0tOFJMOsr57xs1D+QBhIAgELrdhI8cODAhifBrK+c89TpH7jnnnsu6R9IAwkAQOGt7SSdnJzcdP8sr5yrovCrgQJAEs1mkxP+JqgBAEBOijY3gRoAAOSgMzehMzy1MzdBUrCaCjUAAMhBEecmkAAAlFrRmlV6KeLcBBIAgNIq05IPRZybQAIAUFp5NKukVcMo4twEEgCA3KV1Us26WSXNGkYR5yaQAADkKs2TatbNKmnXMHqt+hkKCQBArtI8qWbdrFLEjts0kQAA5CrNk2rWzSpF7LhNEwkAQK7SPqlm2ayStIZRlqGpJAAAuSriaJhektQwyjQ0lQQAIFdFHA2zkX5rGEWc8dsLCQBA7nqdVMvSdLKRMnUcB00AZvaLZvagmT1gZu8zs6eGjAdAOGVqOtlImTqOgyUAM9sl6S2Sxt19t6TtkqZCxQMgrDI1nWykTH0coZuALpN0uZldJqkm6cuB4wGQon6adMrUdLKRMvVxmLuHO7jZfkmzkv5H0kfdfd0nZGbTkqYlqV6vjy0sLCQ61vLysoaHhweItnwocxyKWubFxUUdPHhQFy5ceHLb0NCQDhw40PWWjlNTUzpz5sy67fV6XWv/3xe1zFkapMwTExPH3X183R/cPchD0pWS7pH0TElPkfS3kt6w0WvGxsY8qSNHjiR+bVlR5jgUtcyNRsMlrXs0Go2u+7daLa/VapfsW6vVvNVqrdu3qGXO0iBllnTMu5xTQzYBTUr6ors/4u7flPQBSS8LGA+AFPXbpFOmppOqCHlLyFOSXmJmNa00Ad0g6VjAeACkaGRkREtLS12398KN3PMVrAbg7vdJOizpfkmfaccyFyoeAOkq02iYWAUdBeTub3P369x9t7u/0d0vbP4qAGVAk07xhWwCAlBxNOkUW+h5AACQqyosN5EWagAAotFZbqIz47iz3ISkKGsq1AAARCPt5SbKXpugBgAgGmkuN1GF2gQ1AADR2Gilzn6v5quweB0JAMDAytIU0mtuwitf+cq+l6KuwuJ1JAAAAynTOv695ibcddddfV/Nl2nd/15IAAAGUramkG53I0tyNV+Fmc4kAAADqUJTSJKr+SrMdCYBABhIFZpCkl7N93vD+KIhAQAYSBWaQqpwNZ8E8wAADKRzkpyZmdGpU6c0MjKi2dnZ0p08Y1y3iAQAYGAxnjyrgCYgAIgUCQCouLJM0kL+aAICKqwK69UgO9QAgAor2yQt5IsEAFRY2qtf0pRULSQAoMLSmqRVpvV+sHUkAKBE+r0KTzJJq9sxaEqqJhIAUBK9rsIXFxd7vqbfGa69jrG0tNR1/zKt94P1SABASfS6Cj906NCGr+tnvZpex9i+fXvX/cu03g/WIwEAJdHravvs2bOJ3q9bU0+vYzzxxBOlX+8H6wVNAGb2dDM7bGafNbMTZvbSkPEARdbranvnzp19v1evpp4dO3Z03b/TdBTbYmlVF7oGcJukf3D36yR9t6QTgeMBCqtXh+6+ffv6fq9eTT2d91x7jM7ibmVe+hjrBUsAZvatkr5f0rslyd2/4e5fCxUPUHS9OnQnJyf7fq9eTT3nzp3jSj8i5u5hDmz2PZLmJP27Vq7+j0va7+6PrdlvWtK0JNXr9bGFhYVEx1teXtbw8PAgIZcOZY5DkjJPTU3pzJkz67bX63Ul/T+WJ77n/kxMTBx39/F1f3D3IA9J45Iel/Ti9vPbJP3mRq8ZGxvzpI4cOZL4tWVFmeOQpMytVstrtZpLevJRq9W81WqlH2AG+J77I+mYdzmnhuwDOC3ptLvf135+WNILAsYDlFo/k8RivQMWLhVsNVB3/4qZfcnMnu/un5N0g1aagwD0Kcmqn9zEBaFHAb1Z0ryZfVrS90j6rbDhAOXEUg1IIuj9ANz9U1rpCwAwgDRX/UQ8QtcAAKQgrVU/ERcSAFABSVb9BEgAQAUwqgdJcE9goCIY1YN+UQMAAuI2iwiJGgAQSJKx+0CaqAEAgTB2H6GRAIBAGLuP0EgAQCCM3UdoJAAgEMbuIzQSAJCSfkf0MHYfoTEKCEhB0hE9jN1HSNQAgBQwogdlRAJAKSWZQNXrNWlMxmJED8qIJiCUzkbNLbt27errNffee6/uuOOOgSdjjYyMaGlpqet2oKioAaB0kjS39HrN3NxcKk03jOhBGZEAUDpJmlt6/e2JJ57o+726YUQPyogmIJROkuaWXq/Zvn171ySQpOmGET0oG2oAKJ0kzS29XjM9PZ1L0w2rfqKISAAonY2aWxYXF7ueaHu95p3vfGfmTTedDuilpSW5+5MdzSQBBOfupXmMjY15UkeOHEn82rKKrcytVsuHhoZc0pOPWq3mrVYraFyNRuOSmDqPRqORyvvH9j27U+Z+STrmXc6p1ABQGTMzM7pw4cIl24owGStJpzVNRsgDCQCVUdTJWP2u+kmTEfJCAsCmynI1mtfyyv1+Hv12WrOsBPISPAGY2XYz+6SZfTh0LFivTFejs7OzGhoaumRb2iN6knwe/c4RKGpNBtUTPAFI2i/pROgg0F1eV6Np1DKazaYOHDiQ6YiepJ9Hs9nUyZMndfHiRZ08eXLDmLhRDPISNAGY2dWSXiXpUMg40FseV6Np1jImJye3fKJNIo/Pg2UlkBdbGSEU6OBmhyX9tqSnSTrg7q/uss+0pGlJqtfrYwsLC4mOtby8rOHh4QGiLZ80yjw1NaUzZ86s216v15X0u8jyGFl/z3l8HpK0uLioQ4cO6ezZs9q5c6f27dunycnJrvvybzsOg5R5YmLiuLuPr/tDt7GheTwkvVrSO9u/75H04c1ewzyA/qRR5lar5bVaLdOx9WbWdZy8mfX9Xll/z3l8Hv3i33YcqjYP4OWSXmNmJyUtSNprZq2A8ZRS1iN08ljkrKht3t0+WxZ9Q6V0ywp5P0QNIJHNrkbLUuY0r6rTKnMRr/R7Kcv3nCbK3B8lrQGY2S1mdmU26QeDqMp48SJeVVflswU2spXloJ8l6RNmdr+k2yV9pJ1RUuPuRyUdTfM9Y1Cl8eJFW0q5Sp8t0MumNQB3/3VJ10p6t6Q3SXrIzH7LzL4t49iwiaK2nW+EWcVAcWypE7h9xf+V9uNxSVdKOmxmv5NhbNhE2caLl21WcZk+WyCJrfQBvMXMjkv6HUn3SvpOd/9ZSWOSfjjj+LCBIradbyRJu3qoGkPZPlsgia30AVwl6XXufsn99Nz9opmtm7iFfBWt7Xwj/bard2oMnaTRqTFIyqXMZfpsgSS20gfwG2tP/qv+xho+2LJ+29UZiQNkqwiLwSES/barMxIHyBYJALnpt12dkThAtkgAyFU/yyIzEgfIFgkAhcVIHCBbJIAAyjIZqgj6qTEA6M9WhoEiRaGHNgJABzWAnDG0EUBRkAByxtBGAEVBAsgZQxsBFAUJIGcMbQRQFCSAnDG0EUBRkAACyGNo4/z8vKamphhqCqAnhoFWEENNAWwFNYAKYqgpgK0gAVQQQ00BbAUJoIIYagpgK0gAFcRQUwBbQQKooM5Q03q9zlBTAD0xCqiims2mdu3apT179oQOBUBBBasBmNk1ZnbEzE6Y2YNmtj9ULAAQo5BNQI9Lequ7f4ekl0j6eTO7PmA8hVamewiUKVYgZsESgLs/7O73t3//uqQTknaFiqfIOhO7lpaW5O5PTuwKfWLtdqIvaqwA1jN3Dx2DzGxU0scl7Xb3R9f8bVrStCTV6/WxhYWFRMdYXl7W8PDwgJGGMTU1pTNnzqzbXq/XtdHn0W+ZFxcXdejQIZ09e1Y7d+7Uvn37NDk52XPfgwcP6sKFC09uGxoa0tDQkB599NF1+28Wa1rK/D0nRZnjMEiZJyYmjrv7+Lo/uHvQh6RhScclvW6zfcfGxjypI0eOJH5tnlqtljcaDTczbzQa3mq13Mxc0rqHmW34Xv2UudVqea1Wu+T9a7Wat1qtrvs3Go2uMfV6bBZrWsryPaeJMsdhkDJLOuZdzqlBh4Ga2VMkvV/SvLt/IGQsRdCr+WTHjh1d909zYle/y0f0O6uYSWhA8YQcBWSS3i3phLv/Xqg4iqTXSVhS5hO7+l0+otcJ/RnPeAaT0ICSCFkDeLmkN0raa2afaj9eGTCe4HqdbM+dO5f5PQT6XT6i12zj2267jfsdACURbCKYu/+jJAt1/CIaGRnR0tJS1+3NZjPTk+js7OwlS0hLG1+5d2KZmZnRqVOnNDIyotnZ2Se3c8IHio+lIAok5Bo+Se5UlseNbQBkhwSQgrQmPoW+XSQndCAurAXUh/n5+XVNHpJSvftW1k09ANBBDWCLeg3R3L9/f6nuvsUyDQA6qAFsUa8hmmu3dRTx7lvcKxjAatQAtqgKE5+4VzCA1UgAW1SFiU/cKxjAaiSALarCxCfuFQxgNRLAFm00RLPX8Mmidbhyr2AAq9EJ3Id+hmgWscN1s9m7AOJCDSAjRe1wZbIXgA4SQEbocAVQdCSAjNDhCqDoKp8AOh2xe/fuzbUjlg5XAEVX6QQQ8gbloRd2A4DNVDoBhO6IpcMVQJFVOgHQEQsAvVU6AdARCwC9VToB0BELAL1VOgHQEQsAvVV+KYjO8g1Hjx7Vnj17QocDAIVR6RoAAKA3EgAARIoEAACRCpoAzOxGM/ucmX3ezG4NGQsAxCZYAjCz7ZL+RNIrJF0v6fVmdn2oeAAgNiFrAC+S9Hl3/4K7f0PSgqSbAsYDAFExdw9zYLMfkXSju+9rP3+jpBe7+y1r9puWNC1J9Xp9bGFhIdHxlpeXNTw8PFjQJUOZ40CZ4zBImScmJo67+/ja7SHnAViXbeuykbvPSZqTpPHxcU86lj/GeQCUOQ6UOQ5ZlDlkE9BpSdesen61pC8HigUAohMyAXxC0rVm9lwz+xZJU5LuDBgPAEQlWBOQuz9uZrdI+oik7ZJud/cHQ8UDALEJuhaQu98l6a6QMQBArJgJ3EXnPsLbtm3L9T7CAJCnyq8G2q/OfYQ7t5Ls3EdYEstIA6gUagBrhL6PMADkhQSwBvcRBhALEsAa3EcYQCxIAGtwH2EAsSABrMF9hAHEglFAXXTuIwwAVUYNAAAiRQIAgEiRAAAgUiQAAIgUCQAAIkUCAIBIkQAAIFIkAACIFAkAACJFAgCASJEAACBSJAAAiBQJAAAiRQIAgEiRAAAgUiQAAIhUkARgZr9rZp81s0+b2QfN7Okh4gCAmIWqAdwtabe7f5ek/5D0q4HiAIBoBUkA7v5Rd3+8/fRfJF0dIg4AiJm5e9gAzP5O0l+7e6vH36clTUtSvV4fW1hYSHSc5eVlDQ8PJ46zjChzHChzHAYp88TExHF3H1+7PbMEYGaLkp7V5U8z7v6h9j4zksYlvc63EMj4+LgfO3YsUTxHjx7Vnj17Er22rChzHChzHAYps5l1TQCZNQG5+6S77+7y6Jz8b5b0aknNrZz80zY/P6/R0VFt27ZNo6Ojmp+fzzsEAAjqshAHNbMbJf2KpB9w9/N5H39+fl7T09M6f37l0EtLS5qenpYkNZvNvMMBgCBCjQL6Y0lPk3S3mX3KzP4sz4PPzMw8efLvOH/+vGZmZvIMAwCCClIDcPdvD3HcjlOnTvW1HQCqKMqZwCMjI31tB4AqijIBzM7OqlarXbKtVqtpdnY2UEQAkL8oE0Cz2dTc3JwajYbMTI1GQ3Nzc3QAA4hKkD6AImg2m5zwAUQtyhoAAIAEAADRIgEAQKRIAAAQKRIAAEQq+HLQ/TCzRyQtJXz5VZK+mmI4ZUCZ40CZ4zBImRvu/sy1G0uVAAZhZse6LYdaZZQ5DpQ5DlmUmSYgAIgUCQAAIhVTApgLHUAAlDkOlDkOqZc5mj4AAMClYqoBAABWIQEAQKSiSABmdqOZfc7MPm9mt4aOJ2tmdo2ZHTGzE2b2oJntDx1THsxsu5l90sw+HDqWPJjZ083ssJl9tv1dvzR0TFkzs19s/5t+wMzeZ2ZPDR1T2szsdjM7a2YPrNq2w8zuNrOH2j+vTONYlU8AZrZd0p9IeoWk6yW93syuDxtV5h6X9FZ3/w5JL5H08xGUWZL2SzoROogc3SbpH9z9OknfrYqX3cx2SXqLpHF33y1pu6SpsFFl4i8k3bhm262SPubu10r6WPv5wCqfACS9SNLn3f0L7v4NSQuSbgocU6bc/WF3v7/9+9e1cmLYFTaqbJnZ1ZJeJelQ6FjyYGbfKun7Jb1bktz9G+7+taBB5eMySZeb2WWSapK+HDie1Ln7xyWdW7P5Jkl3tH+/Q9Jr0zhWDAlgl6QvrXp+WhU/Ga5mZqOSvlfSfYFDydofSPplSRcDx5GX50l6RNJ72s1eh8zsitBBZcnd/1PSQUmnJD0s6b/d/aNho8pN3d0fllYu8CTtTONNY0gA1mVbFGNfzWxY0vsl/YK7Pxo6nqyY2aslnXX346FjydFlkl4g6U/d/XslPaaUmgWKqt3ufZOk50p6jqQrzOwNYaMqtxgSwGlJ16x6frUqWG1cy8yeopWT/7y7fyB0PBl7uaTXmNlJrTTx7TWzVtiQMnda0ml379TsDmslIVTZpKQvuvsj7v5NSR+Q9LLAMeXljJk9W5LaP8+m8aYxJIBPSLrWzJ5rZt+ilU6jOwPHlCkzM620DZ9w998LHU/W3P1X3f1qdx/Vyvd7j7tX+srQ3b8i6Utm9vz2phsk/XvAkPJwStJLzKzW/jd+gyre8b3KnZJubv9+s6QPpfGmlb8pvLs/bma3SPqIVkYN3O7uDwYOK2svl/RGSZ8xs0+1t/2au98VLiRk4M2S5tsXNl+Q9JOB48mUu99nZocl3a+VkW6fVAWXhDCz90naI+kqMzst6W2S3iHpb8zsp7SSCH80lWOxFAQAxCmGJiAAQBckAACIFAkAACJFAgCASJEAACBSJAAAiBQJAAAiRQIABmBmLzSzT5vZU83sivZa9btDxwVsBRPBgAGZ2dslPVXS5VpZn+e3A4cEbAkJABhQeymGT0j6X0kvc/cnAocEbAlNQMDgdkgalvQ0rdQEgFKgBgAMyMzu1Moy1M+V9Gx3vyVwSMCWVH41UCBLZvYTkh539/e27z/9T2a2193vCR0bsBlqAAAQKfoAACBSJAAAiBQJAAAiRQIAgEiRAAAgUiQAAIgUCQAAIvV/oxybIWbri+oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(sim_data['x'], sim_data['y'], 'ko') # ko = black points\n",
    "plt.grid(True)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('x vs y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a Model\n",
    "Now, you are ready to build and evaluate the model using Python. Python has libraries that contain linear modeling capabilities.\n",
    "\n",
    "The first python library that is popular for linear regression is `Scikit-learn`.  Scikit-learn is a library that has many different types of machine learning algorithms, but here we concern ourselves with the linear regression capabilities.\n",
    "\n",
    "The second library is called `Stats-models`.  For those that have experience with the programming language, R, this library is the most similar because it provides easy statistical analysis of linear models that are fitted to the data.\n",
    "\n",
    "To start, we will show how to use the library `Scikit-learn` for linear regression. Model creation in Scikit-learn generally has two steps, declaring the model formulation and then fitting the model.\n",
    "\n",
    "In our case we only have one independent variable and one dependent variable. The code in the cell below does the following:  \n",
    "\n",
    "- Compute the Python model object, `linear_model`, using the simple linear regression.\n",
    "- Use the model object to compute scores (predicted values) for the dependent variable `y`. In this case, we just use the data that was originally used to compute the model. In a more general case, you can use other data to make predictions from the model.\n",
    "- The residuals of the model are computed.\n",
    "\n",
    "Execute this code and examine the head of the data frame computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Slope: [[0.9376838]]\n",
      "Intercept: [0.18626013]\n",
      "Root mean squared error:  0.7999411714439316\n",
      "R2 score:  0.9050638541982607\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhoklEQVR4nO3deZRU1dX38e+WQXBEBNGoYRBEiRPaGtQAjyLEAYdHXwccwCkYJ4gaZ0UFjaAGEaPGCYEoyCAKGhLSGF1qRKWxURlkEB6kDUITJ4LMvd8/TuNqoKvprq6qW1X391kri+6iumpXY86+95x99jF3R0RE4meHqAMQEZFoKAGIiMSUEoCISEwpAYiIxJQSgIhITNWNOoCaaNKkibdo0SLqMEREcsqMGTNWunvTrR/PqQTQokULioqKog5DRCSnmNmSyh5P+xSQmQ0zsxVmNqvCY43NrNDMFpT/uUe64xARkS1lYg1gOHDyVo/dBrzp7m2AN8u/FxGRDEp7AnD3d4Bvtnr4TGBE+dcjgLPSHYeIiGwpqiqgZu6+DKD8z70SPdHMeptZkZkVlZaWZixAEZF8l/VloO7+jLsXuHtB06bbLGKLiEiSokoAy81sH4DyP1dEFIeISGxFlQAmAb3Kv+4FTIwoDhGRrFA4Zzn9Js6icM7yjL1nJspARwPTgLZmVmJmVwADga5mtgDoWv69iEgsFc5ZTp/RxYyctoQ+o4szlgTSvhHM3Xsk+Ksu6X5vEZFc8O6CUtZs2ATAmg2beHdBKV3bNUv7+2b9IrCISL7r2KYpDevVAaBhvTp0bFOh4OX772HgQFi9OuXvm1OtIERE8lHXds0Y2qM97y4opWObpuHqf9UqePxxeOQR+PZbRn/XgCY9e6T0zkAJQERyQuGc5VsOkHmma7tm4XOtXg0PPwwPPQQrV1LauSu/bX0GM2hJw9HFDO3RPmWfX1NAIpL1kl0kjaKyJmlr1sCjj0KrVnDLLVBQAB9+yOM3DGZGk5bhKeXrA6miBCAiWa+yRdLtiaqypsbWrYM//QkOOABuvBEOPRTeew/+9jc45piq1wdqSQlARLJeMoNgMkkjo9avh6efhtat4frroU0bePttmDoVjj/+p6dtXh/oeWzzlE7/gNYARCQHVLpIuh0d2zRlXFEJazZsSvmVc61s2AAjR8KAAbBkCRx7LAwfDieeCGaV/shP6wMpZu6e8hdNl4KCAteBMCJSXVm1cLxxI4waBf37wxdfwNFHhyTQrVvCgT9VzGyGuxds/bjuAEQkb6XryrlGNm2CMWPgvvtg/nxo3x4mTYLu3dM+8G+P1gBERNKhrAzGjg2LuhddBA0awKuvwowZcPrpkQ/+oAQgIpJaZWUwYQIcfjicf34Y6MeOheJiOOusrBj4N1MCEBFJBXd4/XU46ig455xQ5fPSS/Dpp3DuubBDGG6zaW+CEoCISG24/1SzzxlnwA8/wIgRMHs2XHgh1Knz01OzbW+CEoCISDLcobAw1OyfeiqUlsLzz8Pnn0PPnlB32xqbbNuboAQgIjkv49Mqb78NnTuHEs6lS+Gpp0KFz+WXQ716CX8snbt6k6EyUBHJaZunVdZs2MS4opKU75bdwnvvQb9+8NZb8LOfhRYOV14JO+5YrR9PZkNbOikBiEhOy8RhKh+Nnkzjhx6g9cz3oVkzGDIEeveGhg1r/FpZsTehnKaARCQSqZq2Seu0SlERpZ1P4pgLT2OPebMYdNKVvPnXadC3b1KDf7bRHYCIZFwqp23SMq0ycybccw9MmsTOu+7OwM6XMvLI0/ixfkN6lqymy1G1f4tsoAQgIhmX6mmblE2rzJoF994Lr7wCjRrBgAF80O0CRry+MPuayqWAEoCIZFzWder8/PMw8I8dC7vsEhZ6b7gBGjXiRGDoLrtmzcJtKqkbqIhEIis6dS5YELpzjhoV5vT79oWbboLGjWv90lnx+col6gaqBCAi8bNoEdx/f+jLX78+XHcd3HwzNE3NnUjFNY6G9eqktzS1GhIlAFUBiUh8LFkSyjfbtoXRo6FPH1i8OBzAnqLBH7Jvx28ikSYAM7vBzGab2SwzG21mDaKMR0Sil5ZdvSUlcM014djFESPg6qvDoSyDB4e6/hTLth2/iUQ2BWRm+wLvAe3cfY2ZjQUmu/vwRD+jKSCR/JbyqZNly2DgwHD2blkZXHEF3HEH7L9/6oJOIBfWAKKuAqoLNDSzDcBOwL8jjkdEIpSy8tAVK2DQIHjyyXAG72WXwZ13QosWqQ24Ctm04zeRyKaA3P0r4BHgS2AZ8L27/2Pr55lZbzMrMrOi0tLsnEcTkcRqMqVT66mTlSvhttugZcvQruH882HePHj22YwO/rkiyimgPYBXgPOB74BxwHh3fzHRz2gKSCS3JDOlk9TUybffhvn8IUNg9Wro0SPs5D3wwNp/iDyQjVNAJwGL3b0UwMwmAMcBCROAiOSWZKZ0ajR18v33YdAfPDgcxHLeeWHgb9eulpHHQ5RVQF8CHcxsJzMzoAswN8J4RCTF0lYNs2oVPPBAmNa5917o0gU++QTGjNHgXwOR3QG4+4dmNh74GNgIFAPPRBWPiKReyhu1rV4NTzwR6vb/8x84/XS47z5o3z41AceMdgKLSPZbswb+/OdQ0rliBZxyShj4jz466shygnYCi0juWbsWHn8cWrWCG2+EQw+Ff/0LJk/W4J8CUe8DEBHZ1rp1MGxYmOf/6qtw/u6YMdCpU9SR5RUlABHJHhs2wPDhoVHbl1/CcceFhm0nnABmKXubbNqlGyVNAYlI9DZuhBdeCE3aeveGvfeGKVPCIewnnpjywb/P6GJGTltCn9HFqe05lGOUAEQkOps2wYsvwsEHw+WXwx57wBtvwAcfQLduKR34N0t1p860NK/LECUAEcm8sjJ4+WU45BC45BLYaSd47TUoKoLTTkvLwL9ZKvcm5PrdhNYARCRzyspgwoSweWv27LBpa9w4OPts2CEz16NV7U2o6dpAqs82zjQlABFJiSoHT3eYNCm0afjkEzjooHAHcO65GRv4K6qs3UTFvkXjikqq1bco6842riElABGptYSDpzv87W/hkPUZM6B1a/jLX0Kztjp1og57C8n2LUrpTucMUwIQkVrbZvCcv4KuJZ+Egf/DD0N75hdegIsvhrrZOewkezWfC33/E8nOfwkRySkVB8/OX83ihgED4OMP4ec/D734e/WCevWiDrNKuX41nwwlABGpta7tmjHywHU0eegBWs6aDvvuG07juvxy2HHHqMOrtly+mk+GEoCI1M60adCvH0dPnRo2cA0dCr/5DTRoEHVksh3aByAiyZk+HU49NbRr+OQTeOQR+OILuP56Df45QncAIlIzxcWhnPP112HPPcPh69deCzvvHHVkUkNKACJSPZ99FjZwTZgAjRqFhm19+sCuu0YdmSRJCUAkBmrV/XLu3DDwjx0Lu+0Wrv5vuAF23z0tsUrmKAGI5LlkdrgCMH8+9O8Po0aF6Z077wyHsjRunP6gJSO0CCyS52rc/XLRIrjsstCh89VX4eabYfHiMOWjwT+vKAGI5Llqd79csiSUbx54YOjT07dvSAaDBkGTJkButz6WbWkKSCTPbXeHa0lJOHrx+edDG+Zrr4XbboN99tniaUlPJUnWUgIQyTHJLOhWusN12TJ48EF4+unQtO3KK+GOO2C//cJ7fDRri/fI9dbHsi1NAYnkkJQcQLJ8eVjMbdUqtGvo2RMWLAhflw/+lb1HKg9SkeygBCCSQ2p1nOHKlXDLLaEz52OPwQUXhEqfZ5+F5s23+x6bp5J6Httc0z95ItIEYGaNzGy8mX1uZnPN7Ngo4xHJdkldhX/zTSjhbNkytGs4++xQ2//CCxSu3XmbRd2q3qNru2b0P/MQDf55wtw9ujc3GwG86+7PmVl9YCd3/y7R8wsKCryoqChj8Ylko2qvAXz3HTz6aPjfqlVw3nlhQ9fBB//0OpsXdRvWq7PFVX2tNo5J1jGzGe5esPXjkS0Cm9luQCfgUgB3Xw+sjyoekVyx3ZbFP/wQOnL+8Y8hCZxzTti9e+ihWzytqkXduLVFjqsop4BaAaXAC2ZWbGbPmdk23aTMrLeZFZlZUWlpDeY7ReLmv/+FgQPDVM/dd0OnTqFx2/jx2wz+oEVdiXAKyMwKgA+A4939QzN7DPjB3e9O9DOaAhKpxI8/hgqeQYPCQu+pp8J990HBNnf829BUTzxk3RQQUAKUuPuH5d+PB26LMB6R3LJmTajhHzgwlHZ26xYG/g4dqv0SmuqJt8imgNz9a2CpmbUtf6gLMCeqeERyxrp18MQT0Lp16MrZrh288w5MmULhbi3VqkGqLeqdwNcDL5VXAC0CLos4HpHstX49DB8emrItXQq/+hW8+CKccAKgVg1Sc5HuA3D3me5e4O6HuftZ7v5tlPGIZKUNG2DYMGjbFq66Khy4PmVKuOovH/yhlpvEJJa0E1gkW23cCCNHhrr9K64IHTknT4b33w/z/WZbPF1VPVJTUU8BicjWNm0Kp2/ddx/MmwdHHAGTJkH37tsM+hVtt+unyFaUAESyRVkZvPJK2K07Z06o3X/lFTjrLNihejfrquqRmtAUkEjU3OG116B9+9CuwR3GjIGZM0PfnmoO/iI1pf+yRKLiDn/9a9iw9b//G+r6X3oJPvssJAIN/JJm+i9MJNPcQxVPhw7QvTs/lv6HCdcPYOr4t+DCC6FOnagjlJjQGoBIprjDP/8J/fqFSp7mzZndfzAXrGnDqjKj4bjPGFq3rubwJWN0ByCSCZtr9k86KRy+/tRTMH8+Yw7ryqqyUNmj2n3JNCUAkXR6//0w6HfuHEo6hw6FhQvht7+F+vVVuy+R0hSQSDp89FHowf/3v8Nee8HgwWHQb9hwi6epdl+ipAQgkkrFxWGO/403YM894aGH4JprYOdtjrr4iWr3JSqaAhJJhU8/DTX7Rx7Jj2+/w4K+t8PixXDzzVUO/iJRUgIQqY3Zs0PN/uGHs6FwKo93uphfXvksZ+zaicKlP0YdnUiVNAUkkox580KvnpdfDlf4d93FwwedzDOffRf+fqszdkWyke4AJGcVzlle48NPEv1MtV9r4ULo1SscwjJxItxyS5jqGTCAo9u3VkWP5JTIzgROhs4Els0qHn7SsF6dah1+kuhnqvVaixeHg1hGjIB69eDaa8Pgv9de27yHKnok2yQ6E1h3AJKTkjn8JNHPVPlaS5eG8s0DDwx9eq67DhYtgkce2Wbwh1DR0//MQzT4S05QApCclMwGqkQ/U+njX30VBvvWrcNpXFddBV98AUOGwD77pOdDiWSYpoAkZyUz3ZLoZzY/3qWR0/nV5+HPfw4Hs1x+Odx5J/z85+n6GCJpl2gKSAlABKC0FAYNgiefDIev9+oFd90FLVum7C20PiBRSZQAVAYqeadGA+1//hPm8x9/PPTjv+giuPtuaNMm5TFtXmgeV1RSrUVrkXTTGoDklc0D7chpS+gzujhxWee334aBvmXLcOV/+ulhU9fIkSkf/CG5RWuRdFMCkLyy3YH2+++hf/8w8N9/P/z616GNw+jRcNBBaYsrmUXrZPY5iNSEpoAkr3Rs05RxRSU/1fT/NNCuWgV/+hM8/HC4+j/rrHD4+uGHZySumnb91JSRZELkCcDM6gBFwFfu3j3qeKRyubKAuc1A23yXMOg/9BCsXAmnnhruAI46qlbvk8zvoyZdPyu7k8nm37vkpmyYAuoLzI06CEms2vPqWaJru2b073YAXaeMglatwo7do46CDz4Ih7CnYPBP9+9DB8VIJkSaAMxsP+A04Lko45CqZWoBMyVz3uvWhameAw6AG2+EQw+F994LB7P88pcpiTMTv4/NdzI9j22u6R9Jm6jvAIYAtwBliZ5gZr3NrMjMikpLVTkRhUxcjdb6qnr9enj66bBz9/rrQyXP22/D1Klw/PEpjTVTV+dqKyHpFtkagJl1B1a4+wwz+59Ez3P3Z4BnIGwEy0x0UlEmji1Mes57w4ZQujlgQDhs/dhjYfhwOPFEMEt5nKBjHCV/RLkIfDxwhpmdCjQAdjOzF9394ghjkgTSfWxhwuqdRDZuhFGjQk/+RYvgmGPCHUC3bmkb+CvSMY6SDyJLAO5+O3A7QPkdwO81+CcnVyp0qlLtq+pNm2DMmDDwz58P7dvD66/DaaelZeDPh9+tSCLbTQBmdh3wkrt/m4F4pIbyqV68yqvqsjIYPz7U7s+dC4cdBq++CmeembYr/nz63YpUpjqLwHsD081srJmdbJb6/7e5+9vaA5CcvG8x4B4G+iOOgPPPD4P92LFQXBw2c6Vxuifvf7cSe9tNAO5+F9AGeB64FFhgZn8wswPSHJtUQ97Wi7uHqZ2jjoKzzw7lnaNGhbYN554LO6S/gC1vf7ci5aq1BuDubmZfA18DG4E9gPFmVujut6QzQKlaLlakVDmv7g5TpkC/fjB9etjINWIEXHgh1M3sklUu/m5FamK75wGYWR+gF7CSsGHrNXffYGY7AAvcPWN3AjoPIPclPH/XHd58Mwz806ZB8+bh60suCWfwikjSanMeQBPgbHdfUvFBdy8rr+UXqbZK6/1XzA2D/bvvwn77hdO4LrsM6tcHVIkjki7VWQPot/XgX+Hv1MNHaqTivPpxX39O3wevhhNOgIULQwuHhQvD+bsVBv9c6kMkkksi7wYq8dK1XTOG/8LZc9AfaD3zfdhrL3j00TDoN2y4zfPVFVMkfaLuBSRxMmMGdO/OL3ucSuuS+aFN86JF8LvfVTr4gypxRNJJdwCSfp98AvfcAxMnQuPG8Ic/hIZtu+yy3R9VJY5I+igBSPrMnh127o4fD7vvHg5i6dsXdtutRi+jvjsi6aEEIKn3+eehV8+YMeEqv18/uOEGaNQo6shEpAIlgIjkZWnjggXhKn/UqDCnf+ut8Pvfw557Rh2ZiFRCCSACeddkbNEiuP/+0Je/fn246Sa4+WZoqgVbkWymKqAI5E2TsSVLoHdvaNs2XPVff31IBg89pMFfJAcoAUQg50sbS0rgmmvCsYsjRoQa/kWLQj3/3ntHHZ2IVJOmgCKQs6WNy5bBwIHh5K1Nm+CKK+COO+DnP486MhFJghJARHKqtHHFChg0CJ58MpzBe+mlcNdd0KJF1JGJSC0oAUhiK1fCI4/A44/D2rVw0UWhpLN166gjE5EUUALIY0mXmn7zDQweDI89BqtXQ48eYeBv2zZ9wYpIxikB5KmkSk2//x6GDAmD/w8/wHnnhRYO7dplJGYRySxVAeWpGpWarloFDzwQ5vTvvRe6dAn9e8aM0eAvkseUAPJUtUpNV68ONfstW4ZF3Y4d4eOPYcIEOOywDEcsIpmmKaA8VWWp6Zo18NRTobJnxQo45ZTQu+foo6MLWEQyTgkgj21Tarp2LTz7bGjH/PXXcNJJYeA/7rjoghSRyCgBxMG6dTBsWJjn/+or6Nw5zO936hR1ZCISocjWAMxsfzN7y8zmmtlsM+sbVSx5a8MGeO45OPDA0LqheXN480146y0N/iIS6SLwRuAmdz8Y6ABca2YqOUmFjRth+PBQt/+b38A++8CUKfDee3DiiWAWdYQikgUiSwDuvszdPy7/ehUwF9g3qnhyQeGc5fSbOIvCOcsrf8KmTfDii3DwwXDZZbDHHvDGGzBtGnTrltGBf7uxikjksqIM1MxaAO2BDyv5u95mVmRmRaWlOdo2OQU2b+waOW0JfUYXbzmwlpWFOf1DDoFLLoGddoLXXoOiIjjttLQO/JUN9FXGKiJZI/IEYGa7AK8Av3P3H7b+e3d/xt0L3L2gaYx7zFe6sausLNTsH344XHAB1KkTzt8tLoYzz0xq4K/JlXuigT5vzjsQyXORJgAzq0cY/F9y9wlRxpJNKhuEt9jYVXcHzl76MRx1FJxzTljsHT0aPv00fL9Dcv+sNb1yTzTQ5/x5ByIxEVkZqJkZ8Dww190HRxVHtknUw6dru2YMveAIvh7zKmdMeJrdZ38SunL+5S+hWVudOrV+78oG9Kr6B3Vs05RxRSWs2bBpi4E+Z887EImZKPcBHA9cAnxmZjPLH7vD3SdHF1L0Kh2ED94Lpk6la79+8MEHoWfPsGFhvr9u6v4JEw3oiVQ10OfUeQciMRVZAnD39wDVI25l60H4jG/mQaerQwnn/vuH07guvTQcvp5iyVy5a6AXyV3m7lHHUG0FBQVeVFQUdRhpVzhnOUsmTuHs156h8Uf/gp/9DO68MxzBuOOOUYcnIjnGzGa4e8HWj6sVRLb54IMw1VNYCM2ahUNZeveGBg2ijkxE8owSQIokffrWZkVF4fCVyZOhSZNwFOPVV4eafhGRNFACqKHKBvqkTt/arLg4HMIyaRI0bgwDB8K118Iuu6TvQ4iIkAUbwXJJSjc+ffZZqNk/8kh45x0YMAAWL4Zbb0374K82DSICSgA1kpKNT3PnwvnnhxO3CgvDtM/ixeFErt12S/tnUJsGEdlMCaAGEg30m8snex7bPPH0z/z5cPHF8ItfwF//CnfcAf/3f2H6p1GjjH0GtWkQkc20BlADSW18WrQI+vcPO3YbNICbb4bf/x4i6mtU081eIpK/tA8gXZYsgfvvhxdegHr1QkXPrbeG0s6I1bpiSURyivYBZEpJSTh68fnnwYwvL7iUMV0u5Ihf/oKuWTD4g3bvikigNYBUWbYM+vSBAw4Ig/8VV/DO5Gn8uuX/44l5a7TgKiJZRwmgtpYvhxtvhFat4MknoWfPsOD71FNM/W99LbiKSNZSAkjWypVwyy3QsmVo13D++TBvHjz7bOjWifrii0h20xpATX3zDfzxjzB0KKxeDRdeCP36wYEHbvNU9cUXkWwWiwSQkqqX776DRx+FIUPghx/CFf8994QD2KugBVcRyVZ5nwBq1acHwmA/dGi46v/uu9C+4Z574NBD0xaziEgm5P0aQNI7X//739CYrWVLuPtu6NQpNG4bP16Dv4jkhbxPADVeiP3xx3C136oV3H47dOgA06fDxIlwxBHpD1hEJEPyfgqo2guxa9eG4xYffDCUdnbtGlo4dOiQ2YBFRDIk7xMAbGchdt26sHHrgQfg3/+GE06AceOgY8fMBikikmGxSACVWr8ehg8PA/+XX8KvfgUvvhgSgIhIDOT9GsA2Nm4MDdratoWrrgoHrv/jH+FQFg3+IhIj8UkAmzaFlswHHwyXXx7O3Z08Gd5/P8z3m0UdoYhIRsUjAbz6ajiIpWdP2HnnUNHz0Udwyika+EUktiJNAGZ2spnNM7OFZnZb2t5o5szQk/+VV+Djj+GMMzTwi0jsRXYgjJnVAeYDXYESYDrQw93nJPqZpA+EWbsW6teHHeJxwyMiUlGiA2GiHBGPARa6+yJ3Xw+8DJyZlndq0ECDv4jIVqIcFfcFllb4vqT8sS2YWW8zKzKzotJS9dMXEUmVKBNAZZPw28xHufsz7l7g7gVNIzpIXUQkH0WZAEqA/St8vx/w74hiERGJnSgTwHSgjZm1NLP6wAXApAjjERGJlchaQbj7RjO7DpgC1AGGufvsqOIREYmbSHsBuftkYHKUMYiIxJVqI0VEYiq+3UC3IyXnCIuIZDHdAVRi8znCI6ctoc/oYgrnLI86JBGRlFMCqETS5wiLiOQQJYBK1PgcYRGRHKQ1gEpU+xxhEZEcpgSQQJXnCIuI5AFNAYmIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMRUJAnAzB42s8/N7FMze9XMGkURh4hInEV1B1AIHOLuhwHzgdsjikNEJLYiSQDu/g9331j+7QfAflHEISISZ9mwBnA58LdEf2lmvc2syMyKSktLMxiWiEh+q5uuFzazqcDelfzVne4+sfw5dwIbgZcSvY67PwM8A1BQUOBpCFVEJJbSlgDc/aSq/t7MegHdgS7uroFdRCTD0pYAqmJmJwO3Ap3d/ccoYgAonLOcdxeU0rFNU7q2axZVGCIikYhqDeBPwK5AoZnNNLM/ZzqAwjnL6TO6mJHTltBndDGFc5ZnOgQRkUhFcgfg7q2jeN+K3l1QypoNmwBYs2ET7y4o1V2AiMRKNlQBRaJjm6Y0rFcHgIb16tCxTdOIIxIRyaxI7gCyQdd2zRjao73WAEQktmKbACAkAQ38IhJXsZ0CEhGJOyUAEZGYUgIQEYkpJQARkZhSAhARiSklABGRmLJc6sNmZqXAkiR/vAmwMoXh5AJ95njQZ46H2nzm5u6+zW7XnEoAtWFmRe5eEHUcmaTPHA/6zPGQjs+sKSARkZhSAhARiak4JYBnog4gAvrM8aDPHA8p/8yxWQMQEZEtxekOQEREKlACEBGJqVgkADM72czmmdlCM7st6njSzcz2N7O3zGyumc02s75Rx5QJZlbHzIrN7I2oY8kEM2tkZuPN7PPyf+tjo44p3czshvL/pmeZ2WgzaxB1TKlmZsPMbIWZzarwWGMzKzSzBeV/7pGK98r7BGBmdYAngFOAdkAPM2sXbVRptxG4yd0PBjoA18bgMwP0BeZGHUQGPQb83d0PAg4nzz+7me0L9AEK3P0QoA5wQbRRpcVw4OStHrsNeNPd2wBvln9fa3mfAIBjgIXuvsjd1wMvA2dGHFNaufsyd/+4/OtVhIFh32ijSi8z2w84DXgu6lgywcx2AzoBzwO4+3p3/y7SoDKjLtDQzOoCOwH/jjielHP3d4Bvtnr4TGBE+dcjgLNS8V5xSAD7AksrfF9Cng+GFZlZC6A98GHEoaTbEOAWoCziODKlFVAKvFA+7fWcme0cdVDp5O5fAY8AXwLLgO/d/R/RRpUxzdx9GYQLPGCvVLxoHBKAVfJYLGpfzWwX4BXgd+7+Q9TxpIuZdQdWuPuMqGPJoLrAkcBT7t4eWE2KpgWyVfm895lAS+BnwM5mdnG0UeW2OCSAEmD/Ct/vRx7eNm7NzOoRBv+X3H1C1PGk2fHAGWb2f4QpvhPN7MVoQ0q7EqDE3Tff2Y0nJIR8dhKw2N1L3X0DMAE4LuKYMmW5me0DUP7nilS8aBwSwHSgjZm1NLP6hEWjSRHHlFZmZoS54bnuPjjqeNLN3W939/3cvQXh3/ef7p7XV4bu/jWw1Mzalj/UBZgTYUiZ8CXQwcx2Kv9vvAt5vvBdwSSgV/nXvYCJqXjRuql4kWzm7hvN7DpgCqFqYJi7z444rHQ7HrgE+MzMZpY/doe7T44uJEmD64GXyi9sFgGXRRxPWrn7h2Y2HviYUOlWTB62hDCz0cD/AE3MrAS4BxgIjDWzKwiJ8NyUvJdaQYiIxFMcpoBERKQSSgAiIjGlBCAiElNKACIiMaUEICISU0oAIiIxpQQgIhJTSgAitWBmR5vZp2bWwMx2Lu9Vf0jUcYlUhzaCidSSmd0PNAAaEvrzPBhxSCLVogQgUkvlrRimA2uB49x9U8QhiVSLpoBEaq8xsAuwK+FOQCQn6A5ApJbMbBKhDXVLYB93vy7ikESqJe+7gYqkk5n1BDa6+6jy86ffN7MT3f2fUccmsj26AxARiSmtAYiIxJQSgIhITCkBiIjElBKAiEhMKQGIiMSUEoCISEwpAYiIxNT/B9nNK3quoCVOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# standard template\n",
    "\n",
    "# imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "# Create numpy arrays for the x (independent) and y (dependent) variables\n",
    "x = sim_data.x.to_numpy().reshape(-1,1)\n",
    "y = sim_data.y.to_numpy().reshape(-1,1)\n",
    "\n",
    "# sckit-learn implementation\n",
    "\n",
    "# Model initialization\n",
    "regression_model = LinearRegression()\n",
    "# Fit the data(train the model)\n",
    "regression_model.fit(x, y)\n",
    "# Predict\n",
    "y_predicted = regression_model.predict(x)\n",
    "\n",
    "# model evaluation\n",
    "rmse = mean_squared_error(y, y_predicted)\n",
    "r2 = r2_score(y, y_predicted)\n",
    "\n",
    "# printing values\n",
    "print('Slope:' ,regression_model.coef_)\n",
    "print('Intercept:', regression_model.intercept_)\n",
    "print('Root mean squared error: ', rmse)\n",
    "print('R2 score: ', r2)\n",
    "\n",
    "# plotting values\n",
    "\n",
    "# data points\n",
    "plt.scatter(x, y, s=10)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "# predicted values\n",
    "plt.plot(x, y_predicted, color='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/anajaved/opt/miniconda3/lib/python3.8/site-packages/seaborn/_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD6CAYAAAC8sMwIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUdElEQVR4nO3df2xV533H8c+3xhScssAWvGXARtiysaTqlMqqkiViqMnUdKuSblqlMq2KtlWHKe1I0cQW1iUdSGhDWyo2UQ0QyahYRFXSoMLGFpJSZEUoqHaabHGBFjlQnIXaqJBAjOUf/e4PG9cGX+x7z8/nnPdLisw5ju99rn3u5z7ne57nOebuAgCE6315NwAAEA9BDgCBI8gBIHAEOQAEjiAHgMAR5AAQuNhBbmZLzOzbZnbczLrM7LEkGgYAmBmLO47czG6VdKu7v2pm8yR1Svqku3+v1s/ccsstvnTp0ljPCwBV09nZed7dF167f1bcB3b3tyW9PfbvS2Z2XNIiSTWDfOnSpero6Ij71ABQKWZ2Zqr9idbIzWyppLskHZvie5GZdZhZR19fX5JPCwCVlliQm9kHJH1D0hfc/d1rv+/uO9y9zd3bFi687swAANCgRILczJo1GuLPuvvzSTwmAGBmkhi1YpKelnTc3b8cv0kAgHok0SO/V9JnJH3UzF4b++93E3hcAMAMJDFq5WVJlkBbKu3IiV5tb+/W2Qv9WrKgRatXLNPK5a15NwtAAJjZWQBHTvTqyf1d6r00oPlzm9V7aUBP7u/SkRO9eTcNQAAI8gLY3t6t5iZTy+xZMhv92txk2t7enXfTAAQgdmkF8Z290K/5c5sn7Zvb3KSeC/05tQhIByXEdNAjL4AlC1p0ZWhk0r4rQyNavKAlpxYByaOEmB6CvABWr1imoRFX/+Cw3Ee/Do24Vq9YlnfTgMRQQkwPQV4AK5e3auNDd6p13hy9c2VIrfPmaONDd3LKiVI5e6Ffc5ubJu2jhJgMauQFsXJ5K8GNUluyoEW9lwbUMvunsUMJMRn0yAFkghJieghyAJmghJgeSisAMkMJMR30yAEgcPTIY2ByA4AioEfeICY3ACgKgrxBTG4AUBQEeYOY3ACgKAjyBrE+CoCiIMgbxOQGAEVBkDeIyQ0AioLhhzGENrmB4ZJAOdEjrwiGSwLlRZBXBMMlgfKitFIR3E4OoaIkOD165BXBcEmEiJLgzBDkFcFwSYSIkuDMEOQVwXBJhIgZ1DNDjbxCQhsuCXB7uJmhRw6gsCgJzgw9cgANS3tEycrlrdqo0Vp5z4V+LWbUypQIcgANuTqipLnJJo0o2SglHuYE941RWgHQEEaUFAdBDqAhjCgpDoIcQEOYZFYcBDmAhjCipDgIcgANYZJZcTBqBUDDGFFSDIn0yM3sGTPrNbM3kng8AMDMJVVa2SXpwYQeCwBQh0SC3N3bJf04iccCANQns4udZhaZWYeZdfT19WX1tABQepkFubvvcPc2d29buHBhVk8LAKXH8EMACBxBDgCBS2QcuZntkbRS0i1m1iPpS+7+dBKPjXLjxrpAfIkEubuvSuJxUC1ZLYOKMPCh3jhmdiI3E5dBlaSW2bPUPzis7e3dvIELJu2Q5UM9HmrkyA3LoIbhasj2XhqYFLJHTvQm9hysbR4PQY7csAxqGLIIWT7U4yHIkRuWQQ1DFiHLh3o8BDlywzKoYcgiZPlQj4eLncgVy6AW3+oVy/Tk/i71Dw5rbnOTrgyNJB6yK5e3aqNGyzg9F/q1mFErdSHIAdxQViHLh3rjCHIA0yJki40aOQAEjiAHgMBVsrTCVGAgXbzHslW5HnkWs9SAKuM9lr3KBTlTgYF08R7LXuVKK2cv9Gv+3OZJ+5gKnBxOqWemzL8n3mPZC6ZHfuREr1bteEX3bT6sVTteafg0janA6eGUembK/nviPZa9III8yQOfqcDpyeqUOqkP9byUvfTAeyx7QQR5kgc+63ukJ4vFlcrQmy37Sn+8x7IXRI086ZpbrVlqRaxbFrFNtSxZ0KLeSwPjN4qQkj+lLsPNKBr9PYV0LDATNFtB9MizqLkVsadXxDbdSBan1GXozTbyewrtWEC2ggjyLAKiiHXLIrbpRrI4pS7DhbRGfk+hHQvIVhCllSxWXyvikKkitmk6aZ9SZ7Gkahbq/T2FeCwgO0EEuZR+QGRR3y1Dm/JW1XWrORZwI8EEedqK2NMrYpuycqMLe1W8kFblYyEkeV2QNndP/Umu1dbW5h0dHZk/73Su/hGK1NMrYpvSdvXCXnOTTQqt6erIIY3qaEQVj4WQNHrc1sPMOt297br9BDmKZtWOV64rI/QPDqt13hztie6e8meyeBMBN9LIcVuvWkEexKgVVEsjQwwZ1YG85Tk0lhp5Smqd5pfl9D/N19HIhT1GdSSjLMdnHvK8IE2PPAW1Jm/8y0vfL8WkjrQnpzQyb6AM48vzxqSjePJcY4YgT0Gt0/ydL79ZitP/tMsYjUyYYaGm+ChPxZPnGjOUVlJQ6zT/vcER/VLg08ulbMoY9Q4xrOr48iRRnoovr6GxBHkKatXKbpo9Opoi9EkdRZ2cUsXx5Ukq6t8V06O0koJap/mfve+2Upz+U8YoJ/6u4aJHnoIbneZ/aPH84E//KWOUE3/XcDEhCAACwYQgACipRILczB40s5NmdsrMHk/iMQEAMxM7yM2sSdJXJH1c0h2SVpnZHXEfFwAwM0n0yD8i6ZS7d7v7oKSvSXo4gccFAMxAEkG+SNLZCds9Y/tqOnPmjA4cOCBJGh4eVhRFOnjwoCRpYGBAURTp0KFDkqTLly8riiIdPnxYknTx4kVFUaT29nZJ0vnz5xVFkY4ePSpJOnfunKIo0rFjx0Yb09OjKIrU2dkpSTp9+rSiKNLrr78uSTp16pSiKFJXV5ck6eTJk4qiSCdPnpQkdXV1KYoinTp1SpL0+uuvK4oinT59WpLU2dmpKIrU09MjSTp27JiiKNK5c+ckSUePHlUURTp//rwkqb29XVEU6eLFi5Kkw4cPK4oiXb58WZJ06NAhRVGkgYEBSdLBgwcVRZGGh4clSQcOHFAUReO/y3379unRRx8d3967d6/WrFkzvr1nzx6tXbt2fHv37t1at27d+PauXbu0fv368e2dO3fqiSeeGN/etm2bNmzYML69detWbdq0aXx7y5Yt2rx58/j2U089paeeemp8e/PmzdqyZcv49qZNm7R169bx7Q0bNmjbtm3j20888YR27tw5vr1+/Xrt2rVrfHvdunXavXv3+PbatWu1Z8+e8e01a9Zo796949uPPvqo9u3bN74dRRHHXkWPvc+t+1v9zmfX677Nh7Vqxyv687WPB3fs1ZLE8EObYt91Q2HMLJIUSdLNN9+cwNMCwMwcOdGrjjMXNOv9N2nB2Doy3Wcvav7PvZd30xIRe/ihmd0j6e/c/WNj2+slyd3/vtbPMPwQRcYKgOWTxVrh00niuKo1/DCJHvl3JN1uZrdJekvSpyX9UQKPC2Ru4g0qJq4AuFGq9N2JQpf3OjKNHlczFbtG7u7Dkj4v6QVJxyV93d274j4ukIdGVgBk+dfiy3uZ47RXlkxkHLm7H3T3X3P3X3H3TdP/RFiOnOjVqh2vjF8k4Q1aXtydqJzyXkcm7bsHMbNzGvS2qqWRnluet/jCzOS5VriU/hkBi2ZNY2JvS5JaZs9S/+Cwtrd3UwMtodUrlunJ/V3qHxyedBPn6e5OlOTyr9Tb05HnMseNHFf1IMinkfdFEmSrkRUAk3yTNnJRrOz3hy2DtFeWZPXDaRRh2FIZlD1Urr6+uG/Seo+3icE/8UPkDz+8SM+9+tZ1+7MsJyB5aQ4/LLW0T4mqIO2hV0WQ1Gl7vWeAtUp/O19+UwvnvX/KkuDVnyvrh2oVcbFzGnlfJCkDRnXMXL0XxWpdaH1vcGTK/T/40btcvC8heuQzwL0g4+E6w8zVewZY7/1hB0dcN3PxvnTokSN1eU/GCEm9Z4D13h929qz3MVQyAUWbW0KPHKnjOkN96jkDrPf+sNvbuxMdKllFRbzmQ5AjddzUN121gr/Wfj5U4yni3BKCHJngOkMx8KEaXxGv+RDkQMXwoRpP0jN5k8DFTgCoQ94LcE2FIAeAOhRxbgmlFQCoU9HKU/TIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOAYflhCZb8bD4DJ6JGXzNWV2bhxAFAdBHnJcDceoHoorZRMEVdmKypKUCgLgrxkirgyWxEV8eYAKJ5QPuwprZRMEVdmK6JGS1BFu8UX0hPS9SaCvGSKuDJbEdW6+/yNSlAhvbGrIO0P1ZCuN1FaKaGircxWRI2UoIp4i6+qyqI0FtL1JnrkqKRGSlCN9OKRjix6y0sWtOjK0MikfUW93kSQo5IaKUGF9MYuuyw+VEO63kRpBaVXa+RBvSWo1SuWcQf6gshidFZIN6o2d8/8Sdva2ryjoyPz50X1TKylTgzfRi8AX/1QKPobu+yS/ruGwsw63b3t2v30yFFqSV+g5EJyMYTUW84CQY5SC2nkAerDh+pPEeSBSnLGWSiz1xrBTFdUQaxRK2b2KTPrMrOfmNl1dRukI8mJKWWf5BLSyAOgUXGHH74h6Q8ktSfQFsxQkmNoQ5q91ghmuqIKYpVW3P24JJlZMq3BjCRZ961CDZlaKsouswlBZhaZWYeZdfT19WX1tKWU5MQUJrkA4Zs2yM3sJTN7Y4r/Hq7nidx9h7u3uXvbwoULG29xHcq6Ul2SdV9qyED4pi2tuPsDWTQkaWVebzrJMbQ3eqwyj2YByqS0ww/LvlJdknXfqR6rzB+EQNnEHX74+2bWI+keSf9pZi8k06z4WKkunrKPZgHKJO6olX2S9iXUlkQxESSeKoxmAcqitMvYchEvHkazAOEobZAzESQePggRqrKOVrsRlrFFTSzZitCUfXlblrFF3ZgRidCUfbRaLaUtrQConqqOViPIAZRGVS/SE+QASqOqF+mpkQMpYpmDbFX1FnAEOZASljnIRxUv0lNaAVLCMgfICkEOpKSqIyiQPYIcSElVR1AgewQ5kJKqjqBA9ghyICWs94OsMGoFSFEVR1CEKPRhovTIAVTa1WGivZcGJg0TDWnVRIIcQKWVYZgoQQ6g0sowTJQgB1BpZRgmSpADqLQyDBMlyAFUWhmGiTL8EEDlhT5MlB45AASOIAeAwBHkABA4ghwAAkeQA0DgCHIACBxBDgCBI8gBIHAEOQAEjiAHgMAR5AAQOIIcAAJHkANA4GIFuZn9o5mdMLP/MbN9ZjY/oXYBAGYobo/8RUkfdPcPSfq+pPXxmwQAqEesIHf3Q+4+PLb5iqTF8ZsEAKhHkjXyP5X0X7W+aWaRmXWYWUdfX1+CTwsA1TbtHYLM7CVJvzDFt77o7t8c+3++KGlY0rO1Hsfdd0jaIUltbW3eUGsBANeZNsjd/YEbfd/MHpH0CUn3uzsBDQAZi3XPTjN7UNJfS/ptd+9PpkkAgHrErZFvlTRP0otm9pqZbUugTQCAOsTqkbv7rybVEABAY5jZCQCBI8gBIHAEOQAEjiAHgMAR5AAQOIIcAAJHkANA4AhyAAgcQQ4AgSPIASBwBDkABI4gB4DAEeQAEDiCHAACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AASOIAeAwMW6+XIRHDnRq+3t3Tp7oV9LFrRo9YplWrm8Ne9mAUBmgu6RHznRqyf3d6n30oDmz21W76UBPbm/S0dO9ObdNADITNBBvr29W81NppbZs2Q2+rW5ybS9vTvvpgFAZoIO8rMX+jW3uWnSvrnNTeq50J9TiwAge0EH+ZIFLboyNDJp35WhES1e0JJTiwAge0EH+eoVyzQ04uofHJb76NehEdfqFcvybhoAZCboIF+5vFUbH7pTrfPm6J0rQ2qdN0cbH7qTUSsAKiX44Ycrl7cS3AAqLegeOQCAIAeA4BHkABA4ghwAAkeQA0DgzN2zf1KzPklnGvzxWySdT7A5IeA1VwOvuRrivOZfdveF1+7MJcjjMLMOd2/Lux1Z4jVXA6+5GtJ4zZRWACBwBDkABC7EIN+RdwNywGuuBl5zNST+moOrkQMAJguxRw4AmIAgB4DABRXkZvagmZ00s1Nm9nje7UmbmS0xs2+b2XEz6zKzx/JuUxbMrMnMvmtm/5F3W7JgZvPN7DkzOzH2t74n7zalzczWjh3Tb5jZHjObk3ebkmZmz5hZr5m9MWHfz5rZi2b2g7GvC5J4rmCC3MyaJH1F0scl3SFplZndkW+rUjcs6S/d/Tck3S3pcxV4zZL0mKTjeTciQ/8s6b/dfbmk31TJX7uZLZK0RlKbu39QUpOkT+fbqlTskvTgNfsel/Qtd79d0rfGtmMLJsglfUTSKXfvdvdBSV+T9HDObUqVu7/t7q+O/fuSRt/gi/JtVbrMbLGk35O0M++2ZMHMfkbSCklPS5K7D7r7xVwblY1Zkuaa2SxJLZL+L+f2JM7d2yX9+JrdD0v66ti/vyrpk0k8V0hBvkjS2QnbPSp5qE1kZksl3SXpWM5NSdsWSX8l6Sc5tyMryyT1Sfq3sXLSTjO7Ke9Gpcnd35L0T5J+KOltSe+4+6F8W5WZn3f3t6XRjpqkRO6KE1KQ2xT7KjF20sw+IOkbkr7g7u/m3Z60mNknJPW6e2febcnQLEkflvSv7n6XpPeU0Ol2UY3VhR+WdJukX5R0k5n9cb6tCltIQd4jacmE7cUq4enYtcysWaMh/qy7P593e1J2r6SHzOy0RktnHzWzf8+3SanrkdTj7lfPtJ7TaLCX2QOS3nT3PncfkvS8pN/KuU1Z+ZGZ3SpJY197k3jQkIL8O5JuN7PbzGy2Ri+O7M+5TakyM9No7fS4u3857/akzd3Xu/tid1+q0b/vYXcvdU/N3c9JOmtmvz62635J38uxSVn4oaS7zaxl7Bi/XyW/wDvBfkmPjP37EUnfTOJBg7n5srsPm9nnJb2g0avcz7h7V87NStu9kj4j6X/N7LWxfX/j7gfzaxJS8BeSnh3roHRL+pOc25Mqdz9mZs9JelWjI7O+qxJO1TezPZJWSrrFzHokfUnSP0j6upn9mUY/0D6VyHMxRR8AwhZSaQUAMAWCHAACR5ADQOAIcgAIHEEOAIEjyAEgcAQ5AATu/wHRry1uj1JMRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "sns.residplot(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept: 0.18626\n",
      "Slope : 0.93768\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArh0lEQVR4nO3deXiU5b3/8fd3skOAsAYIkTXBuitxRTFhEdt6qfV00aq19tdDrVoUW3v0nJ7luk77k/48R8BqF7u3tqXWhdraVgMxWq0bi7ghSQBZwi4EEkjIMvfvj0kwhJmQTGbmyczzeV2XF8mTZ+b53pN4f5/nXs05h4iI+E/A6wBERMQbSgAiIj6lBCAi4lNKACIiPqUEICLiU0oAIiI+pQQgIuJTSgAiIj6lBCAi4lNKAOJ7ZjbZzPaZ2Tnt3481s71mVhrm3HvM7PEux5aY2YPtX3/RzDaaWb2ZbTKz68O8x2gzO2xmwzsdm2Zme8wsI9blE4nEtBSECJjZPwN3AdOAp4C3nXPfCHPeeGAdMNo5d9DM0oBtwKeAt4EdwLnOufVmNgYY5px7N8z7/AX4k3PuB+3fLwLSnXNfi08JRY6nJwARwDn3Y6AaeA0YA/xbhPM2A6uBq9sPzQQOO+debf8+CJxmZjnOuR3hKv92vwRuAGhPItcBv45BUUR6TAlA5CM/Bk4DvuecO9LNeb8lVGEDfL79e5xzh4DPAbcAO8zsGTM7OcJ7/BE4xcwmAXOAA86512NQBpEeUxOQCGBmucBa4Hng48Dpzrl9Ec4dCWwBioB3gAudc+u6nJMDfBs4zzl3SYT3+RGh5qOTgXXOuW/HqDgiPaInAJGQJcAq59yXgWeAH0Y60Tm3B6gEfg5s6qj8zSzfzK40s4HAEaABaOvmmr8CvghcCTwagzKI9IoSgPiemV0FXE6o6QZCncHnhBvB08lvgdnt/3YIAF8HtgP7gEuBWyO9gXPuZUJ9Bqudcx9EG79ItNQEJOIhM6sAfuuc+4nXsYj/KAGIeMTMzgXKgULnXL3X8Yj/qAlIxANm9ktgOXCnKn/xip4ARER8Sk8AIiI+le51AL0xYsQIN2HChKhee+jQIQYOHBjbgPo5ldkfVGZ/6EuZV61atdc5N7Lr8aRKABMmTGDlypVRvbayspLS0tLYBtTPqcz+oDL7Q1/KbGabwx2PexOQmf3MzHab2Tudjg0zs3Izq27/d2i84xARkWMlog/gF4Qm2XR2D7DCOVcErGj/XkREEijuCcA59yKhWZGdXUVoNUTa/7063nGIiMixEjIM1MwmAH92zp3W/n2dcy6v08/3O+fCNgOZ2TxgHkB+fv60pUuXRhVDQ0MDubm5Ub02WanM/qAy+0NfylxWVrbKOVfS9Xi/7wR2zj0CPAJQUlLiou0EUaeRP6jM/qAyx4ZXCWCXmY1xzu1o3zVpt0dxiIh4btmaWu5/dj3b6xoZm5fD3XOncvXZBXG/rlcTwZ4Gbmr/+iZCm2OIiPjOsjW13Pvk29TWNeKA2rpG7n3ybZatqY37tRMxDPR3wCvAVDPbZmb/B1gIzDGzakK7IS2MdxwiIv3R/c+up7Hl2G0jGlvauP/Z9XG/dtybgJxz10X40ax4X1tEpL/bXtd4wuPxGqyjtYBERDw0Ni8n4nHnHN/96/tM/dbf+FvVAaYvrIhp01C/HwUkInJcJ+mZ3e20mVzunjuVe598+5hmoJyMND51dgGzH3iBDXsOAbBqb4Dag6H+ASAmncR6AhCRfi1cJ2nt/sYT3gkvW1PL9IUVTLznmZjfOcfS1WcXcN81p1OQl4MBI3IzGTc0h4eerzla+QNcPi4IxLZ/QAlARPq1cJ2kQee6rQS9HFkTjavPLuDh68/hkuKR7G1opnp3A4Oyj22gyUz76OtI/Qa9pQQgIv1aTzpJu/JyZE1vvb3tAF/6xRtc/fDLvFi1h9ysdObPKuKlf5lJQTf9A7GgPgAR6dfG5uVQG6ay764SjCZpJNq72w+weHk15e/tAmBAZho3XTSBeZdMYujATCBy/8Ddc6fGJAYlABHp18JVggGzbivBaJJGolTtqmdReRV/fWcnANkZAb5w4QS+MmMSw3Ozjjm3o6M39ORST0GMZwkrAYhIv9a5EuwYBVQwtK3bSjDed87RqNndwJIV1fz5re04B5npAW44fzy3lE5i1KDsiK+7+uwCrj67gMrKSr52fWlMY1ICEJF+r6MS7FBZWXnC8wFP1tfpatPeQzy4opo/vllL0EFmWoBrzyvk1tIpjB4SueJPBCUAEUlJXZNGom3dd5gHV1Tz5Jpa2oKO9IBx7XmF3FY2JWLnbqIpAYiIxNC2/Yd5+Pka/rByG61BR1rA+FxJIbfPnMKqzfv57A9f8fyppIMSgIhIDOw40Mj3n9/A0je20NLmCBhcc3YB82cVMWHEwKNzEzr6JTrmJkBsZvVGQwlARKQPdh9s4vuVG/jt61tobg1iBleeOZY7ZhcxeeRHO3h1NzdBCUBEJApebaayt+EIP3phA796ZTNHWkPLNHzy9DHcMbuI4vxBx53fH+cmKAGISNLyolll36FmHnlxI7/8xwdHr3vZKfksmFPMx8YMjvi6/jg3QQlARJJWIppVOp4wausaGZSVTnNb8Ogd/6yTR3Hn7GJOHzfkhO/TH+cmKAGISMLFqtkm3s0qy9bUcs8Tb9HUXuHXH2kF4GNjBnHfNWdwVmFej9+rP81N6KAEICIJFctmm3g2qzQcaeU//vjO0cq/s4ONrb2q/Dt4PTehK60GKiIJFcuVOu+eO5WcjLRjjvW1WeVwcys/fGEDl3y3goNNrWHP6U+LyvWFngBEJKFi2WwTy2aVppY2Hn11Mz98YQN7G5qB0LINzW3HPwH0h0XlYkEJQEQSKtbNNn1tVmlqaWPp61t4uHIDe+qPAHBWYR53zSnmw4Yj/OtT7/S649aroam9pQQgIgnVX0bDNLcGeWzlVh5+voYdB5oAOK1gMHfNKaZs6ijMDAAz61Vl3h9n/EaiBCAiCeX1aJiWtiBPrNrG9ypqjj6JnDx6EHfNKWbOKflHK/7O8fYmtv444zcSJQARSbhIlWo8m05a24I8taaWByuq2bovVPEX5+eyYHYxc08dTSBgJ3iHnumPM34j8TQBmNkC4MuAA94GbnbONXkZk4h4I15NJ21Bx5/WbmfJimo27T0EwKSRA7ljVhFXnDGWtBhV/B3644zfSDxLAGZWAMwHTnHONZrZY8C1wC+8iklEvBPrppNg0PHM2ztYvLyKDXtCFf+E4QO4Y3YRV55ZEPOKv0N/6ePoCa+bgNKBHDNrAQYA2z2OR0RiqDdNOrFqOgkGHc+9t5NF5dWs31UPwLihOcyfVcQ1ZxeQnhbf6U9e93H0hjnnvLu42R3Ad4BG4Dnn3PVhzpkHzAPIz8+ftnTp0qiu1dDQQG5u7olPTCEqsz/01zLXNbZQu7+RYKc6JmBGwdAc8nIyjjt//c76sGPuM9MCTB197Oqa4crsnOPNPW08Vd3ClvrQ+wzLNq6cnMHFBemkx+mOP1H68nsuKytb5Zwr6XrcyyagocBVwESgDviDmd3gnHu083nOuUeARwBKSkpcaWlpVNerrKwk2tcmK5XZH/prmacvrKC2Lu244wV5abx8T+lxx+u69AFAqOnkvmtOp7TL3XPnMjvnqKzaw6LyKt7adhiA/MFZ3F42hc+eW0hW+vExJKN4/J69bAKaDWxyzu0BMLMngYuAR7t9lYgkhd426fS26cQ5x0s1e3mgvIo1W+oAGJGbxa2lk/n8+SeRnZEaFX88eZkAtgAXmNkAQk1As4CVHsYjIjEUzWiYno65X/dhG9//0au8/sE+AIYPzOQrl07ixgsmkJOpir+nPEsAzrnXzOxxYDXQCqyhvalHRJJfPEbDvPHBPh54ropXNjYBTeQNyOArMybzhQvHMzDL6zEtycfTT8w595/Af3oZg4jERyxHw6zesp9F5VX8vXovAAPS4atlxXxx+gQGZR/foSw9o5QpInHT14Xa3tpWx6LyKp5fvweA3Kx0vnTxRKZSyydnFcUqTN9SAhCRfued2gMsXl7N8nW7ABiQmcbN0yfwz5dMIm9AJpWV0U8ZSpaVOhNBCUBE+o33dx5kcXk1f3t3JwDZGQFuunAC82ZMYnhuVp/fP5lW6kwEJQAR8Vz1rnoWr6jmmbd2AJCVHuCGC8Zzy6WTGTmo7xV/h1gvN5HsTxNKACLimY17Gliyopqn127HudCs38+ffxJfLZ1M/uDsmF8vlit1psLThBKAiCTclg8Ps2RFNU+t2UbQQUaa8dmSQm4rmxLXVTO7m5vQ27v5ZFr3PxIlABHps55Wntv2H+ahihr+sGobbUFHWsC4tmQct5VNoXDYgLjHGWluQtnJI3t9N59M6/5HogQgIn3Sk6aQHQcaeaiihsdWbqWlzREw+KdzxjF/1hTGDx+YsFgjzU2I5m4+mdb9j0QJQET6pLvK86LJw/l+5QZ++9oWmtuCmMHVZ41l/qwiJo30ZgXTcHMTFvz+zbDndnc3n0zr/keiBCAifRKpkqyta+SS//c8R1pDFf8VZ4zhztlFTBk1KOz5Xop23SJIjnX/I1ECEJE+iVR5AhxpDfLx00Zzx+wiTh49OMGR9Vy0d/N9nensNSUAEemTu+dO5Z4n3qKp9djNXE4bO5jvfvoMTh07xKPIei4V7uajoQQgIlE70NgS2mi902ZbWekBbiubwvwkW6sn2e/mo6EEICK9Vt/Uwi9e/oAf/30jB5taAbikaAR3zi5m2vihHkcnPaUEICI9dri5lV/+YzM/enEDdYdbALhw0nAWzCnmvInDPI5OeksJQCTFxWK9msbmNh59dTM/fGEDHx5qBuDcCUNZMKeYiyaPiEfYkgBKACIprK/r1TS1tPG717fw/coN7Kk/AsBZhXl8/bJiLp4yAjM7wTtIf6YEIJLCol2v5khrG4+t3MbDFTXsPNgEwBnjhrBgTjGlxSNV8acIJQCRFNbb9Wpa2oI8vmobD1XUHB3b/7Exg7lrTjENTS1866l3fDVMMtUpAYiksJ7OcG1tC/Lkmlq+V1HN1n2h84vzc1kwu5i5p47m6bXb+den3knqpY/leEoAIkkkXIduXjfnn2iGa1vQ8fTaWpYsr+aDDw8DMGpQFi1tQap2NfDtZ9ZxpDWYEksfy/GUAESSRKQO3fsuSov4mkgzXK88cyx/Wrudxcur2LDnEAATRwxk+uThPL5q29FZvR3X6Fr5d0impY/leEoAIkki0l34rgMt3b6u8wzXYNDx7Ls7uXzJi1TtagCgcFgO82cW8amzC7j0/srjlnRobGkjzYw2545772Ra+liOpwQgkiQi3W03twXDHu/MOUf5e7tYtLyadTsOAjB0QAbOwdZ9jSxeXk1GWiDiNdqcIycjLamXPpbjBby8uJnlmdnjZva+ma0zswu9jEekP4t0t52ZFvl/Y+ccFe/v4sqHXmber1exbsdBRg/O5tPTxtHY3EZdY+jpoaOpJ29ARtj3KcjL4b5rTqcgLwfr9L3a/5Ob108AS4C/Oec+bWaZQPz3hBNJUpE6dPOHZB53rnOOF6v38kB5FWu31gEwclAWt5ZO5rrzTmLW/74QtqknKz0Q8U7fj4ulpTrPEoCZDQZmAF8EcM41A81exSPS30Xq0M07UH30HOccr2z4kAfKq1i5eT8Awwdm8tXSyVx//nhyMkMdxpGaeg40trDoc2f5bllkvzIXpmMnIRc2Owt4BHgPOBNYBdzhnDvU5bx5wDyA/Pz8aUuXLo3qeg0NDeTmerMFnVdUZn/oKPP6fW08Wd3M+v2hO/vcDPjExAxmnZRBVvqxM3fX76wP23eQmRZg6uj+t2NXV37+PUejrKxslXOupOtxLxNACfAqMN0595qZLQEOOuf+PdJrSkpK3MqVK6O6XmVlJaWlpVG9NlmpzP7w46dWUPnhQF6u+RCAITkZzJsxiZsumkBuVviH/K5DSiHU1JMs7fp+/D33pcxmFjYBeNkHsA3Y5px7rf37x4F7PIxHJKm8ubWOReVVvFDVBDSRnREgIxDgQGMLv31tCwV5ORErc7/ugCXH8iwBOOd2mtlWM5vqnFsPzCLUHCQi3Xin9gCLyqtY8f5uALLTYMbUfF6s2kP9kdDmLD1ZqkGduuL1KKCvAb9pHwG0EbjZ43hE+q11Ow6yeHkVz767Cwg12dx00QROTdvBwtUHw47q0VIN0h1PE4Bz7k3guHYpEflI1a56liyv5pm3dwChPXdvvGA8t5ROZkRuFpWVO9ledzjsa7VUg3TH6ycAEYlgw54GHlxRzdNrt+McZKYH+Px5J3Fr6WRGDc4+5tyervop0pkSgEg/88HeQzxYUc2yNbUEHWSkGdeeexK3lk1mzJDwFfqJVv0UCUcJQKSf2LrvMA9V1PD46m20BR3pAeOzJeO4feYUxg3tfpK8RvVINJQARDy2va6Rh56v4bE3ttIadKQFjM9MG8fXZhZx0vCer46iUT3SW0oAIh7ZdbCJbzy2lr/X7D16rGT8UO7/zJlMHDHQw8jEL5QARBJsT/0RflC5gV+/+gEtbcfOxH93+0HWbq1TApCEUAIQSZAPG47wyIsb+eUrH9DUEn4Nf43dl0RSAhCJs/2Hmvnx3zfyi398wOHm0CidOafkU/7errDna+y+JIoSgEicHGhs4acvbeJnL22ioX2JhrKpI7lrzlROHzeE6QsrNHZfPKUEIBJj9U0t/PzlD/jx3zdS3xSq+GcUj2TB7CLOPmno0fM0dl+8pgQgEiO/f30L3/nLOg62V/oAF00ezoI5xZw7Ydhx52vsvnhNCUCkjxqb2/iXJ97i6bXbjzmemRbgsyWFYSv/Dhq7L15SAhCJUlNLG795bQs/qNzA3oYjx/28uS2oET3SrykBSFJatqa2100nkV7T2/c60trG79/YysPP17Dr4PEVf2ca0SP9mRKAJJ2u2xl23vwkr5evWbl5H0+sqg37Xl2TQHNrkD+s2spDFTXsONAEwCljBnPXnGL+44/vsL39WGca0SP9mRKAJJ37n11/zMgZ+GgC1XcuCPTqNb97bSttXfbF7joZq6UtyFOra3mwoppt+0N39CePHsSds4u47JTRBAJGw5FWjeiRpKMEIEknUrNK6Hj4JRQivaZr5d/5/Na2IH98czsPVlSz+cPQhitTRuVy5+wiPnHaGAIBO3q+RvRIMlICkKQTzeYnkV6TZhY2CeQNyOCyxS+ycc8hACaNGMgds4u44oyxpHWq+DvTiB5JNuGfl0X6sbvnTiUnI+2YYydqbon0muvOLzzuuAH7D7ewcc8hTho2gP/5zJk8t2AGV51VELHyP5Fla2qZvrCCifc8w/SFFSxbUxvV+4jEkp4AJOl019yy7K/vMX1hxXHHu3vNtJOG8t9/Xse+w80AOKAgL4f5s6ZwzTnjyEjr231Sd53WemIQLykBSFIK19yybE0ttfsbqa0L3dF3rWi7vsY5x4p1u/jJS5uOVv5jhmRz+8wpfGZaIZnpsXlA7q7TWglAvKQEICnj/mfXc21h9yN6IFTxv1C1h0XlVazddgCAUYOyuLV0MteedxLZXZqE+qr7TuvwopnnINJbSgCSMrbXNUJhhOOEKv6Xaz7kgfL1rN5SB8CI3Ey+WjqF68+PfcXfobed1moykkRRApATSpa70VCFWh/2+KsbP+SB8ipe37QPgKEDMrjl0snceOF4BmT27n+D3n4evV31U01GkiieJwAzSwNWArXOuSu8jkeOlUx3o3fPnUrtulXHHMtMCzAgM41rH3kVgCE5GcybMYmbLppAblbv//yj+Tx6O0cgmiYjkWh4ngCAO4B1wGCvA5HjJepuNBZPGVefXcCyne9RkJdGbV0jWekBjrQGqd7dwKDsdL588SRuvngCg7Mzoo4z2s+jN3MEopnnIBINT+cBmNk44JPAT7yMQyJLxN1ox111bV0jjo/uqqMZK1/XHGDq6EEAHGkNkpuVzvyZU3jpmzO5Y3ZRnyp/SMznEc08B5FomIswFT4hFzd7HLgPGAR8I1wTkJnNA+YB5OfnT1u6dGlU12poaCA3N7cP0SafWJR5/c56mtuO38A8M+2jiravYnGNLQfbWFbTwurdbe2vhTknZfDxiRnkZkY3eStesfZEXWMLuw400dwWJDMtQP6QbPJywicv/W37Q1/KXFZWtso5V9L1uGdNQGZ2BbDbObfKzEojneecewR4BKCkpMSVlkY8tVuVlZVE+9pkFYsy13Vp84bQ3eh915xOaYyagG6+5xlcmIdRAzYtLO32tVW76lm8vIq/vL0TgMwAfPHiScybMYkRuVkxia+zRHwevaW/bX+IR5m97AOYDlxpZp8AsoHBZvaoc+4GD2NKOvEeoZOIRc6iafOu2d3Agyuq+dNb23EOMtMD3HD+eM7I3MnVcz8Wk7i6+2yTYVSUyIl4lgCcc/cC9wK0PwF8Q5V/7yRqhE68FznrzTDJD/Ye4sEV1Sx7s5agCzW9XHdeIV8tncLoIdlUVu6OSUwn+mxV4UsqOGECMLPbgd845/YnIB7phVQZL96Tu+qt+w7zvYpqnlhdS1vQkR4wrj2vkNvKplAQh9ExqfLZinSnJ08Ao4E3zGw18DPgWRfjnmPnXCVQGcv39INUGi8e6a66tq6Rhypq+MPKrbQGHWkB43Mlhdw+cwqFwwbELZ5U+mxFIjlhAnDOfcvM/h24DLgZeMjMHgN+6pzbEO8AJbJkHC/e0z6LnQea+H5lDUtf30pzW5CAwTXnFDB/ZhETRoTf9CWWkvGzFemtHvUBOOecme0EdgKtwFDgcTMrd859M54BSmS9XWLAaz3ps9hd38QPKjfwm9e20NwaxAyuPHMsd8wuYvLIxA37S7bPViQaPekDmA/cBOwlNGHrbudci5kFgGpACcAjyTYipbt29YuLRvCjFzbw61c309QSGmf/ydPHcFrBYB59dQuz//eFhJYv2T5bkWj05AlgBHCNc25z54POuWD7WH7xUDKNSInUfl5b18gl333+aHK47JR8FswpZv3Oek/XIUqmz1YkGj3pA/iPbn62LrbhSCqL1K4OoSeBWSeP4s7ZxZw+bggAX/7lSo3EEYkj7QksCXP33Klkh9ll6+TRg1h223R++sVzj1b+oJE4IvHWH1YDFR9oONJKbV0jgU6bqmelB7jl0sksmFMc9jUaiSMSX0oAEleHm1v59Sub+dGLG9l3KLTv7nkTh3HXnGIumDS829dqJI5IfCkBSFw0tbTx6Kub+eELG9jbEKr4p40fytfnFHPh5OGYnXiFTo3EEYkvJQAPJMsWi9Foamnj929s5eHna9hdfwSAMwvzuGtOMTOKRvSo4u9MI3FE4kcJIMGSaYvF3mhuDfLYylDFv+NAEwCnjh3MXXOKmXnyqF5X/CISf0oACZZqi4y1tAV5YtU2vldRc7TD9uTRg1gwp5jLTslXxS/SjykBJFiqDG1sbQvy1JpavldRw5Z9hwEoGpXLgjnFXH7q6GNG+4hI/6QEkGDJPrSxLej409rtLFlRzaa9hwCYNGIgd8wu4oozxpKmil8kaSgBJFiyDm0MBh3PvL2Dxcur2LAnVPGPHz6AO2YVceWZY0lP05xCkWSjBJBgyTa0MRh0PPfeThaVV7N+Vz0A44bmMH9mEZ86p4AMVfwiSUsJwAOJGNq4bE0tu3bWc/M9z0SVZJxzrFi3mwfKq3hvx0EAxgzJ5vaZU/jMtEIywyzpICLJRQkgBXUMNb315CCOQK+GmjrnqKzaw6LyKt7adgCAUYOyuH3mFD53biFZ6Wlxj19EEkMJIAVFM9TUOcdLNXt5oLyKNVvqABiRm8WtpZP5/PknkZ2hil8k1SgBpKDeDjV9ZcOHLCqv4vUP9gEwbGAmt1w6iRsvmEBOpip+kVSlBJCCejrU9I0P9vHAc1W8svFDAPIGZDBvxiRuunACA7P0pyGS6vR/eQrqGGoa2r45pPNQ09Vb9rOovIq/V+8FYFB2OvMumcQXp09gUHaGFyGLiAeUAFJQRzv/rvWrMTg6CmjSyIHc/PPXeX79HgAGZaXzpYsn8qWLJzIkRxW/iN8oAaSoq88uoPJANZsWlvLu9gMsXl5N+Xu7ABiQmcbN0yfwz5dMIm9ApseRiohXPEsAZlYI/AoYDQSBR5xzS7yKJxVtqw/y1UdX8dd3dgKQnRHgpgsnMG/GJIbnZnkcnYh4zcsngFbg68651WY2CFhlZuXOufc8jKnf6s0eAjW761m8vJpn3mrE0UhWeoAbLhjPLZdOZuSg+Ff8qbzfgUgq8SwBOOd2ADvav643s3VAAaAE0EVP9xDYtPcQS5ZX8ce123EO0g2uv2A8t5ZNIX9wdlzi6lrRAym534FIKjLnnNcxYGYTgBeB05xzB7v8bB4wDyA/P3/a0qVLo7pGQ0MDubm5fYzUG+t31tPcFjzueGZagKmjB7H7cJCnN7Twj+2tBB2kGcwYl87M0S0UDu95mesaW9h1oInmtiCZaQHyh2STF6FzuK6xhdr9jQQ7/f0EzAgYtAaP/5vqiDXekvn3HC2V2R/6UuaysrJVzrmSrsc9TwBmlgu8AHzHOfdkd+eWlJS4lStXRnWdyspKSktLo3ptIoW7q17w+zeJ9Fu67rxC/rByG61BR1rA+My0cdxWNoXCYQN6VeauTxkQGjp63zWnh71zn76wIuxcg0gM2LTwkz0+P1rJ8nuOJZXZH/pSZjMLmwA8HQVkZhnAE8BvTlT5+0Gkpp68ARnsP9wS9jW/e30rAYN/Omcc82dNYfzwgVFdu7fLR/R2A5tk2e9AxE+8HAVkwE+Bdc65B7yKoz+JVAlnpQfIyUg77mcGXHXWWObPKmLSyL49Dvd2+YhIs43zcjI40hpMuv0ORPzIyzV9pwM3AjPN7M32/z7hYTyei1TZ1jW2cP7EYcccO6swj+cWzGDxtWf3ufKHyHfokY7fPXcqOV0WiMvJSOO/rjyV+645nYK8HAwoyMuJ2IwkIt7ychTQS4RuYqVdpLtqAyqrQrN3556az4I5xZw8enBMr93bncpOtLGNKnyR/k8zgfuRcJUwgANmfyyfO2cXcVrBkLhcO5qdyhKxsY2IxI8SQAzEauJT2cmjmFE8gufe3XV01M/Hxgxi4TVncGZhXkxjDkcVuoi/KAH0QrwmPjUcaeUXL2/ikRc3crAptILnxVNGsGBOMdPGD41DSURElAB6LNIQzeyMQK933+pw6Egrv3plMz96cQN17cM8z584jLvmFHP+pOFxK4eWaRARUALosUhDNLse69DdOPnG5jYefXUzP3xhAx8eagagZPxQ7rqsmIsmj4hd0F30dEkJEfEHJYAeisXEp6aWNn73+ha+X7mBPfVHgNBwzrvmFHNJ0QhCUyPiJ5q9gkUkdSkB9FBfJj4daW3jsZXbeLiihp0HmwA4vWAIC+YUUTZ1VNwr/g69newlIqlNCaCHIo2T/68rTwXCD59saQvy+KptPFRRczR5fGzMYBbMLmLOKfkJq/g79HSvYBHxByWAHurNxKfWtiCPrdzKd//6/tE2/vSAceMF4/n3K04hEPBm/ltvJ3uJSGpTAuiFE42Tbws6nl5by5Ll1Xzw4eFjftYadCx9YytnFuZ51t4ezWQvEUldSgAxEAw6/vz2DpYsr2LDnkMApAWMti7r4veHDldN9hKRDkoAfRAMOp59dyeLlldRtasBgHFDc5g/q4hvPv5W2Neow1VE+gslgCg45yh/bxeLllezbkdoA7OxQ7K5fWYRn542jsz0AEuWV6vDVUT6NS+Xg06IZWtqmb6wgrdrDzB9YQXL1tRG/V7OOZ5/fzdXPvQy8369inU7DpI/OIv/vupUnr+7lM+ffxKZ6aGPNNJyyepwFZH+IqWfAI6Z+VoY/cxX5xx/r97LA+VVvLm1DoCRg7K4tXQy1513EtldKvrO768OVxHpr1I6AcRi5us/akIV/8rN+wEYPjCTWy6dzA0XjCcn8/iKvzN1uIpIf5bSCaAvM19f37SPB8rX8+rGfQDkDcjgKzMm84ULxzMwK6U/NhHxiZSuyaKZ+bpq834WlVfxUs1eAAZnp/OVSydz00UTyFXFLyIpJKVrtN7MfF27tY5Fy6uoXB/aenFQVjpfungiX7p4IkNyMhIWs4hIoqR0AujcEQv1FITpiH2n9gCLl1exfN1uAAZmpnHz9Il8+ZKJ5A3I9CJsEZGESOkEAB91xFZWVvK160uPHn9/50EWl1fzt3d3AqEngy9cNJ6vzJjMsIGq+EUk9aV8Auiqelc9i1dU88xbOwDISg9w4wXj+cqlkxk5KMvj6EREEsc3CWDnoSB3LF3D02u34xxkpgX4/PkncWvpZEYNzvY6PBGRhEv5BNDY3Ma3lr3Dk6sbcTSSkWZ87txCbiubwpghWpZBRPzL0wRgZpcDS4A04CfOuYWxvkZ2RoANexoIGHy2veIfN3RArC8jIpJ0PEsAZpYGPAzMAbYBb5jZ086592J8Hf7vp07nnTdX8tlPnBHLtxYRSWpeLgZ3HlDjnNvonGsGlgJXxeNCp4wdzKgBKb/unYhIr5hz7sRnxePCZp8GLnfOfbn9+xuB851zt3c5bx4wDyA/P3/a0qVLo7peQ0MDubm5fQs6yajM/qAy+0NfylxWVrbKOVfS9biXfQDhNsY9Lhs55x4BHgEoKSlxpaWlUV2ssrKSaF+brFRmf1CZ/SEeZfayXWQbUNjp+3HAdo9iERHxHS8TwBtAkZlNNLNM4FrgaQ/jERHxFc+agJxzrWZ2O/AsoWGgP3POvetVPCIifuPpPADn3F+Av3gZg4iIX6X8TOBoLFtTq60cRSTlKQF0ccw+wkS/j7CISH+n2VFddLePsIhIKlEC6KIv+wiLiCQTJYAuIu0X3N0+wiIiyUgJoIu7504lJyPtmGOR9hEWEUlm6gTuovM+whoFJCKpTAkgjI59hEVEUpmagEREfEoJQETEp5QARER8SglARMSnlABERHxKCUBExKeUAEREfEoJQETEp5QARER8SglARMSnlABERHxKCUBExKeUAEREfEoJQETEp5QARER8SglARMSnPEkAZna/mb1vZm+Z2VNmludFHCIifubVE0A5cJpz7gygCrjXozhERHzLkwTgnHvOOdfa/u2rwDgv4hAR8TNzznkbgNmfgN875x6N8PN5wDyA/Pz8aUuXLo3qOg0NDeTm5kYdZzJSmf1BZfaHvpS5rKxslXOupOvxuCUAM1sOjA7zo39zzv2x/Zx/A0qAa1wPAikpKXErV66MKp7KykpKS0ujem2yUpn9QWX2h76U2czCJoD0vgYViXNu9gkCugm4ApjVk8o/1patqeX+Z9ezva6RsXk53D13KlefXZDoMEREPBO3BNAdM7sc+BfgUufc4URff9maWu598m0aW9oAqK1r5N4n3wZQEhAR3/BqFNBDwCCg3MzeNLMfJvLi9z+7/mjl36GxpY37n12fyDBERDzlyROAc26KF9ftsL2usVfHRURSkS9nAo/Ny+nVcRGRVOTLBHD33KnkZKQdcywnI4275071KCIRkcTzpAnIax0dvRoFJCJ+5ssEAKEkoApfRPzMl01AIiKiBCAi4ltKACIiPqUEICLiU0oAIiI+5fly0L1hZnuAzVG+fASwN4bhJAOV2R9UZn/oS5nHO+dGdj2YVAmgL8xsZbjlUFOZyuwPKrM/xKPMagISEfEpJQAREZ/yUwJ4xOsAPKAy+4PK7A8xL7Nv+gBERORYfnoCEBGRTpQARER8yhcJwMwuN7P1ZlZjZvd4HU+8mVmhmT1vZuvM7F0zu8PrmBLBzNLMbI2Z/dnrWBLBzPLM7HEze7/9d32h1zHFm5ktaP+bfsfMfmdm2V7HFGtm9jMz221m73Q6NszMys2suv3fobG4VsonADNLAx4GPg6cAlxnZqd4G1XctQJfd859DLgAuM0HZQa4A1jndRAJtAT4m3PuZOBMUrzsZlYAzAdKnHOnAWnAtd5GFRe/AC7vcuweYIVzrghY0f59n6V8AgDOA2qccxudc83AUuAqj2OKK+fcDufc6vav6wlVDCm9+YGZjQM+CfzE61gSwcwGAzOAnwI455qdc3WeBpUY6UCOmaUDA4DtHscTc865F4F9XQ5fBfyy/etfAlfH4lp+SAAFwNZO328jxSvDzsxsAnA28JrHocTbYuCbQNDjOBJlErAH+Hl7s9dPzGyg10HFk3OuFvgfYAuwAzjgnHvO26gSJt85twNCN3jAqFi8qR8SgIU55ouxr2aWCzwB3OmcO+h1PPFiZlcAu51zq7yOJYHSgXOAHzjnzgYOEaNmgf6qvd37KmAiMBYYaGY3eBtVcvNDAtgGFHb6fhwp+NjYlZllEKr8f+Oce9LreOJsOnClmX1AqIlvppk96m1IcbcN2Oac63iye5xQQkhls4FNzrk9zrkW4EngIo9jSpRdZjYGoP3f3bF4Uz8kgDeAIjObaGaZhDqNnvY4prgyMyPUNrzOOfeA1/HEm3PuXufcOOfcBEK/3wrnXErfGTrndgJbzWxq+6FZwHsehpQIW4ALzGxA+9/4LFK847uTp4Gb2r++CfhjLN405TeFd861mtntwLOERg38zDn3rsdhxdt04EbgbTN7s/3Yvzrn/uJdSBIHXwN+035jsxG42eN44so595qZPQ6sJjTSbQ0puCSEmf0OKAVGmNk24D+BhcBjZvZ/CCXCz8TkWloKQkTEn/zQBCQiImEoAYiI+JQSgIiITykBiIj4lBKAiIhPKQGIiPiUEoCIiE8pAYj0gZmda2ZvmVm2mQ1sX6v+NK/jEukJTQQT6SMz+zaQDeQQWp/nPo9DEukRJQCRPmpfiuENoAm4yDnX5nFIIj2iJiCRvhsG5AKDCD0JiCQFPQGI9JGZPU1oGeqJwBjn3O0ehyTSIym/GqhIPJnZF4BW59xv2/ef/oeZzXTOVXgdm8iJ6AlARMSn1AcgIuJTSgAiIj6lBCAi4lNKACIiPqUEICLiU0oAIiI+pQQgIuJT/x+ymQLR1vGJ6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# template on our dataset\n",
    "\n",
    "from sklearn import linear_model\n",
    "# from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# First initialize the model.\n",
    "linear_model = linear_model.LinearRegression()\n",
    "\n",
    "# define x and y\n",
    "x_input = sim_data['x'].values.reshape(n_points, 1)\n",
    "y_output = sim_data['y'].values.reshape(n_points, 1)\n",
    "\n",
    "# fit the model to the existing data with lm.fit(x,y)\n",
    "linear_model.fit(x_input, y_output)\n",
    "\n",
    "# Use the fit model to make predictions\n",
    "y_pred = linear_model.predict(x_input)\n",
    "\n",
    "# Plot output\n",
    "plt.scatter(x_input, y_output)\n",
    "plt.plot(x_input, y_pred, linewidth=2)\n",
    "plt.grid(True)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('x vs y')\n",
    "\n",
    "# Print model parameters by viewing attributes intercept and coef\n",
    "# Get slope (m) and y-intercept (b)\n",
    "print('Intercept: {0:.5f}'.format(linear_model.intercept_[0]))\n",
    "print('Slope : {0:.5f}'.format(linear_model.coef_[0][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stats Model Package\n",
    "Now we show how to fit a line with Stats Models package. The model declaration and fitting are in a different format as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept: 0.18626\n",
      "Slope : 0.93768\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.formula.api as sm\n",
    "\n",
    "# ols = ordinary least squares\n",
    "ols_model = sm.ols(formula = 'y ~ x', data=sim_data)\n",
    "# Alternatively:\n",
    "#ols_model = sm.OLS(y_output, sm.add_constant(x_input))\n",
    "\n",
    "# fit the model\n",
    "results = ols_model.fit()\n",
    "\n",
    "# Get slope (m) and y-intercept (b)\n",
    "print('Intercept: {0:.5f}'.format(results.params.Intercept))\n",
    "print('Slope : {0:.5f}'.format(results.params.x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    0.186260\n",
       "x            0.937684\n",
       "dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>predicted</th>\n",
       "      <th>resids</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.680006</td>\n",
       "      <td>0.186260</td>\n",
       "      <td>0.493746</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.204082</td>\n",
       "      <td>-2.282142</td>\n",
       "      <td>0.377624</td>\n",
       "      <td>-2.659766</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.408163</td>\n",
       "      <td>-0.522395</td>\n",
       "      <td>0.568988</td>\n",
       "      <td>-1.091383</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.612245</td>\n",
       "      <td>1.340155</td>\n",
       "      <td>0.760352</td>\n",
       "      <td>0.579802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.816327</td>\n",
       "      <td>1.379675</td>\n",
       "      <td>0.951716</td>\n",
       "      <td>0.427958</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x         y  predicted    resids\n",
       "0  0.000000  0.680006   0.186260  0.493746\n",
       "1  0.204082 -2.282142   0.377624 -2.659766\n",
       "2  0.408163 -0.522395   0.568988 -1.091383\n",
       "3  0.612245  1.340155   0.760352  0.579802\n",
       "4  0.816327  1.379675   0.951716  0.427958"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add predicted to pandas dataframe\n",
    "sim_data['predicted'] = y_pred\n",
    "# Add residuals to pandas dataframe\n",
    "sim_data['resids'] = y_output - y_pred\n",
    "\n",
    "# View head of data frame\n",
    "sim_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m = 0.9376837979830364\n",
      "b = 0.18626012709053796\n"
     ]
    }
   ],
   "source": [
    "# What are the coefficients?  (m = slope, and b = y-intercept)\n",
    "m = linear_model.coef_[0]\n",
    "b = linear_model.intercept_\n",
    "print('m = {}'.format(m[0]))\n",
    "print('b = {}'.format(b[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting the Parameters!\n",
    "\n",
    "It is very important to know how to interpret the parameters.  Here is an example for the above parameters. (Your results will be different because of the random seed in creating the data.)\n",
    "\n",
    "**Y Intercept**:  When x is zero, y is -0.53...\n",
    "    \n",
    "**Slope**:  When we increase x by 1 then we expect that y will go up by 1.07...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Summary Statistics\n",
    "The code in the cell below is fairly volumonous, but straight forward. In summary, the code computes summary statistics and makes diagnostic plots for ordinary linear models.\n",
    "\n",
    "The Scikit-learn library is primarily made for fitting models and getting predictions. For linear regression, we are also interested in many diagnostic plots.  To get these plots easily, we can use a different library called `statsmodels`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intercept, Slope : Intercept    0.186260\n",
      "x            0.937684\n",
      "dtype: float64\n",
      "\n",
      "SSE, SST, SSR, and RMSE:\n",
      "SSE: 381.3077904336268\n",
      "SST: 421.30484900582337\n",
      "SSR: 39.99705857219658\n",
      "RMSE: 19.527103995053306\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.905\n",
      "Model:                            OLS   Adj. R-squared:                  0.903\n",
      "Method:                 Least Squares   F-statistic:                     457.6\n",
      "Date:                Wed, 13 Jan 2021   Prob (F-statistic):           3.45e-26\n",
      "Time:                        07:39:17   Log-Likelihood:                -65.366\n",
      "No. Observations:                  50   AIC:                             134.7\n",
      "Df Residuals:                      48   BIC:                             138.6\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      0.1863      0.254      0.732      0.468      -0.325       0.698\n",
      "x              0.9377      0.044     21.392      0.000       0.850       1.026\n",
      "==============================================================================\n",
      "Omnibus:                        1.496   Durbin-Watson:                   1.582\n",
      "Prob(Omnibus):                  0.473   Jarque-Bera (JB):                0.788\n",
      "Skew:                          -0.264   Prob(JB):                        0.674\n",
      "Kurtosis:                       3.315   Cond. No.                         11.7\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.formula.api as sm\n",
    "import seaborn as sns\n",
    "\n",
    "ols_model = sm.ols(formula = 'y ~ x', data=sim_data)\n",
    "# Alternatively:\n",
    "#ols_model = sm.OLS(y_output, sm.add_constant(x_input))\n",
    "\n",
    "results = ols_model.fit()\n",
    "\n",
    "# Get slope (m) and y-intercept (b)\n",
    "print('Intercept, Slope : {}'.format(results.params))\n",
    "\n",
    "print('\\nSSE, SST, SSR, and RMSE:')\n",
    "mean_y = np.mean(y_output)\n",
    "sst = np.sum((y_output - mean_y)**2)\n",
    "sse = sst - results.ssr\n",
    "print('SSE: {}'.format(sse))\n",
    "print('SST: {}'.format(sst))\n",
    "print('SSR: {}'.format(results.ssr))\n",
    "print('RMSE: {}'.format(np.sqrt(results.mse_model)))\n",
    "\n",
    "# Get most of the linear regression statistics we are interested in:\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/anajaved/opt/miniconda3/lib/python3.8/site-packages/seaborn/distributions.py:2551: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Residual Histogram')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxR0lEQVR4nO3dd3xc9Znv8c8zoy6rWJZkySqWe8Ed4UIzNWCH0JIQekI2IU6W5O5ushtSliSbzW5yl7vphLC5hAsECC3EAdMx3TaSe5FlbNmWZMlWsyWrl3nuHzMGIWRpZOvozGie9+ull2fOOTrnO7I0z5zz+53fT1QVY4wxkcvjdgBjjDHuskJgjDERzgqBMcZEOCsExhgT4awQGGNMhLNCYIwxEc4KgQlbInKTiLw0wPrXReRLw3CcC0Sk8jT30Swik083izFOsEJgRoSIHBCRtsAb4mEReUBExpzOPlX1T6r6ieHKeKpEREVkap9lPxSRh088V9Uxqlo2yH5Ou+AYcyqsEJiR9ClVHQMsABYC33E3TmQREa/bGUxoskJgRpyqHgZexF8QABCRpSLyrogcE5GtInJBr3VfEJEyETkuIvtF5KZey9/utd2lIrJbRBpF5DeA9Fr3kU/oIlIQ+CQfFXh+m4iUBI5RJiJfGc7X3PusQURWisiuwLEOici3RCQReB6YEDhrahaRCSISKyK/EJGqwNcvRCS2137/RUSqA+u+1Oc4D4jI70RkjYi0ABeKyCdFZLOINIlIhYj8sJ+fyW2BdUdFZJWInCUi2wL/N78Zzp+LCQ1WCMyIE5FcYAWwN/A8B3gO+HcgDfgW8JSIZATeIH8FrFDVJOBsYEs/+0wHngK+D6QD+4BzhhCrBrgCSAZuA34uIotO5fUF4f8CXwm8njnAa6ragv9nUhW4jDRGVauA7wFL8RfN+cBi/K8REbkc+CfgEmAqsLyfY90I/ARIAt4GWoBbgVTgk8BXReTqPt+zBJgGfA74RSDDJcAZwHUi0t9xTBizQmBG0jMichyowP/G+4PA8puBNaq6RlV9qvoyUAysDKz3AXNEJF5Vq1V1Zz/7XgnsUtUnVbUL/xvY4WCDqepzqrpP/d4AXgLOG8Jr2xT4xHxMRI4Bdw6wbRcwW0SSVfWoqm4aYNubgH9T1RpVrQV+BNwSWHcd8EdV3amqrYF1ff1VVd8J/FzbVfV1Vd0eeL4NeJSPF5AfB7Z9CX/heDRw/EPAW/gv65lRxAqBGUlXBz4FXwDMxP/JHWAi8Nk+b6TnAtmBT8qfA1YB1SLynIjM7GffE/AXGADUP5piRT/b9UtEVojIehFpCBx/Za98wVikqqknvoCfDrDtpwP7Pygib4jIsgG2nQAc7PX8YGDZiXW9X2N/r/cjy0RkiYisFZFaEWnE/3Pt+zqP9Hrc1s/z02rkN6HHCoEZcYFP3A8AdwcWVQAP9X4jVdVEVf1pYPsXVfVSIBvYDfxPP7utBvJOPBER6f0c/yfbhF7Ps3ptG4v/stLdwPjAG/kaerUxDCdVLVLVq4BM4Bng8ROr+tm8Cn+hPCE/sAz8rzm317rer/eDw/V5/giwGshT1RTgXhx6nSZ8WCEwbvkFcKmILAAeBj4lIpeJiFdE4gJdKXNFZLyIXBloK+gAmoGefvb3HHCGiFwbaAD+Br3e7PG3K5wvIvkiksJHeyzFALFALdAtIisAR7qlikiM+O9/SAlcwmriw9dzBBgXyHfCo8D3A+0l6cBd+H9e4C8gt4nILBFJCKwbTBLQoKrtIrIYfxuCiXBWCIwrAte7HwT+VVUrgKuA7+J/M64A/hn/76cH+Cb+T8EN+K9nf62f/dUBn8V/SaYef2PnO73Wvwz8GdgGbASe7bXuOP7C8ThwFP+b4+rhfL193AIcEJEm/Jdmbg7k2I3/jb8scIlsAv4G9OJA7u3ApsAyVPV5/A3pa/E3vK8L7L9jgGN/Dfi3QFvNXXx4NmIimNjENMaMDiIyC9gBxKpqt9t5TPiwMwJjwpiIXBO43DQW+BnwNysCZqisEBgT3r6C/3LaPvxtDV91N44JR3ZpyBhjIpydERhjTISLcjvAUKWnp2tBQYHbMYwxJqxs3LixTlUz+lsXdoWgoKCA4uJit2MYY0xYEZGDJ1tnl4aMMSbCWSEwxpgIZ4XAGGMinBUCY4yJcFYIjDEmwjlaCETkchEpFZG9IvKxiToCI0w2isiWwFcwoycaY4wZRo51HxX/RNm/BS4FKoEiEVmtqrv6bPqWql7hVA5jjDEDc/KMYDGwV1XLVLUTeAz/UMPGGGNCiJOFIIePTpNXGVjW1zIR2Soiz4vIGQ7mMcYY0w8n7yzub/q7viPcbQImqmqziKzEP23ftI/tSOR24HaA/Pz8YY5pjPse2VA+ose7cYn9HZkPOXlGUMlH51DN5cO5VgFQ1SZVbQ48XgNEB6bjo89296lqoaoWZmT0O1SGMcaYU+RkISgCponIJBGJAa6nz/R/IpIVmGScwPypHvzTDBpjjBkhjl0aUtVuEbkDeBHwAver6k4RWRVYfy/wGeCrItINtAHXq02QYIwxI8rR0UcDl3vW9Fl2b6/HvwF+42QGY4wxA7M7i40xJsJZITDGmAhnhcAYYyKcFQJjjIlwVgiMMSbCWSEwxpgIZ4XAGGMinBUCY4yJcFYIjDEmwlkhMMaYCGeFwBhjIpwVAmOMiXBWCIwxJsJZITDGmAhnhcAYYyKcFQJjjIlwVgiMMSbCWSEwxpgIZ4XAGGMinBUCY4yJcFYIjDEmwkW5HcAYE7xun4/a4x0ca+0iPtpLVkoccdFet2OZMGeFwJgw0NLRzdt763hvfwNtXT0fLI/yCHNyUvjE7PGkJsS4mNCEMysExoS4kuomnt58iNaObs6YkMwZOSmkJcTQ2tnN7sPH2Vx+jJLqJq5akMOCvFS345owZIXAmBClqrz1fh0v7DxMdkocXzyngOyU+I9sMyMrmfOmZfDkxgoeL66gvauHpZPHuZTYhCtrLDYmRL1ScoQXdh5mbk4KX10+5WNF4IS0xBhuO2cSM7OSWL21iq0Vx0Y2qAl7VgiMCUHv7K1jbWkthRPH8rmz8ojyDvynGu31cOOSfCaOS+DpzZUcbmwfoaRmNLBCYEyI2V3dxHPbqzljQjJXL8zBIxLU90V5PNywOJ+4aC+PFpXT3eNzOKkZLawQGBNC6po7eHxjBRNS47iuMC/oInBCclw01y7MpfZ4B2/sqXUopRltrBAYEyJ6fMqfiyoQhJuWTCR6kMtBJzMjK4l5uSm8XlpL7fGOYU5pRiMrBMaEiDf21HLoWBtXL8xh7GneE3DFvAlEeYUXdx4epnRmNLNCYEwIqG5sY+3uGubmpDA3J+W09zcmNorzpqWzq7qJ8obWYUhoRjMrBMa4rLPbxxPFlcTHeLlq/oRh2+85U9MZExvFS3ZWYAbhaCEQkctFpFRE9orInQNsd5aI9IjIZ5zMY0wo+u3avRxuaueahTkkxA7fPZ6xUV7On55BWV0LFXZWYAbgWCEQES/wW2AFMBu4QURmn2S7nwEvOpXFmFB1sL6F372xj3m5KczKTh72/Z9VMJb4aK/1IDIDcvKMYDGwV1XLVLUTeAy4qp/tvg48BdQ4mMWYkPTjZ3cR7RFWzsl2ZP+xUV6WTk6jpLrJehCZk3KyEOQAFb2eVwaWfUBEcoBrgHsH2pGI3C4ixSJSXFtrn2zM6LC2tIZXSmr4+sXTSI6Pduw4y6ak4/EI68rqHTuGCW9OFoL+7oTRPs9/AXxbVXv62fbDb1K9T1ULVbUwIyNjuPIZ45qO7h7+7W+7mJyeyBfPmeToscbERjE3J4XN5Ufp6B7wT81EKCcLQSWQ1+t5LlDVZ5tC4DEROQB8BrhHRK52MJMxIeH+tw+wv66FH1x5BjFRznfeWzopjY5uH1srGh0/lgk/Tv4GFgHTRGSSiMQA1wOre2+gqpNUtUBVC4Anga+p6jMOZjLGdQ0tndyzdi+XzMpk+fSROcPNS0sgOyWODfvrUe17Ym4inWOFQFW7gTvw9wYqAR5X1Z0iskpEVjl1XGNC3a9fe5+Wzm7uXDFzxI4pIpxVkEZ1YzvVNjKp6cPRiWlUdQ2wps+yfhuGVfULTmYxJhSU17fy8PqDXFeYx9TMpBE99vzcVJ7bXs3G8qMjelwT+uzOYmNG0N0vleL1CP946fQRP3Z8jJfZ2clsrThGZ7cNUW0+ZIXAmBGyvbKR1Vur+LtzJzE+Oc6VDIvyx9La2cNru+22HfMhKwTGjABV5acvlDA2IZqvLJ/iWo6pmWNIjPHy7La+HfhMJLNCYMwIeOv9Ot7ZW88dF00jOc65m8cG4/UIZ+Sk8GpJDW2ddk+B8bNCYIzDfD7lp8/vJndsPDcvzXc7DnNzUmjr6mFtqV0eMn5WCIxx2F+3HmJXdRP/fNkMYqO8bsdhUnoi6WNieG5btdtRTIiwQmCMg9q7erj7xT3MyUnmU/OGb66B0+ERYcWcbF7dfYTWzm6345gQYIXAGAc9vP4gh461cefls/B4hjYRvZNWzs2mvctnvYcMYIXAGMc0tnXxm7V7OW9aOudOS3c7zkcsnpRG+phYuzxkACsExjjmd6/v41hrF9++fOSGkgiW1yOsnJvFa7traOmwy0ORzgqBMQ6obmzjj+/s5+oFE5gzDJPRO+GTc7Pp6Pbxql0einhWCIxxwN0v7kEVvvmJGW5HOanCgjTSx8Tw8q4jbkcxLrNCYMww23Gokac3V3LbOQXkpSW4HeekvB7hopmZvF5aQ1ePjT0UyawQGDOMVJWfPFdCanw0X7twqttxBnXJrPEcb++maH+D21GMi6wQGDOMXi2pYV1ZPf9wyXRSHJyHeLicOy2dmCgPL5fY5aFIZoXAmGHS1ePjP54vYXJGIjcucX8oiWAkxERx7tR0Xik5YjOXRTArBMYMk8feK6estoXvrJhFtDd8/rQunpVJRUMb79c0ux3FuMTRGcqMCVePbCgf0vbtXT3c/VIpk9ITqWlqH/L3u+nimeP5Hjt4peQI08eP7KxpJjSEz8cWY0LY66W1tHX2sHJuNiKhM5REMLJS4pibk8Ir1o00YlkhMOY01Td38M6+OhbkpZKTGu92nFNyyazxbK44Rl1zh9tRjAusEBhzGlSVv22rIsojXDYny+04p+yS2ZmowmsldpdxJLJCYMxpKKk+zp4jzVw8a7yrM4+drtnZyWSnxNlopBHKCoExp6iz28ez26vITIpl2eRxbsc5LSLCBTMyeXtvHZ3ddpdxpLFCYMwpemNPLcdau7hywQS8ITTXwKm6cEYGzR3dFB+0u4wjjRUCY05B7fEO3ny/lvm5KUxOH+N2nGFxztR0Yrwe1trloYhjhcCYIfKp8vTmSmK8HlbOzXY7zrBJjI1iyeQ01pbWuh3FjDArBMYMUdGBBg7Wt7JybhZJYdxA3J8LZmSyt6aZioZWt6OYEWSFwJghaGzr4oUdh5mckcii/LFuxxl2F83MBGBtqV0eiiRBFQIRmeN0EGNCnaqyemsVPlWuWZATdncQB2NSeiIF4xKsnSDCBHtGcK+IvCciXxORVCcDGROqtlQco6S6iUtmjWfcmFi34zjmwpmZvLuvnrbOHrejmBESVCFQ1XOBm4A8oFhEHhGRSx1NZkwIOdbayeqtVUxMS+Ccqelux3HUhTMy6ej2sb6s3u0oZoQE3Uagqu8D3we+DSwHfiUiu0XkWqfCGRMKfKo8uakSVfjMmbl4RuElod4WT0ojPtprdxlHkGDbCOaJyM+BEuAi4FOqOivw+OcDfN/lIlIqIntF5M5+1l8lIttEZIuIFIvIuaf4OoxxzPqyespqW/jk3OxRfUnohLhoL+dMTWdtaY1NVhMhgj0j+A2wCZivqn+vqpsAVLUK/1nCx4iIF/gtsAKYDdwgIrP7bPZqYJ8LgC8CfxjyKzDGQTVN7byw4zAzxidRWDD6egmdzIUzM6g82sa+WpusJhIEWwhWAo+oahuAiHhEJAFAVR86yfcsBvaqapmqdgKPAVf13kBVm/XDjxyJgH38MCGjx6c8sbGSmCgP1y4anb2ETubCGf5upHZ5KDIEWwheAXoPtJ4QWDaQHKCi1/PKwLKPEJFrRGQ38Bz+swJjQsLa0hoOHWvj6gU5o+7GscFMSI1nZlYSa3fbXcaRINhCEKeqH5wjBh4nDPI9/X18+tgnflX9i6rOBK4GftzvjkRuD7QhFNfW2i+mcV55fQuvl9awMC+VOTkpbsdxxQUzMik60EBTe5fbUYzDgi0ELSKy6MQTETkTaBvkeyrxdzc9IReoOtnGqvomMEVEPtY3T1XvU9VCVS3MyMgIMrIxp6a5o5vHN1aSEh/Np+ZPcDuOay6amUm3T3nn/Tq3oxiHBVsI/gF4QkTeEpG3gD8DdwzyPUXANBGZJCIxwPXA6t4biMhUCVx4DRSaGMA6LxtX/Wj1To62dHJdYR5x0V6347hmUX4qyXFRNtxEBIgKZiNVLRKRmcAM/Jd8dqvqgOeLqtotIncALwJe4H5V3SkiqwLr7wU+DdwqIl34zzA+p9Zfzbjo+e3VPLGxkgtnZDBxXKLbcVwV5fVw/vQM1pbW4vMpnlEw54LpX1CFIOAsoCDwPQtFBFV9cKBvUNU1wJo+y+7t9fhnwM+GkMEYxxxubOc7f9nO/NwULpo53u04IeHCGZk8u62aXdVNEdtWEgmCvaHsIeBu4Fz8BeEsoNDBXMaMKJ9P+dYTW+no8vHzzy0YFTOODYflMzIQsW6ko12wZwSFwGy7bGNGq/vf2c/be+v4z2vnMjljDOvLbLpGgPQxsczLTWVtaQ3fuHia23GMQ4ItBDuALKDawSzGnNQjG8od23d1Yxv3vL6PWdnJ+Hzq6LHC0YUzMvjlq+/T0NJJWmKM23GMA4LtNZQO7BKRF0Vk9YkvJ4MZMxK6fT6eKK4kIdrLNQsj6+7hYF00MxNVeGOPXR4arYI9I/ihkyGMccsbpbUcbmrnlqUTGRM7lL4TkWPOhBTSx8Swdnct1yzMdTuOcUCw3UffEJGJwDRVfSUwzlDkdrA2o0J1YxtrS2uYn5vCrOxkt+OELI9HWD49k1dKjtDjU2tIH4WC7TX0ZeBJ4PeBRTnAMw5lMsZxPT7lqU2VxMdE8al5kXv3cLAumplJY1sXm8uPuh3FOCDYNoK/B84BmuCDSWoynQpljNPeer+WqmPtXDV/Agl2SWhQ505Lx+sRu8t4lAq2EHQEhpIGQESisCGjTZhqaOnktd01nDEh2W6SClJKfDRnThzLazYa6agUbCF4Q0S+C8QH5ip+Avibc7GMcc6z26rwiHCFXRIakotmZlJS3cThxna3o5hhFmwhuBOoBbYDX8E/bES/M5MZE8pKqpvYffg4F8/KJCU+suYYOF0nJqt53S4PjTrB9hryAf8T+DImLHX1+Hh2WxWZSbGcPeVjo52bQUwfP4ac1Hhe213D9Yvz3Y5jhlFQhUBE9tP/pDKThz2RMQ55Y08tR1u7+NK5k6wL5CkQES6YkcEzmw/R0d1DbJT1IB8tgr00VMiHg82dB/wKeNipUMYMt2Otnby5p5Z5uSlMzhjjdpywdeGMTFo6eyjab91IR5OgCoGq1vf6OqSqvwAucjaaMcPn1ZIaFLjsjCy3o4S1s6eOIybKY91IR5lgbyhb1OurMDC5TJLD2YwZFocb29lUfpRlk8cxNsEGTTsdCTFRLJ08zgrBKBPsnTT/p9fjbuAAcN2wpzHGAS/uPExstIcLZth818PhohkZ/PBvuzhQ10JBemTP4jZaBHtp6MJeX5eq6pdVtdTpcMacrrLaZkqPHOeC6ZkkxNgdxMPh4ln+2dte2nXY5SRmuATba+ifBlqvqv89PHGMGT6qygs7D5MSH82yKePcjjNq5KUlMDcnhTXbD3P7+VPcjmOGwVB6DX0V/2BzOcAqYDb+dgJrKzAhac+R41QebePimZlEe4P9VTfBuHxOFlsqjlF1rM3tKGYYDGVimkWq+k1V/SZwJpCrqj9S1R85F8+YU6OqvLa7hrEJ0SzMH+t2nFFnxRx/76sXdtjlodEg2EKQD3T2et4JFAx7GmOGyd6aZiqOtrF8eqbdPOaAyRljmJmVxPM7bPba0SDYQvAQ8J6I/FBEfgBsAB50LpYxp05VeXV3DSnx0SzKT3U7zqi1Yk42xQePUtNkg9CFu2B7Df0EuA04ChwDblPV/3AwlzGnrKyuhfKGVpZPzyDK2gYcs3JuFqr+7rkmvA3lryQBaFLVXwKVIjLJoUzGnJbXdteQHBfFmROtbcBJ08YnMSUjkTXbrRCEu2DvLP4B8G3gO4FF0dhYQyYElTe0sr+uhfOmZVhPoRGwcm42G/bXU9/c4XYUcxqC/Uu5BrgSaAFQ1Sqs26gJQW+/X0tctIfCAjsbGAkr5mTjU1iz3RqNw1mwhaBTVZXAUNQiYveVm5DT0NLJzqomFheMsyGSR8is7CRmjE/i6c2H3I5iTkOwheBxEfk9kCoiXwZewSapMSHmnX11iGB3EY8gEeHaRTlsLj/G/roWt+OYUzRoIRARAf4MPAk8BcwA7lLVXzuczZigtXX2sPHAUebnptoUlCPsqgU5iMBfNlW6HcWcokHHGlJVFZFnVPVM4OURyGTMkL23v57OHh/nTrMpKEdaVkoc50xJ5y9bDvGPl07H/9nRhJNgLw2tF5GzHE1izCnq9vlYV1bP1IwxZKfEux0nIl27KIeKhjaKD9rMZeEo2EJwIf5isE9EtonIdhHZ5mQwY4K1s6qJpvZuzp5qbQNuueyMLOKjvTxtl4fC0oCFQETyAw9XAJPxT0/5KeCKwL8DEpHLRaRURPaKyJ39rL8pUFi2ici7IjJ/6C/BRLr1ZfWkJcYwfbz1aHZLYmwUK+Zk8ey2atq7etyOY4ZosDOCZwBU9SDw36p6sPfXQN8oIl7gt/iLyGzgBhGZ3Wez/cByVZ0H/Bi47xReg4lg1Y1tHKxvZcmkNDx2bdpV1yzK4Xh7Ny/tOuJ2FDNEgxWC3n9Zk4e478XAXlUtU9VO4DHgqt4bqOq7qnriouJ6IHeIxzARbn1ZPdFeseEkQsDZU9LJS4vn4fUDfkY0IWiwQqAneRyMHKCi1/PKwLKT+Tvg+f5WiMjtIlIsIsW1tbVDjGFGq7bOHrZUHGN+bqpNQxkCvB7h5iUTeW9/A7sPN7kdxwzBYIVgvog0ichxYF7gcZOIHBeRwf6n+ztP77eYiMiF+AvBt/tbr6r3qWqhqhZmZNgE5MZvU/lRunqUpZOtkThUXFeYR2yUhwfX2VlBOBmwEKiqV1WTVTVJVaMCj088Tx5k35VAXq/nuUBV341EZB7wB+AqVa0f6gswkcmnyvqyevLTEpiQal1GQ8XYxBiunD+BZzYfoqm9y+04JkhODs9YBEwTkUkiEgNcD6zuvUGgV9LTwC2qusfBLGaUKattob6lk6WT09yOYvq4dVkBrZ09PLXRupKGC8cKgap2A3cALwIlwOOqulNEVonIqsBmdwHjgHtEZIuIFDuVx4wuRQcaiI/2csaEFLejmD7m5qawIC+Vh9YdxOcbatOicYOjLWyqugZY02fZvb0efwn4kpMZzOjT3NHNrqomlkxOszkHQtTnz57IP/55K2/vreP86Rk8sqF8RI9/45L8wTcyH7C/IhN2NpcfpUeVswrsslCoWjk3m8ykWO55fa/bUUwQrBCYsKKqFB04Sn5aAuOT49yOY04iNsrLV5ZPYX1ZA+/tb3A7jhmEFQITVg7Ut1LX3GFnA2HgxsX5pI+J4devve92FDMIKwQmrBQfaCA2ysPcHGskDnXxMV6+fN5k3nq/jvKGVrfjmAFYITBho62zh+2HGlmQl0pMlP3qhoObl05kbEI0a3fXuB3FDMD+mkzY2FJxlG6fNRKHk8TYKL503mRKjxzn0NE2t+OYk7BCYMLCiUbinNR4u5M4zNy6bCLx0V6e31mNqt1XEIqsEJiwUHm0jcNN7RQW2Cij4SYpLppLZo+nrLaFkmobjC4UWSEwYaHoQAPRXmF+bqrbUcwpWFyQRmZSLGt2HKa7x+d2HNOHFQIT8po7utlW2ci83FTior1uxzGnwOsRPjkvm4aWTt7ZZ2NLhhorBCbkrd5SRWePzxqJw9y0zCRmZiWxtrTGRiYNMVYITMh7rKic8cmx5I21RuJwt3JuNj09yprt1W5HMb3YtE4mpO2samRbZSNXzMtGbE7iYTPSg8CdkD4mlgtnZvBKSQ1zJjQyx24MDAl2RmBC2mPvVRAT5WFBXqrbUcwwWT49kwmpcfx1yyGaO7rdjmOwQmBCWFtnD89sOcTKOVk2J/Eo4vUIn1mUR3uXj9VbPzZpoXGBFQITsp7bXs3x9m6uX2xjy482WSlxXDwrkx2HGtlaecztOBHPCoEJWY+9V87k9ESWTLLeQqPRedMyyBsbzzObD1Hf3OF2nIhmhcCEpD1HjlN88CjXL86zRuJRyusRrl+cj0eER94rp8tuNHONFQITkh7ZUE6M18NnzsxzO4px0NiEGK4rzKW6sZ2/WXuBa6wQmJDT3tXD05squWxOFmmJMW7HMQ6bkZXMBTMyKD54lI0HbTYzN1ghMCFnzfZqmtq7uWGxnQ1EiktmjWdKRiLPbKniYH2L23EijhUCE3Ie2VDOpPRElk0e53YUM0I8ItywOJ/U+GgeWn+QhpZOtyNFFCsEJqScaCS+wRqJI05CTBSfX1aAKjy47gDtXT1uR4oYVghMSLFG4siWnhTLjUvyqWvu4LGicnp8NpHNSLBCYEKGNRIbgCkZY7hqQQ57jjTb4HQjxO7bNyHjRCPxjXYnccQ7qyCN2uMdvL23jvQxMSybku52pFHNCoEJGScaiZdOtjuJDVw+J4v65g6e3VZNWmIMM7KS3Y40atmlIRMSrJHY9OUR4bqz8shOiePRogqqG9vcjjRqWSEwIeHBdQeIibJGYvNRsVFebllWQFyUhwfXHaSpzWY2c4IVAuO6pvYunt50iCvnT7BGYvMxKfHR3LqsgLbOHh5af5DObhuTaLhZITCue6K4ktbOHj6/rMDtKCZETUiN5/qz8qg61sbjxRX41LqVDicrBMZVPp/y0LoDLMpPZW6uTVtoTm5mdjKfnJfNruomXthx2O04o4qjhUBELheRUhHZKyJ39rN+poisE5EOEfmWk1lMaHpjTy0H6lv5/NkFbkcxYeDsKeksnTyOt/fWsWF/vdtxRg3Huo+KiBf4LXApUAkUichqVd3Va7MG4BvA1U7lMKHtgXcPkJEUy4o52W5HMWHik3OzOdrSyd+2VjE2IYbp45PcjhT2nDwjWAzsVdUyVe0EHgOu6r2BqtaoahFgXQEi0P66Ft7YU8tNS/KJibKrlCY4Xo9w/Vl5ZCbF8eh75dQdt9nNTpeTf305QEWv55WBZUMmIreLSLGIFNfW1g5LOOO++9/eT4zXw41L7E5iMzSx0V5uWTYRr0d4aMNBOmyAutPiZCHo766gU2rqV9X7VLVQVQszMjJOM5YJBfXNHTxeXME1C3PITIpzO44JQ2MTYrhhcT71zR08sbHSehKdBicLQSXQ++6gXMDmojMA/L91B+no9vHl8ye7HcWEsSkZY7h8jr8n0Rt77GrBqXKyEBQB00RkkojEANcDqx08ngkTrZ3dPLjuAJfOHs/UzDFuxzFh7pwp41iQl8oru45QerjJ7ThhybFCoKrdwB3Ai0AJ8Liq7hSRVSKyCkBEskSkEvgn4PsiUikiNrLUKPd4UQXHWrtYtdzOBszpExGuXpBDVkocfy6uoK7ZGo+HytGuGqq6RlWnq+oUVf1JYNm9qnpv4PFhVc1V1WRVTQ08tpI+inX3+Pift/ZTOHEsZ060UUbN8IiJ8nDzkol4RPjThoO0dna7HSmsWJ89M6Ke217NoWNtfGX5FLejmFFmbGIMnyvMo6apg+88vR21xuOgWSEwI6a7x8cvX3mfGeOTuHhmpttxzCg0bXwSl8wez1+3VPHguoNuxwkbNjGNGTHPbKmirK6Fe28+E4/H5hwwzlg+PQNV5cfP7mJOTrJdggyCnRGYEdHV4+OXr+5hTk4yl50x3u04ZhTziPB/rltAzth4vvanTdTanceDskJgRsQTxZVUNLTxzUtn2AxkxnEp8dH87qYzaWzr4uuPbqK7x+YwGIgVAuO49q4efv3a+yzMT+WCGXZnuBkZsyck85Or57K+rIH/erHU7TghzQqBcdwjG8qpbmy3swEz4j59Zi43L83n92+W8fz2arfjhCwrBMZR9c0d/OKVPZw7NZ1zpo5zO46JQP96xWwW5KXyz09uY19ts9txQpIVAuOou18qpbWzhx9eOdvOBowrYqO83HPTImKiPKx6aCMtHXazWV9WCIxjtlc28lhRBV84u4CpmTZ5iHHPhNR4fn3DQvbVNvPtp7bZzWZ92H0E5pQ8sqF8wPU+Ve57s4zEmCgmpMYPur0xTjtnajr/fNlMfvbCbhblj+WL505yO1LIsDMC44gt5ccob2jlsjOyiIv2uh3HGABWLZ/MJ2aP5z/WlFB0oMHtOCHDCoEZdsfbu3huezX5aQkszE91O44xHxAR7r5uPnlpCXztT5uoaWp3O1JIsEJghpWq8tctVXT1+Pj0olw81kBsQkxyXDT33nwmze3d3PHIZrrsZjMrBGZ4FR88yq7qJi6dPZ6MpFi34xjTrxlZSfz003N570ADP3t+t9txXGeFwAybI03tPLutiikZiZwzNd3tOMYM6KoFOXzh7AL+8PZ+nt0W2bPoWiEww6K9q4dHNpQT4/Xw2cI8uyRkwsJ3V87izIlj+Zcnt1FSHblzYlkhMKfNp8rjxRXUt3Rww+J8kuOi3Y5kTFBiojzcc9MikuOi+cIf3+PQsTa3I7nCCoE5LarKs9uq2H34OJ+cm83kDJuM3oSX8clxPPDFs2jt6OEL979HY2uX25FGnBUCc1pe3V3D+rIGzpuazrIp1i5gwtPMrGR+f+uZHKxv5csPFtPe1eN2pBFlhcCcElXl5V2HeW13DYvyx3LZnCy3IxlzWs6eks7d183nvQMNrHp4Ix3dkVMMrBCYIevq8fGXzYdYW1pL4cSxXLsoxxqHzahw5fwJ/Oe1c3m9tJavPrwpYoqBFQIzJNWNbdz0hw0UHzzKhTMyuHqhFQEzutywOJ//uGYur+2u4WsRUgxs0LlRwulB3VSVTeXHWLO9mh6fcl1hHgvyUh09pjFuuXFJPj5Vvv/MDm77YxH33nLmqO4NZ4XADEhV2VfbwislRyhvaCU/LYHPnJlL+hi7a9iMbjcvnUhCjJd/eXIb1927jj/edhbZKfFux3KEFQLTr8a2LnZWNbLx4FGqG9tJiY/m2oU5LJo41i4FmYhx7aJcMpPiWPXwRq757bv84fOFzMlJcTvWsLNCEKa6e3zUt3TS1NZFU3s37x85Tt+pNkQgyuPB6xGiPPKRf70eocendHb76Ozxcay1i/qWTuqbOyhvaKXmeAcA2SlxXL0gh4X5qUR7rUnJRJ5zp6XzxKplfPGBIq69512+f8Usblk6cVTNuGeFIMQ1d3Szu7qJkuomdlUfZ2/NcaqOtXO4qZ0e3/DPspQY42VCajwL81KZlZ1MZnLcsB/DmHAzKzuZ575xHt96Yit3/XUn6/bV89Nr55GSMDraDawQhJi65g7W7atnXVk96/fVU1bX8sG65LgoZmYls2RSGhNS4xmfEkdqfDRJcVGs21dP388nPoUeVbp7lB6fj26f0uPTD/6N8gjRUR5ivB5S4qNJS4yxSWSMOYm0xBj+cGshf3i7jP/9QilFB17nuytncc3CnLA/O7BC4LLG1i7W76/3v/nvq6f0yHEAkmKjWDwpjWsX5TArO5lZ2clkp8Sd9Beu6phNsGGM0zwe4fbzp3D2lHS+/8wO/unxrTxeXMFdV5zB7AnJbsc7ZVYIRlhLRzdFBxpYt6+ed/fVs6OqEVWIi/ZwVkEaVy/MYdmUccyZkEyUXZM3JiTNyUnh6a+ezaNF5fzs+d2s/NVbXDp7PN+4aBpzc8OvMdkKgcMa27ooPtDAe/sbWL+/gR2HGunxKTFeDwvzU/lfF0/j7CnpzM9LITbKLssYEy48HuGmJRO5Yu4E/vjufu5/ez+f2nWEJZPS+GxhHivnZpEQEx5vsaI6/A2OTiosLNTi4mK3Y/Sru8fH+zXN7DjUyI5DjR/M1qUKMV4P8/NSWDJpHEsmp1E4MY34mOF743f6hjJjwsmNS/JH/JjH27t4eH05fy4q50B9K4kxXi6ZPZ4LZmRw/rQMxrl8742IbFTVwv7WOVquRORy4JeAF/iDqv60z3oJrF8JtAJfUNVNTmYaDs0d3Rysb+FAXSv765rZX9fK3tpmdlc30dHtn/80IcbL/Fz/J/7Fk9JYlD/WGmKNGcWS4qL56gVTWLV8MkUHjvLkxgpeLanhr1uqEPGPcDo/N4X5eanMzk5mUkZiyNyt7FghEBEv8FvgUqASKBKR1aq6q9dmK4Bpga8lwO8C/w67Hp/S3tVDV4+Pzm4fHd0+/+PA864eH62dPRxv76aprYvj7d0cb/f30T/W2smRpg6OHG+npqmD5o7uj+w7KzmOyRmJ3LpsInNyUpiTk0LBuES8nvDuSWCMGToRYfGkNBZPSsPnU7YfauT10lo2lh/l+R2Heayo4oNt08fEkDs2gYykWP/XmFjSk2JJjosiPtpLYmwUCTFeEmL8/45NjGFM7PC/bTt5RrAY2KuqZQAi8hhwFdC7EFwFPKj+61PrRSRVRLJVtXq4w6zZXs3XH908pO8RgTGxUaQmRDM+KY5ZWcksnx5LZlIc+WkJTEpPpCA9IWyuAxpjRpbHI8zPS2V+YFwuVaW8oZXSw8fZX9dCWW0LVY1tVDS0sungURpaOxnoav1Xzp/Md1bOGvacTr6D5QAVvZ5X8vFP+/1tkwN8pBCIyO3A7YGnzSJSGmSGdKAu2MAhIJzyhlNWsLxOC6m8Nw28OqSyBuGDvN/9GXz31Pcz8WQrnCwE/V0X+dgoCEFsg6reB9w35AAixSdrHAlF4ZQ3nLKC5XVaOOUNp6wwMnmd7KheCeT1ep4LVJ3CNsYYYxzkZCEoAqaJyCQRiQGuB1b32WY1cKv4LQUanWgfMMYYc3KOXRpS1W4RuQN4EX/30ftVdaeIrAqsvxdYg7/r6F783UdvG+YYQ76c5LJwyhtOWcHyOi2c8oZTVhiBvGF3Q5kxxpjhZYPZGGNMhLNCYIwxES5iCoGIfEtEVETS3c5yMiLyYxHZJiJbROQlEZngdqaBiMh/icjuQOa/iEiq25kGIiKfFZGdIuITkZDsPigil4tIqYjsFZE73c4zGBG5X0RqRGSH21kGIyJ5IrJWREoCvwf/y+1MAxGROBF5T0S2BvL+yKljRUQhEJE8/ENdhPrIbP+lqvNUdQHwLHCXy3kG8zIwR1XnAXuA77icZzA7gGuBN90O0p9ew7KsAGYDN4jIbHdTDeoB4HK3QwSpG/imqs4ClgJ/H+I/3w7gIlWdDywALg/0rhx2EVEIgJ8D/0I/N6uFElVt6vU0kdDP+5Kqnhh4aT3++0BClqqWqGqwd6W74YNhWVS1EzgxLEvIUtU3gQa3cwRDVatPDGqpqseBEvwjGYQk9WsOPI0OfDnynjDqC4GIXAkcUtWtbmcJhoj8REQq8N8lH+pnBL19EXje7RBh7mRDrphhJiIFwEJgg8tRBiQiXhHZAtQAL6uqI3lHxWhpIvIKkNXPqu/hH5rjEyOb6OQGyqqqf1XV7wHfE5HvAHcAPxjRgH0Mljewzffwn3b/aSSz9SeYvCEsqCFXzOkRkTHAU8A/9DkLDzmq2gMsCLS//UVE5qjqsLfHjIpCoKqX9LdcROYCk4Ctgbl+c4FNIrJYVQ+PYMQPnCxrPx4BnsPlQjBYXhH5PHAFcLGGwE0pQ/j5hiIbcsVhIhKNvwj8SVWfdjtPsFT1mIi8jr89ZtgLwai+NKSq21U1U1ULVLUA/x/aIreKwGBEZFqvp1cCu93KEozAxEPfBq5U1Va384wCwQzLYk5RYCKs/wuUqOp/u51nMCKScaInnojEA5fg0HvCqC4EYeinIrJDRLbhv5wV0t3bgN8AScDLgS6v97odaCAico2IVALLgOdE5EW3M/UWaHg/MSxLCfC4qu50N9XARORRYB0wQ0QqReTv3M40gHOAW4CLAr+vW0RkpduhBpANrA28HxThbyN41okD2RATxhgT4eyMwBhjIpwVAmOMiXBWCIwxJsJZITDGmAhnhcAYYyKcFQITsUSkJ9CFcIeI/O1URk8VkUIR+dVJ1h041dFuReSHIvKtU/leY4bKCoGJZG2qukBV5+AfOO3vh7oDVS1W1W8MfzRjRo4VAmP81hEY4E1EpojICyKyUUTeEpGZgeWfDZw9bBWRNwPLLhCRZwOPxwXmkdgsIr8nMHaQiBT0Hq8/MDfGDwOPvywiRYF9PiUiCSP6qo3BCoExJ+YBuJgPh3O4D/i6qp4JfAu4J7D8LuCywPjwV/azqx8Ab6vqwsC+8oM4/NOqelZgnyVAKN+Za0apUTHonDGnKD4wxG8BsBH/UBljgLOBJwIDFQLEBv59B3hARB4H+huw7Hz8E9+gqs+JyNEgMswRkX8HUoEx+IeXMGZE2RmBiWRtgdngJgIx+NsIPMCxQNvBia9ZAKq6Cvg+/hFCt4jIuH722d+YLd189G8trtfjB4A7VHUu8KM+64wZEVYITMRT1UbgG/gvA7UB+0Xks+AfsVJE5gceT1HVDap6F1DHR4eMBv8UmDcFtl0BjA0sPwJkBtoQYvEP231CElAdGB75JkdeoDGDsEJgDKCqm4Gt+Id+vgn4OxHZCuzkw+ki/0tEtgcaft8MbN/bj4DzRWQT/tFjywP77gL+Df9sWM/y0aGE/zWw/GVCfNhxM3rZ6KPGGBPh7IzAGGMinBUCY4yJcFYIjDEmwlkhMMaYCGeFwBhjIpwVAmOMiXBWCIwxJsL9f5MMdRuMhfCvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot a histogram of the residuals\n",
    "sns.distplot(results.resid, hist=True)\n",
    "plt.xlabel('Residual')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Residual Histogram')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Your Turn 1\n",
    "Create a regression model from synthetic data with intercept of 0 and maximum value at ${x = 10, y = 10}$, and with a the error having a standard deviation of 5. Plot the result of your model. How does this slope and intercept of this model compare to the model from the data with a standard deviation of 1? \n",
    "\n",
    "**Hint:** You need need to add columns named `score` and `resids` to the data frame before you create the plots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# New Regression Model sd=5\n",
    "\n",
    "# Replace the ??? with the correct parameters\n",
    "# Paramters of generated data\n",
    "n_points = 50\n",
    "x_start, x_end = 0, ??\n",
    "y_start, y_end = 0, ??\n",
    "y_sd = ??\n",
    "\n",
    "# Generate data columns\n",
    "x_data = np.linspace(x_start, x_end, n_points)\n",
    "y_error = np.random.normal(loc=0, scale=y_sd, size=n_points)\n",
    "y_data = np.linspace(y_start, y_end, n_points) + y_error\n",
    "\n",
    "# Put data in dataframe\n",
    "reg_data_5 = pd.DataFrame({'x':x_data, 'y':y_data})\n",
    "\n",
    "reg_data_5.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace the ??? with the correct parameters\n",
    "# Fit the model and get the linear model summaries/plots.\n",
    "x5_input = reg_data_5['x']\n",
    "y5_output = reg_data_5['y']\n",
    "ols5_model = sm.ols(formula='y ~ x', data = reg_data_5)\n",
    "\n",
    "results5 = ols5_model.fit()\n",
    "\n",
    "# Get slope (m) and y-intercept (b)\n",
    "results5.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace the ??? with the correct parameters\n",
    "# Additional Statistics\n",
    "print('\\nSSE, SST, SSR, and RMSE:')\n",
    "mean5_y = np.mean(y5_output)\n",
    "sst5 = np.sum((y5_output - mean5_y)**2)\n",
    "sse5 = sst5 - results5.ssr\n",
    "print('SSE: {}'.format(sse5))\n",
    "print('SST: {}'.format(sst5))\n",
    "print('SSR: {}'.format(results5.ssr))\n",
    "print('RMSE: {}'.format(np.sqrt(results5.mse_model)))\n",
    "\n",
    "# Get most of the linear regression statistics we are interested in:\n",
    "print(results5.summary())\n",
    "\n",
    "# Plot a histogram of the residuals, using x to predict y\n",
    "y5_pred = ols5_model.fit().predict(x5_input)\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "sns.distplot(results5.resid, hist=True)\n",
    "plt.grid(True)\n",
    "plt.xlabel('Residual')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Residual Histogram')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.scatter(x5_input, y5_output)\n",
    "plt.plot(x5_input, y5_pred, linewidth=2)\n",
    "plt.grid(True)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('x vs y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace the ??? with the correct parameters\n",
    "\n",
    "# Compare slope and intercept\n",
    "print('sd=1 Intercept, Slope : \\n{}'.format(results.params))\n",
    "print('\\nsd=5 Intercept, Slope : \\n{}'.format(results5.params))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluation of regression models\n",
    "\n",
    "Now that you have built a regression model, let's look at how you can quantitatively evaluate the performance of a regression model. The evaluation of regression models is based on measurements of the errors.\n",
    "\n",
    "If we want to know how well we did, we must consider what a good benchmark would be.  A common easy benchmark is to see if we can predict better than just the mean of the target variable.\n",
    "\n",
    "If we define the following terms:\n",
    "\n",
    "$\\bar{Y}$ - The mean of the target values.\n",
    "\n",
    "$Y_{i}$ - Target value of the i-th point.\n",
    "\n",
    "$\\hat{Y_{i}}$ - Predicted value of the i-th point.\n",
    "\n",
    "Then we can defind the following error metrics:\n",
    "\n",
    "**Sum of Squared Residuals (SSR)**: This is the value that the method of least squares was able to minimize to with the best fit line.\n",
    "\n",
    "SSR = $\\sum \\left( Y_{i} - \\hat{Y_{i}} \\right)^{2}$\n",
    "\n",
    "**Sum of Squared Total (SST)**: This is a measure of how bad a prediction the target-mean would be.\n",
    "\n",
    "SST = $\\sum \\left( Y_{i} - \\bar{Y} \\right)^{2}$\n",
    "\n",
    "**Sum of Squared Error (SSE)**: This is a sum of a squared pointwise difference of the (SST - SSR).  You can also think of this as a measurement of how much better we are explaining the variation than the mean.\n",
    "\n",
    "SSE = $\\sum \\left( \\hat{Y_{i}} - \\bar{Y} \\right)^{2}$\n",
    "\n",
    "It is the case that **SST = SSR + SSE**.\n",
    "\n",
    "We also will be calculating something called the **Root-Mean-Squared-Error**, or the **RMSE**.  The reason we do this is because the Sum-squared-error (also known as the mean-squared error) is in units of x-squared.  To get a handle on understanding the error in units of 'x', we take the square root of the mean-squared-error , or the square root of the sum-square-error.  This is the RMSE.\n",
    "\n",
    "---------------\n",
    "\n",
    "The errors of a regression model can be visualized as shown in the figure below. \n",
    "\n",
    "<img src=\"https://library.startlearninglabs.uw.edu/DATASCI410/img/Errors.jpg\" alt=\"Regression_Errors\" style=\"width: 450px;\"/>\n",
    "\n",
    "<center>**Measuring errors for a regression model**\n",
    "$$Where\\\\\n",
    "Y = [y_1, y_2, \\ldots, y_n]\\\\\n",
    "and\\\\\n",
    "y_i = ith\\ data\\ value\\\\\n",
    "\\bar{Y} = mean(Y)\\\\\n",
    "\\\\\\hat{y_i} = regression\\ estimate\\ of\\ y_i\\\\\n",
    "SSE = sum\\ square\\ explained\\ = \\Sigma_i{(\\hat{y_i} - \\bar{Y})^2}\\\\\n",
    "SSR = sum\\ square\\ residual\\ = \\Sigma_i{(y_i - \\hat{y_i})^2}\\\\\n",
    "SST = sum\\ square\\ total\\ = \\Sigma_i(y_i - \\bar{Y})^2$$\n",
    "\n",
    "The goal of regression is to minimize the residual error, $SSR$. Specifically we wish to explain the maximum amount of the variance in the original data as possible with our model. We can quantify this idea with coefficient of determination also known as $R^2$.\n",
    "\n",
    "$$R^2 = 1 - \\frac{SSR}{SST}\\\\\n",
    "so\\ as\\\\\n",
    "SSR \\rightarrow 0\\\\\n",
    "R^2 \\rightarrow 1$$\n",
    "\n",
    "In words, $R^2$ is the fraction of the variance of the original data explained by the model. A model that perfectly explains the data has $R^2 = 1$. A model which does not explain the data at all has $R^2 = 0$.\n",
    "\n",
    "However, there are two problems with $R^2$. </center>\n",
    " - $R^2$ is not bias adjusted for degrees of freedom.\n",
    " - More importantly, there is no adjustment for the number of model parameters. As the number of model parameters increases $SSR$ will generally decrease. Without an adjustment you will get a false sense of model performance.\n",
    " \n",
    "To addresses these related issues, we can use adjusted $R^2$.\n",
    "\n",
    "$$R^2_{adj} = 1 - \\frac{\\frac{SSR}{df_{SSR}}}{\\frac{SST}{df_{SST}}} = 1 - \\frac{var_{residual}}{var_{total}}\\\\\n",
    "where\\\\\n",
    "df_{SSR} = n - k = SSR\\ degrees\\ of\\ freedom\\\\\n",
    "df_{SST} = n - 1 = SST\\ degrees\\ of\\ freedom$$\n",
    "\n",
    "This gives $R^2_{adj}$ as:\n",
    "\n",
    "$$R^2_{adj} = 1 - (1 - R^2) \\frac{n - 1}{n - k}\\\\ \n",
    "where\\\\\n",
    "n = number\\ of\\ data\\ samples\\\\\n",
    "k = number\\ of\\ model\\ coefficients$$\n",
    "\n",
    "Or, we can rewrite $R^2_{adj}$ as:\n",
    "\n",
    "$$R^2_{adj} =  1.0 - \\frac{SSR}{SST}  \\frac{n - 1}{n - 1 - k}$$\n",
    "\n",
    "Another measure of regression performance is root mean square error or $RMSE$:\n",
    "\n",
    "$$RMSE = \\sqrt{ \\frac{\\Sigma^n_{i-1} (y_i - \\hat{y_i})^2}{n}} = \\frac{\\sqrt{SSR}}{n}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your Turn 2\n",
    "\n",
    "Examine the performance metrics for the previous two regressions with sd=1 and sd=5. \n",
    "- How do SSE, SSR, SST, $R^2$, and RMSE compare?\n",
    "\n",
    "**Hint**: Use `results.rsquared` to find the $R^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.rsquared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## fill in ?? \n",
    "print('      sd1                sd5')\n",
    "print('SSE: ', sse, sse5)\n",
    "print('SSR: ', results.ssr, results5.ssr)\n",
    "print('SST: ', sst, sst5)\n",
    "print('Rsq: ', results.??,results5.??)\n",
    "print('RMSE:',np.sqrt(results.mse_model), np.sqrt(results5.mse_model))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Leverage and Cook's Distance\n",
    "\n",
    "Up to now, we have only looked at regression models with Normally distributed noise or errors. But, in the real world there are errors and outliers in data. These errors and outliers can have greater or lesser effect, depending on how extreme they are and their placement with respect to the other data. \n",
    "\n",
    "You can imagine a regression line as a lever. Outliers that occur near the ends of the lever will have a greater influence all other factors being equal. \n",
    "\n",
    "One way to measure influence of a data point is Cook's distance, introduced by Dennis Cook in 1977. The influence for the `ith` data point can be computed as:\n",
    "\n",
    "$$D_i = \\frac{\\Sigma_{j=1}^n (\\hat{Y_j} - \\hat{Y_{j(i)}})^2}{n (p+1)\\hat{\\sigma^2}} \\\\\n",
    "where \\\\\n",
    "p = number\\ of\\ parameters\\\\\n",
    "n = number\\ of\\ data\\ points$$\n",
    "\n",
    "In effect, Cook's distance compares the difference between means with and without a given data point. Computing Cook's distance can be moderately computationally intensive for large data set. Typically, Cook's distance is measured in units of standard deviation.\n",
    "\n",
    "Let's make these concepts concrete with an example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a scale function for a column in a pandas df\n",
    "def scale(col):\n",
    "    mean_col = np.mean(col)\n",
    "    sd_col = np.std(col)\n",
    "    std = (col - mean_col) / sd_col\n",
    "    return std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create new copy of original linear data\n",
    "outlier_data = sim_data[['x', 'y']].copy()\n",
    "# Add an outlier row\n",
    "outlier = pd.DataFrame([[0.0, 20.0]], columns = ['x', 'y'])\n",
    "outlier_data.append(outlier, ignore_index=True)\n",
    "\n",
    "# Scale data\n",
    "outlier_data['x_scaled'] = scale(outlier_data['x'])\n",
    "\n",
    "# Fit linear model\n",
    "outlier_linear = sm.ols(formula = 'y ~ x_scaled', data = outlier_data).fit()\n",
    "\n",
    "# Get most of the linear regression statistics we are interested in:\n",
    "print(outlier_linear.summary())\n",
    "\n",
    "# Plot a histogram of the residuals\n",
    "y_outlier = outlier_linear.predict(outlier_data['x_scaled'])\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "sns.distplot(outlier_linear.resid, hist=True)\n",
    "plt.grid(True)\n",
    "plt.xlabel('Residual')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Residual Histogram (With 1 Outlier)')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.scatter(outlier_data['x_scaled'], outlier_data['y'])\n",
    "plt.plot(outlier_data['x_scaled'], y_outlier, linewidth=2)\n",
    "plt.grid(True)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title('x vs y (With 1 Outlier)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Linear Regression Summary Plots\n",
    "# Summary graphs:\n",
    "import scipy.stats as stats\n",
    "import statsmodels.api as statsmodels\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "\n",
    "# Residuals vs Fitted Values\n",
    "residuals = outlier_linear.resid # outlier_linear = name of linear model on our dataset\n",
    "fitted_vals = outlier_linear.predict(outlier_data['x_scaled']) # making predictions from our fit model\n",
    "plt.plot(fitted_vals, residuals, 'o') # plotting predictions from fit model vs residuals\n",
    "plt.xlabel('Fitted Values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.title('Residuals vs. Fitted Values (w/ 1 outlier)')\n",
    "\n",
    "# Create Q-Q Normal Plot of the Residuals\n",
    "statsmodels.qqplot(residuals, stats.norm, fit=True, line='45')\n",
    "\n",
    "# Fitted Values vs. Square Root of the Standardized Residuals\n",
    "standardized_resid = (outlier_linear.resid - np.min(outlier_linear.resid)) / np.max(outlier_linear.resid)\n",
    "\n",
    "# Leverage Plot (Cook's Distance)\n",
    "influence_plot(outlier_linear)  # from statsmodels.graphics.regressionplots import *\n",
    "\n",
    "# Additional stats models plots\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "fig = statsmodels.graphics.plot_regress_exog(outlier_linear, \"x_scaled\", fig=fig)\n",
    "\n",
    "'''\n",
    "Plot regression results against one regressor.\n",
    "\n",
    "This plots four graphs in a 2 by 2 figure: ‘endog versus exog’, ‘residuals versus exog’, ‘fitted versus exog’ and ‘fitted plus residual versus exog’\n",
    "\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "We looked at a lot of Regression concepts in this lesson so far:\n",
    "\n",
    "- Single regression equation of a line is $y_i = mx_i + b + \\epsilon_i $\n",
    "- Use the Python model object `linear_model.LinearRegression()` to initialize the model\n",
    "- Use the `fit` method to fit the model to the data\n",
    "- Use `predict` method to compute scores(predict values) for the dependent value $y$\n",
    "- Pull the model parameters $m$ (slope) with `coef` and $b$ (y-intercept) with `intercept`\n",
    "- The residuals are the difference between the y_output and the y_predicted\n",
    "\n",
    "- Use the `statsmodels` and `seaborn` libraries to get summary statistics and make diagnostic plots.\n",
    "- `sm.ols(formula = 'y ~ x', data=sim_data)` and `fit()` method to get the results\n",
    "\n",
    "Evaluate the results from the OLS Regression Results with `summary` method\n",
    "- SSR the Sum of the Squared Residuals is the difference in error from the regression line (try to minimize)\n",
    "- SST the Sum of the Squared Total is the difference to the target mean\n",
    "- SSE the Sum of the Squared Error is the variation of the regression line to the mean\n",
    "- RMSE Root Mean Squared Error is the square root of the SSE\n",
    "- $R^2$ is the fraction of the variance of the original data explained by the model, want a number close to 1.\n",
    "\n",
    "We created a function to scale data reconstructing the mean standardization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<div id=\"reminder\" style=\"border-radius: 5px; background-color:#f5f5f5;\" >\n",
    "<h3>Reminder</h3>\n",
    "<p>Use this notebook to answer the quiz questions related to the <b>Your Turn</b> sections.<p>\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}