week_5/5.2_learning_better_pt.2/learning_better_pt.2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 5.2: Learning Better Pt.2\n",
    "\n",
    "This lecture, we are going to regularize a model to improve its generalization.\n",
    "\n",
    "**Learning goals:**\n",
    "\n",
    "- overfit a polynomial regression model\n",
    "- use ridge regularization to prevent overfitting\n",
    "- analyse the effect of the regularization weight on model parameter sizes\n",
    "- train a regularized logistic regression model with polynomial features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This is part 2 of \"Learning Better\". We'll be using the same dataset and trying to solve the same task as notebook 5.1. If you haven't been through notebook 5.1, please start there."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. High Polynomial Regression\n",
    "\n",
    "Let's load our instagram planning dataset from csv into a `DataFrame`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('instagram_planning.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To refresh our minds since last lecture, let's visualize the dataset again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.scatter(df['planned_minutes'], df['actual_minutes'], alpha=0.6)\n",
    "\n",
    "ax.set_xlabel('planned online time (min)')\n",
    "ax.set_ylabel('time spent online (min) ')\n",
    "ax.set_title('Online Procrastination');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember the single feature is the number of minutes we _plan_ to spend online, and the numerical label is the _actual_ amount of time we spend online. 📱\n",
    "\n",
    "In notebook 5.1, we had success adding polynomial features to our linear regression model. The resulting non-linear hypothesis was a better fit to the labels. We are therefore tempted to power up the model even more... let's turn it up all the way to 11! 🎸\n",
    "\n",
    "Let's create our feature matrix and label vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df['planned_minutes'].values.reshape(-1, 1)\n",
    "y = df['actual_minutes'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then add the polynomial features of degree 11, and standardize them using sklearn proprocessors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures, StandardScaler\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "poly = PolynomialFeatures(11, include_bias=False)\n",
    "poly = poly.fit(X)\n",
    "X_poly = poly.transform(X)\n",
    "\n",
    "scaler = StandardScaler()\n",
    "scaler = scaler.fit(X_poly)\n",
    "X_poly_scaled = scaler.transform(X_poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can train our superpowered polynomial regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg = LinearRegression().fit(X_poly_scaled, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like notebook 5.1, we can visualize the optimized hypothesis by using `.predict()`. Let's put this code in a function, so that we may reuse it later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def plot_polynomial_regression(X, y, reg, poly, scaler):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "\n",
    "    # plot the examples\n",
    "    ax.scatter(X, y, alpha=0.6)\n",
    "\n",
    "    # create feature matrix\n",
    "    xmin, xmax = ax.get_xlim()\n",
    "    x_line = np.linspace(xmin, xmax, 30).reshape(-1, 1)\n",
    "\n",
    "    # apply feature transformations\n",
    "    x_line_poly = poly.transform(x_line)\n",
    "    x_line_poly_scaled = scaler.transform(x_line_poly)\n",
    "\n",
    "    # predict\n",
    "    y_line = reg.predict(x_line_poly_scaled)\n",
    "\n",
    "    # plot the hypothesis\n",
    "    ax.plot(x_line, y_line, c='g', linewidth=3)\n",
    "\n",
    "    # formatting\n",
    "    ax.set_xlim(xmin, xmax)\n",
    "    ax.set_xlabel('planned online time (min)')\n",
    "    ax.set_ylabel('time spent online (min)')\n",
    "    ax.set_title('Online Procrastination');\n",
    "    \n",
    "plot_polynomial_regression(X, y, reg, poly, scaler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... what happened? we were expecting awesome modeling power, but this hypothesis looks all crooked! The value are shooting up for small values of `planned_minutes`, and dropping for large values of the feature, both of which make no sense. 😜\n",
    "\n",
    "This is a typical example of _overfitting_. The model is too powerful and tries to fit the noise rather than the underlying pattern.\n",
    "\n",
    "One way of preventing this from happening, is to chill out with the degree polynomial features and maybe keep it to 2 like last lecture. 💆‍♂️ However sometimes it is difficult to know which exact polynomial degree, or which model \"power\", to use. Therefore, a more robust way of preventing overfitting is _regularization_.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Ridge Regression\n",
    "\n",
    "Overfitting typically happens when some of the model parameters become very large. Large $\\theta$ suggests \"abrupt\" variations in the hypothesis, and is a sign that the model is trying to fit noise.\n",
    "\n",
    "We can check the model parameters inside the fitted `LinearRegression` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.append(reg.intercept_, reg.coef_)\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we have added polynomial features of degree 11. As expected, $\\theta$ is a vector of length 12: one $\\theta_{i}$ per feature, plus one intercept term.\n",
    "\n",
    "Some of these $\\theta$ values are _very_ large! We're seeing quite a few values in the tens of millions. There is no universal guideline about an \"appropriate\" range for $\\theta$ values, but keep in mind that our features are standardized, and our labels are in the range of hundreds. Thus, there is no reason for $\\theta_{i} > 1000000$. 🚓\n",
    "\n",
    "Ridge regularization penalizes large $\\theta_{i}$ values by adding a regularization term to the cost function. \n",
    "\n",
    "$$ \n",
    "J(\\theta) = \\frac{1}{m}\\sum_{i=1}^{m}(h_{\\theta}(x^{(i)} - y^{(i)})^{2} + \\alpha \\sum_{j=1}^{n}\\theta_{j}^{2}\n",
    "$$\n",
    "\n",
    "Hopefully this can help keep our model parameters under control!\n",
    "\n",
    "sklearn directly integrates ridge regularization into the model api. Therefore all we have to do is switch the vanilla `LinearRegression` model with a [`Ridge`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html) instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "reg = Ridge(alpha=1, random_state=0).fit(X_poly_scaled, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's that easy! sklearn takes care of choosing the correct regularized cost function, and optimizing it using gradient descent. If you need a refresher on what these steps are, check out lecture 4.2 \n",
    "\n",
    "We can now check if our model parameters are better behaved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.append(reg.intercept_, reg.coef_)\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All $\\theta_{i}$ are in the range of tens, so the ridge regularization worked! Our model is less likely to be overfit. We can verify this graphically: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_polynomial_regression(X, y, reg, poly, scaler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hypothesis no longer seems to model noisy data points. Instead it was \"smoothed\" by our ridge regularization, and we can expect it to generalize better.\n",
    "\n",
    "ℹ️ Note how the hypothesis is very similar to the polynomial regression of degree 4 we have trained in the final task of last lecture. This agrees with the fact that model parameter values beyond $\\theta_{4}$ are quite small compared to the rest. Small $\\theta_{i}$ values mean that the features don't contribute much to the overall hypothesis. This is like saying that we gave extra \"polynomial power\" to the regression model, but the optimization procedure found this extra power unnecessary. 💁‍♂️This is an indication that polynomial features of degree 11 might have been overkill here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We haven't talked about the `Ridge` argument: `alpha`. The regularization term added to the MSE cost function is multiplied by a _regularization weight_ , $\\alpha$. This allows us to control the strength of the regularization penalty. Let's investigate what is the effect of $\\alpha$ on our model parameter values $\\theta$.\n",
    "\n",
    "We'll train 100 different `Ridge` regression models, each with a different value of $\\alpha$. We can then plot the values of $\\theta$ and observe their relationship:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_ridge(alpha, X, y):\n",
    "    reg = Ridge(alpha=alpha).fit(X_poly_scaled, y)\n",
    "    return np.append(reg.intercept_, reg.coef_)\n",
    "    \n",
    "alphas = np.linspace(0.000001, 1, 100)\n",
    "thetas = [train_ridge(alpha, X, y) for alpha in alphas]\n",
    "thetas = np.array(thetas)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "for i in range(0, thetas.shape[1]):\n",
    "    theta = thetas[:, i]\n",
    "    ax.plot(alphas, theta)\n",
    "    \n",
    "ax.set_xlabel('α')\n",
    "ax.set_ylabel('θ')\n",
    "ax.set_title('Model Parameters vs Regularization Weight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each colored line is a parameter $\\theta_{i}$. We can spot a reduction in the absolute value of the model parameters, but it happens all scrunched up for very small values of $\\alpha$. Let's \"zoom in\" a little:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.linspace(0.000001, 0.0001, 100)\n",
    "thetas = [train_ridge(alpha, X, y) for alpha in alphas]\n",
    "thetas = np.array(thetas)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "for i in range(0, thetas.shape[1]):\n",
    "    theta = thetas[:, i]\n",
    "    ax.plot(alphas, theta)\n",
    "    \n",
    "ax.set_xlabel('α')\n",
    "ax.set_ylabel('θ')\n",
    "ax.set_title('Model Parameters vs Regularization Weight');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can clearly see the effect of $\\alpha$ on the model parameters. For big $\\alpha$, the effect of the regularization is stronger, and the large values of $\\theta$ are more penalized. This trend slows down and the $\\theta$ values converge around $\\alpha \\approx 0.0001$. This means that for this regression model, there is no additional regularization provided by $\\alpha > 0.0001$\n",
    "\n",
    "ℹ️ Notice that this $\\alpha$ threshold of 0.0001 is much smaller than the default 1.0 used in the [section above](#3.-Ridge-Regression). If a regularized model is _underfit_ , it can be helpful to lower $\\alpha$ in order to increase modeling power."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Ridge Classification\n",
    "\n",
    "The penalty term added to the cost function in ridge regularization isn't specific to the MSE cost function... it can also be applied to our logistic regression's cross-entropy cost function! \n",
    "\n",
    "You are going to learn how to apply ridge regularization in classification task _all on your own_. No need to worry, this is very similar to the sections above since sklearn has a model api class especially dedicated to ridge regularization for logistic regression: [`RidgeClassifier`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeClassifier.html).\n",
    "\n",
    "💪💪💪 Train a ridge regularized polynomial logistic regression classifier. And then tell your friends, because that's a badass model name. 😎Here's a few steps to help out:\n",
    "- Load the `wine_quality.csv` dataset\n",
    "- use the `sulphates` and `alcohol` columns as features, and the `tasty` column as categorical label\n",
    "- store the feature matrix and label vector in two NumPy `ndarray`s, `X` and `y`\n",
    "- use a polynomial feature preprocessor named `poly` and a feature scaler called `scaler`\n",
    "- name your `RidgeClassifier` instance `clf` (for classifier)\n",
    "- use polynomial degree 5, and default regularization weight 1.0\n",
    "\n",
    "You are only expected to _fit_ your model. You can then run the two cells below to unit test your code, and visualize the non-linear decision boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INSERT CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_decision_boundary(ax, X, poly, scaler, clf):\n",
    "    h = .02  # step size in the mesh\n",
    "    # create a mesh to plot in\n",
    "    x1_min, x1_max = X.iloc[:, 0].min() - 1, X.iloc[:, 0].max() + 1\n",
    "    x2_min, x2_max = X.iloc[:, 1].min() - 1, X.iloc[:, 1].max() + 1\n",
    "    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, h),\n",
    "                         np.arange(x2_min, x2_max, h))\n",
    "\n",
    "\n",
    "    # Plot the decision boundary. For that, we will assign a color to each\n",
    "    # point in the mesh [x_min, m_max]x[y_min, y_max].\n",
    "    X_mesh = np.c_[xx1.ravel(), xx2.ravel()]\n",
    "    X_mesh_poly = poly.transform(X_mesh)\n",
    "    X_mesh_poly_scaled = scaler.transform(X_mesh_poly)\n",
    "    Z = clf.predict(X_mesh_poly_scaled)\n",
    "\n",
    "    # Put the result into a color plot\n",
    "    Z = Z.reshape(xx1.shape)\n",
    "    ax.contour(xx1, xx2, Z, cmap=plt.cm.Greens)\n",
    "\n",
    "def plot_classification(X, y, poly, scaler, clf):\n",
    "    colors = sns.color_palette('husl').as_hex()\n",
    "\n",
    "    fig = plt.figure(dpi=120)\n",
    "    ax = fig.add_subplot(111)\n",
    "\n",
    "    # plot points\n",
    "    positives = X[y == 1]\n",
    "    negatives = X[y == 0]\n",
    "    ax.scatter(positives.iloc[:, 0], positives.iloc[:, 1], alpha=0.4, s=8,c=colors[0], label='tasty')\n",
    "    ax.scatter(negatives.iloc[:, 0], negatives.iloc[:, 1], alpha=0.4, s=8, c=colors[4], label='not tasty')\n",
    "\n",
    "\n",
    "    # plot decision boundary\n",
    "    plot_decision_boundary(ax, X, poly, scaler, clf)\n",
    "\n",
    "    # formatting\n",
    "    ax.set_ylim((-2, 4))\n",
    "    ax.set_xlim((-2, 4))\n",
    "    ax.set_xlabel('Sulphates (norm)')\n",
    "    ax.set_ylabel('Alcohol (norm)')\n",
    "    ax.legend()\n",
    "    ax.set_title('Wine Tastiness vs Sulphates & Alcohol Contents');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def test_classification():\n",
    "    X_predict = np.array([1, -1]).reshape(1, 2)\n",
    "    X_predict_poly = poly.transform(X_predict)\n",
    "    X_predict_poly_scaled = scaler.transform(X_predict_poly)\n",
    "    y_predict = clf.decision_function(X_predict_poly_scaled)\n",
    "    assert math.isclose(y_predict, -0.05693, rel_tol=1e-4)\n",
    "    print('Success! 🎉')\n",
    "\n",
    "test_classification()\n",
    "plot_classification(X, y, poly, scaler, clf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congrats on training your first model from scratch! 🎊\n",
    "\n",
    "🧠 Can you tell which sides of the decision boundary are positive vs negative predictions?  \n",
    "\n",
    "🧠 Are there graphical signs that the regularization worked?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Summary\n",
    "\n",
    "Today, we learned about **regularization**. First we defined **overfitting**, and how avoiding it promises a better model **generalization**. We then explained how to avoid overfitting, by preventing model parameters values from getting too big. We understood that this can be implemented by adding an extra term to the model's **cost-function**. We also highlighted how a data scientist must find the right **balance** between underfit and overfit models. We tested the **ridge** regularization technique on a polynomial regression model of degree 11. We then analysed the influence of the **regularization weight**, $\\alpha$, before implementing our own **regularized classifier**.\n",
    "\n",
    "# Resources\n",
    "\n",
    "## Core Resources\n",
    "\n",
    "- [Regularization in ML](https://towardsdatascience.com/regularization-in-machine-learning-76441ddcf99a)  \n",
    "Comprehensive list of cost-function regularization methods\n",
    "\n",
    "## Additional Resources\n",
    "\n",
    "- [Visualizing Ridge Regression](https://xavierbourretsicotte.github.io/animation_ridge.html)"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}