R-集成学习实用指南-二-R 集成学习实用指南（二）第五章. 简洁的增强算法我们所说的简洁的增强算法是什么意思？增强

R 集成学习实用指南（二）

原文：annas-archive.org/md5/80cdff39151999bd78d34d1c7cef3132

译者：飞龙

协议：CC BY-NC-SA 4.0

第五章. 简洁的增强算法

我们所说的简洁的增强算法是什么意思？增强算法（及其变体）可以说是机器学习工具箱中最重要的算法之一。任何数据分析师都需要了解这个算法，并且最终追求更高的准确率不可避免地会推动对增强技术的需求。据报道，在www.kaggle.org论坛上，复杂和大量数据的增强算法运行了数周，而且大多数获奖方案都是基于这个算法。此外，这些算法在现代图形设备机器上运行。

考虑到其重要性，我们将在下面详细研究增强算法。简洁的当然不是增强算法的变体。由于增强算法是非常重要且关键的算法之一，我们将首先陈述该算法并以最基本的方式实现它，这将展示算法的每个步骤是如何运作的。

我们将从自适应增强算法开始——通常被称为AdaBoost算法——并使用非常简单和原始的代码，我们将展示其在分类问题中的应用。该说明是在一个toy数据集上进行的，以便读者可以清楚地跟随步骤。

在下一节中，我们将扩展分类增强算法到回归问题。对于这个问题，增强变体是著名的梯度增强算法。通过一系列具有单个分割的基本决策树，即所谓的树桩，将创建一个有趣的非线性回归问题。我们将通过平方误差损失函数的选择来展示梯度增强算法。对于增强方法，我们将阐明变量重要性计算。本章的倒数第二节将讨论gbm包的细节。结论部分将对垃圾邮件数据集的袋装、随机森林和增强方法进行比较。本章包括以下主题：

通用增强算法
自适应提升
梯度提升
- 基于树桩的梯度提升
- 基于平方误差损失的梯度提升
增强技术中的变量重要性
使用 gbm 包
袋装、随机森林和增强算法的比较

技术要求

在本章中，我们将使用以下库：

rpart
gbm

通用增强算法

在前几章中介绍的基于树的集成方法，Bagging 和 随机森林，是决策树的重要扩展。然而，虽然 Bagging 通过平均多个决策树提供了更大的稳定性，但偏差仍然存在。这种局限性促使 Breiman 在每个分割点采样协变量以生成一组“独立”的树，并为随机森林奠定了基础。随机森林中的树可以像 Bagging 一样并行开发。在多个树上进行平均的想法是为了确保偏差和方差权衡之间的平衡。提升是决策树的第三大扩展，也可能是最有效的一种。它同样基于集成同质基学习器（在这种情况下，是树），就像 Bagging 和随机森林一样。提升算法的设计完全不同。它是一种顺序集成方法，因为在算法的下一个运行中，前一个学习者的残差/误分类点被改进。

通用提升技术包括一系列算法，这些算法将弱学习器转换为最终强学习器。在分类问题的背景下，弱学习器是一种比随机猜测更好的技术。该方法将弱学习算法转换为更好的算法，这就是它被称为提升的原因。提升技术旨在提供一种接近完美分类器的强学习器。

子空间	子空间	子空间	子空间	准确率
R	R	R	Q	0.75
R	R	Q	R	0.75
R	Q	R	R	0.75
Q	R	R	R	0.75

表 1 简单分类器场景

通过一个简单的例子可以更广泛地理解提升的动机。假设从样本空间通用提升算法中抽取大小为 n 的随机样本，作为独立同分布（IID）的样本。随机样本的分布假设为 D。假设在中有 T=4 个分类器，分类器用于真值函数 f。

考虑一个假设场景，其中样本空间通用提升算法在中由四个部分组成，四个分类器按前表所示执行。提升方法背后的思想是以顺序方式改进分类器。也就是说，分类器是一个接一个地组合，而不是同时全部组合。现在，的错误将被纠正到一个新的分布 D'，其中分类器在区域通用提升算法的错误将给予更多权重，导致分布 D''。提升方法将继续对剩余的分类器进行过程，并给出一个总的组合/集成。Zhou（2012）的伪提升算法（见第二章）总结如下：

步骤 0：初始样本分布是 D，并设置 D1 = D
步骤：
- Dt
- Dt+1 = 改进分布 (Dt,)
最后一步：

算法中的两个步骤改进分布和组合输出显然需要可执行的操作。在下一节中，我们将通过清晰的数值说明来开发自适应提升方法。

自适应提升

Schapire 和 Freund 发明了自适应提升方法。Adaboost是这个技术的流行缩写。

通用自适应提升算法如下：

均匀初始化观测权重：
对于 m，分类器 hm 从 1 到 m 对数据进行多次遍历，执行以下任务：
- 使用权重对训练数据拟合分类器 hm。
- 计算每个分类器的误差如下：
- 计算分类器 hm 的投票权重：
- 设置
输出：

简而言之，算法展开如下：

初始时，我们对所有观测值使用均匀权重。
在下一步中，我们计算考虑的每个分类器的加权误差。
分类器（通常是树桩，或单分支的决策树）需要被选择，通常的做法是选择具有最大准确率的分类器。
在改进分布和组合输出的情况下，如果有准确率相同的分类器，则选择任意一个。
接下来，错误分类的观察值被赋予更多的权重，而正确分类的值被降低权重。这里需要记录一个重要点：

注意

在权重更新步骤中，正确分类为观察值的权重之和将等于错误分类的观察值的权重之和。

从计算分类器的误差到权重更新步骤的步骤会重复 M 次，从而获得每个分类器的投票权重。对于任何给定的观察值，我们然后使用 M 个分类器的预测，这些预测根据各自的投票权重加权，并使用算法中指定的符号函数进行预测。

尽管算法可能很简单，但通过玩具数据集来执行自适应提升方法的工作是一个有用的练习。数据和计算方法取自 Jessica Noss 的视频，可在www.youtube.com/watch?v=gmok1h8wG-Q找到。自适应提升算法的说明现在开始。

考虑一个包含五个三元点的玩具数据集：两个解释变量和一个二元输出值。变量和数据可以用自适应提升来总结，这里的数据点有、、、和。数据首先输入到 R 中，然后作为初步步骤进行可视化：

> # ADAPTIVE BOOSTING with a Toy Dataset
> # https://www.youtube.com/watch?v=gmok1h8wG-Q 
> # Jessica Noss
> # The Toy Data
> x1 <- c(1,5,3,1,5)
> x2 <- c(5,5,3,1,1)
> y <- c(1,1,-1,1,1)
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="The TOY Data Depiction")
> text(x1,x2,labels=names(y),pos=1)

自适应提升

图 1：玩具数据集的简单表示

树桩是决策树的一种特殊情况，已在讨论中提到。在这里，我们将使用树桩作为基学习器。简单看一下前面的图可以帮助我们轻松找到准确率高于随机猜测的树桩。

例如，我们可以在自适应提升处放置一个树桩，并将左侧的所有观察值标记为正，右侧的观察值标记为负。在以下程序中，绿色阴影区域的点被树桩预测为正，而红色阴影区域的点被预测为负。同样，我们可以在和处使用额外的树桩。通过symmetry()函数，我们还可以交换相同树桩的预测。因此，我们之前将绿色阴影区域放在自适应提升的左侧，并预测值为正，通过反转顺序，树桩右侧的区域将被标记为正。对负值也进行类似的分类。在和处重复此任务。使用par、plot、text和rect图形函数，我们在以下内容中展示这些基学习器的可视化表示：

> # Visualizing the stump models
> windows(height=200,width=300)
> par(mfrow=c(2,3))
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1<2")
> text(x1,x2,labels=names(y),pos=1)
> plim <- par("usr")
> rect(xleft=2,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "red",col="red",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 2,ytop = plim[4],
+      border = "green",col="green",density=20 )
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1<4")
> text(x1,x2,labels=names(y),pos=1)
> rect(xleft=4,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "red",col="red",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 4,ytop = plim[4],
+      border = "green",col="green",density=20 )
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1<6")
> text(x1,x2,labels=names(y),pos=1)
> rect(xleft=6,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "red",col="red",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 6,ytop = plim[4],
+      border = "green",col="green",density=20 )
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1>2")
> text(x1,x2,labels=names(y),pos=1)
> rect(xleft=2,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "green",col="green",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 2,ytop = plim[4],
+      border = "red",col="red",density=20 )
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1>4")
> text(x1,x2,labels=names(y),pos=1)
> rect(xleft=4,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "green",col="green",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 4,ytop = plim[4],
+      border = "red",col="red",density=20 )
> plot(x1,x2,pch=c("+","+","-","+","+"),cex=2,
+      xlim=c(0,6),ylim=c(0,6),
+      xlab=expression(x[1]),ylab=expression(x[2]),
+      main="Classification with Stump X1>6")
> text(x1,x2,labels=names(y),pos=1)
> rect(xleft=6,ybottom = plim[3],xright = plim[2],ytop = plim[4],
+      border = "green",col="green",density=20 )
> rect(xleft=plim[1],ybottom = plim[3],xright = 6,ytop = plim[4],
+      border = "red",col="red",density=20 )

前面 R 程序的结果如下所示：

自适应提升

图 2：基于 X1 的树桩分类器

注意，在点 2、4 和 6 处，对于变量自适应提升可以获得类似的分类。尽管不需要给出基于的树桩的完整 R 程序，我们只需在以下图表中生成输出。程序可在代码包中找到。在接下来的讨论中，将忽略基于的树桩：

自适应提升

图 3：基于 X2 的树桩分类器

基于自适应提升选择的树桩导致了一些误分类，我们可以看到观察点 P1、P4 和 P3 被正确分类，而 P2 和 P5 被误分类。基于这个树桩的预测可以表示为（1，-1，-1，1，-1）。基于的树桩正确分类了点 P1 和 P4，而 P2、P3 和 P5 被误分类，这里的向量预测为（1，-1，1，1，-1）。这里考虑的六个模型在 R 程序中用 M1、M2、…、M6 表示，根据前面指定的算法，我们有自适应提升。同样，我们还有其他四个树桩的预测，并将它们输入到 R 中，如下所示：

> # The Simple Stump Models
> M1 <- c(1,-1,-1,1,-1)   # M1 = X1<2 predicts 1, else -1
> M2 <- c(1,-1,1,1,-1)    # M2 = X1<4 predicts 1, else -1
> M3 <- c(1,1,1,1,1)      # M3 = X1<6 predicts 1, else -1
> M4 <- c(-1,1,1,-1,1)    # M4 = X1>2 predicts 1, else -1;M4=-1*M1
> M5 <- c(-1,1,-1,-1,1)   # M5 = X1>4 predicts 1, else -1;M5=-1*M2
> M6 <- c(-1,-1,-1,-1,-1) # M6 = X1>6 predicts 1, else -1;M6=-1*M3

使用六个模型M1-M6给出的预测，我们可以将它们与y中的真实标签进行比较，以查看在每个模型中哪些观察值被误分类：

> # Stem Model Errors
> Err_M1 <- M1!=y
> Err_M2 <- M2!=y
> Err_M3 <- M3!=y
> Err_M4 <- M4!=y
> Err_M5 <- M5!=y
> Err_M6 <- M6!=y
> # Their Misclassifications
> rbind(Err_M1,Err_M2,Err_M3,Err_M4,Err_M5,Err_M6)
          P1    P2    P3    P4    P5
Err_M1 FALSE  TRUE FALSE FALSE  TRUE
Err_M2 FALSE  TRUE  TRUE FALSE  TRUE
Err_M3 FALSE FALSE  TRUE FALSE FALSE
Err_M4  TRUE FALSE  TRUE  TRUE FALSE
Err_M5  TRUE FALSE FALSE  TRUE FALSE
Err_M6  TRUE  TRUE FALSE  TRUE  TRUE

因此，TRUE的值表示在名为模型的行中，名为点的列被误分类。初始化权重自适应提升，并在以下 R 代码块中计算每个模型的加权误差：

> # ROUND 1
> # Weighted Error Computation
> weights_R1 <- rep(1/length(y),length(y)) #Initializaing the weights
> Err_R1 <- rbind(Err_M1,Err_M2,Err_M3,Err_M4,Err_M5,Err_M6)%*%
+   weights_R1
> Err_R1 # Error rate
       [,1]
Err_M1  0.4
Err_M2  0.6
Err_M3  0.2/
Err_M4  0.6
Err_M5  0.4
Err_M6  0.8

由于对应模型 3，或自适应提升的误差是最小的，我们首先选择它，并按照以下方式计算分配给它的投票权重：

> # The best classifier error rate
> err_rate_r1 <- min(Err_R1)
> alpha_3 <- 0.5*log((1-err_rate_r1)/err_rate_r1)
> alpha_3
[1] 0.6931472

因此，提升算法步骤表明自适应提升给出了所需的预测：

> alpha_3*M3
[1] 0.6931472 0.6931472 0.6931472 0.6931472 0.6931472
> sign(alpha_3*M3)
[1] 1 1 1 1 1

中心观察点 P3 仍然被错误分类，所以我们继续到下一步。

现在我们需要更新权重自适应提升，对于分类问题，简化形式的规则如下：

自适应提升

因此，我们需要一个函数，它将接受前一次运行的权重、错误率和模型错误分类作为输入，然后返回包含先前公式的更新权重。我们定义这样的函数如下：


> # Weights Update Formula and Function
> Weights_update <- function(weights,error,error_rate){
+   weights_new <- NULL
+   for(i in 1:length(weights)){
+     if(error[i]==FALSE) weights_new[i] <- 0.5*weights[i]/(1-error_rate)
+     if(error[i]==TRUE) weights_new[i] <- 0.5*weights[i]/error_rate
+   }
+   return(weights_new)
+ }

现在，我们将更新权重并计算六个模型中的每个模型的误差：

> # ROUND 2
> # Update the weights and redo the analyses
> weights_R2 <- Weights_update(weights=weights_R1,error=Err_M3,
+                              error_rate=err_rate_r1)
> Err_R2 <- rbind(Err_M1,Err_M2,Err_M3,Err_M4,Err_M5,Err_M6)%*%
+   weights_R2
> Err_R2 # Error rates
       [,1]
Err_M1 0.25
Err_M2 0.75
Err_M3 0.50
Err_M4 0.75
Err_M5 0.25
Err_M6 0.50

在这里，模型 M1 和 M5 使用新的权重具有相等的错误率，我们简单地选择模型 1，计算其投票权重，并基于更新后的模型进行预测：

> err_rate_r2 <- min(Err_R2)
> alpha_1 <- 0.5*log((1-err_rate_r2)/err_rate_r2)
> alpha_1
[1] 0.5493061
> alpha_3*M3+alpha_1*M1
[1] 1.242453 0.143841 0.143841 1.242453 0.143841
> sign(alpha_3*M3+alpha_1*M1)
[1] 1 1 1 1 1

由于点 P3 仍然被错误分类，我们继续迭代并再次应用循环：

> # ROUND 3
> # Update the weights and redo the analyses
> weights_R3 <- Weights_update(weights=weights_R2,error=Err_M1,
+                              error_rate=err_rate_r2)
> Err_R3 <- rbind(Err_M1,Err_M2,Err_M3,Err_M4,Err_M5,Err_M6)%*%
+   weights_R3
> Err_R3 # Error rates
            [,1]
Err_M1 0.5000000
Err_M2 0.8333333
Err_M3 0.3333333
Err_M4 0.5000000
Err_M5 0.1666667
Err_M6 0.6666667
> err_rate_r3 <- min(Err_R3)
> alpha_5 <- 0.5*log((1-err_rate_r3)/err_rate_r3)
> alpha_5
[1] 0.804719
> alpha_3*M3+alpha_1*M1+alpha_5*M5
[1]  0.4377344  0.9485600 -0.6608779  0.4377344  0.9485600
> sign(alpha_3*M3+alpha_1*M1+alpha_5*M5)
[1]  1  1 -1  1  1

现在分类是完美的，经过三次迭代后，我们没有任何错误分类或错误。本节编程的目的是以基本的方式展示自适应提升算法的步骤。在下一节中，我们将探讨 梯度提升 技术。

梯度提升

自适应提升方法不能应用于回归问题，因为它旨在解决分类问题。梯度提升方法可以使用适当的损失函数来解决分类和回归问题。实际上，梯度提升方法的应用不仅限于这两个标准问题。这项技术起源于 Breiman 的一些观察，并由 Friedman（2000）发展成回归问题。在下一节中，我们将对基本代码进行解释，而无需展示算法。在设置清楚之后，我们将在下一个小节中正式陈述针对平方误差损失函数的提升算法，并创建一个新的函数来实现该算法。

下面的图示是标准正弦波函数的描述。很明显，这是一个非线性关系。在不显式使用正弦变换的情况下，我们将看到使用提升算法来学习这个函数。当然，我们需要简单的回归树桩，我们从一个简单的函数 getNode 开始，它将给出我们想要的分割：

梯度提升

图 4：提升算法能否用于非线性正弦数据？

从头开始构建

在上一节中，我们使用了简单的分类树桩。在那个例子中，简单的视觉检查就足以识别树桩，我们很快获得了 12 个分类树桩。对于回归问题，我们首先定义一个getNode函数，这是对 Tattar（2017）第九章中定义的函数的轻微修改。首先设置所需的符号。

假设我们有一对 n 个数据点从头开始构建，我们正在尝试学习之间的关系，其中f的形式对我们来说完全未知。

对于回归树，分割标准相当直接。对于按 x 值分割的数据，我们计算每个分割部分的ys的平均差平方和，然后将它们加起来。分割标准被选为那个 x 值。这最大化了感兴趣变量中的平均差平方和。R 函数getNode实现了这种思考：

> getNode <- function(x,y)	{
+   xu <- sort(unique(x),decreasing=TRUE)
+   ss <- numeric(length(xu)-1)
+   for(i in 1:length(ss))	{
+     partR <- y[x>xu[i]]
+     partL <- y[x<=xu[i]]
+     partRSS <- sum((partR-mean(partR))²)
+     partLSS <- sum((partL-mean(partL))²)
+     ss[i] <- partRSS + partLSS
+   }
+   xnode <- xu[which.min(ss)]
+   minss <- min(ss)
+   pR <- mean(y[x>xnode])
+   pL <- mean(y[x<=xnode])
+   return(list(xnode=xnode,yR=pR,yL=pL))
+ }

