XGBoost详解

Prerequisite: CART回归树

CART回归树是假设树为二叉树，通过不断将特征进行分裂。比如当前树结点是基于第j个特征值进行分裂的，设该特征值小于s的样本划分为左子树，大于s的样本划分为右子树。
$R _ { 1 } ( j , s ) = \{ x | x ^ { ( j ) } \leq s \} \text { and } R _ { 2 } ( j , s ) = \{ x | x ^ { ( j ) } > s \}$
而CART回归树实质上就是在该特征维度对样本空间进行划分, 而这种空间划分的优化是一种NP难问题，因此在决策树模型中是使用启发式方法解决。典型CART回归树产生的目标函数为:
$\sum _ { x _ { i } \in R _ { m } } \left( y _ { i } - f \left( x _ { i } \right) \right) ^ { 2 }$
当我们为了求解最优的切分特征j和最优的切分点s，就转化为求解这么一个目标函数:
$\min _ { j , s } \left[ \min _ { c _ { 1 } } \sum _ { x _ { i } \in R _ { 1 } ( j , s ) } \left( y _ { i } - c _ { 1 } \right) ^ { 2 } + \min _ { c _ { 2 } } \sum _ { x _ { i } \in R _ { 2 } ( j , s ) } \left( y _ { i } - c _ { 2 } \right) ^ { 2 } \right]$

XGBoost数学原理推导

该算法思想就是不断地添加树，不断地进行特征分裂来生长一棵树，每次添加一个树，其实是学习一个新函数，去拟合上次预测的残差。当我们训练完成得到k棵树，我们要预测一个样本的分数，其实就是根据这个样本的特征，在每棵树中会落到对应的一个叶子节点，每个叶子节点就对应一个分数，最后只需要将每棵树对应的分数加起来就是该样本的预测值。

如下图例子，训练出了2棵决策树，小孩的预测分数就是两棵树中小孩所落到的结点的分数相加。爷爷的预测分数同理。

爷爷的预测分数

XGBoost考虑正则化项，目标函数定义如下:
$L ( \phi ) = \sum _ { i } l \left( y _ { i } , \hat { y } _ { i } \right) + \sum _ { k } \Omega \left( f _ { k } \right) ,\quad 其中 \quad \Omega \left( f _ { k } \right) = \gamma T + \frac { 1 } { 2 } \lambda | | w | | ^ { 2 }$
很显然, $\sum _ { i } l \left( y _ { i } , \hat { y } _ { i } \right)$ 代表损失函数, 而 $\sum _ { k } \Omega \left( f _ { k } \right)$ 代表正则化项

其中 $\hat{y}_i$ 为预测输出， $y_i$ 为label值， $f_k$ 为第k树模型， $T$ 为树叶子节点数， $w$ 为叶子权重值， $\gamma$ 为叶子树惩罚正则项，具有剪枝作用， $\lambda$ 为叶子权重惩罚正则项，防止过拟合。XGBoost也支持一阶正则化，容易优化叶子节点权重为0，不过不常用。
$\hat { y } = \phi \left( x _ { i } \right) = \sum _ { k = 1 } ^ { K } f _ { k } \left( x _ { i } \right)$
$\text {where } F = \left\{ f ( x ) = w _ { q ( x ) } \right\} \left( q : R ^ { m } \rightarrow T , w \in R ^ { T }\right)$
$w_{q(x)}$ 是叶子节点q的分数, $q(x)$ 是叶子节点的编号, $f(x)$ 是其中一颗回归树。也就是说对于任意一个样本 $x$ , 其最后会落在树的某个叶子节点上, 其值为 $w_{q(x)}$

正如上文所说，新生成的树是要拟合上次预测的残差的，即当生成t棵树后，预测分数可以写成:
$\hat{y}_i^{(t)}=\hat{y}_i^{(t-1)}+f_t(x)$
同时，可以将目标函数改写成(其中 $f_t(x_i)$ 代表当前树上某个叶子节点上的值):
$L^{(t)}=\sum_{i=1} l\ (\ y_i\ ,\ (\ \hat{y}_i^{(t-1)} + f_t(x_i)\ ) \ )+ \Omega(f_t)$
很明显，我们接下来就是要去找到一个 $f_t$ 能够最小化目标函数。XGBoost的想法是利用其在 $f_t=0$ 处的 $\color{blue}{泰勒二阶}$ 展开近似它。所以，目标函数近似为：
$\mathcal { L } ^ { ( t ) } \simeq \sum _ { i = 1 } ^ { n } \left[ l \left( y _ { i } , \hat { y } ^ { ( t - 1 ) } \right) + g _ { i } f _ { t } \left( \mathbf { x } _ { i } \right) + \frac { 1 } { 2 } h _ { i } f _ { t } ^ { 2 } \left( \mathbf { x } _ { i } \right) \right] + \Omega \left( f _ { t } \right)$
其中 $g_i$ 为一阶导数， $h_i$ 为二阶导数, 其中 $h_i$ 前面的 $\frac{1}{2}$ 为泰勒二阶展开后的二阶导的系数 $\frac{1}{2!}$ :
$g_i=\partial_{\hat{y}_i^{(t-1)}}\;l(y_i, \hat{y}_i^{(t-1)}) \quad , \quad h_i=\partial_{\hat{y}_i^{(t-1)}}^2 \; l(y_i, \hat{y}_i^{(t-1)})$
为什么此处可以用二阶导呢?

