10 - Gradient Boosting/02 - XGBoost and ensemble models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/80x15.png\" /></a>&nbsp;| [Emmanuel Rachelson](https://personnel.isae-supaero.fr/emmanuel-rachelson?lang=en) | <a href=\"https://supaerodatascience.github.io/machine-learning/\">https://supaerodatascience.github.io/machine-learning/</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"font-size:22pt; line-height:25pt; font-weight:bold; text-align:center;\">XGBoost<br>Introduction to XGBoost</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Practice Course is composed of 3 parts - each part is meant to be done in about 1 hour :\n",
    "* In the **first notebook**, you learned the **basic of XGBoost**, how to apply it on a dataset and tune it to obtain the best performances.\n",
    "* In the **second notebook**, we will focus on **ensemble methods** and explain what makes XGBoost different from other models.\n",
    "* Finally in the **last notebook** you will see how the choice of a method (such as XGBoost) is a key element of a tradeoff between **Bias and Variance**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# About Ensemble methods\n",
    "\n",
    "Ensemble method are based on the hypothesis that combining multiple models together can often produce a much more powerful model.\n",
    "\n",
    "The goal of this notebook is to understand and manipulate the model behind XGBoost, to better understand the various parameters related to this model.\n",
    "We will have 3 exercices :\n",
    "* A short exercice focused on Weak Learners\n",
    "* Then we will see the difference between Boosting vs Bagging vs stacking\n",
    "* A specific focus on Ensemble Trees, where you will learn to tune the parameters of XGBoost that are specific to trees.\n",
    "\n",
    "# Weak learners\n",
    "\n",
    "**Weak learners** are basic models that do not perform so well by themselves, either because they have a high bias (low degree of freedom models, for example) or because they have too much variance to be robust (high degree of freedom models, for example). The idea of ensemble methods is to try reducing bias and/or variance of such weak learners by combining several of them together in order to create a strong learner (or ensemble model) that achieves better performances.\n",
    "\n",
    "In the following exemple, we will use a \"decision stumps.\" A decision stump is simply a decision tree where the whole tree is just one node. \n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Exercice 1</b><br>\n",
    "      Based on the following code, produce a weak learner on the given dataset.\n",
    "      <br>\n",
    "      Subsidiary question : what can be parallelized in the following code ?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split, cross_val_predict\n",
    "from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, accuracy_score\n",
    "from sklearn.tree import DecisionTreeClassifier \n",
    "\n",
    "class Stump:\n",
    "    def __init__(self):\n",
    "        self.gtLabel = None\n",
    "        self.ltLabel = None\n",
    "        self.splitThreshold = None\n",
    "        self.splitFeature = None\n",
    "\n",
    "    def predict(self, listpoint):\n",
    "        return np.array([self.predict_single(x) for x in listpoint])\n",
    "        \n",
    "    def predict_single(self, point):\n",
    "        if point[self.splitFeature] >= self.splitThreshold:\n",
    "            return self.gtLabel\n",
    "        else:\n",
    "            return self.ltLabel\n",
    "\n",
    "    def __call__(self, point):\n",
    "        return self.predict(point)\n",
    "\n",
    "    def errorFunction(self, data,index, t):\n",
    "        left_side = sum([1.0 for x,y in data if ... != y]) / len(data)\n",
    "        right_side = 1. - left_side\n",
    "        return min(left_side, right_side)\n",
    "\n",
    "\n",
    "    def makeThreshold(self,x, t,index):\n",
    "        return  1 if x[index] <= t else 0\n",
    "\n",
    "    def bestThreshold(self,data, index):\n",
    "        '''Compute best threshold for a given feature. Returns (threshold, error)'''\n",
    "\n",
    "        thresholds = [x[index] for x,y in data]\n",
    "        errors = [(threshold, self.errorFunction(data, index, threshold)) for threshold in thresholds]\n",
    "        return min(errors, key=lambda p: p[1])\n",
    "\n",
    "\n",
    "    def majorityVote(self,labels):\n",
    "        try:\n",
    "            return max(set(labels), key=labels.count)\n",
    "        except:\n",
    "            return 0\n",
    "\n",
    "    def fit(self,X, Y):\n",
    "        data = list(zip(X,Y))\n",
    "        # find the index of the best feature to split on, and the best threshold for that index.\n",
    "\n",
    "        bestThresholds = [(i,) + self.bestThreshold(data, i) for i in range(len(X[0]))]\n",
    "        feature, thresh, _ = min(bestThresholds, key = lambda p: p[2])\n",
    "\n",
    "        self.splitFeature = feature\n",
    "        self.splitThreshold = thresh\n",
    "        self.gtLabel = self.majorityVote([y for x,y in data if x[self.splitFeature] >= self.splitThreshold])\n",
    "        self.ltLabel = self.majorityVote([y for x,y in data if x[self.splitFeature] < self.splitThreshold])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "breast = datasets.load_breast_cancer()\n",
    "#print(breast.DESCR)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(breast.data, breast.target, test_size=0.2, random_state=42)\n",
    "\n",
    "stump = Stump()\n",
    "stump.fit(X_train,y_train)\n",
    "print(accuracy_score(y_train, stump.predict(X_train)))\n",
    "print(accuracy_score(y_test, stump.predict(X_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# you can use this code as a comparison point\n",
    "tree = DecisionTreeClassifier(max_depth=1, random_state=0)\n",
    "tree.fit(X_train,y_train)\n",
    "print(accuracy_score(y_test, tree.predict(X_test)))\n",
    "\n",
    "\n",
    "n_nodes = tree.tree_.node_count\n",
    "children_left = tree.tree_.children_left\n",
    "children_right = tree.tree_.children_right\n",
    "feature = tree.tree_.feature\n",
    "threshold = tree.tree_.threshold\n",
    "\n",
    "node_depth = np.zeros(shape=n_nodes, dtype=np.int64)\n",
    "is_leaves = np.zeros(shape=n_nodes, dtype=bool)\n",
    "stack = [(0, -1)]  # seed is the root node id and its parent depth\n",
    "while len(stack) > 0:\n",
    "    node_id, parent_depth = stack.pop()\n",
    "    node_depth[node_id] = parent_depth + 1\n",
    "\n",
    "    # If we have a test node\n",
    "    if (children_left[node_id] != children_right[node_id]):\n",
    "        stack.append((children_left[node_id], parent_depth + 1))\n",
    "        stack.append((children_right[node_id], parent_depth + 1))\n",
    "    else:\n",
    "        is_leaves[node_id] = True\n",
    "\n",
    "print(\"The binary tree structure has %s nodes and has \"\n",
    "      \"the following tree structure:\"\n",
    "      % n_nodes)\n",
    "for i in range(n_nodes):\n",
    "    if is_leaves[i]:\n",
    "        print(\"%snode=%s leaf node.\" % (node_depth[i] * \"\\t\", i))\n",
    "    else:\n",
    "        print(\"%snode=%s test node: go to node %s if X[:, %s] <= %s else to \"\n",
    "              \"node %s.\"\n",
    "              % (node_depth[i] * \"\\t\",\n",
    "                 i,\n",
    "                 children_left[i],\n",
    "                 feature[i],\n",
    "                 threshold[i],\n",
    "                 children_right[i],\n",
    "                 ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Boosting vs Bagging vs stacking\n",
    "The main hypothesis is that when weak models are correctly combined we can obtain more accurate and/or robust models.\n",
    "\n",
    "* **Bagging** considers homogeneous weak learners, learns them independently from each other in parallel and combines them following some kind of deterministic averaging process\n",
    "* **Boosting** considers homogeneous weak learners, learns them sequentially in a very adaptative way (a base model depends on the previous ones) and combines them following a deterministic strategy\n",
    "* **Stacking** considers heterogeneous weak learners, learns them in parallel and combines them by training a meta-model to output a prediction based on the different weak models predictions\n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Question 2</b><br>\n",
    "      Which kind is XGBoost ? Which kind is Random Forest ?\n",
    "</div>\n",
    "\n",
    "Each of these technics aims at improving either the *bias* or the *variance* of the individual weak learners.\n",
    "<img src=\"img/1 5pA6iY-qDP2JIsLoyfje-Q@2x.png\">\n",
    "\n",
    "## Focus on Boosting\n",
    "The intuition behind boosting is that each new model focus its efforts on the most difficult observations to fit up to now, so that we obtain, at the end of the process, a strong learner with lower bias. Boosting is a gerneral technic that is used with success by the the [top solutions](https://www.kaggle.com/c/porto-seguro-safe-driver-prediction/discussion/44629#250927) in Kaggle competition - it is a crucial competence to master by going beyond using black-box library.\n",
    "\n",
    "Being mainly focused at reducing bias, the base models that are often considered for boosting are models with low variance but high bias. For instance, if we want to use trees as our base models, we will choose most of the time shallow decision trees with only a few depths.\n",
    "\n",
    "There are two emblematic Boosting algorithms : **Adaboost** and **Gradient Boosting**. Since Gradient Boosting is the algorithm used in XGboost, we will try to understand how it works.\n",
    "\n",
    "Here are the steps of a simpler version of Gradient Boosting :\n",
    "* 1 - Fit a simple  decision tree on data\n",
    "* 2 - Calculate error residuals. Actual target value, minus predicted target value\n",
    "* 3 - Fit a new model on error residuals as target variable with same input variables\n",
    "* 4 - Add the predicted residuals to the previous predictions\n",
    "* 5 - Fit another model on residuals that is still left. Repeat steps 2 to 5 until overfitting or the sum of residuals becomes constant\n",
    "\n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Exercice 3</b><br>\n",
    "      Use the following code to improve the weak learner by one step of boosting. Consider it as a regression problem and use directly DecisionTreeRegressor as the stump. Evaluate the result with the mean_squared_error function\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "stump = DecisionTreeRegressor(max_depth = 1)\n",
    "stump.fit(X_train,y_train)\n",
    "prediction = stump.predict(X_train)\n",
    "print(1-mean_squared_error(y_train, prediction))\n",
    "print(1-mean_squared_error(y_test, stump.predict(X_test)))\n",
    "\n",
    "...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ensemble Tree\n",
    "XGBoost use Trees as the Weak Learners in the Gradient Boosting algorithm. A nice overview can be found on the [XGBoost documentation](https://xgboost.readthedocs.io/en/latest/tutorials/model.html).\n",
    "\n",
    "In XGBoost we can draw the constructed tree after fitting with the following code :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import style\n",
    "import seaborn as sns\n",
    "sns.set_style('whitegrid')\n",
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split, cross_val_predict\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.linear_model import SGDClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import LabelEncoder, MinMaxScaler\n",
    "from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score\n",
    "from sklearn.model_selection import cross_val_score\n",
    "import time\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning)\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "from xgboost import XGBClassifier\n",
    "import xgboost as xgb\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "X_train, X_test, y_train, y_test = train_test_split(iris.data, iris.target, test_size=0.2, random_state=42)\n",
    "\n",
    "dtrain = xgb.DMatrix(X_train, label=y_train)\n",
    "dtest = xgb.DMatrix(X_test, label=y_test)\n",
    "\n",
    "param = {'max_depth': 3,  'eta': 0.3, 'objective': 'multi:softprob', 'num_class': 3}\n",
    "num_round = 20 \n",
    "\n",
    "bst = xgb.train(param, dtrain, num_round)\n",
    "\n",
    "bst.dump_model('dump.raw.txt') # dump for model explanation\n",
    "\n",
    "with open('dump.raw.txt', 'r') as file:\n",
    "    data = file.read()\n",
    "    print (data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we have the *boosting* logic of having Decision Trees one after the other, while each single Decision Tree is very shallow.\n",
    "\n",
    "We will now focus on some parameters specific to trees. \n",
    "\n",
    "* n_estimators: number of trees you want to build.\n",
    "* max_depth: determines how deeply each tree is allowed to grow during any boosting round.\n",
    "* subsample: percentage of samples used per tree.\n",
    "* gamma: minimum loss reduction required to make a further partition on a leaf node of the tree.\n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Exercice 4</b><br>\n",
    "      Among these 4 parameters and using the following code, determine for each parameter if a low or a high value leads to Overfitting/Underfitting.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def ground_truth(x):\n",
    "    \"\"\"Ground truth -- function to approximate\"\"\"\n",
    "    return x * np.sin(x) + np.sin(2 * x)\n",
    "\n",
    "def gen_data(n_samples=200):\n",
    "    \"\"\"generate training and testing data\"\"\"\n",
    "    np.random.seed(13)\n",
    "    x = np.random.uniform(0, 10, size=n_samples)\n",
    "    x.sort()\n",
    "    y = ground_truth(x) + 0.75 * np.random.normal(size=n_samples)\n",
    "    train_mask = np.random.randint(0, 2, size=n_samples).astype(np.bool)\n",
    "    x_train, y_train = x[train_mask, np.newaxis], y[train_mask]\n",
    "    x_test, y_test = x[~train_mask, np.newaxis], y[~train_mask]\n",
    "    return x_train, x_test, y_train, y_test\n",
    "\n",
    "X_train, X_test, y_train, y_test = gen_data(200)\n",
    "\n",
    "# plot ground truth\n",
    "x_plot = np.linspace(0, 10, 500)\n",
    "def plot_data(figsize=(8, 5)):\n",
    "    fig = plt.figure(figsize=figsize)\n",
    "    gt = plt.plot(x_plot, ground_truth(x_plot), alpha=0.4, label='ground truth')\n",
    "\n",
    "    # plot training and testing data\n",
    "    plt.scatter(X_train, y_train, s=10, alpha=0.4)\n",
    "    plt.scatter(X_test, y_test, s=10, alpha=0.4, color='red')\n",
    "    plt.xlim((0, 10))\n",
    "    plt.ylabel('y')\n",
    "    plt.xlabel('x')\n",
    "    \n",
    "plot_data(figsize=(8, 5))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import islice\n",
    "from xgboost import XGBRegressor\n",
    "\n",
    "est = XGBRegressor(max_depth=3,subsample=1,gamma=0,n_estimators=100)\n",
    "est.fit(X_train, y_train)\n",
    "\n",
    "ax = plt.gca()\n",
    "\n",
    "# step over prediction as we added 20 more trees.\n",
    "plot_data()\n",
    "plt.plot(x_plot, est.predict(x_plot[:, np.newaxis]), color='r', alpha=0.8)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">**Takeaway Questions**<br>\n",
    "a) What is the model behind XGBoost?\n",
    "<br>\n",
    "b) What is the difference between a random forest and boosted trees?\n",
    "<br>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-danger\"><a href=\"#answers1\" data-toggle=\"collapse\">**Answers (click to unhide):**</a><br>\n",
    "<div id=\"answers1\" class=\"collapse\">\n",
    "a) Boosted trees, which are ensemble trees, i.e. a finite set of classification and regression trees. Like other ensemble methods, each tree participates in the prediction made by the model which is a sum of the outputs of all the trees.\n",
    "<br>\n",
    "b) Random forest and boosted trees use the same model but are trained in a different way. Random forests are trained using the technique of bootstrap aggregating, or bagging, while boosted trees are trained with gradient descent.\n",
    "<br>\n",
    "</div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion of the notebook\n",
    "\n",
    "In the notebook, we have seen the algorithm and the model behind XGBoost. \n",
    "The next and last notebook will continue the tuning of the parameters, with regards to the tradeoff between biais and variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sources :\n",
    "* https://medium.com/@aravanshad/ensemble-methods-95533944783f\n",
    "* https://www.datacamp.com/community/tutorials/xgboost-in-python\n",
    "* https://towardsdatascience.com/ensemble-methods-bagging-boosting-and-stacking-c9214a10a205\n",
    "* https://jeremykun.com/2015/05/18/boosting-census/\n",
    "* https://medium.com/mlreview/gradient-boosting-from-scratch-1e317ae4587d"
   ]
  }
 ],
 "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.8.10"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": false,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}