机器学习-Sklearn（第三版）Day11 XGBoost

sklearn与XGBoost

momokofly 于2022-03-14 14:50:34

• sklearn

1 在学习XGBoost之前

1.1 xgboost库与XGB的sklearn API

陈天奇创造了XGBoost算法后，很快和一群机器学习爱好者建立了专门调用XGBoost库，名为xgboost。xgboost是一个独立的、开源的，并且专门提供梯度提升树以及XGBoost算法应用的算法库。它和sklearn类似，有一个详细的官方网站可以提供学习资料，并且可以与C、Python、R、Julia等语言连用，但需要单独安装和下载。
xgboost documents：XGBoost Documentation — xgboost 2.0.0-dev documentation
本文基于Pythonl来运行。xgboost库要求必须提供适合的Scipy环境。接下来提供在windows中使用pip安装xgboost的代码：

pip install xgboost  #安装xgboost库
pip install --upgrade xgboost   #更新xgboost库

安装完毕之后，就能够使用这个库中所带的XGB相关的类了。

import xgboost as xgb

现在有两种方式来使用xgboost库。第一种方式，是直接使用xgboost库自己的建模流程。

其中最核心的是DMatrix这个读取数据的类，以及train()这个用于训练的类。与sklearn把所有的参数都写在类中的方式不同，xgboost库中必须先使用字典设定参数集，再使用train来将参数集输入，然后进行训练。会这样设计的原因，是因为XGB所涉及到的参数实在太多，全部写在xgb.train()太长也容易出错。下面给出params可能的取值以及xgboost.train的列表。
params {eta, gamma, max_depth, min_child_weight, max_delta_step, subsample, colsample_bytree, colsample_bylevel, colsample_bynode, lambda, alpha, tree_method_string, sketch_eps, scale_pos_weight, updater, refresh_leaf, process_type, grow_policy, max_leaves, max_bin, predictor, num_parallel_tree}
xgboost.train(params, dtrain, num_boost_round = 10, evals = (), obj = None, feval = None, maximize = False, early_stopping_rounds = None, evals_result = None, verbose_eval = True, xgb_model = None, callbacks = None, learning_rates = None)
或者，也可以选择第二种方法，使用xgboost库中的sklearn中的API。这就是说，可以调用如下的类，并用sklearn当中的惯例的实例化，fit和predict的流程来运行XGB，并且也可以调用属性（比如coef_）等等。当然，这是回归的类，也有用于分类、排序的类。他们与回归的类非常相似，因此了解一个类即可。
xgboost.XGBRegressor(max_depth = 3, learning_rate = 0.1, n_estimators = 100, silent = True, objective = ‘reg:linear’, booster = ‘gbtree’, n_jobs = 1, nthread = None, gamma = 0, min_child_weight = 1, max_delta_step = 0, subsample = 1, colsample_bytree = 1, colsample_bylevel = 1, reg_alpha = 0, reg_lambda = 1, scale_pos_weight = 1, base_score = 0.5, random_state = 0, seed = None, missing = None, importance_type = ‘gain’, ∗ ∗ ** ∗∗kwargs)
调用xgboost.train和调用sklearn API中的类XGBRegressor，需要输入的参数是不同的，而且看起来相当的不同。但其实，这些参数只是写法不同，功能是相同的。如params字典中的第一个参数eta，其实就是XGBRegressor里面的参数learning_rate，他们的含义和实现的功能是一模一样的。只不过在sklearn API中，开发团队友好的将参数的名称调节成了与sklearn中其他算法类更相似的样子。
所以，使用xgboost中设定的建模流程来建模，和使用sklearn API中的类来建模，模型效果是比较相似的，但是xgboost库本身的运算速度（尤其是交叉验证）以及调参手段比sklearn要简单。

1.2 XGBoost的三大板块

