集成学习（GBDT）

参考了统计学习方法，西瓜书，Machine Learnig with python做的总结，还包含自己用sklearn做的一些对比实验，原文是写在jupyter上的，这里是直接转为.md导过来的，所以格式有些问题，有些东西还待完善…
参考博客
https://www.cnblogs.com/pinard/p/6140514.html，
https://blog.csdn.net/qq_22238533/article/details/79185969

建议先看 决策树,回归树,随机森林,Adaboosting,boosting tree

（三） 梯度提升(gradient boosting decision tree)

当损失函数是平方损失和指数函数损失函数，每一步优化很简单，但对于一般的损失函数而言往往不易这就引入了梯度提升，利用最速下降法的近似方法，关键是使用损失函数的负梯度在当前模型的值

梯度提升树算法（GBDT）-回归算法

输入：训练数据集$T=(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$;损失函数$L(y,f(x))$;
输出：回归树$\hat{f}(x)$
(1)初始化$$f_0(x)=argmi{n}_c\sum _{i=1}^{N}L(y_i,c)$$

(2)对$m=1,2,...,M$:
\quad\quad a.对$i=1,2,...,N$计算$$r_{mi}=-[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}{]}_{f(x)=f_{m-1}(x)}$$
\quad\quad b.对新数据集$(x_i,r_{mi})$拟合一个CART回归树，注意里标签值是$r_{mi}$，得到第$m$棵树的叶节点区域$R_{mj},j=1,2,..,J$
\quad\quad c.对$j=1,2,...,J$,计算$$c_{mj}=argmi{n}_c\sum _{x_i\in R_{mj}}L(y_i,f_{m-1}(x_i)+c)（CART回归树中是各个区域标签均值）$$
\quad\quad 得到回归树$$T_m=\sum _{m=1}^J{c}_{mj}I(x\in R_{mj})$$
\quad\quad d.更新$f_m(x)=f_{m-1}(x)+{\sum }_{j=1}^J{c}_{mj}I(x\in R_{mj})$
(3)得到回归树$$\hat{f}(x)=f_M(x)=\sum _{m=1}^M\sum _{j=1}^J{c}_{mj}I(x\in R_{mj})$$
注意:
1.算法第一步初始化，估计使损失函数极小化的常数值，它是只有根节点的树
2.第2(a)步中计算损失函数的负梯度在当前模型的值，将它作为残差估计，对于平方损失函数，它就是通常所说的残差，对于一般损失函数，它就是残差近似值
3.第2(b)步估计回归树叶节点区域，以拟合残差近似值
4.第2©步利用线性搜索估计叶节点区域的值，使损失函数极小化

例子

import pandas as pd
import numpy as np
data0={"x":[1,2,3,4,5,6,7,8,9,10],
     "y":[5.56,5.70,5.91,6.40,6.80,7.05,8.90,8.70,9.00,9.05]}#来自统计学习方法
data0=pd.DataFrame(data0)
data0

	x	y
0	1	5.56
1	2	5.70
2	3	5.91
3	4	6.40
4	5	6.80
5	6	7.05
6	7	8.90
7	8	8.70
8	9	9.00
9	10	9.05

现在我们来计算其GBRT模型,选取的最小化整体损失函数是$$\sum _{i=1}^N(y_i-f(x_i){)}^2$$,其单个样本损失函数为$$(y_i-f(x_i){)}^2$$梯度为$$-[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}]=2(y_i-f(x_i))$$这里直接使用$y_i-f(x_i)$

第一步先求$f_0(x)$即回归树$T_0(x)$，也就是求$f(x)=argmi{n}_{c_0}{\sum }_{i=0}^9(y_i-c_0{)}^2$,对此式求导可得$c_0=\frac{{\sum }_{i=0}^N{y}_i}{N}$

data_y=data0["y"]
data_L0=[]
c_0=np.sum(data0["y"])/len(data0)
print(c_0)
data_L0.append(np.sum((data_y-c_0)**2))
print(data_L0)

7.307
[19.114210000000003]

根据上面计算当我们选取$f_0(x)=7.05$，损失函数有最小值

第二步根据第2a步求出残差表

f_0=7.307
data1={"x":[1,2,3,4,5,6,7,8,9,10],
     "r":list(data0["y"]-f_0)}
