lectures/Lecture_21.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Deep Learning: I\n",
    "\n",
    "<font color='red'>CHANGE BOSTON HOUSING TO CALIFORNIA HOUSING</font>\n",
    "\n",
    "*S. R. Taylor (2022)*\n",
    "\n",
    "This lecture and notebook are based on the \"Classification3\", \"NeuralNetworks\", and \"TensorFloeCodeOnly\" lectures of of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), which in turn are based on materials from Andy Connolly, and Ivezic et al. Chapter 9, Andy Connolly's [blog](http://connolly.github.io/introAstroML/blog/regression.html), and Aurelien Geron's [book](https://github.com/ageron/handson-ml2). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 9.\n",
    "- Many blogs and videos.\n",
    "- Free online book! http://neuralnetworksanddeeplearning.com/index.html\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [Loss functions](#one)\n",
    "* [Gradient descent](#two)\n",
    "* [AdaBoost](#three)\n",
    "* [Neural networks](#four)\n",
    "    \n",
    "---\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---\n",
    "\n",
    "We'll be talking about neural networks and deep learning today. [This](https://www.youtube.com/watch?v=bxe2T-V8XRs&list=PLiaHhY2iBX9hdHaRr6b7XevZtgZRa1PoU&ab_channel=WelchLabs) series of videos provides a particularly helpful intro that simplifies the explanation. Before we delve into neural networks, let me introduce a few useful concepts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Loss Functions <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "A **[loss function](https://en.wikipedia.org/wiki/Loss_functions_for_classification)** is like a cost function, or a likelihood, or optimization function, or an objective function. We're trying to minimize the offset between a model and some data. The difference with a loss function is that we evaluate it on a single training example rather than the full data set.\n",
    "\n",
    "- Whether you realize it or not, you are typically working with **$L2$ loss functions**:\n",
    "$$(y-f(x))^2,$$\n",
    "where the corresponding ***cost function*** is the mean of those for all $x_i$ or the **Mean Squared Error (MSE)** (just another name for $\\chi^2$, which of course is related to the probability of Gaussian uncertainties).\n",
    "\n",
    "\n",
    "- We also talked about **$L1$ loss functions** way back when we looked at the ***Huber loss*** which was more robust to outliers:\n",
    "$$|y-f(x)|,$$\n",
    "and the L1 penalty was important in LASSO regression.\n",
    "\n",
    "For classification we plot the **loss vs. $(y\\times f(x))$**, where the latter corresponds to the known class (either $+1$ or $-1$) times the predicted value. If that product is positive, we predict $+1$. If it is negative, we predict $-1$. We define the loss to be zero for $y*f(x)=1$, that is, when we have gotten the right answer.\n",
    "\n",
    "So what do the mean squared error (MSE or $L2$) and mean absolute error (MAE or $L1$) look like for classification?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Mathematical formulas for various loss functions\n",
    "def log_loss(raw_model_output):\n",
    "    return np.log(1+np.exp(-raw_model_output))\n",
    "\n",
    "def hinge_loss(raw_model_output):\n",
    "    return np.maximum(0,1-raw_model_output)\n",
    " \n",
    "def l2(raw_model_output):\n",
    "    return (raw_model_output-1)**2  \n",
    "\n",
    "def l1(raw_model_output):\n",
    "    return np.abs(raw_model_output-1)  \n",
    " \n",
    "def zero_one(raw_model_output):\n",
    "    return np.where(raw_model_output < 0, 1, 0)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Create a grid of values and plot\n",
    "grid = np.linspace(-3,3,1000)\n",
    "plt.plot(grid, l2(grid), \"g\", label='L2')\n",
    "plt.plot(grid, l1(grid), \"brown\", label='L1')\n",
    "\n",
    "plt.fill_between([0,3], y1=-0.02, y2=5,\n",
    "                 color=\"b\", alpha=0.2)\n",
    "plt.fill_between([-3,0], y1=-0.02, y2=5,\n",
    "                 color=\"r\", alpha=0.2)\n",
    "plt.xlim([-3,3])\n",
    "plt.ylim([-0.02,5])\n",
    "plt.xlabel(\"y*f(x)\", fontsize=12)\n",
    "plt.ylabel(\"loss\", fontsize=12)\n",
    "plt.title(\"predict -1 (incorrect)         predict +1 (correct)\", fontsize=12)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This does something reasonable for $y*f(x)\\le1$.  However, look what happens at larger values (where we are even more confident that $y*f(x)$ is positive and that our class should be $+1$.  The loss goes **up**.  That's bad.\n",
    "\n",
    "Now, you may be wondering how $y*f(x)$ can be larger than 1. Here's how this works.\n",
    "\n",
    "$f(x)$ isn't just a value between $-1$ and $1$-- it is a function. For example, let's say that our training data looks like this (where `xx` are features and `yy` are class labels):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "xx = np.array([0.1, 0.2, 0.4, 0.6,\n",
    "               0.8, 1.1, 1.4, 1.5])\n",
    "yy = np.array([-1, -1, -1, -1,\n",
    "               1, 1, 1, 1])\n",
    "\n",
    "if 1:\n",
    "    \n",
    "    fig = plt.figure(1)\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "    # Move left y-axis and bottim x-axis to centre, passing through (0,0)\n",
    "    ax.spines['left'].set_position(('axes',0.045))\n",
    "    ax.spines['bottom'].set_position('center')\n",
    "\n",
    "    # Eliminate upper and right axes\n",
    "    ax.spines['right'].set_color('none')\n",
    "    ax.spines['top'].set_color('none')\n",
    "\n",
    "    # Show ticks in the left and lower axes only\n",
    "    ax.xaxis.set_ticks_position('bottom')\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "    \n",
    "    ax.plot(1, 0, \">k\", transform=ax.get_yaxis_transform(), clip_on=False)\n",
    "    ax.plot(0, 1, \"^k\", transform=ax.get_xaxis_transform(), clip_on=False)\n",
    "        \n",
    "    ax.scatter(xx,yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now let's fit a linear model to the training data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(xx[:,None],yy)\n",
    "ypred = linreg.predict(grid[:,None])\n",
    "\n",
    "if 1:\n",
    "    \n",
    "    fig = plt.figure(1)\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "    # Move left y-axis and bottim x-axis to centre, passing through (0,0)\n",
    "    ax.spines['left'].set_position(('axes',0.25))\n",
    "    ax.spines['bottom'].set_position('center')\n",
    "\n",
    "    # Eliminate upper and right axes\n",
    "    ax.spines['right'].set_color('none')\n",
    "    ax.spines['top'].set_color('none')\n",
    "\n",
    "    # Show ticks in the left and lower axes only\n",
    "    ax.xaxis.set_ticks_position('bottom')\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "    \n",
    "    ax.plot(1, 0, \">k\", transform=ax.get_yaxis_transform(), clip_on=False)\n",
    "    ax.plot(0, 1, \"^k\", transform=ax.get_xaxis_transform(), clip_on=False)\n",
    "    \n",
    "    ax.plot(grid,ypred)\n",
    "    ax.set_xlim([-1,3])\n",
    "    ax.set_ylim([-2,2])\n",
    "        \n",
    "    ax.scatter(xx,yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We see that for $x$ greater than about 1.3, $f(x)$ can indeed be larger than 1 and so can $y*f(x)$, which is indicating an increased certainty of the $+1$ class. \n",
    "\n",
    "OK, so now we can understand the plot, but we still need a loss function that makes sense for classification.\n",
    "\n",
    "- The first we'll try is the so-called **[\"Zero-One Loss\"](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.zero_one_loss.html)** shown in **black** below. It is 1 for $yf(x)<0$ and 0 for $yf(x)>0$; thus the name. You increment the loss function by 1 every time you make a wrong prediction. It is just a count of the total number of mistakes.\n",
    "\n",
    "However, the Zero-One loss is hard to minimize, so instead we can try something that allows the loss to be continuous function in $y*f(x)$.  \n",
    "\n",
    "- The **[Hinge Loss](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.hinge_loss.html)**, which looks like\n",
    "$${\\rm max}(0,1-y*f(x)),$$\n",
    "is plotted in **orange**. Here there is no contribution to the loss for values $\\ge 1$, but there is a linearly increasing loss for smaller values. So, it penalizes both wrong predictions and also correct predictions that have low confidence.\n",
    "\n",
    "- A **[Logistic Loss](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.log_loss.html#sklearn.metrics.log_loss)** (also called the *log loss* and *cross entropy loss*) function has similar properties as shown in **blue**, but is smoother and has slightly less and less penalty for more and more confident $+1$ predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Create a grid of values and plot\n",
    "grid = np.linspace(-3,3,1000)\n",
    "\n",
    "plt.plot(grid, zero_one(grid), \"k\", label='0-1')\n",
    "plt.plot(grid, hinge_loss(grid), \"orange\", label='hinge')\n",
    "plt.plot(grid, log_loss(grid), \"b\", label='logistic')\n",
    "\n",
    "plt.fill_between([0,3], y1=-0.02, y2=5,\n",
    "                 color=\"b\", alpha=0.2)\n",
    "plt.fill_between([-3,0], y1=-0.02, y2=5,\n",
    "                 color=\"r\", alpha=0.2)\n",
    "plt.xlim([-3,3])\n",
    "plt.ylim([-0.02,5])\n",
    "plt.xlabel(\"y*f(x)\",fontsize=12)\n",
    "plt.ylabel(\"loss\",fontsize=12)\n",
    "plt.title(\"predict -1 (incorrect)         predict +1 (correct)\",fontsize=12)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For more see\n",
    "- [Linear Classsifiers in Python course](https://learn.datacamp.com/courses/linear-classifiers-in-python).\n",
    "- https://datascience103579984.wordpress.com/2019/09/18/linear-classifiers-in-python-from-datacamp/\n",
    "- https://towardsdatascience.com/common-loss-functions-in-machine-learning-46af0ffc4d23\n",
    "- https://www.analyticsvidhya.com/blog/2019/08/detailed-guide-7-loss-functions-machine-learning-python-code\n",
    "- http://www.datasciencecourse.org/notes/linear_classification/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gradient Descent <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "That brings us to the topic of **[Gradient Descent](https://en.wikipedia.org/wiki/Gradient_descent)** for finding the optimal extremal position of a cost function.\n",
    "\n",
    "Throughout this couse we have been trying to determine model parameters, $\\theta$, that minimize either the regression error or the classification error when fitting our training data (and not overfitting!). \n",
    "\n",
    "- Sometimes we have been able to write an analytic solution for $\\theta$.\n",
    "- In MCMC we semi-randomly sampled the multi-dimensional $\\theta$ space to find the best answer (and map the full parameter posterior distribution along the way). \n",
    "- But what happens if you are stuck on top of a freezing cold mountain and you have no map and can't magically jump from place to place? You start walking **down**. That's the basic idea of ***gradient descent***--take a look around you, figure out which way is sloping downward the most and go *that* way.\n",
    "\n",
    "We are going to determine the local gradient of the loss function with respect to $\\theta$ and go in the steepest direction, until the gradient is zero (and we have arrived at our destination).\n",
    "\n",
    "Mathematically, for a simple linear model the MSE is defined as\n",
    "\n",
    "$$ {\\rm MSE}({\\mathbf \\theta}) = \\frac{1}{N}(Y - X\\theta)^T (Y - X\\theta).$$\n",
    "\n",
    "Therefore the partial derivative with respect to the vector of parameters $\\theta$ is\n",
    "\n",
    "$$\\nabla_{\\theta}{\\rm MSE}({\\mathbf \\theta}) = \\frac{2}{N}X^T(X\\theta - Y).$$\n",
    "\n",
    "That gives the uphill direction, but we want to move against the gradient, so we compute the next step as\n",
    "\n",
    "$$\\theta^{\\rm next step} = \\theta - \\eta\\nabla_{\\theta}{\\rm MSE}({\\mathbf \\theta}),$$\n",
    "\n",
    "where $\\eta$ is the **\"learning rate\"** and the rest are all matrices or vectors. \n",
    "\n",
    "- The initial values for $\\theta$ are chosen randomly.\n",
    "- The **[learning rate](https://en.wikipedia.org/wiki/Learning_rate)** controls how big your steps \"down\" are. \n",
    "- If your step size is too small, it will take too long to converge. \n",
    "- If it is too big, you might miss the bottom completely (and possibly end up diverging from the solution).  \n",
    "\n",
    "\n",
    "![https://miro.medium.com/max/1400/0*GaO7X6j3coh3oNwf.png](https://miro.medium.com/max/1400/0*GaO7X6j3coh3oNwf.png)\n",
    "\n",
    "We also have to be careful that we don't end up in a local minimum instead of the global minimum. One of the nice things about the $L2$ cost function is that it is guaranteed to have just a single global minimum. Gradient descent is also useful for regression where there are too many training points (or too many features) to fit into memory at a given time.\n",
    "\n",
    "<font color='red'>Here's an example from Geron where we apply gradient descent to a simple linear regression problem.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Run the next 4 cells\n",
    "# N points randomly drawn from a linear distribution\n",
    "N = 100\n",
    "X = 2 * np.random.rand(N, 1)\n",
    "y = 4 + 3 * X + np.random.randn(N, 1)\n",
    "\n",
    "# Turn X into a matrix\n",
    "X_b = np.c_[np.ones((100, 1)), X]  # add x0 = 1 to each \n",
    "\n",
    "# Grid for plotting\n",
    "X_new = np.array([[0], [2]])\n",
    "X_new_b = np.c_[np.ones((2, 1)), X_new]  # add x0 = 1 to each instance\n",
    "\n",
    "X_b = np.c_[np.ones((100, 1)), X]  # add x0 = 1 to each instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eta = 0.1  # learning rate\n",
    "n_iterations = 100\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "\n",
    "for iteration in range(n_iterations):\n",
    "    gradients = 2/N * X_b.T.dot(X_b.dot(theta) - y)\n",
    "    theta = theta - eta * gradients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n_lines=10\n",
    "color_idx = np.linspace(0, 1, n_lines)\n",
    "\n",
    "# Helper function\n",
    "theta_path_bgd = []\n",
    "def plot_gradient_descent(theta, eta, theta_path=None):\n",
    "    m = len(X_b)\n",
    "    plt.plot(X, y, \"b.\")\n",
    "    n_iterations = 1000\n",
    "    \n",
    "    for iteration in range(n_iterations):\n",
    "        if iteration < 10:\n",
    "            y_predict = X_new_b.dot(theta)\n",
    "            plt.plot(X_new, y_predict, color=plt.cm.cool(color_idx[iteration]))\n",
    "            \n",
    "        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)\n",
    "        theta = theta - eta * gradients\n",
    "        if theta_path is not None:\n",
    "            theta_path.append(theta)\n",
    "            \n",
    "    plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "    plt.axis([0, 2, 0, 15])\n",
    "    plt.title(r\"$\\eta = {}$\".format(eta), fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "plt.figure(figsize=(12,4))\n",
    "\n",
    "plt.subplot(131); plot_gradient_descent(theta, eta=0.02)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "\n",
    "plt.subplot(132); plot_gradient_descent(theta, eta=0.1, \n",
    "                                        theta_path=theta_path_bgd)\n",
    "plt.subplot(133); plot_gradient_descent(theta, eta=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We've got a \"three bears\" scenario.\n",
    "\n",
    "- On the left our learning rate is too low. We'll eventually get to the solution, but it will take a long time.  \n",
    "- On the right it is too high and we have completely missed the solution.  \n",
    "- In the middle the learning rate is just right.\n",
    "\n",
    "<font color='red'>Try $\\eta = 0.3$ and $0.4$ to see if you can understand what is going on in the right panel.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## AdaBoost <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "When I introduced Boosting in *Classification II*, I didn't mention **[AdaBoost](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.AdaBoostClassifier.html)** or **[Gradient Boosting](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html?highlight=gradient%20boosting#sklearn.ensemble.GradientBoostingClassifier)**, because it makes more sense to know about the other items above first. But ***AdaBoost = \"Adaptive Boosting\"***.\n",
    "\n",
    "Below is an example of AdaBoost from Geron. <font color='red'>See what happens when you change the learning rate, where half the learning rate means that weights are boosted half as much for each iteration.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Run the next 3 cells\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.datasets import make_moons\n",
    "\n",
    "X, y = make_moons(n_samples=500, noise=0.30, random_state=42)\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Helper function for plotting\n",
    "from matplotlib.colors import ListedColormap\n",
    "def plot_decision_boundary(clf, X, y, \n",
    "                           axes=[-1.5, 2.45, -1, 1.5], \n",
    "                           alpha=0.5, contour=True):\n",
    "    x1s = np.linspace(axes[0], axes[1], 100)\n",
    "    x2s = np.linspace(axes[2], axes[3], 100)\n",
    "    x1, x2 = np.meshgrid(x1s, x2s)\n",
    "    \n",
    "    X_new = np.c_[x1.ravel(), x2.ravel()]\n",
    "    y_pred = clf.predict(X_new).reshape(x1.shape)\n",
    "    \n",
    "    custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])\n",
    "    plt.contourf(x1, x2, y_pred, alpha=0.3, cmap=custom_cmap)\n",
    "    if contour:\n",
    "        custom_cmap2 = ListedColormap(['#7d7d58','#4c4c7f','#507d50'])\n",
    "        plt.contour(x1, x2, y_pred, cmap=custom_cmap2, alpha=0.8)\n",
    "    plt.plot(X[:, 0][y==0], X[:, 1][y==0], \"yo\", alpha=alpha)\n",
    "    plt.plot(X[:, 0][y==1], X[:, 1][y==1], \"bs\", alpha=alpha)\n",
    "    \n",
    "    plt.axis(axes)\n",
    "    plt.xlabel(r\"$x_1$\", fontsize=18)\n",
    "    plt.ylabel(r\"$x_2$\", fontsize=18, rotation=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "learning_rate=1.0\n",
    "m = len(X_train)\n",
    "\n",
    "fix, axes = plt.subplots(ncols=1, figsize=(10,6))\n",
    "sample_weights = np.ones(m)\n",
    "for i in range(5):\n",
    "    svm_clf = SVC(kernel=\"rbf\", C=0.05, gamma=\"scale\", random_state=42)\n",
    "    svm_clf.fit(X_train, y_train, sample_weight=sample_weights)\n",
    "    \n",
    "    y_pred = svm_clf.predict(X_train)\n",
    "    sample_weights[y_pred != y_train] *= (1 + learning_rate)\n",
    "    \n",
    "    plot_decision_boundary(svm_clf, X, y, alpha=0.2)\n",
    "    plt.title(\"learning_rate = {}\".format(learning_rate), fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<font color='red'>Try some learning rates between $0.1$ and $1.0$.  It is important to note the order of these decision boundaries; that is, are they diverging?</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Neural Networks <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "**[Artificial Neural Networks](https://en.wikipedia.org/wiki/Artificial_neural_network)** are a simplified computational architecture based loosely on the real neural networks found in brains.\n",
    "\n",
    "![Neuron example](https://4.bp.blogspot.com/-Z5LfY6yoIcE/U-OFKWHoAbI/AAAAAAAAAKo/ytH6BzDLeo4/s1600/Picture-533.png)\n",
    "\n",
    "In reality, what we are going to explore is a **[multi-layer perceptron](https://en.wikipedia.org/wiki/Multilayer_perceptron)**.\n",
    "\n",
    "In the image below, \n",
    "- the circles on the ***left*** represent the **features/attributes** of our input data, $X$, which here is 3 dimensional.  \n",
    "- the circles in the ***middle*** represent the **neurons**. They take in the information from the input and, based on some criterion decide whether or not to \"fire\". These middle layers are called \"**hidden layers**\".\n",
    "- the collective results of the neurons in the hidden layer produce the **output**, $y$, which is represented by the circles on the ***right***, which here is 2 dimensional result.  \n",
    "- the lines connecting the circles represent the synapses.  \n",
    "\n",
    "This is a simple example with just one layer of neurons; however, there can be many layers of neurons.\n",
    "\n",
    "![Cartoon of Neural Network](https://upload.wikimedia.org/wikipedia/commons/thumb/e/e4/Artificial_neural_network.svg/500px-Artificial_neural_network.svg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In more detail\n",
    "\n",
    "The job of a synapse is to take input values and multiply them by some **weight**, $w$, and add a **bias**, $b$, before passing them to the neuron (hidden layer):\n",
    "\n",
    "$$z = \\sum_i w x_i + b$$\n",
    "\n",
    "The bias determines the input level at which the neuron \"fires\". It is always present, and unique to each neuron, but we'll set that to zero for the sake of simplicity here.\n",
    "\n",
    "The neuron then sums up the inputs from all of the synapses connected to it and applies an \"**activation function**\", e.g., a **[sigmoid](https://en.wikipedia.org/wiki/Sigmoid_function) activation function**.\n",
    "\n",
    "$$a = \\frac{1}{1+e^{-z}}.$$\n",
    "\n",
    "![Sigmoid Function](https://upload.wikimedia.org/wikipedia/commons/thumb/8/88/Logistic-curve.svg/500px-Logistic-curve.svg.png)\n",
    "\n",
    "> ***The neural network learns the weights and biases of the synapses that are needed to produce an accurate model of $y_{\\rm train}$.***\n",
    "\n",
    "![Ivezic Figure 9.17](https://www.astroml.org/_images/fig_neural_network_1.png)\n",
    "\n",
    "- Rather than think about the inputs individually, we can write this process in matrix form as\n",
    "$$Z^{(2)} = X W^{(1)}$$\n",
    "\n",
    "\n",
    "- If $D$ is the number of attributes and $H$ is the number of neurons in the hidden layer, then $X$ is an $N\\times D$ matrix, while $W^{(1)}$ is a $D\\times H$ matrix.  The result, $Z^{(2)}$, is then an $N\\times H$ matrix.\n",
    "\n",
    "\n",
    "- We then apply the activation function to each entry of $Z^{(2)}$ independently: \n",
    "$$A^{(2)} = f(Z^{(2)}),$$\n",
    "where $A^{(2)}$ is the output of the neurons in the hidden layer and is also $N\\times H$.\n",
    "\n",
    "\n",
    "- These values are then the inputs for the next set of synapses, where we multiply the inputs by another set of weights, $W^{(2)}:$\n",
    "$$Z^{(3)} = A^{(2)} W^{(2)},$$\n",
    "where $W^{(2)}$ is an $H\\times O$ matrix and $Z^{(3)}$ is an $N\\times O$ matrix with $O$-dimensional output.\n",
    "\n",
    "\n",
    "- Another (potentially different type of) activation function is then applied to $Z^{(3)}$ to give\n",
    "$$\\hat{y} = g(Z^{(3)}),$$\n",
    "which is our estimator of $y$.\n",
    "\n",
    "---\n",
    "\n",
    "For example we might have $N=100$ people with known height and weight, for whom we have measured\n",
    "\n",
    "1. shoe size\n",
    "2. belt size\n",
    "3. hat size \n",
    "\n",
    "We are going to use this to predict the height and weight for people where we only know shoe size, belt size, and hat size.\n",
    "\n",
    "The neural network essentially boils down to determining the weights and biases of the synapses, which are usually initialized randomly. We do that by minimizing the cost function (which compares the true values of $y$ to our predicted values).  Typically an MSE cost:\n",
    "\n",
    "$$ {\\rm Cost} = J = \\sum\\frac{1}{2}(y - \\hat{y})^2.$$\n",
    "\n",
    "If we just had 1 weight and we wanted to check 1000 possible values, that wouldn't be so bad. ***But...*** if we have 20 weights, that means checking $20^{1000}$ possible combinations. Remember the curse of dimensionality?  That might take a while. Indeed, far, far longer than the age of the Universe.\n",
    "\n",
    "### Backpropagation\n",
    "\n",
    "However some clever techniques have evolved out of realizing that we can write an **analytic formula for the *gradient* going backwards through the network, and use that to update our weights (and biases, of course)**.\n",
    "\n",
    "For example, how about just checking 3 points for each weight and see if we can at least figure out which way is \"down hill\"? That's a start. We can rewrite $J$ as\n",
    "\n",
    "$$ J = \\sum\\frac{1}{2}\\left[y - g\\left( f(X W^{(1)}) W^{(2)} \\right) \\right]^2$$\n",
    "\n",
    "and then compute\n",
    "\n",
    "$$\\frac{\\partial J}{\\partial W}$$\n",
    "\n",
    "in order to determine the slope of the cost function for each weight.  This is the **gradient descent** method, which we encountered before.  Your choice of cost function is important here; specifically you want it to be differentiable.\n",
    "\n",
    "We'll want $\\partial J/\\partial W^{(1)}$ and $\\partial J/\\partial W^{(2)}$ separately. This allows us to ***[backpropagate](https://en.wikipedia.org/wiki/Backpropagation)*** the error contributions along each neuron and to change the weights where they most need to be changed.  It is like each observation gets a vote on which way is \"down hill\".  We compute the vector sum to decide the ultimate down hill direction.\n",
    "\n",
    "Once we know the down hill direction from the derivative, we update the weights by subtracting a scalar (the *learning rate*) times that derivative from the original weights-- this is just gradient descent! This is obviously much faster than randomly sampling all the possible combinations of weights.  \n",
    "\n",
    "<font color='red'>Watch the following short video to understand more of the mathematics of **backpropagation (= \"backward propagation of errors\")**.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('https://www.youtube.com/watch?v=GlcnxUlrtek&list=PLiaHhY2iBX9hdHaRr6b7XevZtgZRa1PoU&index=4&ab_channel=WelchLabs', width=1000, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Once the weights and biases are set through training, you've got your Neural Network classifier/regressor.***\n",
    "\n",
    "![Ivezic Figure 9.17](https://www.astroml.org/_images/fig_neural_network_1.png)\n",
    "\n",
    "Scikit-Learn has both **[unsupervised Neural Network](http://scikit-learn.org/stable/modules/neural_networks_unsupervised.html#neural-networks-unsupervised)** and **[supervised Neural Network](http://scikit-learn.org/stable/modules/neural_networks_supervised.html#neural-networks-supervised)** examples. \n",
    "\n",
    "<font color='red'>Let's try to use the **multi-layer perceptron classifier** on the Boston House Price dataset (using 75% of the data for training and 25% for testing).</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import preprocessing\n",
    "\n",
    "boston = load_boston()\n",
    "\n",
    "X = boston.data\n",
    "y = boston.target\n",
    "\n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.25, random_state=42)\n",
    "\n",
    "Xtrain_scaled = preprocessing.scale(Xtrain)\n",
    "Xtest_scaled = preprocessing.scale(Xtest)\n",
    "Xscaled = preprocessing.scale(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPRegressor\n",
    "\n",
    "clf = MLPRegressor(solver='lbfgs', alpha=1e-5, \n",
    "                   hidden_layer_sizes=(5,2), \n",
    "                   random_state=1, max_iter=5000)\n",
    "clf.fit(Xtrain_scaled, ytrain)\n",
    "\n",
    "# Look at the weights\n",
    "print([coef.shape for coef in clf.coefs_])\n",
    "\n",
    "ypred = clf.predict(Xtest_scaled)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "plt.scatter(ytest,ypred)\n",
    "plt.xlabel(\"Actual Value [x$1000]\")\n",
    "plt.ylabel(\"Predicted Value [x$1000]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, that only predicts the value for a fraction of the data set. <font color='red'>Again, we can use Scikit-Learn's [cross_val_predict](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_predict.html#sklearn.model_selection.cross_val_predict) to make predictions for the full data set.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_predict\n",
    "\n",
    "yCVpred = cross_val_predict(clf, ____, ____, cv=____) # Complete\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "plt.scatter(y,yCVpred)\n",
    "plt.xlabel(\"Actual Value [x$1000]\")\n",
    "plt.ylabel(\"Predicted Value [x$1000]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as with other model hyperparameters, we can use cross-validation to determine the optimal number of layers, neurons per layer, etc. But for now, let's talk about some guidelines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Number of Layers\n",
    "\n",
    "*Ivezic*: \n",
    "\n",
    "> \"*For data that can be represented by a linear model, no layers are required (McCullagh & Nelder 1989).  A single layer network can approximate any continuous function.  Two layers can represent arbitrary decision boundaries for smooth functions (Lippmann 1987).  More layers can represent non-continuous or complex structure within the data.*\"  \n",
    "\n",
    "Think about why no layers are needed for linear regression.  We just connect our input to the output where the synapses are the weights (slopes) and the output neurons add the constant (intecept).\n",
    "\n",
    "So you might start with a single layer, then add more layers and use cross-validation to determine when you are overfitting.\n",
    "\n",
    "Say you had to draw a whole forest, but you couldn't cut and paste anything.  That would be very tedious.  But if you could draw just one leaf and copy that to make a small branch, then scale up a small branch to a big branch and copy that, then attach the branches to a tree trunk and then make copies of the full tree, you wouldn't have to draw all that much. **Each of the layers in a neural network handles more and more detailed aspects of the problem.**  \n",
    "\n",
    "For image recognition, you might need dozens of layers, but also a huge training set to populate those layers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Number of Neurons\n",
    "\n",
    "The number of neurons in each layer is also a free parameter. \n",
    "\n",
    "- ***Typically choose somewhere between twice the number of input nodes and a number between the number of input and output nodes.***\n",
    "\n",
    "\n",
    "- If there are lots hidden layers (where \"lots\" is not clearly defined) then we call that a **deep neural network or [deep learning](https://en.wikipedia.org/wiki/Deep_learning)**. \n",
    "\n",
    "\n",
    "- Sometimes the number of neurons in each layer goes down.  But it can also be useful to have the same number in each layer so that there is only one hyperparameter (the number of neurons) and not one per layer.\n",
    "\n",
    "\n",
    "- In practice a reasonable approach is to simply **specify many more layers and neurons than you need and perform regularization**. This can be as simple as just stopping the training when the cross-validation error reaches a minimum, which appropriately (for once) is called **[early stopping](https://en.wikipedia.org/wiki/Early_stopping)**. Basically, you put your `fit` method into a loop and instantiate with `max_iter=1` and `warm_start=True`.\n",
    "\n",
    "\n",
    "While the number of neurons in the hidden layers are free parameters the number of input and output nodes are constrained by the data and the desired output.  For example, the MNIST digits data requires 784 input neurons (one for each pixel in the 28x28 images) and 10 output neurons (one for each class [digit]).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Activation Functions\n",
    "\n",
    "The **[Activation function](https://en.wikipedia.org/wiki/Activation_function)** controls how much \"signal\" it takes for a neuron to \"fire\".  The more signal, the more likely the neuron will fire. See https://mlfromscratch.com/activation-functions-explained/#/\n",
    "\n",
    "The cells below show different activation functions, using the same visualization as we used for loss functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mathematical formulas for activation functions\n",
    "  \n",
    "def binary(raw_model_output):\n",
    "    return np.where(raw_model_output < 0, \n",
    "                    0, \n",
    "                    1)  \n",
    "\n",
    "def sigmoid(raw_model_output):\n",
    "    return 1.0 / (1 + np.exp(-raw_model_output))\n",
    "\n",
    "def ReLU(raw_model_output):\n",
    "    return np.where(raw_model_output < 0, \n",
    "                    0, \n",
    "                    raw_model_output)\n",
    "\n",
    "def LReLU(raw_model_output):\n",
    "    alpha=0.1\n",
    "    return np.where(raw_model_output < 0, \n",
    "                    alpha*raw_model_output, \n",
    "                    raw_model_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid of values and plot\n",
    "grid = np.linspace(-3,3,1000)\n",
    "#plt.plot(grid, log_loss(grid), label='logistic')\n",
    "plt.plot(grid, binary(grid), \"k\", label='binary')\n",
    "plt.plot(grid, sigmoid(grid), \"b\", label='sigmoid')\n",
    "#plt.plot(grid, ReLU(grid), label='ReLU')\n",
    "#plt.plot(grid, LReLU(grid), label='LReLU')\n",
    "#plt.plot(grid, l2(grid), label='L2')\n",
    "#plt.plot(grid, l1(grid), label='L1')\n",
    "\n",
    "plt.fill_between([0,3], y1=-1, y2=3, \n",
    "                 color=\"b\", alpha=0.2)\n",
    "plt.fill_between([-3,0], y1=-1, y2=3, \n",
    "                 color=\"r\", alpha=0.2)\n",
    "plt.xlim([-3,3])\n",
    "plt.ylim([-1,3])\n",
    "plt.xlabel(\"z\",fontsize=12)\n",
    "plt.title(\"don't fire                 fire\",fontsize=12)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"*On vs. Off*\" activation actually isn't quite true or what we want. The way your eyes work is that you need 1-10 photons to trigger a rod, but several rods must be triggered to send a signal to the brain. The sigmoid activation function captures the probabilistic nature of neuron firing. More importantly it is differentiable, so can be used for backpropagation.\n",
    "\n",
    "Another important aspect of non-linear activation functions is that they are what allow neural networks to solve non-linear problems.  If we used a strictly linear activation function, then we could only solve linear problems.  That is, you could fit a straight line, but not an exponential.\n",
    "\n",
    "\n",
    "#### Vanishing and Exploding Gradients\n",
    "\n",
    "However, neural network research suffered significant limitations and problems at the hands of the [vanishing and exploding gradients](https://en.wikipedia.org/wiki/Vanishing_gradient_problem) problem.  We won't go into detail there except to say that ***around 2010 there were suggestions for different activation functions that were published.***  \n",
    "\n",
    "- For example, the **[Rectified Linear Unit (ReLU) activation function](https://www.wikiwand.com/en/Rectifier_(neural_networks))** is another commonly used activation function as it solves the vanishing gradient problem (since the gradient is only 0 or 1).\n",
    "$${\\rm ReLU}(z) = max(0,z)$$\n",
    "\n",
    "\n",
    "- It isn't ideal both because the derivative is 0 for $z<0$ and it can end up producing dead nodes. However, it is fast and the resulting sparsity can be good (sort of like regularization).\n",
    "\n",
    "\n",
    "- A number of papers since 2015 describe various improvements, including the **Leaky ReLU**, the **exponential linear unit (ELU)**, and **scaled exponential linear unit (SELU)**.\n",
    "\n",
    "\n",
    "Note that the activation functions can be different in different layers.  For example, for regression, one typically doesn't use any activation function in the output layer as including one would restrict the range of possible outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid of values and plot\n",
    "grid = np.linspace(-3,3,1000)\n",
    "#plt.plot(grid, log_loss(grid), label='logistic')\n",
    "#plt.plot(grid, binary(grid), \"k\", label='binary')\n",
    "#plt.plot(grid, sigmoid(grid), \"b\", label='sigmoid')\n",
    "plt.plot(grid, ReLU(grid), \"b--\", label='ReLU')\n",
    "plt.plot(grid, LReLU(grid), \"orange\", label='LReLU')\n",
    "#plt.plot(grid, l2(grid), label='L2')\n",
    "#plt.plot(grid, l1(grid), label='L1')\n",
    "\n",
    "plt.fill_between([0,3], y1=-1, y2=3, \n",
    "                 color=\"b\", alpha=0.2)\n",
    "plt.fill_between([-3,0], y1=-1, y2=3, \n",
    "                 color=\"r\", alpha=0.2)\n",
    "plt.xlim([-3,3])\n",
    "plt.ylim([-1,3])\n",
    "plt.xlabel(\"z\",fontsize=12)\n",
    "plt.title(\"don't fire                       fire\",fontsize=12)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some general guidance on activation functions:\n",
    "    \n",
    "* **Use sigmoid for output of binary classification (with binary cross entropy loss)**\n",
    "\n",
    "\n",
    "* **Use ReLU for layers other than output (at least to start with because it is faster)**\n",
    "\n",
    "\n",
    "* **Use softmax for output with more than 2 classes (with categorical cross entropy loss)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularization\n",
    "\n",
    "Just as we can use regularization for standard regression and classification tasks, so too can we with neural networks.  \n",
    "\n",
    "Not only can we apply the usual $L1$ (LASSO) or $L2$ (Ridge) regularization techniques, we can also use **dropout** which, as the name indicates, causes some neurons to be temporarily \"dropped out\" during training (usually by setting some probability for that to happen, typically 10-50%). After training, all of the neurons are used.\n",
    "\n",
    "*Geron* explains this in terms of a company needing to try to figure out how to adapt to a crucial employee being out sick for a period of time.  In the end, it can make the company stronger as more people (neurons) are able to handle certain parts of the process.\n",
    "\n",
    "\n",
    "### Batch normalization\n",
    "\n",
    "Just as it is often necessary to normalize or standardize our features, sometimes it is helpful to do the same to the output of the hidden layers.  This is called **[batch normalization](https://en.wikipedia.org/wiki/Batch_normalization)** and is done before passing the data to the activation function.  It make the process more stable and can also make it faster.  We'll do an example next time. \n",
    "\n",
    "\n",
    "### Faster Optimizers\n",
    "\n",
    "We aren't going to talk about optimizers, but it might be useful for you to have some options to feed into a search for the best parameter using cross validation.  For example **`['mse', 'adam', 'sgd', 'adagrad']`**."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}