首先我们需要明确一个概念, 我们的boosting每一轮迭代是在优化什么呢? 换句话说我们在用损失函数 $\mathcal { l }$ 在干什么呢? 其实我们看Boosting的框架, 无论是GBDT还是Adaboost, 其在每一轮迭代中, 根本没有理会损失函数具体是什么, 仅仅用到了损失函数的一阶导数。仅仅用一阶导数的问题就是, 我们根本无法保证我们找到的是全局最优。除非问题本身 $f(x)$ 是强凸的 (convex)而且最好是smooth的。每次迭代相当于就是最优化负梯度。即下式中的 $\nabla f(x^{(i)})$ , 因为是负梯度所以前面要加负号, $c$ 代表学习率。
$x_c^{(i+1)} = x^{(i)} - c \nabla f(x^{(i)})$
有没有感觉这个公式形式很熟悉, 是不是就是一般常见的linear regression 的 stochastic 梯度更新公式。既然GBDT用的是Stochastic Gradient Descent, 我们回想一下, 是不是还有别的梯度更新的方式, 这时, 牛顿法 Newton's method就出现了。可以说, XGBoost是用了牛顿法进行的梯度更新。

我们先来看一下泰勒展开式:
$\sum_{n=0} ^\infty \frac{1}{n!} \frac{\partial^n f(x)}{\partial x^n}\Delta x^n=f(x) + \frac{\partial f(x)}{\partial x}\Delta x + \frac{1}{2!}\frac{\partial^2 f(x)}{\partial x^2}\Delta x^2 + \cdots$
而对于上面这个公式值, 其与 $f(x+\Delta x)$ 之间的误差, 可以用如下公式表示:
$f(x+\Delta x)- \sum_{n=0} ^\infty \frac{1}{n!} \frac{\partial^n f(x)}{\partial x^n}\Delta x^n=h_k(x) a^k$
我们保留二阶泰勒展开:
$\begin{split} f(x+\Delta x)\simeq f(x)+f'(x)\Delta x+\frac12f''(x)\Delta x^2 \end{split}$
这个式子是成立的, 当且仅当 $\Delta x$ 趋近于0, 我们对上式求导 (对 $\Delta x$ 求导) 并令导数为0。
$f'(x)+f''(x)\Delta x=0$
即得到:
$\Delta x=-\frac{f'(x)}{f''(x)}$
所以得出迭代公式:
$x_{n+1}=x_{n}-\frac{f'(x)}{f''(x)}$
将损失函数与 $f(x)$ 对应起来:
$\color{blue}{\begin{split} f(x)&=l(y_i, \hat{y}_i^{(t-1)} )&=(y_i - \hat{y}_i^{(t-1)})^2 \\ f(x+\Delta x)&=l(y_i, \hat{y}_i^{(t-1)} + f_t(x_i))& =(y_i - (\ \hat{y}_i^{(t-1)} + f_t(x_i)\ )\ )^2 \end{split}}$
所以实际上 x 即为 $\hat{y}^{t-1}_i$ ，而 $\Delta x$ 即为 $f_t(x_i)$ 。 $\color{red}{故对f(x)求导数时即对\hat{y}^{t-1}_i求偏导}$ ，故根据二阶泰勒展开, $f(x+\Delta x)$ 可以表示为:
$f ( x + \Delta x )=\left( y _ { i } - \hat { y } _ { i } ^ { ( t - 1 ) } \right) ^ { 2 }+\left[\left( y _ { i } - \hat { y } _ { i } ^ { ( t - 1 ) } \right) ^ { 2 }\right]'*f_t(x_i)+\frac{1}{2}\left[\left( y _ { i } - \hat { y } _ { i } ^ { ( t - 1 ) } \right) ^ { 2 }\right]''*f_t^2(x_i)$
用 $g_i\ 和\ h_i$ 替代上式中的 $f'(x)\ 和\ f''(x)$ , 即可得到 (此处的 $obj^{(t)}$ 就是上面的 $\mathcal { L } ^ { ( t ) }$ , 不同写法而已):
$\begin{split} \text{obj}^{(t)} = \sum_{i=1}^n [l(y_i, \hat{y}_i^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i)] + \Omega(f_t) + constant \end{split}$
这里有必要再明确一下， $g_i$ 和 $h_i$ 的含义。gi怎么理解呢？假设现有 $t-1$ 棵树, 这 $t-1$ 棵树组成的模型对第i个训练样本有一个预测值 $\hat{y}_i$ 。这个 $\hat{y}_i$ 与第i个样本的真实标签 $y_i$ 肯定有差距, 这个差距可以用 $l(yi,\hat{y}_i)$ 这个损失函数来衡量。所以此处的 $g_i$ 和 $h_i$ 就是对于该损失函数的一阶导和二阶导。

我们来看一个具体的例子，假设我们正在优化第11棵CART树，也就是说前10棵 CART树已经确定了。这10棵树对样本 $(x_i,y_i=1)$ 的预测值是 $\hat{y}_i=-1$ ，假设我们现在是做分类，我们的损失函数是:
$\large L ( \theta ) = \sum _ { i } \left[ y _ { i } \ln \left( 1 + e ^ { - \hat { y } _ { i } } \right) + \left( 1 - y _ { i } \right) \ln \left( 1 + e ^ { \hat { y } _ { i } } \right) \right]$
在类别标签 $在y_i=1$ 时, $(1-y_i)=0$ 故此时损失函数变成了:
$\large \ln \left( 1 + e ^ { - \hat { y } _ { i } } \right)$
我们可以求出这个损失函数对于 $\hat{y}_i$ 的梯度，如下所示:
$\large (-1)*\frac{1}{1+e^{-\hat{y}_i}}*e^{-\hat{y}_i}=\frac{-e^{-\hat{y}_i}}{1+e^{-\hat{y}_i}}$
将 $\hat{y}_i =-1$ 代入上面的式子，计算得到-0.27。这个-0.27就是 $g_i$ 。该值是负的，也就是说，如果我们想要减小这10棵树在该样本点上的预测损失，我们应该沿着梯度的反方向去走，也就是要增大 $\hat{y}_i =-1$ 的值,使其趋向于正，因为我们的 $y_i=1$ 就是正的。

那么在优化第t棵树时，有多少个 $g_i$ 和 $h_i$ 要计算？嗯，没错就是各有N个，N是训练样本的数量。如果有10万样本，在优化第t棵树时，就需要计算出个10万个 $g_i$ 和 $h_i$ 。很显然, 这10万个 $g_i$ 之间并没有什么关系, 那么是不是可以并行计算呢？这就是XGBoost速度如此快的原因。

再总结一下, $g_i$ 和 $h_i$ 是可以并行地求出来的。而且，在并行计算的时候 $g_i$ 和 $h_i$ 是不依赖于损失函数的形式的，只要这个损失函数二次可微就可以了。这有什么好处呢？好处就是 $\color{red}{xgboost可以支持自定义损失函数，只需满足二次可微即可。}$

不过此处所说的不依赖于损失函数形式, 意思不是说在自定义损失函数的时候, 给xgboost传入一个损失函数, xgboost就会自动计算其一阶和二阶导。相反的, $\color{red}{我们需要给xgboost传入一个损失函数, 同时还需要给xgboost传入这个损失函数的一阶导和二阶导的形式。}$ 一阶二阶导和损失函数之间其实算是独立的。XGBoost在每次并行计算的时候独立调用损失函数和一阶二阶导。

我们来看一下python接口自定义损失函数调用XGBoost的方式(代码地址https://github.com/dmlc/xgboost/blob/master/demo/guide-python/custom_objective.py):

#coding:utf-8
import numpy as np
import xgboost as xgb

print('start running example to used customized objective function')

dtrain = xgb.DMatrix('../data/agaricus.txt.train')
dtest = xgb.DMatrix('../data/agaricus.txt.test')

param = {'max_depth': 2, 'eta': 1, 'silent': 1}
watchlist = [(dtest, 'eval'), (dtrain, 'train')]
num_round = 2

#下面这部分就是自定义损失函数, 可以看到, 函数接收两个向量, 一个是预测值pred, 一个是训练数据, 训练数据中包含标签值。在这个函数中, 计算了当前预测值和标签值的损失函数的一阶梯度和二阶导, 并将计算结果返回
#可以看出, 返回的是一阶导和二阶导, 二阶导相当于是Hessian矩阵, 所以此处用了hess作为变量名。很显然, 既然这么设计函数的话, logregobj是可以并行的, 只要每个并行的程序自己调用logregobj即可, 他们之间完全可以互不影响。
def logregobj(preds, dtrain):
    labels = dtrain.get_label()
    preds = 1.0 / (1.0 + np.exp(-preds))
    grad = preds - labels
    hess = preds * (1.0 - preds)
    return grad, hess

#在上面定义完一阶导和二阶导后, 我们此时要计算损失函数, 需要注意的是, 上面的一阶导和二阶导都是矩阵形式的, 不过此处的损失函数的返回值需要是单一的值, 又可称作margin值。回想一下树形结构的公式, 在每个分裂的节点, 该叶节点都对应一个值, 这个值是落在这个叶节点下所有样本的综合值。
def evalerror(preds, dtrain):
    labels = dtrain.get_label()
    # return a pair metric_name, result. The metric name must not contain a colon (:) or a space
    # since preds are margin(before logistic transformation, cutoff at 0)
    return 'my-error', float(sum(labels != (preds > 0.0))) / len(labels)

# training with customized objective, we can also do step by step training
# simply look at xgboost.py's implementation of train
#在obj中传入返回一阶二阶导的函数, 在feval中传入损失函数。
bst = xgb.train(param, dtrain, num_round, watchlist, obj=logregobj, feval=evalerror)

我们回到之前的公式推导环节, 我们已经理解了一阶二阶导, 现在来看正则化项:
$\Omega(f_t)=\gamma T+\frac{1}{2}\lambda\sum\limits_{j=1}^{T}w_j^2$
其实正则为什么可以控制模型复杂度呢？有很多角度可以看这个问题，最直观就是，我们为了使得目标函数最小，自然正则项也要小，正则项要小，叶子节点个数T要小(叶子节点个数少，树就简单)。

而为什么要对叶子节点的值进行L2正则，这个可以参考一下LR里面进行正则的原因，简单的说就是 $\color{red}{LR没有加正则，整个w的参数空间是无限大的，只有加了正则之后，才会把w的解规范在一个范围内。}$ (对此困惑的话可以跑一个不带正则的LR，每次出来的权重w都不一样，但是loss都是一样的，加了L2正则后，每次得到的w都是一样的)

由于前t-1棵树的预测分数与y的残差对目标函数优化不影响, 可以直接去掉, 所以有:
$\hat { \mathcal { L } } ^ { ( t ) } = \sum _ { i = 1 } ^ { n } \left[ g _ { i } f _ { t } \left( \mathbf { x } _ { i } \right) + \frac { 1 } { 2 } h _ { i } f _ { t } ^ { 2 } \left( \mathbf { x } _ { i } \right) \right] + \Omega \left( f _ { t } \right)$
上式是将每个样本的损失函数值加起来，我们知道，每个样本都最终会落到一个叶子结点中，所以我们可以将所有同一个叶子结点的样本重组起来，过程如下:
$\begin{align} \hat{ \mathcal { L }}^{(t)} &=\sum_{i=1}^n [g_if_t(x_i) + \frac12h_tf_t^2(x_i)] + \Omega(f_t) \\ & =\sum_{i=1}^n [g_if_t(x_i) + \frac12h_tf_t^2(x_i)] + \gamma T+\frac{1}{2}\lambda\sum\limits_{j=1}^{T}w_j^2 \\ & = \sum\limits_{j=1}^{T} [(\sum\limits_{i\in I_j}g_i)w_j+\frac{1}{2}(\sum\limits_{i\in I_j}h_i+\lambda) w_j^2]+\gamma T \\ & = \frac{1}{2}\sum\limits_{j=1}^{T} (H_j+\lambda)(w_j + \frac{G_j}{H_j+\lambda})^2+\gamma T -\frac{1}{2}\sum\limits_{j=1}^{T}\frac{G_j^2}{H_j+\lambda} \end{align}$
$\color{red}{其中G_j = \sum\limits_{i\in I_j}g_i为落入叶子 i 所有样本一阶梯度统计值总和, H_j = \sum\limits_{i\in I_j}h_i为落入叶子 i 所有样本二阶梯度统计值总和。}$

因为XGBoost的基本框架是boosting, 也就是一棵树接着一棵树的形式, 所以此处的 $\hat{ \mathcal { L }}^{(t)}$ 中的t代表第t颗树, 注意上式中 j 代表树的一个叶子节点, i 也代表一个叶子节点, 根据树的判别条件, n 个样本点被分到了 T 个叶子节点中

因此通过上式的改写，我们可以将目标函数改写成关于叶子结点分数w的一个一元二次函数，求解最优的w和目标函数值就变得很简单了，直接使用顶点公式即可。因此， $\color{blue}{对于第t颗树而言, 每个叶子节点最优的w为}$ (还记得上面推导出来的牛顿法吗 $\color{red}{\large \frac{f'(x)}{f''(x)}}$ , 是一阶导除以二阶导的形式, $\color{red}{跟此处的 w_j^*的公式联系起来了}$ ):
$w_j^* = -\frac{G_j}{H_j+\lambda}$
当w取上式最优解时, 之前 $\hat{\mathcal { L }}^{t}$ 公式中的第一项等于0, 此时 $\hat{\mathcal { L }}^{*}$ 剩下的项即为最终的目标函数公式:
$\hat{\mathcal { L }}^{*} =-\frac{1}{2}\sum\limits_{j=1}^{T}\frac{G_j^2}{H_j+\lambda}+\gamma T\tag{1}$

我们花了这么大的功夫，得到了叶子结点取值的表达式。在xgboost里，叶子节点取值的表达式很简洁，推导起来也比GBDT的要简便许多。

到这里，我们一直都是在围绕目标函数进行分析，这个到底是为什么呢？这个主要是为了后面我们寻找 $f_k(x)$ ，也就是建树的过程。

具体来说，我们回忆一下建树的时候需要做什么， $\color{red}{建树的时候最关键的一步就是选择一个分裂的准则，也就如何评价分裂的质量}$ 。比如在详解GBDT里，我们可以选择MSE，MAE来评价我们的分裂的质量，但是，我们所选择的分裂准则似乎不总是和我们的损失函数有关，因为这种选择是启发式的。
比如，在分类任务里面，损失函数可以选择logloss，分裂准确选择MSE，这样看来， $\color{red}{似乎分裂的好坏和我们的损失并没有直接挂钩}$ 。

但是， $\color{red}{在xgboost里面，我们的分裂准则是直接与损失函数挂钩的准则}$ ，这个也是xgboost和GBDT一个很不一样的地方。

具体来说，xgboost选择这个准则，计算增益Gain (其中 $G_L$ 代表如果分裂的话左叶子节点中样本点的集合的一阶梯度和, $G_R$ 代表有边界点的):
$Gain=\frac{1}{2}[\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}]-\gamma \tag{2}$

其实选择这个作为准则的原因很简单也很直观。我们这样考虑。由上面的(1)式知道，对于一个结点，假设我们不分裂的话, 那么对于左叶子节点和右叶子节点而言其包含样本点的一阶二阶导的和就不是单独来算了, 而是合在一起。

因为此时所有样本都存在于未分裂的节点中, 所以此时该节点中一阶导的和就相当于 $(G_L+G_R)$ , 未分裂节点上的二阶导为 $(H_L+H_R)$ , 将其代入公式(1)中, 即可得到不分裂情况下的损失值:
$\color{blue}{ \frac{(G_L+G_R)^2}{H_L+H_R+\lambda}}$
假设在这个节点分裂的话，分裂之后左右叶子节点的损失分别为:
$\frac{G_L^2}{H_L+\lambda}\quad和\quad\frac{G_R^2}{H_R+\lambda}$
既然要分裂的时候，我们当然是选择分裂成左右子节点后，损失减少的最多(或者看回(1)式，由于前面的负号，所以欲求(1)的最小值，其实就是求(2)的最大值), 也就是找到一种分裂有:
$\large max\left(\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\right)$

那么 $\gamma$ 的作用是什么呢？利用 $\gamma$ 可以控制树的复杂度，进一步来说，利用 $\gamma$ 来作为阈值，只有大于 $\gamma$ 时候才选择分裂。这个其实起到预剪枝的作用。

最后就是如何得到左右子节点的样本集合？这个和普通的树一样，都是通过遍历特征所有取值，逐个尝试。

至此，我们把xgboost的基本原理阐述了一遍。我们总结一下，其实就是做了以下几件事情。

1.在损失函数的基础上加入了正则项。
2.对目标函数进行二阶泰勒展开。
3.利用推导得到的表达式作为分裂准则, 来构建每一颗树。

手动还原XGBoost实例过程

下面这个例子由https://blog.csdn.net/qq_22238533/article/details/79477547给出。

上面，我只是简单的阐述整个流程，有一些细节的地方可能都说的不太清楚，我以一个简单的UCI数据集，一步一步的演算整个xgboost的过程。数据集如下:

$\begin{bmatrix} \hline ID & x_1 & x_2 & y \\ \hline 1 & 1 & -5 & 0 \\ 2 & 2 & 5 & 0 \\ 3 & 3 & -2 & 1 \\ 4 & 1 & 2 & 1 \\ 5 & 2 & 0 & 1 \\ 6 & 6 & -5 & 1 \\ 7 & 7 & 5 & 1 \\ 8 & 6 & -2 & 0 \\ 9 & 7 & 2 & 0 \\ 10 & 6 & 0 & 1 \\ 11 & 8 & -5 & 1 \\ 12 & 9 & 5 & 1 \\ 13 & 10 & -2 & 0 \\ 14 & 8 & 2 & 0 \\ 15 & 9 & 0 & 1\\ \hline \end{bmatrix}$

这里为了简单起见，树的深度设置为3(max_depth=3)，树的颗数设置为2(num_boost_round=2),学习率为0.1(eta=0.1)。另外再设置两个正则的参数, $\lambda=1,\gamma=0$ 。
损失函数选择logloss。
由于后面需要用到logloss的一阶导数以及二阶导数，这里先简单推导一下:
$\large L(y_i,\hat{y}_i)=y_iln(1+e^{-\hat{y}_i})+(1-y_i)ln(1+e^{\hat{y}_i})$
两边对 $\hat{y}_i$ 求一阶导数, 其中: $\large y_{i,pred}=\frac{1}{1+e^{-\hat{y}_i}}$ :
$\begin{align} \large L'(y_i,\hat{y}_i) &= \large y_i\frac{-e^{-\hat{y}_i}}{1+e^{-\hat{y}_i}}+(1-y_i)\frac{e^{\hat{y}_i}}{1+e^{\hat{y}_i}}\\ &=\large y_i\frac{-1}{1+e^{\hat{y}_i}}+(1-y_i)\frac{1}{1+e^{-\hat{y}_i}}\\ &=\large y_i*(y_{i,pred}-1)+(1-y_i)*y_{i,pred}\\ &=\large y_{i,pred}-y_i \end{align}$
在一阶导的基础上再求一次可得二阶导(其实就是sigmod函数求导):
$\large L’'(y_i,\hat{y}_i)=y_{i,pred}*(1-y_{i,pred})$
所以样本的一阶导数值为:
$\large g_i=y_{i,pred}-y_i\quad\tag{3}$
二阶导数值为:
$\quad\large h_i=y_{i,pred}*(1-y_{i,pred})\tag{4}$
下面就是xgboost利用上面的参数拟合这个数据集的过程:

$\LARGE 建立第一颗树(k=1)$

建树的时候从根节点开始(Top-down greedy)，在根节点上的样本有ID1～ID15。那么回顾xgboost的算法流程，我们需要在根节点处进行分裂，分裂的时候需要计算之前的式子(2), 即如下:
$\operatorname { Gain } = \frac { 1 } { 2 } \left[ \frac { G _ { L } ^ { 2 } } { H _ { L } + \lambda } + \frac { G _ { R } ^ { 2 } } { H _ { R } + \lambda } - \frac { \left( G _ { L } + G _ { R } \right) ^ { 2 } } { H _ { L } + H _ { R } + \lambda } \right] - \gamma$
那么式子(2)表达是：在结点处把样本分成左子结点和右子结点两个集合。分别求两个集合的 $G_L、H_L、G_R、H_R$ ，然后计算增益Gain。

而这里，我其实可以先计算每个样本的一阶导数值和二阶导数值，即按照式子(3)和(4)计算，但是这里你可能碰到了一个问题，那就是第一颗树的时候每个样本的预测的概率 $\large y_{i,pred}$ 是多少？这里和GBDT一样，应该说和所有的Boosting算法一样，都需要一个初始值。而 $\color{red}{在xgboost里面，对于分类任务只需要初始化为(0,1)中的任意一个数都可以}$ 。具体来说就是参数base_score。（其默认值是0.5）

(值得注意的是 $\color{red}{base\_score是一个经过sigmod映射的值}$ ，可以理解为一个概率值，提这个原因是后面建第二颗树会用到，需要注意这个细节）

这里我们也设base_score=0.5。然后我们就可以计算每个样本的一阶导数值和二阶导数值了。具体如下表:

一阶导和二阶导

比如说对于ID=1样本, 由上面得到的一阶二阶导公式(3)(4)可得:

$g_1=y_{1,pred}-y_1=0.5-0=0.5$
$h_1=y_{1,pred}*(1-y_{1,pred})=0.5*(1-0.5)=0.25$

那么把样本如何分成两个集合呢？这里就是上面说到的选取一个最佳的特征以及分裂点使得Gain最大。
比如说对于特征x1，一共有[1,2,3,6,7,8,9,10]8种取值(注意上表的ID不代表取值)。可以得到以下这么多划分方式。

1)以1为划分时 $(x_1<1$ ):
左子节点包含的样本: $I_L=[]$
右子节点包含的样本: $I_R=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]$
那么左子节点为空， $G_L=0\ 和\ H_L=0$
右子节点的一阶导数和: $G_R=\sum_{(i \in I_R)}g_i=(0.5+0.5+....-0.5)=-1.5$
右子节点的二阶导数和: $H_R=\sum_{(i \in I_R)}h_i=(0.25+0.25...0.25)=3.75$
最后我们在计算一下增益Gain，可以有得到Gain=0。
计算出来发现Gain=0,不用惊讶,因为左子结点为空，也就是这次分裂把全部样本都归到右子结点，这个和没分裂有啥区别？所以，分裂的时候每个结点至少有一个样本。

$\color{red}{在上面以1为划分标准时, 我们注意到一阶和二阶导是直接用的上表中提前计算好的, 而不用重新计算, 样本被分到哪部分, 其已经计算好的一阶二阶导就直接累积到那部分上即可。}$

2)以2为划分时 $(x_1<2$ ):
左子节点包含的样本: $I_L=[1,4]$
右子节点包含的样本: $I_R=[2,3,5,6,7,8,9,10,11,12,13,14,15]$
左子节点的一阶导数和: $G_L=\sum_{(i \in I_L)}g_i=(0.5-0.5)=0$
左子节点的二阶导数和: $H_L=\sum_{(i \in I_L)}h_i=(0.25+0.25)=0.5$
右子节点的一阶导数和: $G_R=\sum_{(i \in I_R)}g_i=-1.5$
右子节点的二阶导数和: $H_R=\sum_{(i \in I_R)}h_i=3.25$
$\operatorname { Gain } = \frac { 1 } { 2 } \left[ \frac { G _ { L } ^ { 2 } } { H _ { L } + \lambda } + \frac { G _ { R } ^ { 2 } } { H _ { R } + \lambda } - \frac { \left( G _ { L } + G _ { R } \right) ^ { 2 } } { H _ { L } + H _ { R } + \lambda } \right] - \gamma=0.0557275541796$

3)以3为划分时 $(x_1<3$ ):
左子节点包含的样本: $I_L=[1,2,4,5]$
右子节点包含的样本: $I_R=[3,6,7,8,9,10,11,12,13,14,15]$
左子节点的一阶导数和: $G_L=\sum_{(i \in I_L)}g_i=0$
左子节点的二阶导数和: $H_L=\sum_{(i \in I_L)}h_i=1$
右子节点的一阶导数和: $G_R=\sum_{(i \in I_R)}g_i=-1.5$
右子节点的二阶导数和: $H_R=\sum_{(i \in I_R)}h_i=2.75$
$\operatorname { Gain } = \frac { 1 } { 2 } \left[ \frac { G _ { L } ^ { 2 } } { H _ { L } + \lambda } + \frac { G _ { R } ^ { 2 } } { H _ { R } + \lambda } - \frac { \left( G _ { L } + G _ { R } \right) ^ { 2 } } { H _ { L } + H _ { R } + \lambda } \right] - \gamma=0.126315789474$

其他的值[6,7,8,9,10]类似，计算归总到下面表:

计算汇总

从上表我们可以到，如果特征x1以x1<10分裂时可以得到最大的增益0.615205。

按照算法的流程，这个时候需要遍历下一个特征x2，对于特征x2也是重复上面这些步骤，可以得到类似的表如下:

特征x2

可以看到，以x2特征来分裂时，最大的增益是0.2186<0.615205。所以在根节点处，我们以x1<10来进行分裂。

由于我设置的最大深度是3，此时只有1层，所以还需要继续往下分裂。
左子节点的样本集合 $I_L=[1,2,3,4,5,6,7,8,9,10,11,12,14,15]$
右子节点的样本集合 $I_R=[13]$
右子节点此时只剩一个样本，不需要分裂了，也就是已经是叶子结点。 $\color{red}{可以计算其对应的叶子结点值了}$ ，按照如下公式:
$w^*=-\frac{G_j}{H_j+\lambda}$
有:
$w_1=-\frac{G_R}{H_R+\lambda}=-\frac{g_{13}}{h_{13}+1}=-\frac{0.5}{1+0.25}=-0.4$
那么下面就是对左子结点 $I_L$ 进行分裂。分裂的时候把此时的结点看成根节点，其实就是循环上面的过程，同样也是需要遍历所有特征 $(x_1,x_2)$ 的所有取值作为分裂点，选取增益最大的点。
这里为了说的比较清晰，也重复一下上面的过程:

此时样本有: $I=[1,2,3,4,5,6,7,8,9,10,11,12,14,15]$ , 注意, $\color{red}{此时已经没有ID为13的样本了, 因为已经被上一层树的右子节点分出去了。}$

先考虑特征 $x_1$ ，此时 $x_1$ 的取值有[1, 2, 3, 6, 7, 8, 9]

注意, 现在用的还是之前求出来的一阶二阶导, 并没有重新计算, 因为此时是在构建当前这棵树的第二层, 而不是构建新树, 在构建新树时才会。

1)以1为分裂Gain=0。
2)以2为划分时 $(x_1<2$ ):
左子节点包含的样本: $I_L=[1,4]$
右子节点包含的样本: $I_R=[2,3,5,6,7,8,9,10,11,12,14,15]$
左子节点的一阶导数和: $G_L=\sum_{(i \in I_L)}g_i=(0.5-0.5)=0$
左子节点的二阶导数和: $H_L=\sum_{(i \in I_L)}h_i=(0.25+0.25)=0.5$
右子节点的一阶导数和: $G_R=\sum_{(i \in I_R)}g_i=-2$
右子节点的二阶导数和: $H_R=\sum_{(i \in I_R)}h_i=3$
$\operatorname { Gain } = \frac { 1 } { 2 } \left[ \frac { G _ { L } ^ { 2 } } { H _ { L } + \lambda } + \frac { G _ { R } ^ { 2 } } { H _ { R } + \lambda } - \frac { \left( G _ { L } + G _ { R } \right) ^ { 2 } } { H _ { L } + H _ { R } + \lambda } \right] - \gamma=0.111111111111$

