52ML系列(3)--Xgboost详解

本来觉得xgboost已经弄懂了，但听了AI Lab陈博士的讲座之后，又有了更深入的认知，本文将详细解释一些细节，帮助大家理解。顺便表示对陈博的膜拜，讲的太清楚了👍。

首先呢，Xgboost是一个究极进化体，它的祖先其实是棵树。我们来看下它的进化史：
$\begin{align} X\!G\!Boost&=eXtreme+GBDT\\ &=eXtreme+(Gradient+BDT) \\ &=eXtreme+Gradient+(Boosting+DecisionTree) \end{align}$

$Boosting \to BDT \to GBDT \to X\!G\!Boost$

本文为了让你真正懂Xgboost，将追根溯源，从他的祖先开始一层一层的讲解Xgboost是怎么进化来的😝

But！本文默认你是懂DecisionTree的，不清楚决策树是啥的请自行百度😌

But，but！还是要提一嘴决策树的表示形态，这个很重要，因为Xgboost的进化与决策树的形态变化密不可分。

决策树有哪些表示形态啊？

树形结构
if-else判断结构
图像区域划分
。。。

前三种是我们常见的形式了，点点点我们后面再说😝

好！我们正式开始研究Xgboost的进化史，首先是Boosting。

Boosting

提升方法(boosting)就是集成方法中的一员(另一员是Bagging，其中差别不在本文讨论范围内，请自行百度)，它将所有臭皮匠加起来，最后赛过一个诸葛亮。所以boosting本质是一个加法模型。

加法模型：
$f\left(x\right)=\sum_{m=1}^M\beta_m b\left(x;\gamma_m\right) \tag{1.1}$
其中， $b\left(x;\gamma_m\right)$ 为基函数， $\gamma_m$ 为基函数的参数， $\beta_m$ 为基函数的系数。 $b\left(x;\gamma_m\right)$ 基函数是什么？只要保证函数可加性就可以，比如决策树。

有了这个加法模型，那么我们的目标就是：在给定训练数据 $\{\left(x_i,y_i\right)\}_{i=1}^N$ 及损失函数 $L\left(y,f\left(x\right)\right)$ 的条件下，学习加法模型 $f\left(x\right)$ 成为经验风险极小化问题：
$\min_{\beta_m,\gamma_m}\sum_{i=1}^N L\left(y_i,\sum_{m=1}^M\beta_m b\left(x_i;\gamma_m\right)\right)\tag{1.2}$
稍微解释下经验风险极小化是什么鬼，其实就是损失函数。还有一个结构风险极小化，指的是正则化项。

式（1.2）这个模型看上去有点麻烦，但我们可以从前向后，每一步只学习一个基函数及其系数，逐步逼近优化目标函数式（1.2），这样就可以简化优化复杂度。具体地，每步只需优化如下损失函数：
$\min_{\beta,\gamma}\sum_{i=1}^N L\left(y_i,\beta b\left(x_i;\gamma\right)\right)\tag{1.3}$

于是有了前向分布算法：

算法1：前向分步算法

输入：训练数据集 $T=\{\left(x_1,y_1\right),\left(x_2,y_2\right),\dots,\left(x_N,y_N\right)\}$ ；损失函数 $L\left(y,f\left(x\right)\right)$ ；基函数集合 $\{b\left(x;\gamma\right)\}$ ；
输出：加法模型 $f\left(x\right)$
（1）初始化 $f_0\left(x\right)=0$
（2）对 $m=1,2,\dots,M$
（a）极小化损失函数
$\left(\beta_m,\gamma_m\right)=\mathop{\arg\min}_{\beta,\gamma} \sum_{i=1}^N L\left(y_i, f_{m-1}\left(x_i\right)+\beta b\left(x_i;\gamma\right)\right) \tag{1.4}$
得到参数 $\beta_m$ ， $\gamma_m$
（b）更新
$f_m\left(x\right)=f_{m-1}\left(x\right)+\beta_m b\left(x;\gamma_m\right) \tag{1.5}$
（3）得到加法模型
$f\left(x\right)=f_M\left(x\right)=\sum_{m=1}^M\beta_m b\left(x;\gamma_m\right) \tag{1.6}$

前向分步算法将同时求解从 $m=1$ 到 $M$ 所有参数 $\beta_m,\gamma_m$ 的优化问题简化为逐次求解各个 $\beta_m, \gamma_m$ 的优化问题。

