十一、XGBoost

11.1 梯度提升树

class 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)

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

XGBoost就是由梯度提升树发展而来的，既可以回归，也可以分类。XGBoost中所有的树都是二叉的；
接下来以回归树为例来进行介绍梯度提升：

• 首先，梯度提升回归树是专注于回归的树模型的提升集成模型，其建模过程大致如下：最开始先建立一棵树，然后逐渐迭代，每次迭代过程中都增加一棵树，逐渐形成众多树模型集成的强评估器；

每一棵树之间不是相互独立，而是相互递进的关系

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

对于梯度提升回归树来说，每个样本的预测结果可以表示为所有树上的结果的加权求和：

其中，K是树的总数量，k代表第k棵树，r_k是这棵树的权重，h_k表示这棵树上的预测结果，是均值；

• 值得注意的是，XGB作为GBDT的改进，在y^ 上却有所不同：对于XGB来说，每个叶子节点上会有一个预测分数/叶子权重。这个叶子权重就是所有在这个叶子节点上的样本在这一棵树上的回归取值，用f_k(x_i)或者w来表示（要另外计算），其中f_k表示第k棵决策树，x_i表示样本i对应的特征向量；
• 当只有一棵树的时候，就是 提升集成算法返回的结果，但这个结果往往非常糟糕。当有多棵树的时候，集成模型的回归结果就是 所有树的预测分数之和，假设这个集成模型中总共有 棵决策树，则整个模型在这个样本 上给出的预测结果为：
  

从式子中不难看出，需要确定的是K，即集成模型中决策树的个数，可以作为超参数来人为控制

• 常用参数：
  

若直接用xgb的话记得把10改大点。在XGB中，我们也期待相似的表现，虽然XGB的集成方式与随机森林不同，但使用更多的弱分类器来增强模型整体的学习能力这件事是一致的

• 先建个模试试：

1. 导入需要的库，模块以及数据

#直接从xgb导入回归
from xgboost import XGBRegressor as XGBR
#从sklearn导入回归来进行对比
from sklearn.ensemble import RandomForestRegressor as RFR
from sklearn.linear_model import LinearRegression as LinearR
from sklearn.datasets import load_boston
#KFold交叉验证
from sklearn.model_selection import KFold, cross_val_score as CVS, train_test_split as TTS
#MSE均方误差，作为模型评估指标
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()
#波士顿数据集非常简单，但它所涉及到的问题却很多

X = data.data
y = data.target

2. 建模，查看其他接口和属性

Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420)

reg = XGBR(n_estimators=100).fit(Xtrain,Ytrain)
#传统接口predict
reg.predict(Xtest)
#你能想出这里应该返回什么模型评估指标么？
#返回的是r²
reg.score(Xtest,Ytest)
#0.9050988968414799

MSE(Ytest,reg.predict(Xtest))
#8.830916343629323

#树模型的优势之一：能够查看模型的重要性分数，可以使用嵌入法进行特征选择，分数越大说明这个特征贡献越大
reg.feature_importances_
#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)

3. 交叉验证，与线性回归&随机森林回归进行对比

#交叉验证要导入没有训练过的模型
reg = XGBR(n_estimators=100)
CVS(reg,Xtrain,Ytrain,cv=5).mean()
#0.7995062821902295
#这里应该返回什么模型评估指标，还记得么？和模型本身的评估接口返回相同的指标奥，这里是r²
#严谨的交叉验证与不严谨的交叉验证之间的讨论：训练集or全数据？

#严谨与否？
#在CVS时若输入全数据会被认作不严谨，因为其中也有数据被作为测试集使用，这样就会让模型提前知道了
#但实际上用不用全数据无所谓，都可以

#在均方误差上的结果，使用scoring来要求其用指定指标来进行评估
CVS(reg,Xtrain,Ytrain,cv=5,scoring='neg_mean_squared_error').mean()
#-16.215644229762717

#来查看一下sklearn中所有的模型评估指标,可以加到scoring中来指定评估指标
import sklearn
sorted(sklearn.metrics.SCORERS.keys())

#使用随机森林和线性回归进行一个对比
rfr = RFR(n_estimators=100)
CVS(rfr,Xtrain,Ytrain,cv=5).mean()
#0.8061118332352855

CVS(rfr,Xtrain,Ytrain,cv=5,scoring='neg_mean_squared_error').mean()
#-16.669141208008043

#线性回归
lr = LinearR()
CVS(lr,Xtrain,Ytrain,cv=5).mean()
#0.6835070597278081

CVS(lr,Xtrain,Ytrain,cv=5,scoring='neg_mean_squared_error').mean()
#-25.349507493648456
#线性回归，不太行

#开启参数slient：在数据巨大，预料到算法运行会非常缓慢的时候可以使用这个参数来监控模型的训练进度
reg = XGBR(n_estimators=10,silent=False)
CVS(reg,Xtrain,Ytrain,cv=5,scoring='neg_mean_squared_error').mean()

设置为False后就会展示每一棵树的状况，此过程可以帮助理解模型建模到了什么程度

4. 定义绘制以训练样本数为横坐标的学习曲线的函数

#4. 定义绘制以训练样本数为横坐标的学习曲线的函数

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
                                                           ,n_jobs=n_jobs)
    if ax == None:
        ax = plt.gca()
    else:
        ax = plt.figure()
    ax.set_title(tile)
    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 examples")
    ax.plot(train_sizes, np.mean(test_scores, axis=1),'o-'
           , color="g",label="Test examples")
    ax.legend(loc="best")
    return ax

5. 使用学习曲线观察XGB在波士顿数据集上的潜力

#5. 使用学习曲线观察XGB在波士顿数据集上的潜力

#定义一个交叉验证模式，n_split表示分几份，shuffle表示是否打乱数据集，
cv = KFold(n_splits=5, shuffle = True, random_state=4)
plot_learning_curve(XGBR(n_estimators=100,random_states=420)
                   ,"XGB",Xtrain,Ytrain,ax=None,cv=cv)
plt.show()

多次运行上面的代码会发现每次的图像都不一样，这是因为在自定义函数中没有固定随机取样策略；