其他的同理，最后所有值都遍历完后可以得到下表:

对所有值遍历后

可以看到，x1选择x<3时能够获得最大的增益0.2539。

同理，我们对x2再次遍历可以得到下表：

对x2再次遍历

可以看到 $x_2$ 在 $x_2<2$ 时分裂可以获得最大的增益0.4444。
比较知道，应该选择 $x_2<2$ 作为分裂点。
分裂后左右叶子结点的集合如下：
左子节点的样本集合 $I_L=[1,3,5,6,8,10,11,15]$
右子节点的样本集合 $I_R=[2,4,7,9,12,14]$
然后继续对 $I_L$ 和 $I_R$ 进行分裂。这里就不在啰嗦了。这里直接给出第一个树的结构。

第一颗树

这里有的人可能对叶子结点取值感到困惑。为何我算出来的是-0.4，可图上却是-0.04？
这里以最左边的一个叶子结点为例，展示一下计算的过程。
落在最左边这个叶子结点上的样本有 $I=[1]$

所以由公式:
$w^*=-\frac{G_j}{H_j+\lambda}$
可得:
$w_2=-\frac{G_R}{H_R+\lambda}=-\frac{g_{1}}{h_{1}+1}=-\frac{0.5}{1+0.25}=-0.4$
落在从左边数起第二个叶子结点上的样本有 $I=[3,5,6,8,10,11,15]$
$w_3=-\frac{G_R}{H_R+\lambda}=-\frac{g_{3}+g_{5}+g_{6}+g_{8}+g_{10}+g_{11}+g_{15}}{h_{3}+h_{5}+h_{6}+h_{8}+h_{10}+h_{11}+h_{15}+1}=-\frac{-2.5}{1+1.75}=0.909$

