03_decision_trees.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decision trees\n",
    "\n",
    "Welcome to the first jupyter notebook! In this session, we won't go into too many minute details yet, but will cover the basics of machine learning. We've tried to keep it as simple as possible, also because many might only have little experience with python and/or programming a 'learning machine'. If it looks like you're gonna be through the content of this notebook in ten minutes or so, because you're already familiar with all of its concepts, then feel free to challenge yourself a little more by trying out more advanced code snippets or doing some last-minute preparations for the group presentations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "To allow the next code blocks to run smoothly, this section sets a couple of settings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some imports that we will be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the random seed to a fixed number. This will guarantee that the notebook output is generated the same way for every run, otherwise the seed would be – random, as the name suggests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the numpy output precision for some prettier numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some figure plotting settings: increase the axis labels of our figures a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we start, we'll define the following function to show decision boundaries in our plots. Ignore the details of this implementation ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "def plot_decision_boundary(clf, X, y, axes=[-2.0, 3.0, -1.5, 2.0]):\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(['#FFC38D','#767676','#A7D1F2'])\n",
    "    plt.contourf(x1, x2, y_pred, alpha=0.3, cmap=custom_cmap)\n",
    "    plt.plot(X[:, 0][y==1], X[:, 1][y==1], \"s\")\n",
    "    plt.plot(X[:, 0][y==0], X[:, 1][y==0], \"o\")\n",
    "    plt.axis(axes)\n",
    "    plt.xlabel(r\"$x_1$\", fontsize=18)\n",
    "    plt.ylabel(r\"$x_2$\", fontsize=18, rotation=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with an easy classification problem and see what can be done about it with decision trees. The first step is to generate some random data. The Scikit Learn package provides a couple of functions for these, so we won't have to worry about it ourselves. The function used is [make_moons](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_moons.html). Run the following cell to import the function and generate 100 samples of two classes in half-moon shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons\n",
    "X, y = make_moons(n_samples=100, noise=0.25, random_state=53)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can examine the structure of our data as we did in the previous notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember: this command outputs the feature vectors of the first three instances. From what we can tell above, every of these instances has two features (the square brackets mark the beginning/the end of an entry in the overarching vector). Let's also have a look at the target values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look's like a typical classification problem: our entries are either 0's, meaning that this instance belongs to \"class 0\", or 1's for instances of \"class 1\". Since our dataset is reasonably small with only 100 instances, we should maybe just plot it to get an idea what we're talking about (ignore the details of these plotting commands):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X[:, 0][y==0], X[:, 1][y==0], \"o\")\n",
    "plt.plot(X[:, 0][y==1], X[:, 1][y==1], \"s\")\n",
    "plt.axis([-2.0, 3.0, -1.5, 2.0])\n",
    "plt.xlabel(r\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(r\"$x_2$\", fontsize=18, rotation=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, this indeed looks like two half moons, or something like a yin-and-yang structure. If we used something like a linear regression, we would definitely not be able to fit this dataset very well, but maybe a decision tree could help us with that?\n",
    "\n",
    "Let's define a decision tree! For that, we use the [DecisionTreeClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html) class of Scikit Learn, that comes with lots of functionalities useful for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeClassifier\n",
    "tree = DecisionTreeClassifier(random_state=42)\n",
    "tree.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the output, you can already see that the model comes with a bunch of parameters that we could have set by hand. Many of these have very resonable default values, such as `criterion='gini'`. Just a quick reminder of the Gini impurity index:\n",
    "\n",
    "$$ G_i = 1 - \\sum_{k=1}^n p_{i,k}^2 $$\n",
    "\n",
    "for classes $k=1,\\dots,n$ and leaf $i$ of a tree. We could also switch to using entropy as a measure when/where/how often to split the decision tree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we discuss all of these parameters in a bit more detail, let's just use the \"default\" decision-tree classifier and see how well it performs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_decision_boundary(tree, X, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm. The model seems to produce quite a few small \"boxes\" with only very few, sometimes even just one instance in it. Would you say this model looks particularly great? Would it generalise well on a new dataset generated from the same shape?\n",
    "\n",
    "We can also visualise the tree decisions with the graphviz package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import export_graphviz\n",
    "import graphviz\n",
    "dot_data = export_graphviz(tree, out_file=None,\n",
    "                           feature_names=['x_1', 'x_2'],\n",
    "                           filled=True, rounded=True, \n",
    "                           special_characters=True)\n",
    "graphviz.Source(dot_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe it would be useful to apply some sort of regularisation to the model. As we're extremely sensitive to outliers at the moment, maybe we could require a minimum number of instances for each leaf of the tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = DecisionTreeClassifier(min_samples_leaf=8, random_state=42)\n",
    "tree.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create the same plot as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_decision_boundary(tree, X, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualise the decisions once again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dot_data = export_graphviz(tree, out_file=None,\n",
    "                           feature_names=['x_1', 'x_2'],\n",
    "                           filled=True, rounded=True, \n",
    "                           special_characters=True)\n",
    "graphviz.Source(dot_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This definitely does look a lot better, doesn't it? There are a few outliers that the decision boundary doesn't catch, but probably – on a completely new dataset – this model would perform a lot better than the one above. We usually say: _it generalises well with unseen data_. There are many more regularisation parameters, feel free to try them out:\n",
    "* Maximum number of features considered: in our case we only have two features anyways, but in many cases this makes sense.\n",
    "* Maximum number of leaves in the entire tree.\n",
    "* Minimum decrease of impurity (in our case, the Gini impurity) when moving from one level to the next.\n",
    "* Minimum impurity value to still allow further splits.\n",
    "* Minimum number of instances required for a split.\n",
    "* ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bagging trees\n",
    "\n",
    "We've learnt about techniques to combine multiple decision trees into an _ensemble_ of models. Usually, this includes sampling instances from the total number of instances. Or in normal people's language: you make sure to not consider all, but just a _subset_ of datapoints when building each tree. This ensures a bit of variety in what the trees are sensitive to, what they focus on, etc. This might make them less powerful by themselves (we like to call them \"weak learners\"), but they combine into a \"strong learner\".\n",
    "\n",
    "As a quick reminder, there are two common techniques to sample our data: bagging (bootstrap aggregation), which is sampling _with_ replacement, and pasting, which is sampling _without_ replacement. Quick example:\n",
    "\n",
    "```\n",
    "Total dataset: \n",
    "[0, 1, 2, 3, 4, 5, 6, 7, 8]\n",
    "\n",
    "Bagging:\n",
    "Subset 1: [0, 4, 7]\n",
    "Subset 2: [2, 5, 7]\n",
    "Subset 3: [1, 1, 6]\n",
    "Not sampled: [3, 8]\n",
    "\n",
    "Pasting:\n",
    "Subset 1: [0, 4, 7]\n",
    "Subset 2: [2, 5, 8]\n",
    "Subset 3: [1, 3, 6]\n",
    "Not sampled: []\n",
    "```\n",
    "Whereas pasting makes sure that all instances are sampled and each instance only occurs once in all subsets, bagging allows duplicates of the same instance in the subsets. On the other hand, some instances are not sampled at all. Usually, bagging results in a greater variety among the subsets, but also a larger bias. When combining them, however, the ensemble usually performs quite well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we implement a bagging ensemble ourselves, we should make sure we can quantify how good a model actually is. For that, we should evaluate the model's performance on some test data it has never seen before. We've already generated the half-moon dataset above, but let's generate it again with a larger number of instances, and then use the Scikit Learn function [train_test_split](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html) to split it into two subsets: one for training, one for testing. With no parameters provided, the function splits the dataset into 3/4 and 1/4 for training and testing, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X, y = make_moons(n_samples=400, 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": [
    "Let's just test if it that splitting actually worked fine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks good: 300 instances of train data, 100 instances of test data, and each of them have two entries. Now, let's quickly train a simple decision tree on the training dataset again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_clf = DecisionTreeClassifier(random_state=42)\n",
    "tree_clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's predict the target values for our test dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = tree_clf.predict(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This we need to compare to the _true_ target values now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range (0, len(X_train)):\n",
    "    if i > 5: break\n",
    "    print(\"Predicted class: %s \\t True class: %s\" % (y_test[i], y_pred[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, we can do it like that, but of course there are automated ways to compare the two. Scikit Learn comes with a function to calculate the _accuracy score_, that is, the fraction of correctly predicted target values ([accuracy_score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html)). Let's try it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "print(accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so we're correct in about 90% of the cases. Let's compare this to a bagging ensemble. We will be using the [BaggingClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.BaggingClassifier.html) of Scikit Learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import BaggingClassifier\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, this comes with a bunch of parameters, most importantly: `n_estimators`, the number of models in the ensemble, and `max_samples`, the maximum number of instances that is fed into the individual models. `bootstrap=True` enables bootstrap aggregation, otherwise we would be using pasting. Let's calculate the predicted target values for the bagging ensemble, too, and then have a look at the accuracy score:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = bag_clf.predict(X_test)\n",
    "print(accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, looks better than before! To understand why, maybe we should just plot the two models overlaid with the training and the test data. The top row shows training data, the bottom row shows test data. What's your conclusion when looking at the plots? Which model probably has the higher accuracy score on the training data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(11,8))\n",
    "plt.subplot(221)\n",
    "plot_decision_boundary(tree_clf, X_train, y_train)\n",
    "plt.title(\"Decision Tree (train/test)\", fontsize=14)\n",
    "plt.subplot(222)\n",
    "plot_decision_boundary(bag_clf, X_train, y_train)\n",
    "plt.title(\"Bagging Ensemble (train/test)\", fontsize=14)\n",
    "plt.subplot(223)\n",
    "plot_decision_boundary(tree_clf, X_test, y_test)\n",
    "plt.subplot(224)\n",
    "plot_decision_boundary(bag_clf, X_test, y_test)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive boosting\n",
    "\n",
    "As a last step, let's look at one of the two common boosting techniques: AdaBoost, or adaptive boosting. The idea behind the boosting is to start off with a model that performs reasonably well, but then iteratively build improved models that \"correct\" the mistakes done before. In the particular case of AdaBoost, this includes reweighting instances that were misclassified in the previous step. This makes the model in the next iteration step a lot more sensitive to these \"difficult\" instances, and they are quite likely to be predicted better. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you're interested in the actual maths behind it, let me give a quick summary:\n",
    "1. Calculate an error rate for the current model $j$. $w^{(i)}$ is the weight of instance $i$:\n",
    "\n",
    "$$ r_j = \\frac{\\sum_{i=1}^{m} w^{(i)} ~ \\text{with $\\hat{y}^{(i)} \\neq y^{(i)}$}}{\\sum_{i=1}^{m} w^{(i)}}$$\n",
    "\n",
    "2. Calculate a weight for the current model $j$ based on its error rate $r_j$ ($\\eta$ is our usual learning rate):\n",
    "\n",
    "$$ \\alpha_j = \\eta \\log \\frac{1 - r_j}{r_j} $$\n",
    "\n",
    "3. Update the weights of the misclassified instances. If an instance was predicted correctly ($\\hat{y}^{(i)} = y^{(i)}$), do not make any changes to the weight:\n",
    "\n",
    "$$ w^{(i)} \\leftarrow w^{(i)} \\exp \\alpha_j \\quad \\text{if} ~\\hat{y}^{(i)} \\neq y^{(i)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, at some point you'll stop building up the ensemble (maybe just based on a maximum number of models you'd like to have). Then, the prediction of the ensemble is calculated according to:\n",
    "\n",
    "$$ \\hat{y} = \\text{argmax}_k \\sum_{j=1}^N \\alpha_j \\quad \\text{with} ~ \\hat{y}_j = k $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means you'll still consider _all_ models of your ensemble in the prediction, but their \"voting rights\" are weighted according to their model weight $\\alpha_j$, based on their error rates $r_j$. This ensures to \"trust\" those models more that generally perform better on the training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try this on our half-moon dataset. We'll be using the [AdaBoostClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.AdaBoostClassifier.html) of Scikit Learn."
   ]
  },
  {
   "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",
    "    learning_rate=0.2, random_state=42)\n",
    "ada_clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's make our predictions on the test data once again and calculate the accuracy score:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = ada_clf.predict(X_test)\n",
    "print(accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, this number is lower than that of the bagging ensemble. Do you have any idea why? Let's also have a look at the decision boundaries once again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(11,4))\n",
    "plt.subplot(121)\n",
    "plt.title(\"AdaBoost (train)\", fontsize=14)\n",
    "plot_decision_boundary(ada_clf, X_train, y_train)\n",
    "plt.subplot(122)\n",
    "plt.title(\"AdaBoost (test)\", fontsize=14)\n",
    "plot_decision_boundary(ada_clf, X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.14"
  },
  "nav_menu": {
   "height": "309px",
   "width": "468px"
  },
  "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
}