XGBoost本身的核心是基于梯度提升树实现的集成算法，整体来说可以有三个核心部分：集成算法本身，用于集成的弱评估器，以及应用中的其他过程。三个部分中，前两个部分包含了XGBoost的核心原理以及数学过程，最后的部分主要是在XGBoost应用中占有一席之地。

参数	集成算法	弱评估器	其他过程
n_estimators	∗ \ast ∗
learning_rate	∗ \ast ∗
silent(verbosity)	∗ \ast ∗
subsample	∗ \ast ∗
max_depth		∗ \ast ∗
objective		∗ \ast ∗
booster		∗ \ast ∗
gamma		∗ \ast ∗
min_child_weight		∗ \ast ∗
max_delta_step		∗ \ast ∗
colsample_bytree		∗ \ast ∗
colsample_bylevel		∗ \ast ∗
reg_alpha		∗ \ast ∗
reg_lambda		∗ \ast ∗
nthread			∗ \ast ∗
n_jobs			∗ \ast ∗
scale_pos_weight			∗ \ast ∗
base_score			∗ \ast ∗
seed			∗ \ast ∗
random_state			∗ \ast ∗
missing			∗ \ast ∗
importance_type			∗ \ast ∗

2 梯度提升树

2.1 提升集成算法：重要参数n_estimators

XGBoost的基础是梯度提升算法，因此必须先从了解梯度提升算法开始。梯度提升（Gradient boosting）算法是构建预测模型的最强大技术之一，它是集成算法中提升法（Boosting）的代表算法。集成算法通过在数据集上建立多个弱评估器，汇总所有弱评估器的预测模型，以获取比单个模型更好的回归或分类表现。弱评估器被要求至少比随机猜测更好的模型，即预测准确度不低于50%的任意模型。
集成不同弱评估器的方法有很多种，有一次性建立多个平行独立的弱评估器的装袋法（Bagging），也有逐一构建弱评估器，经过多次迭代逐渐累积多个弱评估器的提升法（Boosting）。提升法中最著名的算法包括Adaboost和梯度提升树，XGBoost就是由梯度提升树发展而来的。梯度提升树中可以有回归树也可以有分类树，两者都是以CART树算法作为主流，XGBoost背后也是CART树，这意味着XGBoost中所有的树都是二叉的。以梯度提升回归树为例，了解一下Boosting算法的工作原理。
首先梯度提升回归树是专注于回归的树模型的提升集成模型，其建模过程大致如下：最开始先建立一棵树，然后逐渐迭代，每次迭代过程中都增加一棵树，逐渐形成众多树模型集成的强评估器。

对于决策树而言，每个被放入模型的任意样本 i i i最终一定都会落到一个叶子节点上。而对于回归树，每个叶子节点上的值是这个叶子节点上所有样本的均值。

对于梯度提升树回归树来说，每个样本的预测结果可以表示为所有树上的结果的加权求和： y ^ i ( k ) = ∑ k K γ k h k ( x i ) \hat y_i^{(k)}=\sum_k^K\gamma_kh_k(x_i) y^​i(k)​=k∑K​γk​hk​(xi​)
其中， K K K是树的总数量， k k k代表第 k k k棵树， γ k \gamma_k γk​是这棵树的权重， h k h_k hk​表示这棵树上的预测结果。
值得注意的是，XGB作为GBDT的改进，在 y ^ \hat y y^​上却有所不同。对于XGB来说，每个叶子节点上会有一个预测分数（prediction score），也被称为叶子权重。这个叶子权重就是所有在这个叶子节点上的样本在这一棵树上的回归取值，用 f k ( x i ) f_k(x_i) fk​(xi​)或者 w w w来表示，其中 f k f_k fk​表示第 k k k棵决策树， x i x_i xi​表示样本 i i i对应的特征向量。当只有一棵树的时候， f 1 ( x i ) f_1(x_i) f1​(xi​)就是提升集成算法返回的结果，但这个结果往往非常糟糕。当有多棵树的时候，集成模型的回归结果就是所有树的预测分数之和，假设这个集成模型中总共有 K K K棵决策树，则整个模型在这个样本 i i i上给出的预测结果为： y ^ i ( k ) = ∑ k K f k ( x i ) \hat y_i^{(k)}=\sum_k^Kf_k(x_i) y^​i(k)​=k∑K​fk​(xi​)