训练集上的表现展示了模型的学习能力，测试集上的表现展示了模型的泛化能力，通常模型在测试集上的表现不太可能超过训练集，因此我们希望我们的测试集的学习曲线能够努力逼近我们的训练集的学习曲线；

我们希望两条曲线要尽可能的靠近，如下图中的中间和右边的状况：

中间这条线是让训练集分数下降，测试集分数上升从而达到消除过拟合；
右边的则是不断提升测试的分数；

还是右边的那种好，但是并非所有模型都可以调成这样。大多数时候还是中间这样，

6. 使用参数学习曲线观察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值呢？

可以看到，这个数在图像中的位置偏大，且效果提升不明显

7. 进化的学习曲线：方差与泛化误差

• 泛化误差：衡量模型在未知数据上的准确率的指标。一个集成模型(f)在未知数据集(D)上的泛化误差E(f;D),由方差(var)，偏差(bais)和噪声(ε)共同决定：
  

其中偏差就是训练集上的拟合程度决定，方差是模型的稳定性决定，噪音是不
可控的。而泛化误差越小，模型就越理想

现在我们的数据量非常少，模型会相对不稳定，因此我们应当将方差也纳入考虑的范围。在绘制学习曲线时，我们不仅要考虑偏差的大小，还要考虑方差的大小，更要考虑泛化误差中我们可控的部分。用这种思路来改进学习曲线：

#7. 进化的学习曲线：方差与泛化误差

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-偏差，即r²
    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这个点的泛化误差最小

8. 细化学习曲线，找出最佳n_estimators