data1=pd.DataFrame(data1)
data1

	x	r
0	1	-1.747
1	2	-1.607
2	3	-1.397
3	4	-0.907
4	5	-0.507
5	6	-0.257
6	7	1.593
7	8	1.393
8	9	1.693
9	10	1.743

第三步对残差进行拟合回归树(见CART回归树的建立):
首先计算搜索最优切分特征和最优切分点，由于特征就一个所以就是$x$，现在搜索最优切分点

data_L1=[]
C11=[]
C12=[]
# S=[(data1["x"][i]+data1["x"][i+1])/2 for i in  range(len(data1["x"])-1)]#备选切分点,这里用的是x间隔均值
# print(S)
# for s in S:
#     sum=0
#     c1=[]
#     c2=[]
#     for x,r in zip(data1["x"],data1["r"]):
#         if x<=s:
#             c1.append(r)
#         else:
#             c2.append(r)
#     c1=np.mean(c1)
#     c2=np.mean(c2)
#     for x,r in zip(data1["x"],data1["r"]):
#         if x<=s:
#             sum+=(r-c1)**2
#         else:
#             sum+=(r-c2)**2
#     data_L1.append(sum)

##直接选取x作为备选切分点
for s in data1["x"]:
    sum=0
    c11=[]
    c12=[]
    for x,r in zip(data1["x"],data1["r"]):
        if x<=s:
            c11.append(r)
        else:
            c12.append(r)
    c11=np.mean(c1)
    c12=np.mean(c2)
    C11.append(c1)
    C12.append(c2)
    for x,r in zip(data1["x"],data1["r"]):
        if x<=s:
            sum+=(r-c11)**2
        else:
            sum+=(r-c12)**2
    data_L1.append(sum)
print(data_L1)
print("c11:",C11[5])
print("c12:",C12[5])

[15.723088888888887, 12.0833875, 8.365638095238095, 5.775475000000003, 3.9113200000000017, 1.9300083333333335, 8.009809523809526, 11.735399999999997, 15.738599999999998, 19.11421]
c1: -1.0703333333333338
c2: 1.6054999999999997

可以得知最优切分点是$x=6$,个各区域输出值为$c_{11}=-1.0703,c_{12}=1.6055$,所以我们可以得到回归树$f_1(x_i)=7.307+{\sum }_{k=1}^2{c}_{1k}I(x_i\in R_{1k})$

第一轮迭代已经完成

进行第二轮迭代，首先还是计算新的需要拟合的数据

def F1(x):
    if x<=6:
        return 7.307-1.0703
    else:
        return 7.307+1.6055

求新残差表

data2={"x":[1,2,3,4,5,6,7,8,9,10],
     "r":[data0["y"][i]-F1(data0["x"][i]) for i in range(10)]}
data2=pd.DataFrame(data2)
data2

	x	r
0	1	-0.6767
1	2	-0.5367
2	3	-0.3267
3	4	0.1633
4	5	0.5633
5	6	0.8133
6	7	-0.0125
7	8	-0.2125
8	9	0.0875
9	10	0.1375

第三步对残差进行拟合回归树(见CART回归树的建立):
首先计算搜索最优切分特征和最优切分点，由于特征就一个所以就是?，现在搜索最优切分点

data_L2=[]
C21=[]
C22=[]
for s in data2["x"]:
    sum=0
    c21=[]
    c22=[]
    for x,r in zip(data2["x"],data2["r"]):
        if x<=s:
            c21.append(r)
        else:
            c22.append(r)
    c21=np.mean(c21)
    c22=np.mean(c22)
    C21.append(c21)
    C22.append(c22)
    for x,r in zip(data2["x"],data2["r"]):
        if x<=s:
            sum+=(r-c21)**2
        else:
            sum+=(r-c22)**2
    data_L2.append(sum)
print(data_L2)
print("c21:",C21[2])
print("c22:",C22[2])

[1.421235199999999, 1.0098567799999985, 0.8006163352380933, 1.1402758533333317, 1.665360511999999, 1.930008333333334, 1.9299332152380952, 1.89835646, 1.9089952799999999, 1.9300083360000004]
c21: -0.5133666666666675
c22: 0.2199857142857142