XGB和GBDT核心区别1：求解预测值 y ^ \hat y y^​的方式不同
GBDT中预测值是由所有弱分类器上的预测结果的加权求和，其中每个样本上的预测结果就是样本所在的叶子节点的均值。而GBDT中的预测值是所有弱分类器上的叶子权重直接求和得到，计算叶子权重是一个复杂的过程。

从上式来看，在集成中需要考虑的第一件事是超参数 K K K，究竟要建多少棵树？

参数含义	xgb.train()	xgb.XGBRegressor()
集成中弱评估器的数量	num_round，默认为10	n_estimators，默认为100
训练者是否打印每次训练的结果	silent，默认为False	silent，默认True

在最新版本的xgboost中没有参数silent，改为verbosity，可以取值0-3，其中取0时表示silent。
在随机森林中理解n_estimators：n_estimators越大，模型的学习能力就会越强，模型也越容易过拟合。在随机森林中，调整的第一个参数就是n_estimators，这个参数非常强大，常常能够一次性将模型调整到极限。在XGB中，也期待相似的表现，虽然XGB的集成方式与随机森林不同，但使用更多的弱分类器来增强模型整体的学习能力这件事是一致的。

#导入需要的库，模块以及数据
from xgboost import XGBRegressor as XGBR
from sklearn.ensemble import RandomForestRegressor as RFR
from sklearn.linear_model import LinearRegression as LinearR
from sklearn.datasets import load_boston
from sklearn.model_selection import KFold,cross_val_score as CVS,train_test_split as TTS
from sklearn.metrics import mean_squared_error as MSE
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from time import time
import datetime

data = load_boston()#波士顿数据集非常简单，但它涉及到的问题却很多
#data
x = data.data
y = data.target
x.shape
#结果：(506, 13)
y.shape
#结果：(506,)

