【机器学习-监督学习】集成学习与梯度提升决策树

【作者主页】Francek Chen
【专栏介绍】$\lceil$Python机器学习$\rfloor$ 机器学习是一门人工智能的分支学科，通过算法和模型让计算机从数据中学习，进行模型训练和优化，做出预测、分类和决策支持。Python成为机器学习的首选语言，依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。
【GitCode】专栏资源保存在我的GitCode仓库：https://gitcode.com/Morse_Chen/Python_machine_learning。

---

  本文将会首先介绍集成学习的思路以及一些常用的集成学习方法，然后介绍梯度提升决策树模型。在前面的文章中，我们讲解了许多不同的机器学习算法，每个算法都有其独特的优缺点。同时，对于同一个任务，往往有多种算法可以将其解决。例如我们要将平面上的一些点分类，假设算法一和算法二的正确率是75%，算法三是50%。3种算法都已经通过调参达到了其最优表现，已经无法再进一步了。那么，我们是否能通过组合这些算法，得到比75%更高的正确率呢？看上去组合之后，算法三会拖算法一和二的后腿，反而会拉低整体表现，更别说提升了。然而，我们考虑表1中的例子。

表1 通过较差算法组合出更好的算法

样本	算法一	算法二	算法三	多数投票
A	$\surd$	$\surd$	$\times$	$\surd$
B	$\times$	$\surd$	$\surd$	$\surd$
C	$\surd$	$\times$	$\surd$	$\surd$
D	$\surd$	$\surd$	$\times$	$\surd$

  该例中共有4个样本和3个算法，每个样本只有两种可能的类别，$\surd$表示算法的预测结果与真实结果相符，$\times$表示预测结果与真实结果不符。各个算法的准确率也符合我们的假设。对每个样本，我们用多数投票（majority voting）制集成3个算法，即3个算法分别给出预测的分类，再取较多者作为最终预测。令人惊讶的是，虽然3个算法单独最高的正确率仅有75%，它们组合之后竟然达到了100%的准确率。仔细观察可以发现，算法一预测错误的B样本，算法二和三都预测正确了。其他样本也类似，当一个算法在此有缺陷时，另外两个算法就可以将其弥补。这一现象启发我们，我们可以通过将不同算法得到的模型按某些方式进行组合，取长补短，从而获得比任一单个模型都要好的表现，这就是集成学习（ensemble learning）的思想。上面的例子中，我们采用了最简单的投票组合方法。接下来，我们来介绍几种实践中常用的集成学习方法。

一、自举聚合与随机森林

  在机器学习的基本思想一文中，我们曾讲解过交叉验证方法。该方法将数据集随机划分成$k$个部分，进行$k$次独立的训练。每次训练时选其中 $k-1$ 份作为训练集，剩下的1份作为验证集。最后，用模型$k$次训练得到的验证误差的平均值来选择模型，由此降低模型过拟合到某一部分数据集上的风险。既然这一过程中会得到多个模型，将其与集成学习的思想结合起来，就得到了自举聚合（bootstrap aggregation，bagging）算法。

  顾名思义，bagging算法由自举和聚合两部分构成。其中，自举指的是自举采样（bootstrap sampling）。与交叉验证中的无重复均匀划分不同，自举采样为了保证随机性，尽可能降低不同子数据集之间的相关性，采用了允许重复的有放回采样。假设数据集的大小为$N$，每次从中有放回地采出$N$个样本。那么，某个样本没有被采样到的概率为$$P(\mathbf{x}没有被采样到)={\left(1-\frac{1}{N}\right)}^N\approx {\mathrm{e}}^{-1}=36.8\%$$ 也就是说，每次按自举采样法采出的大小与原数据集大小相同的样本集合，平均都包含了$63.2\%$的样本。这一比例用来训练和防止过拟合都是较为合适的。当然，我们也可以根据情况调整采样的样本数量。

  假设我们共进行了$B$次自举采样，得到了$B$个子数据集。接下来，我们分别在第$i$个数据集上独立地训练出模型，记为${\hat{f}}_i(\mathbf{x})$。这时，我们可以用某种集成方法将这些模型组合起来，得到更优的模型。例如，我们可以取所有模型预测的平均值作为新的聚合模型：$$\hat{f}(\mathbf{x})=\frac{1}{B}\sum _{i=1}^B{\hat{f}}_i(\mathbf{x})$$

  在交叉验证中，我们的目标是通过多个模型选出最优的超参数，因此我们为每个模型单独计算损失后再取平均。而在bagging算法中，我们希望将这些模型聚合成表现更好的新模型，因此应当衡量最后集成模型的整体效果。我们既可以像普通的机器学习算法那样，提前将数据集中第一部分划分出来作为测试集，不参与聚合采样；也可以计算bagging算法的out-of-bag（OOB）误差。对于数据集中的每个样本$\mathbf{x}$，我们选择那些训练集中不包含$\mathbf{x}$的模型进行测试，将它们的输出用与集成模型相同的集成方式组合起来，得到$\mathbf{x}$的预测结果，并用该结果与真实值计算误差。最后，再和交叉验证一样，将所有样本的OOB误差取平均，作为模型整体的误差。实践表明，由于自举采样保持了合适的采样比例，用OOB误差进行评估与单独划分测试集进行评估的结果没有显著区别。

  为了分析bagging算法的作用原理，我们来考虑bagging相比于单个模型的提升部分。首先对机器学习的预测误差做一个普适性的偏差-方差分解（bias-variance decomposition）分析。设数据的分布为 $y=f(\mathbf{x})+ϵ$，其中$ϵ$是随机噪声，满足 $E(ϵ)=0$，$\mathrm{Var}(ϵ)={\sigma }_{ϵ}^2$。对于某个机器学习模型$\hat{f}(\mathbf{x})$，其学习具有一定随机因素，但与数据集的噪声$ϵ$独立。那么其在样本$\mathbf{x}$上的期望误差平方损失为
