顺序集成之Boosting--Gradient Boosting梯度提升

  和上一期讲的AdaBoost一样，Gradient Boosting 也是通过按顺序添加predictors 到集成中来工作。但是，并不像AdaBoost那样，在每次iteration的时候调整样本的权重，Gradient Boosting这个方法是使用新的predictor去拟合旧的predictor产生的的残差。也就是说在残差的基础上进行拟合，拟合完成后剩下的残差又可以用新的predictor来拟合，步骤如下：

  第一步：使用DecisionTreeRegressor 来拟合训练集;
  第二步：对于第一个predictor产生的残差用第二个DecisionTreeRegressor 来训练；
  第三步：在第二个predictor产生的残差上面训练第三个regressor；
  …
  最后一步：们就会有一个包含n棵树的集合，它通过把所有的树的预测结果相加，从而对新的样本进行预测。

使用代码一步步实现上面的步骤

from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import numpy as np 
import pandas as pd
from sklearn.ensemble import AdaBoostClassifier
from sklearn.tree import DecisionTreeClassifier

  产生一些数据样本点：

np.random.seed(42)
X=np.random.rand(100,1)-0.5
y=3*X[:,0]**2+0.05*np.random.randn(100)

  数据点大致长下面这样子：

plt.scatter(X,y,c=y)

  使用单棵最大深度为2的决策树：

from sklearn.tree import DecisionTreeRegressor
tree_reg1=DecisionTreeRegressor(max_depth=2)
y1=tree_reg1.fit(X,y).predict(X)
resid1=y-y1

resid1_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid1).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid1_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
#ax1.title('h(x)=h0(x)+h1(x)')
ax2.scatter(X,resid1,label='residual1')
ax2.plot(df2.X,df2.y,label='h1(x)')
ax2.legend()
#ax2.title('h1(x)')

  使用决策树，设置max_depth=2，然后拟合预测得到单棵决策树的预测值y1,并且使用y-y1得到resid1，得到单棵决策树拟合一次后的残差，然后使用剩下的残差再进行拟合，拟合得到的结果用resid1_pred来表示，做出如图所示的结果。

  下图中左图散点为原始的数据集，线段为为一棵决策树拟合的结果，右图散点为拟合一次残差的结果，线段为残差用单棵决策树拟合的结果。

  然后我们再用最原始的y减去第一次预测得到的y1，再减去第一次残差拟合得到的预测值resid1_pred得到第二个残差resid2,然后再用第二个残差训练单棵决策树再进行预测，得到第二个残差的预测值resid2_pred。

resid2=y-y1-resid1_pred
resid2_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid2).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred+resid2_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid2_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)+h2(x)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(X,resid2,label='residual2')
ax2.plot(df2.X,df2.y,label='h2(x)')
ax2.legend()

  下图中左图线段为第一次拟合和拟合两次残差相加的结果，右图散点为第二次残差，线段为第二次残差的预测值。

  第三次：

resid3=y-y1-resid1_pred-resid2_pred
resid3_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid3).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred+resid2_pred+resid3_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid3_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))

ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)+h2(x)+h3(x)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(X,resid3,label='residual3')
ax2.plot(df2.X,df2.y,label='h3(x)')
ax2.legend()

  第四次：

resid4=y-y1-resid1_pred-resid2_pred-resid3_pred
resid4_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid4).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred+resid2_pred+resid3_pred+resid4_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid4_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)+h2(x)+h3(x)+h4(x)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(X,resid4,label='residual4')
ax2.plot(df2.X,df2.y,label='h4(x)')
ax2.legend()

  第五次：

resid5=y-y1-resid1_pred-resid2_pred-resid3_pred-resid4_pred
resid5_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid5).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred+resid2_pred+resid3_pred+resid4_pred+resid5_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid5_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)+h2(x)+h3(x)+h4(x)+h5(x)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(X,resid5,label='residual5')
ax2.plot(df2.X,df2.y,label='h5(x)')
ax2.legend()

  第六次：

resid6=y-y1-resid1_pred-resid2_pred-resid3_pred-resid4_pred-resid5_pred
resid6_pred=DecisionTreeRegressor(max_depth=2).fit(X,resid6).predict(X)

df1=pd.DataFrame({"X":X[:,0],"y":y1+resid1_pred+resid2_pred+resid3_pred+resid4_pred+resid5_pred+resid6_pred})
df1=df1.sort_values(by=['X'])
df2=pd.DataFrame({"X":X[:,0],"y":resid6_pred})
df2=df2.sort_values(by=['X'])

fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df1.X,df1.y,label='h(x)=h0(x)+h1(x)+h2(x)+h3(x)+h4(x)+h5(x)+h(6)')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(X,resid5,label='residual6')
ax2.plot(df2.X,df2.y,label='h6(x)')
ax2.legend()

  以此类推，可以发现最后剩下的残差越来越小，曲线也越来越接近于y=3*x^2。上面的分步骤也可以通过下面的代码来一次性查看：

通过全局变量实现重复部分累加

np.random.seed(42)
X=np.random.rand(100,1)-0.5
y=3*X[:,0]**2+0.05*np.random.randn(100)
plt.scatter(X,y,c=y)

from sklearn.tree import DecisionTreeRegressor
y_sum=0
y_resid=y
y1=DecisionTreeRegressor(max_depth=2).fit(X,y).predict(X)
y_resid=y_resid-y1
y_sum=y_sum+y1

def calculate(X,dy):
    global y_resid, y_sum
    dy=DecisionTreeRegressor(max_depth=2).fit(X,dy).predict(X)
    y_resid=y_resid-dy
   # y_resid_p=DecisionTreeRegressor(max_depth=2).fit(X,y_resid).predict(X)
    y_sum=y_sum+dy
    df1=pd.DataFrame({"X":X[:,0],"y":y_sum}).sort_values(by=['X'])
    df2=pd.DataFrame({"X":X[:,0],"y":dy}).sort_values(by=["X"])
    
    fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
    ax1.plot(df1.X,df1.y)
    ax1.scatter(X,y,c=y)
    
    ax2.scatter(X,y_resid)
    ax2.plot(df2.X,df2.y)

for i in range(25):
    calculate(X,y_resid)

GrandientBoostingRegressor

  下面使用GrandientBoostingRegressor来拟合曲线：

from sklearn.ensemble import GradientBoostingRegressor

import matplotlib.patches as mpatches
gbrt=GradientBoostingRegressor(max_depth=2,n_estimators=3,learning_rate=1)
gbrt.fit(X,y)
df=pd.DataFrame({"X":X[:,0],"y":gbrt.predict(X)})
df=df.sort_values(by=['X'])
plt.scatter(X,y,c=y)
line1,=plt.plot(df.X,df.y,c='red',label='Ensemble predictions')
plt.legend(handles=[line1],loc='upper center')

  可以看到这个图和h(x)=h0(x)+h1(x)+h2(x)这个图像是一样的，即和使用单棵决策树拟合三次的结果是一样的。因为图片尺寸问题可能不是很明显，查看具体的数据可以看到，其y值是一模一样的。

  上面是通过GradientBoostingRegressor使用3棵决策树的结果，下面看看使用200棵决策树的结果：

gbrt=GradientBoostingRegressor(max_depth=2,n_estimators=200,learning_rate=0.1)
gbrt.fit(X,y)
df=pd.DataFrame({"X":X[:,0],"y":gbrt.predict(X)})
df=df.sort_values(by=['X'])
plt.scatter(X,y,c=y)
line1,=plt.plot(df.X,df.y,c='red',label='Ensemble predictions')
plt.legend(handles=[line1],loc='upper center')

resid=y-gbrt.predict(X)

df2=pd.DataFrame({"X":X[:,0],"y":resid})
df2=df2.sort_values(by=['X'])
fig,[ax1,ax2]=plt.subplots(1,2,figsize=(12,6))
ax1.plot(df.X,df.y,c='red',label='Ensemble predictions')
ax1.scatter(X,y,c=y,label='Training set')
ax1.legend()
ax2.scatter(df2.X,df2.y,label='residual200')
ax2.plot(df2.X,gbrt.fit(X,df2.y).predict(X),label='h200(x)')

  可以发现左图中的红色拟合线越来越接近于y=3*x^2并且残差已经减少到正负0.04左右。

简单说说GradientBoostingRegressor

Sklearn.ensemble.GradientBoostingRegressor(*,loss=’ls’,learning_rate=0.1,n_estimator=100,subsample=1.0,criterion=’fridman_mse’,min_samples_split=2,min_samples_leaf=1,min_weight_fraction_leaf=0.0,warm_start=False,…etc.)

  好多参数其实在单棵决策树里面有的。

Loss:{‘ls’,’lad’,’huber’,’quantile’},default=’ls’

  Ls: least squares regression，最小二乘回归，也就是通常说的L2-norm,L2正则，通过最小化下式S来达到优化的目的，

为目标值，

为估计值：

  Lad:least absolute deviation，最小绝对偏差，通常说的L1-norm，L1正则，通过最小化下式S来达到优化的目的，

为目标值，

为估计值：