以上就是Boosting

好啦，接下来要开始进化了：boosting->boosting+decision tree

BDT

我们把boosting的基函数设置成决策树，那么boosting就变成了BDT，提升决策树。
提升决策树模型可以表示为决策树的加法模型：
$f_M=\sum_{m=1}^M T\left(x;\Theta_m\right) \tag{2.1}$
其中， $T\left(x;\Theta_m\right)$ 表示决策树； $\Theta_m$ 为决策树的参数； $M$ 为树的个数。

还记得前文说的决策树形式吗？这下又多了一种表示形式 $T\left(x;\Theta_m\right)$ 。下面是对这种表示形式的解释：

已知训练数据集 $T=\{\left(x_1,y_1\right),\left(x_2,y_2\right),\dots\left(x_N,y_N\right)\}$ ， $x_i\in\mathcal{X}\subseteq\mathbb{R}^n$ ， $\mathcal{X}$ 为输入空间， $y_i\in\mathcal{Y}\subseteq\mathbb{R}$ ， $\mathcal{Y}$ 为输出空间。如果将输入空间 $\mathcal{X}$ 划分为 $J$ 个互不相交的区域 $R_1,R_2,\dots,R_J$ ，并且在每个区域上确定输出的常量 $c_j$ ，那么决策树可表示为
$T\left(x;\Theta\right)=\sum_{j=1}^J c_j I\left(x\in R_j\right) \tag{2.4}$
其中，参数 $\Theta=\{\left(R_1,c_1\right),\left(R_2,c_2\right),\dots,\left(R_J,c_J\right)\}$ 表示决策树的区域划分和各区域上的常量值。 $J$ 是决策树的复杂度即叶子结点个数。

我们回到BDT上来，提升决策树使用以下前向分步算法：
$\begin{align} f_0\left(x\right)&=0 \\ f_m\left(x\right)&=f_{m-1}\left(x\right)+T\left(x;\Theta_m\right),\quad m=1,2,\dots,M \\ f_M\left(x\right)&=\sum_{m=1}^M T\left(x;\Theta_m\right) \end{align}$
在前向分步算法的第 $m$ 步，给定当前模型 $f_{m-1}\left(x\right)$ ，需要求解
$\hat{\Theta}_m=\mathop{\arg\min}_{\Theta_m}\sum_{i=1}^N L\left(y_i,f_{m-1}\left(x_i\right)+T\left(x_i;\Theta_m\right)\right)$
得到 $\hat{\Theta}_m$ ，即第 $m$ 棵树的参数。

当采用平方误差损失函数时，
$L\left(y,f\left(x\right)\right)=\left(y-f\left(x\right)\right)^2$
其损失变为
$\begin{align} L\left(y,f_{m-1}\left(x\right)+T\left(x;\Theta_m\right)\right) &=\left[y-f_{m-1}\left(x\right)-T\left(x;\Theta_m\right)\right]^2 \\ &=\left[r-T\left(x;\Theta_m\right)\right]^2 \end{align}$
其中，
$r=y-f_{m-1}\left(x\right) \tag{2.5}$
是当前模型拟合数据的残差（residual）。对回归问题的提升决策树，只需要简单地拟合当前模型的残差。

算法2：回归问题的提升决策树算法

输入：训练数据集 $T=\{\left(x_1,y_1\right),\left(x_2,y_2\right),\dots,\left(x_N,y_N\right)\}$ ；
输出：提升决策树 $f_M\left(x\right)$
（1）初始化 $f_0\left(x\right)=0$
（2）对 $m=1,2,\dots,M$
（a）按照式（2.5）计算残差
$r_{mi}=y_i-f_{m-1}\left(x_i\right), \quad i=1,2,\dots,N$
（b)拟合残差 $r_{mi}$ 学习一个回归树，得到 $T\left(x;\Theta_m\right)$
（c）更新 $f_m\left(x\right)=f_{m-1}\left(x\right)+T\left(x;\Theta_m\right)$
（3）得到回归提升决策树
$f_M\left(x\right)=\sum_{m=1}^M T\left(x;\Theta_m\right)$

简单总结下，BDT其实就是把boosting的基函数设定成了决策树而已。

形象的说说Boosting和BDT：
诸葛亮的智力值为100，有个智力值为40的臭皮匠A想打败他。两者差了60点智力值，实力太悬殊，于是A又找了个智力值为30的臭皮匠B，这下A、B加起来智力值就是70了，只差诸葛亮30点了。于是他们又找到了智力值为20的臭皮匠C。。。好吧，你们慢慢找，总有一天能够打败诸葛亮。