$$\begin{array}{rl}\mathcal{L}(\mathbf{x})& =E((f(\mathbf{x})+ϵ-\hat{f}(\mathbf{x}){)}^2)\\ & =E(((f(\mathbf{x})-\hat{f}(\mathbf{x}))+ϵ{)}^2)\\ & =f^2(\mathbf{x})+E({\hat{f}}^2(\mathbf{x}))-2f(\mathbf{x})E(\hat{f}(\mathbf{x}))+E({ϵ}^2)+2E(f(\mathbf{x})-\hat{f}(\mathbf{x}))E(ϵ)\\ & =f^2(\mathbf{x})-2f(\mathbf{x})E(\hat{f}(\mathbf{x}))+(E(\hat{f}(\mathbf{x})){)}^2+E({\hat{f}}^2(\mathbf{x}))-(E(\hat{f}(\mathbf{x})){)}^2+{\sigma }_{ϵ}^2\\ & =(f(\mathbf{x})-E(\hat{f}(\mathbf{x})){)}^2+(E({\hat{f}}^2(\mathbf{x}))-(E(\hat{f}(\mathbf{x})){)}^2)+{\sigma }_{ϵ}^2\\ & ={\mathrm{Bias}}^2(\hat{f}(\mathbf{x}))+\mathrm{Var}(\hat{f}(\mathbf{x}))+{\sigma }_{ϵ}^2\end{array}$$ 其中，${\mathrm{Bias}}^2(\hat{f}(\mathbf{x}))$表示模型$\hat{f}(\mathbf{x})$与数据去除噪声后的真实分布$f(\mathbf{x})$之间的期望偏差，$\mathrm{Var}(\hat{f}(\mathbf{x}))$表示模型在此处的方差。将这一损失函数应用到bagging算法得到的模型上，设第$i$个模型表示的函数为${\hat{f}}_i(\mathbf{x})$，其在样本$\mathbf{x}$上的期望MSE误差为$${\mathcal{L}}_i(\mathbf{x})={\mathrm{Bias}}^2({\hat{f}}_i(\mathbf{x}))+\mathrm{Var}({\hat{f}}_i(\mathbf{x}))+{\sigma }_{ϵ}^2$$ 由于bagging算法中得到每个模型的过程都是独立且相同的，其期望误差和方差可以认为相同，我们将其分别记为$\mu (\mathbf{x})$和${\sigma }^2(\mathbf{x})$。假设不同模型之间的相关系数为$\rho$，即模型之间的协方差 $\mathrm{Cov}({\hat{f}}_i(\mathbf{x}),{\hat{f}}_j(\mathbf{x}))=\rho \sqrt{\mathrm{Var}({\hat{f}}_i(\mathbf{x}))\mathrm{Var}({\hat{f}}_j(\mathbf{x}))}$。以聚合模型取各个模型的平均值为例，其期望MSE误差中的期望偏差和方差分别为$${\mathrm{Bias}}^2(\hat{f}(\mathbf{x}))=f(\mathbf{x})-E\left(\frac{1}{B}\sum _{i=1}^B{\hat{f}}_i(\mathbf{x})\right)=\frac{1}{B}\sum _{i=1}^B(f(\mathbf{x})-E({\hat{f}}_i(\mathbf{x})))=\frac{1}{B}\sum _{i=1}^B\mu (\mathbf{x})=\mu (\mathbf{x})$$ $$\begin{array}{rl}\mathrm{Var}(\hat{f}(\mathbf{x}))& =\mathrm{Var}\left(\frac{1}{B}\sum _{i=1}^B{\hat{f}}_i(\mathbf{x})\right)\\ & =\frac{1}{B^2}\left(\sum _{i=1}^B\mathrm{Var}({\hat{f}}_i(\mathbf{x}))+\sum _{i=1}^B\sum _{j=1,j\ne i}^B\mathrm{Cov}({\hat{f}}_i(\mathbf{x}),{\hat{f}}_j(\mathbf{x}))\right)\\ & =\frac{1}{B^2}\left(B{\sigma }^2(\mathbf{x})+\rho \sum _{i=1}^B\sum _{j=1,j\ne i}^B\sqrt{\mathrm{Var}({\hat{f}}_i(\mathbf{x}))\mathrm{Var}({\hat{f}}_j(\mathbf{x}))}\right)\\ & =\frac{1}{B^2}\left(B{\sigma }^2(\mathbf{x})+\rho (B^2-B){\sigma }^2(\mathbf{x})\right)\\ & =\frac{1-\rho }{B}{\sigma }^2(\mathbf{x})+\rho {\sigma }^2(\mathbf{x})\end{array}$$

  考虑到不同模型所用的数据集是按相同方式采样的，其相关系数应满足 $0\le \rho \le 1$。从上式中可以看出，聚合模型并不改变单一模型的期望偏差，但是可以缩小模型预测的方差。当bagging采样的单模型数量$B$趋于无穷大时，其方差可以降低为单模型方差$\rho$倍，因为 $0\le \rho \le 1$，所有相对于单模型，方差是降低了。因此，bagging算法对于低偏差、高方差的模型（如神经网络和决策树等）的稳定性有较大提升。

  对于决策树模型，其bagging算法的改进版本又称为随机森林（random forest），如图1所示。除了用bagging方法为每个决策树随机采样进行训练之外，上面的推导说明，bagging在不同数据集上训练出的模型相关性较强的时候，降低方差的能力会被削弱。因此，为了进一步降低模型的相关性关系，在决策树每次分裂节点前，我们都从全部的$M$个特征中采样$m$个特征，只在这$m$个特征中选择最优划分特征。例如图1的第二棵树中，每次我们都先从全部的4个特征中采出2个，再从中进一步筛选。这样，由采样方式引起的模型之间的相关性就被进一步削弱了。通常来说，我们选择 $m=\sqrt{M}$，甚至 $m=1$ 来增加模型的随机性。

图1 随机森林和bagging树的区别在于每个分类节点先会采样部分特征再选择划分特征

  下面，我们来动手实现决策树的bagging算法和随机森林算法，并测试其在分类任务上的表现。简单起见，决策树的部分我们直接采用sklearn中的DecisionTreeClassifier。为了体现出随机森林采样特征的特点，我们用sklearn中的工具生成一个高维的点集用作分类数据集。

from tqdm import tqdm
import numpy as np
from matplotlib import pyplot as plt
from sklearn.datasets import make_classification
from sklearn.tree import DecisionTreeClassifier as DTC
from sklearn.model_selection import train_test_split

# 创建随机数据集
X, y = make_classification(
    n_samples=1000, # 数据集大小
    n_features=16, # 特征数，即数据维度
    n_informative=5, # 有效特征个数
    n_redundant=2, # 冗余特征个数，为有效特征的随机线性组合
    n_classes=2, # 类别数
    flip_y=0.1, # 类别随机的样本个数，该值越大，分类越困难
    random_state=0 # 随机种子
)

print(X.shape)

  接下来，我们实现bagging和随机森林的部分。由于这两个算法只有决策树分裂节点时不同，而数据采样部分相同，我们将其写成一个类，用参数来控制执行哪个算法。

class RandomForest():

    def __init__(self, n_trees=10, max_features='sqrt'):
        # max_features是DTC的参数，表示结点分裂时随机采样的特征个数
        # sqrt代表取全部特征的平方根，None代表取全部特征，log2代表取全部特征的对数
        self.n_trees = n_trees
        self.oob_score = 0
        self.trees = [DTC(max_features=max_features) for _ in range(n_trees)]

    # 用X和y训练模型
    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.n_classes = np.unique(y).shape[0]   
        # 集成模型的预测，累加单个模型预测的分类概率，再取较大值作为最终分类
        ensemble = np.zeros((n_samples, self.n_classes))
            
        for tree in self.trees:
            # 自举采样，该采样允许重复
            idx = np.random.randint(0, n_samples, n_samples)
            # 没有被采到的样本
            unsampled_mask = np.bincount(idx, minlength=n_samples) == 0
            unsampled_idx = np.arange(n_samples)[unsampled_mask]
            # 训练当前决策树
            tree.fit(X[idx], y[idx])
            # 累加决策树对OOB样本的预测
            ensemble[unsampled_idx] += tree.predict_proba(X[unsampled_idx])
        # 计算OOB分数，由于是分类任务，我们用正确率来衡量
        self.oob_score = np.mean(y == np.argmax(ensemble, axis=1))
    
    # 预测类别
    def predict(self, X):
        proba = self.predict_proba(X)
        return np.argmax(proba, axis=1)
    
    def predict_proba(self, X):
        # 取所有决策树预测概率的平均
        ensemble = np.mean([tree.predict_proba(X) for tree in self.trees], axis=0)
        return ensemble
    
    # 计算正确率
    def score(self, X, y):
        return np.mean(y == self.predict(X))

  我们通过调整max_feature参数的值来测试两种算法的表现，其中None代表bagging算法，sqrt代表随机森林算法。最后，我们再将两种算法的训练准确率和OOB分数随其中包含的决策树个数的变化曲线画在一张图中，方便对比两种算法的表现。

# 算法测试与可视化
num_trees = np.arange(1, 101, 5)
np.random.seed(0)
plt.figure()

# bagging算法
oob_score = []
train_score = []
with tqdm(num_trees) as pbar:
    for n_tree in pbar:
        rf = RandomForest(n_trees=n_tree, max_features=None)
        rf.fit(X, y)
        train_score.append(rf.score(X, y))
        oob_score.append(rf.oob_score)
        pbar.set_postfix({
            'n_tree': n_tree, 
            'train_score': train_score[-1], 
            'oob_score': oob_score[-1]
        })
plt.plot(num_trees, train_score, color='blue', label='bagging_train_score')
plt.plot(num_trees, oob_score, color='blue', linestyle='-.', label='bagging_oob_score')

# 随机森林算法
oob_score = []
train_score = []
with tqdm(num_trees) as pbar:
    for n_tree in pbar:
        rf = RandomForest(n_trees=n_tree, max_features='sqrt')
        rf.fit(X, y)
        train_score.append(rf.score(X, y))
        oob_score.append(rf.oob_score)
        pbar.set_postfix({
            'n_tree': n_tree, 
            'train_score': train_score[-1], 
            'oob_score': oob_score[-1]
        })
plt.plot(num_trees, train_score, color='red', linestyle='--', label='random_forest_train_score')
plt.plot(num_trees, oob_score, color='red', linestyle=':', label='random_forest_oob_score')

plt.ylabel('Score')
plt.xlabel('Number of trees')
plt.legend()
plt.show()

  从上面的训练结果可以看出，当数据集比较小的时候，随机森林算法相比于简单的bagging算法，可以更有效地防止过拟合。此外，由于随机森林对特征进行了采样，在选择最优特征进行划分时需要的时间也更少，当包含的决策树数量较多时，其训练时间显著小于bagging算法。读者可以通过调整最开始生成数据集的参数观察两种算法效果的变化，当数据集中的样本个数或者数据维度越来越多时，过拟合现象会被自然削弱，随机森林的优势也随之下降。当然，bagging算法只是一个集成学习的框架，除了决策树之外，还可以应用到神经网络等其他不同的模型上。

  最后，我们用sklearn中的bagging和随机森林算法在同样的数据集上进行测试，与我们上面的结果进行比较，验证实现的正确性。