getNode函数的第一步是找到x的唯一值，然后按降序排序。对于唯一值，我们通过 for 循环计算平方和。循环的第一步是将数据分割成左右两部分。

对于每个特定唯一值，在每个分区中计算平均差平方和，然后将它们加起来以得到总的残差平方和。

然后，我们获取x的值，它导致最小的残差平方和。在分割区域中的预测是该区域 y 值的平均值。

getNode函数通过返回x的分割值和左右分区的预测值来结束。我们现在可以创建回归树桩。

正弦波数据首先很容易创建，我们允许 x 值在从头开始构建区间内变化。y 值是简单地将正弦函数应用于 x 向量：

> # Can Boosting Learn the Sine Wave!
> x <- seq(0,2*pi,pi/20)
> y <- sin(x)
> windows(height=300,width=100)
> par(mfrow=c(3,1))
> plot(x,y,"l",col="red",main="Oh My Waves!")

前面的显示结果将是图 1。我们继续获取数据的第一次分割，并在图上显示左右分区的平均值。残差将来自正弦波，它们也将放在同一显示中，如下所示：

> first_split <- getNode(x,y)
> first_split
$xnode
[1] 3.141593
$yR
[1] -0.6353102
$yL
[1] 0.6050574

现在，我们的第一个分割点发生在x值为，从头开始构建这里，3.141593。分割点右侧的预测值为-0.6353102，左侧的预测值为0.6050574。预测值使用segments函数在同一显示上绘制：

> segments(x0=min(x),y0=first_split$yL,
+          x1=first_split$xnode,y1=first_split$yL)
> segments(x0=first_split$xnode,y0=first_split$yR,
+          x1=max(x),y1=first_split$yR)

现在，预测很容易获得，简单的ifelse函数有助于计算它们。与正弦波之间的偏差是残差，我们计算第一组残差和summary函数给出了残差值的简要概述：

> yfit1 <- ifelse(x<first_split$xnode,first_split$yL,first_split$yR)
> GBFit <- yfit1
> segments(x0=x,x1=x,y0=y,y1=yfit1)
> first_residuals <- y-yfit1
> summary(first_residuals)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.60506 -0.25570  0.04752  0.03025  0.32629  0.63531

预测的第一步被保存在GBFit对象中，拟合与预测之间的差异在first_residuals向量中找到。这完成了梯度提升算法的第一轮迭代。第一轮迭代的残差将成为第二轮迭代的回归因变量/输出变量。使用getNode函数，我们执行第二轮迭代，这模仿了早期的代码集：

> second_split <- getNode(x,first_residuals)
> plot(x,first_residuals,"l",col="red",main="The Second Wave!")
> segments(x0=min(x),y0=second_split$yL,
+          x1=second_split$xnode,y1=second_split$yL)
> segments(x0=second_split$xnode,y0=second_split$yR,
+          x1=max(x),y1=second_split$yR)
> yfit2 <- ifelse(x<second_split$xnode,second_split$yL,second_split$yR)
> GBFit <- GBFit+yfit2
> segments(x0=x,x1=x,y0=first_residuals,y1=yfit2)
> second_residuals <- first_residuals-yfit2
> summary(second_residuals)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.51678 -0.24187 -0.02064 -0.01264  0.25813  0.56715

这里的一个重要区别是我们通过累加而不是平均来更新预测。请注意，我们正在模拟第一步的残差，因此下一个拟合所解释的残差剩余部分需要累加而不是平均。残差的范围是多少？建议读者将残差值与早期迭代进行比较。对第三次迭代执行类似的扩展：

> third_split <- getNode(x,second_residuals)
> plot(x,second_residuals,"l",col="red",main="The Third Wave!")
> segments(x0=min(x),y0=third_split$yL,
+          x1=third_split$xnode,y1=third_split$yL)
> segments(x0=third_split$xnode,y0=third_split$yR,
+          x1=max(x),y1=third_split$yR)
> yfit3 <- ifelse(x<third_split$xnode,third_split$yL,third_split$yR)
> GBFit <- GBFit+yfit3
> segments(x0=x,x1=x,y0=second_residuals,y1=yfit3)
> third_residuals <- second_residuals-yfit3
> summary(third_residuals)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.47062 -0.27770 -0.03927 -0.01117  0.18196  0.61331

所有的可视化显示在以下图中：

从头开始构建

图 5：梯度提升算法的三次迭代

显然，我们不可能每次都进行详细的执行迭代，循环是非常重要的。代码被保存在一个块中，并执行了 22 次更多迭代。每次迭代的输出在图中展示，我们将它们全部放入一个外部文件，Sine_Wave_25_Iterations.pdf：

> pdf("Sine_Wave_25_Iterations.pdf")
> curr_residuals <- third_residuals
> for(j in 4:25){
+   jth_split <- getNode(x,curr_residuals)
+   plot(x,curr_residuals,"l",col="red",main=paste0(c("The ", j, "th Wave!")))
+   segments(x0=min(x),y0=jth_split$yL,
+            x1=jth_split$xnode,y1=jth_split$yL)
+   segments(x0=jth_split$xnode,y0=jth_split$yR,
+            x1=max(x),y1=jth_split$yR)
+   yfit_next <- ifelse(x<jth_split$xnode,jth_split$yL,jth_split$yR)
+   GBFit <- GBFit+yfit_next
+   segments(x0=x,x1=x,y0=curr_residuals,y1=yfit_next)
+   curr_residuals <- curr_residuals-yfit_next
+ }
> dev.off()
> summary(curr_residuals)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.733811 -0.093432  0.008481 -0.001632  0.085192  0.350122

经过 25 次迭代后，我们在GBFit中有一个整体拟合，我们可以将其与实际的 y 值进行比较，以查看梯度提升算法的表现如何：

> plot(y,GBFit,xlab="True Y",ylab="Gradient Boosting Fit")

从头开始构建

图 6：梯度拟合与实际正弦数据

对于非线性模型来说，拟合是相当好的。这种方法是为了清楚地理解梯度提升算法。在下一小节中讨论和开发了提升算法的更一般形式。

平方误差损失函数

用表示数据，并将迭代/树的数量固定为一个数字 B。选择一个收缩因子和树深度 d。基于平方误差损失函数的梯度提升算法在此简要说明。参见 Efron 和 Hastie（2016）的第 17.2 算法，如下：

初始化残差和梯度提升预测
对于：
- 为数据拟合深度为 d 的回归树
- 获得预测值
- 通过更新提升预测
- 更新残差
返回函数序列

现在，我们将定义一个名为GB_SqEL的函数，该函数将实现由平方误差损失函数驱动的梯度提升算法。该函数必须提供五个参数：y和x将构成数据，depth将指定树的深度（即回归树中的分割数），iter表示迭代次数，shrinkage是因子。GB_SqEL函数的设置如下：

> # Gradiend Boosting Using the Squared-error Loss Function
> GB_SqEL <- function(y,X,depth,iter,shrinkage){
+   curr_res <- y
+   GB_Hat <- data.frame(matrix(0,nrow=length(y),ncol=iter))
+   fit <- y*0
+   for(i in 1:iter){
+     tdf <- cbind(curr_res,X)
+     tpart <- rpart(curr_res~.,data=tdf,maxdepth=depth)
+     gb_tilda <- predict(tpart)
+     gb_hat <- shrinkage*gb_tilda
+     fit <- fit+gb_hat
+     curr_res <- curr_res-gb_hat
+     GB_Hat[,i] <- fit
+   }
+   return(list(GB_Hat = GB_Hat))
+ }

初始化发生在参数设置中，行fit <- y*0。算法的深度参数在行maxdepth=depth中指定，使用rpart函数创建所需深度的树。predict函数在每次迭代时提供平方误差损失函数的值，而fit+gb_hat执行必要的更新。请注意，GB_Hat[,i]包含每个迭代结束时的预测值。

我们将以 Efron 和 Hastie（2016）的例子来说明算法。考虑的数据与 Lu Gerig 的疾病，或肌萎缩侧索硬化症（ALS）有关。数据集包含有关 1,822 名患有 ALS 疾病的人的信息。目标是预测功能评分的进展率dFRS。研究有关于 369 个预测因子/协变量的信息。在这里，我们将使用GB_SqEL函数来拟合梯度提升技术，并随着迭代次数的增加分析均方误差。详细信息和数据可以从web.stanford.edu/~hastie/CASI/data.html的源文件中获得。现在，我们将平方误差损失函数驱动的提升方法付诸实践：

> als <- read.table("../Data/ALS.txt",header=TRUE)
> alst <- als[als$testset==FALSE,-1]
> temp <- GB_SqEL(y=alst$dFRS,X=alst[,-1],depth=4,
+                 iter=500,shrinkage = 0.02)
> MSE_Train <- 0
> for(i in 1:500){
+   MSE_Train[i] <- mean(temp$GB_Hat[,i]-alst$dFRS)²
+ }
> windows(height=100,width=100)
> plot.ts(MSE_Train)

使用read.table函数，我们将代码包中的数据导入到als对象中。数据以.txt格式从源文件中提供。列testset表示观察值是为了训练目的还是为了测试。我们选择了训练观察值，并删除了第一个变量testset，将其存储在对象alst中。对alst对象应用了GB_SqEL函数，并指定了适当的参数。

每次迭代之后，我们计算均方误差并将其存储在GB_Hat中，如前所述。从以下图中我们可以看出，随着迭代次数的增加，均方误差逐渐减小。在这里，算法在接近 200 次迭代后稳定下来：

平方误差损失函数

图 7：梯度提升和迭代均方误差

在下一节中，我们将看到两个强大 R 包的使用。

使用 adabag 和 gbm 包

将提升方法作为集成技术确实非常有效。算法从头开始展示了分类和回归问题。一旦我们清楚地理解了算法，我们就可以使用 R 包来提供未来的结果。有许多包可用于实现提升技术。然而，在本节中，我们将使用两个最受欢迎的包adabag和gbm。首先，我们需要查看这两个函数的选项。名称很明显，adabag实现了自适应提升方法，而gbm处理梯度提升方法。首先，我们查看以下代码中这两个函数可用的选项：

使用 adabag 和 gbm 包

提升和 gbm 函数

公式是常规的参数。在adabag中的mfinal参数和gbm中的n.trees参数允许指定树的数量或迭代次数。提升函数提供了boos选项，这是使用每个观测值的权重在该迭代中抽取的培训集的 bootstrap 样本。梯度提升是一个更通用的算法，能够处理比回归结构更多的内容。它可以用于分类问题。gbm函数中的distribution选项提供了这些选项。同样，在这里可以看到gbm函数提供了许多其他选项。我们既不会承担解释所有这些选项的艰巨任务，也不会将它们应用于复杂的数据集。用于解释和阐述自适应和梯度提升算法的两个数据集将使用boosting和gbm函数继续进行。

需要更改玩具数据集，我们将多次复制它们，以便我们有足够的观测值来运行boosting和gbm函数：

> # The adabag and gbm Packages
> x1 <- c(1,5,3,1,5)
> x1 <- rep(x1,times=10)
> x2 <- c(5,5,3,1,1)
> x2 <- rep(x2,times=10)
> y <- c(1,1,0,1,1)
> y <- rep(y,times=10)
> toy <- data.frame(x1=x1,x2=x2,y=y)
> toy$y <- as.factor(toy$y)
> AB1 <- boosting(y~.,data=toy,boos=TRUE,mfinal = 10,
+                 maxdepth=1,minsplit=1,minbucket=1)
> predict.boosting(AB1,newdata=toy[,1:2])$class
 [1] "1" "1" "0" "1" "1" "1" "1" "0" "1" "1" "1" "1" "0" "1" "1" "1" "1" "0"
[19] "1" "1" "1" "1" "0" "1" "1" "1" "1" "0" "1" "1" "1" "1" "0" "1" "1" "1"
[37] "1" "0" "1" "1" "1" "1" "0" "1" "1" "1" "1" "0" "1" "1"

maxdepth=1函数确保我们只使用树桩作为基础分类器。很容易看出提升函数工作得非常完美，因为所有观测都被正确分类。

与boosting函数一样，我们需要更多的数据点。我们通过seq函数增加这些数据点，并使用distribution="gaussian"选项，要求gbm函数拟合回归提升技术：

> x <- seq(0,2*pi,pi/200)
> y <- sin(x)
> sindata <- data.frame(cbind(x,y))
> sin_gbm <- gbm(y~x,distribution="gaussian",data=sindata,
+                n.trees=250,bag.fraction = 0.8,shrinkage = 0.1)
> par(mfrow=c(1,2))
> plot.ts(sin_gbm$fit, main="The gbm Sine Predictions")
> plot(y,sin_gbm$fit,main="Actual vs gbm Predict")

使用绘图函数，我们比较了梯度提升方法的拟合情况。以下图表表明拟合是适当的。然而，图表也显示故事中有些地方不太对劲。提升方法在使用 adabag 和 gbm 包和处的函数近似留下了很多遗憾，实际与预测的图表表明在 0 处存在不连续/性能不佳的问题。然而，我们不会对这些问题深入探讨：

使用 adabag 和 gbm 包

图 8：使用 gbm 函数进行的正弦波近似

接下来，我们将讨论变量重要性的概念。

变量重要性

提升方法本质上使用树作为基学习器，因此变量重要性的概念在这里与树、袋装和随机森林相同。我们只需像在袋装或随机森林中那样，将变量在树之间的重要性相加。

对于来自adabag包的提升拟合对象，变量重要性提取如下：

> AB1$importance
 x1  x2 
100   0

这意味着提升方法根本未使用x2变量。对于梯度提升对象，变量重要性由summary函数给出：

> summary(sin_gbm)
  var rel.inf
x   x     100

现在很明显，我们只有一个变量，因此解释回归因变量很重要，我们当然不需要某些软件来告诉我们。当然，在复杂情况下是有用的。比较是基于树的集成方法的。让我们继续下一节。

比较袋装、随机森林和提升

在前一章中，我们进行了袋装和随机森林方法的比较。现在，我们使用gbm函数将提升准确率添加到早期分析中：

> data("spam")
> set.seed(12345)
> Train_Test <- sample(c("Train","Test"),nrow(spam),replace = TRUE,
+ prob = c(0.7,0.3))
> head(Train_Test)
[1] "Test"  "Test"  "Test"  "Test"  "Train" "Train"
> spam_Train <- spam[Train_Test=="Train",]
> spam_TestX <- within(spam[Train_Test=="Test",],
+                      rm(type))
> spam_TestY <- spam[Train_Test=="Test","type"]
> spam_Formula <- as.formula("type~.")
> spam_rf <- randomForest(spam_Formula,data=spam_Train,coob=TRUE,
+                         ntree=500,keepX=TRUE,mtry=5)
> spam_rf_predict <- predict(spam_rf,newdata=spam_TestX,type="class")
> rf_accuracy <- sum(spam_rf_predict==spam_TestY)/nrow(spam_TestX)
> rf_accuracy
[1] 0.9436117
> spam_bag <- randomForest(spam_Formula,data=spam_Train,coob=TRUE,
+                          ntree=500,keepX=TRUE,mtry=ncol(spam_TestX))
> spam_bag_predict <- predict(spam_bag,newdata=spam_TestX,type="class")
> bag_accuracy <- sum(spam_bag_predict==spam_TestY)/nrow(spam_TestX)
> bag_accuracy
[1] 0.9350464
> spam_Train2 <- spam_Train
> spam_Train2$type <- ifelse(spam_Train2$type=="spam",1,0)
> spam_gbm <- gbm(spam_Formula,distribution="bernoulli",data=spam_Train2,
+                 n.trees=500,bag.fraction = 0.8,shrinkage = 0.1)
> spam_gbm_predict <- predict(spam_gbm,newdata=spam_TestX,
+                             n.trees=500,type="response")
> spam_gbm_predict_class <- ifelse(spam_gbm_predict>0.5,"spam","nonspam")
> gbm_accuracy <- sum(spam_gbm_predict_class==spam_TestY)/nrow(spam_TestX)
> gbm_accuracy
[1] 0.945753
> summary(spam_gbm)
                                var      rel.inf
charExclamation     charExclamation 21.985502703
charDollar               charDollar 18.665385239
remove                       remove 11.990552362
free                           free  8.191491706
hp                               hp  7.304531600

num415                       num415  0.000000000
direct                       direct  0.000000000
cs                               cs  0.000000000
original                   original  0.000000000
table                         table  0.000000000
charHash                   charHash  0.000000000

提升准确率0.9457高于随机森林准确率0.9436。进一步的微调，将在下一章中探讨，将有助于提高准确率。变量重要性也可以通过summary函数轻松获得。

摘要

提升是决策树的另一种推论。它是一种迭代技术，通过更无顾忌地针对前一次迭代的误差。我们从重要的自适应提升算法开始，并使用非常简单的玩具数据来说明其基础。然后，该方法被扩展到回归问题，我们通过两种不同的方法说明了梯度提升方法。对adabag和gbm这两个包进行了简要阐述，并再次强调了变量重要性的概念。对于垃圾邮件数据集，我们通过提升方法获得了更高的准确率，因此提升算法的讨论特别有用。

本章考虑了提升算法的不同变体。然而，我们没有讨论它为什么有效。在下一章中，这些方面将更详细地介绍。

第六章：提升改进

在上一章中，我们学习了提升算法。我们研究了算法的结构形式，用一个数值示例进行了说明，然后将算法应用于回归和分类问题。在本章简要介绍中，我们将涵盖提升算法及其基础的一些理论方面。提升理论在这里也很重要。

本章，我们还将从几个不同的角度探讨提升算法为什么有效。不同类型的问题需要不同类型的损失函数，以便提升技术能够有效。在下一节中，我们将探讨我们可以选择的不同类型的损失函数。极端梯度提升方法在专门用于处理xgboost包的章节中概述。此外，h2o包将在最后一节中最终讨论，这可能对其他集成方法也很有用。本章将涵盖以下主题：

提升为什么有效？
gbm包
xgboost包
h2o包

技术要求

本章我们将使用以下 R 库：

adabag
gbm
h2o
kernlab
rpart
xgboost

提升为什么有效？