可以得知最优切分点是$x=3$,个各区域输出值为$c_{21}=-0.5134,c_{22}=0.22$,所以我们可以得到回归树$f_2(x_i)=7.307+{\sum }_{k=1}^2{c}_{1k}I(x_i\in R_{1}k)+{\sum }_{k=1}^2{c}_{2k}I(x_i\in R_{2k})$

进行第三轮迭代，首先还是计算新的需要拟合的数据

def F2(x):
    f1=F1(x)
    if x<=3:
        return f1-0.5134
    else:
        return f1+0.22
    

求新残差表

data3={"x":[1,2,3,4,5,6,7,8,9,10],
     "r":[data0["y"][i]-F2(data0["x"][i]) for i in range(10)]}
data3=pd.DataFrame(data3)
data3

	x	r
0	1	-0.1633
1	2	-0.0233
2	3	0.1867
3	4	-0.0567
4	5	0.3433
5	6	0.5933
6	7	-0.2325
7	8	-0.4325
8	9	-0.1325
9	10	-0.0825

第三步对残差进行拟合回归树(见CART回归树的建立):
首先计算搜索最优切分特征和最优切分点，由于特征就一个所以就是?，现在搜索最优切分点

data_L3=[]
C31=[]
C32=[]
for s in data3["x"]:
    sum=0
    c31=[]
    c32=[]
    for x,r in zip(data3["x"],data3["r"]):
        if x<=s:
            c31.append(r)
        else:
            c32.append(r)
    c31=np.mean(c31)
    c32=np.mean(c32)
    C31.append(c31)
    C32.append(c32)
    for x,r in zip(data3["x"],data3["r"]):
        if x<=s:
            sum+=(r-c31)**2
        else:
            sum+=(r-c32)**2
    data_L3.append(sum)
print(data_L3)
print("c31:",C31[5])
print("c32:",C32[5])

[0.7709864622222217, 0.7788541149999993, 0.8006163352380948, 0.7992815233333329, 0.767737584, 0.47794967333333394, 0.6009705066666668, 0.771725715, 0.7930538399999999, 0.8006163399999996]
c31: 0.14666666666666592
c32: -0.2200000000000002

可以得知最优切分点是$x=6$,个各区域输出值为$c_{31}=0.1467,c_{32}=-0.22$,所以我们可以得到回归树$f_3(x_i)=7.307+{\sum }_{k=1}^2{c}_{1k}I(x_i\in R_{1}k)+{\sum }_{k=1}^2{c}_{2k}I(x_i\in R_{2k})+{\sum }_{k=1}^2{c}_{3k}I(x_i\in R_{3k})$

这样继续迭代下去直至满足停止条件

梯度提升树算法-分类

GBDT的分类算法从思想上和GBDT的回归算法没有区别，但是由于样本标签不是连续的值，而是离散的类别，导致我们无法直接从输出类别去拟合类别输出的误差。为了解决这个问题，主要有两个方法，一个是用指数损失函数，此时GBDT退化为Adaboost算法。另一种方法是用类似于逻辑回归的对数似然损失函数的方法。也就是说，我们用的是类别的预测概率值和真实概率值的差来拟合损失。本文仅讨论用对数似然损失函数的GBDT分类。而对于对数似然损失函数，我们又有二元分类和多元分类的区别。

梯度提升树算法（GBDT）-二分类

对于二元分类GBDT,这里使用类似于逻辑回归的对数似然损失函数，则单样本损失函数为：$$L(y,f(x))=-logP(Y=y|x)=log(1+exp(-yf(x)))(1)$$

这个损失函数这么来的：
对于$x_i$,其标签是$y_i$,那么我们当然希望模型判断$x_i$为$y_i$的概率$P(Y=y_i|x_i)$越大越好，也就是希望$-P(Y=y_i|x_i)$越小越好。在二分类中 y ∈ { − 1 , + 1 } y\in \{-1,+1\} y∈{−1,+1}，$$P(Y=y|x)=\frac{exp(yf(x))}{1+exp(yf(x))}$$所以有：
$$P(Y=1|x)=\frac{exp(f(x))}{1+exp(f(x))}$$，$$P(Y=-1|x)=1-P(Y=1|x)=\frac{exp(-f(x))}{1+exp(-f(x))}=\frac{1}{1+exp(f(x))}$$.