from sklearn.ensemble import BaggingClassifier, RandomForestClassifier

bc = BaggingClassifier(n_estimators=100, oob_score=True, random_state=0)
bc.fit(X, y)
print('bagging：', bc.oob_score_)

rfc = RandomForestClassifier(n_estimators=100, max_features='sqrt', oob_score=True, random_state=0)
rfc.fit(X, y)
print('随机森林：', rfc.oob_score_)

Scikit-learn库中ensemble模块的BaggingClassifier类可以建立自举聚合模型，其基本使用格和常用参数说明式如下。

sklearn.ensemble.BaggingClassifier(self, base_estimator=None, n_estimators=10, max_samples=1.0, max_features=1.0, bootstrap=True, bootstrap_features=False, oob_score=False, n_jobs=None, random_state=None, verbose=0)

参数名称	参数说明
base_estimator	基础分类器或回归器。默认为None，即使用DecisionTreeClassifier（分类）或DecisionTreeRegressor（回归）。
n_estimators	基础模型的数量。默认为10。
max_samples	用于训练每个基础模型的样本数量或比例。默认为1.0（即使用所有样本）。
max_features	用于训练每个基础模型的特征数量或比例。默认为1.0（即使用所有特征）。
bootstrap	是否使用自助法（bootstrap sampling）。默认为True。
bootstrap_features	是否对特征进行自助法采样。默认为False。
oob_score	是否使用袋外样本来估计泛化误差。默认为False。
n_jobs	并行运行的作业数量。默认为None（即不并行）。
random_state	控制随机性。默认为None。
verbose	控制是否打印详细信息。默认为0（即不打印）。

Scikit-learn库中ensemble模块的RandomForestClassifier类可以建立随机森林模型，其基本使用格和常用参数说明式如下。

sklearn.ensemble.andomForestClassifier(self, n_estimators=100, criterion='gini', max_depth=None, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features='auto', max_leaf_nodes=None, min_impurity_decrease=0.0, bootstrap=True, oob_score=False, n_jobs=None, random_state=None, verbose=0, warm_start=False, class_weight=None, ccp_alpha=0.0)

参数名称	参数说明
n_estimators	树的数量。默认为100。
criterion	用于度量分裂质量的标准。分类器使用"gini"或"entropy"，回归器使用"squared_error"、“absolute_error"或"poisson”。默认为"gini"（分类）或"squared_error"（回归）。
max_depth	树的最大深度。默认为None（即树会扩展到所有叶子节点都纯净或包含少于min_samples_split个样本为止）。
min_samples_split	拆分内部节点所需的最小样本数。默认为2。
min_samples_leaf	叶节点所需的最小样本数。默认为1。
min_weight_fraction_leaf	叶节点所需的最小样本权重的比例。默认为0.0。
max_features	每个树训练时使用的最大特征数。默认为"auto"（即使用sqrt(n_features)）。
max_leaf_nodes	最大叶节点数。默认为None。
min_impurity_decrease	节点分裂的最小不纯度减少。默认为0.0。
bootstrap	是否使用自助法采样。默认为True。
oob_score	是否使用袋外样本来估计泛化误差。默认为False。
n_jobs	并行运行的作业数量。默认为None（即不并行）。
random_state	控制随机性。默认为None。
verbose	控制是否打印详细信息。默认为0（即不打印）。
warm_start	是否使用前一次调用的结果来初始化下一次调用。默认为False。
class_weight	分类任务中的类别权重。默认为None（即所有类别权重相同）。
ccp_alpha	表示最小的剪枝复杂度参数。默认为0.0，表示不进行剪枝，默认情况下，剪枝是不启用的。

二、集成学习器

  Bagging算法虽然适用范围较广，但其底层的基础模型需要是同一种类的，例如都是决策树或都是神经网络。否则，其理论推导中假设的各个模型的平均偏差和方差相同的假设将不再成立，其提升模型稳定性的效果也难以保证。假设我们训练好了$n$个模型，分别为 $f_1(\mathbf{x}),\cdots ,f_n(\mathbf{x})$。这些模型的种类可以不同。为了将它们集成起来，我们可以训练一个新模型 $F(\mathbf{x})=g(f_1(\mathbf{x}),f_2(\mathbf{x}),\cdots ,f_n(\mathbf{x}))$，也即将这$n$个模型的预测结果作为输入，由新模型给出最终的输出。我们通常将底层的模型$f$称为基学习器（base learner），而将$g$称为元学习器（meta learner），图2展示了这一过程。

图2 集成学习器的结构

  事实上，该思路在本文最开始的例子中已经有所体现，其中使用的元学习器的作用是选择 $f_1(\mathbf{x}),\cdots ,f_n(\mathbf{x})$ 中的多数分类。此外，我们还很容易想到用取平均值等简单的方式作为新模型的输出。而对于更一般的情况，如果我们把 $f_1(\mathbf{x}),\cdots ,f_n(\mathbf{x})$ 看作一个$n$维输入，那么，我们的目标其实是寻找该输入与样本$\mathbf{x}$对应的输出$y$之间的函数$g$。因此，我们可以使用任意的模型来拟合该函数，如决策树或神经网络等模型，从而提升元学习器的表达能力。

  在由基学习器构建元学习器的训练数据集时，如果让基学习器直接在原本的训练集上进行预测，那么有极大可能发生过拟合现象，得到的数据质量不高。如果让其在划分出的验证集上预测，那么训练集的部分就无法用来训练元学习器，造成数据浪费。为了尽可能提升数据利用效率，同时保证数据质量，我们通常采用堆垛（stacking）方法，如图3所示。假设目前训练的模型是$f_i$，我们采用与交叉验证中相同的做法，把整个数据集均匀划分为$k$份，让$f_i$分别在其中 $k-1$ 份上训练，再在剩余的一份上进行测试，给出其预测结果。这样，当整个训练过程完成时，数据集中的每一份数据都有$f_i$的预测结果，且进行预测的模型必定没有将其用于训练。

图3 集成学习器的堆垛训练

  对每个基学习器重复这一步骤，我们就得到了完整数据集上的 $f_1(\mathbf{x}),\cdots ,f_n(\mathbf{x})$。对于每个原始样本$\mathbf{x}$，我们构建新样本 $\mathbf{s}=(f_1(\mathbf{x}),\cdots ,f_n(\mathbf{x}){)}^{\mathrm{T}}$，并保留原始的标签$y$作为新样本的标签，将$(\mathbf{s},y)$作为训练集来训练元学习器$g$。对于元学习器的测试集，我们不再像训练集那样做划分，而是将同一个基学习器在数据集的$k$个不同部分训练出的结果取平均值，作为新的训练数据。例如，我们在交叉验证中得到了模型$f_i$的$k$个版本 $f_i^{(1)},\cdots ,f_i^{(k)}$，对测试集中的样本$\mathbf{x}$，我们将 $\begin{array}{r}\frac{1}{k}\sum _{j=1}^k{f}_i^{(j)}(\mathbf{x})\end{array}$ 作为新数据的第$i$维，完整的新数据为$$\mathbf{s}={\left(\frac{1}{k}\sum _{j=1}^k{f}_1^{(j)}(\mathbf{x}),\cdots ,\frac{1}{k}\sum _{j=1}^k{f}_n^{(j)}(\mathbf{x})\right)}^{\mathrm{T}}$$

  此外，有时我们会将原始样本$\mathbf{x}$也拼接在新数据$\mathbf{s}$上，一起作为元学习器的输入。这一做法虽然可以防止信息损失，但也可能引入原始数据中的很多干扰因素。因此，是否进行拼接要根据数据集的特点和各个基础学习器的特点而定。理论上来说，堆垛方法中引入了新的模型和参数，理论上只要元学习器的模型合适，就可以替代其他任何集成学习算法。但是考虑到训练难度等问题，我们通常采用逻辑斯谛回归作为元学习器。

  下面，我们来动手实现堆垛方法。简单起见，我们直接使用sklearn库中的一些模型作为基学习器和元学习器。

from sklearn.model_selection import KFold
from sklearn.base import clone