axisx = range(100,300,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()

#看看泛化误差的可控部分如何？
plt.figure(figsize=(20,5))
plt.plot(axisx,ge,c="gray",linestyle='-.')
plt.show()

9. 检测模型效果

#验证模型效果是否提高了？
time0 = time()
print(XGBR(n_estimators=100,random_state=420).fit(Xtrain,Ytrain).score(Xtest,Ytest))
print(time()-time0)
#0.9050988968414799
#0.059462547302246094

time0 = time()
print(XGBR(n_estimators=110,random_state=420).fit(Xtrain,Ytrain).score(Xtest,Ytest))
print(time()-time0)
#0.9050738056164094
#0.07250857353210449

time0 = time()
print(XGBR(n_estimators=120,random_state=420).fit(Xtrain,Ytrain).score(Xtest,Ytest))
print(time()-time0)
#0.905061100288288
#0.06450581550598145

这个结果和教程不太一样，不过这种曲线也不一定永远有效。若数据集很大则不太好用这个方法，先用其他方法(网格搜索等)后再调整树的数量；

• 结论：

1. XGB中的树的数量决定了模型的学习能力，树的数量越多，模型的学习能力越强。在数量很多时，哪怕数据量很少，也可以很好的学习；
   2.XGB中树的数量很少的时候，对模型的影响较大，当树的数量已经很多的时候，对模型的影响比较小，只能有微弱的变化。当唯一指标R²或者准确率给出的n_estimators看起来不太可靠的时候，我们可以改造学习曲线来帮助我们；
   3.树的数量提升对模型的影响有极限，最开始，模型的表现会随着XGB的树的数量一起提升，但到达某个点之后，树的数量越多，模型的效果会逐步下降，这也说明了暴力增加n_estimators不一定有效果；

建议优先调整n_estimators，一般都不会建议一个太大的数目，300以下为佳；

11.1.2 有放回随机抽样：重要参数subsample

我们训练模型之前，必然会有一个巨大的数据集。我们都知道树模型是天生过拟合的模型，并且如果数据量太过巨大，树模型的计算会非常缓慢，因此，我们要对我们的原始数据集进行有放回抽样bootstrap：有放回的抽样每次只能抽取一个样本，若我们需要总共N个样本，就需要抽取N次。每次抽取一个样本的过程是独立的，这一次被抽到的样本会被放回数据集中，下一次还可能被抽到，因此抽出的数据集中，可能有一些重复的数据；

在抽样时可能出现相同的数据，这样就可以有效防止过拟合。因为我们人为打乱了数据信息

在梯度提升树中，我们每一次迭代都要建立一棵新的树，因此我们每次迭代中，都要有放回抽取一个新的训练样本。但这并不能保证每次建模的效果都比之前那次要好。因此我们规定，在梯度提升树中，每构建一个评估器，都让模型更加集中于数据集中容易被判错的那些样本：

首先我们有一个巨大的数据集，在建第一棵树时，我们对数据进行初次又放回抽样，然后建模。建模完毕后，我们对模型进行一个评估，然后将模型预测错误的样本反馈给我们的数据集，一次迭代就算完成。紧接着，我们要建立第二棵决策树，于是开始进行第二次又放回抽样。但这次有放回抽样，和初次的随机有放回抽样就不同了，在这次的抽样中，我们加大了被第一棵树判断错误的样本的权重。也就是说，被第一棵树判断错误的样本，更有可能被我们抽中
好像就是adaboost的原理？

我们相信，只要弱分类器足够强大，随着模型整体不断在被判错的样本上发力，这些样本会渐渐被判断正确；

在sklearn中，我们使用参数subsample来控制我们的随机抽样，代表随机抽样样本的比例。在xgb和sklearn中，这个参数都默认为1且不能取到0，这说明我们无法控制模型是否进行随机有放回抽样，只能控制抽样抽出来的样本量大概是多少：

正常来说样本量越大，模型才不容易过拟合，现在展现出来的效果，是由于我们的样本量太小造成的一个巧合。从这个角度来看，我们的subsample参数对模型的影响应该会非常不稳定，大概率应该是无法提升模型的泛化能力的，但也不乏提升模型的可能性：

axisx = np.linspace(0,1,20)
rs = []
for i in axisx:
    reg = XGBR(n_estimators=110,subsample=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="green",label="XGB")
plt.legend()
plt.show()

从上图无法看出什么，所以继续细化画：

#细化学习曲线
axisx = np.linspace(0.5,1,20)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators=110,subsample=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()

将上述区域进行放大后可以看出不再是直线了，而且更可以看出哪些点比较好，那么进行进一步细化

#继续细化学习曲线
axisx = np.linspace(0.7,1,25)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators=110,subsample=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()

说明80%会有最好的表现，导入模型看看结果

#不要盲目找寻泛化误差可控部分的最低值，注意观察结果
#看看泛化误差的情况如何
reg = XGBR(n_estimators=110
           ,subsample=0.7999999999999999
           ,random_state=420).fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#0.910883452569069
MSE(Ytest,reg.predict(Xtest))
#8.292640959937806
#这样的结果说明了什么？

并没有提升多少（虽然教程里头是降低）
总体来说这个参数并没有对波士顿房价数据集上的结果造成太大的影响，由于我们的数据集过少，降低抽样的比例反而让数据的效果更低，不如就让它保持默认

11.1.3 迭代决策树：重要参数eta

从数据的角度而言，我们让模型更加倾向于努力攻克那些难以判断的样本。但是，困难样本被加重权重是因为前面的树没能把它判断正确，所以对于下一棵树来说，它要判断的测试集的难度，是比之前的树所遇到的数据的难度都要高的；

除了保证模型逐渐倾向于困难样本的方向，我们必须保证，每次新添加的树一定得是对这个新数据集预测效果最优的那一棵树；

枚举是指不断随机抽取特征，就可以有不同的组合，对于每个组合都建模，从而找到其中最好的，但是很耗时耗力；
sklearn的那种无法保证每次建的树中一定有最好的树，而这也很难实现；

我们希望能找出一种方法，直接帮我们求解出最优的集成算法结果——>最优化问题？

思路：我们首先找到一个损失函数 ，这个损失函数应该可以通过带入我们的预测结果 来衡量我们的梯度提升树在样本的预测效果。然后，我们利用梯度下降来迭代我们的集成算法；

我们让这个过程持续下去，直到找到能够让损失函数最小化的 ，这个 就是我们模型的预测结果。在XGB中，我们完整的迭代决策树的公式应该写作：

其中η读作"eta"，是迭代决策树时的步长（shrinkage），又叫做学习率（learning rate）

和逻辑回归中的α类似，η越大，迭代的速度越快，算法的极限很快被达到，有可能无法收敛到真正的最佳。η越小，越有可能找到更精确的最佳值，更多的空间被留给了后面建立的树，但迭代速度会比较缓慢：

在sklearn中，我们使用参数learning_rate来干涉我们的学习速率：

eta取0表示不进行迭代

在实际应用中，往往使用网格搜索来同时调节两个参数。

#参数eta

def regassess(reg,Xtrain,Ytrain,cv,scoring = ["r2"], show=True):
    score = []
    for i in range(len(scoring)):
        if show:
            print("{}:{:.2f}".format(scoring[i]
                                    ,csv(reg
                                        ,Xtrain,Ytrain
                                        ,cv=cv,scoring=scoring[i]).mean()))
        score.append(CVS(reg,Xtrain,Ytrain,cv=cv,scoring=scoring[i]).mean())
    return score

#运行一下函数来看看效果
regassess(reg,Xtrain,Ytrain,cv,scoring = ["r2","neg_mean_squared_error"])
#r2:0.83
#neg_mean_squared_error:-12.91
#[0.8305252185799183, -12.906703363062721]
                  
#关闭打印功能试试看？
regassess(reg,Xtrain,Ytrain,cv,scoring = ["r2","neg_mean_squared_error"],show=False) 
#[0.8305252185799183, -12.906703363062721]

#观察一下eta如何影响我们的模型：
from time import time
import datetime
                  
for i in [0,0.2,0.5,1]:
    time0=time()
    reg = XGBR(n_estimators=110,random_state=420,learning_rate=i)
    print("learning_rate = {}".format(i))
    regassess(reg,Xtrain,Ytrain,cv,scoring = ["r2","neg_mean_squared_error"])
    print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
    print("\t")

可以看到为0时，R²的值为负数；而当其为0.2的时候，R²也能有个0.83，这说明只要不是为0，那么至少有个不错的结果。而当为1时，R²就下降了

可以看出，eta对于模型影响不是特别直接，在这个很小的数据集没啥特别影响；

除了运行时间，步长还是一个对模型效果影响巨大的参数，如果设置太大模型就无法收敛（可能导致R²很小或者MSE很大的情况），如果设置太小模型速度就会非常缓慢，但它最后究竟会收敛到何处很难由经验来判定，在训练集上表现出来的模样和在测试集上相差甚远，很难直接探索出一个泛化误差很低的步长：

#调整步长

axisx = np.arange(0.05,1,0.05)
rs = []
te = []
for i in axisx:
    reg = XGBR(n_estimators=110,random_state=420,learning_rate=i)
    score = regassess(reg,Xtrain,Ytrain,cv,scoring = ["r2","neg_mean_squared_error"],show=False)
    test = reg.fit(Xtrain,Ytrain).score(Xtest,Ytest)
    rs.append(score[0])
    te.append(test)
print(axisx[rs.index(max(rs))],max(rs))
plt.figure(figsize=(20,5))
plt.plot(axisx,te,c="gray",label="XGB")
plt.plot(axisx,rs,c="green",label="XGB")
plt.legend()
plt.show()

两者都看不出什么特别的趋势，很难判断eta对于泛化误差的影响

所以通常，我们不调整η，即便调整，一般它也会在[0.01,0.2]之间变动。如果我们希望模型的效果更好，更多的可能是从树本身的角度来说，对树进行剪枝，而不会寄希望于调整η。

• 现在来看，我们的梯度提升树可是说是由三个重要的部分组成：

1. 一个能够衡量集成算法效果的，能够被最优化的损失函数Obj；
2. 一个能够实现预测的弱评估器f_k(x)；
3. 一种能够让弱评估器集成的手段，包括我们讲解的迭代方法，抽样手段，样本加权等等过程；

XGBoost是在梯度提升树的这三个核心要素上运行。并且，XGBoost将原本的梯度提升树拓展开来，让XGBoost不再是单纯的树的集成模型，也不只是单单的回归模型。只要我们调节参数，我们可以选择任何我们希望集成的算法，以及任何我们希望实现的功能；

11.2 XGBoost的智慧

class xgboost.XGBRegressor (kwargs，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')

11.2.1 选择弱评估器：重要参数booster

在XGB中，除了树模型，我们还可以选用线性模型，比如线性回归，来进行集成。基于XGB的这种性质，我们有参数“booster"来控制我们究竟使用怎样的弱评估器：

两个参数都默认为"gbtree"，如果不想使用树模型，则可以自行调整。简单看看现有的数据集上，是什么样的弱评估器表现更好：

#XGBoost
#参数booster

for booster in ["gbtree","gblinear","dart"]:
    reg = XGBR(n_estimators=110
               ,random_state=420
               ,learning_rate=0.1
               ,booster=booster).fit(Xtrain,Ytrain)
    print(booster)
    print(reg.score(Xtest,Ytest))
    
#gbtree
#0.9261386578074982
#gblinear
#0.6028097421971284
#dart
#0.9261386707423716

#自己找线性数据试试看"gblinear"的效果吧~

11.2.2 XGB的目标函数：重要参数objective

集成算法中的损失函数是可选的，要选用什么损失函数取决于我们希望解决什么问题，以及希望使用怎样的模型：

• 回归预测：选择调节后的均方误差RMSE作为我们的损失函数；
• 分类预测：以选择错误率error或者对数损失log_loss；

XGB的目标函数被写作：传统损失函数 + 模型复杂度：

中 代表数据集中的第 个样本， 表示导入第 棵树的数据总量， 代表建立的所有树(n_estimators)。第一项代表传统的损失函数，衡量真实标签 与预测值 之间的差异，通常是RMSE，调节后的均方误差。第二项代表模型的复杂度，使用树模型的某种变换 表示，这个变化代表了一个从树的结构来衡量树模型的复杂度的式子，可以有多种定义

第一项传统损失函数也是与已经建好的所有树相关的，相关在这里：

y^_i 中已经包含了所有树的迭代结果，因此整个目标函数都与K棵树相关

• 泛化误差：用来衡量模型在未知数据上的准确率的指标。一个集成模型(f)在未知数据集(D)上的泛化误差 ，由方差(var)，偏差(bais)和噪声(ε)共同决定，而泛化误差越小，模型就越理想；

从上面的图可以看出来，方差和偏差是此消彼长的，并且模型的复杂度越高，方差越大，偏差越小。方差可以被简单地解释为模型在不同数据集上表现出来地稳定性，而偏差是模型预测的准确度

那方差-偏差困境就可以对应到我们的Obj中了：

第一项是衡量我们的偏差，模型越不准确，第一项就会越大。第二项是衡量我们的方差，模型越复杂，模型的学习就会越具体，到不同数据集上的表现就会差异巨大，方差就会越大

所以我们求解Obj的最小值，其实是在求解方差与偏差的平衡点，以求模型的泛化误差最小，运行速度最快；

在应用中，我们使用参数“objective"来确定我们目标函数的第一部分中的l(y_i,y^_i )，也就是衡量损失的部分：

常用的选择有：

hinge是支持向量机的损失函数；softmax是多分类时常用的损失函数

在xgboost中，我们被允许自定义损失函数，但通常我们还是使用类已经为我们设置好的损失函数。注意：分类型的目标函数导入回归类中会直接报错；

XGB调用方式

在train中导入设置好的树

试试看在其他参数相同的情况下，我们xgboost库本身和sklearn比起来，效果如何：

#重参数objective

#默认reg:linear
reg = XGBR(n_estimators=110,random_state=420).fit(Xtrain,Ytrain)
reg.score(Xtest,Ytest)
#0.9050738056164094
MSE(Ytest,reg.predict(Xtest))
#8.83325117960259

#xgb实现法
import xgboost as xgb
#使用类Dmatrix读取数据
#特征矩阵和标签都要传入
dtrain = xgb.DMatrix(Xtrain,Ytrain)
dtest = xgb.DMatrix(Xtest,Ytest) 

#非常遗憾无法打开来查看，所以通常都是先读到pandas里面查看之后再放到DMatrix中
dtrain
#<xgboost.core.DMatrix at 0x255db93bd00>

#写明参数，silent默认为False，通常需要手动将它关闭
param = {'silent':False,'objective':'reg:linear',"eta":0.1}
num_round = 180 #即n_estimators
#类train，可以直接导入的参数是训练数据，树的数量，其他参数都需要通过params来导入
bst = xgb.train(param, dtrain, num_round) 

#接口predict
from sklearn.metrics import r2_score
r2_score(Ytest,bst.predict(dtest))
#0.9261386578074982
MSE(Ytest,bst.predict(dtest))
#6.873084845395747

随着样本量的逐渐上升，sklearnAPI中调用的结果与xgboost中直接训练的结果会比较相似，但是xgboost会快很多

11.2.3 求解XGB的目标函数

在XGB中我们无法使用梯度下降，原因是XGB的损失函数没有需要求解的参数。我们在传统梯度下降中迭代的是参数，而我们在XGB中迭代的是树，树f_k不是数字组成的向量；

在求解XGB的目标函数的过程中，我们考虑的是如何能够将目标函数转化成更简单的，与树的结构直接相关的写法，以此来建立树的结构与模型的效果（包括泛化能力与运行速度）之间的直接联系。

剩下的看网课+课件吧

11.2.4 参数化决策树f_k(x)：参数alpha，lambda

class 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) *

对于回归树，通常来说每个叶子节点上的预测值是这个叶子节点上所有样本的标签的均值。但值得注意的是，XGB作为普通回归树的改进算法，在y^ 上却有所不同；
对于XGB来说，每个叶子节点上会有一个预测分数（prediction score），也被称为叶子权重。这个叶子权重就是所有在这个叶子节点上的样本在这一棵树上的回归取值，用f_k(x_i)或者w来表示；

当有多棵树的时候，集成模型的回归结果就是所有树的预测分数之和，假设这个集成模型中总共有 棵决策树，则整个模型在这个样本 上给出的预测结果为：

来看正则化系数分别对应的参数：

我们往往认为两种正则化达到的效果是相似的，只不过细节不同。

对于两种正则化如何选择的问题，从XGB的默认参数来看，我们优先选择的是L2正则化。当然，如果想尝试L1也不是不可。两种正则项还可以交互，因此这两个参数的使用其实比较复杂。在实际应用中，正则化参数往往不是我们调参的最优选择；

#求解XGB的目标函数

#使用网格搜索来查找最佳的参数组合
from sklearn.model_selection import GridSearchCV

param = {"reg_alpha":np.arange(0,5,0.05),"reg_lambda":np.arange(0,2,0.05)}

gscv = GridSearchCV(reg,param_grid = param,scoring = "neg_mean_squared_error",cv=cv)

#======【TIME WARNING：10~20 mins】======#
time0=time()
gscv.fit(Xtrain,Ytrain)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))