到这里完全没有问题，但是为什么和图上的不一样呢？这里其实和我们在GBDT中的处理一样， $\color{red}{我们会以一个学习率来乘这个值}$ ，当完全取-0.4时说明学习率取1，这个时候很容易过拟合。所以每次得到叶子结点的值后需要乘上学习率eta，在前面我们已经设置了学习率是0.1。这里也是GBDT和xgboost一个共同点，大家都是通过学习率来进行Shrinkage， $\color{red}{以减少过拟合的风险}$ 。

至此，我们学习完了第一颗树。

$\LARGE建立第2颗树(k=2)$

这里，我们开始拟合我们第二颗树。其实过程和第一颗树完全一样。只不过 $\color{red}{对于y_{i,pred}需要进行更新}$ ，也就是拟合第二颗树是在第一颗树预测的结果基础上。这和GBDT一样，因为大家都是Boosting思想的算法。

在第一颗树里面由于前面没有树，所以初始yi,pred=0.5 (用户自己设置的)

假设此时，模型只有这一颗树(K=1)，那么模型对样例xi进行预测时，预测的结果表达是什么呢？
由加法模型
$y_i^K=\sum_{k=0}^Kf_k(x_i)$
有:
$y_i^1=f_0(x_i)+f_1(x_i)$
其中 $f_1(x_i)$ 的值是样例 $x_i\color{red}{落在第一棵树上的叶子结点值}$ 。那 $f_0(x_i)$ 是什么呢？这里就涉及前面提到一个问题 $\color{red}{base\_score是一个经过sigmod映射后的值}$ (因为选择使用Logloss做损失函数，概率 $p=\large\frac{1}{1+e^{−x}}$ )