；
  Huber:结合了ls和lad；与平方误差相比，Huber损失对数据离群值的敏感度低。当误差较小时，它将变成二次方，那具体有多小，取决于可调整的超参数
（delta)。Huber的公式：

当

小于

，更接近与MSE(mean square error)：

  Quantile:分位数，取值范围（0，1）。基于分位数的回归旨在给定预测变量的某些值得情况下估计变量的条件“分位数”。之所以选择quantile，是因为它对于heteroscedastic data（异方差数据）有比较好的惩罚机制，从而使得模型得到更好的拟合结果。从以下公式可以发现，通过设定

值的大小来调整预测比实际值偏大或者偏小的惩罚程度。分位数损失实际上只是MAE的扩展，当quantile(

)为0.5时，下式就和MAE很接近了。又比如，当

=0.25时，那么对于预测值大于目标值的预测结果惩罚更大。


Learning_rate:float, default=0.1;
学习速率，这个参数和单棵决策树学习速率是的一样。

N_estimators:boosting阶段的数据量，即决策树的数量。

Subsample:float，default=1.0。用来拟合单棵决策树的样本比例，即子样本比例。

Criterion：{‘friedman_mse’,‘mse’,‘mae’}，default = ‘friedman_mse’
  The function to measure the quality of a split. Supported criteria are ‘mse’ for the mean squared error, which is equal to variance reduction as feature selection criterion, and ‘mae’ for the mean absolute error.

  普通的MSE：

is the predicted values;

  The least-squares improvement criterion used to evaluate potential splits of a currently terminal region R into two subregions
is:

  Friedman Criterion:

Where

are the left and right daughter response means respectively, and

are the corresponding sums of the weights.

参考文献：https://statweb.stanford.edu/~jhf/ftp/stobst.pdf

Warm_start:bool, default=False.
  When set to True, reuse the solution of the previous call to fit and add more estimators to the ensemble, otherwise, just erase the previous solution.

  由上文知道，在Gradient Boosting中，后面的决策树的结果是在前面的决策树的残差上进行拟合的，一步步对残差拟合得到。那么当warm_start这个参数设置为True的时候，使用之前的模型拟合的结果，对于新的决策树，则直接添加到集成中。设置为False的时候，擦除之前的结果，并重新开始计算，如果不设置这个参数，那么就默认为False。

Early stop 早停

  实际上也可以提前停止训练。而不是先训练大量的树然后再回去查找到最优estimator的数量。

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
X_train,X_val,y_train,y_val=train_test_split(X,y)
gbrt=GradientBoostingRegressor(max_depth=2,n_estimators=120)
gbrt.fit(X_train,y_train)
errors=[mean_squared_error(y_val,y_pred) for y_pred in gbrt.staged_predict(X_val)]
bst_n_estimators=np.argmin(errors)
gbrt_best=GradientBoostingRegressor(max_depth=2,n_estimators=bst_n_estimators)
gbrt_best.fit(X_train,y_train)

axes=plt.gca()
axes.set_ylim([0,0.01])
axes.set_xlim([0,120])
plt.plot(np.arange(120),errors,label='Validation error')
plt.plot((0,120),(errors[bst_n_estimators],errors[bst_n_estimators]),linestyle='dashed')
plt.plot((bst_n_estimators,bst_n_estimators),(0,errors[bst_n_estimators]),linestyle='dashed')
plt.scatter(bst_n_estimators,errors[bst_n_estimators],color='red',label='Minimum')
plt.legend()
plt.xlabel('Number of trees')
plt.ylabel('Validation error')

  上图中，可以看到随着树的数量增加，error的值是越来越小的，但是经过某一个点之后，也就是说当决策树的数量大于某值的时候，error的值并不是逐渐减小的，此时增加树的数量并不能减少误差，那么这个数值就是当前状况下，最优的决策树数量。

  通过计算在不同estimator的数量下，预测的结果和实际结果的误差，得到的误差结果最小的模型中对应的estimator的数量为74。然后将n_estimators的数量设置为74，然后再使用GradientBoostingRegressor进行拟合，得到的结果如下：

gbrt=GradientBoostingRegressor(max_depth=2,n_estimators=bst_n_estimators,learning_rate=0.1)
gbrt.fit(X,y)
df=pd.DataFrame({"X":X[:,0],"y":gbrt.predict(X)})
df=df.sort_values(by=['X'])
plt.scatter(X,y,c=y)
line1,=plt.plot(df.X,df.y,c='red',label='Ensemble predictions')
plt.legend(handles=[line1],loc='upper center')

本文地址：https://blog.csdn.net/weixin_44746091/article/details/110872048