# 堆垛分类器，继承sklearn中的集成分类器基类EnsembleClassifier
class StackingClassifier():

    def __init__(
        self, 
        classifiers, # 基础分类器
        meta_classifier, # 元分类器
        concat_feature=False, # 是否将原始数据拼接在新数据上
        kfold=5 # K折交叉验证
    ):
        self.classifiers = classifiers
        self.meta_classifier = meta_classifier
        self.concat_feature = concat_feature
        self.kf = KFold(n_splits=kfold)
        # 为了在测试时计算平均，我们需要保留每个分类器
        self.k_fold_classifiers = []
        
    def fit(self, X, y):
        # 用X和y训练基础分类器和元分类器
        n_samples, n_features = X.shape
        self.n_classes = np.unique(y).shape[0]
        
        if self.concat_feature:
            features = X
        else:
            features = np.zeros((n_samples, 0))
        for classifier in self.classifiers:
            self.k_fold_classifiers.append([])
            # 训练每个基础分类器
            predict_proba = np.zeros((n_samples, self.n_classes))
            for train_idx, test_idx in self.kf.split(X):
                # 交叉验证
                clf = clone(classifier)
                clf.fit(X[train_idx], y[train_idx])
                predict_proba[test_idx] = clf.predict_proba(X[test_idx])
                self.k_fold_classifiers[-1].append(clf)
            features = np.concatenate([features, predict_proba], axis=-1)
        # 训练元分类器
        self.meta_classifier.fit(features, y)
    
    def _get_features(self, X):
        # 计算输入X的特征
        if self.concat_feature:
            features = X
        else:
            features = np.zeros((X.shape[0], 0))
        for k_classifiers in self.k_fold_classifiers:
            k_feat = np.mean([clf.predict_proba(X) for clf in k_classifiers], axis=0)
            features = np.concatenate([features, k_feat], axis=-1)
        return features
    
    def predict(self, X):
        return self.meta_classifier.predict(self._get_features(X))
        
    def score(self, X, y):
        return self.meta_classifier.score(self._get_features(X), y)

  我们选择KNN、随机森林和逻辑斯谛回归作为基础分类器，用逻辑斯谛回归作为元分类器。下面展示了单个基础分类器与最后的元分类器在数据集上的效果，可以看出，将不同分类器组合后得到的集成分类器要由于任何一个单个分类器。同时，将原始数据拼接到新数据上会使数据中的干扰信息增加，降低集成分类器的效果。大家可以进一步尝试非线性的元分类器，如神经网络和决策树，再次观察原始数据拼接到新数据上的模型预测性能的改变。

from sklearn.linear_model import LogisticRegression as LR
from sklearn.ensemble import RandomForestClassifier as RFC
from sklearn.neighbors import KNeighborsClassifier as KNC

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# 基础分类器
rf = RFC(n_estimators=10, max_features='sqrt', random_state=0).fit(X_train, y_train)
knc = KNC().fit(X_train, y_train)
# multi_class='ovr'表示二分类问题
lr = LR(solver='liblinear', multi_class='ovr', random_state=0).fit(X_train, y_train)
print('随机森林：', rf.score(X_test, y_test))
print('KNN：', knc.score(X_test, y_test))
print('逻辑斯谛回归：', lr.score(X_test, y_test))
# 元分类器
meta_lr = LR(solver='liblinear', multi_class='ovr', random_state=0)

sc = StackingClassifier([rf, knc, lr], meta_lr, concat_feature=False)
sc.fit(X_train, y_train)
print('Stacking分类器：', sc.score(X_test, y_test))

# 带原始特征的stacking分类器
sc_concat = StackingClassifier([rf, knc, lr], meta_lr, concat_feature=True)
sc_concat.fit(X_train, y_train)
print('带原始特征的Stacking分类器：', sc_concat.score(X_test, y_test))

Scikit-learn库中ensemble模块的StackingClassifier类可以建立Stacking分类模型，其基本使用格式和常用参数说明如下。

sklearn.ensemble.StackingClassifier(estimators, final_estimator=None, *, cv=None, stack_method='auto', n_jobs=None, passthrough=False, verbose=0)

参数名称	参数说明
estimators	接收list。表示指定初级学习器。无默认值。
final_estimator	接收str。表示将用于组合初级学习器结果的分类器，即次级分类器。默认为None。
cv	接收int或交叉验证生成器。表示用于训练final_estimator时使用的交叉验证策略。默认为None。
stack_method	接收str。表示每个初级学习器所调用的方法，可选auto、predict_proba、decision_function、predict。默认为“auto”。
n_jobs	并行运行的作业数量，-1表示使用所有可用处理器。默认为None（即不并行）。
passthrough	是否将原始输入数据传递给最终分类器。默认为False。
verbose	控制是否打印详细信息。默认为0（即不打印）。

三、提升算法

  提升（boosting）算法是另一种集成学习的框架，其基本思路是利用当前模型的误差来调整训练数据的权重，使下一个模型更多关注目前误差较大的部分。假如我们在某一数据集上训练出了模型$f_1$，并计算了该模型在各个数据点上预测值与真实值的偏差。既然$f_1$对某些样本的预测已经较为准确，而在另一部分样本上的预测偏差还较大，那么下一个模型$f_2$就应当更多地关注偏差较大的部分。为了达到这一目标，我们可以提高这些样本的损失在总损失中的权重，从而让$f_2$对这些样本的预测更加准确。如图4所示，上排方框内每个蓝色的点代表数据集中的一个样本，点的大小代表样本的权重；下排样本的颜色代表模型对该样本的预测偏差，越红代表偏差越大，越绿色代表偏差越小。接下来，我们不断重复这一过程，用$f_2$的损失调整$f_3$训练数据的权重，以此类推，直到模型数量达到我们预设的上限为止。最后我们得到了$n$个模型，记为 $f_1,\cdots ,f_n$，这些模型称为弱学习器，它们在数据集上各有侧重。因此，我们将这些模型的加权平均，使不同弱学习器的优势都能发挥出来，这样，最终得到的强学习器相比于弱学习器有就更好的效果。

图4 提升算法的原理示意

  与上面介绍的bagging算法和堆垛算法不同，提升算法用到的各个模型在训练时是串行的。这是因为前两种框架中的子模型互相独立，其训练互不干扰，可以并行训练。而提升算法则在一系列模型中用前一个模型的偏差来调整下一个模型训练集的权重，因此这些模型的训练必须由前到后依次进行。

  提升算法的框架中有几个关键问题没有给出确定的答案，分别是如何计算偏差、如何通过偏差计算样本权重、以及如何将得到的各个较弱的学习器结合。这些问题的不同解决方式就导出了不同的提升算法，如适应提升（adaptive boosting，AdaBoost）、逻辑提升（logit boosting）和梯度提升（gradient boosting）等。下面，我们来详细讲解AdaBoost和梯度提升算法，以及梯度提升在决策树上的应用——梯度提升决策树（gradient boost decision tree，GBDT）。

（一）适应提升

  我们先来考虑提升框架的基本模型。由于各个弱学习器的组合方式是加权平均，我们可以用加性模型（additive model）来表示强学习器：$$F(\mathbf{x})=\sum _{i=1}^M{\alpha }_i{f}_i(\mathbf{x})$$ 其中，$M$是弱学习器的数量，$f_i(\mathbf{x})$是弱学习器，${\alpha }_i$是权重。采用不同的损失函数$\mathcal{L}(F)$作为优化目标就可以导出不同的算法。对于标签 y ∈ { − 1 , 1 } y\in\{-1,1\} y∈{−1,1} 的二分类问题，罗伯特·夏皮尔（Robert Schapire）和约拉姆·辛格（Yoram Singer）在1998年的论文中提出，引入了一个新的损失函数$$\mathcal{L}(F)=E_{\mathbf{x}}({\mathrm{e}}^{-yF(\mathbf{x})})$$ 来优化$F(\mathbf{x})$，就得到了AdaBoost算法。该损失函数中的期望代表对数据集中的所有样本取平均。这一目标函数的形式与逻辑斯谛回归中的交叉熵十分相似。设样本$\mathbf{x}$的类别为$y$，$\sigma$为逻辑斯谛函数，那么用$F$预测样本$\mathbf{x}$的交叉熵损失可以写为