所以需要将0.5逆运算 $x=ln\frac{y}{1-y}$ 后可以得到 $f_0(x_i)=0$ 。

所以第一颗树预测的结果就是 $y_i^1=f_0(x_i)+f_1(x_i)=0+\large\color{red}{ w_{q(x_i)}}$ (其实当训练次数K足够多的时候，初始化这个值几乎不起作用的，这个在官网文档上有说明)

比如对于ID=1的样本，其落在-0.04这个节点。 $\color{red}{那么经过sigmod映射后的值为}$ (预测后将其映射成概率):
$\color{red}{p_{1,pred}=\frac{1}{1+e^{(0+0.04)}}=0.490001}$

所以，我们可以得到第一棵树预测的结果为下表, $\color{red}{注意此处得到的是y_{i,pred}的值, 其和p_{1,pred}是等同的}$ :
$\begin{bmatrix} \hline ID & y_{i,pred}\\ \hline 1 & 0.490001 \\ 2 & 0.494445 \\ 3 & 0.522712 \\ 4 & 0.494445 \\ 5 & 0.522712 \\ 6 & 0.522712 \\ 7 & 0.494445 \\ 8 & 0.522712 \\ 9 & 0.494445 \\ 10 & 0.522712 \\ 11 & 0.522712 \\ 12 & 0.509999 \\ 13 & 0.490001 \\ 14 & 0.494445 \\ 15 & 0.522712 \\ \hline \end{bmatrix}$
由之前计算的一阶导和二阶导的公式, 注意此处的 $y_i$ 是最开始的UCI数据集中的样本的标签值:
$\large g_i=y_{i,pred}-y_i\quad \quad\large h_i=y_{i,pred}*(1-y_{i,pred})$
我们可以接着得到新一轮所有样本的一阶和二阶导, 具体如下表, 注意此处之前在第一棵树中被分出去的譬如ID为13的样本, $\color{red}{此时又重新放回训练数据中}$ , 因为此时是第二颗树的训练:
$\begin{bmatrix} \hline ID & g_i & h_i\\ \hline 1 & 0.490001320839 & 0.249900026415 \\ 2 & 0.494444668293 & 0.24996913829 \\ 3 & -0.477288365364 & 0.249484181652 \\ 4 & -0.505555331707 & 0.24996913829 \\ 5 & -0.477288365364 & 0.249484181652 \\ 6 & -0.477288365364 & 0.249484181652 \\ 7 & -0.505555331707 & 0.24996913829 \\ 8 & 0.522711634636 & 0.249484181652 \\ 9 & 0.494444668293 & 0.24996913829 \\ 10 & -0.477288365364 & 0.249484181652 \\ 11 & -0.477288365364 & 0.249484181652 \\ 12 & -0.490001320839 & 0.249900026415 \\ 13 & 0.490001320839 & 0.249900026415 \\ 14 & 0.494444668293 & 0.24996913829 \\ 15 & -0.477288365364 & 0.249484181652 \\ \hline \end{bmatrix}$
之后，我们和第一颗树建立的时候一样以公式(2), 即下式去建树。
$Gain=\frac{1}{2}[\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{(G_L+G_R)^2}{H_L+H_R+\lambda}]-\gamma$
拟合完后第二颗树如下图:

第二颗树

后面的所有过程都是重复这个过程，这里就不再啰嗦了。

至此，整个xgboost的训练过程已经完了，但是其实里面还有一些细节的东西，下面已单独一个部分来说明这个部分。

训练过程的细节-参数min_child_weight

在选择分裂的时候，我们是选取分裂点为一个最大的增益Gain。但是其实在xgboost里面有一个参数叫min_child_weight。

我先来看看官网对这个参数的解释：

min_child_weight

可能看完大概知道这个权重指的应该是二阶导数，但是具体是怎么一回事呢。
其实是这样的：
在前面建立第一颗的根结点的时候，我们得到特征x1每个分裂点的这个表：

x1分裂点

我们当时选取了 $x_1<10$ 作为分裂特征，但是！这个是有一个前提的，那就是 $\color{red}{参数 min\_child\_weight<min(H_L,H_R)}$
如果我们设置min_child_weight=0.26的时候，分裂点就不是选择10，而是 $\color{red}{放弃这个最大增益}$ ，考虑次最大增益。如果次最大增益也不满足min_child_weight< $min(H_L,H_R)$ ，则继续往下找，如果没有找到一个满足的，则不进行分裂。（在上面中min_child_weight取的0，所以只要是最大的增益就选为分裂点）

$\LARGE 训练过程细节-参数\gamma$
在前面例子里，我们把 $\gamma$ 设成了0，如果我们设置成其他值比如1的话， $\color{red}{在考虑最大增益的同时，也要考虑这个最大的增益是否比\gamma大}$ , 如果小于 $\gamma$ 则不进行分裂(预剪枝)。