所以可以得到$$L(y,f(x))=-logP(Y=y|x)=-log(\frac{exp(yf(x))}{1+exp(yf(x))})=log(\frac{1+exp(yf(x))}{exp(yf(x))})$$$$=log(1+exp(-yf(x)))$$这就是式（1）

则此时的负梯度误差为$$r_{mi}=-[\frac{\partial L(y,f(x_i))}{\partial f(x_i)}{]}_{f(x)=f_{m-1}(x)}=\frac{y_i}{1+exp(y_{i}f(x_i))}$$

对于生成的决策树，各个叶子节点的最佳负梯度拟合值为$$c_{mk}=argmi{n}_c\sum _{x_i\in R_{mk}}log(1+exp(-y_i(f_{t-1}(x_i)+c)))$$

上式难解，一般使用近似值代替$$c_{mk}=\frac{{\sum }_{x_i\in R_{mk}}{r}_{mi}}{{\sum }_{x_i\in R_{mi}}}|r_{mk}|(1-|r_{mk}|)$$

也就是除了负梯度计算和叶子结点的最佳负梯度拟合，二元GBDT分类同GBDT回归算法过程相同

注意

在上面推导中若 y ∈ { 1 , 0 } y\in \{1,0\} y∈{1,0}，则单样本损失函数可以写成逻辑回归的对数似然函数$$L(y,f(x))=ylog(P(Y=1|x))+(1-y)log(1-P(Y=1|x))$$

其中$$P(Y=1|x)=\frac{exp(f(x))}{1+exp(f(x))}$$，$$P(Y=0|x)=1-P(Y=1|x)=\frac{1}{1+exp(f(x))}$$.就不能写为$P(Y=y|x)=\frac{exp(yf(x))}{1+exp(yf(x))}$的合体形式

梯度提升树算法（GBDT）-多分类

注意：
1.对于多分类任务，GDBT的做法是采用一对多的策略,即，对每个类别训练$M$个分类器。假设有$K$个类别，那么训练完之后总共有$M\times K$颗树。
2.$K$个类别都拟合完第一颗树之后才开始拟合第二颗树，不允许先把某一个类别的M颗树学习完，再学习另外一个类别

对于多元分类，我们假设有$K$类，则此时的单样本对数似然损失函数函数$$L(y,f(x))=-\sum _{k=1}^K{y}_{k}log(P(Y=k|x))$$

整体的就是$$L(y,f(x))=-\sum _{i=1}^N{y}_{i}log(P(Y=y_i|x_i))$$

其中如若模型对样本$x$输出类别为$y_k$则$y_k=1$,否则$y_k=0$，也就是实际上只算了其标签类别对应的预测值得损失函数$$P(Y=k|x)=\frac{exp(f_k(x))}{{\sum }_{l=1}^{K}exp(f_l(x))}$$

也就是实际上模型输出的是一个K维向量，表示样例属于各个类别的相对可能性，然后通过上式进行归一化处理.

此时，在第$m$轮迭代，第$i$个样本对应类别$l$的负梯度误差是$$r_{mil}=-[\frac{\partial L(y,f(x_i))}{\partial f(x_i)}{]}_{f(x)=f_{m-1}l(x)}=y_{il}-P(Y=l|x)$$

对于各个类别生成的决策树，各个叶子节点的最佳负梯度拟合值为$$c_{mlk}=argmi{n}_c\sum _{x_i\in R_{m}lk}L(l,f_{m-1,l}(x)+c)$$$c_{mlk}$就表示$l$类别的第$m$轮迭代的决策树的第$k$个分区的拟合值.

由于上式难优化，所以用近似值代替$$c_{mlk}=\frac{K-1}{K}\frac{{\sum }_{x_i\in R_{mlk}}{r}_{mlk}}{{\sum }_{x_i\in R_{mlk}}|r_{mlk}|(1-|r_{mlk}|)}$$这里K表示类别数

那么就可以得到$l$类别第$m$轮迭代的决策树$$T_{ml}=\sum _{i=1}^N{c}_{mlk}I(x_i\in R_{mlk})$$

那么当所有类别的第$m$轮迭代分类器都生成了之后再继续下一轮迭代