Python实现梯度提升树(GDBT)

一 概念介绍

加法模型

下面公式代表的是一个加法模型
其中 b(x;γ_m）为基函数，γ_m 是基函数的参数，β_m 是基函数的系数。

前向分布算法

如果给定了训练数据集以及损失函数L(y,f(x)) ,学习加法模型的过程就是极小化损失函数的过程。即求解如下问题:

因为加法模型是由多个基模型组合而成的，很难直接去求解上述问题。 所以我们选择一个叫做前向分布算法的方法来求解。它的基本思想是: 如果能从前向后，每一步只学习一个基函数及其系数，逐步逼近目标函数式（上图中的式子),就可以简化优化的复杂度。

所以每步只需优化如下损失函数:

前向分布算法流程:

二 梯度提升树 (GDBT)

在前面的博客统计学习方法学习笔记-提升方法(一)之Adaboost（内含Python实现) 中学习了AdaBoost算法,它是利用前一轮迭代学习的弱学习器的误差率来更新训练集的权重，然后一轮轮迭代下去。 梯度提升树也是一种迭代学习的模型，它使用了前向学习算法，弱学习器限定为CART回归树 (在博客统计学习方法学习笔记-决策树(四)之Python实现CART算法中讲述)。

梯度提升树的学习

梯度提升树的学习过程使用前向分步算法 。初始时梯度提升树为f₀(x)=0,第m步模型为
其中f_m-1(x) 为当前模型，通过极小化损失函数确定下一棵梯度提升树的参数。
其中损失函数用的是平方损失函数
所以梯度提升树的损失函数可以表示为

其中 r = y-f_m-1(x) 我们称其为当前模型拟合数据的残差，我们在下一轮模型进行拟合的时候只需要拟合残差即可。

举一个通俗的例子
有天考试成绩出来了，小明只考了0分，回家妈妈打了他一顿，说三天以后我要重新让你考这张试卷，你必须得满分。小明就开始疯狂的学习这张试卷，第一天他学习前面30分的题，第二天他学习后40分的题，最后一天他学习最后30分的题，第四天，妈妈对他进行考试，满分。