好啦，又要进化了：BDT->GBDT

GBDT

前文的加法模型中，新的基模型都是在拟合真实残差，既 $r_{mi}=y_i-f_{m-1}\left(x_i\right), \quad i=1,2,\dots,N$
如何加速这个拟合过程呢？我们知道负导数是函数下降最快的方向，那么对损失函数求偏导就能得到极值。

GBDT使用损失函数的负梯度在当前模型的值
$-\left[\frac{\partial L\left(y,f\left(x_i\right)\right)}{\partial f\left(x_i\right)}\right]_{f\left(x\right)=f_{m-1}\left(x\right)} \tag{3.1}$
作为回归问题提升决策树算法中残差的近似值，拟合一个回归树。

算法3：梯度提升算法

输入：训练数据集 $T=\{\left(x_1,y_1\right),\left(x_2,y_2\right),\dots,\left(x_N,y_N\right)\}$ ；损失函数 $L\left(y,f\left(x\right)\right)$
输出：梯度提升决策树 $\hat{f}\left(x\right)$
（1）初始化
$f_0\left(x\right)=\mathop{\arg\min}_c\sum_{i=1}^N L\left(y_i,c\right)$
（2）对 $m=1,2,\dots,M$
（a）对 $i=1,2,\dots,N$ ，计算
$r_{mi}=-\left[\frac{\partial L\left(y,f\left(x_i\right)\right)}{\partial f\left(x_i\right)}\right]_{f\left(x\right)=f_{m-1}\left(x\right)}$
（b)对 $r_{mi}$ 拟合一个回归树，得到第 $m$ 棵树的叶结点区域 $R_{mj},j=1,2,\dots,J$
（c）对 $j=1,2,\dots,J$ ，计算
$c_{mj}=\mathop{\arg\min}_c\sum_{x_i\in R_{mj}} L\left(y_i, f_{m-1}\left(x_i\right)+c\right)$
（d）更新 $f_m\left(x\right)=f_{m-1}\left(x\right)+\sum_{j=1}^J c_{mj} I\left(x\in R_{mj}\right)$
（3）得到回归梯度提升决策树
$\hat{f}\left(x\right)=f_M\left(x\right)=\sum_{m=1}^M \sum_{j=1}^J c_{mj} I\left(x\in R_{mj}\right)$

这里和BDT有几处不同：

第一步的基分类器不再是粗暴的设置为 $f_0\left(x\right)=0$ ，而是拟合了一个当前最优解c。
拟合的残差是损失函数的负梯度。
我们观察式（3.1），可以发现它是一个二元函数。与常见的逻辑回归、svm中的梯度下降不太一样的是，前者求得是参数 $w$ ，而这里我们直接把 $f\left(x\right)$ 作为变量求解。
决策树的表示形式又不一样了，这次变成了 $\sum_{j=1}^J c_{mj} I\left(x\in R_{mj}\right)$ ，这样表示的意思是以树的叶节点作为考量，样本放入树模型后，被划分到哪个叶节点？而该叶节点又将输出何值？

终于进化到究极体了！GBDT->Xgboost

Xgboost

在Xgboost里，我们又双叒叕重新改变了决策树的表示形式：
训练数据集 $\mathcal{D}=\{\left(\mathbf{x}_i,y_i\right)\}$ ，其中 $\mathbf{x}_i\in\mathbb{R}^m,y_i\in\mathbb{R},\left|\mathcal{D}\right|=n$ ，则决策树表示为：
$f\left(\mathbf{x}\right)=w_{q\left(\mathbf{x}\right)} \tag{4.1}$
其中， $q:\mathbb{R}^m\to \{1,\dots,T \},w\in\mathbb{R}^T$ , $T$ 为决策树叶子节点数。

这里是第一个重点，需要解释下。
首先我们对决策树来个功能上的认知：给树模型一个样本，树模型就输出一个预测值，对吧？这个预测值实际上就是叶子节点的分数，也就是说一个叶子节点对应一个预测值。
那么这里的 $q$ 函数表示的是叶子节点的序号， $q\left(\mathbf{x}\right)$ 的意思就是样本 $x$ 给到 $q$ 函数之后， $q$ 函数将样本分配到某个叶子节点上，并输出这个叶子节点的序号。
而 $w$ 表示的是叶子节点的分数，用下标（就是 $q(x)$ ）的形式来表示是哪一个叶子节点的分数。

