Xgboost 算法.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Xgboost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 在之前的算法中，梯度提升树都是使用一阶导求解来拟合残差，xgboost提出了一种使用二阶导，也就是海森矩阵来拟合\n",
    "\n",
    "1. [知乎资料](https://www.zhihu.com/question/41354392)\n",
    "2. [52cs](http://www.52cs.org/?p=429)\n",
    "3. [xgboost](http://www.kdd.org/kdd2016/papers/files/rfp0697-chenAemb.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn import datasets\n",
    "from sklearn.metrics import accuracy_score,mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 定义Sigmoid函数用来做对数几率回归\n",
    "\n",
    "class Sigmoid():\n",
    "    def __call__(self, x):\n",
    "        return 1 / (1 + np.exp(-x))\n",
    "\n",
    "    def gradient(self, x):\n",
    "        return self.__call__(x) * (1 - self.__call__(x))\n",
    "\n",
    "class LogisticLoss():\n",
    "    def __init__(self):\n",
    "        sigmoid = Sigmoid()\n",
    "        self.log_func = sigmoid\n",
    "        self.log_grad = sigmoid.gradient\n",
    "\n",
    "    def loss(self, y, y_pred):\n",
    "        y_pred = np.clip(y_pred, 1e-15, 1 - 1e-15)\n",
    "        p = self.log_func(y_pred)\n",
    "        return y * np.log(p) + (1 - y) * np.log(1 - p)\n",
    "\n",
    "    # gradient w.r.t y_pred\n",
    "    def gradient(self, y, y_pred):\n",
    "        p = self.log_func(y_pred)\n",
    "        return -(y - p)\n",
    "\n",
    "    # w.r.t y_pred\n",
    "    def hess(self, y, y_pred):\n",
    "        p = self.log_func(y_pred)\n",
    "        return p * (1 - p)\n",
    "def to_categorical(x, n_col=None):\n",
    "    \"\"\" One-hot encoding of nominal values \"\"\"\n",
    "    if not n_col:\n",
    "        n_col = np.amax(x) + 1\n",
    "    one_hot = np.zeros((x.shape[0], n_col))\n",
    "    one_hot[np.arange(x.shape[0]), x] = 1\n",
    "    return one_hot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def divide_on_feature(X, feature_i, threshold):\n",
    "    \"\"\" Divide dataset based on if sample value on feature index is larger than\n",
    "        the given threshold \"\"\"\n",
    "    split_func = None\n",
    "    if isinstance(threshold, int) or isinstance(threshold, float):\n",
    "        split_func = lambda sample: sample[feature_i] >= threshold\n",
    "    else:\n",
    "        split_func = lambda sample: sample[feature_i] == threshold\n",
    "\n",
    "    X_1 = np.array([sample for sample in X if split_func(sample)])\n",
    "    X_2 = np.array([sample for sample in X if not split_func(sample)])\n",
    "\n",
    "    return np.array([X_1, X_2])\n",
    "class DecisionNode():\n",
    "    \"\"\"Class that represents a decision node or leaf in the decision tree\n",
    "    Parameters:\n",
    "    -----------\n",
    "    feature_i: int\n",
    "        Feature index which we want to use as the threshold measure.\n",
    "    threshold: float\n",
    "        The value that we will compare feature values at feature_i against to \n",
    "        determine the prediction.\n",
    "    value: float\n",
    "        The class prediction if classification tree, or float value if regression tree.\n",
    "    true_branch: DecisionNode\n",
    "        Next decision node for samples where features value met the threshold.\n",
    "    false_branch: DecisionNode\n",
    "        Next decision node for samples where features value did not meet the threshold.\n",
    "    \"\"\"\n",
    "    def __init__(self, feature_i=None, threshold=None,\n",
    "                 value=None, true_branch=None, false_branch=None):\n",
    "        self.feature_i = feature_i          # Index for the feature that is tested\n",
    "        self.threshold = threshold          # Threshold value for feature\n",
    "        self.value = value                  # Value if the node is a leaf in the tree\n",
    "        self.true_branch = true_branch      # 'Left' subtree\n",
    "        self.false_branch = false_branch    # 'Right' subtree\n",
    "\n",
    "\n",
    "# Super class of RegressionTree and ClassificationTree\n",
    "class DecisionTree(object):\n",
    "    \"\"\"Super class of RegressionTree and ClassificationTree.\n",
    "    Parameters:\n",
    "    -----------\n",
    "    min_samples_split: int\n",
    "        The minimum number of samples needed to make a split when building a tree.\n",
    "    min_impurity: float\n",
    "        The minimum impurity required to split the tree further. \n",
    "    max_depth: int\n",
    "        The maximum depth of a tree.\n",
    "    loss: function\n",
    "        Loss function that is used for Gradient Boosting models to calculate impurity.\n",
    "    \"\"\"\n",
    "    def __init__(self, min_samples_split=2, min_impurity=1e-7,\n",
    "                 max_depth=float(\"inf\"), loss=None):\n",
    "        self.root = None  # Root node in dec. tree\n",
    "        # Minimum n of samples to justify split\n",
    "        self.min_samples_split = min_samples_split\n",
    "        # The minimum impurity to justify split\n",
    "        self.min_impurity = min_impurity\n",
    "        # The maximum depth to grow the tree to\n",
    "        self.max_depth = max_depth\n",
    "        # Function to calculate impurity (classif.=>info gain, regr=>variance reduct.)\n",
    "        self._impurity_calculation = None\n",
    "        # Function to determine prediction of y at leaf\n",
    "        self._leaf_value_calculation = None\n",
    "        # If y is one-hot encoded (multi-dim) or not (one-dim)\n",
    "        self.one_dim = None\n",
    "        # If Gradient Boost\n",
    "        self.loss = loss\n",
    "\n",
    "    def fit(self, X, y, loss=None):\n",
    "        \"\"\" Build decision tree \"\"\"\n",
    "        self.one_dim = len(np.shape(y)) == 1\n",
    "        self.root = self._build_tree(X, y)\n",
    "        self.loss=None\n",
    "\n",
    "    def _build_tree(self, X, y, current_depth=0):\n",
    "        \"\"\" Recursive method which builds out the decision tree and splits X and respective y\n",
    "        on the feature of X which (based on impurity) best separates the data\"\"\"\n",
    "\n",
    "        largest_impurity = 0\n",
    "        best_criteria = None    # Feature index and threshold\n",
    "        best_sets = None        # Subsets of the data\n",
    "\n",
    "        # Check if expansion of y is needed\n",
    "        if len(np.shape(y)) == 1:\n",
    "            y = np.expand_dims(y, axis=1)\n",
    "\n",
    "        # Add y as last column of X\n",
    "        Xy = np.concatenate((X, y), axis=1)\n",
    "\n",
    "        n_samples, n_features = np.shape(X)\n",
    "\n",
    "        if n_samples >= self.min_samples_split and current_depth <= self.max_depth:\n",
    "            # Calculate the impurity for each feature\n",
    "            for feature_i in range(n_features):\n",
    "                # All values of feature_i\n",
    "                feature_values = np.expand_dims(X[:, feature_i], axis=1)\n",
    "                unique_values = np.unique(feature_values)\n",
    "\n",
    "                # Iterate through all unique values of feature column i and\n",
    "                # calculate the impurity\n",
    "                for threshold in unique_values:\n",
    "                    # Divide X and y depending on if the feature value of X at index feature_i\n",
    "                    # meets the threshold\n",
    "                    Xy1, Xy2 = divide_on_feature(Xy, feature_i, threshold)\n",
    "                    \n",
    "                    if len(Xy1) > 0 and len(Xy2) > 0:\n",
    "                        # Select the y-values of the two sets\n",
    "                        y1 = Xy1[:, n_features:]\n",
    "                        y2 = Xy2[:, n_features:]\n",
    "\n",
    "                        # Calculate impurity\n",
    "                        impurity = self._impurity_calculation(y, y1, y2)\n",
    "\n",
    "                        # If this threshold resulted in a higher information gain than previously\n",
    "                        # recorded save the threshold value and the feature\n",
    "                        # index\n",
    "                        if impurity > largest_impurity:\n",
    "                            largest_impurity = impurity\n",
    "                            best_criteria = {\"feature_i\": feature_i, \"threshold\": threshold}\n",
    "                            best_sets = {\n",
    "                                \"leftX\": Xy1[:, :n_features],   # X of left subtree\n",
    "                                \"lefty\": Xy1[:, n_features:],   # y of left subtree\n",
    "                                \"rightX\": Xy2[:, :n_features],  # X of right subtree\n",
    "                                \"righty\": Xy2[:, n_features:]   # y of right subtree\n",
    "                                }\n",
    "\n",
    "        if largest_impurity > self.min_impurity:\n",
    "            # Build subtrees for the right and left branches\n",
    "            true_branch = self._build_tree(best_sets[\"leftX\"], best_sets[\"lefty\"], current_depth + 1)\n",
    "            false_branch = self._build_tree(best_sets[\"rightX\"], best_sets[\"righty\"], current_depth + 1)\n",
    "            return DecisionNode(feature_i=best_criteria[\"feature_i\"], threshold=best_criteria[\n",
    "                                \"threshold\"], true_branch=true_branch, false_branch=false_branch)\n",
    "\n",
    "        # We're at leaf => determine value\n",
    "        leaf_value = self._leaf_value_calculation(y)\n",
    "\n",
    "        return DecisionNode(value=leaf_value)\n",
    "\n",
    "\n",
    "    def predict_value(self, x, tree=None):\n",
    "        \"\"\" Do a recursive search down the tree and make a prediction of the data sample by the\n",
    "            value of the leaf that we end up at \"\"\"\n",
    "\n",
    "        if tree is None:\n",
    "            tree = self.root\n",
    "\n",
    "        # If we have a value (i.e we're at a leaf) => return value as the prediction\n",
    "        if tree.value is not None:\n",
    "            return tree.value\n",
    "\n",
    "        # Choose the feature that we will test\n",
    "        feature_value = x[tree.feature_i]\n",
    "\n",
    "        # Determine if we will follow left or right branch\n",
    "        branch = tree.false_branch\n",
    "        if isinstance(feature_value, int) or isinstance(feature_value, float):\n",
    "            if feature_value >= tree.threshold:\n",
    "                branch = tree.true_branch\n",
    "        elif feature_value == tree.threshold:\n",
    "            branch = tree.true_branch\n",
    "\n",
    "        # Test subtree\n",
    "        return self.predict_value(x, branch)\n",
    "\n",
    "    def predict(self, X):\n",
    "        \"\"\" Classify samples one by one and return the set of labels \"\"\"\n",
    "        y_pred = []\n",
    "        for x in X:\n",
    "            y_pred.append(self.predict_value(x))\n",
    "        return y_pred\n",
    "\n",
    "    def print_tree(self, tree=None, indent=\" \"):\n",
    "        \"\"\" Recursively print the decision tree \"\"\"\n",
    "        if not tree:\n",
    "            tree = self.root\n",
    "\n",
    "        # If we're at leaf => print the label\n",
    "        if tree.value is not None:\n",
    "            print (tree.value)\n",
    "        # Go deeper down the tree\n",
    "        else:\n",
    "            # Print test\n",
    "            print (\"%s:%s? \" % (tree.feature_i, tree.threshold))\n",
    "            # Print the true scenario\n",
    "            print (\"%sT->\" % (indent), end=\"\")\n",
    "            self.print_tree(tree.true_branch, indent + indent)\n",
    "            # Print the false scenario\n",
    "            print (\"%sF->\" % (indent), end=\"\")\n",
    "            self.print_tree(tree.false_branch, indent + indent)\n",
    "\n",
    "\n",
    "\n",
    "class XGBoostRegressionTree(DecisionTree):\n",
    "    \"\"\"\n",
    "    Regression tree for XGBoost\n",
    "    - Reference -\n",
    "    http://xgboost.readthedocs.io/en/latest/model.html\n",
    "    \"\"\"\n",
    "\n",
    "    def _split(self, y):\n",
    "        \"\"\" y contains y_true in left half of the middle column and \n",
    "        y_pred in the right half. Split and return the two matrices \"\"\"\n",
    "        col = int(np.shape(y)[1]/2)\n",
    "        y, y_pred = y[:, :col], y[:, col:]\n",
    "        return y, y_pred\n",
    "\n",
    "    def _gain(self, y, y_pred):\n",
    "        nominator = np.power((y * self.loss.gradient(y, y_pred)).sum(), 2)\n",
    "        denominator = self.loss.hess(y, y_pred).sum()\n",
    "        return 0.5 * (nominator / denominator)\n",
    "\n",
    "    def _gain_by_taylor(self, y, y1, y2):\n",
    "        # Split\n",
    "        y, y_pred = self._split(y)\n",
    "        y1, y1_pred = self._split(y1)\n",
    "        y2, y2_pred = self._split(y2)\n",
    "\n",
    "        true_gain = self._gain(y1, y1_pred)\n",
    "        false_gain = self._gain(y2, y2_pred)\n",
    "        gain = self._gain(y, y_pred)\n",
    "        return true_gain + false_gain - gain\n",
    "\n",
    "    def _approximate_update(self, y):\n",
    "        # y split into y, y_pred\n",
    "        y, y_pred = self._split(y)\n",
    "        # Newton's Method\n",
    "        gradient = np.sum(y * self.loss.gradient(y, y_pred), axis=0)\n",
    "        hessian = np.sum(self.loss.hess(y, y_pred), axis=0)\n",
    "        update_approximation =  gradient / hessian \n",
    "\n",
    "        return update_approximation\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        self._impurity_calculation = self._gain_by_taylor\n",
    "        self._leaf_value_calculation = self._approximate_update\n",
    "        super(XGBoostRegressionTree, self).fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class XGBoost(object):\n",
    "    \"\"\"The XGBoost classifier.\n",
    "    Reference: http://xgboost.readthedocs.io/en/latest/model.html\n",
    "    Parameters:\n",
    "    -----------\n",
    "    n_estimators: int\n",
    "        The number of classification trees that are used.\n",
    "    learning_rate: float\n",
    "        The step length that will be taken when following the negative gradient during\n",
    "        training.\n",
    "    min_samples_split: int\n",
    "        The minimum number of samples needed to make a split when building a tree.\n",
    "    min_impurity: float\n",
    "        The minimum impurity required to split the tree further. \n",
    "    max_depth: int\n",
    "        The maximum depth of a tree.\n",
    "    \"\"\"\n",
    "    def __init__(self, n_estimators=200, learning_rate=0.001, min_samples_split=2,\n",
    "                 min_impurity=1e-7, max_depth=2):\n",
    "        self.n_estimators = n_estimators            # Number of trees\n",
    "        self.learning_rate = learning_rate          # Step size for weight update\n",
    "        self.min_samples_split = min_samples_split  # The minimum n of sampels to justify split\n",
    "        self.min_impurity = min_impurity              # Minimum variance reduction to continue\n",
    "        self.max_depth = max_depth                  # Maximum depth for tree\n",
    "\n",
    "        \n",
    "        # Log loss for classification\n",
    "        self.loss = LogisticLoss()\n",
    "\n",
    "        # Initialize regression trees\n",
    "        self.trees = []\n",
    "        for _ in range(n_estimators):\n",
    "            tree = XGBoostRegressionTree(\n",
    "                    min_samples_split=self.min_samples_split,\n",
    "                    min_impurity=min_impurity,\n",
    "                    max_depth=self.max_depth,\n",
    "                    loss=self.loss)\n",
    "\n",
    "            self.trees.append(tree)\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        y = to_categorical(y)\n",
    "\n",
    "        y_pred = np.zeros(np.shape(y))\n",
    "        for i in range(self.n_estimators):\n",
    "            tree = self.trees[i]\n",
    "            y_and_pred = np.concatenate((y, y_pred), axis=1)\n",
    "            tree.fit(X, y_and_pred)\n",
    "            update_pred = tree.predict(X)\n",
    "\n",
    "            y_pred -= np.multiply(self.learning_rate, update_pred)\n",
    "\n",
    "    def predict(self, X):\n",
    "        y_pred = None\n",
    "        # Make predictions\n",
    "        for tree in self.trees:\n",
    "            # Estimate gradient and update prediction\n",
    "            update_pred = tree.predict(X)\n",
    "            if y_pred is None:\n",
    "                y_pred = np.zeros_like(update_pred)\n",
    "            y_pred -= np.multiply(self.learning_rate, update_pred)\n",
    "\n",
    "        # Turn into probability distribution (Softmax)\n",
    "        y_pred = np.exp(y_pred) / np.sum(np.exp(y_pred), axis=1, keepdims=True)\n",
    "        # Set label to the value that maximizes probability\n",
    "        y_pred = np.argmax(y_pred, axis=1)\n",
    "        return y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = datasets.load_iris()\n",
    "X = data.data\n",
    "y = data.target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = XGBoost()\n",
    "clf.fit(X_train, y_train)\n",
    "y_pred = clf.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy = accuracy_score(y_test,y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.94666666666666666"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}