#建模，查看其他接口和属性
xtrain,xtest,ytrain,ytest = TTS(x,y,test_size = 0.3, random_state = 420)
reg = XGBR(n_estimators = 100).fit(xtrain,ytrain)#训练
reg.predict(xtest)#传统接口predict
'''
结果：
array([ 6.6689262, 22.34918  , 31.052807 , 13.911595 ,  9.467966 ,
       22.658588 , 14.514282 , 15.092699 , 15.293644 , 12.680115 ,
       24.140797 , 35.890083 , 21.573483 , 27.07066  , 19.052658 ,
        9.89033  , 23.386076 , 23.588493 , 23.311466 , 22.401644 ,
       18.98444  , 15.766946 , 25.8352   , 20.193802 , 19.982517 ,
       15.611423 , 22.883228 , 29.838228 , 22.815304 , 16.779037 ,
       37.13194  , 20.133305 , 19.67352  , 23.525528 , 22.845137 ,
       23.87397  , 15.17887  , 23.45934  , 16.685331 , 31.761686 ,
       18.525843 , 22.441063 , 38.48728  , 17.93719  , 15.10122  ,
       28.980541 , 46.363487 , 12.842797 ,  9.618281 , 35.40579  ,
       25.657566 , 20.605602 , 20.800055 , 49.228447 , 31.355848 ,
       29.382515 , 18.911947 , 21.049877 , 16.165182 , 18.098577 ,
       14.659002 , 21.720213 , 19.413454 , 28.932102 , 30.573524 ,
       19.228426 , 20.531511 , 15.666289 , 23.52929  , 19.30554  ,
       28.384985 , 42.83562  , 29.429724 , 23.306015 , 19.741224 ,
       24.202463 , 38.735516 , 26.782232 , 22.168324 , 20.67139  ,
       19.534992 , 47.360317 , 24.008802 ,  8.184814 , 25.613977 ,
       20.691843 , 17.127483 , 21.10027  , 22.279167 ,  7.755469 ,
       20.0476   , 15.406112 , 42.38165  , 33.828186 , 22.788246 ,
        9.302419 , 10.416187 , 20.033014 ,  8.241165 , 12.816478 ,
       30.793974 , 10.078776 , 25.383692 , 21.933594 , 30.53568  ,
       42.866497 , 19.598145 ,  8.321976 , 23.194368 , 12.547767 ,
       46.838867 , 22.961082 , 20.244467 , 23.170694 , 18.948856 ,
       29.682056 , 24.363865 , 19.715958 , 44.975193 , 17.64582  ,
       24.3169   , 24.946495 , 18.23235  , 16.922577 , 14.77714  ,
       21.757038 , 33.83876  , 10.938419 , 20.035763 , 21.085218 ,
       19.331562 , 20.505526 ,  8.285518 , 28.01946  , 29.631227 ,
       19.908175 , 18.30325  , 14.204149 , 10.795732 , 18.232307 ,
       42.266888 , 17.304502 , 22.974121 , 20.946724 , 30.724297 ,
       20.072989 , 12.772859 , 41.472908 , 27.652838 , 22.20476  ,
       14.235871 , 25.88964  ], dtype=float32)
'''
reg.score(xtest,ytest)#R^2，shift+tab可以看函数具体原理
#结果：0.9050988968414799
MSE(ytest,reg.predict(xtest))
#结果：8.830916343629323
reg.feature_importances_#树模型的优势之一：能够查看模型的重要性分数，可以使用嵌入法(SelectFromModel)进行特征选择
'''
结果：
array([0.01902167, 0.0042109 , 0.01478316, 0.00553537, 0.02222196,
       0.37914088, 0.01679686, 0.0469872 , 0.04073574, 0.05491759,
       0.06684221, 0.00869464, 0.3201119 ], dtype=float32)
'''
#交叉验证，与线性回归&随机森林回归进行对比
reg = XGBR(n_estimators = 100)#交叉验证中导入没有经过训练的模型
CVS(reg,xtrain,ytrain,cv = 5).mean()#这里返回的模型评估指标是R^2
#严谨的交叉验证与不严谨的交叉验证之间的讨论：训练集 or 全数据？
#结果：0.7995062821902295
CVS(reg,xtrain,ytrain,cv = 5,scoring = "neg_mean_squared_error").mean()
#结果：-16.215644229762717
#来查看一下sklearn中所有的模型评估指标
import sklearn
sorted(sklearn.metrics.SCORERS.keys())
'''
结果：
['accuracy',
 'adjusted_mutual_info_score',
 'adjusted_rand_score',
 'average_precision',
 'balanced_accuracy',
 'completeness_score',
 'explained_variance',
 'f1',
 'f1_macro',
 'f1_micro',
 'f1_samples',
 'f1_weighted',
 'fowlkes_mallows_score',
 'homogeneity_score',
 'jaccard',
 'jaccard_macro',
 'jaccard_micro',
 'jaccard_samples',
 'jaccard_weighted',
 'max_error',
 'mutual_info_score',
 'neg_brier_score',
 'neg_log_loss',
 'neg_mean_absolute_error',
 'neg_mean_absolute_percentage_error',
 'neg_mean_gamma_deviance',
 'neg_mean_poisson_deviance',
 'neg_mean_squared_error',
 'neg_mean_squared_log_error',
 'neg_median_absolute_error',
 'neg_root_mean_squared_error',
 'normalized_mutual_info_score',
 'precision',
 'precision_macro',
 'precision_micro',
 'precision_samples',
 'precision_weighted',
 'r2',
 'rand_score',
 'recall',
 'recall_macro',
 'recall_micro',
 'recall_samples',
 'recall_weighted',
 'roc_auc',
 'roc_auc_ovo',
 'roc_auc_ovo_weighted',
 'roc_auc_ovr',
 'roc_auc_ovr_weighted',
 'top_k_accuracy',
 'v_measure_score']
'''
#使用随机森林和线性回归进行一个对比
rfr = RFR(n_estimators = 100)
CVS(rfr,xtrain,ytrain,cv = 5).mean()
#结果：0.7922215794897552
CVS(rfr,xtrain,ytrain,cv = 5,scoring = "neg_mean_squared_error").mean()
#结果：-16.71229234925553
lr = LinearR()
CVS(lr,xtrain,ytrain,cv = 5).mean()
#结果：0.6835070597278079
CVS(lr,xtrain,ytrain,cv = 5,scoring = "neg_mean_squared_error").mean()
#结果：-25.349507493648463
#如果开启参数silent：在数据巨大，预料到算法运行会非常缓慢的时候可以使用这个参数来监控模型的训练进度
#最新版本的sklearn中的XGBR的参数中没有silent这一参数，改为verbosity
#verbosity : Optional[int]        The degree of verbosity. Valid values are 0 (silent) - 3 (debug).
reg = XGBR(n_estimators = 10, verbosity = 0)
CVS(reg,xtrain,ytrain,cv = 5,scoring = 'neg_mean_squared_error').mean()
#结果为：-18.633733952333067