$\LARGE 训练过程细节-缺失值的处理$
我前面放了许多树的结构图，上面有一个方向是missing，其实这个是xgboost自动学习的一个对缺失值处理的方向。可以看到，所有的missing方向都是指向右子结点，这是因为我们数据集中没有缺失值。

xgboost对缺失值的处理思想很简单，具体看下面的算法流程:

缺失值处理

首先需要注意到两个集合一个是 $I$ ,其包含所有的样本（包括含空缺值的样本）。
$I_k$ 是不包含空缺值样本的集合。
$\color{red}{在计算总的G和H时用的是\ I,\ 也就说空缺值的一阶导数和二阶导数已经包含进去了。}$

可以看到内层循环里面有两个for，第一个for是从把特征取值从小到大排序，然后从小到大进行扫描，这个时候 $\color{red}{在计算G_R的时候是用总的G减去G_L, H_R也是同样用总的H减去H_L,这意味着把空缺样本归到了右子结点。}$

$\color{red}{第二个for相反过来，把空缺样本归到了左子结点。}$

$\color{red}{只要比较这两次最大增益出现在第一个for中还是第二个for中就可以知道对于空缺值的分裂方向}$ , 这就是xgboost如何学习空缺值的思想。

$\LARGE 训练过程细节-近似贪婪算法$