$$\begin{array}{rl}{\mathcal{L}}_{\mathrm{CE}}(y,\mathbf{x})& =-\mathbb{I}(y=1)\log \sigma (F(\mathbf{x}))-\mathbb{I}(y=-1)\log (1-\sigma (F(\mathbf{x})))\\ & =-\frac{1+y}{2}\log \frac{{\mathrm{e}}^{F(\mathbf{x})}}{1+{\mathrm{e}}^{F(\mathbf{x})}}-\frac{1-y}{2}\log \frac{1}{1+{\mathrm{e}}^{F(\mathbf{x})}}\\ & =-\frac{1+y}{2}F(\mathbf{x})+\frac{1+y}{2}\log (1+{\mathrm{e}}^{F(\mathbf{x})})+\frac{1-y}{2}\log (1+{\mathrm{e}}^{F(\mathbf{x})})\\ & =-\log {\mathrm{e}}^{\frac{1+y}{2}F(\mathbf{x})}+\log (1+{\mathrm{e}}^{F(\mathbf{x})})\\ & =\log \frac{1+{\mathrm{e}}^{F(\mathbf{x})}}{{\mathrm{e}}^{\frac{1+y}{2}F(\mathbf{x})}}=\{\begin{array}{ll}\log (1+{\mathrm{e}}^{F(\mathbf{x})}),& y=-1\\ \log (1+{\mathrm{e}}^{-F(\mathbf{x})}),& y=1\end{array}\\ & =\log (1+{\mathrm{e}}^{-yF(\mathbf{x})})\end{array}$$ 可以看出，上述两个公式只差了一个变换 $\log (1+z)$。

  回到AdaBoost上，我们再来推导如何从分类器的预测$F(\mathbf{x})$得到样本$\mathbf{x}$类别为1的概率。由于总的损失等于整个数据集上每个样本的损失的平均值，我们先考虑对单个样本$\mathbf{x}$的损失$E({\mathrm{e}}^{-yF(\mathbf{x})}|\mathbf{x})$。记$p$为样本$\mathbf{x}$的类别 $y=1$ 的概率，损失函数可以写为$$E({\mathrm{e}}^{-yF(\mathbf{x})}|\mathbf{x})=p{\mathrm{e}}^{-F(\mathbf{x})}+(1-p){\mathrm{e}}^{F(\mathbf{x})}$$

  当损失函数最小时，其对$F(\mathbf{x})$的偏导数应该为零，即$$0=\frac{\partial E({\mathrm{e}}^{-yF(\mathbf{x})}|\mathbf{x})}{\partial F(\mathbf{x})}=-p{\mathrm{e}}^{-F(\mathbf{x})}+(1-p){\mathrm{e}}^{F(\mathbf{x})}$$ 解得$$p=\frac{1}{1+{\mathrm{e}}^{-2F(\mathbf{x})}}=\sigma (2F(\mathbf{x}))$$ 其中，$\sigma$是逻辑斯谛函数。从这里可以看出，AdaBoost与逻辑斯谛回归在预测上只差了一个常系数2，本质上是等价的。

  要优化诸如AdaBoost的加性模型，我们有许多方法。记${\gamma }_i$是模型$f_i$的参数，那么整个强学习器的参数是 { α 1 , γ 1 , ⋯   , α M , γ M } \{\alpha_1,\gamma_1,\cdots,\alpha_M,\gamma_M\} {α1​,γ1​,⋯,αM​,γM​}。同时优化这些参数复杂度通常很高，而像求解支持向量机时用到的SMO算法那样，每次固定除${\alpha }_i$和${\gamma }_i$以外的所有参数，只优化弱学习器$f_i$的参数和权重，在一般情况下并不能保证收敛。因此，我们可以采用前向分步（forward stagewise）算法，向模型中不断添加弱学习器。记 $F_m(\mathbf{x})=\sum _{i=1}^m{\alpha }_i{f}_i(\mathbf{x})$。假设前 $m-1$ 步的优化已经完成，我们得到了模型 $F_{m-1}$ 并将其固定。在第$m$步的模型 $F_m=F_{m-1}+{\alpha }_m{f}_m$ 中，我们只需要优化${\alpha }_m$和${\gamma }_m$就可以了。这样，整个优化问题的复杂度就降低了许多。

  对于二分类问题，我们假设每个弱学习器$f_i(\mathbf{x})$的输出都是具体的类别-1或1。将上面的优化目标函数 $\mathcal{L}(F)=E_{\mathbf{x}}({\mathrm{e}}^{-yF(\mathbf{x})})$ 应用到第$m$步上，就得到$$\mathcal{L}(F_m)=E_{\mathbf{x}}({\mathrm{e}}^{-y{F}_{m-1}(\mathbf{x})}{\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})=E_{\mathbf{x}}(w_{\mathbf{x}}{\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})=E_{\mathbf{w}}({\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})$$ 其中，$E_{\mathbf{w}}$表示以$w_{\mathbf{x}}$为权重取期望。由于我们固定了$F_{m-1}$，$w_{\mathbf{x}}={\mathrm{e}}^{-y{F}_{m-1}(\mathbf{x})}$与当前待优化的${\alpha }_m$和$f_m$都无关。为了求$\mathcal{L}(F_m)$的最小值，我们令其对${\alpha }_m$和$f_m(\mathbf{x})$的偏导数分别等于零。

  对于${\alpha }_m$，有$$0=\frac{\partial E_{\mathbf{w}}({\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})}{\partial {\alpha }_m}=E_{\mathbf{w}}(-y{f}_m(\mathbf{x}){\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})$$ 注意，$y$和$f_m(\mathbf{x})$都属于 { − 1 , 1 } \{-1,1\} {−1,1}，其乘积当 $y=f_m(\mathbf{x})$ 时为1，当 $y\ne f_m(\mathbf{x})$ 时为-1，上式可以化为
$$\begin{array}{rl}0& =E_{\mathbf{w}}(-y{f}_m(\mathbf{x}){\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})})\\ & =E_{\mathbf{w}}(P(y\ne f_m(\mathbf{x})){\mathrm{e}}^{{\alpha }_m}+(1-P(y\ne f_m(\mathbf{x})))(-{\mathrm{e}}^{-{\alpha }_m}))\\ & =\text{err}\cdot {\mathrm{e}}^{{\alpha }_m}-(1-\text{err})\cdot {\mathrm{e}}^{-{\alpha }_m}\\ \Rightarrow {\alpha }_m& =\frac{1}{2}\log \frac{1-\text{err}}{\text{err}}\end{array}$$ 式中，$\text{err}=E_{\mathbf{w}}(\mathbb{I}(y\ne f_m(\mathbf{x})))$ 表示以$\mathbf{w}$加权的分类错误率。权重${\alpha }_m$的图像如图5所示，可以看出，自变量$\text{err}$的取值范围是$(0,1)$，且函数单调递减。分类错误率越小，我们就应当更看重$f_m$的判断，因此赋予其更大的权重。注意，对于二分类问题，错误率高于0.5时，我们只需要将原本的分类反过来，也即是对应负数权重，就可以达到错误率低于0.5的分类器。所以，在图5中，当 $\text{err}>0.5$ 时，分类器$f_m$的权重是负值。

图5 权重的图像

  对于$f_m$，直接计算梯度求解比较困难，我们利用${\mathrm{e}}^t$函数在 $t=1$ 处的二阶泰勒展开 $\begin{array}{r}{\mathrm{e}}^t\approx 1+t+\frac{1}{2}{t}^2\end{array}$ 以及 $f_m^2(\mathbf{x})=1$，得到$\mathcal{L}(F_m)$的近似形式为$$\mathcal{L}(F_m)\approx E_{\mathbf{w}}\left(1-y{\alpha }_m{f}_m(\mathbf{x})+\frac{1}{2}{\alpha }_m^2\right)$$

  假设分类错误率小于0.5，那么 ${\alpha }_m>0$。令损失函数最小，可以得到
$$\begin{array}{rl}{f}_m(\mathbf{x})& =\arg \underset{f}{\min }\mathcal{L}(F_m(\mathbf{x}))\\ & \approx \arg \underset{f}{\min }{E}_{\mathbf{w}}\left(1-y{\alpha }_m{f}_{(}\mathbf{x})+\frac{1}{2}{\alpha }_m^2|\mathbf{x}\right)\\ & =\arg \underset{f}{\min }{E}_{\mathbf{w}}(-y{\alpha }_{m}f(\mathbf{x})|\mathbf{x})\\ & =\arg \underset{f}{\max }{E}_{\mathbf{w}}(yf(\mathbf{x})|\mathbf{x})\end{array}$$

  由于我们考虑的是$f_m$对某个样本的作用，因此上式右端的期望中限定了$\mathbf{x}$。这一结果提示我们，我们在训练$f_m$时，应当为数据集中的样本添加权重 $w_{\mathbf{x}}={\mathrm{e}}^{-y{F}_{m-1}(\mathbf{x})}$。考虑到 f m ( x ) ∈ { − 1 , 1 } f_m(\boldsymbol x)\in\{-1,1\} fm​(x)∈{−1,1}，上式的解为
