完整实现代码请参考本人的 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_tree import RegressionTree
- from ..tree.regression_tree import RegressionTree
2.2 创建 GradientBoostingBase 类
初始化, 存储回归树, 学习率, 初始预测值和变换函数.(注: 回归不需要做变换, 因此函数的返回值等于参数)
- class GradientBoostingBase(object):
- def __init__(self):
- self.trees = None
- self.lr = None
- self.init_val = None
- self.fn = lambda x: x
- class GradientBoostingBase(object): def __init__(self): self.trees = None self.lr = None self.init_val = None self.fn = lambda x: x
2.3 计算初始预测值
初始预测值即 y 的平均值.
- def _get_init_val(self, y):
- return sum(y) / len(y)
- def _get_init_val(self, y): return sum(y) / len(y)
2.4 计算残差
- def _get_residuals(self, y, y_hat):
- return [yi - self.fn(y_hat_i) for yi, y_hat_i in zip(y, y_hat)]
- def _get_residuals(self, y, y_hat): return [yi - self.fn(y_hat_i) for yi, y_hat_i in zip(y, y_hat)]
2.5 训练模型
训练模型的时候需要注意以下几点: 1. 控制树的最大深度 max_depth; 2. 控制分裂时最少的样本量 min_samples_split; 3. 训练每一棵回归树的时候要乘以一个学习率 lr, 防止模型过拟合; 4. 对样本进行抽样的时候要采用有放回的抽样方式.
- def fit(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 _ in range(n_estimators):
- idx = range(n)
- if subsample is not None:
- k = int(subsample * n)
- idx = choices(population=idx, k=k)
- X_sub = [X[i] for i in idx]
- residuals_sub = [residuals[i] for i in idx]
- y_hat_sub = [y_hat[i] for i in idx]
- 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_i for y_hat_i,
- res_hat_i in zip(y_hat, tree.predict(X))]
- residuals = self._get_residuals(y, y_hat)
- self.trees.append(tree)
- def fit(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 _ in range(n_estimators): idx = range(n) if subsample is not None: k = int(subsample * n) idx = choices(population=idx, k=k) X_sub = [X[i] for i in idx] residuals_sub = [residuals[i] for i in idx] y_hat_sub = [y_hat[i] for i in idx] 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_i for y_hat_i, res_hat_i in zip(y_hat, tree.predict(X))] residuals = self._get_residuals(y, y_hat) self.trees.append(tree)
2.6 预测一个样本
- def _predict(self, Xi):
- return self.fn(self.init_val + sum(self.lr * tree._predict(Xi) for tree in self.trees))
- def _predict(self, Xi): return self.fn(self.init_val + sum(self.lr * tree._predict(Xi) for tree in self.trees))
2.7 预测多个样本
- def predict(self, X):
- return [self._predict(Xi) for Xi in X]
- def predict(self, X): return [self._predict(Xi) for Xi in X]
3 效果评估
3.1 main 函数
使用著名的波士顿房价数据集, 按照 7:3 的比例拆分为训练集和测试集, 训练模型, 并统计准确度.
- @run_time
- def main():
- 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_timedef main(): 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/