gscv.best_params_

gscv.best_score_

preds = gscv.predict(Xtest)

from sklearn.metrics import r2_score，mean_squared_error as MSE
r2_score(Ytest,preds)

MSE(Ytest,preds) 
#网格搜索的结果有什么样的含义呢？为什么会出现这样的结果？你相信网格搜索得出的结果吗？试着用数学和你对XGB的理解来解释一下吧

11.2.5 寻找最佳树结构：求解w与T

树我们使用叶子节点上的预测分数来表达，而树的复杂度则是叶子数目加上正则项：

XGBoost默认使用L2正则化

剩下的和11.2.6都是数学原理，先跳过；

11.2.7 让树停止生长：重要参数gamma

从目标函数和结构分数之差Gain的式子中来看，γ是我们每增加一片叶子就会被剪去的惩罚项。增加的叶子越多，结构分数之差Gain会被惩罚越重，所以γ又被称之为是“复杂性控制”（complexity control），所以γ是我们用来防止过拟合的重要参数，还能让树停止生长；

在逻辑回归中，我们使用参数 来设定阈值，并规定如果梯度下降时损失函数减小量小于 下降就会停止。在XGB中，我们规定，只要结构分数之差Gain是大于0的，即只要目标函数还能够继续减小，我们就允许树继续进行分枝——>Gain大于0就说明还可以继续优化，那么久继续减少；

