Lab 02 - Linear Regression with MLE and MAP.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Mini-lab 02 - Linear Regression with MLE and MAP.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yDW9MxvfHY7i"
      },
      "source": [
        "# Lab 2: Linear Regression with MLE and MAP\n",
        "\n",
        "In this lab, you will experiment with a linear regression model. The data and optimization method can be changed. Try them out and see what happens. Which optimization methods are better? Why? Could they be adapted for different scenarios?\n",
        "\n",
        "As mentioned in the lecture, least square linear regression can often be done analytically. The main pedagogical point of this lab is to illustrate finding model parameters using interative optimization. Regression models are primarily used as illustration.\n",
        "\n",
        "There are several regression models where the parameters cannot be found analytically. Most of these models fall outside the scope of this course.\n",
        "\n",
        "As usual, we start with imports."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OXLcCtULdPeb"
      },
      "source": [
        "import numpy as np\n",
        "import time\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import ipywidgets\n",
        "# from IPython import display\n",
        "\n",
        "%matplotlib inline"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ABBhnkNEL-hX"
      },
      "source": [
        "## Data\n",
        "\n",
        "Do you remember Anscombe's Quartet? There are four sets of points here. All have the same means, variances, and best parameters (for a least square line)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TUblOQoQczp5"
      },
      "source": [
        "def anscombe(n=1):\n",
        "  \"\"\"Load one set of Anscombe's Quartet\"\"\"\n",
        "  assert 1 <= n and n <= 4\n",
        "  X = np.vstack([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5])\n",
        "  if n == 1:\n",
        "    y = [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]\n",
        "  elif n == 2:\n",
        "    y = [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]\n",
        "  elif n == 3:\n",
        "    y = [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]\n",
        "  elif n == 4:\n",
        "    X = [8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8]\n",
        "    y = [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]\n",
        "  X = np.vstack(X)\n",
        "  y = np.asarray(y)\n",
        "  return X, y\n",
        "\n",
        "X, y = anscombe()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LxPcKLvaNOh7"
      },
      "source": [
        "## The first model\n",
        "\n",
        "For simplicity, we can define a model with intercept and slope as explicit arguments."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xDjzsIgJaTdl"
      },
      "source": [
        "def simple_model(x, intercept, slope):\n",
        "  \"\"\"Simple linear model\"\"\"\n",
        "  ret = intercept + x*slope\n",
        "  return ret.ravel() # Sklearn output data is flat, like this"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "As_9k9zGg5hm"
      },
      "source": [
        "Each set in Anscombe's Quartet looks different but has the same parameters. You can change these to get a different linear model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tVuOM4DfhjGw"
      },
      "source": [
        "intercept = 3\n",
        "slope = .5"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vFtoBhbWtCtF"
      },
      "source": [
        "Now for plotting the loaded data and the model given predefined parameter values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tVhL_7mPdN1Y"
      },
      "source": [
        "plt.figure(figsize=(5, 5), dpi=100)\n",
        "\n",
        "x_plot = np.asarray([X.min(), X.max()])\n",
        "plt.plot(x_plot, simple_model(x_plot, intercept, slope), color='C1', linewidth=3, zorder=0)\n",
        "\n",
        "for i in range(X.shape[0]):\n",
        "  plt.plot([X[i], X[i]], [y[i], simple_model(X[i], intercept, slope)], 'k--', alpha=.3, zorder=-1)\n",
        "\n",
        "plt.scatter(X, y, color='C0', edgecolors='k', zorder=1)\n",
        "\n",
        "plt.axis('equal')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GPvwb9UAhoC-"
      },
      "source": [
        "By changing the parameters in a widget, we can see different models and their loss right away without rerunning cells."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qA5hf0pg_wZc"
      },
      "source": [
        "def interactive_residual_plot(intercept, slope):\n",
        "  # Create the plot\n",
        "  plt.figure(figsize=(5, 5), dpi=100)\n",
        "  # Plot regression line\n",
        "  x_plot = np.asarray([X.min(), X.max()])\n",
        "  plt.plot(x_plot, intercept + slope*x_plot, color='C1', linewidth=3, zorder=0, label=\"Trendline\")\n",
        "  # Plot training data\n",
        "  plt.scatter(X, y, color='C0', edgecolors='k', zorder=1, label=\"Data\")\n",
        "  # Error 'path'\n",
        "  for i in range(X.shape[0]):\n",
        "    plt.plot([X[i], X[i]], [y[i], intercept + slope*X[i]], 'k--', alpha=.3, zorder=-1)\n",
        "  # Fix look\n",
        "  plt.axis('equal')\n",
        "  plt.legend(loc='lower right')\n",
        "  # Find and print the residual\n",
        "  residual = np.sum((y - (intercept + slope*X.ravel()))**2)\n",
        "  plt.title(\"Residual (squared error): %.1f\" % residual)\n",
        "  # Show the plot\n",
        "  plt.show()\n",
        "\n",
        "interact_plot = ipywidgets.interact(interactive_residual_plot,\n",
        "                                    intercept=ipywidgets.FloatSlider(value=5, min=0, max=10, step=0.1, description=\"Intercept\"),\n",
        "                                    slope=ipywidgets.FloatSlider(value=0, min=-1, max=1, step=0.05, description=\"Slope\"));\n",
        "output = interact_plot.widget.children[-1] # This should prevent flickering\n",
        "output.layout.height = '500px'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OUKp_MjXXSIs"
      },
      "source": [
        "## Parameter optimization\n",
        " \n",
        "What are the best parameter values? Any best/good set of parameters are always evaluated using some quality metric, often called \"loss.\" For simple linear regression, the loss function is squared loss, as in the code above. This is defined as $\\sum_i (y_i - \\hat y_i)^2$. We could image other loss functions for simple regression, like using the $l_1$ distance instead of $l_2$ (squared euclidean distance). The loss function then becomes $\\sum_i |y_i - \\hat y_i|$.\n",
        "\n",
        "Maybe the most simple way of finding the best parameters would be to randomly choose the values for the parameters and only keep the set with the lowest loss (i.e., the lowest residual). Let's see what this might look like."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BPuEjg8Wmuv6"
      },
      "source": [
        "def model(x, theta):\n",
        "  \"\"\"Basic linear model\"\"\"\n",
        "  ret = theta[0] + x*theta[1]\n",
        "  return ret.ravel()\n",
        "  \n",
        "def lnorm_loss(y, y_hat, l=2):\n",
        "  \"\"\"Finds the residual given the target data (y) and the predictions (y_hat).\"\"\"\n",
        "  if l==2:\n",
        "    return np.sum((y - y_hat)**2) # l2 loss\n",
        "  elif l==1:\n",
        "    return np.sum(np.abs(y - y_hat)) # l1 loss"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kbnCxM6mmu56"
      },
      "source": [
        "for i in range(10):\n",
        "  theta = np.random.uniform(-10, 10, size=2)  # Random guess (within constraints) for parameters.\n",
        "  print(\"theta: %.1f, %.1f \\tl2 loss: %.1f   \\tl1 loss: %.1f\" % (theta[0], theta[1], lnorm_loss(y, model(X, theta), l=2), lnorm_loss(y, model(X, theta), l=1)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MRP1Pv6JnSBQ"
      },
      "source": [
        "This will take some time, even when the distribution is constrained. There must be a better way!\n",
        " \n",
        "We could imagine other ways. Below, I have defined a couple of ways to do optimization. Go through them all (not necessarily in order) and try to understand what is happening. Each algorithm uses the model, data, and residual defined above. When you want to change something from above, simply load the new data or rewrite (and reload) the loss function to use some other metric. Each algorithm will use your `model` and `loss` functions and return a history of parameters that can be used with the plotting functions further down in this notebook."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yrwhNV67ouVj"
      },
      "source": [
        "def coordinate_descent(function_to_minimize, initial_theta, n_iterations=50, step=.25):\n",
        "  \"\"\"Coordinate descent (look at wikipedia for definition)\"\"\"\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    theta_changed_during_iter = False\n",
        "    for dim in range(len(initial_theta)):\n",
        "      for direction in [-1, 1]:\n",
        "        new_theta = theta.copy()\n",
        "        new_theta[dim] += direction*step\n",
        "        while function_to_minimize(new_theta) < function_to_minimize(theta):\n",
        "          theta = new_theta\n",
        "          theta_changed_during_iter = True\n",
        "          history.append(theta)\n",
        "          function_values.append(function_to_minimize(new_theta))\n",
        "          new_theta = theta.copy()\n",
        "          new_theta[dim] += direction*step\n",
        "    if not theta_changed_during_iter:\n",
        "      step /= 2\n",
        "      # break\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def finite_difference_gradient_descent(function_to_minimize, initial_theta, n_iterations=20, gamma=.01):\n",
        "  \"\"\"Gradient descent optimization, where the gradient is approximated using\n",
        "  finite difference. Look at scipy.optimize.approx_fprime to learn more.\"\"\"\n",
        "  from scipy.optimize import approx_fprime\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    gradient = approx_fprime(theta, function_to_minimize, epsilon=.001)\n",
        "    theta = theta - gamma*gradient\n",
        "    history.append(theta)\n",
        "    function_values.append(function_to_minimize(theta))\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def random_changes(function_to_minimize, initial_theta, n_iterations=200, noise_sd=.1):\n",
        "  \"\"\"Adds random noise to the parameter vector in the hope of \n",
        "  improving the loss.\"\"\"\n",
        "  theta = initial_theta.copy()\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  for iteration in range(n_iterations):\n",
        "    new_theta = theta + np.random.normal(scale=noise_sd, size=theta.shape)\n",
        "    if function_to_minimize(new_theta) < function_to_minimize(theta):\n",
        "      theta = new_theta\n",
        "    history.append(theta.copy())\n",
        "    function_values.append(function_to_minimize(theta))\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])\n",
        "\n",
        "def scipy_minimize(function_to_minimize, initial_theta):\n",
        "  \"\"\"Optimization using external library. Look at the manual of \n",
        "  scipy.optimize.minimize for specifics.\"\"\"\n",
        "  from scipy.optimize import minimize\n",
        "  history = [initial_theta.copy()]\n",
        "  function_values = [function_to_minimize(initial_theta)]\n",
        "  def store_theta_during_optimization(current_theta):\n",
        "    history.append(current_theta)\n",
        "    function_values.append(function_to_minimize(current_theta))\n",
        "  res = minimize(lambda t: function_to_minimize(t), initial_theta, \n",
        "                 callback=store_theta_during_optimization)\n",
        "  return theta, function_values, np.concatenate([np.vstack(h).T for h in history])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "njcuFRAzo8nR"
      },
      "source": [
        "def plot_loss_landscape(history, loss_function):\n",
        "  xx, yy = np.meshgrid(np.linspace(min(history[:, 0].min(), -5), max(history[:, 0].max(), 10), 100), # Intercept\n",
        "                      np.linspace(min(history[:, 1].min(), -1), max(history[:, 1].max(), 2), 100)) # Slope\n",
        "  Z = np.apply_along_axis(lambda theta: loss_function(theta), 1, \n",
        "                          np.concatenate([np.vstack(xx.ravel()), np.vstack(yy.ravel())], axis=1))\n",
        "  Z = Z.reshape(xx.shape)\n",
        "\n",
        "  fig = plt.figure(figsize=(6, 6))\n",
        "  ax = fig.subplots(1, 1)\n",
        "  ax.contourf(xx, yy, np.log(Z))\n",
        "  ax.set_xlabel(\"Intercept\")\n",
        "  ax.set_ylabel(\"Slope\")\n",
        "  ax.plot([h[0] for h in history], [h[1] for h in history], 'r.-', alpha=.7, zorder=1, label=\"Path to optimum\")\n",
        "  ax.scatter(history[-1, 0], history[-1, 1], s=100, color='r', edgecolors='k', zorder=2, label=\"Optimum\")\n",
        "  ax.legend()\n",
        "  fig.tight_layout()\n",
        "  fig.show()\n",
        "\n",
        "def plot_loss(loss):\n",
        "  fig = plt.figure(figsize=(8, 4))\n",
        "  ax = fig.subplots(1, 1)\n",
        "  ax.plot(loss, linewidth=2)\n",
        "  # ax.set_yscale('log')\n",
        "  ax.set_ylabel(\"Loss\")\n",
        "  ax.set_xlabel(\"Step\")\n",
        "  fig.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SlhJbP1kCN34"
      },
      "source": [
        "# Be carefull when defining functions like this. It convinient here but might\n",
        "# create confusion due to it using global data.\n",
        "loss_function = lambda theta: lnorm_loss(y, model(X, theta), l=2)\n",
        "\n",
        "# theta, loss, history = coordinate_descent(function_to_minimize=loss_function,\n",
        "#                                           initial_theta=np.random.uniform(-5, 5, size=2), \n",
        "#                                           n_iterations=5000, step=.1)\n",
        "\n",
        "# theta, loss, history = finite_difference_gradient_descent(function_to_minimize=loss_function,\n",
        "#                                           initial_theta=np.random.uniform(-5, 5, size=2), \n",
        "#                                           n_iterations=20, gamma=.0001)\n",
        "\n",
        "theta, loss, history = random_changes(function_to_minimize=loss_function,\n",
        "                                      initial_theta=np.random.uniform(-5, 5, size=2),\n",
        "                                      n_iterations=200, noise_sd=.1)\n",
        "\n",
        "# theta, loss, history = scipy_minimize(function_to_minimize=loss_function,\n",
        "#                                       initial_theta=np.random.uniform(-5, 5, size=2))\n",
        "\n",
        "plot_loss(loss)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Vv4kcgzj9ZPj"
      },
      "source": [
        "plot_loss_landscape(history, loss_function=loss_function)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1QKSRDEGezMY"
      },
      "source": [
        "def plot_wrapper(iteration):\n",
        "  intercept, slope = history[iteration]\n",
        "  print(\"Intercept: %.1f, Slope: %.1f\" % (intercept, slope))\n",
        "  interactive_residual_plot(intercept, slope)\n",
        "interact_plot = ipywidgets.interact(plot_wrapper,\n",
        "                                    iteration = ipywidgets.Play(min=0, max=len(history)-1, step=max(1, len(history)//100), value=0));\n",
        "output = interact_plot.widget.children[-1] # This should prevent flickering\n",
        "output.layout.height = '500px'\n",
        "# Note that the play button doesn't always show up in colab. The clickable area is still there though (leftmost gray area)."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "52MK_jKi79cv"
      },
      "source": [
        "from mpl_toolkits.mplot3d.axes3d import Axes3D\n",
        "\n",
        "xx, yy = np.meshgrid(np.linspace(-5, 10, 100), # Intercept\n",
        "                    np.linspace(-1, 2, 100)) # Slope\n",
        "\n",
        "fig = plt.figure(figsize=(16, 8))\n",
        "for order in [1, 2]:\n",
        "  ax = fig.add_subplot(1, 2, order, projection='3d')\n",
        "  Z = np.apply_along_axis(lambda theta: lnorm_loss(y, model(X, theta), l=order), 1, \n",
        "                          np.concatenate([np.vstack(xx.ravel()), np.vstack(yy.ravel())], axis=1))\n",
        "  Z = Z.reshape(xx.shape)\n",
        "  ax.plot_surface(xx, yy, Z, cmap='viridis', rstride=5, cstride=5)\n",
        "  ax.set_title(\"$l_%i$ loss surface\" % order)\n",
        "  ax.set_xlabel(\"Intercept\")\n",
        "  ax.set_ylabel(\"Slope\")\n",
        "  ax.set_zlabel(\"Loss\")\n",
        "  ax.view_init(30, 120)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KNjpWKQY7PpW"
      },
      "source": [
        "## Probabilistic model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fnfkT7A4-cXc"
      },
      "source": [
        "This is an implementation of the slide on Bayesian linear regression. To get a simple model and also being true to the lecture slides, the hack for the variance is also implemented. Try to understand the model, using both the code and the slides. Minimizing the negative of this loss function is the same as maximizing the probability (Maximum A Posteriori, MAP). See this as both an exercise in modelling and parsing code."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_F2gb4k__WYi"
      },
      "source": [
        "def bayesian_model(x, theta):\n",
        "  \"\"\"Bayesian linear model\"\"\"\n",
        "  ret = theta[0] + x*theta[1]\n",
        "  return ret.ravel()\n",
        "\n",
        "from scipy.stats import norm\n",
        "\n",
        "def model_log_likelihood(X, y, theta):\n",
        "  \"\"\"log likelihood of the model\n",
        "  Note that the bayesian_model function name is hardcoded here. This is bad design.\"\"\"\n",
        "  return np.sum(norm(scale=theta[2]).logpdf(y-bayesian_model(X, theta)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GU4aB4tuEO45"
      },
      "source": [
        "print(\"Intercept\\tSlope\\tNoise\\t\\tLog loss\")\n",
        "for i in range(10):\n",
        "  theta = np.random.uniform(-10, 10, size=3)\n",
        "  theta[2] = np.abs(theta[2]) # Theta_2 must not be negative, why?\n",
        "  print(\"%.1f\\t\\t%.1f\\t%.1f\\t->\\t%.1f\" % (theta[0], theta[1], theta[2], model_log_likelihood(X, y, theta)))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GczkrLnFBUbc"
      },
      "source": [
        "from mpl_toolkits.mplot3d.axes3d import Axes3D\n",
        "\n",
        "xx, yy = np.meshgrid(np.linspace(-5, 10, 100), # Intercept\n",
        "                    np.linspace(-1, 2, 100)) # Slope\n",
        "\n",
        "fig = plt.figure(figsize=(16, 8))\n",
        "for i, noise_level in enumerate([2, 6]):\n",
        "  ax = fig.add_subplot(1, 2, i+1, projection='3d')\n",
        "  thetas_on_grid = np.concatenate([np.vstack(xx.ravel()), np.vstack(yy.ravel())], axis=1)\n",
        "  thetas_on_grid = np.append(thetas_on_grid, noise_level*np.ones((thetas_on_grid.shape[0], 1)), axis=1)\n",
        "  Z = np.apply_along_axis(lambda theta: model_log_likelihood(X, y, theta), 1, thetas_on_grid)\n",
        "  Z = Z.reshape(xx.shape)\n",
        "  ax.plot_surface(xx, yy, np.exp(Z), cmap='viridis')\n",
        "  ax.set_title(\"Distribution over intercept and slope with fixed noise=%.1f\" % noise_level)\n",
        "  ax.set_xlabel(\"Intercept\")\n",
        "  ax.set_ylabel(\"Slope\")\n",
        "  ax.set_zlabel(\"Log Likelihood\")\n",
        "  ax.view_init(30, 120)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7nkKXg5bCa59"
      },
      "source": [
        "from scipy.optimize import minimize\n",
        "\n",
        "history = [np.asarray([3, .5, .2])]\n",
        "\n",
        "def store_theta_during_optimization(current_theta):\n",
        "  history.append(current_theta)\n",
        "res = minimize(lambda theta: -model_log_likelihood(X, y, theta), history[0], \n",
        "               bounds=[(None, None), (None, None), (1e-3, None)],\n",
        "               callback=store_theta_during_optimization)\n",
        "theta = res.x\n",
        "\n",
        "print(\"MAP estimate of theta: %s\" % theta)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IcN_p6k3FRME"
      },
      "source": [
        "intercept, slope = theta[:2]\n",
        "plt.figure(figsize=(5, 5), dpi=100)\n",
        "# Plot regression line\n",
        "x_plot = np.asarray([X.min(), X.max()])\n",
        "plt.plot(x_plot, intercept + slope*x_plot, color='C1', linewidth=3, zorder=0, label=\"Trendline\")\n",
        "plt.plot(x_plot, theta[2] + intercept + slope*x_plot, '--', color='C1', linewidth=1, zorder=0, label=\"1st sd\")\n",
        "plt.plot(x_plot, -theta[2] + intercept + slope*x_plot, '--', color='C1', linewidth=1, zorder=0)\n",
        "plt.plot(x_plot, 2*theta[2] + intercept + slope*x_plot, '-.', color='C1', linewidth=1, zorder=0, label=\"2nd sd\")\n",
        "plt.plot(x_plot, -2*theta[2] + intercept + slope*x_plot, '-.', color='C1', linewidth=1, zorder=0)\n",
        "# Plot training data\n",
        "plt.scatter(X, y, color='C0', edgecolors='k', zorder=1, label=\"Data\")\n",
        "# Error 'path'\n",
        "for i in range(X.shape[0]):\n",
        "  plt.plot([X[i], X[i]], [y[i], intercept + slope*X[i]], 'k--', alpha=.3, zorder=-1)\n",
        "# Fix look\n",
        "plt.axis('equal')\n",
        "plt.legend(loc='lower right')\n",
        "# Find and print the residual\n",
        "residual = np.sum((y - (intercept + slope*X.ravel()))**2)\n",
        "plt.title(\"Residual (squared error): %.1f\" % residual)\n",
        "# Show the plot\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}