上一章中关于自适应提升算法的部分包含了m个模型、分类器提升为什么有效？、n个观察值和权重，以及一个按顺序确定的投票权。通过一个玩具示例说明了自适应提升方法的适应，然后使用专用函数应用了该方法。与袋装法和随机森林方法相比，我们发现提升提供了最高的准确率，这你可能还记得上一章前面提到的那个部分的结果。然而，算法的实现并没有告诉我们为什么它预期会表现得更好。

我们没有一个普遍接受的答案来解释提升为什么有效，但根据 Berk（2016）的第 6.2.2 小节，有三种可能的解释：

提升是边缘最大化器
提升是统计优化器
提升是一个插值器

但这些实际上意味着什么呢？现在我们将逐一介绍这些观点。在提升算法中，一个观察值的边缘计算如下：

提升为什么有效？

我们可以看到，边界值是正确分类和错误分类总和之间的差异。在先前的公式中，为什么提升工作？表示投票权重。提升算法之所以工作得如此之好，尤其是在分类问题上，是因为它是一个边界最大化器。在他们的开创性论文中，提升算法的发明者 Schapire 等人声称，提升特别擅长找到具有大边界的分类器，因为它专注于那些边界值小（或负）的例子，并迫使基本学习算法为这些例子生成良好的分类。下面引文中的粗体部分将在下一个 R 代码块中展示。

将使用来自kernlab包的垃圾邮件数据集来说明这一关键思想。gbm包中的boosting函数将拟合数据，以区分垃圾邮件和正常邮件。我们将从一个仅包含一次迭代的初始模型开始，获取准确率，获取边界值，然后在下面的代码块中生成摘要：