如此，我们可以直接通过设定γ的大小来让XGB中的树停止生长。γ因此被定义为，在树的叶节点上进行进一步分枝所需的最小目标函数减少量。γ设定越大，算法就越保守，树的叶子数量就越少，模型的复杂度就越低：

如果我们希望从代码中来观察 的作用，使用sklearn中传统的学习曲线等工具就比较困难了：

#参数gamma

#======【TIME WARNING: 1 min】=======#
axisx = np.arange(0,5,0.05)
rs = []
var = []
ge = []
for i in axisx:
    reg = XGBR(n_estimators=110,random_state=420,gamma=i)
    result = CVS(reg,Xtrain,Ytrain,cv=cv)
    rs.append(result.mean())
    var.append(result.var())
    ge.append((1 - result.mean())**2+result.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.1
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()

可以看出线段是完全在波动的，没有趋于平衡的趋势

在sklearn下XGBoost太不稳定，如果这样来调整参数的话，效果就很难保证。因此，为了调整 ，我们需要来引入新的工具，xgboost库中的类xgboost.cv：

xgboost.cv (params, dtrain, num_boost_round=10, nfold=3, stratified=False, folds=None, metrics=(), obj=None,feval=None, maximize=False, early_stopping_rounds=None, fpreproc=None, as_pandas=True, verbose_eval=None,show_stdv=True, seed=0, callbacks=None, shuffle=True)

这个工具会同时返回训练集和测试集上的结果

#XGB中自带的交叉验证工具xgboost.cv

import xgboost as xgb

#为了便捷，使用全数据
dfull = xgb.DMatrix(X,y) 

#设定参数
param1 = {'silent':True,'obj':'reg:linear',"gamma":0}
num_round = 110
n_fold=5 

#使用类xgb.cv
time0 = time()
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:393787

#看看类xgb.cv生成了什么结果？
#可以看到随着树不断增加，模型的效果如何变化
cvresult1

这四列分别是：训练样本的调节后的均值/方差、测试样本的调节后的均值/方差；

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,111),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,111),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.legend()
plt.show()

#xgboost中回归模型的默认模型评估指标是什么？
#是调整后的均方误差

可以看到，随着数量的不断增多，学习曲线先是下降然后趋于平稳，平衡后的差值空间就是过拟合的空间；

来看用于回归和分类的评估指标都有哪些：

param1 = {'silent':True,'obj':'reg:linear',"gamma":0,"eval_metric":"mae"}
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,111),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,111),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.legend()
plt.show()

#从这个图中，我们可以看出什么？
#怎样从图中观察模型的泛化能力？
#从这个图的角度来说，模型的调参目标是什么？

图形看起来大致和之前的图一样，但是纵坐标变小了

1. 首先可以看出，这个模型仍然处于过拟合的状态，这说明树的数量还是太多了（也就是说该曲线可以帮助我们找到最佳的n_estimators；
2. 泛化能力即过拟合程度（两条线地差距）；
3. 调参目标就是让两条线尽可能地靠近；

来看看如果我们调整γ，会发生怎样的变化：

param1 = {'silent':True,'obj':'reg:linear',"gamma":0}
param2 = {'silent':True,'obj':'reg:linear',"gamma":20}
num_round = 110
n_fold=5

time0 = time()
cvresult1 = xgb.cv(param1, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:396075

time0 = time()
cvresult2 = xgb.cv(param2, dfull, num_round,n_fold)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:412829

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,181),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,181),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.plot(range(1,181),cvresult2.iloc[:,0],c="green",label="train,gamma=20")
plt.plot(range(1,181),cvresult2.iloc[:,2],c="blue",label="test,gamma=20")
plt.legend()
plt.show()

#从这里，你看出gamma是如何控制过拟合了吗？
#通过控制训练集上的训练，即降低训练集上的表现来控制。但他不一定能影响测试集的表现

γ调大后时间变长说明这是个后剪枝

可以看到，设置γ为20后确实控制了模型的过拟合程度；但是训练集集的方差就要大很多

试一个分类的例子：

#分类的例子

from sklearn.datasets import load_breast_cancer
data2 = load_breast_cancer()

x2 = data2.data
y2 = data2.target

dfull2 = xgb.DMatrix(x2,y2)