#定义绘制以训练样本数为横坐标的学习曲线的函数
def plot_learning_curve(estimator, title, x, y
                       ,ax = None  #选择子图
                       ,ylim = None  #设置纵坐标的取值范围
                       ,cv = None  #交叉验证
                       ,n_jobs = None  #设定索要使用的线程
                       ):
    from sklearn.model_selection import learning_curve
    import matplotlib.pyplot as plt
    import numpy as np
    train_sizes, train_scores, test_scores = learning_curve(estimator, x, y
                                                           ,shuffle = True
                                                           ,cv = cv
                                                           #,random_state = 420
                                                            #random_state固定时绘制的图像不会每次运行都变化
                                                           ,n_jobs = n_jobs)
    if ax == None:
        ax = plt.gca()
    else:
        ax = plt.figure()
    ax.set_title(title)
    if ylim is not None:
        ax.set_ylim(*ylim)
    ax.set_xlabel("Training examples")
    ax.set_ylabel("Score")
    ax.grid()#绘制网格，不是必须
    ax.plot(train_sizes, np.mean(train_scores, axis = 1),"o-", color = "r", label = "Training score")
    ax.plot(train_sizes, np.mean(test_scores,axis = 1), "o-", color = "g", label = "Test score")
    ax.legend(loc = "best")
    return ax
#使用学习曲线观察XGB在波士顿数据集上的潜力
cv = KFold(n_splits = 5, shuffle = True, random_state = 42)#交叉验证模式
plot_learning_curve(XGBR(n_estimators = 100, random_state = 420)
                   ,"XGB", xtrain,ytrain,ax = None,cv = cv)
plt.show()
#多次运行，观察结果，这是怎么造成的？——random_state
#在现在的状况下，如何看数据的潜力，还能调上去吗？

训练集上的表现展示了模型的学习能力，测试集上的表现展示了模型的泛化能力，通常模型在测试集上的表现不太可能超过训练集，因此希望测试集的学习曲线能够努力逼近训练集的学习曲线。观察三种学习曲线组合：希望将模型调整成什么样？能够将模型调整成什么样？

#使用参数学习曲线观察n_estimators对模型的影响
axisx = range(10,1010,50)
rs = []
for i in axisx:
    reg = XGBR(n_estimators = i, random_state = 420)
    rs.append(CVS(reg,xtrain,ytrain,cv = cv).mean())
