GBDT 回归的原理与 Python 实现
完整实现代码请参考本人的 p... 哦不是...github:gbdt_base.pygithub.com https://github.com/tushushu/imylu/blob/master/imylu/ensemble/gbdt_base.py
gbdt_regressor.pygithub.com https://github.com/tushushu/imylu/blob/master/imylu/ensemble/gbdt_regressor.py
gbdt_regressor_example.pygithub.com https://github.com/tushushu/imylu/blob/master/examples/gbdt_regressor_example.py
1. 原理篇
我们用人话而不是大段的数学公式来讲讲 GBDT 回归是怎么一回事.
1.1 温故知新
回归树是 GBDT 的基础, 之前的一篇文章曾经讲过回归树的原理和实现. 链接如下:
李小文: 回归树的原理与 Python 实现 zhuanlan.zhihu.com
1.2 预测年龄
仍然以预测同事年龄来举例, 从回归树那篇文章中我们可以知道, 如果需要通过一个常量来预测同事的年龄, 平均值是最佳的选择之一.
1.3 年龄的残差
我们不妨假设同事的年龄分别为 5 岁, 6 岁, 7 岁, 那么同事的平均年龄就是 6 岁. 所以我们用 6 岁这个常量来预测同事的年龄, 即[6, 6, 6]. 每个同事年龄的残差 = 年龄 - 预测值 = [5, 6, 7] - [6, 6, 6], 所以残差为[-1, 0, 1]
1.4 预测年龄的残差
为了让模型更加准确, 其中一个思路是让残差变小. 如何减少残差呢? 我们不妨对残差建立一颗回归树, 然后预测出准确的残差. 假设这棵树预测的残差是 [-0.9, 0, 0.9], 将上一轮的预测值和这一轮的预测值求和, 每个同事的年龄 = [6, 6, 6] + [-0.9, 0, 0.9] = [5.1, 6, 6.9], 显然与真实值[5, 6, 7] 更加接近了, 年龄的残差此时变为[-0.1, 0, 0.1]. 显然, 预测的准确性得到了提升.
1.5 GBDT
重新整理一下思路, 假设我们的预测一共迭代 3 轮 年龄:[5, 6, 7]
第 1 轮预测:[6, 6, 6] (平均值)
第 1 轮残差:[-1, 0, 1]
第 2 轮预测:[6, 6, 6] (平均值) + [-0.9, 0, 0.9] (第 1 颗回归树) = [5.1, 6, 6.9]
第 2 轮残差:[-0.1, 0, 0.1]
第 3 轮预测:[6, 6, 6] (平均值) + [-0.9, 0, 0.9] (第 1 颗回归树) + [-0.08, 0, 0.07] (第 2 颗回归树) = [5.02, 6, 6.97]
第 3 轮残差:[-0.08, 0, 0.03]
看上去残差越来越小, 而这种预测方式就是 GBDT 算法.
1.6 公式推导
看到这里, 相信您对 GBDT 已经有了直观的认识. 这么做有什么科学依据么, 为什么残差可以越来越小呢? 前方小段数学公式低能预警.
假设要做 m 轮预测, 预测函数为 Fm, 初始常量或每一轮的回归树为 fm, 输入变量为 X, 有:
设要预测的变量为 y, 采用 MSE 作为损失函数:
我们知道泰勒公式的一阶展开式是长成这个样子滴:
如果:
那么, 根据式 3 和式 4 可以得出:
根据式 2 可以知道, 损失函数的一阶偏导数为:
根据式 6 可以知道, 损失函数的二阶偏导数为:
蓄力结束, 开始放大招. 根据式 1, 损失函数的一阶导数为:
根据式 5, 将式 8 进一步展开为:
令式 9, 即损失函数的一阶偏导数为 0, 那么:
将式 6, 式 7 代入式 9 得到:
因此, 我们需要通过用第 m-1 轮残差的均值来得到函数 fm, 进而优化函数 Fm. 而回归树的原理就是通过最佳划分区域的均值来进行预测. 所以 fm 可以选用回归树作为基础模型, 将初始值, m-1 颗回归树的预测值相加便可以预测 y.
2. 实现篇
本人用全宇宙最简单的编程语言 --Python 实现了 GBDT 回归算法, 没有依赖任何第三方库, 便于学习和使用. 简单说明一下实现过程, 更详细的注释请参考本人 github 上的代码.
2.1 导入回归树类
回归树是我之前已经写好的一个类, 在之前的文章详细介绍过, 代码请参考:regression_tree.pygithub.com https://github.com/tushushu/Imylu/blob/master/regression_tree.py
from..tree.regression_treeimportRegressionTree
from..tree.regression_treeimportRegressionTree
2.2 创建 GradientBoostingBase 类
初始化, 存储回归树, 学习率, 初始预测值和变换函数.(注: 回归不需要做变换, 因此函数的返回值等于参数)classGradientBoostingBase(object):
def__init__(self):
self.trees=None
self.lr=None
self.init_val=None
self.fn=lambdax:x
classGradientBoostingBase(object):def__init__(self):self.trees=Noneself.lr=Noneself.init_val=Noneself.fn=lambdax:x
2.3 计算初始预测值
初始预测值即 y 的平均值.def_get_init_val(self,y):
returnsum(y)/len(y)
def_get_init_val(self,y):returnsum(y)/len(y)
2.4 计算残差def_get_residuals(self,y,y_hat):
return[yi-self.fn(y_hat_i)foryi,y_hat_iinzip(y,y_hat)]
def_get_residuals(self,y,y_hat):return[yi-self.fn(y_hat_i)foryi,y_hat_iinzip(y,y_hat)]
2.5 训练模型
训练模型的时候需要注意以下几点: 1. 控制树的最大深度 max_depth; 2. 控制分裂时最少的样本量 min_samples_split; 3. 训练每一棵回归树的时候要乘以一个学习率 lr, 防止模型过拟合; 4. 对样本进行抽样的时候要采用有放回的抽样方式.deffit(self,X,y,n_estimators,lr,max_depth,min_samples_split,subsample=None):
self.init_val=self._get_init_val(y)
n=len(y)
y_hat=[self.init_val]*n
residuals=self._get_residuals(y,y_hat)
self.trees=[]
self.lr=lr
for_inrange(n_estimators):
idx=range(n)
ifsubsampleisnotNone:
k=int(subsample*n)
idx=choices(population=idx,k=k)
X_sub=[X[i]foriinidx]
residuals_sub=[residuals[i]foriinidx]
y_hat_sub=[y_hat[i]foriinidx]
tree=RegressionTree()
tree.fit(X_sub,residuals_sub,max_depth,min_samples_split)
self._update_score(tree,X_sub,y_hat_sub,residuals_sub)
y_hat=[y_hat_i+lr*res_hat_ifory_hat_i,
res_hat_iinzip(y_hat,tree.predict(X))]
residuals=self._get_residuals(y,y_hat)
self.trees.append(tree)
deffit(self,X,y,n_estimators,lr,max_depth,min_samples_split,subsample=None):self.init_val=self._get_init_val(y)n=len(y)y_hat=[self.init_val]*n residuals=self._get_residuals(y,y_hat)self.trees=[]self.lr=lrfor_inrange(n_estimators):idx=range(n)ifsubsampleisnotNone:k=int(subsample*n)idx=choices(population=idx,k=k)X_sub=[X[i]foriinidx]residuals_sub=[residuals[i]foriinidx]y_hat_sub=[y_hat[i]foriinidx]tree=RegressionTree()tree.fit(X_sub,residuals_sub,max_depth,min_samples_split)self._update_score(tree,X_sub,y_hat_sub,residuals_sub)y_hat=[y_hat_i+lr*res_hat_ifory_hat_i,res_hat_iinzip(y_hat,tree.predict(X))]residuals=self._get_residuals(y,y_hat)self.trees.append(tree)
2.6 预测一个样本def_predict(self,Xi):
returnself.fn(self.init_val+sum(self.lr*tree._predict(Xi)fortreeinself.trees))
def_predict(self,Xi):returnself.fn(self.init_val+sum(self.lr*tree._predict(Xi)fortreeinself.trees))
2.7 预测多个样本defpredict(self,X):
return[self._predict(Xi)forXiinX]
defpredict(self,X):return[self._predict(Xi)forXiinX]
3 效果评估
3.1 main 函数
使用著名的波士顿房价数据集, 按照 7:3 的比例拆分为训练集和测试集, 训练模型, 并统计准确度.@run_time
defmain():
print("Tesing the accuracy of GBDT regressor...")
X,y=load_boston_house_prices()
X_train,X_test,y_train,y_test=train_test_split(
X,y,random_state=10)
reg=GradientBoostingRegressor()
reg.fit(X=X_train,y=y_train,n_estimators=4,
lr=0.5,max_depth=2,min_samples_split=2)
get_r2(reg,X_test,y_test)
@run_timedefmain():print("Tesing the accuracy of GBDT regressor...")X,y=load_boston_house_prices()X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=10)reg=GradientBoostingRegressor()reg.fit(X=X_train,y=y_train,n_estimators=4,lr=0.5,max_depth=2,min_samples_split=2)get_r2(reg,X_test,y_test)
3.2 效果展示
最终拟合优度 0.851, 运行时间 2.2 秒, 效果还算不错~
3.3 工具函数
本人自定义了一些工具函数, 可以在 github 上查看
utils.pygithub.com https://github.com/tushushu/imylu/blob/master/imylu/utils.py
1. run_time - 测试函数运行时间
2. load_boston_house_prices - 加载波士顿房价数据
3. train_test_split - 拆分训练集, 测试集
4. get_r2 - 计算拟合优度
总结
GBDT 回归的原理: 平均值加回归树
GBDT 回归的实现: 加加减减 for 循环
来源: http://blog.jobbole.com/114351/