其实我们的梯度提升树拟合模型就像小明学习试卷一样。初始时小明成绩为0，但是他的目标是100分，在第一轮迭代中他朝着100分的目标努力，只学习了30分，第二天他已经可以考30分了，他只用再拟合剩下的70分(这个70分就是经过第一轮学习的残差 目标值减去第一轮学习器学习到的值) 就可以了，他又拟合了40分，最后一天他只用再拟合最后30分(这个30分就是第二轮学习后的残差 目标值减去第一轮学习器加第二轮学习器学习到的值） 将三轮弱学习器组合成我们的加法模型。 我们可以发现加法模型的预测值为100，和我们的目标值是一样的，这时损失值为0， 这也是我们学习一个模型最希望看到的结果。

回归问题的梯度提升树算法

一个例子

我们拿一个预测房价的数据来讲解一下
(注意数据没有给完整，所以大家看看思路就可以了 不要尝试计算呀)

我们利用之前建立CART回归决策树的方法建立一个单层决策树，建立的规则是找出一个特征s，并且找一个切点将数据分成两部分。

比如说我们找出的特征是 每栋住宅的房间数，将这个特征的所有取值按大小进行排序 {6.421,6.575,6.998,7.147,7.185…} ，那它的切点可以这样取{(6.421+6.575)/2=6.498, (6.575+6.998)/2 =6.7865,(6.998+7.147)/2=7.0725 (7.147+7.185)/2=14.332…}

假设我们的切点是 6.498 那我们就可以把数据分成两个单元 房间数<6.498 与房间数>6.498 对于这两个单元来讲 c1 是第一个单元数据的房屋价格的均值 c2 是第二个单元数据的房屋价格的均值 根据这些我们遍历所有特征的所有切点使选取的特征以及切点可以满足以下公式
如果看了以上讲解仍然不是很清楚可以参考博客统计学习方法学习笔记-决策树(四)之Python实现CART算法

我们在第一轮迭代时创建出决策树

然后更新残差

生成新数据
根据拟合该数据生成第二棵梯度提升树
计算新残差

拟合该数据 生成第三棵决策树

给定了数据应该如何预测？

三 梯度提升树代码实现

模型代码

class Gx:
    def __init__(self,feature,split,left,right):
        #一个分布 这个分布是以哪个feature 作为划分的 划分点 split是 这个分布的系数是coef
         self.feature=feature
         self.split= split
         self.left =left#小于切分点的值为多少
         self.right = right #大于切分点的值为多少
class GDBT:
    def __init__(self,Train,Test,e):
        self.Train=Train #训练集
        self.Test=Test #测试集
        self.feature_num = self.Train.shape[1] - 2  # 求特征个数 因为训练数据的最后一列是标签 所以特征数等于维度数减一
        self.label = self.Train.shape[1] - 2  # 标签所在列
        self.r = self.Train.shape[1] - 1  # 残差所在列
        self.Gx=[] #基分类器
        #self.M = M #迭代次数
        self.e=e # 阈值 在训练数据的损失函数值小于e时停止训练
        self.split = self.cal_split()

    def cal_split(self):
        # 因为数据的特征都是连续值 计算每一个特征的切分点 每两个数据的中点集合作为这个特征的切点集合
        split = {}
        for i in range(self.feature_num):
            split[i] = []  # 将某一个特征的所有切分点放到一个列表中
            d = self.Train[np.argsort(self.Train[:, i])]  # 以这个特征的大小来进行排序
            for j in range(d.shape[0] - 1):
                sp = (d[j][i] + d[j + 1][i]) / 2
                if sp not in split[i]:  # 可能有的切分点是重复的 这一步去重
                    split[i].append(sp)
        return split

    def cost(self, data): #求损失值
        X = data[:, self.r:self.r + 1]
        avg = np.sum(X) / X.shape[0]  # 平均值
        sum = 0
        for i in range(X.shape[0]):
            sum += (X[i][0] - avg) ** 2
        return sum,avg

    def update_r(self,G): #更新残差值
        feature=G.feature
        split=G.split
        left=G.left
        right=G.right
        for i in range(self.Train.shape[0]):
            if self.Train[i][feature]  <= split:
                self.Train[i][self.r]=self.Train[i][self.r]-left
            else:
                self.Train[i][self.r] = self.Train[i][self.r] - right

    def cal_best_split(self):
        min_cost=sys.maxsize
        min_split=0
        min_feature=0
        min_left=0
        min_right=0
        #选择最佳特征，以及最佳切分点
        for i in range(self.feature_num):  # 遍历所有特征取值
             for j in self.split[i]:#遍历每一个切分点
                 left = self.Train[(self.Train[:, i] <= j), :]
                 right = self.Train[(self.Train[:, i] > j), :]
                 cost1,avg1=self.cost(left)
                 cost2,avg2=self.cost(right)
                 # print(cost1)
                 # print(cost2)
                 # print("******")
                 if cost1+cost2 < min_cost:
                     min_cost= cost1+cost2
                     min_feature=i
                     min_split=j
                     min_left=avg1
                     min_right=avg2
        G=Gx(min_feature,min_split,min_left,min_right)
        self.update_r(G)
        self.Gx.append(G)

    def fit(self):
        L=sys.maxsize
        while L>self.e:
            self.cal_best_split()
            print("当前损失值:")
            L=self.cal_cost(self.Train)
            print(L)
            print("***************")

    def print_Gx(self):
        for gx in self.Gx:
            print("特征:")
            print(gx.feature)
            print("切分点")
            print(gx.split)
            print("左值")
            print(gx.left)
            print("右值")
            print(gx.right)
            print("******************************")

    def cal_cost(self,data):
        L=0#预测损失值
        for i in range(data.shape[0]):
            sum=0# 预测值
            for gx in self.Gx:
                feature=gx.feature
                split=gx.split
                left=gx.left
                right=gx.right
                if data[i][feature] <=split:
                    sum+=left
                else:
                    sum+=right
            L+=(data[i][self.label]-sum)**2
        return L

调用模型

boston=datasets.load_boston()
data1=boston.data
print(data1.shape)
data2=boston.target
data2.resize([data2.shape[0],1])
data3=np.concatenate((data1,data2), axis=1)
data=np.concatenate((data3,data2), axis=1)
print(data.shape)
print(data)
train_size = int(len(data) * 0.7)#划分训练集与测试集
train = np.array(data[:train_size])
# # print(X_train.shape)
test=np.array(data[train_size:])
model=GDBT(train,test,100)
model.fit()