GBDT原理分析以及XGBoost代码实践

简介

GBDT中文译为梯度提升决策树。GBDT是以分类树或者回归树作为基本分类器的提升方法，它被认为是统计学习中性能最好的方法之一。在深度学习兴起和流行之前，GBDT是公认效果最出色的几个模型之一。虽然现在已经号称进入了深度学习以及人工智能时代，但是GBDT也没有落伍，它依然在很多的场景和公司当中被广泛使用。

GBDT基础概念

GBDT的英文原文是Gradient Boosting Decision Tree，即梯度提升决策树。从它的英文表述我们可以看出来，GBDT的基础还是决策树。关于决策树，读者可以参考其他文章，本文不再赘述。在GBDT当中用到的主要是决策树的CART算法，在CART算法当中，我们每次都会选择一个特征并且寻找一个阈值进行二分。将样本根据阈值分成小于等于阈值的以及大于阈值的两个部分，在CART树当中，同一个特征可以重复使用，其他类似的ID3和C4.5都没有这个性质。

另外一个关键词是Boosting，Boosting表示一种集成模型的训练方法，具体可搜索AdaBoost模型相关内容。它最大的特点就是会训练多个模型，通过不断地迭代来降低整体模型的偏差。比如在Adaboost模型当中，会设置多个弱分类器，根据这些分类器的表现我们会给与它们不同的权值。通过这种设计尽可能让效果好的分类器拥有高权重，从而保证模型的拟合能力。但GBDT的Boosting方法与众不同，它是一个由多棵CART决策回归树构成的加法模型。我们可以简单理解成最后整个模型的预测结果是所有回归树预测结果的和，理解了这一点对于后面理解梯度和残差非常重要。

GBDT概述

GBDT的思想可以用一个通俗的例子解释，假如有个人30岁，我们首先用20岁去拟合，发现损失有10岁，这时我们用6岁去拟合剩下的损失，发现差距还有4岁，第三轮我们用3岁拟合剩下的差距，差距就只有一岁了。如果我们的迭代轮数还没有完，可以继续迭代下面，每一轮迭代，拟合的岁数误差都会减小。
对应到GBDT算法里面，描述如下，我们每一轮都会训练一棵决策树，每棵树都会给出对应的预测值，我们将所有决策树的预测值加起来就是我们最终的预测值 $y$ ，上述这个小例子可以用以下图例来表示：

注意第一轮的决策树只有一个叶子节点，文章后续会详细解释。

GBDT算法

下面给出GDBT算法的具体步骤：
（1）输入数据为训练数据集 $T=\{ (x_1,y_1),(x_2,y_2),...,(x_n.,y_n)\}$ ，损失函数为 $L(y_i, F(x))$ 。
（2）使用一个常量来初始化模型：

（3）从

m=1

到

M

（M是决策树的个数）:
　　(A) 对于

i=1,2,...,n

，计算负梯度：

　　(B) 将上步得到的残差作为样本新的真实值，并将数据 $(x_i,r_{im})， i=1,2,..N$ 作为下一棵树的训练数据，得到一棵新的回归树 $f_m(x)$ ，其对应的叶子节点区域为 $R_{jm},j=1,2,..,J$ ，其中 $J$ 为回归树的叶子节点个数。
　　(C)对于叶子区域

j=1,2,...,J_m

，计算：

　　(D) 更新

f_m(x)

：

（4）得到最终强学习器：

GBDT实例讲解

上述算法看起来比较晦涩难懂，接下来举一个具体的例子，并且按照上述算法步骤一步一步地展示计算过程。

第一步：
（1）输入数据为训练数据集 $T=\{ (x_1,y_1),(x_2,y_2),...,(x_n.,y_n)\}$ ，损失函数为 $L(y_i, F(x))$ 。

为了方便演示，我们选取一个极其简单的训练数据集，它只包含3条数据。每条数据包含3个特征，分别是“身高”、“喜欢的颜色”、“性别”，标签值为“体重”。训练数据集如下：

训练数据集

接着是选取损失函数，在回归问题中，最常使用的是MSE损失函数（当然也可以选取其他的损失函数），定义如下：

