task_9_random_forests.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 9: Random Forests\n",
    "\n",
    "_All credit for the code examples of this notebook goes to the book \"Hands-On Machine Learning with Scikit-Learn & TensorFlow\" by A. Geron. Modifications were made and text was added by K. Zoch in preparation for the hands-on sessions._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup\n",
    "\n",
    "First, import a few common modules, ensure MatplotLib plots figures inline and prepare a function to save the figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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",
    "# Function to save a figure. This also decides that all output files \n",
    "# should stored in the subdirectorz 'classification'.\n",
    "PROJECT_ROOT_DIR = \".\"\n",
    "EXERCISE = \"forests\"\n",
    "\n",
    "def save_fig(fig_id, tight_layout=True):\n",
    "    path = os.path.join(PROJECT_ROOT_DIR, \"output\", EXERCISE, fig_id + \".png\")\n",
    "    print(\"Saving figure\", fig_id)\n",
    "    if tight_layout:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig(path, format='png', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bagging decision trees\n",
    "\n",
    "First, let's create some half-moon data (as done in one of the earlier tasks)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "markdown",
   "metadata": {},
   "source": [
    "This code example shows how \"bagging\" multiple decision trees can improve the classification performance, compared to a single decision tree. Notice how bias and variance change when combining 500 trees as in the example below (it can be seen very nicely in the plot). Please try the following:\n",
    "\n",
    "1. How does the number of samples affect the performance of the ensemble classifier? Try changing it to the training size (m = 500), or go even higher.\n",
    "2. How is the performance different when pasting is used instead of bagging (_no_ replacement of instances)?\n",
    "3. How relevant is the number of trees in the ensemble?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import BaggingClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "# Create an instance of a bagging classifier, composed of\n",
    "# 500 decision tree classifiers. bootstrap=True activates\n",
    "# replacement when picking the random instances, i.e.\n",
    "# turning it off will switch from bagging to pasting.\n",
    "bag_clf = BaggingClassifier(\n",
    "    DecisionTreeClassifier(random_state=42), n_estimators=500,\n",
    "    max_samples=100, bootstrap=True, n_jobs=-1, random_state=42)\n",
    "bag_clf.fit(X_train, y_train)\n",
    "y_pred = bag_clf.predict(X_test)\n",
    "\n",
    "# Create an instance of a single decision tree to compare with.\n",
    "tree_clf = DecisionTreeClassifier(random_state=42)\n",
    "tree_clf.fit(X_train, y_train)\n",
    "y_pred_tree = tree_clf.predict(X_test)\n",
    "\n",
    "# Now do the plotting of the two.\n",
    "from matplotlib.colors import ListedColormap\n",
    "def plot_decision_boundary(clf, X, y, axes=[-1.5, 2.5, -1, 1.5], 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",
    "    X_new = np.c_[x1.ravel(), x2.ravel()]\n",
    "    y_pred = clf.predict(X_new).reshape(x1.shape)\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",
    "    plt.axis(axes)\n",
    "    plt.xlabel(r\"$x_1$\", fontsize=18)\n",
    "    plt.ylabel(r\"$x_2$\", fontsize=18, rotation=0)\n",
    "\n",
    "plt.figure(figsize=(11,4))\n",
    "plt.subplot(121)\n",
    "plot_decision_boundary(tree_clf, X, y)\n",
    "plt.title(\"Decision Tree\", fontsize=14)\n",
    "plt.subplot(122)\n",
    "plot_decision_boundary(bag_clf, X, y)\n",
    "plt.title(\"Decision Trees with Bagging\", fontsize=14)\n",
    "save_fig(\"decision_tree_without_and_with_bagging_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you need an additional performance measure, you can use the accuracy score:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "print(\"Bagging ensemble: %s\" % accuracy_score(y_test, y_pred))\n",
    "print(\"Single tree:      %s\" % accuracy_score(y_test, y_pred_tree))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Out-of-Bag evaluation\n",
    "\n",
    "When a bagging classifier is used, evaluation of the performance can be performed _out-of-bag_. Remember what bagging does, and how many instances (on average) are picked from all training instances if the bag is chosen to be the same size as the number of training instances. The fraction of chosen instances should converge against something like\n",
    "\n",
    "$$1 - \\exp(-1) \\approx 63.212\\%$$\n",
    "\n",
    "But that also means, that more than 35% of instances are _not seen_ in training. The [BaggingClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.BaggingClassifier.html) allows to set evaluation on out-of-bag instances automatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bag_clf = BaggingClassifier(\n",
    "    DecisionTreeClassifier(random_state=42), n_estimators=500,\n",
    "    bootstrap=True, n_jobs=-1, oob_score=True, random_state=40)\n",
    "bag_clf.fit(X_train, y_train)\n",
    "bag_clf.oob_score_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boosting via AdaBoost\n",
    "\n",
    "The performance of decision trees can be much improved through the procedure of _hypothesis boosting_. AdaBoost, probably the most popular algorithm, uses a very common technique: models are trained _sequentially_, where each model tries to correct for mistakes the previous model made. AdaBoost in particular _boosts_ the weights of those instances that were classified incorrectly. The next classifier will then be more sensitive to these instances and probably do an overall better job. In the end, the outputs of all sequential classifiers are combined into a prediction value. Each classifier enters this global value weighted according to its error rate. Please check/answer the following questions to familiarise yourself with AdaBoost:\n",
    "1. What is the error rate of a predictor? \n",
    "2. How is the weight for each predictor calculated? \n",
    "3. How are weights of instances updated if they were classified correctly? How are they updated if classified incorrectly?\n",
    "4. How is the final prediction made from an AdaBoost ensemble?\n",
    "5. The [AdaBoostClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.AdaBoostClassifier.html) implements the AdaBoost algorithm in Scikit-Learn. The following bit of code implements AdaBoost with decision tree classifiers. Make yourself familiar with the class and its arguments, then try to tweak it to achieve better performance than in the example below!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import AdaBoostClassifier\n",
    "\n",
    "ada_clf = AdaBoostClassifier(\n",
    "    DecisionTreeClassifier(max_depth=1), n_estimators=200,\n",
    "    algorithm=\"SAMME.R\", learning_rate=0.5, random_state=42)\n",
    "ada_clf.fit(X_train, y_train)\n",
    "\n",
    "plot_decision_boundary(ada_clf, X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following bit of code is a visualisation of how the weight adjustment in AdaBoost works. While not relying on the above AdaBoostClassifier class, this implements a support vector machine classifier ([SVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html)) and boosts the weights of incorrectly classified instances by hand. With different learning rates, the \"amount\" of boosting can be controlled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "m = len(X_train)\n",
    "\n",
    "plt.figure(figsize=(11, 4))\n",
    "for subplot, learning_rate in ((121, 1), (122, 0.5)):\n",
    "    # Start with equal weights for all instances.\n",
    "    sample_weights = np.ones(m)\n",
    "    plt.subplot(subplot)\n",
    "    \n",
    "    # Now let's go through five iterations where the\n",
    "    # weights get adjusted based on the previous step.\n",
    "    for i in range(5):\n",
    "        # As an example, use SVM classifier with Gaussian kernel.\n",
    "        svm_clf = SVC(kernel=\"rbf\", C=0.05, gamma=\"auto\", random_state=42)\n",
    "        svm_clf.fit(X_train, y_train, sample_weight=sample_weights)\n",
    "        y_pred = svm_clf.predict(X_train)\n",
    "        # The most important step: increase the weights of\n",
    "        # incorrectly predicted instances according to the\n",
    "        # learning_rate parameter.\n",
    "        sample_weights[y_pred != y_train] *= (1 + learning_rate)\n",
    "        # And do the plotting.\n",
    "        plot_decision_boundary(svm_clf, X, y, alpha=0.2)\n",
    "        plt.title(\"learning_rate = {}\".format(learning_rate), fontsize=16)\n",
    "    if subplot == 121:\n",
    "        plt.text(-0.7, -0.65, \"1\", fontsize=14)\n",
    "        plt.text(-0.6, -0.10, \"2\", fontsize=14)\n",
    "        plt.text(-0.5,  0.10, \"3\", fontsize=14)\n",
    "        plt.text(-0.4,  0.55, \"4\", fontsize=14)\n",
    "        plt.text(-0.3,  0.90, \"5\", fontsize=14)\n",
    "\n",
    "save_fig(\"boosting_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient Boosting\n",
    "\n",
    "An alternative to AdaBoost is gradient boosting. Again, gradient boosting sequentially trains multiple predictors which are then combined for a global prediction in the end. Gradient boosting fits the new predictor to the _residual errors_ made by the previous predictor, but doesn't touch instance weights. This can be visualised very well with a regression problem (of course, classification can also be performed. Scikit-Learn comes with the two classes [GradientBoostingRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html) and [GradientBoostingClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html) for these tasks. As a first step, the following example implements regression with decision trees by hand.\n",
    "\n",
    "First, generate our random data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "X = np.random.rand(100, 1) - 0.5\n",
    "y = 3*X[:, 0]**2 + 0.05 * np.random.randn(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "# Start with the first tree and fit it to X, y.\n",
    "tree_reg1 = DecisionTreeRegressor(max_depth=2, random_state=42)\n",
    "tree_reg1.fit(X, y)\n",
    "\n",
    "# Calculate the residual errors the previous tree\n",
    "# has made and fit a second tree to these.\n",
    "y2 = y - tree_reg1.predict(X)\n",
    "tree_reg2 = DecisionTreeRegressor(max_depth=2, random_state=42)\n",
    "tree_reg2.fit(X, y2)\n",
    "\n",
    "# Again, calculate the residual errors of the previous\n",
    "# tree and fit a third tree.\n",
    "y3 = y2 - tree_reg2.predict(X)\n",
    "tree_reg3 = DecisionTreeRegressor(max_depth=2, random_state=42)\n",
    "tree_reg3.fit(X, y3)\n",
    "\n",
    "# And the rest is just plotting ...\n",
    "def plot_predictions(regressors, X, y, axes, label=None, style=\"r-\", data_style=\"b.\", data_label=None):\n",
    "    x1 = np.linspace(axes[0], axes[1], 500)\n",
    "    y_pred = sum(regressor.predict(x1.reshape(-1, 1)) for regressor in regressors)\n",
    "    plt.plot(X[:, 0], y, data_style, label=data_label)\n",
    "    plt.plot(x1, y_pred, style, linewidth=2, label=label)\n",
    "    if label or data_label:\n",
    "        plt.legend(loc=\"upper center\", fontsize=16)\n",
    "    plt.axis(axes)\n",
    "\n",
    "plt.figure(figsize=(11,11))\n",
    "\n",
    "plt.subplot(321)\n",
    "plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label=\"$h_1(x_1)$\", style=\"g-\", data_label=\"Training set\")\n",
    "plt.ylabel(\"$y$\", fontsize=16, rotation=0)\n",
    "plt.title(\"Residuals and tree predictions\", fontsize=16)\n",
    "\n",
    "plt.subplot(322)\n",
    "plot_predictions([tree_reg1], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label=\"$h(x_1) = h_1(x_1)$\", data_label=\"Training set\")\n",
    "plt.ylabel(\"$y$\", fontsize=16, rotation=0)\n",
    "plt.title(\"Ensemble predictions\", fontsize=16)\n",
    "\n",
    "plt.subplot(323)\n",
    "plot_predictions([tree_reg2], X, y2, axes=[-0.5, 0.5, -0.5, 0.5], label=\"$h_2(x_1)$\", style=\"g-\", data_style=\"k+\", data_label=\"Residuals\")\n",
    "plt.ylabel(\"$y - h_1(x_1)$\", fontsize=16)\n",
    "\n",
    "plt.subplot(324)\n",
    "plot_predictions([tree_reg1, tree_reg2], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label=\"$h(x_1) = h_1(x_1) + h_2(x_1)$\")\n",
    "plt.ylabel(\"$y$\", fontsize=16, rotation=0)\n",
    "\n",
    "plt.subplot(325)\n",
    "plot_predictions([tree_reg3], X, y3, axes=[-0.5, 0.5, -0.5, 0.5], label=\"$h_3(x_1)$\", style=\"g-\", data_style=\"k+\")\n",
    "plt.ylabel(\"$y - h_1(x_1) - h_2(x_1)$\", fontsize=16)\n",
    "plt.xlabel(\"$x_1$\", fontsize=16)\n",
    "\n",
    "plt.subplot(326)\n",
    "plot_predictions([tree_reg1, tree_reg2, tree_reg3], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label=\"$h(x_1) = h_1(x_1) + h_2(x_1) + h_3(x_1)$\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=16)\n",
    "plt.ylabel(\"$y$\", fontsize=16, rotation=0)\n",
    "\n",
    "save_fig(\"gradient_boosting_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following piece of code now uses the Scikit-Learn class for regression with gradient boosting. Two examples are given: (1) with a fast learning rate, but only very few predictors, (2) with a slower learning rate, but a high number of predictors. Clearly, the second ensemble overfits the problem. Can you try to tweak the parameters to get a model that generalises better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "\n",
    "# First regression instance with only three estimaters,\n",
    "# but a fast learning rate. The max_depth parameter\n",
    "# controls the number of 'layers' in the decision\n",
    "# tree estimators of the ensemble. Increase for stronger\n",
    "# bias of the individual trees.\n",
    "gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=3, learning_rate=1.0, random_state=42)\n",
    "gbrt.fit(X, y)\n",
    "\n",
    "# Second instance with many estimators and slower\n",
    "# learning rate.\n",
    "gbrt_slow = GradientBoostingRegressor(max_depth=2, n_estimators=200, learning_rate=0.5, random_state=42)\n",
    "gbrt_slow.fit(X, y)\n",
    "\n",
    "plt.figure(figsize=(11,4))\n",
    "\n",
    "plt.subplot(121)\n",
    "plot_predictions([gbrt], X, y, axes=[-0.5, 0.5, -0.1, 0.8], label=\"Ensemble predictions\")\n",
    "plt.title(\"learning_rate={}, n_estimators={}\".format(gbrt.learning_rate, gbrt.n_estimators), fontsize=14)\n",
    "\n",
    "plt.subplot(122)\n",
    "plot_predictions([gbrt_slow], X, y, axes=[-0.5, 0.5, -0.1, 0.8])\n",
    "plt.title(\"learning_rate={}, n_estimators={}\".format(gbrt_slow.learning_rate, gbrt_slow.n_estimators), fontsize=14)\n",
    "\n",
    "save_fig(\"gbrt_learning_rate_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to solve this overfitting is to use _early stopping_ to find the optimal number of iterations/predictors for this problem. For that, we first need to split the dataset into training and validation set, because of course we cannot evaluate performance on instances the predictor used in training. The following code implements the already known [model_selection.train_test_split](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html) function. Then, train another ensemble (with a fixed number of 120 predictors), but this time only on the training set. Errors are calculated based on the validation set and the optimal number of iterations is extracted. The code also creates a plot of the performance on the validation set to point out the optimal iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "# Split dataset into training and validation set.\n",
    "X_train, X_val, y_train, y_val = train_test_split(X, y, random_state=49)\n",
    "\n",
    "# Fit an ensemble. Let's start with 120 estimators, which\n",
    "# is probably too much (as we saw above).\n",
    "gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=120, random_state=42)\n",
    "gbrt.fit(X_train, y_train)\n",
    "\n",
    "# Calculate the errors for each iteration (on the validation set)\n",
    "# and find the optimal iteration step.\n",
    "errors = [mean_squared_error(y_val, y_pred)\n",
    "          for y_pred in gbrt.staged_predict(X_val)]\n",
    "bst_n_estimators = np.argmin(errors)\n",
    "min_error = np.min(errors)\n",
    "\n",
    "# Retrain a new ensemble with those settings.\n",
    "gbrt_best = GradientBoostingRegressor(max_depth=2,n_estimators=bst_n_estimators, random_state=42)\n",
    "gbrt_best.fit(X_train, y_train)\n",
    "\n",
    "# And do the plotting of validation error as well\n",
    "# as the optimised ensemble.\n",
    "plt.figure(figsize=(11, 4))\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(errors, \"b.-\")\n",
    "plt.plot([bst_n_estimators, bst_n_estimators], [0, min_error], \"k--\")\n",
    "plt.plot([0, 120], [min_error, min_error], \"k--\")\n",
    "plt.plot(bst_n_estimators, min_error, \"ko\")\n",
    "plt.text(bst_n_estimators, min_error*1.2, \"Minimum\", ha=\"center\", fontsize=14)\n",
    "plt.axis([0, 120, 0, 0.01])\n",
    "plt.xlabel(\"Number of trees\")\n",
    "plt.title(\"Validation error\", fontsize=14)\n",
    "\n",
    "plt.subplot(122)\n",
    "plot_predictions([gbrt_best], X, y, axes=[-0.5, 0.5, -0.1, 0.8])\n",
    "plt.title(\"Best model (%d trees)\" % bst_n_estimators, fontsize=14)\n",
    "\n",
    "save_fig(\"early_stopping_gbrt_plot\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.3"
  },
  "nav_menu": {
   "height": "252px",
   "width": "333px"
  },
  "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": 1
}