print(axisx[rs.index(max(rs))],max(rs))
plt.figure(figsize = (20,5))
plt.plot(axisx,rs,c = "red",label = "XGB")
plt.legend()
plt.show()
#选出来的n_estimators非常不寻常，是否要选择准确率最高的n_estimators值呢？
#结果为：160 0.8320776498992342

进化的学习曲线：方差与泛化误差
在随机森林中有方差-偏差困境。在机器学习中，用来衡量模型在未知数据上的准确率的指标，叫做泛化误差（Genelization error）。一个集成模型(f)在未知数据集(D)上的泛化误差 E ( f ; D ) E(f;D) E(f;D)，由方差(var)，偏差(bias)和噪声(e)共同决定。其中偏差就是训练集上的拟合程度决定，方差是模型的稳定性决定，噪音是不可控的。而泛化误差越小，模型就越理想。
E ( f ; D ) = b i a s 2 + v a r + e 2 E(f;D)=bias^2+var+e^2 E(f;D)=bias2+var+e2
过去往往直接取学习曲线获得的分数的最高点，即考虑偏差最小的点，是因为模型极度不稳定，方差很大的情况其实比较少见。但是现在的数据量非常少，模型会相对不稳定，因此应当方差也纳入考虑的范围。在绘制学习曲线时，不仅要考虑偏差的大小，还要考虑方差的大小，更要考虑泛化误差中可控的部分。当然，并不是说可控的部分比较小，整体的泛化误差就一定小，因为误差有时候可能占主导。基于这种思路改进学习曲线：

axisx = range(50,1050,50)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators = i, random_state = 420)
    cvresult = CVS(reg, xtrain, ytrain, cv = cv)
    #记录  1-偏差
    rs.append(cvresult.mean())
    #记录  方差
    var.append(cvresult.var())
    #计算泛化误差的可控部分
    ge.append((1-cvresult.mean())**2 + cvresult.var())
#打印R2最高所对应的参数取值，并打印这个参数下的方差
print(axisx[rs.index(max(rs))], max(rs), var[rs.index(max(rs))])
#打印方差最低时对应的参数取值，并打印这个参数下的R2
print(axisx[var.index(min(var))], rs[var.index(min(var))], min(var))
#打印泛化误差可控部分的参数取值，并打印这个参数下的R2，方差以及泛化误差的可控部分
print(axisx[ge.index(min(ge))], rs[ge.index(min(ge))], var[ge.index(min(ge))], min(ge))
plt.figure(figsize = (20,5))
plt.plot(axisx, rs, c = "red", label = "XGB")
plt.legend()
plt.show()
'''
结果：
100 0.8320924293483107 0.005344212126112929
100 0.8320924293483107 0.005344212126112929
100 0.8320924293483107 0.005344212126112929 0.03353716440826495
'''

#细化学习曲线，找出最佳n_estimators
axisx = range(50,250,10)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators = i, random_state = 420)
    cvresult = CVS(reg, xtrain, ytrain, cv = cv)
    rs.append(cvresult.mean())
    var.append(cvresult.var())
    ge.append((1-cvresult.mean())**2 + cvresult.var())
print(axisx[rs.index(max(rs))], max(rs), var[rs.index(max(rs))])
print(axisx[var.index(min(var))], rs[var.index(min(var))], min(var))
print(axisx[ge.index(min(ge))], rs[ge.index(min(ge))], var[ge.index(min(ge))], min(ge))
rs = np.array(rs)
var = np.array(var)*0.01
plt.figure(figsize = (20,5))
plt.plot(axisx,rs,c = "black",label = "XGB")
#添加方差线
plt.plot(axisx, rs+var, c = "red", linestyle = "-.")
plt.plot(axisx, rs-var, c = "red", linestyle = "-.")
plt.legend()
plt.show()
'''
结果：
100 0.8320924293483107 0.005344212126112929
100 0.8320924293483107 0.005344212126112929
100 0.8320924293483107 0.005344212126112929 0.03353716440826495
'''