$$f_m(\mathbf{x})=\{\begin{array}{ll}1,& E_{\mathbf{w}}(y|\mathbf{x})>0\\ -1,& 其他\end{array}$$

  综上所述，AdaBoost算法的过程如下。

  （1）初始化，设数据集大小为$N$，为所有样本赋予相同的权重 $w_{\mathbf{x}}=1/N$。

  （2）开始迭代，$m=1,2,\cdots ,M$

  a. 在加权的数据集上训练弱分类器$f_m$。

  b. 计算分类器的加权误差 $\text{err}=E_{\mathbf{w}}(\mathbb{I}(y\ne f_m(\mathbf{x})))$。

  c. 计算分类器的权重 $\begin{array}{r}{\alpha }_m=\frac{1}{2}\log \frac{1-\text{err}}{\text{err}}\end{array}$。

  d. 更新数据集中的样本权重 $w_{\mathbf{x}}\leftarrow w_{\mathbf{x}}{\mathrm{e}}^{-y{\alpha }_m{f}_m(\mathbf{x})}$。

  （3）迭代完成，得到强学习器 $\begin{array}{r}F(\mathbf{x})=\text{sgn}\left(\sum _{m=1}^M{\alpha }_m{f}_m(\mathbf{x})\right)\end{array}$。

  上述AdaBoost算法由于弱分类器$f_m$的输出只有-1或1两个值，因此又称为离散适应提升（discrete AdaBoost）。当弱分类器$f_m$的输出是连续的实数时，我们可以将其看作是离散的弱分类器乘上分类器权重${\alpha }_m$的结果，仍然可以用上面的AdaBoost算法进行优化和组合。具体来说，只需要将上面的${\alpha }_m{f}_m(\mathbf{x})$合成输出连续值的弱分类器$f_m^{\prime }(\mathbf{x})$，再将第b~d行合并为：

  b’. 在加权的数据集上估计样本类别为1的概率 $p_m(\mathbf{x})=E_{\mathbf{w}}(y=1|\mathbf{x})\in [0,1]$。

  c’. 令弱分类器 $\begin{array}{r}{f}_m^{\prime }(\mathbf{x})=\frac{1}{2}\log \frac{p_m(\mathbf{x})}{1-p_m(\mathbf{x})}\end{array}$。

  这一算法的可扩展性相比于离散AdaBoost有了进一步的提升，称为实适应提升（real AdaBoost）算法。

  这里，我们不再手动实现AdaBoost算法，而是直接使用sklearn库中的工具，比较其与bagging算法和随机森林算法的效果。其中，所有算法的弱分类器都采用深度为1的决策树，即只从根节点开始做一次分裂，称为桩（stump）。由于桩结构简单、训练快速，但自身效果并不好，因此常用来测试集成学习算法的效果。

from sklearn.ensemble import AdaBoostClassifier

# 初始化stump
stump = DTC(max_depth=1, min_samples_leaf=1, random_state=0)

# 弱分类器个数
M = np.arange(1, 101, 5)
bg_score = []
rf_score = []
dsc_ada_score = []
real_ada_score = []
plt.figure()

with tqdm(M) as pbar:
    for m in pbar:
        # bagging算法
        bc = BaggingClassifier(base_estimator=stump, n_estimators=m, random_state=0)
        bc.fit(X_train, y_train)
        bg_score.append(bc.score(X_test, y_test))
        # 随机森林算法
        rfc = RandomForestClassifier(n_estimators=m, max_depth=1, min_samples_leaf=1, random_state=0)
        rfc.fit(X_train, y_train)
        rf_score.append(rfc.score(X_test, y_test))
        # 离散 AdaBoost，SAMME是分步加性模型（stepwise additive model）的缩写
        dsc_adaboost = AdaBoostClassifier(base_estimator=stump, n_estimators=m, algorithm='SAMME', random_state=0)
        dsc_adaboost.fit(X_train, y_train)
        dsc_ada_score.append(dsc_adaboost.score(X_test, y_test))
        # 实 AdaBoost，SAMME.R表示弱分类器输出实数
        real_adaboost = AdaBoostClassifier(base_estimator=stump, n_estimators=m, algorithm='SAMME.R', random_state=0)
        real_adaboost.fit(X_train, y_train)
        real_ada_score.append(real_adaboost.score(X_test, y_test))

# 绘图
plt.plot(M, bg_score, color='blue', label='Bagging')
plt.plot(M, rf_score, color='red', ls='--', label='Random Forest')
plt.plot(M, dsc_ada_score, color='green', ls='-.', label='Discrete AdaBoost')
plt.plot(M, real_ada_score, color='purple', ls=':', label='Real AdaBoost')
plt.xlabel('Number of trees')
plt.ylabel('Test score')
plt.legend()
plt.tight_layout()
plt.savefig('output_26_1.png')
plt.savefig('output_26_1.pdf')
plt.show()

Scikit-learn库中ensemble模块的AdaBoostClassifier类可以建立适应提升决策树模型，其基本使用格式和常用参数说明如下。

sklearn.ensemble.AdaBoostClassifier(base_estimator='DecisionTreeClassifier', n_estimators=50, learning_rate=1.0, algorithm='SAMME.R', random_state=None, loss=deviance, warm_start=False, verbose=0, class_weight=None)

参数名称	参数说明
base_estimator	基础分类器，如决策树。这是AdaBoost构建过程中每次迭代使用的模型。默认为DecisionTreeClassifier。
n_estimators	迭代次数，即训练了多少个基础分类器并加权求平均的过程。默认为50。
learning_rate	学习速率，影响每个弱分类器对最终预测的影响程度。较小的值可以使模型更加稳定，但也可能减慢收敛速度。默认为1.0。
algorithm	有两种算法选择：SAMME和SAMME.R。SAMME.R更适合连续特征，而SAMME仅支持离散标签。默认为SAMME.R。
random_state	随机种子，用于保证结果的可重复性。默认为None。
loss	接收str。表示指定使用的损失函数。deviance表示使用对数损失函数；exponential表示使用指数损失函数，此时模型只能用于二分类问题。默认为deviance。
warm_start	如果设置为True，可以在之前训练的基础上继续训练，加速多次训练。默认为False。
verbose	控制输出信息的详细程度。默认为0。
class_weight	可以平衡类别权重，防止多数类过拟合。默认为None。

小故事
  集成学习的思想来自于1989年哈佛大学的莱斯利·瓦利安特（Leslie Valiant）和迈克尔·卡恩斯（Michael Kearns）提出的问题，其大意是：如果在某一问题上已有一个比随机方法好、但是仍然较弱的学习算法，能否保证一定存在很强的学习算法？这一反直觉的问题在1990年由夏皮尔给出了回答，他证明总可以通过不断让新的弱学习器关注当前模型的弱点，将弱学习器组合起来得到强学习器。1995年，约夫·弗雷德（Yoav Freund）改进了夏皮尔的算法，并在1996年和夏皮尔共同提出了里程碑式的集成学习算法AdaBoost。由AdaBoost开始衍生出的一大类提升算法，目前是机器学习领域集成学习方法的绝对主流。与此同时，利奥·布雷曼（Leo Breiman）在1994年提出了bagging算法，何天琴（Tin Kam Ho）于1995年开创了随机森林，而布雷曼在2001年进一步发展和完善了随机森林算法。
  由于AdaBoost的非凡表现，1996年开始，不断有学者尝试从数学上分析这一算法，给出泛化误差的理论解释和上界。最初，大家认为与支持向量机相同，AdaBoost的关键在于最大化了最小间隔。但是1999年，布雷曼从最小间隔出发，提出的Arc-gv算法却在比AdaBoost最小间隔更大的情况下泛化误差也更大，这一实验结果无疑给了最小间隔理论沉重的打击。2006年，夏皮尔和列夫·雷津（Lev Reyzin）发现了布雷曼实验中的漏洞，得出了与其相反的结论。随后，高尉和周志华在2013年证明，AdaBoost的关键在于使平均间隔最大的同时使间隔方差最小。直到2019年，艾伦·格隆德（Allan Grønlund）等人最终给出了与下界相匹配的泛化误差上界，解决了AdaBoost的理论问题。

（二）梯度提升

  我们可以换一个视角来考虑加性模型的优化过程。考虑连续的回归问题，和实AdaBoost一样将弱学习器的权重与学习器本身合并，记 $F_{m-1}=\sum _{i=1}^{m-1}{f}_i(\mathbf{x})$ 是到第 $m-1$ 步为止的学习器，这些弱学习器已经固定，不再训练。设总损失函数为$\mathcal{L}(F)$，单个样本的损失函数为$l(y,\hat{y})$，其中$y$和$\hat{y}$分别是样本$\mathbf{x}$对应的真实值和模型的预测值。当我们添加第$m$个弱学习器$f_m$后，模型在数据集上的总损失为$$\mathcal{L}(F_m)=E_{\mathbf{x}}\left(l\left(y,\sum _{i=1}^m{f}_i(\mathbf{x})\right)\right)=E_{\mathbf{x}}(l(y,F_{m-1}(\mathbf{x})+f_m(\mathbf{x})))$$

  简单起见，记 ${\hat{y}}_{m-1}=F_{m-1}(\mathbf{x})$。新添加的$f_m$应当使上面的总损失尽可能小，即$$f_m=\arg \underset{f}{\min }{E}_{\mathbf{x}}(l(y,{\hat{y}}_{m-1}+f(\mathbf{x})))$$

  这里，我们以误差平方损失为例来进一步说明。对于误差平方损失，$l$的形式为$$l(y,{\hat{y}}_{m-1}+f(\mathbf{x}))=\frac{1}{2}(y-{\hat{y}}_{m-1}-f(\mathbf{x}){)}^2$$ 如果把 $y-{\hat{y}}_{m-1}$ 看作是样本$\mathbf{x}$的标签，那么$f_m$就相当于在拟合前 $m-1$ 轮训练后模型$F_{m-1}$的残差。而对于更一般的损失函数形式，这一结论并不一定成立。因此，我们应当寻找其他的方式来寻找最优的$f_m(\mathbf{x})$。在梯度下降中我们讲过，沿损失函数梯度的反方向更新参数是让损失函数的值下降最快的办法。仿照这一思想，如果我们将上面的$F(\mathbf{x})$看成函数$\mathcal{L}(F)$的参数，那么新添加的$f_m$就应当沿着$\mathcal{L}$在$F_{m-1}$处的梯度方向。对函数求“梯度”有些不够严谨，但我们这里可以分别求其对每个样本$\mathbf{x}$上预测值$F(\mathbf{x})$的梯度，为$${\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})={\nabla }_{F(\mathbf{x})}l(y,{\hat{y}}_{m-1})$$ 做一步梯度下降得到$F_m$，为$$F_m(\mathbf{x})={\hat{y}}_{m-1}-{\eta }_m{\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})$$ 其中，${\eta }_m>0$ 是第$m$步的学习率。而 $F_m=F_{m-1}+f_m$，于是新的弱学习器$f_m$为$$f_m(\mathbf{x})=-{\eta }_m{\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})$$

  因此，我们可以通过求损失函数在第 $m-1$ 步时的梯度来得到$f_m$。由于这一方法采用了梯度作为优化目标，我们将其称为梯度提升算法。我们也可以从函数近似的角度来理解梯度提升算法。对于损失函数 $\mathcal{L}(F_m(\mathbf{x}))=\mathcal{L}({\hat{y}}_{m-1}+f_m(\mathbf{x}))$，将其在已经固定的${\hat{y}}_{m-1}$处做泰勒展开，得到$$\mathcal{L}({\hat{y}}_{m-1}+f_m(\mathbf{x}))=\mathcal{L}({\hat{y}}_{m-1})+{\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})\cdot f_m(\mathbf{x})+\frac{1}{2}{\nabla }_{F(\mathbf{x})}^2\mathcal{L}({\hat{y}}_{m-1})f_m^2(\mathbf{x})+\cdots$$如果我们保留一阶近似，就得到$$\mathcal{L}({\hat{y}}_{m-1}+f_m(\mathbf{x}))\approx \mathcal{L}({\hat{y}}_{m-1})+{\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})\cdot f_m(\mathbf{x})$$ 这时，如果寻找使 $\mathcal{L}({\hat{y}}_{m-1}+f_m(\mathbf{x}))$ 最小的$f_m(\mathbf{x})$，考虑到在整个数据集上$f_m(\mathbf{x})$和梯度${\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})$组合起来会变成向量，因而$f_m$应当和梯度共线反向的来使其乘积最小，也即$$f_m(\mathbf{x})=-{\eta }_m{\nabla }_{F(\mathbf{x})}\mathcal{L}({\hat{y}}_{m-1})$$ 这样，我们就得到了和上面仿照梯度下降思路写出的相同的求解公式。并且泰勒展开的思路也提示我们，${\eta }_m$不应当太大，否则被我们省略的高阶项的将会显著影响求解的质量。和上面的AdaBoost相比，学习率${\eta }_m$并不是学习器$f_m$的权重，而是手动设置的超参数，同时可以起到防止过拟合的作用。因此，这一过程又称为梯度提升的收缩（shrinkage），其本质上是通过减小每次迭代中对残差的收敛程度，并增加学习树模型的总数量。

  梯度提升算法一般要求各个弱学习器的模型一致，实践中通常与决策树相结合，组成梯度提升决策树算法。除了简单的直接组合之外，GBDT算法还可以针对决策树的特点进行各种优化，其中极限梯度提升（extreme gradient boosting，XGBoost）是应用最广泛的优化之一。相比于普通的GBDT算法，XGBoost首先在损失函数中添加与决策树复杂度有关的正则化约束，防止单个弱学习器发生过拟合现象。设第$m$棵树为$f_m$，其复杂度为$\Omega (f_m)$，那么损失函数为$$\mathcal{L}(F_m)=\sum _{\mathbf{x}}l(y,{\hat{y}}_{m-1}+f_m(\mathbf{x}))+\Omega (f_m)$$ 此式把每个样本的损失函数从期望改为相加，本质上没有区别。此外，由于前 $m-1$ 棵树已经固定，它们的复杂度也不会计入损失函数中。我们依然对$\mathcal{L}(F_m)$的第一个求和项在$F_{m-1}$处泰勒展开，只不过为了更高的精度，我们保留到二阶项