在前面，我们可以看到选取特征的时候是遍历了每个特征所有可能的取值，但是实际上，我们可以使用近似的方法来进行特征分裂的选择，具体算法流程:

近似贪婪算法

xgboost如何用于特征选择？

相信很多做过数据挖掘比赛的人都利用xgboost来做特征选择。
一般我们调用xgb库的get_fscore()。但其实xgboost里面有三个指标用于对特征进行评价，而get_fscore()只是其中一个指标weight。

这个指标大部分人都很熟悉，其代表着某个特征被选作分裂的次数。
比如在前面举的例子里，我们得到这两颗树:

树1

树2

可以看到特征x1被选作分裂点的次数为6 (在同一棵树中可以反复选择一个特征)，x2被选做分裂点的次数为2。
get_fscore()就是返回这个指标。

而xgboost还提供了另外两个指标，一个叫gain，一个叫cover。可以利用get_score()来选择。

那么gain是指什么呢？其代表着某个特征的平均增益。

比如，特征 $x_1$ 被选了6次作为分裂的特征，每次的增益假如为 $Gain_1,Gain_2,…Gain_6$ ，那么其平均增益为 $(Gain_1+Gain_2+...Gain_3)/6$

最后一个cover指的是什么呢？其代表着每个特征在分裂时结点处的平均二阶导数。
这个为了加深理解，我举个例子，这个例子用的还是UCI数据集，不过此时只有max_depth=2,num_boost_round=1

最大深度为2的树

建树完之后，其结构如上。
在第一个结点分裂时，x1的特征增益表如下：

x1的特征增益表

第二个结点分裂时，x2的特征增益表如下：

x2的特征增益表

那么特征x1的cover是如下计算的：
在x1在第一个结点处被选择分裂特征，此时结点处的总二阶导数 G=3.75

故x1的cover值为3.75。
这里x1只是被分裂了一次，如果后续还有就是求平均。

同样地, x2的cover值为3.5。

举这个例子就是先给大家说一下何为平均二阶导数。

其实为什么需要选择二阶导数？这个二阶导数背后有什么意义吗？我们先看看官网对cover的定义：
‘cover’ - the average coverage of the feature when it is used in trees。大概的意义就是特征覆盖了多少个样本。

这里，我们不妨假设损失函数是mse，也就是 $0.5*(y_i−y_{pred})^2$ ,我们求其二阶导数，很容易得到为常数1。也就是每个样本对应的二阶导数都是1，那么这个cover指标不就是意味着，在某个特征在某个结点进行分裂时所覆盖的样本个数吗？

xgboost与传统GBDT的区别与联系？

至此，我们来简单的总结一下xgboost和GBDT的区别以及联系。

区别：
1.xgboost和GBDT的一个区别在于目标函数上。
在xgboost中，损失函数+正则项。 GBDT中一般只有损失函数。
2.xgboost中利用二阶导数的信息，而GBDT只利用了一阶导数, 即在GBDT回归中利用了残差的概念。
3.xgboost在建树的时候利用的准则来源于目标函数推导，即可以理解为牛顿法。而GBDT建树利用的是启发式准则。
4.xgboost中可以自动处理空缺值，自动学习空缺值的分裂方向，GBDT(sklearn版本)不允许包含空缺值。

相似点：
1.xgboost和GBDT的学习过程都是一样的，都是基于Boosting的思想，先学习前n-1个学习器，然后基于前n-1个学习器学习第n个学习器。而且其都是将损失函数和分裂点评估函数分开了。
2.建树过程都利用了损失函数的导数信息,。
3.都使用了学习率来进行Shrinkage，从前面我们能看到不管是GBDT还是xgboost，我们都会利用学习率对拟合结果做缩减以减少过拟合的风险。

参考资料:
https://blog.csdn.net/qq_22238533/article/details/79477547,

https://homes.cs.washington.edu/~tqchen/pdf/BoostedTree.pdf

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 158,736评论 4赞 362
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 67,167评论 1赞 291
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 108,442评论 0赞 243
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 43,902评论 0赞 204
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 52,302评论 3赞 287
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 40,573评论 1赞 216
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 31,847评论 2赞 312
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 30,562评论 0赞 197
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 34,260评论 1赞 241
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 30,531评论 2赞 245
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 32,021评论 1赞 258
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 28,367评论 2赞 253
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 33,016评论 3赞 235
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 26,068评论 0赞 8
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 26,827评论 0赞 194
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 35,610评论 2赞 274
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 35,514评论 2赞 269