理解了究极体的树模型表示形式，接下来我们来看看Xgboost模型：
$\hat{y}_i=\phi\left(\mathbf{x}_i\right)=\sum_{k=1}^K f_k\left(\mathbf{x}_i\right) \tag{4.2}$
其中， $f_k\left(\mathbf{x}\right)$ 为第 $k$ 棵决策树。

这种模型跟前文的BDT、GBDT看起来并没有什么区别，但是基于全新的树模型表示方式，我们可以把正则化项加到目标函数中去：
$\mathcal{L}\left(\phi\right)=\sum_i l\left(\hat{y}_i,y_i\right)+\sum_k \Omega\left(f_k\right) \tag{4.3}$
其中， $\Omega\left(f\right)=\gamma T+\frac{1}{2}\lambda\|w\|^2=\gamma T+\frac{1}{2}\lambda\sum_{j=1}^T w_j^2$ 。

这样第 $t$ 轮目标函数就是这样：
$\mathcal{L}^{\left(t\right)}=\sum_{i=1}^n l\left(y_i,\hat{y}^{\left(t-1\right)}_i+f_t\left(\mathbf{x}_i\right)\right)+\Omega\left(f_t\right) \tag{4.4}$

这里是第二个重点，在目标函数里加入了GBDT没有考虑的正则项，降低了模型过拟合风险。

接下来是第三个重点，GBDT采用损失函数负梯度（一阶求导）作为残差，Xgboost嫌他慢了，采用二阶求导。为了让损失函数具有一般性，这里要应用泰勒展开公式，什么是泰勒公式自行百度😈。

这里的损失函数是二元的，我们只对 $\hat{y}$ 做变换。
第 $t$ 轮目标函数在 $\hat{y}^{\left(t-1\right)}$ 处的二阶泰勒展开得到：
$\mathcal{L}^{\left(t\right)}\simeq\sum_{i=1}^n\left[l\left(y_i,\hat{y}^{\left(t-1\right)}\right)+g_i f_t\left(\mathbf{x}_i\right)+\frac{1}{2}h_i f^2_t\left(\mathbf{x}_i\right)\right]+\Omega\left(f_t\right) \tag{4.5}$
其中， $g_i=\partial_{\hat{y}^{\left(t-1\right)}}l\left(y_i,\hat{y}^{\left(t-1\right)}\right),h_i=\partial^2_{\hat{y}^{\left(t-1\right)}}l\left(y_i,\hat{y}^{\left(t-1\right)}\right)$ 。

第 $t$ 轮目标函数的二阶泰勒展开并移除关于 $f_t\left(\mathbf{x}_i\right)$ 常数项，得到：
$\begin{align} \tilde{\mathcal{L}}^{\left(t\right)}&=\sum_{i=1}^n\left[g_if_t\left(\mathbf{x}_i\right)+\frac{1}{2}h_i f^2_t\left(\mathbf{x}_i\right)\right]+\Omega\left(f_t\right) \\ &=\sum_{i=1}^n\left[g_i f_t\left(\mathbf{x}_i\right)+\frac{1}{2}h_i f^2_t\left(\mathbf{x}_i\right)\right]+\gamma T+\frac{1}{2}\lambda\sum_{j=1}^T w_j^2 \end{align} \tag{4.6} \\$

这里稍微解释下，我们的目标是求出让损失函数最小的 $f_t(x)$ ，而 $\hat{y}^{(t-1)}$ 是在上一轮已经求出的已知项，因此 $l\left(y_i,\hat{y}^{\left(t-1\right)}\right)$ 在这里是常数项。

式（4.6）左边是用样本点的标号来遍历，而右边是用叶子节点的标号遍历，能不能统一一下呢？可以啊，就全部以叶子标号来遍历好了：

定义叶结点 $j$ 上的样本的下标集合 $I_j=\{i|q\left(\mathbf{x}_i\right)=j\}$ ，则目标函数可表示为按叶结点累加的形式
$\tilde{\mathcal{L}}^{\left(t\right)}=\sum_{j=1}^T\left[\left(\sum_{i\in I_j}g_i\right)w_j+\frac{1}{2}\left(\sum_{i\in I_j}h_i+\lambda\right)w_j^2\right]+\gamma T \tag{4.7}$