$$\begin{array}{rl}\mathcal{L}(F_m)& =\sum _{\mathbf{x}}l(y,{\hat{y}}_{m-1}+f_m(\mathbf{x}))+\Omega (f_m)\\ & \approx \sum _{\mathbf{x}}\left(l(y,{\hat{y}}_{m-1})+g(\mathbf{x})f_m(\mathbf{x})+\frac{1}{2}h(\mathbf{x})f_m^2(\mathbf{x})\right)+\Omega (f_m)\end{array}$$ 其中，$g(\mathbf{x})$和$h(\mathbf{x})$分别是损失函数$l$在$F_{m-1}$处的一阶和二阶梯度：
$$\begin{array}{r}g(\mathbf{x})={\nabla }_{F(\mathbf{x})}l(y,{\hat{y}}_{m-1})\\ h(\mathbf{x})={\nabla }_{F(\mathbf{x})}^{2}l(y,{\hat{y}}_{m-1})\end{array}$$

  为了进一步简化损失函数，使其可以求解，我们考虑决策树模型$f_m$应当具有怎样的形式。设决策树的叶节点数目为$T$，$\mathbf{w}\in {\mathbb{R}}^T$ 为每个叶节点上的预测值组成的向量，函数$q(\mathbf{x})$表示样本$\mathbf{x}$被决策树分到的叶节点的编号，那么决策树对样本$\mathbf{x}$的预测为 $f({\mathbf{x}}_m)=w_{q(\mathbf{x})}$。仿照决策树一文中的正则化约束，我们希望决策树的叶节点数目和每个叶节点上的预测值都不要太大，否则容易产生过拟合。因此，复杂度导出的正则化约束可以写为$$\Omega (f_m)=\gamma T+\frac{1}{2}\lambda \Vert \mathbf{w}{\Vert }^2$$ 式中，$\gamma$和$\lambda$都是正则化系数。

  我们将$f_m$和复杂度函数的表达式代入损失函数，并用 I j = { x ∣ q ( x ) = j } \mathcal I_j=\{x|q(\boldsymbol x)=j\} Ij​={x∣q(x)=j} 表示被决策树分到叶节点$j$上的样本的集合，就得到