param1 = {'silent':True,'obj':'binary:logistic',"gamma":0,"nfold":5}
param2 = {'silent':True,'obj':'binary:logistic',"gamma":2,"nfold":5}
num_round = 100

time0 = time()
#返回error和log_loss区别不大
cvresult1 = xgb.cv(param1, dfull2, num_round,metrics=("error"))
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:164005

time0 = time()
cvresult2 = xgb.cv(param2, dfull2, num_round,metrics=("error")) 
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))
#00:00:228206

plt.figure(figsize=(20,5))
plt.grid()
plt.plot(range(1,101),cvresult1.iloc[:,0],c="red",label="train,gamma=0")
plt.plot(range(1,101),cvresult1.iloc[:,2],c="orange",label="test,gamma=0")
plt.plot(range(1,101),cvresult2.iloc[:,0],c="green",label="train,gamma=2")
plt.plot(range(1,101),cvresult2.iloc[:,2],c="blue",label="test,gamma=2")
plt.legend()
plt.show()

有了xgboost.cv这个工具，我们的参数调整就容易多了。这个工具可以让我们直接看到参数如何影响了模型的泛化能力。这个目标函数及这个目标函数所衍生出来的各种数学过程是XGB原理的重中之重，大部分XGB中基于原理的参数都集中在这个模块之中，到这里大家应该已经基本掌握；

11.3 XGBoost应用中的其他问题

11.3.1 过拟合：剪枝参数与回归模型调参

class 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中，最大深度的功能与参数γ相似，因此如果先调节了 ，则最大深度可能无法展示出巨大的效果。当然，如果先调整了最大深度，则 也有可能无法显示明显的效果。通常来说，这两个参数中我们只使用一个，不过两个都试试也没有坏处。
• 三个随机抽样特征的参数中，前两个比较常用，且有较大影响。但是在XGBoost之前，这种方法并没有被使用到boosting算法当中过。Boosting算法一直以抽取样本（横向抽样）来调整模型过拟合的程度，而实践证明其实纵向抽样（抽取特征）更能够防止过拟合；
• 参数min_child_weight不太常用，它是一篇叶子上的二阶导数h_i之和。本质上来说，这个参数其实是在控制叶子上所需的最小样本量，但是效果没有直接调整深度/叶子来的便捷；

通常当我们获得了一个数据集后，我们先使用网格搜索找出比较合适的n_estimators和eta组合，然后使用gamma或者max_depth观察模型处于什么样的状态（过拟合还是欠拟合，处于方差-偏差图像的左边还是右边？），最后再决定是否要进行剪枝；

们先从最原始的，设定默认参数开始，先观察一下默认参数下，我们的交叉验证曲线长什么样：

#XGBoost应用中的其他问题
#过拟合：剪枝参数与回归模型调参

dfull = xgb.DMatrix(X,y)

param1 = {'silent':True
          ,'obj':'reg:linear'
          ,"subsample":1
          ,"max_depth":6
          ,"eta":0.3
          ,"gamma":0
          ,"lambda":1
          ,"alpha":0
          ,"colsample_bytree":1
          ,"colsample_bylevel":1
          ,"colsample_bynode":1
          ,"nfold":5}
num_round = 200

time0 = time()
cvresult1 = xgb.cv(param1, dfull, num_round)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))

fig,ax = plt.subplots(1,figsize=(15,10))
#ax.set_ylim(top=5)
ax.grid()
ax.plot(range(1,201),cvresult1.iloc[:,0],c="red",label="train,original")
ax.plot(range(1,201),cvresult1.iloc[:,2],c="orange",label="test,original")
ax.legend(fontsize="xx-large")
plt.show()

从曲线上可以看出，模型现在处于过拟合的状态。我们决定要进行剪枝。在这里，我们要使用三组曲线。一组用于展示原始数据上的结果，一组用于展示上一个参数调节完毕后的结果，最后一组用于展示现在我们在调节的参数的结果：

fig,ax = plt.subplots(1,figsize=(15,8))
ax.set_ylim(top=5)
ax.grid()
ax.plot(range(1,201),cvresult1.iloc[:,0],c="red",label="train,original")
ax.plot(range(1,201),cvresult1.iloc[:,2],c="orange",label="test,original")

#调参结果1
param2 = {'silent':True
          ,'obj':'reg:linear'
          ,"subsample":1
          ,"eta":0.05
          ,"gamma":20
          ,"lambda":3.5
          ,"alpha":0.2
          ,"max_depth":4
          ,"colsample_bytree":0.4
          ,"colsample_bylevel":0.6
          ,"colsample_bynode":1
          ,"nfold":5} 
#调参结果2
param3 = {'silent':True
          ,'obj':'reg:linear'
          ,"max_depth":2
          ,"eta":0.05
          ,"gamma":0
          ,"lambda":1
          ,"alpha":0
          ,"colsample_bytree":1
          ,"colsample_bylevel":0.4
          ,"colsample_bynode":1
          ,"nfold":5}

time0 = time()
cvresult2 = xgb.cv(param2, dfull, num_round)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))

time0 = time()
cvresult3 = xgb.cv(param3, dfull, num_round)
print(datetime.datetime.fromtimestamp(time()-time0).strftime("%M:%S:%f"))


ax.plot(range(1,201),cvresult2.iloc[:,0],c="green",label="train,last")
ax.plot(range(1,201),cvresult2.iloc[:,2],c="blue",label="test,last")
ax.plot(range(1,201),cvresult3.iloc[:,0],c="gray",label="train,this")
ax.plot(range(1,201),cvresult3.iloc[:,2],c="pink",label="test,this")
ax.legend(fontsize="xx-large")#调节字体大小
plt.show()

在max_depth、gamma、eta确定下来后，剩余的参数对于模型的影响很小
按param1的顺序逐个进行调参，最终确定了最优参数

先调gamma/最大深度对于正则化项有着巨大影响，不过正则化那两个一般都是放到最后调整。优先调整最大深度/gamma（顺序随意）