MSE损失函数

其实上式前面的1/2不是必须的，这里加上的原因是为了在计算负梯度的时候，刚好和指数部分的2抵消掉。

其对应的负梯度为：

第二步：
（2）使用一个常量来初始化模型：

我们首先需要找到一个常量来初始化模型，即上式中的 $c$ 。我们将训练数据集中的数据带入上式得：
$f_0(x) = \frac {1}{2}(88-c)^2+ \frac {1}{2}(76-c)^2+ \frac {1}{2}(56-c)^2$ 现在的问题是我们如何求出上式的最小值呢？可以使用梯度下降法，用下面的图来展示：

左边的3个绿点自上而下分别代表88、76、56，我们的目标是要找到一个

c

，使得它与这些数据的残差的平方和最小。右边的曲线代表的是

c

和残差和的关系图。可以看到谷底在黑色箭头所指的地方。计算方式也很简单，首先要计算上式的导数，然后设为0，即可求出。计算过程如下，首先计算

f_0(x)

的导数：

f_0'(x) = (88-c)+ (76-c)+ (56-c)

然后令其为0，得：

f_0'(x) = (88-c)+ (76-c)+ (56-c) = 0

求出

c = \frac {88+76+56}{3} = 73.3

，其实就是平均值啦。

这里要注意一点，很多文章中是直接说 $f_0(x)$ 就等于训练数据集标签的平均值，实际上这只有在损失函数使用的是MSE的情况下才成立。我们要知道，其实这是通过梯度下降法，令导数等于0求出的，只是这里很巧，刚好就等于平均值。

okay，至此，我们其实已经计算出来了第一棵回归树，它就只有一个叶子节点，值为73.3。即这棵回归树把所有的样本的体重都预测为73.3。看起来很蠢不是嘛，但是别急，我们还会继续生成其他的回归树来减小误差。

第三步：
（3）从 $m=1$ 到 $M$ （M是决策树的个数）:
　　(A) 对于 $i=1,2,...,n$ ，计算：

这一步是个循环，即对GBDT中的 $M$ 个决策树都执行相同的过程。步骤（A）是要计算每一个样本数据的负梯度，在第一步中我们其实已经计算好了，当前 $m=1$ ，那么 $f_{m-1}(x) = f_0(x)$ ，又 $f_0(x) = 73.3$ ，故负梯度为 $(y - 73.3)$ 。然后我们把所有样本的标签值带入这个表达式中，可以得到所有样本的残差 $r_{im}$ 如下：

残差

以第一行为例，

r_{i,1} = 88 - 73.3 = 14.7

，其他以此类推。

第三步：
　　(B) 将上步得到的残差作为样本新的真实值，并将数据 $(x_i,r_{im})， i=1,2,..N$ 作为下一棵树的训练数据，得到一棵新的回归树 $f_m(x)$ ，其对应的叶子节点区域为 $R_{jm},j=1,2,..,J$ ，其中 $J$ 为回归树的叶子节点个数。

这一步很重要，在这一步中我们构造的决策树是要预测上一步中计算出来的参数 $r_{im}$ ，而不再是体重数据了。见下图，最后一列才是我们要预测的：

我们以“身高<1.55”为划分节点，拟合出来一棵回归树如下：

R_{1,1}、R_{2,1}

分别是对应的叶子节点，残差为对应叶子节点的值。
此时有同学可能会问了，为什么是“身高小于1.55”为划分节点，用其他特征不行么？这里解释一下，实际上我们可以根据“身高”、“喜欢的颜色”、“性别”划分出若干棵决策树，但是究竟使用哪一棵呢？我们选择均方损失最小的哪一棵。比如，针对上述的问题，我们其实可以还构造以下2棵决策树，如下：

以“喜爱颜色”和“性别”进行划分得到的两棵决策树

以这三种特征划分的3棵决策树对应的均方损失如下：

可以看到，使用“身高”特征划分的决策树的损失最小，因此我们选用以“身高”特征进行划分的决策树。

