02_Training_linear_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Chapter 4 – Training Linear Models**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You only need to put code in cells just below where a task is stated.\n",
    "The rest is setup and plotting. Don't worry if you don't understand what those big chunks of code are doing!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/UAPH4582/PH482_582/blob/main/02_Training_linear_models.ipynb\"><img src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import a few common modules, ensure MatplotLib plots figures inline and prepare a function to save the figures. We also check that Python 3.5 or later is installed (although Python 2.x may work, it is deprecated so we strongly recommend you use Python 3 instead), as well as Scikit-Learn ≥0.20."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python ≥3.5 is required\n",
    "import sys\n",
    "assert sys.version_info >= (3, 5)\n",
    "\n",
    "# Scikit-Learn ≥0.20 is required\n",
    "import sklearn\n",
    "assert sklearn.__version__ >= \"0.20\"\n",
    "\n",
    "# Common imports\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# to make this notebook's output stable across runs\n",
    "np.random.seed(42)\n",
    "\n",
    "# To plot pretty figures\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "\n",
    "# Where to save the figures\n",
    "PROJECT_ROOT_DIR = \".\"\n",
    "CHAPTER_ID = \"training_linear_models\"\n",
    "IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, \"images\", CHAPTER_ID)\n",
    "os.makedirs(IMAGES_PATH, exist_ok=True)\n",
    "\n",
    "def save_fig(fig_id, tight_layout=True, fig_extension=\"png\", resolution=300):\n",
    "    path = os.path.join(IMAGES_PATH, fig_id + \".\" + fig_extension)\n",
    "    print(\"Saving figure\", fig_id)\n",
    "    if tight_layout:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig(path, format=fig_extension, dpi=resolution)\n",
    "\n",
    "# Ignore useless warnings (see SciPy issue #5998)\n",
    "import warnings\n",
    "warnings.filterwarnings(action=\"ignore\", message=\"^internal gelsd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = 2 * np.random.rand(100, 1)\n",
    "y = 4 + 3 * X + np.random.randn(100, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X, y, \"b.\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "save_fig(\"generated_data_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1\n",
    "Fit a linear regression model and print out its intercept and coefficient. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lin_reg = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2 \n",
    "With the model you created above, predict the new values `X_new`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new = np.array([[0], [2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# y_new = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the data (blue) and our predictions (red)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(X, y, \"b.\")\n",
    "plt.plot(X_new, y_new, \"r.\", markersize = 20)\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "save_fig(\"generated_data_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3\n",
    "The `LinearRegression` class is based on the `scipy.linalg.lstsq()` function (the name stands for \"least squares\").\n",
    "\n",
    "Basically we fit to the form y = 1 theta0 + theta1 x = x_b.theta with x_b = (1, x)^T and theta = (theta0, theta1)^T. So we need to add a 1 to each instance of `X`. See slide 13 of Lecture 4.\n",
    "\n",
    "Look up `scipy.linalg.lstsq()`, call it directly and print out the output. Compare the output with the output from your `LinearRegression` model above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_b = np.c_[np.ones((100, 1)), X]  # add x0 = 1 to each instance\n",
    "X_new_b = np.c_[np.ones((2, 1)), X_new]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# theta_best_svd, residuals, rank, s = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression using batch gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just click through this section.\n",
    "\n",
    "The first part shows how to manually do gradient descent. The formula for the gradients has been calculated by hand and we just plug in `X` and `theta`.\n",
    "\n",
    "The second part (function `plot_gradient_descent`) does the gradient descent and plots the first 10 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta = 0.1  # learning rate\n",
    "n_iterations = 1000\n",
    "m = 100\n",
    "\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "\n",
    "for iteration in range(n_iterations):\n",
    "    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)\n",
    "    theta = theta - eta * gradients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new_b.dot(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "    for iteration in range(n_iterations):\n",
    "        if iteration < 10:\n",
    "            y_predict = X_new_b.dot(theta)\n",
    "            style = \"b-\" if iteration > 0 else \"r--\"\n",
    "            plt.plot(X_new, y_predict, style)\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",
    "    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": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "theta_path_bgd = []\n",
    "\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.subplot(131); plot_gradient_descent(theta, eta=0.02)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)\n",
    "plt.subplot(133); plot_gradient_descent(theta, eta=0.5)\n",
    "\n",
    "save_fig(\"gradient_descent_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section shows how one can implement Stochastic Gradient Descent by hand.\n",
    "\n",
    "Try to understand the part that is not marked with `#not shown`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_sgd = []\n",
    "m = len(X_b)\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 50\n",
    "t0, t1 = 5, 50  # learning schedule hyperparameters\n",
    "\n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)\n",
    "\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    for i in range(m):\n",
    "        if epoch == 0 and i < 20:                    # not shown in the book\n",
    "            y_predict = X_new_b.dot(theta)           # not shown\n",
    "            style = \"b-\" if i > 0 else \"r--\"         # not shown\n",
    "            plt.plot(X_new, y_predict, style)        # not shown\n",
    "        random_index = np.random.randint(m)\n",
    "        xi = X_b[random_index:random_index+1]\n",
    "        yi = y[random_index:random_index+1]\n",
    "        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(epoch * m + i)\n",
    "        theta = theta - eta * gradients\n",
    "        theta_path_sgd.append(theta)                 # not shown\n",
    "\n",
    "plt.plot(X, y, \"b.\")                                 # not shown\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)                     # not shown\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)           # not shown\n",
    "plt.axis([0, 2, 0, 15])                              # not shown\n",
    "save_fig(\"sgd_plot\")                                 # not shown\n",
    "plt.show()                                           # not shown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4\n",
    "\n",
    "Fit an SGD Regressor to `X` and `y` and print out its intercept and coefficients. Compare the values you find to the ones from the stochastic gradient descent above.\n",
    "\n",
    "Use `max_iter=50, tol=-np.infty, penalty=None, eta0=0.1, random_state=42` as parameters for the SGD Regressor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import SGDRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sgd_reg = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mini-batch gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code shows how one could implement mini-batch gradient descent by hand. Can you find the difference to Stochastic gradient descent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_mgd = []\n",
    "\n",
    "n_iterations = 50\n",
    "minibatch_size = 20\n",
    "\n",
    "np.random.seed(42)\n",
    "theta = np.random.randn(2,1)  # random initialization\n",
    "\n",
    "t0, t1 = 200, 1000\n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)\n",
    "\n",
    "t = 0\n",
    "for epoch in range(n_iterations):\n",
    "    shuffled_indices = np.random.permutation(m)\n",
    "    X_b_shuffled = X_b[shuffled_indices]\n",
    "    y_shuffled = y[shuffled_indices]\n",
    "    for i in range(0, m, minibatch_size):\n",
    "        t += 1\n",
    "        xi = X_b_shuffled[i:i+minibatch_size]\n",
    "        yi = y_shuffled[i:i+minibatch_size]\n",
    "        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(t)\n",
    "        theta = theta - eta * gradients\n",
    "        theta_path_mgd.append(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_bgd = np.array(theta_path_bgd)\n",
    "theta_path_sgd = np.array(theta_path_sgd)\n",
    "theta_path_mgd = np.array(theta_path_mgd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], \"r-s\", linewidth=1, label=\"Stochastic\")\n",
    "plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], \"g-+\", linewidth=2, label=\"Mini-batch\")\n",
    "plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], \"b-o\", linewidth=3, label=\"Batch\")\n",
    "plt.legend(loc=\"upper left\", fontsize=16)\n",
    "plt.xlabel(r\"$\\theta_0$\", fontsize=20)\n",
    "plt.ylabel(r\"$\\theta_1$   \", fontsize=20, rotation=0)\n",
    "plt.axis([2.5, 4.5, 2.3, 3.9])\n",
    "save_fig(\"gradient_descent_paths_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 5\n",
    "Explain which Linear Regression training algorithm you can use if you have a training set with millions of features? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a little trick Linear Regression turns out to be rather powerful also for data that is not linear!\n",
    "The trick is to add \"new features\" to `X`. In our case we need a second feature that is just x^2. Then we fit again, but because our feature vector is now (1, x, x^2), the dot product with (theta0, theta1, theta2) gives the form for the parabola (check this if you don't see it!).\n",
    "\n",
    "Here we show an example of how to fit a parabola with Linear Regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "m = 100\n",
    "X = 6 * np.random.rand(m, 1) - 3\n",
    "y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)\n",
    "\n",
    "poly_features = PolynomialFeatures(degree=2, include_bias=False)\n",
    "X_poly = poly_features.fit_transform(X)\n",
    "\n",
    "lin_reg = LinearRegression()\n",
    "lin_reg.fit(X_poly, y)\n",
    "\n",
    "X_new=np.linspace(-3, 3, 100).reshape(100, 1)\n",
    "X_new_poly = poly_features.transform(X_new)\n",
    "y_new = lin_reg.predict(X_new_poly)\n",
    "\n",
    "for style, width, degree in ((\"g-\", 1, 300), (\"b--\", 2, 2), (\"r-+\", 2, 1)):\n",
    "    polybig_features = PolynomialFeatures(degree=degree, include_bias=False)\n",
    "    std_scaler = StandardScaler()\n",
    "    lin_reg = LinearRegression()\n",
    "    polynomial_regression = Pipeline([\n",
    "            (\"poly_features\", polybig_features),\n",
    "            (\"std_scaler\", std_scaler),\n",
    "            (\"lin_reg\", lin_reg),\n",
    "        ])\n",
    "    polynomial_regression.fit(X, y)\n",
    "    y_newbig = polynomial_regression.predict(X_new)\n",
    "    plt.plot(X_new, y_newbig, style, label=str(degree), linewidth=width)\n",
    "\n",
    "plt.plot(X, y, \"b.\", linewidth=3)\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([-3, 3, 0, 10])\n",
    "save_fig(\"high_degree_polynomials_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training and Validation Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "def plot_learning_curves(model, X, y):\n",
    "    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=10)\n",
    "    train_errors, val_errors = [], []\n",
    "    for m in range(1, len(X_train)):\n",
    "        model.fit(X_train[:m], y_train[:m])\n",
    "        y_train_predict = model.predict(X_train[:m])\n",
    "        y_val_predict = model.predict(X_val)\n",
    "        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))\n",
    "        val_errors.append(mean_squared_error(y_val, y_val_predict))\n",
    "\n",
    "    plt.plot(np.sqrt(train_errors), \"r-+\", linewidth=2, label=\"train\")\n",
    "    plt.plot(np.sqrt(val_errors), \"b-\", linewidth=3, label=\"val\")\n",
    "    plt.legend(loc=\"upper right\", fontsize=14)   \n",
    "    plt.xlabel(\"Training set size\", fontsize=14) \n",
    "    plt.ylabel(\"RMSE\", fontsize=14)              "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare different degrees of polynomials for fitting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linear Regression\n",
    "lin_reg = LinearRegression()\n",
    "plot_learning_curves(lin_reg, X, y)\n",
    "plt.axis([0, 80, 0, 3])                        \n",
    "plt.show()                                     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "# Degree 2\n",
    "polynomial_regression = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=2, include_bias=False)),\n",
    "        (\"lin_reg\", LinearRegression()),\n",
    "    ])\n",
    "\n",
    "plot_learning_curves(polynomial_regression, X, y)\n",
    "plt.axis([0, 80, 0, 3])           \n",
    "plt.show()                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "# Degree 10\n",
    "polynomial_regression = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=10, include_bias=False)),\n",
    "        (\"lin_reg\", LinearRegression()),\n",
    "    ])\n",
    "\n",
    "plot_learning_curves(polynomial_regression, X, y)\n",
    "plt.axis([0, 80, 0, 3])           \n",
    "plt.show()                        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 6\n",
    "Look at the three plots above that show the training and validation errors for linear regression, quadratic regression and order 10. Compare them to the plot above (`high_degree_polynomials_plot`).\n",
    "\n",
    "Describe which one is likely overfitting and which one is likely underfitting and why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularized models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "m = 20\n",
    "X = 3 * np.random.rand(m, 1)\n",
    "y = 1 + 0.5 * X + np.random.randn(m, 1) / 1.5\n",
    "X_new = np.linspace(0, 3, 100).reshape(100, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 7\n",
    "Fit a ridge linear model with `alpha=1` and Stochastic Average Gradient descent (`\"sag\"`) solver. Check which value it predicts for `X=[[1.5]]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ridge_reg_sag = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we fit a linear model (left plot) and a degree 10 polynomial model (right plot) with ridge regression for different values of alpha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "def plot_model(model_class, polynomial, alphas, **model_kargs):\n",
    "    for alpha, style in zip(alphas, (\"b-\", \"g--\", \"r:\")):\n",
    "        model = model_class(alpha, **model_kargs) if alpha > 0 else LinearRegression()\n",
    "        if polynomial:\n",
    "            model = Pipeline([\n",
    "                    (\"poly_features\", PolynomialFeatures(degree=10, include_bias=False)),\n",
    "                    (\"std_scaler\", StandardScaler()),\n",
    "                    (\"regul_reg\", model),\n",
    "                ])\n",
    "        model.fit(X, y)\n",
    "        y_new_regul = model.predict(X_new)\n",
    "        lw = 2 if alpha > 0 else 1\n",
    "        plt.plot(X_new, y_new_regul, style, linewidth=lw, label=r\"$\\alpha = {}$\".format(alpha))\n",
    "    plt.plot(X, y, \"b.\", linewidth=3)\n",
    "    plt.legend(loc=\"upper left\", fontsize=15)\n",
    "    plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "    plt.axis([0, 3, 0, 4])\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "plot_model(Ridge, polynomial=False, alphas=(0, 10, 100), random_state=42)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.subplot(122)\n",
    "plot_model(Ridge, polynomial=True, alphas=(0, 10**-5, 1), random_state=42)\n",
    "\n",
    "save_fig(\"ridge_regression_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: to be future-proof, we set `max_iter=1000` and `tol=1e-3` because these will be the default values in Scikit-Learn 0.21."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sgd_reg = SGDRegressor(penalty=\"l2\", max_iter=1000, tol=1e-3, random_state=42)\n",
    "sgd_reg.fit(X, y.ravel())\n",
    "sgd_reg.predict([[1.5]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same, but with Lasso:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "plot_model(Lasso, polynomial=False, alphas=(0, 0.1, 1), random_state=42)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.subplot(122)\n",
    "plot_model(Lasso, polynomial=True, alphas=(0, 10**-7, 1), random_state=42)\n",
    "\n",
    "save_fig(\"lasso_regression_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 8 \n",
    "Fit a Lasso linear model with alpha=1. Check which value it predicts for `X=[[1.5]]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lasso_reg = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Early stopping example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "m = 100\n",
    "X = 6 * np.random.rand(m, 1) - 3\n",
    "y = 2 + X + 0.5 * X**2 + np.random.randn(m, 1)\n",
    "\n",
    "X_train, X_val, y_train, y_val = train_test_split(X[:50], y[:50].ravel(), test_size=0.5, random_state=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "\n",
    "poly_scaler = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=90, include_bias=False)),\n",
    "        (\"std_scaler\", StandardScaler())\n",
    "    ])\n",
    "\n",
    "X_train_poly_scaled = poly_scaler.fit_transform(X_train)\n",
    "X_val_poly_scaled = poly_scaler.transform(X_val)\n",
    "\n",
    "sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True,\n",
    "                       penalty=None, learning_rate=\"constant\", eta0=0.0005, random_state=42)\n",
    "\n",
    "minimum_val_error = float(\"inf\")\n",
    "best_epoch = None\n",
    "best_model = None\n",
    "for epoch in range(1000):\n",
    "    sgd_reg.fit(X_train_poly_scaled, y_train)  # continues where it left off\n",
    "    y_val_predict = sgd_reg.predict(X_val_poly_scaled)\n",
    "    val_error = mean_squared_error(y_val, y_val_predict)\n",
    "    if val_error < minimum_val_error:\n",
    "        minimum_val_error = val_error\n",
    "        best_epoch = epoch\n",
    "        best_model = deepcopy(sgd_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True,\n",
    "                       penalty=None, learning_rate=\"constant\", eta0=0.0005, random_state=42)\n",
    "\n",
    "n_epochs = 500\n",
    "train_errors, val_errors = [], []\n",
    "for epoch in range(n_epochs):\n",
    "    sgd_reg.fit(X_train_poly_scaled, y_train)\n",
    "    y_train_predict = sgd_reg.predict(X_train_poly_scaled)\n",
    "    y_val_predict = sgd_reg.predict(X_val_poly_scaled)\n",
    "    train_errors.append(mean_squared_error(y_train, y_train_predict))\n",
    "    val_errors.append(mean_squared_error(y_val, y_val_predict))\n",
    "\n",
    "best_epoch = np.argmin(val_errors)\n",
    "best_val_rmse = np.sqrt(val_errors[best_epoch])\n",
    "\n",
    "plt.annotate('Best model',\n",
    "             xy=(best_epoch, best_val_rmse),\n",
    "             xytext=(best_epoch, best_val_rmse + 1),\n",
    "             ha=\"center\",\n",
    "             arrowprops=dict(facecolor='black', shrink=0.05),\n",
    "             fontsize=16,\n",
    "            )\n",
    "\n",
    "best_val_rmse -= 0.03  # just to make the graph look better\n",
    "plt.plot([0, n_epochs], [best_val_rmse, best_val_rmse], \"k:\", linewidth=2)\n",
    "plt.plot(np.sqrt(val_errors), \"b-\", linewidth=3, label=\"Validation set\")\n",
    "plt.plot(np.sqrt(train_errors), \"r--\", linewidth=2, label=\"Training set\")\n",
    "plt.legend(loc=\"upper right\", fontsize=14)\n",
    "plt.xlabel(\"Epoch\", fontsize=14)\n",
    "plt.ylabel(\"RMSE\", fontsize=14)\n",
    "save_fig(\"early_stopping_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 9\n",
    "Suppose you use Batch Gradient Descent and you plot the validation error at every epoch. If you notice that the validation error consistently goes up, what is likely going on? How can you fix this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 10\n",
    "We want to find at which exact epoch our model is the best.\n",
    "For that we go through 1000 \"epochs\".\n",
    "\n",
    "We have already prepared part of the code for this task. You only have to fill in the last part which should do the following:\n",
    "\n",
    "- if the validation error (of the current epoch) is smaller than the smalles validation error until now (`minimum_val_error`), then:\n",
    "    - set the `minimum_val_error` to be the current validation error\n",
    "    - set the `best_epoch` to be the current epoch\n",
    "    - set the `best_model` to be the current model (`clone(cgd_reg)`)\n",
    "- else: next epoch (you can also just leave out `else`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import clone\n",
    "sgd_reg = SGDRegressor(max_iter=1, tol=-np.infty, warm_start=True, penalty=None,\n",
    "                       learning_rate=\"constant\", eta0=0.0005, random_state=42)\n",
    "\n",
    "minimum_val_error = float(\"inf\")\n",
    "best_epoch = None\n",
    "best_model = None\n",
    "for epoch in range(1000):\n",
    "    sgd_reg.fit(X_train_poly_scaled, y_train)  # continues where it left off\n",
    "    y_val_predict = sgd_reg.predict(X_val_poly_scaled)\n",
    "    val_error = mean_squared_error(y_val, y_val_predict)\n",
    "    #### put your code here \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    ######################################## \n",
    "    \n",
    "print(\"Best epoch: \", best_epoch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(-10, 10, 100)\n",
    "sig = 1 / (1 + np.exp(-t))\n",
    "plt.figure(figsize=(9, 3))\n",
    "plt.plot([-10, 10], [0, 0], \"k-\")\n",
    "plt.plot([-10, 10], [0.5, 0.5], \"k:\")\n",
    "plt.plot([-10, 10], [1, 1], \"k:\")\n",
    "plt.plot([0, 0], [-1.1, 1.1], \"k-\")\n",
    "plt.plot(t, sig, \"b-\", linewidth=2, label=r\"$\\sigma(t) = \\frac{1}{1 + e^{-t}}$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.legend(loc=\"upper left\", fontsize=20)\n",
    "plt.axis([-10, 10, -0.1, 1.1])\n",
    "save_fig(\"logistic_function_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()\n",
    "list(iris.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(iris.DESCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = iris[\"data\"][:, 3:]  # petal width\n",
    "y = (iris[\"target\"] == 2).astype(np.int)  # 1 if Iris virginica, else 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: To be future-proof we set `solver=\"lbfgs\"` since this will be the default value in Scikit-Learn 0.22."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "log_reg = LogisticRegression(solver=\"lbfgs\", random_state=42)\n",
    "log_reg.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new = np.linspace(0, 3, 1000).reshape(-1, 1)\n",
    "y_proba = log_reg.predict_proba(X_new)\n",
    "decision_boundary = X_new[y_proba[:, 1] >= 0.5][0]\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(X[y==0], y[y==0], \"bs\")\n",
    "plt.plot(X[y==1], y[y==1], \"g^\")\n",
    "plt.plot([decision_boundary, decision_boundary], [-1, 2], \"k:\", linewidth=2)\n",
    "plt.plot(X_new, y_proba[:, 1], \"g-\", linewidth=2, label=\"Iris virginica\")\n",
    "plt.plot(X_new, y_proba[:, 0], \"b--\", linewidth=2, label=\"Not Iris virginica\")\n",
    "plt.text(decision_boundary+0.02, 0.15, \"Decision  boundary\", fontsize=14, color=\"k\", ha=\"center\")\n",
    "plt.arrow(decision_boundary, 0.08, -0.3, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')\n",
    "plt.arrow(decision_boundary, 0.92, 0.3, 0, head_width=0.05, head_length=0.1, fc='g', ec='g')\n",
    "plt.xlabel(\"Petal width (cm)\", fontsize=14)\n",
    "plt.ylabel(\"Probability\", fontsize=14)\n",
    "plt.legend(loc=\"center left\", fontsize=14)\n",
    "plt.axis([0, 3, -0.02, 1.02])\n",
    "save_fig(\"logistic_regression_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decision_boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_reg.predict([[1.7], [1.5]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = iris[\"data\"][:, (2, 3)]  # petal length, petal width\n",
    "y = iris[\"target\"]\n",
    "\n",
    "softmax_reg = LogisticRegression(multi_class=\"multinomial\",solver=\"lbfgs\", C=10, random_state=42)\n",
    "softmax_reg.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0, x1 = np.meshgrid(\n",
    "        np.linspace(0, 8, 500).reshape(-1, 1),\n",
    "        np.linspace(0, 3.5, 200).reshape(-1, 1),\n",
    "    )\n",
    "X_new = np.c_[x0.ravel(), x1.ravel()]\n",
    "\n",
    "\n",
    "y_proba = softmax_reg.predict_proba(X_new)\n",
    "y_predict = softmax_reg.predict(X_new)\n",
    "\n",
    "zz1 = y_proba[:, 1].reshape(x0.shape)\n",
    "zz = y_predict.reshape(x0.shape)\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(X[y==2, 0], X[y==2, 1], \"g^\", label=\"Iris virginica\")\n",
    "plt.plot(X[y==1, 0], X[y==1, 1], \"bs\", label=\"Iris versicolor\")\n",
    "plt.plot(X[y==0, 0], X[y==0, 1], \"yo\", label=\"Iris setosa\")\n",
    "\n",
    "from matplotlib.colors import ListedColormap\n",
    "custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])\n",
    "\n",
    "plt.contourf(x0, x1, zz, cmap=custom_cmap)\n",
    "contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)\n",
    "plt.clabel(contour, inline=1, fontsize=12)\n",
    "plt.xlabel(\"Petal length\", fontsize=14)\n",
    "plt.ylabel(\"Petal width\", fontsize=14)\n",
    "plt.legend(loc=\"center left\", fontsize=14)\n",
    "plt.axis([0, 7, 0, 3.5])\n",
    "save_fig(\"softmax_regression_contour_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "softmax_reg.predict([[5, 2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "softmax_reg.predict_proba([[5, 2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 11\n",
    "Can Gradient Descent get stuck in a local minimum when training a Logistic Regression model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 12\n",
    "Which Gradient Descent algorithm (among those we discussed) will reach the vicinity of the optimal solution the fastest? Which will actually converge? How can you make the others converge as well?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}