• 问题：
• 1. 一个个参数调整太麻烦，可不可以使用网格搜索呢？
     没问题，只是使用的时候要注意，首先XGB的参数非常多，参数可取的范围也很广，究竟是使用np.linspace或者np.arange作为参数的备选值也会影响结果，而且网格搜索的运行速度往往不容乐观(还不可解释)，因此建议至少先使用xgboost.cv来确认参数的范围，否则很可能花很长的时间做了无用功。并且，在使用网格搜索的时候，最好不要一次性将所有的参数都放入进行搜索，最多一次两三个。有一些互相影响的参数需要放在一起使用，比如学习率eta和树的数量n_estimators；
• 2. 调参的时候参数的顺序会影响调参结果吗？
     会影响，因此在现实中，我们会优先调整那些对模型影响巨大的参数。在这里，我建议的剪枝上的调参顺序是：n_estimators与eta共同调节，gamma或者max_depth，采样和抽样参数（纵向抽样影响更大），最后才是正则化的两个参数。当然，可以根据自己的需求来进行调整；
• 3. 调参之后测试集上的效果还没有原始设定上的效果好怎么办？

11.3.2 XGBoost模型的保存和调用

在使用Python进行编程时，我们可能会需要编写较为复杂的程序或者建立复杂的模型。模
型的参数复杂繁多，并且调参过程不是太容易，一旦训练完毕，我们往往希望将训练完毕后的模型保存下来，以便日后用于新的数据集：

(1) 使用Pickle保存和调用模型

pickle是python编程中比较标准的一个保存和调用模型的库，我们可以使用pickle和open函数的连用，来将我们的模型保存到本地：

#XGBoost模型的保存和调用
#1.使用Pickle保存和调用模型

import pickle

dtrain = xgb.DMatrix(Xtrain,Ytrain)

param2 = {'silent':True
          ,'obj':'reg:linear'
          ,"subsample":1
          ,"eta":0.05
          ,"gamma":20
          ,"lambda":3.5
          ,"alpha":0.2
          ,"max_depth":4
          ,"colsample_bytree":0.4
          ,"colsample_bylevel":0.6
          ,"colsample_bynode":1
          ,"nfold":5} 
num_round = 200

bst = xgb.train(param, dtrain, num_round) 

#保存模型
pickle.dump(bst, open("xgboostonboston.dat","wb"))
#注意，open中我们往往使用w或者r作为读取的模式，但其实w与r只能用于文本文件，当我们希望导入的不是文本文件，而是模型本身的时候，我们使用"wb"和"rb"作为读取的模式。其中wb表示以二进制写入，rb表示以二进制读入

#看看模型被保存到了哪里？
import sys
sys.path

#重新打开jupyter
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split as TTS
from sklearn.metrics import mean_squared_error as MSE
import pickle
import xgboost as xgb
data = load_boston()

X = data.data
y = data.target

Xtrain,Xtest,Ytrain,Ytest = TTS(X,y,test_size=0.3,random_state=420) 

#注意，如果我们保存的模型是xgboost库中建立的模型，则导入的数据类型也必须是xgboost库中的数据类型
dtest = xgb.DMatrix(Xtest,Ytest) 

#导入模型
loaded_model = pickle.load(open("xgboostonboston.dat", "rb"))
print("Loaded model from: xgboostonboston.dat") 

#做预测
ypreds = loaded_model.predict(dtest)

from sklearn.metrics import mean_squared_error as MSE, r2_score
MSE(Ytest,ypreds)
#9.18818551539731

r2_score(Ytest,ypreds)
#0.9012595174151566

(2) 使用Joblib保存和调用模型

Joblib是SciPy生态系统中的一部分，它为Python提供保存和调用管道和对象的功能，处理NumPy结构的数据尤其高效，对于很大的数据集和巨大的模型非常有用

#(2) 使用Joblib保存和调用模型

bst = xgb.train(param, dtrain, num_round)

import joblib

#同样可以看看模型被保存到了哪里
joblib.dump(bst,"xgboost-boston.dat")

#不用open
loaded_model = joblib.load("xgboost-boston.dat")

ypreds = loaded_model.predict(dtest)
MSE(Ytest, ypreds)
#9.18818551539731

r2_score(Ytest,ypreds) 
#0.9012595174151566

#使用sklearn中的模型
from xgboost import XGBRegressor as XGBR

bst = XGBR(n_estimators=200
           ,eta=0.05,gamma=20
           ,reg_lambda=3.5
           ,reg_alpha=0.2
           ,max_depth=4
           ,colsample_bytree=0.4
           ,colsample_bylevel=0.6).fit(Xtrain,Ytrain)

joblib.dump(bst,"xgboost-boston.dat")
loaded_model = joblib.load("xgboost-boston.dat") 

#则这里可以直接导入Xtest
ypreds = loaded_model.predict(Xtest)
MSE(Ytest, ypreds)
#9.18818551539731

r2_score(Ytest,ypreds) 
#0.9012595174151566

注意，模型的保存调用与自写函数的保存调用是两回事，大家要注意区分

11.3.3 分类案例：XGB中的样本不均衡问题

存在分类，就会存在样本不平衡问题带来的影响，XGB中存在着调节样本不平衡的参数scale_pos_weight，这个参数非常类似于之前随机森林和支持向量机中我们都使用到过的class_weight参数，通常我们在参数中输入的是负样本量与正样本量之比：

这里默认正样本比例为1，只用调节负样本比例即可

1. 导库，创建样本不均衡的数据集

#1. 导库，创建样本不均衡的数据集

import numpy as np
import xgboost as xgb
import matplotlib.pyplot as plt
from xgboost import XGBClassifier as XGBC
from sklearn.datasets import make_blobs #自创数据集
from sklearn.model_selection import train_test_split as TTS
from sklearn.metrics import confusion_matrix as cm, recall_score as recall, roc_auc_score as auc

class_1 = 500 #类别1有500个样本
class_2 = 50 #类别2只有50个
centers = [[0.0, 0.0], [2.0, 2.0]] #设定两个类别的中心
clusters_std = [1.5, 0.5] #设定两个类别的方差，通常来说，样本量比较大的类别会更加松散
X, y = make_blobs(n_samples=[class_1, class_2],
                  centers=centers,
                  cluster_std=clusters_std,
                  random_state=0, shuffle=False)

Xtrain, Xtest, Ytrain, Ytest = TTS(X,y,test_size=0.3,random_state=420) 