第三步：
　　(C)对于叶子区域 $j=1,2,...,J_m$ ，计算：

根据上一步产生的回归树，我们可以按照上式来计算 $\Upsilon_{jm}$ 。记得此时 $m=1，j = \{1,2\}$ 。则有：
$\Upsilon_{11} = argmin( \frac {1}{2}(56 - (77.3 + \Upsilon))^2 = argmin \frac {1}{2}(- 17.3 - \Upsilon))^2)$ 依然利用梯度下降法，求导令为0，得到 $\Upsilon_{11} =$ - $17.3$ 。同理有：
$\Upsilon_{21} = argmin( \ \frac {1}{2}(14.7 - \Upsilon))^2 + \frac {1}{2}(2.7 - \Upsilon))^2 )$ 利用梯度下降法，求导令为0，得到 $\Upsilon_{21} = 8.7$ 。

其实，由于我们在第一步已经证明了，对于MSE损失，最小值就是平均值，因此我们可以很快地计算出每个叶子节点对应的标签 $\Upsilon$ 。

回归树更新如下：

其实这里可以这么理解，由于我们每一棵决策树的叶子节点其实都是需要输出一个label的，但是有些叶子节点包含了不止一个元素，因此我们需要找到一个令均方损失最小的值来作为这个叶子节点的整体标签。

第三步：
　　(D) 更新 $f_m(x)$ ：

根据前面几步，我们得到了 $f_0(x)=77.3$ ，以及一棵决策树，那么更新就很简单了。

其中上图中的0.1代表的是学习率，它的作用是给我们拟合的每一棵回归树乘上一个介于0~1之间的参数，可以通过增大回归树的个数来避免GBDT陷入过拟合。接下来我们就可以用新的模型来对训练数据进行预测。
对第一条数据的预测：

对第二条数据的预测：

对第三条数据的预测：

对训练数据集进行新一轮的预测值之后，我们可以得到新一轮的残差，分别是

88-74.2=13.8

，

76-74.2=1.8

，

56-71.6=-15.6

。我们重复第三步的操作，可以得到新的一棵回归树如下：

由此，我们就可以继续更新强学习器，得到：

假设，我们现在有一个测试数据，我们就可以使用此前训练好的GBDT来进行体重预测了，测试数据如下：

测试数据

预测结果为：

到这里，我们就完整地演示了一遍GBDT算法的执行过程，并且产生了一棵梯度提升决策树用于预测体重。

代码实践

XGBoost

XGBoost是陈天奇等人开发的一个开源机器学习项目，高效地实现了GBDT算法并进行了算法和工程上的许多改进，被广泛应用在Kaggle竞赛及其他许多机器学习竞赛中并取得了不错的成绩。

房价预测

本文使用XGBoost库来进行房价预测，问题来自于Kaggle的一个比赛项目房价预测，给出房子的众多特征，要求建立数值回归模型，预测房子的价格。
数据集可到此处下载。
首先加载数据：

import pandas as pd

trainDF = pd.read_csv("train.csv")
testDF = pd.read_csv("test.csv")

先看一下整体的价格分布：

import seaborn as sns
import matplotlib.pyplot as plt

sns.distplot(trainDF['SalePrice'])

# 因为seaborn是matplotlib的高级库，所以可以用plt.show()来调动绘图
plt.show()

价格分布图

查看相关系数并绘制热力图：

corrmat = trainDF.corr()
# 得到saleprice相对系数前十大的数据 
cols = corrmat.nlargest(10,'SalePrice').index
largest_price = trainDF[cols].corr()

# 绘制这前十的相关系数的热点图,其中annot=True表示将数值写入格子
sns.heatmap(largest_price,annot=True, xticklabels=largest_price.columns,yticklabels=largest_price.index)
plt.show()

print(corrmat.nlargest(10,'SalePrice'))

热点图

先查看所有特征相对于价格的散点图，查看哪些是异常值，然后去除异常值：

for var in trainDF.select_dtypes(include=[np.number]).columns:
    plt.scatter(trainDF[var],trainDF['SalePrice'])
    plt.xlabel(var)
    plt.ylabel('SalePrice')
    plt.show()