这一步是第四个重点，一个很有技巧的操作，也是很多人看不明白的地方，我们来理解下：
我们的本来目的是要遍历所有的样本点对吧？而样本点又无重复的分布在所有叶子节点上对吧？那我们遍历所有的叶子节点不也能遍历到所有的样本吗？
想像上课时老师点名，可以按学号从1到N来点名，也可以按小组序号来点名啊，那么这里的学号就对应这样本点标号，小组序号就对应着叶子节点标号。
有了上面的解释，相信大家应该能理解式（4.7）是怎么得来的了吧？😝

好啦，最困难的一步理解之后，就是喜闻乐见的求导了：
由于 $w_j^*=\mathop{\arg\min}_{w_j}\tilde{\mathcal{L}}^{\left(t\right)}$
可令 $\frac{\partial\tilde{\mathcal{L}}^{\left(t\right)}}{\partial w_j}=0$
得到每个叶结点 $j$ 的最优分数为
$w_j^*=-\frac{\sum_{i\in I_j}g_i}{\sum_{i\in I_j} h_i+\lambda} \tag{4.8}$

至此，我们求出了每个叶子节点的分数，但树长什么样还没有完全得出。

最开始有个假设被我们忽略了，就是叶子节点个数 $T$ ，这个条件一直被我们当做常量来计算了，也就是说我们已经人为设定好了树的形状就是 $T$ 个叶子的树，然后通过各种神操作得到了每个叶子的输出值。

那么接下来的问题是，哪些样本放在1号叶子里，哪些放2号叶子里呢？

我们回想下决策树是怎么做的？依照信息增益、信息增益比、gini系数对数据集进行分裂对吧？那我们是不是可以创造一种新的分裂方法呢？

代入每个叶结点 $j$ 的最优分数，得到最优化目标函数值：
$\tilde{\mathcal{L}}^{\left(t\right)}\left(q\right)=-\frac{1}{2}\sum_{j=1}^T \frac{\left(\sum_{i\in I_j} g_i\right)^2}{\sum_{i\in I_j} h_i+\lambda}+\gamma T \tag{4.9}$

我们发现这个目标函数仅与叶节点分数的表达式（4.8）非常像呀。是不是可以认为，叶节点的输出衡量了最终的模型性能？

ok，那我们假设 $I_L$ 和 $I_R$ 分别为分裂后左右结点的实例集，令 $I=I_L\cup I_R$ ，则分裂后损失减少量由下式得出：
$\mathcal{L}_{split}=\frac{1}{2}\left[\frac{\left(\sum_{i\in I_L} g_i\right)^2}{\sum_{i\in I_L}h_i+\lambda}+\frac{\left(\sum_{i\in I_R} g_i\right)^2}{\sum_{i\in I_R}h_i+\lambda}-\frac{\left(\sum_{i\in I} g_i\right)^2}{\sum_{i\in I}h_i+\lambda}\right]-\gamma \tag{4.10}$
式（4.10）表示的是数据集分裂前后，目标函数究竟提升了多少，这种操作是不是起到了跟gini系数一样的作用呢？显然答案是肯定的呀。

那我们就可以用（4.10）贪婪的遍历数据集的每一个特征，每一个特征下的数据类别：

算法4 分裂查找的精确贪婪算法

输入：当前结点实例集 $I$ ;特征维度 $d$
输出：根据最大分值分裂
（1） $gain\leftarrow 0$
（2） $G\leftarrow\sum_{i\in I}g_i$ ， $H\leftarrow\sum_{i\in I}h_i$
（3）for $k=1$ to $d$ do
（3.1） $G_L \leftarrow 0$ ， $H_L \leftarrow 0$
（3.2）for $j$ in sorted( $I$ , by $\mathbf{x}_{jk}$ ) do
（3.2.1） $G_L \leftarrow G_L+g_j$ ， $H_L \leftarrow H_L+h_j$
（3.2.2） $G_R \leftarrow G-G_L$ ， $H_R=H-H_L$
（3.2.3） $score \leftarrow \max\left(score,\frac{G_L^2}{H_L+\lambda}+\frac{G_R^2}{H_R+\lambda}-\frac{G^2}{H+\lambda}\right)$
（3.3）end
（4）end

至此，Xgboost打完收功😽

附送陈天奇老师的图片

xgboost的生成

参考
七月在线课程
《XGBoost A Scalable Tree Boosting System》-陈天奇
《统计学习方法》-李航

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