(y == 1).sum() / y.shape[0]
#0.09090909090909091

2. 在数据集上建模：sklearn模式

#2. 在数据集上建模：sklearn模式

#在sklearn下建模#

clf = XGBC().fit(Xtrain,Ytrain)
ypred = clf.predict(Xtest)
clf.score(Xtest,Ytest)
#0.9272727272727272

cm(Ytest,ypred,labels=[1,0])
#array([[  9,   4],
#       [  8, 144]], dtype=int64)

recall(Ytest,ypred)
#0.6923076923076923

auc(Ytest,clf.predict_proba(Xtest)[:,1])
#0.9701417004048585

#负/正样本比例
clf_ = XGBC(scale_pos_weight=10).fit(Xtrain,Ytrain)
ypred_ = clf_.predict(Xtest)
clf_.score(Xtest,Ytest) #分类的默认评估指标：准确率
#0.9333333333333333

cm(Ytest,ypred_,labels=[1,0])
#array([[ 10,   3],
#       [  8, 144]], dtype=int64)

recall(Ytest,ypred_)
#0.7692307692307693

auc(Ytest,clf_.predict_proba(Xtest)[:,1])
#0.9696356275303644

#随着样本权重逐渐增加，模型的recall,auc和准确率如何变化？
for i in [1,5,10,20,30]:
    clf_ = XGBC(scale_pos_weight=i).fit(Xtrain,Ytrain)
    ypred_ = clf_.predict(Xtest)
    print(i)
    print("\tAccuracy:{}".format(clf_.score(Xtest,Ytest)))
    print("\tRecall:{}".format(recall(Ytest,ypred_)))
    print("\tAUC:{}".format(auc(Ytest,clf_.predict_proba(Xtest)[:,1])))

image.png

3. 在数据集上建模：xgboost模式

#3. 在数据集上建模：xgboost模式

dtrain = xgb.DMatrix(Xtrain,Ytrain)
dtest = xgb.DMatrix(Xtest,Ytest) 

#看看xgboost库自带的predict接口
param= {'silent':True,'objective':'binary:logistic',"eta":0.1,"scale_pos_weight":1}
num_round = 100

bst = xgb.train(param, dtrain, num_round)

preds = bst.predict(dtest) 

#看看preds返回了什么？
preds

返回的是每个样本分类的概率

#自己设定阈值
ypred = preds.copy()
ypred[preds > 0.5] = 1
ypred[ypred != 1] = 0

此时XGBoost就知道是要做分类了

#写明参数
scale_pos_weight = [1,5,10]
names = ["negative vs positive: 1"
         ,"negative vs positive: 5"
         ,"negative vs positive: 10"] 

#导入模型评估指标
from sklearn.metrics import accuracy_score as accuracy, recall_score as recall, roc_auc_score as auc

for name,i in zip(names,scale_pos_weight):
    param= {'silent':True,'objective':'binary:logistic'
            ,"eta":0.1,"scale_pos_weight":i}
    clf = xgb.train(param, dtrain, num_round)
    preds = clf.predict(dtest)
    ypred = preds.copy()
    ypred[preds > 0.5] = 1
    ypred[ypred != 1] = 0
    print(name)
    print("\tAccuracy:{}".format(accuracy(Ytest,ypred)))
    print("\tRecall:{}".format(recall(Ytest,ypred)))
    print("\tAUC:{}".format(auc(Ytest,preds)))

提示说silent这个参数没有，去掉之后就简洁了好多；

#当然我们也可以尝试不同的阈值
for name,i in zip(names,scale_pos_weight):
    for thres in [0.3,0.5,0.7,0.9]:
        param= {'silent':True,'objective':'binary:logistic'
                ,"eta":0.1,"scale_pos_weight":i}
        clf = xgb.train(param, dtrain, num_round)
        preds = clf.predict(dtest)
        ypred = preds.copy()
        ypred[preds > thres] = 1
        ypred[ypred != 1] = 0
        print("{},thresholds:{}".format(name,thres))
        print("\tAccuracy:{}".format(accuracy(Ytest,ypred)))
        print("\tRecall:{}".format(recall(Ytest,ypred)))
        print("\tAUC:{}".format(auc(Ytest,preds)))

当阈值越来越大，recall就不断下降，准确率不断上升；但是一般都是一次调一个，不然最好还是用网格搜索的好。不然看着太乱了；

xgboost除了可以做分类和回归，还有其他的多种功能，在一些需要使用精确概率的领域（比如排序ranking），我们希望能够保持概率原有的模样，而提升模型的效果。这种时候，我们就无法使用scale_pos_weight来帮助我们；这里得再深入了解，以后再看吧

11.3.4 XGBoost类中的其他参数和功能

上面的参数和功能主要覆盖了XGBoost中的梯度提升树的原理以及XGBoost自身所带的一些特性。还有一些其他的参数和用法，是算法实际应用时需要考虑的问题。接下来，我们就来看看这些参数：

集成算法就是梯度提升树的部分；弱评估器就是决策树的部分；

• 更多计算资源：n_jobs：输入整数表示使用的线程，输入-1表示使用计算机全部的计算资源。如果我们的数据量很大，则我们可能需要这个参数来为我们调用更多线程；
• 降低学习难度：base_score：它被叫做全局偏差，在分类问题中，它是我们希望关注的分类的先验概率。对于回归来说，这个分数默认0.5，但其实这个分数在这种情况下并不有效；
• 生成树的随机模式：random_state：如果希望控制这种随机性，可以在random_state参数中输入固定整数。需要注意的是，xgb库和sklearn库中，在random_state参数中输入同一个整数未必表示同一个随机模式，不一定会得到相同的结果，因此导致模型的feature_importances也会不一致；
• 自动处理缺失值：missing：XGBoost被设计成是能够自动处理缺失值的模型。我们可
  以在参数missing中输入一个对象，比如np.nan，或数据的任意取值，表示将所有含有这个对象的数据作为空值处理。常用于处理稀疏矩阵；

11.4 总结

XGBoost的难点在于它是一个集大成的模型。它所涉及到的知识点和模型流程，多半在其他常用的机器学习模型中出现过。还可以继续探索下去，慢慢来吧