其中部分异常值筛选如下：

汇总起来，以下特征的取值范围中price的值异常，需要去除该数据：

trainDF = trainDF.drop(trainDF[(trainDF['GrLivArea']>4000) & (trainDF['SalePrice']<300000)].index) # 默认axis=0也就是删除行
trainDF = trainDF.drop(trainDF[(trainDF['LotFrontage']>300) & (trainDF['SalePrice']<300000)].index)
trainDF = trainDF.drop(trainDF[(trainDF['BsmtFinSF1']>5000) & (trainDF['SalePrice']<200000)].index)
trainDF = trainDF.drop(trainDF[(trainDF['TotalBsmtSF']>6000) & (trainDF['SalePrice']<200000)].index)
trainDF = trainDF.drop(trainDF[(trainDF['1stFlrSF']>4000) & (trainDF['SalePrice']<200000)].index)

将测试集和训练集数据进行合并，然后看一下哪些值是缺失的，再删除缺失值超过80%的特征:

allDF = pd.concat([trainDF.drop(["SalePrice"],axis=1),testDF],axis=0)
ratio = (allDF.isnull().sum()/len(allDF)).sort_values(ascending=False)
print(ratio)
allDF.drop(ratio[ratio>0.8].index,axis=1)

对未删除的特征，进行缺失值处理：

allDF["LotFrontage"] = allDF.groupby("Neighborhood")["LotFrontage"].transform(
lambda x: x.fillna(x.median()))

for col in ('FireplaceQu',\
    'GarageType', 'GarageFinish', 'GarageQual', 'GarageCond', \
    'BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2',\
    'MasVnrType',\
    'MSSubClass'):
    allDF[col] = allDF[col].fillna('None')

for col in ('GarageYrBlt', 'GarageArea', 'GarageCars',\
    'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath',\
    'MasVnrArea'):
    allDF[col] = allDF[col].fillna(0)

# 因为数值都一样，没有存在意义
allDF = allDF.drop(['Utilities'], axis=1)

for col in ('MSZoning','Functional','Electrical','KitchenQual',\
    'Exterior1st','Exterior2nd','SaleType'):
    allDF[col] = allDF[col].fillna(allDF[col].mode()[0])

有些数值看起来是连续值，但是上它们的取值只是在有限数值上变动，因此可以变成离散值。方法就是将数值类型转换为字符串类型。
但是注意，最终我们需要将字符串通过哑变量转换为数值类型，以便给XGBoost计算:

allDF['MSSubClass'] = allDF['MSSubClass'].apply(str)   # 应用到每一个元素
allDF['OverallCond'] = allDF['OverallCond'].astype(str)
allDF['YrSold'] = allDF['YrSold'].astype(str)
allDF['MoSold'] = allDF['MoSold'].astype(str)

创造新特征：

allDF['TotalSF'] = allDF['TotalBsmtSF'] + allDF['1stFlrSF'] + allDF['2ndFlrSF']

转换为哑变量:

dm_allDF = pd.get_dummies(allDF)

使用XGBoost训练和预测：

from xgboost import XGBRegressor
from sklearn.model_selection import train_test_split

# 重新分割训练和测试数据
dm_trainDF = dm_allDF[:len(trainDF)]
dm_testDF = dm_allDF[len(trainDF):]

# 去掉id号
train_data = dm_trainDF.drop(['Id'],axis=1).values
train_label = trainDF["SalePrice"]

X_test_ids = dm_testDF['Id'].values
X_test = dm_testDF.drop(['Id'],axis=1).values

# 分割训练集和验证集
X_train, X_valid, Y_train, Y_valid = train_test_split(train_data, train_label, test_size=0.2)

xgb = XGBRegressor(n_estimators=500, learning_rate=0.05, min_child_weight=5, max_depth=4)
xgb.fit(X_train,Y_train)
print("Validation:", xgb.score(X_valid,Y_valid))

predict = xgb.predict(X_test)
print("predict:", predict)

输出结果：

测试结果

参考

最后编辑于：2021.01.14 12:41:01

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