> data("spam")
> set.seed(12345)
> Train_Test <- sample(c("Train","Test"),nrow(spam),replace = TRUE,prob = c(0.7,0.3))
> spam_Train <- spam[Train_Test=="Train",]
> spam_Formula <- as.formula("type~.")
> spam_b0 <- boosting(spam_Formula,data=spam_Train,mfinal=1)
> sum(predict(spam_b0,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 0.905
> mb0 <- margins(spam_b0,newdata=spam_Train)$margins
> mb0[1:20]
 [1]  1 -1  1 -1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
> summary(mb0)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -1.000   1.000   1.000   0.809   1.000   1.000

接下来，我们将迭代次数增加到 5、10、20、50 和 200。读者应跟踪以下结果的准确率和边界值摘要：

> spam_b1 <- boosting(spam_Formula,data=spam_Train,mfinal=5)
> sum(predict(spam_b1,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 0.948
> mb1 <- margins(spam_b1,newdata=spam_Train)$margins
> mb1[1:20]
 [1]  1.0000 -0.2375  0.7479 -0.2375  0.1771  0.5702  0.6069  0.7479  1.0000
[10]  0.7479  1.0000  1.0000 -0.7479  1.0000  0.7479  1.0000  0.7479  1.0000
[19] -0.0146  1.0000
> summary(mb1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -1.000   0.631   1.000   0.783   1.000   1.000 
> spam_b2 <- boosting(spam_Formula,data=spam_Train,mfinal=10)
> sum(predict(spam_b2,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 0.969		
> mb2 <- margins(spam_b2,newdata=spam_Train)$margins
> mb2[1:20]
 [1]  0.852  0.304  0.245  0.304  0.288  0.629  0.478  0.678  0.827  0.678
[11]  1.000  1.000 -0.272  0.517  0.700  0.517  0.700  0.478  0.529  0.852
> summary(mb2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -0.517   0.529   0.807   0.708   1.000   1.000 
> spam_b3 <- boosting(spam_Formula,data=spam_Train,mfinal=20)
> sum(predict(spam_b3,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 0.996
> mb3 <- margins(spam_b3,newdata=spam_Train)$margins
> mb3[1:20]
 [1] 0.5702 0.3419 0.3212 0.3419 0.3612 0.6665 0.4549 0.7926 0.7687 0.6814
[11] 0.8958 0.5916 0.0729 0.6694 0.6828 0.6694 0.6828 0.6130 0.6813 0.7467
> summary(mb3)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -0.178   0.537   0.719   0.676   0.869   1.000 
> spam_b4<- boosting(spam_Formula,data=spam_Train,mfinal=50)
> sum(predict(spam_b4,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 1
> mb4 <- margins(spam_b4,newdata=spam_Train)$margins
> mb4[1:20]
 [1] 0.407 0.333 0.386 0.333 0.379 0.518 0.486 0.536 0.579 0.647 0.695 0.544
[13] 0.261 0.586 0.426 0.586 0.426 0.488 0.572 0.677
> summary(mb4)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.098   0.444   0.590   0.586   0.729   1.000 
> spam_b5<- boosting(spam_Formula,data=spam_Train,mfinal=200)
> sum(predict(spam_b5,newdata=spam_Train)$class==
+       spam_Train$type)/nrow(spam_Train)
[1] 1
> mb5 <- margins(spam_b5,newdata=spam_Train)$margins
> mb5[1:20]
 [1] 0.386 0.400 0.362 0.368 0.355 0.396 0.368 0.462 0.489 0.491 0.700 0.486
[13] 0.317 0.426 0.393 0.426 0.393 0.385 0.624 0.581
> summary(mb5)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.272   0.387   0.482   0.510   0.607   0.916

第一个重要的区别是，边界值完全摆脱了负数，并且在迭代次数达到 50 或更多之后，每个边界值都是非负的。为了更清楚地了解这一点，我们将边界值进行列绑定，然后查看初始化时具有负边界的所有观察值：

> View(cbind(mb1,mb2,mb3,mb4,mb5)[mb1<0,])

以下快照显示了前面代码的结果：

为什么提升工作？

图 1：迭代过程中误分类观察值的边界

因此，我们可以清楚地看到，随着迭代次数的增加，边界值也在增加。

关于提升，特别是自适应提升，需要记住的第二点是它是一个统计优化器。Berk（2016）第 264-5 页和 Zhou（2012）第 25-6 页显示，提升集成达到了贝叶斯误差率。这意味着由于指数损失被最小化，分类错误率也被最小化。

关于提升作为插值器的第三点非常直接。很明显，提升的迭代可以看作是随机森林的加权平均。

到目前为止，提升方法仅解决了分类和回归问题。损失函数对于机器学习算法至关重要，下一节将讨论各种损失函数，这些函数将有助于为不同格式的数据设置提升算法。

`gbm`包

由 Greg Ridgeway 创建的 R gbm 包是一个非常通用的包。该包的详细信息可以在www.saedsayad.com/docs/gbm2.pdf找到。该文档详细介绍了梯度提升的理论方面，并说明了 gbm 函数的各种其他参数。首先，我们将考虑 gbm 函数中可用的收缩因子。

收缩参数非常重要，也有助于解决过拟合问题。通过此选项实现惩罚。对于垃圾邮件数据集，我们将收缩选项设置为 0.1（非常大）和 0.0001（非常小），并观察性能如何受到影响：

> spam_Train2 <- spam_Train
> spam_Train2$type <- as.numeric(spam_Train2$type)-1
> spam_gbm <- gbm(spam_Formula,distribution="bernoulli",
+       data=spam_Train2, n.trees=500,bag.fraction = 0.8,
+       shrinkage = 0.1)
> plot(spam_gbm) # output suppressed
> summary(spam_gbm)
                                var      rel.inf
charExclamation     charExclamation 21.740302311
charDollar               charDollar 18.505561199
remove                       remove 11.722965305
your                           your  8.282553567
free                           free  8.142952834
hp                               hp  7.399617456

num415                       num415  0.000000000
direct                       direct  0.000000000
cs                               cs  0.000000000
original                   original  0.000000000
table                         table  0.000000000
charSquarebracket charSquarebracket  0.000000000

摘要函数还绘制了相对变量重要性图。如下截图所示：

gbm 包

图 2：相对变量影响图

接下来获取拟合对象的详细信息：

> spam_gbm
gbm(formula = spam_Formula, distribution = "bernoulli", data = spam_Train2, 
    n.trees = 500, shrinkage = 0.1, bag.fraction = 0.8)
A gradient boosted model with bernoulli loss function.
500 iterations were performed.
There were 57 predictors, of which 43 had nonzero influence.

Here, the choice of shrinkage = 0.1 leads to 43 nonzero influential variables. We can now reduce the shrinkage factor drastically and observe the impact: 

> spam_gbm2 <- gbm(spam_Formula,distribution="bernoulli",
+       data=spam_Train2,n.trees=500,bag.fraction = 0.8,
+       shrinkage = 0.0001)
> spam_gbm2
gbm(formula = spam_Formula, distribution = "bernoulli", data = spam_Train2, 
    n.trees = 500, shrinkage = 1e-04, bag.fraction = 0.8)
A gradient boosted model with Bernoulli loss function.
500 iterations were performed.
There were 57 predictors of which 2 had nonzero influence.

收缩参数太低，因此几乎没有变量具有影响力。接下来，我们将生成以下图表：

> windows(height=100,width=200)
> par(mfrow=c(1,2))
> gbm.perf(spam_gbm,plot.it=TRUE)
Using OOB method...
[1] 151
Warning message:
In gbm.perf(spam_gbm, plot.it = TRUE) :
  OOB generally underestimates the optimal number of iterations although predictive performance is reasonably competitive. Using cv.folds>0 when calling gbm usually results in improved predictive performance.
> gbm.perf(spam_gbm2,plot.it=TRUE)
Using OOB method...
[1] 500
Warning message:
In gbm.perf(spam_gbm2, plot.it = TRUE) :
  OOB generally underestimates the optimal number of iterations although predictive performance is reasonably competitive. Using cv.folds>0 when calling gbm usually results in improved predictive performance.

gbm 包

图 3：收缩因子作为收敛的一个因素

我们对于极小的收缩因子没有清晰的收敛。

gbm 函数的更多详细信息可以从包文档或之前提供的源中获取。提升是一种非常通用的技术，Ridgeway 为各种数据结构实现了这一技术。下表列出了四种最常见的几种数据结构，展示了统计模型、偏差（与损失函数相关）、初始值、梯度和每个终端节点的输出估计：

输出类型	统计模型	梯度	终端节点估计
数值	高斯
二元	伯努利
计数	泊松
生存数据	Cox 比例风险模型	0	牛顿-拉夫森算法

表 1：GBM 提升选项

我们将应用 gbm 函数对计数数据和生存数据进行处理。

计数数据的提升

事故、错误/打字错误、出生等数量数据是流行的例子。在这里，我们计算特定时间、地点和/或空间内事件的数量。泊松分布非常适合用于建模数量数据。当我们有以协变量和独立变量形式提供的额外信息时，相关的回归问题通常很有兴趣。广义线性模型是建模数量数据的一种流行技术。

让我们看看可从 stats.idre.ucla.edu/r/dae/poisson-regression/ 获取的模拟数据集。必要的更改如本源所示进行。首先，我们使用 glm 函数拟合泊松回归模型。接下来，我们拟合提升模型，如下所示：

> # Poisson regression and boosting
> # https://stats.idre.ucla.edu/r/dae/poisson-regression/
> pregnancy <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
> pregnancy <- within(pregnancy, {
+   prog <- factor(prog, levels=1:3, 
+   labels=c("General", "Academic","Vocational"))
+   id <- factor(id)
+ })
> summary(pregnancy)
       id        num_awards           prog          math     
 1      :  1   Min.   :0.00   General   : 45   Min.   :33.0  
 2      :  1   1st Qu.:0.00   Academic  :105   1st Qu.:45.0  
 3      :  1   Median :0.00   Vocational: 50   Median :52.0  
 4      :  1   Mean   :0.63                    Mean   :52.6  
 5      :  1   3rd Qu.:1.00                    3rd Qu.:59.0  
 6      :  1   Max.   :6.00                    Max.   :75.0  
 (Other):194                                                 
> pregnancy_Poisson <- glm(num_awards ~ prog + math, 
+                     family="poisson", data=pregnancy)
> summary(pregnancy_Poisson)

Call:
glm(formula = num_awards ~ prog + math, family = "poisson", data = pregnancy)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.204  -0.844  -0.511   0.256   2.680  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -5.2471     0.6585   -7.97  1.6e-15 ***
progAcademic     1.0839     0.3583    3.03   0.0025 ** 
progVocational   0.3698     0.4411    0.84   0.4018    
math             0.0702     0.0106    6.62  3.6e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 287.67  on 199  degrees of freedom
Residual deviance: 189.45  on 196  degrees of freedom
AIC: 373.5

Number of Fisher Scoring iterations: 6

> pregnancy_boost <- gbm(num_awards ~ prog+math,dist="poisson",
+ n.trees=100,interaction.depth = 2,shrinkage=0.1,data=pregnancy)
> cbind(pregnancy$num_awards,predict(m1,type="response"),
+       predict(pboost,n.trees=100,type="response"))
    [,1]   [,2]   [,3]
1      0 0.1352 0.1240
2      0 0.0934 0.1072
3      0 0.1669 0.3375
4      0 0.1450 0.0850
5      0 0.1260 0.0257
6      0 0.1002 0.0735

195    1 1.0469 1.4832
196    2 2.2650 2.0241
197    2 1.4683 0.4047
198    1 2.2650 2.0241
199    0 2.4296 2.0241
200    3 2.6061 2.0241
> sum((pregnancy$num_awards-predict(m1,type="response"))²)
[1] 151
> sum((pregnancy$num_awards-predict(pboost,n.trees=100,
+    type="response"))²)
[1] 141
> summary(pregnancy_boost)
      var rel.inf
math math    89.7
prog prog    10.3

这是一个非常简单的例子，只有两个协变量。然而，在实践中，我们多次看到数量数据被当作回归问题处理。这是不幸的，一般回归技术并不是任何形式的炼金术。我们必须尊重数据结构，并且需要对数量数据进行分析。请注意，这里拟合的模型是非线性的。尽管这里的好处并不明显，但随着变量数量的增加，数量数据框架变得更加合适。我们将通过两个树的变量重要性图以及表格显示来结束对数量数据分析的讨论：

用于数量数据的提升

图 4：提升数量数据 – 变量重要性

> pretty.gbm.tree(pregnancy_boost,i.tree=18)
  SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight
0        1      64.50000        1         5           6          10.41    100
1        1      57.50000        2         3           4           1.14     90
2       -1      -0.01146       -1        -1          -1           0.00     71
3       -1       0.02450       -1        -1          -1           0.00     19
4       -1      -0.00387       -1        -1          -1           0.00     90
5       -1       0.05485       -1        -1          -1           0.00     10
6       -1       0.00200       -1        -1          -1           0.00    100
  Prediction
0    0.00200
1   -0.00387
2   -0.01146
3    0.02450
4   -0.00387
5    0.05485
6    0.00200
> pretty.gbm.tree(pregnancy_boost,i.tree=63)
  SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight
0        1      60.50000        1         5           6          3.837    100
1        0      20.00000        2         3           4          0.407     79
2       -1      -0.00803       -1        -1          -1          0.000     40
3       -1       0.05499       -1        -1          -1          0.000     39
4       -1       0.02308       -1        -1          -1          0.000     79
5       -1       0.02999       -1        -1          -1          0.000     21
6       -1       0.02453       -1        -1          -1          0.000    100
  Prediction
0    0.02453
1    0.02308
2   -0.00803
3    0.05499
4    0.02308
5    0.02999
6    0.02453

pretty.gbm.tree 函数有助于提取 gbm 对象的隐藏树。在下一节中，我们将处理生存数据的梯度提升技术。

生存数据的提升

pbc 数据集已在第一章和 集成技术简介 以及第二章和 自助法 中介绍。如前所述，生存数据有缺失观测值，我们需要专门的技术来处理这种情况。在表 1 中，我们看到偏差函数相当复杂。多亏了 Ridgeway，我们不必过多担心这类计算。相反，我们只需使用带有 dist="coxph" 选项的 gbm 函数，并按以下方式进行数据分析：

> # Survival data
> pbc_boost <- gbm(Surv(time,status==2)~trt + age + sex+ascites +
+                    hepato + spiders + edema + bili + chol + 
+                    albumin + copper + alk.phos + ast + trig + 
+                    platelet + protime + stage,
+                  n.trees=100,interaction.depth = 2,
+                  shrinkage=0.01,dist="coxph",data=pbc)
> summary(pbc_boost)
              var rel.inf
bili         bili  54.220
age           age  10.318
protime   protime   9.780
stage       stage   7.364
albumin   albumin   6.648
copper     copper   5.899
ascites   ascites   2.361
edema       edema   2.111
ast           ast   0.674
platelet platelet   0.246
alk.phos alk.phos   0.203
trig         trig   0.178
trt           trt   0.000
sex           sex   0.000
hepato     hepato   0.000
spiders   spiders   0.000
chol         chol   0.000
> pretty.gbm.tree(pbc_boost,i.tree=2)  # output suppressed
> pretty.gbm.tree(pbc_boost,i.tree=72) # output suppressed

因此，使用多功能的 gbm 函数，我们可以轻松地对各种数据结构执行梯度提升技术。

xgboost 软件包

xgboost R 包是梯度提升方法的优化、分布式实现。这是一个已知效率高、灵活且可移植的工程优化——有关更多详细信息及常规更新，请参阅github.com/dmlc/xgboost。这提供了并行树提升，因此在数据科学社区中被发现非常有用。特别是在www.kaggle.org的比赛中，许多获胜者使用了xgboost技术。Kaggle 获胜者的部分列表可在github.com/dmlc/xgboost/tree/master/demo#machine-learning-challenge-winning-solutions找到。

极端梯度提升实现的优点如下所示：

并行计算: 此包通过 OpenMP 启用并行处理，然后使用计算机的所有核心
正则化: 这通过结合正则化思想来避免过拟合问题
交叉验证: 执行交叉验证不需要额外的编码
剪枝: 这将树增长到最大深度，然后反向剪枝
缺失值: 缺失值在内部被处理
保存和重新加载: 这不仅有助于保存现有模型，还可以从上次停止的地方继续迭代
跨平台: 这适用于 Python、Scala 等

我们将使用书中早些时候看到的垃圾邮件数据集来说明这些想法。xgboost包的函数要求所有变量都是数值型，输出也应标记为0和1。此外，协变量矩阵和输出需要分别提供给xgboost R 包。因此，我们首先加载垃圾邮件数据集，然后创建分区和公式，如下所示：

> ## The xgboost Package
> data("spam")
> spam2 <- spam
> spam2$type <- as.numeric(spam2$type)-1
> head(data.frame(spam2$type,spam$type))
  spam2.type spam.type
1          1      spam
2          1      spam
3          1      spam
4          1      spam
5          1      spam
6          1      spam
> # 1 denotes spam, and 0 - nonspam
> set.seed(12345)
> Train_Test <- sample(c("Train","Test"),nrow(spam2),replace = TRUE,
+                      prob = c(0.7,0.3))
> spam2_Train <- spam2[Train_Test=="Train",]
> spam_Formula <- as.formula("type~.")
> spam_train <- list()

xgboost包还要求指定训练回归数据为特殊的dgCMatrix矩阵。因此，我们可以使用as函数将其转换：

> spam_train$data <- as(spam_train$data,"dgCMatrix")
> spam_train$label <- spam2_Train$type
> class(spam_train$data)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

我们现在已准备好应用xgboost函数的基础设施。选择nrounds=100和逻辑函数，结果如下所示：

> # Simple XGBoosting
> spam_xgb<- xgboost(data=spam_train$data,label=spam_train$label,+                     nrounds = 100,objective="binary:logistic")
[1]	train-error:0.064062 
[2]	train-error:0.063437 
[3]	train-error:0.053438 
[4]	train-error:0.050313 
[5]	train-error:0.047812 
[6]	train-error:0.045313 

[95]	train-error:0.002188 
[96]	train-error:0.001875 
[97]	train-error:0.001875 
[98]	train-error:0.001875 
[99]	train-error:0.000937 
[100]	train-error:0.000937

使用拟合的增强模型，我们现在应用predict函数并评估准确性：

> xgb_predict <- predict(spam_xgb,spam_train$data)
> sum(xgb_predict>0.5)
[1] 1226
> sum(spam_train$label)
[1] 1229
> table(spam_train$label,c(xgb_predict>0.5))

    FALSE TRUE
  0  1971    0
  1     3 1226

如果观察到的概率（在这种情况下是响应）被标记为 1 的概率超过 0.5，那么我们可以将观察到的标签标记为 1，否则为 0。通过使用 table R 函数，我们可以获得预测标签和实际标签的列联表。显然，我们具有非常好的准确性，并且只有三个错误分类。

我们声称xgboost包不需要额外的编码来进行交叉验证分析。xgb.cv函数在这里很有用，它使用与xgboost函数相同的参数，并通过nfold选项指定的交叉验证折数来工作。在这里，我们选择nfold=10。现在，使用xgb.cv函数，我们进行该分析并评估预测准确率：

> # XGBoosting with cross-validation
> spam_xgb_cv <- xgb.cv(data=spam_train$data,
+           label=spam_train$label,nfold=10,nrounds = 100,
+           objective="binary:logistic",prediction = TRUE)
[1]	train-error:0.064410+0.001426	test-error:0.091246+0.01697
[2]	train-error:0.058715+0.001862	test-error:0.082805+0.01421
[3]	train-error:0.052986+0.003389	test-error:0.077186+0.01472
[4]	train-error:0.049826+0.002210	test-error:0.073123+0.01544
[5]	train-error:0.046910+0.001412	test-error:0.070937+0.01340
[6]	train-error:0.043958+0.001841	test-error:0.066249+0.01346

[95]	train-error:0.001667+0.000340	test-error:0.048119+0.00926
[96]	train-error:0.001528+0.000318	test-error:0.047181+0.01008
[97]	train-error:0.001458+0.000260	test-error:0.046868+0.00974
[98]	train-error:0.001389+0.000269	test-error:0.047181+0.00979
[99]	train-error:0.001215+0.000233	test-error:0.047182+0.00969
[100]	train-error:0.001111+0.000260	test-error:0.045932+0.01115
> xgb_cv_predict <- spam_xgb_cv$pred
> sum(xgb_cv_predict>0.5)
[1] 1206
> table(spam_train$label,c(xgb_cv_predict>0.5))

    FALSE TRUE
  0  1909   62
  1    85 1144

交叉验证分析显示准确率有所下降。这表明我们在xgboost函数中遇到了过拟合问题。现在，我们将查看xgboost包的其他特性。在本节的开头，我们声称该技术通过提前停止和恢复早期拟合的模型对象提供了灵活性。

然而，一个重要的问题是*何时需要提前停止迭代？*我们没有关于迭代次数的任何基础理论，这个次数是作为变量数量和观测数量的函数。因此，我们将以指定的迭代次数开始这个过程。如果误差减少的收敛性没有低于阈值水平，那么我们将继续进行更多的迭代，这项任务将在下一部分进行。然而，如果指定的迭代次数过高，并且提升方法的性能正在变差，那么我们不得不停止迭代。这是通过指定early_stopping_rounds选项来实现的，我们将在下面的代码中将其付诸实践：

> # Stop early
> spam_xgb_cv2 <- xgb.cv(data=spam_train$data,label=
+             spam_train$label, early_stopping_rounds = 5,
+             nfold=10,nrounds = 100,objective="binary:logistic",
+             prediction = TRUE)
[1]	train-error:0.064271+0.002371	test-error:0.090294+0.02304
Multiple eval metrics are present. Will use test_error for early stopping.
Will train until test_error hasn't improved in 5 rounds.

[2]	train-error:0.059028+0.003370	test-error:0.085614+0.01808
[3]	train-error:0.052048+0.002049	test-error:0.075930+0.01388
[4]	train-error:0.049236+0.002544	test-error:0.072811+0.01333
[5]	train-error:0.046007+0.002775	test-error:0.070622+0.01419
[6]	train-error:0.042882+0.003065	test-error:0.066559+0.01670

[38]	train-error:0.010382+0.001237	test-error:0.048121+0.01153[39]	train-error:0.010069+0.001432	test-error:0.048434+0.01162
[40]	train-error:0.009653+0.001387	test-error:0.048435+0.01154
[41]	train-error:0.009236+0.001283	test-error:0.048435+0.01179
[42]	train-error:0.008924+0.001173	test-error:0.048748+0.01154
Stopping. Best iteration:
[37]	train-error:0.010625+0.001391	test-error:0.048121+0.01162

在这里，最佳迭代已经发生在编号37，这是通过early_stopping_rounds = 5选项获得的确认。现在我们已经找到了最佳迭代，我们停止这个过程。

现在，我们将查看如何添加更多的迭代。以下代码仅用于说明。使用nrounds = 10选项，以及早期拟合的spam_xgb，以及数据和标签的选项，我们将要求xgboost函数进行十次额外的迭代：

> # Continue training
> xgboost(xgb_model=spam_xgb,
+         data=spam_train$data,label=spam_train$label,
+         nrounds = 10)
[101]	train-error:0.000937 
[102]	train-error:0.000937 
[103]	train-error:0.000937 
[104]	train-error:0.000937 
[105]	train-error:0.000937 
[106]	train-error:0.000937 
[107]	train-error:0.000937 
[108]	train-error:0.000937 
[109]	train-error:0.000937 
[110]	train-error:0.000937 
##### xgb.Booster
raw: 136 Kb 
call:
  xgb.train(params = params, data = dtrain, nrounds = nrounds, 
    watchlist = watchlist, verbose = verbose, print_every_n = print_every_n, early_stopping_rounds = early_stopping_rounds, maximize = maximize, save_period = save_period, save_name = save_name, xgb_model = xgb_model, callbacks = callbacks)
params (as set within xgb.train):
  silent = "1"
xgb.attributes:
  niter
callbacks:
  cb.print.evaluation(period = print_every_n)
  cb.evaluation.log()
  cb.save.model(save_period = save_period, save_name = save_name)
niter: 110
evaluation_log:
    iter train_error
       1    0.064062
       2    0.063437
---                 
     109    0.000937
     110    0.000937

迭代次数的粗体和粗大字体并不是 R 软件输出的格式。这种变化是为了强调从早期拟合的spam_xgb对象开始的迭代次数现在将从101继续，并增加到110。使用xgboost函数很容易增加额外的迭代次数。

xgb.plot.importance函数与xgb.importance函数一起使用，可以用来提取并显示由提升方法确定的最重要的变量：

> # Variable Importance
> xgb.plot.importance(xgb.importance(names(spam_train$data),
+                                    spam_xgb)[1:10,])

结果是以下图表：

xgboost 包

图 5：xgboost 函数的变量重要性图

我们现在已经看到了xgboost包的力量。接下来，我们将概述h2o包的功能。

h2o 包

R 软件的.exe文件大小为 75 MB（版本 3.4.1）。h2o R 包的大小为 125 MB。这可能会让你意识到h2o包的重要性。本书中使用的所有数据集大小都非常有限，观测数不超过 10,000。在大多数情况下，文件大小最大仅为几 MB。然而，数据科学界的工作非常努力，经常处理 GB 级别甚至更高格式的文件。因此，我们需要更多的能力，而h2o包正好提供了这些能力。我们只需简单地加载h2o包，就可以一窥究竟：

> library(h2o)

----------------------------------------------------------------------

Your next step is to start H2O:
    > h2o.init()

For H2O package documentation, ask for help:
    > ??h2o

After starting H2O, you can use the Web UI at http://localhost:54321
For more information visit http://docs.h2o.ai

----------------------------------------------------------------------

Attaching package: 'h2o'

The following objects are masked from 'package:stats':

    cor, sd, var

The following objects are masked from 'package:base':

    %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
    colnames<-, ifelse, is.character, is.factor, is.numeric, log,
    log10, log1p, log2, round, signif, trunc

Warning message:
package 'h2o' was built under R version 3.4.4 
> h2o.init()

H2O is not running yet, starting it now...

Note:  In case of errors look at the following log files:
    C:\Users\tprabhan\AppData\Local\Temp\Rtmpu6f0fO/h2o_TPRABHAN_started_from_r.out
    C:\Users\tprabhan\AppData\Local\Temp\Rtmpu6f0fO/h2o_TPRABHAN_started_from_r.err

java version "1.7.0_67"
Java(TM) SE Runtime Environment (build 1.7.0_67-b01)
Java HotSpot(TM) 64-Bit Server VM (build 24.65-b04, mixed mode)

Starting H2O JVM and connecting: .. Connection successful!

R is connected to the H2O cluster: 
    H2O cluster uptime:         4 seconds 449 milliseconds 
    H2O cluster version:        3.16.0.2 
    H2O cluster version age:    7 months and 7 days !!! 
    H2O cluster name:           H2O_started_from_R_TPRABHAN_saz680 
    H2O cluster total nodes:    1 
    H2O cluster total memory:   1.76 GB 
    H2O cluster total cores:    4 
    H2O cluster allowed cores:  4 
    H2O cluster healthy:        TRUE 
    H2O Connection ip:          localhost 
    H2O Connection port:        54321 
    H2O Connection proxy:       NA 
    H2O Internal Security:      FALSE 
    H2O API Extensions:         AutoML, Algos, Core V3, Core V4 
    R Version:                  R version 3.4.0 (2017-04-21) 

Warning message:
In h2o.clusterInfo() : 
Your H2O cluster version is too old (7 months and 7 days)!
Please download and install the latest version from http://h2o.ai/download/

集群和线程有助于扩展计算。对于更热情的读者，以下资源将帮助使用h2o包：

使用gbm、xgboost和h2o包，读者可以分析复杂和大型数据集。

摘要

我们在本章开始时简要思考了为什么提升（boosting）有效。有三个视角可能解释了提升的成功，这些内容在我们深入探讨这个主题之前已经讨论过了。gbm包非常强大，它提供了调整梯度提升算法的不同选项，该算法可以处理多种数据结构。我们通过收缩选项展示了其能力，并将其应用于计数和生存数据结构。xgboost包是梯度提升方法的更高效实现。它更快，还提供了其他灵活性。我们通过使用xgboost函数进行交叉验证、早期停止以及根据需要继续迭代来展示其使用方法。h2o包/平台有助于在大规模上实现集成机器学习技术。

在下一章中，我们将探讨为什么集成（ensembling）有效。特别是，我们将看到为什么将多个模型组合在一起通常是一种有用的实践，我们还将探讨我们可以这样做的情况。

第七章. 通用集成技术

前四章讨论了决策树的集成技术。在这些章节中讨论的每个主题中，基础学习器都是决策树，因此我们深入研究了同质集成技术。在本章中，我们将展示基础学习器可以是任何统计或机器学习技术，并且它们的集成将提高预测的精度。一个重要要求是基础学习器应该比随机猜测更好。通过 R 程序，我们将讨论和阐明集成将有效的情况。投票是分类器的一个重要特性——我们将陈述两种不同的方法，并在 bagging 和随机森林集成器的情况下进行说明。平均技术是回归变量的集成器，它将遵循分类方法的讨论。本章将以对第一章中非正式介绍的堆叠方法的详细讨论结束，集成技术简介。主题流程如下：

集成为什么有效？
通过投票进行集成
通过平均进行集成
栈集成

技术要求

本章将使用以下库：

rpart
randomForest

集成为什么有效？

当使用 bagging 方法时，我们将许多决策树的结果结合起来，通过多数计数来生成单个输出/预测。在不同的采样机制下，随机森林的结果被结合起来生成单个预测。在决策树的顺序误差减少方法下，提升方法也提供了改进的答案。尽管我们处理的是不确定数据，这涉及到概率，但我们不打算有那种给出黑盒结果且行为不一致的方法论。一个理论应该解释其工作原理，我们需要保证结果的一致性，并且对此没有神秘之处。任意和不确定的答案是完全不受欢迎的。在本节中，我们将探讨集成解决方案如何工作以及它们不工作的场景。

集成方法有强大的数学和统计基础，解释了为什么它们能给出这样的解决方案。我们首先考虑分类问题。我们将从一个简化的设置开始，并假设我们有T个相互独立的分类器，并且每个分类器相关的准确率与为什么集成工作？相同。这是一个最简单的情况，我们稍后会推广这个场景。现在，如果我们有T个分类器，并且每个分类器对+1 或-1 等观察进行投票，这就引出了一个问题，整体准确率会是什么？由于T个分类器的正确分类数量必须超过错误分类数量，我们需要至少为什么集成工作？个分类器来投票正确的结果。在这里，表示小于给定分数的最大整数。只要或更多数量的分类器为正确的类别投票，多数分类就是正确的。

为了澄清，重要的是要注意，当我们说一个分类器的准确率是p时，我们并不是指分类器将观察标记为+1 的概率是p。相反，我们在这里的意思是，如果分类器做出 100 次预测，预测可以是+1 和-1 的任何组合；100p次预测被分类器正确识别。准确率与+1 和-1 在人群中的分布无关。

在这种设置下，分类器标记正确观察的概率遵循参数为n = T和概率p的二项分布。因此，多数投票得到正确预测的概率如下：

为什么集成工作？

既然我们已经提到分类器必须比随机猜测更好，我们就需要分类器的准确率超过 0.5。然后我们将逐步增加准确率，并观察分类器数量的增加如何影响多数投票的概率：

> source("Utilities.R")
> windows(height=100,width=100)
> # Ensembling in Classification
> # Illustrating the ensemble accuracy with same accuracy for each classifier
> # Different p's and T's with p > 0.5
> classifiers <- seq(9,45,2) # Number of classifiers 
> accuracy <- seq(0.55,0.85,.05)
> plot(0,type='n',xlim=range(classifiers),ylim=c(0.6,1),
+      xlab="Number of Classifiers",ylab="Probability of Majority Voting")
> for(i in 1:length(accuracy)){
+   Prob_MV <- NULL
+   for(j in 1:length(classifiers)){
+     Prob_MV[j] <- sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+        prob=accuracy[i],size=classifiers[j]))
+   }
+   points(classifiers,Prob_MV,col=i,"l")
+ }
> title("Classifiers with Accuracy Better Than Random Guess")

seq函数设置了一个奇数序列，表示classifiers R数值向量中分类器的数量。准确率百分比在accuracy向量中从0.55到0.85不等。为了启动这个过程，我们设置了一个带有适当的x和y轴标签的空plot。现在，对于每个准确率值，我们将计算范围floor(classifiers[j]/2+1):classifiers[j]内多数投票的概率。floor(./2+1)确保我们选择了正确的起点。例如，如果分类器的数量是九个，那么floor(./2+1)的值是5。此外，当我们有九个分类器时，我们需要至少五个赞成事件发生的投票。另一方面，对于偶数个分类器（例如，八个）floor(./2+1)的值也是5。dbinom函数计算给定大小和概率的特定值的概率。在floor(classifiers[j]/2+1): classifiers[j]的范围内，它给出了多数投票的概率，或者多数投票的准确率。前面代码的输出在图 1中展示。我们可以从结果中看到，随着分类器数量的增加（每个分类器具有相同的准确率且优于随机猜测），多数投票的准确率也在增加：

为什么集成方法有效？

图 1：为什么集成方法应该有效？

这将帮助我们看到某个准确率选择下的Prob_MV值——例如，0.65。我们将分别对prob=0.65运行循环，观察随着分类器数量的增加，多数投票的准确率是如何提高的：

> Prob_MV <- NULL
> for(j in 1:length(classifiers)){
+   Prob_MV[j] <- sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+                            prob=0.65,size=classifiers[j]))
+ }
> Prob_MV
 [1] 0.8282807 0.8513163 0.8705318 0.8867689 0.9006211 0.9125264 0.9228185
 [8] 0.9317586 0.9395551 0.9463770 0.9523633 0.9576292 0.9622714 0.9663716
[15] 0.9699991 0.9732133 0.9760651 0.9785984 0.9808513

因此，随着具有相同准确率的分类器数量的增加，我们可以看到多数投票的准确率也在增加。值得注意的是，尽管我们每个分类器的准确率仅为0.65，但集成方法的准确率要高得多，几乎成为了一个完美的分类器。这就是集成方法的主要优势。

集成方法是否对任何类型的分类器都有帮助？如果我们有准确率低于随机猜测（即小于0.5）的分类器，那么我们将以与之前相同的方式进行搜索。对于许多准确率低于0.5的分类器，我们将计算多数投票分类器的准确率：

> # When p < 0.5, ensemble accuracy goes to zero
> classifiers <- seq(6,50,2)
> accuracy <- seq(0.45,0.05,-0.05)
> plot(0,type='n',xlim=range(classifiers),ylim=c(0,0.3),
+      xlab="Number of Classifiers",ylab="Probability of Majority Voting")
> for(i in 1:length(accuracy)){
+   Prob_MV <- NULL
+   for(j in 1:length(classifiers)){
+     Prob_MV[j] <- sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+                              prob=accuracy[i],size=classifiers[j]))
+   }
+   points(classifiers,Prob_MV,col=i,"l")
+   }
> title("Classifiers with Accuracy Worse Than Random Guess")

前面 R 程序的结果显示在 图 2 中。现在，第一个观察结果是，无论准确率是接近 0.5 还是接近 0，多数投票分类器的概率/准确率都在下降，这不利于性能。在所有情况下，我们都看到准确率最终会接近零。R 代码块中的变化是分类器序列 seq(6,50,2)，准确率水平从 0.45 下降到 0.05，在 seq(0.45,0.05,-0.05) 中。现在，考虑准确率略小于 0.5 的情况。例如，让我们将其保持在 0.4999。我们现在会幸运地看到性能改进吗？

为什么集成有效？

图 2：集成不是炼金术！

> classifiers <- seq(10,200,10)
> Prob_MV <- NULL
> for(j in 1:length(classifiers)){
+   Prob_MV[j] <- sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+                            prob=0.4999,size=classifiers[j]))
+ }
> Prob_MV
 [1] 0.3767071 0.4115491 0.4273344 0.4368132 0.4433011 0.4480955 0.4518222
 [8] 0.4548247 0.4573097 0.4594096 0.4612139 0.4627854 0.4641698 0.4654011
[15] 0.4665053 0.4675025 0.4684088 0.4692370 0.4699975 0.4706989

再次，我们发现我们无法匹配单个分类器的准确率。因此，我们有一个重要且关键的条件，即分类器必须比随机猜测更好。那么随机猜测本身呢？假装我们有一系列都是随机猜测的分类器并不困难。如果集成随着随机猜测的性能提高，我们通常不需要构建任何统计或机器学习技术。给定一组随机猜测，我们总能提高准确率。让我们来看看。

有两种情况——分类器数量为奇数和偶数——我们为这两种情况都提供了程序：

> accuracy <- 0.5
> classifiers <- seq(5,45,2)
> Prob_MV <- NULL
> for(j in 1:length(classifiers)){
+   Prob_MV[j] <- sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+                            prob=accuracy,size=classifiers[j]))
+   }
> Prob_MV
 [1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
[19] 0.5 0.5 0.5
> classifiers <- seq(10,50,2)
> Prob_MV <- NULL
> for(j in 1:length(classifiers)){
+   Prob_MV[j] <- (sum(dbinom(floor(classifiers[j]/2):classifiers[j],
+                             prob=accuracy,size=classifiers[j]))+
+                    sum(dbinom(floor(classifiers[j]/2+1):classifiers[j],
+                               prob=accuracy,size=classifiers[j])))/2
+   }
> Prob_MV
 [1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
[19] 0.5 0.5 0.5

这很有趣！无论分类器的数量如何，随机猜测的集成保持不变。在这里，既没有改进也没有恶化。因此，为了集成目的，我们总是需要比随机猜测更好的分类器。

在理解集成如何工作时，遵循你的直觉是很好的。我们从一个所有模型都具有相同准确率的过度简化假设开始，但如果我们处理具有不同准确率的模型，这样的假设就不适用了。因此，我们需要考虑可能对不同分类器有不同的准确率的情况。我们首先考虑每个分类器的准确率都高于 0.5，或者每个分类器都比随机猜测更好的情况。找到多数投票准确率的方法是评估分类器结果的每种可能组合的概率。我们考虑分类器数量为奇数时的简单情况。

假设我们拥有 T 个分类器，每个分类器的准确率如为什么集成有效？所示。请注意，因为这些对应于不同的度量。

评估具有不等准确率的多数投票概率所涉及步骤如下：

列出所有可能的基本事件。如果每个分类器对一个给定案例投票为真或假，这意味着它有两种可能的结果，以及 T 个分类器。列出可能的结果：
- 示例：如果我们有三个分类器，那么会有八种可能的情况，如下所示：
  
  分类器 1 分类器 2 分类器 3
  真真真
  假真真
  真假真
  假假真
  真真假
  假真假
  真假假
  假假假
计算每个可能事件的概率。由于每个分类器的准确度不同，因此每个可能结果的概率也会不同：
- 示例：如果三个分类器（对于真）的准确度分别为 0.6、0.7 和 0.8，那么假的概率分别为 0.4、0.3 和 0.2，前表中的概率如下：
  
  分类器 1 分类器 2 分类器 3
  0.6 0.7 0.8
  0.4 0.7 0.8
  0.6 0.3 0.8
  0.4 0.3 0.8 0.336
  0.6 0.7 0.2
  0.4 0.7 0.2
  0.6 0.3 0.2
  0.4 0.3 0.2
在下一步中，获取基本事件的概率，这将是每列数字的乘积：

分类器 1 分类器 2 分类器 3 概率
0.6 0.7 0.8 0.336
0.4 0.7 0.8 0.224
0.6 0.3 0.8 0.144
0.4 0.3 0.8 0.096
0.6 0.7 0.2 0.084
0.4 0.7 0.2 0.056
0.6 0.3 0.2 0.036
0.4 0.3 0.2 0.024
找出具有多数计数的事件。在这种情况下，这指的是大于或等于 2 的总和：

分类器 1 分类器 2 分类器 3 投票数
真真真 3
假真真 2
真假真 2
假假真 1
真真假 2
假真假 1
真假假 1
假假假 0
多数投票的概率然后就是投票数大于或等于 2 的情况的概率之和。这是概率列中第 1、2、3 和 5 行条目的总和，即 0.336 + 0.224 + 0.144 + 0.084 = 0.788。

分类器 1	分类器 2	分类器 3
真	真	真
假	真	真
真	假	真
假	假	真
真	真	假
假	真	假
真	假	假
假	假	假

分类器 1	分类器 2	分类器 3
0.6	0.7	0.8
0.4	0.7	0.8
0.6	0.3	0.8
0.4	0.3	0.8	0.336
0.6	0.7	0.2
0.4	0.7	0.2
0.6	0.3	0.2
0.4	0.3	0.2

分类器 1	分类器 2	分类器 3	概率
0.6	0.7	0.8	0.336
0.4	0.7	0.8	0.224
0.6	0.3	0.8	0.144
0.4	0.3	0.8	0.096
0.6	0.7	0.2	0.084
0.4	0.7	0.2	0.056
0.6	0.3	0.2	0.036
0.4	0.3	0.2	0.024

分类器 1	分类器 2	分类器 3	投票数
真	真	真	3
假	真	真	2
真	假	真	2
假	假	真	1
真	真	假	2
假	真	假	1
真	假	假	1
假	假	假	0

我们需要在这里定义一个名为 Get_Prob 的函数，如下所示：

> Get_Prob <- function(Logical,Probability){
+   return(t(ifelse(Logical,Probability,1-Probability)))
+ }

给定一个逻辑向量和相应的概率向量，Get_Prob 函数将返回一个向量，该向量包含逻辑条件为 真 的概率。如果逻辑值为 假，则返回补数（1 – 概率）。

上述步骤被放入 R 程序中，如下所示：

> # Different accuracies T's illustration
> # For simplicity, we set the number of classifiers at odd number
> # Each p_i's greater than 0.5
> accuracy <- c(0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9)
> NT <- length(accuracy) # Number of classifiers 
> APC <- expand.grid(rep(list(c(TRUE,FALSE)),NT)) # All possible combinations
> head(APC)
   Var1  Var2  Var3 Var4 Var5 Var6 Var7 Var8 Var9
1  TRUE  TRUE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
2 FALSE  TRUE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
3  TRUE FALSE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
4 FALSE FALSE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
5  TRUE  TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
6 FALSE  TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
> Elements_Prob <- t(apply(APC,1,Get_Prob,Probability=accuracy))
> head(Elements_Prob)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,]  0.5 0.55  0.6 0.65  0.7 0.75  0.8 0.85  0.9
[2,]  0.5 0.55  0.6 0.65  0.7 0.75  0.8 0.85  0.9
[3,]  0.5 0.45  0.6 0.65  0.7 0.75  0.8 0.85  0.9
[4,]  0.5 0.45  0.6 0.65  0.7 0.75  0.8 0.85  0.9
[5,]  0.5 0.55  0.4 0.65  0.7 0.75  0.8 0.85  0.9
[6,]  0.5 0.55  0.4 0.65  0.7 0.75  0.8 0.85  0.9
> Events_Prob <- apply(Elements_Prob,1,prod)
> Majority_Events <- (rowSums(APC)>NT/2)
> sum(Events_Prob*Majority_Events)
[1] 0.9112646

给定一个名为 accuracy 的包含准确率的数值向量，其中分类器的数量为奇数，我们首先使用 length 函数找到其中的分类器数量，并将其存储在 NT 中。然后使用 expand.grid 函数生成所有可能的 APC 组合，其中 rep 函数将向量 (TRUE, FALSE) NT 重复 NT 次。APC 对象的每一列将生成一个列，其中 TRUE 和 FALSE 条件将使用 Get_Prob 函数替换为相应的分类器准确率以及适当的补数。由于我们考虑的是奇数个分类器，当该基本事件中的 TRUE 数量大于分类器数量的 50%（即大于 NT/2）时，才会进行多数投票。其余的计算比较容易理解。如果九个分类器的准确率分别为 0.5、0.55、0.6、0.65、0.7、0.75、0.8、0.85 和 0.9，那么计算表明集成准确率为 0.9113，高于这里最准确的分类器，即 0.9。然而，我们必须记住，八个分类器中的每一个的准确率都低于 0.9。尽管如此，集成准确率仍然高于我们手头上的最高分类器。为了验证计算是否正常工作，我们将这种方法应用于周（2012）第 74 页上的示例，并确认最终多数投票概率为 0.933：

> accuracy <- c(0.7,0.7,0.7,0.9,0.9)
> NT <- length(accuracy) # Number of classifiers
> APC <- expand.grid(rep(list(c(TRUE,FALSE)),NT)) # All possible combinations
> Elements_Prob <- t(apply(APC,1,Get_Prob,Probability=accuracy))
> Events_Prob <- apply(Elements_Prob,1,prod)
> Majority_Events <- (rowSums(APC)>NT/2)
> sum(Events_Prob*Majority_Events)
[1] 0.93268

当每个分类器都不如随机猜测时会发生什么？我们将简单地输出九个分类器场景的准确率，并重复程序以获得以下答案：

> # Each p_i's lesser than 0.5
> accuracy <- 1-c(0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9)
> NT <- length(accuracy) # Number of classifiers
> APC <- expand.grid(rep(list(c(TRUE,FALSE)),NT)) # All possible combinations
> head(APC)
   Var1  Var2  Var3 Var4 Var5 Var6 Var7 Var8 Var9
1  TRUE  TRUE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
2 FALSE  TRUE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
3  TRUE FALSE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
4 FALSE FALSE  TRUE TRUE TRUE TRUE TRUE TRUE TRUE
5  TRUE  TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
6 FALSE  TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
> Elements_Prob <- t(apply(APC,1,Get_Prob,Probability=accuracy))
> head(Elements_Prob)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,]  0.5 0.45  0.4 0.35  0.3 0.25  0.2 0.15  0.1
[2,]  0.5 0.45  0.4 0.35  0.3 0.25  0.2 0.15  0.1
[3,]  0.5 0.55  0.4 0.35  0.3 0.25  0.2 0.15  0.1
[4,]  0.5 0.55  0.4 0.35  0.3 0.25  0.2 0.15  0.1
[5,]  0.5 0.45  0.6 0.35  0.3 0.25  0.2 0.15  0.1
[6,]  0.5 0.45  0.6 0.35  0.3 0.25  0.2 0.15  0.1
> Events_Prob <- apply(Elements_Prob,1,prod)
> Majority_Events <- (rowSums(APC)>NT/2)
> sum(Events_Prob*Majority_Events)
[1] 0.08873544

当每个分类器都不如随机猜测时，在集成的情况下，多数投票分类器会给出可怕的结果。这让我们面临最后一个情况。如果我们有一组分类器，其中一些比随机猜测分类器好，而一些比随机猜测分类器差呢？我们将计算代码块放入一个名为 Random_Accuracy 的函数中。然后，分类器中的准确率变成了单位区间内随机生成的数字。Random_Accuracy 函数随后运行十次，生成以下输出：

> # Mixture of p_i's, some > 0.5, and some < 0.5
> Random_Accuracy <- function() {
+   accuracy <- runif(9)
+   NT <- length(accuracy) 
+   APC <- expand.grid(rep(list(c(TRUE,FALSE)),NT)) 
+   Elements_Prob <- t(apply(APC,1,Get_Prob,Probability=accuracy))
+   Events_Prob <- apply(Elements_Prob,1,prod)
+   Majority_Events <- (rowSums(APC)>NT/2)
+   return(sum(Events_Prob*Majority_Events))
+ }
> Random_Accuracy()
[1] 0.3423631
> Random_Accuracy()
[1] 0.3927145
> Random_Accuracy()
[1] 0.5341844
> Random_Accuracy()
[1] 0.1624876
> Random_Accuracy()
[1] 0.4065803
> Random_Accuracy()
[1] 0.4687087
> Random_Accuracy()
[1] 0.7819835
> Random_Accuracy()
[1] 0.3124515
> Random_Accuracy()
[1] 0.6842173
> Random_Accuracy()
[1] 0.2531727

结果参差不齐。因此，如果我们需要从集成方法中获得合理的准确性和性能，确保每个分类器都比随机猜测要好是至关重要的。到目前为止，我们分析的一个核心假设是分类器之间是相互独立的。在实际设置中，这个假设很少成立，因为分类器是使用相同的训练集构建的。然而，这个话题将在下一章中讨论。

我们现在将转向投票集成的问题。

投票集成

通过投票集成可以有效地用于分类问题。我们现在有一组分类器，我们需要使用它们来预测未知案例的类别。分类器预测的组合可以以多种方式进行。我们将考虑的两个选项是多数投票和加权投票。

多数投票

通过使用决策树作为基础学习器构建的集成，我们可以通过投票相关的想法进行说明，正如在开发袋装和随机森林时使用的那样。首先，我们将使用 randomForest 函数创建 500 个基础学习器，并重复第一个块中的程序，如第四章中所示，随机森林。集成已经在那一章中完成，我们将在这里详细说明那些步骤。首先，给出设置随机森林的代码块：

> load("../Data/GC2.RData")
> set.seed(12345)
> Train_Test <- sample(c("Train","Test"),nrow(GC2),
+ replace = TRUE,prob = c(0.7,0.3))
> GC2_Train <- GC2[Train_Test=="Train",]
> GC2_TestX <- within(GC2[Train_Test=="Test",],rm(good_bad))
> GC2_TestY <- GC2[Train_Test=="Test","good_bad"]
> GC2_Formula <- as.formula("good_bad~.")
> # RANDOM FOREST ANALYSIS
> GC2_RF <- randomForest(GC2_Formula,data=GC2_Train,keep.inbag=TRUE,
+                        ntree=500)

接下来，我们将使用标准的 predict 函数来预测 GC2_TestX 数据的类别，然后，使用 predict.all=TRUE 选项，获取随机森林中生成的每个树的预测结果：

> # New data voting
> GC2_RF_Test_Margin <- predict(GC2_RF,newdata = GC2_TestX,
+                          type="class")
> GC2_RF_Test_Predict <- predict(GC2_RF,newdata=GC2_TestX,
+                           type="class",predict.all=TRUE
+                           )

预测的 GC2_RF_Test_Predict 对象将包含进一步的 individual 对象，这些对象将包含每个决策树的预测。我们首先定义一个名为 Row_Count_Max 的函数，该函数将返回森林中计数最大的预测。然后，基本的投票方法将在以下代码块中与 predict 函数的结果进行比较：

> Row_Count_Max <- function(x) names(which.max(table(x))) 
> # Majority Voting
> Voting_Predict <- apply(GC2_RF_Test_Predict$individual,1,
+ Row_Count_Max)
> head(Voting_Predict);tail(Voting_Predict)
     1      2      3      4      9     10 
"good"  "bad" "good"  "bad" "good"  "bad" 
   974    980    983    984    988    996 
 "bad"  "bad" "good" "good" "good" "good" 
> all(Voting_Predict==GC2_RF_Test_Predict$aggregate)
[1] TRUE
> all(Voting_Predict==GC2_RF_Test_Margin)
[1] TRUE
> sum(Voting_Predict==GC2_TestY)/313
[1] 0.7795527

因此，我们可以看到 predict 函数实现了多数计数技术。接下来，我们将快速说明加权投票背后的思想和思考。

加权投票

在简单投票的使用中存在的一个隐含假设是所有分类器都是同样准确的，或者所有分类器都有相同的投票权。考虑一个更简单的情况，我们有五个分类器，其中三个的准确率为 0.51，剩下的两个准确率为 0.99。如果准确度较低的分类器将观察结果投票为负案例(-1)，而两个更准确的分类器将其投票为正案例(+1)，那么简单的投票方法将把观察结果标记为(-1)。在这种投票模式中，观察结果为-1 的概率是，而为+1 的概率是。因此，我们不能假装所有分类器都应该有相同的投票权。这就是我们将充分利用加权投票方法的地方。

在这次分析中，我们将使用训练数据集上分类器的准确率作为权重。我们将加权投票视为与相关的权重。权重的一个重要特征是它们应该是非负的，并且它们的总和应该是 1，即。我们将归一化分类器的准确率以满足这一约束。

我们将继续使用德国信用数据集进行分析。首先，我们将获得训练数据集上 500 棵树的预测，然后获得准确率：

> # Analyzing Accuracy of Trees of the Fitted Forest
> GC2_RF_Train_Predict <- predict(GC2_RF,newdata=GC2_Train[,-20],
+                                 type="class",predict.all=TRUE)
> head(GC2_RF_Train_Predict$individual[,c(1:5,496:500)])  
   [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]   [,10] 
5  "bad"  "bad"  "bad"  "bad"  "good" "bad"  "bad"  "bad"  "bad"  "bad" 
6  "good" "good" "good" "good" "good" "good" "bad"  "bad"  "bad"  "good"
7  "good" "good" "good" "good" "good" "good" "good" "good" "good" "good"
8  "good" "good" "good" "good" "good" "bad"  "good" "bad"  "good" "good"
11 "bad"  "bad"  "bad"  "bad"  "bad"  "bad"  "bad"  "bad"  "bad"  "bad" 
12 "good" "bad"  "bad"  "bad"  "bad"  "good" "bad"  "bad"  "bad"  "bad" 
> RF_Tree_Train_Accuracy <- NULL
> for(i in 1:GC2_RF$ntree){
+   RF_Tree_Train_Accuracy[i] <- sum(GC2_RF_Train_Predict$individual[,i]==
+                                   GC2_Train$good_bad)/nrow(GC2_Train)
+ }
> headtail(sort(RF_Tree_Train_Accuracy),10)
 [1] 0.8340611 0.8369723 0.8384279 0.8398836 0.8398836 0.8413392 0.8413392
 [8] 0.8413392 0.8413392 0.8427948 0.8908297 0.8908297 0.8908297 0.8908297
[15] 0.8922853 0.8922853 0.8937409 0.8937409 0.8966521 0.8981077

headtail函数是什么？它在Utilities.R文件中可用。以下是对bagging集成器进行重复分析：

> # Bagging ANALYSIS
> GC2_Bagg <- randomForest(GC2_Formula,data=GC2_Train,keep.inbag=TRUE,
+                          mtry=ncol(GC2_TestX),ntree=500)
> GC2_Bagg_Test_Predict <- predict(GC2_Bagg,newdata=GC2_TestX,
+                                 type="class",predict.all=TRUE)
> GC2_Bagg_Train_Predict <- predict(GC2_Bagg,newdata=GC2_Train[,-20],
+                                 type="class",predict.all=TRUE)
> Bagg_Tree_Train_Accuracy <- NULL
> for(i in 1:GC2_Bagg$ntree){
+   Bagg_Tree_Train_Accuracy[i] <- sum(GC2_Bagg_Train_Predict$individual[,i]==
+                                   GC2_Train$good_bad)/nrow(GC2_Train)
+ }
> headtail(sort(Bagg_Tree_Train_Accuracy),10)
 [1] 0.8369723 0.8384279 0.8413392 0.8457060 0.8457060 0.8471616 0.8471616
 [8] 0.8471616 0.8471616 0.8486172 0.8966521 0.8966521 0.8966521 0.8966521
[15] 0.8966521 0.8981077 0.8995633 0.8995633 0.9024745 0.9097525

接下来，我们将归一化权重并计算测试样本中观测值的加权投票，如下所示：

> # Weighted Voting with Random Forest
> RF_Weights <- RF_Tree_Train_Accuracy/sum(RF_Tree_Train_Accuracy)
> Bagg_Weights <- Bagg_Tree_Train_Accuracy/sum(Bagg_Tree_Train_Accuracy)
> RF_Weighted_Vote <- data.frame(matrix(0,nrow(GC2_TestX),ncol=3))
> names(RF_Weighted_Vote) <- c("Good_Weight","Bad_Weight","Prediction")
> for(i in 1:nrow(RF_Weighted_Vote)){
+   RF_Weighted_Vote$Good_Weight[i] <- 
+     sum((GC2_RF_Test_Predict$individual[i,]=="good")*RF_Weights)
+   RF_Weighted_Vote$Bad_Weight[i] <- 
+     sum((GC2_RF_Test_Predict$individual[i,]=="bad")*RF_Weights)
+   RF_Weighted_Vote$Prediction[i] <- c("good","bad")[which.max(RF_Weighted_Vote[i,1:2])]
+ }
> head(RF_Weighted_Vote,10)
   Good_Weight Bad_Weight Prediction
1    0.8301541 0.16984588       good
2    0.3260033 0.67399668        bad
3    0.8397035 0.16029651       good
4    0.4422527 0.55774733        bad
5    0.9420565 0.05794355       good
6    0.2378956 0.76210442        bad
7    0.4759756 0.52402435        bad
8    0.7443038 0.25569624       good
9    0.8120180 0.18798195       good
10   0.7799587 0.22004126       good

如下所示，对bagging对象重复进行加权投票分析：

> # Weighted Voting with Bagging
> Bagg_Weights <- Bagg_Tree_Train_Accuracy/sum(Bagg_Tree_Train_Accuracy)
> Bagg_Weights <- Bagg_Tree_Train_Accuracy/sum(Bagg_Tree_Train_Accuracy)
> Bagg_Weighted_Vote <- data.frame(matrix(0,nrow(GC2_TestX),ncol=3))
> names(Bagg_Weighted_Vote) <- c("Good_Weight","Bad_Weight","Prediction")
> for(i in 1:nrow(Bagg_Weighted_Vote)){
+   Bagg_Weighted_Vote$Good_Weight[i] <- 
+     sum((GC2_Bagg_Test_Predict$individual[i,]=="good")*Bagg_Weights)
+   Bagg_Weighted_Vote$Bad_Weight[i] <- 
+     sum((GC2_Bagg_Test_Predict$individual[i,]=="bad")*Bagg_Weights)
+   Bagg_Weighted_Vote$Prediction[i] <- c("good","bad")[which.max(Bagg_Weighted_Vote[i,1:2])]
+ }
> head(Bagg_Weighted_Vote,10)
   Good_Weight Bad_Weight Prediction
1    0.9279982 0.07200181       good
2    0.1634505 0.83654949        bad
3    0.8219618 0.17803818       good
4    0.4724477 0.52755226        bad
5    0.9619528 0.03804725       good
6    0.1698628 0.83013718        bad
7    0.4540574 0.54594265        bad
8    0.7883772 0.21162281       good
9    0.8301772 0.16982283       good
10   0.7585720 0.24142804       good

现在，随着投票机制的问题解决，我们将注意力转向回归问题。

通过平均进行集成

在回归模型的背景下，预测是感兴趣变量的数值。结合由于各种集成器导致的输出预测相对简单；由于集成机制，我们只需将集成器之间预测值的平均值解释为预测值。在分类问题的背景下，我们可以进行简单的平均和加权平均。在前一节中，集成器具有同质的基础学习器。然而，在本节中，我们将处理异质的基础学习器。

现在，我们将考虑一个在第八章《集成诊断》中详细处理的回归问题。问题是基于超过 60 个解释变量的房价预测。我们拥有训练和测试数据集，并将它们加载以启动过程：

> # Averaging for Regression Problems
> load("../Data/ht_imp_author.Rdata") # returns ht_imp object
> load("../Data/htest_imp_author.Rdata") # returns htest_imp
> names(ht_imp)[69] <- "SalePrice"
> dim(ht_imp)
[1] 1460   69
> dim(htest_imp)
[1] 1459   68

因此，我们有大量观测数据来构建我们的模型。SalePrice是这里感兴趣的变量。首先，我们创建一个公式并构建一个线性模型；四个不同深度的回归树；四个具有不同隐藏神经元的神经网络；以及以下代码块中的支持向量机模型：

> hf <- as.formula("SalePrice~.")
> SP_lm <- lm(hf,data=ht_imp)
> SP_rpart2 <- rpart(hf,data=ht_imp,maxdepth=2)
> SP_rpart4 <- rpart(hf,data=ht_imp,maxdepth=4)
> SP_rpart6 <- rpart(hf,data=ht_imp,maxdepth=6)
> SP_rpart8 <- rpart(hf,data=ht_imp,maxdepth=8)
> SP_nn2 <- nnet(hf,data=ht_imp,size=2,linout=TRUE)
# weights:  267
initial  value 56996872361441.906250 
final  value 9207911334609.976562 
converged
> SP_nn3 <- nnet(hf,data=ht_imp,size=3,linout=TRUE)
# weights:  400
initial  value 56997125121706.257812 
final  value 9207911334609.960938 
converged
> SP_nn4 <- nnet(hf,data=ht_imp,size=4,linout=TRUE)
# weights:  533
initial  value 56996951452602.304687 
iter  10 value 19328028546738.226562
iter  20 value 19324281941793.617187
final  value 9080312934601.205078 
converged
> SP_nn5 <- nnet(hf,data=ht_imp,size=5,linout=TRUE)
# weights:  666
initial  value 56997435951836.507812 
final  value 9196060713131.609375 
converged
> SP_svm <- svm(hf,data=ht_imp)

我们已经有了考虑异质集成的所需设置。

简单平均

我们使用训练数据集构建了十个模型，现在我们将使用predict函数在这些模型上对训练数据集进行预测，如下所示：

> # Simple Averaging
> SP_lm_pred <- predict(SP_lm,newdata=htest_imp)
Warning message:
In predict.lm(SP_lm, newdata = htest_imp) :
  prediction from a rank-deficient fit may be misleading
> SP_rpart2_pred <- predict(SP_rpart2,newdata=htest_imp)
> SP_rpart4_pred <- predict(SP_rpart4,newdata=htest_imp)
> SP_rpart6_pred <- predict(SP_rpart6,newdata=htest_imp)
> SP_rpart8_pred <- predict(SP_rpart8,newdata=htest_imp)
> SP_nn2_pred <- predict(SP_nn2,newdata=htest_imp)
> SP_nn3_pred <- predict(SP_nn3,newdata=htest_imp)
> SP_nn4_pred <- predict(SP_nn4,newdata=htest_imp)
> SP_nn5_pred <- predict(SP_nn5,newdata=htest_imp)
> SP_svm_pred <- predict(SP_svm,newdata=htest_imp)

当涉及到分类问题时，预测要么基于类别标签，要么基于感兴趣类别的概率。因此，在预测的幅度方面，我们不会遇到不良预测，尽管我们至少需要检查预测是否给出+1 和-1 的混合。如果分类器只预测+1 或-1，那么这样的分类器可以被从进一步的分析中丢弃。对于回归问题，我们需要看看模型是否能在幅度上做出合理的预测，我们将简单地获得预测幅度的图，如下所示：

> windows(height=300,width=400)
> par(mfrow=c(2,5))
> plot.ts(SP_lm_pred,col=1)
> plot.ts(SP_rpart2_pred,col=2)
> plot.ts(SP_rpart4_pred,col=3)
> plot.ts(SP_rpart6_pred,col=4)
> plot.ts(SP_rpart8_pred,col=5)
> plot.ts(SP_nn2_pred,col=6)
> plot.ts(SP_nn3_pred,col=7)
> plot.ts(SP_nn4_pred,col=8)
> plot.ts(SP_nn5_pred,col=9)
> plot.ts(SP_svm_pred,col=10)

前一个代码块的结果如下所示：

简单平均

图 3：十个异构基学习器的预测简单图

我们可以看到，与具有两个或三个隐藏神经元的神经网络模型相关的预测没有产生预测上的变化。因此，我们将这两个模型从进一步的分析中删除。集成预测只是剩余八个模型预测的平均值：

> Avg_Ensemble_Prediction <- rowMeans(cbind(SP_lm_pred,SP_rpart2_pred,
+     SP_rpart4_pred,SP_rpart6_pred,
+                SP_rpart8_pred,SP_nn4_pred,SP_nn5_pred,SP_svm_pred))
> plot.ts(Avg_Ensemble_Prediction)

简单平均

图 4：住房数据集的集成预测

正如将简单投票扩展到加权投票一样，我们现在将探讨加权平均。

权重平均

在分类器的情况下，权重是从训练数据集的分类器的准确性中选择的。在这种情况下，我们需要像这样的统一度量。如果回归模型具有更小的残差方差，则更倾向于选择回归模型，我们将选择方差作为准确性的度量。假设弱基模型i的估计残差方差为权重平均。在集成神经网络的情况下，Perrone 和 Cooper（1993）声称可以使用以下方程获得i个弱基模型的最佳权重：

权重平均

由于比例常数无关紧要，我们将简单地用残差平方的平均值权重平均来代替。在这个方向上，我们将首先通过简单地计算简单平均情况下考虑的八个模型的mean(residuals(model)²)来获得（加上一个常数），如下所示：

> # Weighted Averaging
> SP_lm_sigma <- mean(residuals(SP_lm)²)
> SP_rp2_sigma <- mean(residuals(SP_rpart2)²)
> SP_rp4_sigma <- mean(residuals(SP_rpart4)²)
> SP_rp6_sigma <- mean(residuals(SP_rpart6)²)
> SP_rp8_sigma <- mean(residuals(SP_rpart8)²)
> SP_nn4_sigma <- mean(residuals(SP_nn4)²)
> SP_nn5_sigma <- mean(residuals(SP_nn5)²)
> SP_svm_sigma <- mean(residuals(SP_svm)²)

接下来，我们简单地实现权重公式权重平均，如下所示：

> sigma_sum <- SP_lm_sigma + SP_rp2_sigma + SP_rp4_sigma +
+   SP_rp6_sigma + SP_rp8_sigma + SP_nn4_sigma +
+   SP_nn5_sigma + SP_svm_sigma
> sigma_sum
[1] 20727111061
> SP_lm_wts <- SP_lm_sigma/sigma_sum
> SP_rp2_wts <- SP_rp2_sigma/sigma_sum
> SP_rp4_wts <- SP_rp4_sigma/sigma_sum
> SP_rp6_wts <- SP_rp6_sigma/sigma_sum
> SP_rp8_wts <- SP_rp8_sigma/sigma_sum
> SP_nn4_wts <- SP_nn4_sigma/sigma_sum
> SP_nn5_wts <- SP_nn5_sigma/sigma_sum
> SP_svm_wts <- SP_svm_sigma/sigma_sum

rowMeans和cbind函数简单地给出了加权平均预测：

> Weighted_Ensemble_Prediction <- rowMeans(cbind(SP_lm_wts*SP_lm_pred,
+                                           SP_rp2_wts*SP_rpart2_pred,
+                                           SP_rp4_wts*SP_rpart4_pred,
+                                           SP_rp6_wts*SP_rpart6_pred,
+                                           SP_rp8_wts*SP_rpart8_pred,
+                                           SP_nn4_wts*SP_nn4_pred,
+                                           SP_nn5_wts*SP_nn5_pred,
+                                           SP_svm_wts*SP_svm_pred))
> plot.ts(Weighted_Ensemble_Prediction)

前一个代码的输出如下所示：

权重平均

图 5：住房价格的加权平均预测

堆叠集成

在第一章中提供了一个堆叠回归的入门和激励示例，集成技术简介。在这里，我们将继续讨论一个尚未开发的回归问题的堆叠集成。

在堆叠集成中，几个弱模型的输出作为输入变量，以及用于构建早期模型的协变量，来构建一个堆叠模型。堆叠模型的形式可能是以下之一，或者可以是不同的模型。在这里，我们将简单地使用前几节中使用的八个回归模型作为弱模型。堆叠回归模型被选为梯度提升模型，并将给出原始输入变量和新模型的预测，如下所示：

> SP_lm_train <- predict(SP_lm,newdata=ht_imp)
Warning message:
In predict.lm(SP_lm, newdata = ht_imp) :
  prediction from a rank-deficient fit may be misleading
> SP_rpart2_train <- predict(SP_rpart2,newdata=ht_imp)
> SP_rpart4_train <- predict(SP_rpart4,newdata=ht_imp)
> SP_rpart6_train <- predict(SP_rpart6,newdata=ht_imp)
> SP_rpart8_train <- predict(SP_rpart8,newdata=ht_imp)
> SP_nn4_train <- predict(SP_nn4,newdata=ht_imp)
> SP_nn5_train <- predict(SP_nn5,newdata=ht_imp)
> SP_svm_train <- predict(SP_svm,newdata=ht_imp)
> 
> ht_imp2 <- cbind(ht_imp[,-69],SP_lm_train,SP_rpart2_train,SP_rpart4_train,
+                           SP_rpart6_train,SP_rpart8_train,SP_nn4_train,SP_nn5_train,
+                           SP_svm_train,ht_imp[,69])
> names(ht_imp2)[77] <- "SalePrice"
> SP_gbm <- gbm(hf,data=ht_imp2,distribution = "gaussian",n.trees=200)
> headtail(predict(SP_gbm,n.trees=100),20)
 [1] 180260.6 177793.3 181836.9 177793.3 191927.7 177793.3 191927.7 182237.3
 [9] 177793.3 177793.3 177793.3 191927.7 177793.3 187520.7 177793.3 177793.3
[17] 177793.3 177793.3 177793.3 177793.3 177908.2 177793.3 191927.7 177793.3
[25] 177793.3 177793.3 177793.3 191927.7 177793.3 177793.3 177793.3 191927.7
[33] 177793.3 177793.3 177793.3 177793.3 179501.7 191927.7 177793.3 177793.3

这结束了我们对堆叠集成回归的简单讨论。

摘要

在本章中，我们探讨了在分类问题背景下集成为什么有效。一系列详细程序说明了每个分类器必须比随机猜测更好的观点。我们考虑了所有分类器具有相同准确度、不同准确度和最后是完全任意准确度的场景。在随机森林和袋装方法背景下说明了多数和加权投票。对于回归问题，我们使用了不同的基学习器选择，并允许它们是异质的。在住房销售价格数据相关方面，我们展示了简单和加权平均方法。堆叠回归的简单说明最终结束了本章的技术部分。

在下一章中，我们将探讨集成诊断。

第八章. 集合诊断

在前面的章节中，我们发现集合方法非常有效。在前一章中，我们探讨了集合方法如何提高预测的整体准确性的场景。之前一直假设不同的基学习器之间是相互独立的。然而，除非我们有非常大的样本，并且基模型是使用一组不同观察值的学习者，否则这样的假设是非常不切实际的。即使我们有足够大的样本，相信分区是非重叠的，每个基模型都是建立在不同的分区上，每个分区都携带与任何其他分区相同的信息。然而，测试此类验证是困难的，因此我们需要采用各种技术来验证同一数据集上基模型的独立性。为此，我们将探讨各种不同的方法。本章将简要讨论集合诊断的必要性，并在下一节中介绍基模型多样性的重要性。对于分类问题，可以将分类器相互比较。然后我们可以进一步评估集合的相似性和准确性。在第三部分将介绍实现这一任务的统计测试。最初，将比较一个基学习器与另一个基学习器，然后我们将一次性查看集合中的所有模型。

本章将涵盖以下主题：

集合诊断
集合多样性
配对比较
评分者间一致性

技术要求

我们将在本章中使用以下库：

rpart

什么是集合诊断？

集成方法的力量在前几章中得到了展示。由决策树组成的集成是一个同质集成，这是第三章Bagging到第六章Boosting Refinements的主要内容。在第一章集成技术介绍和第七章通用集成技术中，我们简要介绍了堆叠集成。集成的一个中心假设是模型之间相互独立。然而，这个假设很少成立，我们知道相同的数据分区被反复使用。这并不意味着集成是坏的；我们有充分的理由在使用集成的同时预览集成应用中的担忧。因此，我们需要了解基础模型在预测上的相似程度以及整体上的相似程度。如果预测彼此很接近，那么我们可能需要在集成中使用这些基础模型。在这里，我们将为德国信用数据集构建逻辑回归、朴素贝叶斯、SVM 和决策树作为基础模型。分析和程序在这里略有重复，因为它是从早期章节继承过来的：

> load("../Data/GC2.RData")
> table(GC2$good_bad)
 bad good 
 300  700 
> set.seed(12345)
> Train_Test <- sample(c("Train","Test"),nrow(GC2),replace = 
+ TRUE,prob = c(0.7,0.3))
> head(Train_Test)
[1] "Test"  "Test"  "Test"  "Test"  "Train" "Train"
> GC2_Train <- GC2[Train_Test=="Train",]
> GC2_TestX <- within(GC2[Train_Test=="Test",],rm(good_bad))
> GC2_TestY <- GC2[Train_Test=="Test","good_bad"]
> GC2_TestY_numeric <- as.numeric(GC2_TestY)
> GC2_Formula <- as.formula("good_bad~.")
> p <- ncol(GC2_TestX)
> ntr <- nrow(GC2_Train) 
> nte <- nrow(GC2_TestX) 
> # Logistic Regression
> LR_fit <- glm(GC2_Formula,data=GC2_Train,family = binomial())
> LR_Predict_Train <- predict(LR_fit,newdata=GC2_Train,
+ type="response")
> LR_Predict_Train <- as.factor(ifelse(LR_Predict_Train>0.5,
+ "good","bad"))
> LR_Accuracy_Train <- sum(LR_Predict_Train==GC2_Train$good_bad)/
+ ntr
> LR_Accuracy_Train
[1] 0.78
> LR_Predict_Test <- predict(LR_fit,newdata=GC2_TestX,
+ type="response")
> LR_Predict_Test_Bin <- ifelse(LR_Predict_Test>0.5,2,1)
> LR_Accuracy_Test <- sum(LR_Predict_Test_Bin==
+ GC2_TestY_numeric)/nte
> LR_Accuracy_Test
[1] 0.757
> # Naive Bayes
> NB_fit <- naiveBayes(GC2_Formula,data=GC2_Train)
> NB_Predict_Train <- predict(NB_fit,newdata=GC2_Train)
> NB_Accuracy_Train <- sum(NB_Predict_Train==
+ GC2_Train$good_bad)/ntr
> NB_Accuracy_Train
[1] 0.767
> NB_Predict_Test <- predict(NB_fit,newdata=GC2_TestX)
> NB_Accuracy_Test <- sum(NB_Predict_Test==GC2_TestY)/nte
> NB_Accuracy_Test
[1] 0.808
> # Decision Tree
> CT_fit <- rpart(GC2_Formula,data=GC2_Train)
> CT_Predict_Train <- predict(CT_fit,newdata=GC2_Train,
+ type="class")
> CT_Accuracy_Train <- sum(CT_Predict_Train==
+ GC2_Train$good_bad)/ntr
> CT_Accuracy_Train
[1] 0.83
> CT_Predict_Test <- predict(CT_fit,newdata=GC2_TestX,
+ type="class")
> CT_Accuracy_Test <- sum(CT_Predict_Test==GC2_TestY)/nte
> CT_Accuracy_Test
[1] 0.706
> # Support Vector Machine
> SVM_fit <- svm(GC2_Formula,data=GC2_Train)
> SVM_Predict_Train <- predict(SVM_fit,newdata=GC2_Train,
+ type="class")
> SVM_Accuracy_Train <- sum(SVM_Predict_Train==
+ GC2_Train$good_bad)/ntr
> SVM_Accuracy_Train
[1] 0.77
> SVM_Predict_Test <- predict(SVM_fit,newdata=GC2_TestX,
+ type="class")
> SVM_Accuracy_Test <- sum(SVM_Predict_Test==GC2_TestY)/nte
> SVM_Accuracy_Test
[1] 0.754

在下一节中，我们将强调在集成中多样性的必要性。

集成多样性

在一个集成中，我们有许多基础模型——比如说有L个。对于分类问题，我们有基础模型作为分类器。如果我们有一个回归问题，我们有基础模型作为学习器。由于诊断仅在训练数据集上执行，我们将放弃训练和验证分区这一惯例。为了简单起见，在接下来的讨论中，我们将假设我们有N个观测值。L个模型意味着对于N个观测值中的每一个，我们都有L个预测，因此预测的数量是集成多样性。我们正是在这些预测中试图找到集成的多样性。集成的多样性取决于我们处理的问题类型。首先，我们将考虑回归问题。

数值预测

在回归问题的情况下，观测值的预测值可以直接与它们的实际值进行比较。我们可以很容易地看到哪个基础模型的预测值更接近观测的实际值，哪个离它更远。如果所有预测值都彼此接近，则基础模型不具有多样性。在这种情况下，任何一个预测可能就足够了。如果预测值表现出一些变化，通过使用平均值组合它们可能提供稳定性。在评估多样性时，了解集成预测与真实观测值有多接近也是非常重要的。

让我们考虑一个假设场景，在这个场景中我们有六个观测值，它们的实际值，三个基础学习器，学习器的预测，以及集成预测。以下表格提供了一个示例数据集，有助于您理解集成多样性的复杂性：

观测编号	实际值	E1	E2	E3	EP
1	30	15	20	25	20
2	30	40	50	60	50
3	30	25	30	35	30
4	30	28	30	32	30
5	30	20	30	40	30
6	30	10	15	65	30

表 1：六个观测值、三个基础学习器和集成

为了便于比较，表 1中所有观测值的真实值都保持在 30。六个观测/案例的集成预测范围从 10 到 65，而集成预测——基础学习器预测的平均值——范围从 20 到 50。为了理解特定观测值及其相关预测的集成多样性，我们将使用以下程序块可视化数据：

> DN <- read.csv("../Data/Diverse_Numeric.csv")
> windows(height=100,width=100)
> plot(NULL,xlim=c(5,70),ylim=c(0,7),yaxt='n',xlab="X-values",ylab="")
> points(DN[1,2:6],rep(1,5),pch=c(19,1,1,1,0),cex=2)
> points(DN[2,2:6],rep(2,5),pch=c(19,1,1,1,0),cex=2)
> points(DN[3,2:6],rep(3,5),pch=c(19,1,1,1,0),cex=2)
> points(DN[4,2:6],rep(4,5),pch=c(19,1,1,1,0),cex=2)
> points(DN[5,2:6],rep(5,5),pch=c(19,1,1,1,0),cex=2)
> points(DN[6,2:6],rep(6,5),pch=c(19,1,1,1,0),cex=2)
> legend(x=45,y=7,c("Actual","Model","Ensemble"),pch=c(19,1,0))
> axis(2,at=1:6,labels=paste("Case",1:6),las=1)

程序的解释在这里。代码的第一行从代码包文件夹中导入Diverse_Numeric.csv数据。windows (X11)函数在 Windows（Ubuntu）操作系统中设置一个新的图形设备。然后plot函数设置一个空白的绘图，并通过xlim和ylim指定坐标轴的范围。使用绘图和points函数将表 1中的每一行数据突出显示。选择pch需要进一步说明。如果我们选择pch为，例如，19、1和0，那么这意味着我们选择了填充的圆圈、圆圈和正方形。这三种形状分别表示实际值、模型预测和集成预测。轴命令帮助我们获得正确的标签显示。前面 R 代码块的结果是以下绘图：

数值预测

图 1：理解回归问题中的集成多样性

我们有六个观察结果，每个都被标记为案例。首先考虑案例 1。每个观察结果的实心圆圈表示实际值——在这种情况下是30，这在整个数据集中都是相同的。对于这个观察结果，集成预测是20。空白圆圈中描绘的是三个基础模型的预测值15、20和25。集成预测——基础学习器预测的平均值——是20，并用空白正方形表示。现在，这三个值分布得不太分散，这被解释为表明对于这个观察结果，集成多样性较低。因此，这是一个低多样性-低估计的例子。

在表 1的第二种情况下，三个预测值分布得很好，具有很高的多样性。然而，50的集成估计值与实际值30相差太远，我们将这种情况称为高多样性-低估计。案例 3和案例 4因此被视为低多样性-高估计，因为集成预测与实际值相符，且三个集成预测值彼此接近。案例 5在多样性和准确性之间取得了良好的平衡，因此我们可以将其标记为高多样性-高估计的例子。最后一种情况准确性良好，尽管多样性过高，以至于集成预测并不好。你可以参考 Kuncheva (2014)了解更多关于集成学习器多样性-准确性困境的细节。

我们将在下一节考虑分类问题的多样性-准确性问题。

类别预测

上一节讨论了回归问题中的多样性-准确性问题。在分类问题的案例中，我们可以清楚地标记分类器的预测是否与实际输出/标签匹配。此外，我们只有两种可能的预测：0 或 1。因此，我们可以比较两个分类器在所有观察结果上的接近程度。例如，对于分类器类别预测的两个可能结果和的两个可能结果，对于给定的观察结果有四种可能的场景：

预测标签为 1；预测为 1
预测标签为 1；预测为 0
预测标签为 0；预测为 1
预测标签为 0；预测为 0

在情景 1 和 4 中，两个分类器同意彼此，而在 2 和 3 中，它们不同意。如果我们有N个观测，每个观测用两个模型预测将落入上述四种情景之一。在我们考虑两个或更多模型的正式一致性或不一致性度量之前，我们将在接下来的讨论中考虑两个更简单的案例。

有一个流行的说法，如果两个人总是同意，其中一个人是不必要的。这与分类器工作的方式相似。同样，假设一对鹅被认为是非常忠诚的；它们彼此相伴，共同面对问题。现在，如果我们有两个模型在所有观测中都以与这些鹅相同的方式表现，那么多样性就会永远丧失。因此，在任何给定的集成情景中，我们需要消除鹅对，只保留其中之一。假设我们有一个L预测的矩阵，其中列对应分类器，行对应N个观测。在这种情况下，我们将定义一个名为GP的函数，以及鹅对的缩写，它将告诉我们哪些分类器有一个鹅对分类器在所有观测中与它们一致：

> # Drop the Geese Pair
>GP<- function(Pred_Matrix) {
+   L<- ncol(Pred_Matrix) # Number of classifiers
+   N<- nrow(Pred_Matrix)
+   GP_Matrix <- matrix(TRUE,nrow=L,ncol=L)
+   for(i in 1:(L-1)){
+     for(j in (i+1):L){
+       GP_Matrix[i,j] <- ifelse(sum(Pred_Matrix[,i]==Pred_Matrix[,j])==N,
+                                TRUE,FALSE)
+       GP_Matrix[j,i] <- GP_Matrix[i,j]
+     }
+   }
+   return(GP_Matrix)
+ }

鹅对GP函数是如何工作的？我们向这个函数输入一个矩阵预测作为输入，其中列对应分类器，行对应观测。这个函数首先创建一个逻辑矩阵，其阶数为类别预测，默认逻辑值为TRUE。由于一个分类器显然会与自己一致，我们接受默认值。此外，由于分类器与在相同的方式上同意/不同意，我们利用这一事实通过对称关系计算下矩阵。在两个嵌套循环中，我们比较一个分类器的预测与另一个分类器的预测。ifelse函数检查一个分类器的所有预测是否与另一个分类器匹配，如果对于单个观测条件不成立，我们说正在考虑的两个分类器不是鹅对，或者它们至少在某个场合上不同意：

接下来，将GP函数应用于为分类问题设置的 500 个分类器。CART_Dummy数据集来自RSADBE包。CART_DUMMY数据集和相关问题描述可以在 Tattar（2017）的第九章中找到。我们从这个相同的来源改编了代码和结果输出：

> data(CART_Dummy)
> CART_Dummy$Y <- as.factor(CART_Dummy$Y)
> attach(CART_Dummy)
> windows(height=100,width=200)
> par(mfrow=c(1,2))
> plot(c(0,12),c(0,10),type="n",xlab="X1",ylab="X2")
> points(X1[Y==0],X2[Y==0],pch=15,col="red")
> points(X1[Y==1],X2[Y==1],pch=19,col="green")
> title(main="A Difficult Classification Problem")
> plot(c(0,12),c(0,10),type="n",xlab="X1",ylab="X2")
> points(X1[Y==0],X2[Y==0],pch=15,col="red")
> points(X1[Y==1],X2[Y==1],pch=19,col="green")
> segments(x0=c(0,0,6,6),y0=c(3.75,6.25,2.25,5),
+          x1=c(6,6,12,12),y1=c(3.75,6.25,2.25,5),lwd=2)
> abline(v=6,lwd=2)
> title(main="Looks a Solvable Problem Under Partitions")

如程序所示，我们这里有三个变量：X1、X2和Y。由Y表示的变量是一个二元变量——一个类别用绿色表示，另一个用红色表示。使用X1和X2变量提供的信息，目标是预测Y的类别。红色和绿色点交织在一起，因此单一线性分类器不足以在这里将红色和绿色分开。然而，如果我们通过X1和X2递归地划分数据空间，如图 2 右侧的最终图所示，那么红色和绿色看起来是可分离的。前面的 R 代码块产生以下图表：

类别预测

图 2：一个典型的分类问题

为CART_DUMMY数据集设置了包含500棵树的随机森林。固定的种子确保了在任何执行中输出的可重复性。使用拟合的随机森林，我们接下来使用500棵树预测所有观测值的输出。type="class"和predict.all=TRUE选项是此代码块的核心。然后，将GP函数应用于500棵树的预测矩阵。请注意，GP矩阵的对角元素始终为TRUE。因此，如果有任何分类器与所有观测值完全一致，该单元格的值将是TRUE。如果行和超过计数 2，则该分类器有鹅分类器。以下代码捕获了整个计算过程：

> CD <- CART_Dummy 
> CD$Y <- as.factor(CD$Y)
> set.seed(1234567)
> CD_RF <- randomForest(Y~.,data=CD,ntree=500)
> CD_RF_Predict <- predict(CD_RF,newdata=CD,
+                           type="class",predict.all=TRUE)
> CD_RF_Predict_Matrix <- CD_RF_Predict$individual
> CD_GP <- GP(CD_RF_Predict_Matrix)
> CD_GP[1:8,1:8]
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[2,] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[3,] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
[4,] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
[5,] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
[6,] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
[7,] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
[8,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
> rowSums(CD_ST)
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1
[186] 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
[223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
[260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
[297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1
[371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1

[482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

读者应注意的是，前面输出中 2 的粗体和较大字体不是由 R 给出的。这是由处理文本内容的软件修改的。因此，我们有很多具有鹅分类器的分类器，它们与各自的预测相匹配。使用which函数，我们首先找到所有满足条件的分类器索引，然后，通过应用which函数到CD_GP矩阵的行，我们得到相关的鹅分类器：

> which(rowSums(CD_GP)>1)
 [1]  21  42 176 188 206 221 256 278 290 363 365 385 424 442
> which(CD_GP[21,]==TRUE)
[1]  21 188
> which(CD_GP[42,]==TRUE)
[1]  42 290
> which(CD_GP[176,]==TRUE)
[1] 176 363
> which(CD_GP[206,]==TRUE)
[1] 206 256
> which(CD_GP[221,]==TRUE)
[1] 221 278
> which(CD_GP[365,]==TRUE)
[1] 365 424
> which(CD_GP[385,]==TRUE)
[1] 385 442

运行前面的代码后，我们能够识别与分类器关联的鹅分类器。我们可以选择移除鹅对中的任何一个成员。在下一个示例中，我们将应用此方法到德国信用数据。程序尝试以下方式识别鹅分类器：

> set.seed(12345)
> GC2_RF3 <- randomForest(GC2_Formula,data=GC2_Train,mtry=10,
+                         parms = list(split="information",
+                                      loss=matrix(c(0,1,1000,0),byrow = TRUE,nrow=2)),
+                         ntree=1000)
> GC2_RF_Train_Predict <- predict(GC2_RF3,newdata=GC2_Train,
+                                 type="class",predict.all=TRUE)
> GC2_RF_Train_Predict_Matrix <- GC2_RF_Train_Predict$individual
> GC2_GP <- GP(GC2_RF_Train_Predict_Matrix)
> rowSums(GC2_GP)
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [37] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[973] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> which(rowSums(GC2_GP)>1)
integer(0)

由于没有分类器有对应的鹅分类器，我们不需要消除任何树。

在 Kuncheva (2014)，第 112 页，有一个有用的度量标准，称为神谕输出。接下来，我们正式定义这个量。记住，我们有L个分类器和N个观测值。标签的原始/实际值用类别预测表示。我们将使用分类器 j 表示第 i 个预测值，用表示。

Oracle 输出：如果预测值 Class prediction 等于，则 Oracle 输出被定义为1；否则，它被定义为0。用数学术语来说，Oracle 输出使用以下数学表达式给出：

Class prediction

那么，Oracle 输出和预测之间有什么区别？预测包括数据的标签，标签可能是 1/0、GOOD/BAD、+1/-1、YES/NO 或其他二进制标签对。此外，在二进制标签的情况下，预测为 1 并不一定意味着原始值是 1；它也可能是 0。如果预测 1 为 1 或 0 为 0，Oracle 输出取值为 1；否则，它取值为 0。使用 Oracle 输出的一个后果是，分类器中 1 的比例将给我们提供分类器的准确度。

我们现在将创建一个名为Oracle的 R 函数，该函数在预测矩阵和实际标签作为输入时将给出 Oracle 输出。之后，我们将计算分类器的准确度：

> # Oracle Output
> Oracle <- function(PM,Actual){
+   # PM = Prediction Matrix, Actual = the true Y's
+   OM <- matrix(0,nrow=nrow(PM),ncol=ncol(PM))
+   for(i in 1:ncol(OM)) {
+     OM[,i] <- as.numeric(PM[,i]==Actual)
+   }
+   return(OM)
+ }
> GC_Oracle <- Oracle(PM=GC2_RF_Train_Predict$individual,
+                     Actual=GC2_Train$good_bad)
> colSums(GC_Oracle)/nrow(GC_Oracle)
   [1] 0.872 0.884 0.859 0.869 0.866 0.878 0.888 0.872 0.869 0.875 0.885 0.869
  [13] 0.881 0.866 0.879 0.856 0.870 0.869 0.857 0.870 0.878 0.868 0.886 0.892
  [25] 0.881 0.863 0.866 0.856 0.886 0.876 0.873 0.879 0.875 0.885 0.872 0.872

[973] 0.860 0.873 0.869 0.888 0.863 0.879 0.882 0.865 0.891 0.863 0.878 0.879
 [985] 0.878 0.869 0.856 0.872 0.889 0.881 0.868 0.881 0.884 0.854 0.882 0.882
 [997] 0.862 0.884 0.873 0.885

Oracle 矩阵帮助我们获得分类器的准确度。在下一节中，我们将讨论一些有助于我们了解分类器之间距离的度量。

配对度量

在本节中，我们将提出一些度量两个分类器之间一致性的方法。目的是固定两个分类器之间的一致性/不一致性概念，然后在下一节中将该概念应用到集成分类器的整体分类器中。如果配对度量和是具有预测的分类器模型，那么我们可以获得一个表格，它给出以下内容：

预测为 1；预测它为 1
预测为 1；预测它为 0
预测为 0；预测它为 1
预测为 0；预测它为 0

在N个观测值中的信息可以以表格形式表示，如下所示：

	M1 预测 1	M1 预测 0
M2 预测 1	n11	n10
M2 预测 0	n01	n00

表 2：两个分类器/评分者的列联表

前一个表的对角线元素显示两个模型/分类器之间的一致性，而离对角线元素显示不一致性。这些模型有时被称为评分者。频率表也被称为列联表。使用这种设置，我们现在将讨论一些有用的一致性度量。比较被称为成对度量，因为我们只分析了一对分类器。

不一致性度量

两个分类器/评分者之间的一致性度量定义为以下公式：

不一致性度量

我们现在将定义一个DM函数，该函数接受两个分类器的预测。该函数首先为预测准备列联表。不一致性度量的计算是直接的，如下代码块所示：

> # Disagreement Measure
> DM <- function(prediction1,prediction2){
+   tp <- table(prediction1,prediction2)
+   Diss <- (tp[1,2]+tp[2,1])/length(prediction1)
+   return(Diss)
+ }

在第一部分，我们根据逻辑回归模型、朴素贝叶斯、SVM 和分类树对德国信用数据进行了预测。现在我们将 DM 函数应用于这些预测，看看这些分类器之间有多大的不一致性：

> DM(LR_Predict_Train,NB_Predict_Train)
[1] 0.121
> DM(LR_Predict_Train,CT_Predict_Train)
[1] 0.154
> DM(LR_Predict_Train,SVM_Predict_Train)
[1] 0.153
> DM(NB_Predict_Train,CT_Predict_Train)
[1] 0.179
> DM(NB_Predict_Train,SVM_Predict_Train)
[1] 0.154
> DM(CT_Predict_Train,SVM_Predict_Train)
[1] 0.167

由于我们有四个分类器，因此将会有 3 + 2 + 1 = 6 个成对比较。朴素贝叶斯和分类树之间的不一致性最大，而逻辑回归和朴素贝叶斯分类器之间的一致性最小。可以使用 DM 度量来轻松地获得两个模型的不一致性。

Yule 或 Q 统计量

Yule 系数是一致性的度量，当其值接近于零时，它将给出两个评分者之间的一致性。Yule 度量使用以下公式给出：

Yule 的或 Q 统计量

Q 统计量的值在相关系数的范围内——即 Yule 的或 Q 统计量。因此，如果 Q 值接近于 1，这意味着两个度量几乎总是相互一致，而接近于-1 的值意味着两个模型预测的是相反的。当 Q 值接近于 0 时，这意味着两个评分者之间有非常弱的相关性。以下代码块创建并应用了一个Yule函数，用于不同的模型预测：

> # Q-statistic 
> Yule <- function(prediction1,prediction2){
+   tp <- table(prediction1,prediction2)
+   Yu <- (tp[1,1]*tp[2,2]-tp[1,2]*tp[2,1])/(tp[1,1]*tp[2,2]+tp[1,2]*tp[2,1])
+   return(Yu)
+ }
> Yule(LR_Predict_Train,NB_Predict_Train)
[1] 0.949
> Yule(LR_Predict_Train,CT_Predict_Train)
[1] 0.906
> Yule(LR_Predict_Train,SVM_Predict_Train)
[1] 0.98
> Yule(NB_Predict_Train,CT_Predict_Train)
[1] 0.865
> Yule(NB_Predict_Train,SVM_Predict_Train)
[1] 0.985
> Yule(CT_Predict_Train,SVM_Predict_Train)
[1] 0.912

朴素贝叶斯预测和 SVM 预测之间的一致性最高。请注意，如果我们取不一致性度量的补数并使用以下代码轻松执行，我们得到以下一致性的度量：

> 1-DM(LR_Predict_Train,NB_Predict_Train)
[1] 0.879
> 1-DM(LR_Predict_Train,CT_Predict_Train)
[1] 0.846
> 1-DM(LR_Predict_Train,SVM_Predict_Train)
[1] 0.847
> 1-DM(NB_Predict_Train,CT_Predict_Train)
[1] 0.821
> 1-DM(NB_Predict_Train,SVM_Predict_Train)
[1] 0.846
> 1-DM(CT_Predict_Train,SVM_Predict_Train)
[1] 0.833

然而，这项分析表明，逻辑回归和朴素贝叶斯评分者之间的一致性最高。因此，我们注意到输出和比较可能会导致不同的结论。也可以计算两个评分者之间的相关系数；我们将在下一部分介绍。

科亨统计量

科亨统计量首次出现在 1960 年。它基于两个评分者因偶然或巧合而达成一致的概率。两个评分者达成一致的概率如下所示：

科亨统计量

然而，随机或偶然达成一致的概率如下所示：

科亨统计量

使用和的定义，科亨统计量如下定义：

科亨统计量

科亨的 Kappa 值也可以是负数。如果其值为 1，这意味着评分者完全一致。0 的值表示一致仅是偶然的，负值表示偶然的一致性低于预期。首先，在以下代码中创建了 R 函数Kappa：

> # Kappa-statistic 
> # Cohen's Statistic
> Kappa <- function(prediction1, prediction2){
+   tp <- table(prediction1,prediction2)
+   a <- tp[1,1]; b <- tp[2,1]; c <- tp[1,2]; d <- tp[2,2]
+   n <- length(prediction1)
+   theta1 <- (a+d)/n
+   theta2 <- (((a+b)*(a+c))+((c+d)*(b+d)))/n²
+   kappa <- (theta1-theta2)/(1-theta2)
+   return(kappa)
+ }

编码部分是公式的清晰实现，对a、b、c、d、theta1和theta2的选择已经做出，以便代码易于理解和遵循。接下来，我们将预测应用于德语训练数据集：

> Kappa(LR_Predict_Train,NB_Predict_Train)
[1] 0.69
> Kappa(LR_Predict_Train,CT_Predict_Train)
[1] 0.592
> Kappa(LR_Predict_Train,SVM_Predict_Train)
[1] 0.524
> Kappa(NB_Predict_Train,CT_Predict_Train)
[1] 0.53
> Kappa(NB_Predict_Train,SVM_Predict_Train)
[1] 0.525
> Kappa(CT_Predict_Train,SVM_Predict_Train)
[1] 0.453

再次，逻辑回归和朴素贝叶斯预测的一致性是最高的。我们现在转向最终的分歧度量。

双误度量

在网球中，双误指的是发球失败的情况。发球者有两个机会发球正确，如果他们没有做到，则该分被判给对手。双误度量发生在两个分类器都做出错误预测时：

双误度量

显然，我们需要 DF 尽可能低，接近 0。这个函数易于理解，因此这将被留给读者作为练习来跟随。以下代码给出了双误度量的 R 函数及其应用：

> # Double-fault Measure
> Double_Fault <- function(prediction1,prediction2,actual){
+   DF <- sum((prediction1!=actual)*(prediction2!=actual))/
+         length(actual)
+   return(DF)
+ }
> Double_Fault(LR_Predict_Train,NB_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.166
> Double_Fault(LR_Predict_Train,CT_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.118
> Double_Fault(LR_Predict_Train,SVM_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.148
> Double_Fault(NB_Predict_Train,CT_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.709
> Double_Fault(NB_Predict_Train,SVM_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.154
> Double_Fault(CT_Predict_Train,SVM_Predict_Train,
+ GC2_Train$good_bad)
[1] 0.116

读者应使用双重错误度量来识别最佳一致性。

练习：在多标签（超过两个类别）的情况下，本节讨论的度量扩展变得繁琐。相反，可以使用 Oracle 矩阵并重复这些度量。读者应将这些措施应用于 Oracle 输出。

到目前为止讨论的方法仅适用于一对分类器。在下一节中，我们将测量集成中所有分类器的多样性。

交互一致性

在前节讨论的集成分类器度量简单扩展是计算所有可能的集成对度量，然后简单地平均这些值。这项任务构成了下一个练习。

练习：对于所有可能的集成对组合，计算不一致度量、Yule 统计量、相关系数、Cohen's kappa 和双重错误度量。完成这些后，获得比较的平均值，并将它们报告为集成多样性。

在这里，我们将提出多样性的替代度量，并使用熵度量启动讨论。在本节的所有讨论中，我们将使用 Oracle 输出。

熵度量

你可能记得，我们根据来表示 Oracle 输出。对于特定实例，如果错误分类该实例的分类器数量为，则集成具有最大多样性。这意味着的值为 0，其余的、、的值为 1。然后，集成熵度量定义为以下：

熵度量

熵度量 E 的值位于单位区间内。如果 E 值接近 0，这意味着集成中没有多样性，而接近 1 的值意味着多样性达到最高可能水平。给定 Oracle 矩阵，我们可以轻松地计算熵度量，如下所示：

> # Entropy Measure
> # Page 250 of Kuncheva (2014)
> Entropy_Measure <- function(OM){
+   # OM = Oracle Matrix
+   N <- nrow(OM); L <- ncol(OM)
+   E <- 0
+   for(i in 1:N){
+     E <- E+min(sum(OM[i,]),L-sum(OM[i,]))
+   }
+   E <- 2*E/(N*(L-1))
+   return(E)
+ }
> Entropy_Measure(GC_Oracle)
[1] 0.255

通过在德国信用数据集的集成上应用Entropy_Measure，我们可以看到熵度量值为0.255。由于熵度量值没有接近 0，随机森林集成表现出多样性。然而，它也远离 1，这表明存在多样性。然而，没有临界值或测试来解释多样性是否太低，甚至太高。

Kohavi-Wolpert 度量

Kohavi-Wolpert 度量基于预测的方差为 1 或 0。它基于分类器错误率的分解公式。对于二元问题或使用 Oracle 输入时，方差与 Gini 指数相同。这如下公式给出：

Kohavi-Wolpert 度量

Kohavi-Wolpert 度量是所有观察到的方差的平均值。通过使用 Oracle 矩阵给出的预测概率，或者作为拟合对象的副产品，我们可以获得方差，然后对观察到的方差进行平均。现在创建了一个 R 函数，并将其应用于德国信用数据的一些预测，如下所示：

> # Kohavi-Wolpert variance 
> # Using the predicted probability
> KW <- function(Prob){
+   N <- nrow(Prob)
+   kw <- mean(1-Prob[,1]²-Prob[,2]²)/2
+   return(kw)
+ }
> GC2_RF_Train_Predict_Prob <- predict(GC2_RF3,newdata=GC2_Train,
+                                 type="prob",predict.all=TRUE)
> GC2_RF_Train_Prob <- GC2_RF_Train_Predict_Prob$aggregate
> KW(GC2_RF_Train_Prob)
[1] 0.104

Kohavi-Wolpert 度量也可以通过 Oracle 输出获得。我们定义一个数学实体，用来计算正确分类观察到的分类器的数量如下：

Kohavi-Wolpert 度量

被正确预测的概率如下：

Kohavi-Wolpert 度量

使用这些概率，我们可以获得方差如下：

Kohavi-Wolpert 度量

此方法通过以下代码使用KW_OM函数实现：

> # Using the Oracle matrix
> KW_OM<- function(OM){
+   # OM is the oracle matrix
+   N <- nrow(OM); L <- ncol(OM)
+   kw <- 0
+   for(i in 1:N){
+     lz <- sum(OM[i,])
+     kw <- kw + lz*(L-lz)
+   }
+   kw <- kw/(N*L²)
+   return(kw)
+ }
> KW_OM(GC_Oracle)
[1] 0.104

从这里我们可以看出，这两种方法给出了相同的结果。很明显，在随机森林构建之后，我们没有看到太多的多样性。

集成不一致度量

两个分类器之间的不一致度量可以定义为以下：

集成不一致度量

集成的不一致度量如下所示：

集成不一致度量

Kohavi-Wolpert 度量与不一致度量之间的关系如下：

集成不一致度量

下一个 R 代码块展示了使用 Oracle 输出实现 Kohavi-Wolper 度量的方法如下：

> # Disagreement Measure OVerall on Oracle Matrix
> DMO <- function(OM){
+   # OM is the oracle matrix
+   N <- nrow(OM); L <- ncol(OM)
+   dmo <- 0
+   for(i in 1:L){
+     for(j in c(c(1:L)[c(1:L)!=i])){
+       dmo <- dmo + sum((OM[,i]-OM[,j])²)/N
+     }
+   }
+   dmo <- dmo/(L*(L-1))
+   return(dmo)
+ }
> DM_GC <- DMO(OM=GC_Oracle)
> DM_GC
[1] 0.208
> KW(GC_Oracle)
[1] 0.104
> DM_GC*999/2000
[1] 0.104

再次，我们没有看到在集成中表现出太多的多样性。现在我们将继续探讨集成多样性的最终度量。

评分者间一致性度量

在对 Oracle 输出介绍的讨论中，我们展示了如何轻松地使用它来获取分类器的准确率。分类器准确率的平均值定义为平均个体分类准确率，并用评分者间一致性度量表示。评分者间一致性的度量由以下定义：

评分者间一致性度量

此度量与 Kohavi-Wolpert 度量相关，如下所示：

评分者间一致性度量

通过以下代码块，我们可以理解前面关系的实现：

> Avg_Ensemble_Acc <- function(Oracle){
+   return(mean(colSums(GC_Oracle)/nrow(GC_Oracle)))
+ }
> Avg_Ensemble_Acc(GC_Oracle)
[1] 0.872
> Kappa <- function(Oracle){
+   pbar <- Avg_Ensemble_Acc(Oracle)
+   AvgL <- 0
+   N <- nrow(Oracle); L <- ncol(Oracle)
+   for(i in 1:N){
+     lz <- sum(Oracle[i,])
+     AvgL <- AvgL + lz*(L-lz)
+   }
+   Avgl <- AvgL/L
+   kappa <- 1-Avgl/(N*(L-1)*pbar*(1-pbar))
+   return(kappa)
+ }
> Kappa(GC_Oracle)
[1] 0.0657
> 1-DM_GC/(2*Avg_Ensemble_Acc(GC_Oracle)*(1-
+ Avg_Ensemble_Acc(GC_Oracle)))
[1] 0.0657

这就结束了我们对集成一致性的讨论。

摘要

集成方法被发现对于分类、回归和其他相关问题是十分有效的。任何统计和机器学习方法都必须始终伴随着适当的诊断。所有基础模型之间相互独立这一假设是集成方法成功的关键。然而，这种独立性条件很少得到满足，尤其是因为基础模型是建立在相同的数据集之上的。我们以最简单的度量方法——鹅对法——开始了这一章节。通过这种方法，我们实际上是在寻找在所有时刻都达成一致的模型。如果这样的模型存在于集成中，那么移除其中之一会更安全。在拥有大量数据集和高数量变量的情况下，确实可能没有任何基础模型与其他模型使用相同的语言。然而，我们仍然需要检查它们是否相等。考虑到这一点，我们首先提出了仅比较两个基础模型的措施。不同的措施可能导致相互矛盾的结论。然而，这通常并不是问题。随后，我们将成对比较的概念扩展到了整个集成基础模型。虽然我们发现我们的基础模型并不太多样化，但在此也要注意，大多数值都远离边界值 0。当我们对集成进行诊断并发现值等于零时，那么很明显，基础模型并没有提供任何形式的多样性。在下一章中，我们将探讨回归数据的专门主题。

分类器 1	分类器 2	分类器 3
真	真	真
假	真	真
真	假	真
假	假	真
真	真	假
假	真	假
真	假	假
假	假	假

分类器 1	分类器 2	分类器 3	投票数
真	真	真	3
假	真	真	2
真	假	真	2
假	假	真	1
真	真	假	2
假	真	假	1
真	假	假	1
假	假	假	0

观测编号	实际值	E1	E2	E3	EP
1	30	15	20	25	20
2	30	40	50	60	50
3	30	25	30	35	30
4	30	28	30	32	30
5	30	20	30	40	30
6	30	10	15	65	30

分类器 1	分类器 2	分类器 3
真	真	真
假	真	真
真	假	真
假	假	真
真	真	假
假	真	假
真	假	假
假	假	假

分类器 1	分类器 2	分类器 3	投票数
真	真	真	3
假	真	真	2
真	假	真	2
假	假	真	1
真	真	假	2
假	真	假	1
真	假	假	1
假	假	假	0

观测编号	实际值	E1	E2	E3	EP
1	30	15	20	25	20
2	30	40	50	60	50
3	30	25	30	35	30
4	30	28	30	32	30
5	30	20	30	40	30
6	30	10	15	65	30

R-集成学习实用指南-二-