$$\begin{array}{rl}\mathcal{L}(F_m)& \approx \sum _{\mathbf{x}}\left(l(y,{\hat{y}}_{m-1})+g(\mathbf{x})f_m(\mathbf{x})+\frac{1}{2}h(\mathbf{x})f_m^2(\mathbf{x})\right)+\Omega (f_m)\\ & =\sum _{\mathbf{x}}\left(g(\mathbf{x})w_{q(\mathbf{x})}+\frac{1}{2}h(\mathbf{x})w_{q(\mathbf{x})}^2\right)+\gamma T+\frac{1}{2}\lambda \Vert \mathbf{w}{\Vert }^2+C\\ & =\sum _{j=1}^T\left(\left(\sum _{\mathbf{x}\in {\mathcal{I}}_j}g(\mathbf{x})\right)w_j+\frac{1}{2}\left(\sum _{\mathbf{x}\in {\mathcal{I}}_j}h(\mathbf{x})+\lambda \right)w_j^2\right)+\gamma T+C\end{array}$$ 其中，$C$是与$f_m$无关的常数。上式的最后一个等号从对数据集中的样本求和变为对决策树$f_m$的每个叶节点求和，但求和都会遍历数据集中的所有样本，因此没有本质区别。简记 $G_j=\sum _{\mathbf{x}\in {\mathcal{I}}_j}g(\mathbf{x})$，$H_j=\sum _{\mathbf{x}\in {\mathcal{I}}_j}h(\mathbf{x})$，并舍去常数项，上式可以写为$$\mathcal{L}(F_m)=\sum _{j=1}^T\left(G_j{w}_j+\frac{1}{2}(H_j+\lambda )w_j^2\right)+\gamma T$$

  当决策树对样本的分类方式$q(\mathbf{x})$和叶节点个数$T$固定时，也即决策树的内部结构固定时，上式可以变化的参数只有决策树叶节点上的值$\mathbf{w}$。上式关于$\mathbf{w}$是一个简单的二次函数优化问题，我们省略具体过程，直接给出其最优解为$$w_j^{\ast }={\left(\arg \underset{w}{\min }\mathcal{L}(F_m)\right)}_j=-\frac{G_j}{H_j+\lambda }$$ 其对应的损失函数的最小值为$$\mathcal{L}(F_m;\mathbf{w})=-\frac{1}{2}\sum _{j=1}^T\frac{G_j^2}{H_j+\lambda }+\gamma T$$

  上式是在我们固定决策树结构的前提下求出的损失函数的最小值，那么对于某个结构的决策树，无论其叶节点的值怎样优化，其损失函数都无法比上式更低了。因此，上式可以用来评价一个决策树结构的好坏。与决策树一文中的做法相同，在分裂节点前，我们首先计算上式在分裂前后的值，只有当分裂后的最小损失比分裂前的最小损失更低时，我们才会执行分裂操作。XGBoost算法通过如上的步骤设计的损失函数，往往能使决策树的结构比用其他损失函数构造出的更优，因此取得了较好的均衡表现。在实践中，XGBoost还进行了许多算法和工程上的优化，例如随机森林的特征采样、缓存优化、并行优化等，大大提升了它在大规模数据集上的性能。

  本文最后，我们用xgboost库来展示XGBoost算法的表现，并和上文介绍的其他算法进行对比。为了体现各个算法的水平，我们用sklearn中的make_friedman1工具生成一个更复杂的非线性回归数据集进行测试。设输入样本为$\mathbf{x}$，其每一维都在$[0,1]$内，那么其标签$y$为$$y=10\sin (x_1{x}_2\pi )+20(x_3-0.5{)}^2+10x_4+5x_5+ϵ,ϵ\sim \mathcal{N}(0,{\sigma }^2)$$ 式中的$ϵ$是均值为0、方差为${\sigma }^2$的高斯噪声。可以看出，$\mathbf{x}$只有前5个维度参与计算$y$，剩余维度都是干扰信息。这一数据集最初在中由统计学家杰尔姆·弗里德曼（Jerome H. Friedman）提出。结果表明，XGBoost算法无论是预测精度还是运行效率，都比其他用来比较的算法更加优秀。

# 安装并导入xgboost库
!pip install xgboost
import xgboost as xgb
from sklearn.datasets import make_friedman1
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import BaggingRegressor, RandomForestRegressor, StackingRegressor, AdaBoostRegressor

# 生成回归数据集
reg_X, reg_y = make_friedman1(
    n_samples=2000, # 样本数目
    n_features=100, # 特征数目
    noise=0.5, # 噪声的标准差
    random_state=0 # 随机种子
)

# 划分训练集与测试集
reg_X_train, reg_X_test, reg_y_train, reg_y_test = train_test_split(reg_X, reg_y, test_size=0.2, random_state=0)

def rmse(regressor):
    # 计算regressor在测试集上的RMSE
    y_pred = regressor.predict(reg_X_test)
    return np.sqrt(np.mean((y_pred - reg_y_test) ** 2))

# XGBoost回归树
xgbr = xgb.XGBRegressor(
    n_estimators=100, # 弱分类器数目
    max_depth=1, # 决策树最大深度
    learning_rate=0.5, # 学习率
    gamma=0.0, # 对决策树叶结点数目的惩罚系数，当弱分类器为stump时不起作用
    reg_lambda=0.1, # L2正则化系数
    subsample=0.5, # 与随机森林类似，表示采样特征的比例
    objective='reg:squarederror', # MSE损失函数
    eval_metric='rmse', # 用RMSE作为评价指标
    random_state=0 # 随机种子
)

xgbr.fit(reg_X_train, reg_y_train)
print(f'XGBoost：{rmse(xgbr):.3f}')

# KNN回归
knnr = KNeighborsRegressor(n_neighbors=5).fit(reg_X_train, reg_y_train)
print(f'KNN：{rmse(knnr):.3f}')

# 线性回归
lnr = LinearRegression().fit(reg_X_train, reg_y_train)
print(f'线性回归：{rmse(lnr):.3f}')

# bagging
stump_reg = DecisionTreeRegressor(max_depth=1, min_samples_leaf=1, random_state=0)
bcr = BaggingRegressor(base_estimator=stump_reg, n_estimators=100, random_state=0)
bcr.fit(reg_X_train, reg_y_train)
print(f'Bagging：{rmse(bcr):.3f}')

# 随机森林
rfr = RandomForestRegressor(n_estimators=100, max_depth=1, max_features='sqrt', random_state=0)
rfr.fit(reg_X_train, reg_y_train)
print(f'随机森林：{rmse(rfr):.3f}')

# 堆垛，默认元学习器为带L2正则化约束的线性回归
stkr = StackingRegressor(estimators=[
    ('knn', knnr), 
    ('ln', lnr), 
    ('rf', rfr)
])
stkr.fit(reg_X_train, reg_y_train)
print(f'Stacking：{rmse(stkr):.3f}')

# 带有输入特征的堆垛
stkr_pt = StackingRegressor(estimators=[
    ('knn', knnr), 
    ('ln', lnr), 
    ('rf', rfr)
], passthrough=True)
stkr_pt.fit(reg_X_train, reg_y_train)
print(f'带输入特征的Stacking：{rmse(stkr_pt):.3f}')

# AdaBoost，回归型AdaBoost只有连续型，没有离散型
abr = AdaBoostRegressor(base_estimator=stump_reg, n_estimators=100, learning_rate=1.5, loss='square', random_state=0)
abr.fit(reg_X_train, reg_y_train)
print(f'AdaBoost：{rmse(abr):.3f}')

Scikit-learn库中ensemble模块的GradientBoostingClassifier类可以建立梯度提升决策树模型，其基本使用格式和常用参数说明如下。

sklearn.ensemble.GradientBoostingClassifier(loss='deviance', learning_rate=0.1, n_estimators=100, max_depth=None, subsample=1.0, random_state=None, min_samples_split=2, min_samples_leaf=1,max_depth=None)

参数名称	参数说明
loss	接收str。表示指定使用的损失函数。deviance表示使用对数损失函数；exponential表示使用指数损失函数，此时模型只能用于二分类问题。默认为deviance。
learning_rate	接收float。表示每一棵树的学习率，该参数设定的越小，所需要的基础决策树的数量就越多。默认为0.1。
n_estimators	接收int。表示基础决策树数量。默认为100。
max_depth	每个决策树的最大深度。如果设置为None，则树会尽可能深，直到所有的叶子节点都纯。默认为None。
subsample	接收float。表示用于训练基础决策树的子集占样本集的比例。默认为1.0。
random_state	随机种子，用于保证结果的可重复性。默认为None。
min_samples_split	接收int或float。表示每个基础决策树拆分内部结点所需的最小样本数，若是float，则表示拆分所需的最小样本数占样本数的百分比。默认为2。
min_samples_leaf	接收int或float。表示每个基础决策树模型叶节点所包含的最小样本数，若是float，则表示叶节点最小样本数占样本数的百分比。默认为1。