集成学习

结合R实例

Posted by T.L on December 20, 2017

融合方法

有人把它称为机器学习中的“屠龙刀”,非常万能且有效,集成模型是一种能在各种的机器学习任务上提高准确率的强有力技术,集成算法往往是很多数据竞赛关键的一步,能够很好地提升算法的性能。哲学思想为“三个臭皮匠赛过诸葛亮

现实生活中,我们经常会通过投票,开会等方式,以做出更加可靠的决策。集成学习就与此类似。集成学习就是有策略的生成一些基础模型,然后有策略地把它们都结合起来以做出最终的决策。集成学习又叫多分类器系统。

集成方法是由多个较弱的模型集成模型组,一般的弱分类器可以是DT, SVM, NN, KNN等构成。其中的模型可以单独进行训练,并且它们的预测能以某种方式结合起来去做出一个总体预测。该算法主要的问题是要找出哪些较弱的模型可以结合起来,以及如何结合的方法。这是一个非常强大的技术集,因此广受欢迎。

目前,有三种常见的集成学习框架:bagging,boosting和stacking。

算法鸟瞰

Bagging:基于数据随机重抽样的分类器构建方法。从训练集抽样组成每个基模型所需要的子训练集,对所有基模型预测的结果进行综合产生最终的预测结果 : bagging boosting:训练过程为阶梯状,基模型按次序一一进行训练(实现上可以做到并行),基模型的训练集按照某种策略每次都进行一定的转化。对所有基模型预测的结果进行线性综合产生最终的预测结果:) boosting stacking:将训练好的所有基模型对训练基进行预测,第j个基模型对第i个训练样本的预测值将作为新的训练集中第i个样本的第j个特征值,最后基于新的训练集进行训练。同理,预测的过程也要先经过所有基模型的预测形成新的测试集,最后再对测试集进行预测 stacking

关于基础分类器结果整合的主要方式

关于回归

  1. 简单平均(uniformly),就是取各个分类器结果的平均值。 \(G(x)=\frac{1}{T}\sum g_t(x)\)
  2. 加权平均(non-uniform) \(G(x)=\frac{1}{T}\sum \alpha _tg_t(x),\alpha _t \ge 0\)
  3. 条件平均(non-uniform) \(G(x)=\frac{1}{T}\sum q_t(x) g_t(x) ,q_t(x)\ge 0\)

关于分类

和上面相同,公式同步变为: \(G(x)=sign(\sum 1g_t(x))\)

集成框架

集成学习的关键有两点:

  1. 如何构建具有差异性的分类器
  2. 选择多少个模型
  3. 如何多这些分类器的结果进行整合。

多样性

基学习器之间的性能要有较大的差别,否则集成效果会很不好,也就是要保证多样性(Diversity),基学习器可以是DT, SVM, NN, KNN等,也可以是训练过程不同的同一种模型,例如不同的参数,不同的训练集,或者不同的特征选择等。

为了说明多样性可以提升,以下面三个模型为例:

1
2
3
1111111100 = 80% accuracy 
1111111100 = 80% accuracy 
1011111100 = 70% accuracy. 

下面我们比较3个性能相对较差,但是高度不相关的模型:

1
2
3
1111111100 = 80% accuracy 
0111011101 = 70% accuracy 
1000101111 = 60% accuracy 

在比赛中挑选相关度低的模型进行组合。

模型数量

我们有3个正确率为60%的二分类器记为A,B,C。你可以将这些分类器视为伪随机数产生器,以60%的概率产生”1”,40%的概率产生”0”。 通过集成正确率多少?

集成框架

bagging

bagging实际上是一种把多个在同一数据集上训练的模型组合起来的直观过程,它对分类模型采用多数票,对回归模型使用平均值。 以二元分类为例: 输入: ·data:包含输入特征和带有二元输出标签的一个列的输入数据框 ·M:一个代表我们要训练的模型数量的整数

输入: ·data:包含输入特征和带有二元输出标签的一个列的输入数据框
·M:一个代表我们要训练的模型数量的整数 输出: models: 一个训练好的二元分类模型 方法:

1)创建一个大小为n的有放回的随机样本,其中的n是原始数据集里的观测数据数量。这意味着原始训练集里的某些观测数据会重复使用,而某些则根本不会被选取。这个过程被称为自助法(bootstrapping)、自助抽样(bootstrap sampling)或自助重抽样(bootstrap resampling)。
2)利用这些抽样的数据集训练一个分类模型。典型情况下,我们在该步骤不会选择用于约简过拟合的正则化或收缩方法,因为最终采用的聚合过程会用来消减过拟合。
3)对抽样数据集里的每条观测数据,记录模型给它分配的类别。
4)重复这个过程M次,训练M个模型。
5)对于原始训练数据集里的每条观测数据,通过不同模型的多数票计算预测出的类别。例如,假设M=61,通过自主抽样,某条观测数据出现在其中50个模型的训练数据中。如果这些模型中有37个对这条观测数据预测了类别1,13个预测了类别-1,根据多数票原则,整体的预测就是类别1。
6)利用训练集提供的标签来计算模型的精确度。

有个值得探讨的问题是:“每次进行有放回的抽样时,我们能得到多少不同的观测数据?

在步骤1里,我们构建了一个训练数据的样本,用来训练模型。在原始数据集里,我们把在该过程的某次迭代没有选中的观测数据称为袋外(out-of-bag,OOB)观测数据。那么,这些观测数据在那一次迭代中就没有用在模型的训练中。因此,我们实际上可以利用OOB观测数据来记录某个模型的精确度,而不必每一步都依靠那些用来训练模型的观测数据。

由于是有放回的抽样,我们可以认为这些数据是独立同分布的,这些观测数据均值的方差是$\sigma^2/n$。

代表算法: 随机森林

boosting

对于分类来说,增强起作用的原理是通过在训练数据上构建一个模型,然后衡量在那些训练数据上的分类精确度。被模型错误分类的每条观测数据会得到一个比正确分类的观测数据更大的权值,然后这个模型会用这些新的权值重新训练。这个过程会重复很多次,每次都会根据上一次迭代时是否正确分类的结果来调整每条观测数据的权值。

为了应对过拟合问题,集成分类器的构建方式是取在这个序列里训练的所有模型的加权平均,其中的权值通常是和每个模型的分类精确度成比例的。因为我们用到了全部训练数据,所以就没有袋外观测数据,因此每个模型的精确度都是用训练数据本身来衡量的。而用增强进行回归的做法,通常是根据预测值和标签值之间距离的某种衡量指标来调整观测数据的权值。

用于二元分类的自适应增强
输入:
·data:包含了输入特征和一个二元输出标签列的输入数据框
·M:表示我们需要训练的模型数量的一个整数
输出:
·model:由M个训练好的模型组成的序列
·alpha:包含M个模型权值的向量

方法: 1)初始化一条观测数据权值的向量w,它的长度为n,初始元素为wi=1/n。该向量会在每次迭代中更新。
2)使用观测数据权值的当前值和训练集里的所有数据训练一个分类器模型Gm。
3)计算权值的误差率,用所有错误分类观测数据乘以其权值然后求和,再除以权值向量的总和。根据用xi作为一条观测数据并用yi作为其标签的惯例,我们可以利用下面的公式来表示:
\(err_m=\frac{\sum \omega_iI(y_i\ne G_m(x_i)) }{\sum \omega_i}\) 4)设置该模型的模型权值$\alpha_m$为精确度和误差率之比的对数。用公式表达为: \(\alpha_m=\frac{1}{2}ln(\frac{1-err_m}{err_m})\) 5)为下一次迭代更新观测数据权值向量w。错误分类的观测数据的权值会乘以$e^{\alpha m}$,从而增大它们在下一次迭代时的权值。正确分类的观测数据的权值会乘以$e^{-\alpha m}$,从而减小它们在下一次迭代时的权值。 6)对权值向量重新进行归一化,使权值的总和等于1。
7)把步骤2到6重复M次,产生M个模型。
8)定义我们的集成分类器为所有增强模型输出的加权总和:
\(G(x)=sign(\sum \alpha_m G_m(x))\)

stacking

参照CCF大赛搜狗用户画像总结

代码实践

bagging

针对Statlog心脏病数据集的逻辑回归分类器

心脏病数据集
该数据包含了270个可能心脏有问题的病人的观测数据。其中有120个病人表现出心脏有问题,因此两个类别的划分是比较平均的。我们的任务是根据病人的基本情况和一系列的医学检查来预测病人是否患有心脏病。

1
2
3
4
5
6
7
8
library(shiny)
heart <- read.table("./data/heart.dat", quote="\"")
names(heart) <- c("AGE", "SEX", "CHESTPAIN", "RESTBP", "CHOL", "SUGAR", "ECG", "MAXHR", "ANGINA", "DEP", "EXERCISE", "FLUOR", "THAL", "OUTPUT")
numericInput("n",
 "How many records?", 5)
renderTable({
 head(heart, input$n)
}) 
列名 类型 定义
AGE 数值 岁为单位
SEX 二元类型 性别
CHESTPAIN 类别 4种胸痛级别
RESTBP 数值 静息血压
CHOL 数值 血清胆固醇
SUGAR 二元类型 空腹血糖浓度
ECG 类别 心电图结果,三个类别
MAXHR 数值 每分钟最大心率
ANGINA 二元类别 心绞痛来自运动?
DEP 数值 ST诱发比例,不懂
EXERCISE 有序类别 运动ST曲线峰值斜率
FLUOR 数值 透视着色的主要血管数量
THAL 类别 地中海贫血,三类
OUTPUT 二元类别 是否有心脏病(输出)

CHESTPAIN、THAL和ECG特征都是分类特征。EXERCISE变量虽然是一个有序的分类变量,但它还是属于分类变量,所以它也必须编码为一个因子。
具体而言,R语言会把k个因子水平中的1个作为参考水平,然后根据其他因子水平创建k-1个二元特征,类似OneHot。 最后把结果变换到0和1即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#因子变换
heart$CHESTPAIN = factor(heart$CHESTPAIN)
heart$ECG = factor(heart$ECG)
heart$THAL = factor(heart$THAL)
heart$EXERCISE = factor(heart$EXERCISE)
heart$OUTPUT = heart$OUTPUT-1
library(caret)
#拆分训练集和测试集
heart_sampling_vector <- createDataPartition(heart$OUTPUT, p = 0.85, list = FALSE)
heart_train <- heart[heart_sampling_vector,]
heart_train_labels <- heart$OUTPUT[heart_sampling_vector]
heart_test <- heart[-heart_sampling_vector,]
heart_test_labels <- heart$OUTPUT[-heart_sampling_vector]
#glm-->广义线性模型
heart_model = glm(OUTPUT~.,data=heart_train,family=binomial("logit"))
summary(heart_model)
heart_model2 = glm(OUTPUT~AGE,data=heart_train,family=binomial("logit"))
summary(heart_model2)

准确性评估

1
2
3
4
5
6
train_predictions = predict(heart_model, newdata=heart_train, type="response")
train_class_predictions = as.numeric(train_predictions>0.5)
mean(train_class_predictions==heart_train$OUTPUT)
test_predictions = predict(heart_model, newdata=heart_test, type="response")
test_class_predictions = as.numeric(test_predictions>0.5)
mean(test_class_predictions==heart_test$OUTPUT)

注:利用了glmnet包来通过岭回归和lasso进行正则化可以进一步提高正确率到90% 混淆矩阵

1
2
3
4
(confusion_matrix <- table(predicted = train_class_predictions, actual = heart_train$OUTPUT))
(precision <- confusion_matrix[2,2]/sum(confusion_matrix[2,]))
(recall <- confusion_matrix[2,2]/sum(confusion_matrix[,2]))
(f = 2 * precision * recall / (precision + recall))

用bagging预测心脏病

用有放回的方式抽取样本,并用它们来训练我们的模型,假设我们训练11个模型。

M <- 11
n <- nrow(heart_train)
sample_vectors<-sapply(seq(1,11),function(x) { return(sample(n,n,replace=T)) })

train_1glm <- function(sample_indices) { 
	data <- heart_train[sample_indices,]; 
	model <- glm(OUTPUT~., data=data, family=binomial("logit")); 
	return(model)
}

它现在包含了11个训练好的模型。 作为评估我们模型的第一个方法,我们要利用每个模型训练用到的数据。下面构造bag变量,代表每个模型具体用到哪些不重复的数据。

1
2
3
4
5
6
7
8
9
models <- apply(sample_vectors, 2, train_1glm)

get_1bag <- function(sample_indices) {
	unique_sample <- unique(sample_indices); 
	df <- heart_train[unique_sample, ]; 
	df$ID <- unique_sample; 
	return(df)
}
bags<-apply(sample_vectors, 2, get_1bag)

现在我们有了一个模型的列表和一个训练时用到的数据框的列表,且后者中没有重复的观测数据。根据这两个列表,我们可以得出一批预测结果。对于每个训练数据框,我们会附加一个叫PREDICTIONS{m}的新列,其中的{m}会是作出这个预测所用到的模型编号。这样,bags列表里的第一个数据框就会有一个名为PREDICTIONS 1的预测结果列,第二个数据框会有一个名为PREDICTIONS 2的预测结果列,第三个的名为PREDICTIONS 3,以此类推。

1
2
3
4
5
6
7
glm_predictions <- function(model, data, model_index) {
	colname <- paste("PREDICTIONS",model_index);
	data[colname] <- as.numeric(predict(model,data,type="response") > 0.5); 
	return(data[,c("ID",colname), drop=FALSE])
}
training_predictions <- mapply(glm_predictions,models,bags,1:M,SIMPLIFY=F)
#head(training_predictions[1],10)

下一步,我们要把所有这些数据框合并为一个数据框,其中的行是原始数据框的行(因而对应了数据集里的观测数据),列则是每个模型对该观测数据的预测结果。对于在抽样过程中没有被选取来训练特定模型的特定行(即观测数据),它在对应此模型作出预测结果的列里会有一个NA值。

1
2
3
4
train_pred_df <- Reduce(function(x, y) merge(x, y, by = "ID", all = T), training_predictions)
renderTable({
 head(train_pred_df, 10)
}) 

根据这个数据框每一行的多数票,我们现在就可以产生整个装袋模型对训练数据的预测结果。一旦我们得到了这些结果,就只需要用预测结果和原始的heart_train数据框里对应的标签值进行匹配,并计算出我们的精确度:

1
2
train_pred_vote<-apply(train_pred_df[,-1],1,function(x) as.numeric(mean(x,na.rm=TRUE)>0.5))
(training_accuracy<-mean(train_pred_vote==heart_train$OUTPUT[as.numeric(train_pred_df$ID)]))

现在我们有了装袋模型的第一个精确度衡量指标。这和衡量在训练数据上的精确度是类似的。现在对每个模型采用袋外观测数据重复这个过程,计算出袋外精确度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# Out of Bag Accuracy

get_1oo_bag <- function(sample_indices) {
	unique_sample <- setdiff(1:n,unique(sample_indices)); 
	df <- heart_train[unique_sample,]; 
	df$ID <- unique_sample; 
	if (length(unique(heart_train[sample_indices,]$ECG)) < 3) df[df$ECG == 1,"ECG"] = NA; 
	return(df)
}
oo_bags <- apply(sample_vectors,2, get_1oo_bag)
oob_predictions<-mapply(glm_predictions, models, oo_bags, 1:M,SIMPLIFY=F)
oob_pred_df <- Reduce(function(x, y) merge(x, y, by="ID", all=T), oob_predictions)
oob_pred_vote<-apply(oob_pred_df[,-1], 1, function(x) as.numeric(mean(x,na.rm=TRUE)>0.5))
(oob_accuracy<-mean(oob_pred_vote==heart_train$OUTPUT[as.numeric(oob_pred_df$ID)],na.rm=TRUE))

最后,利用heart_test数据框第三次重复这个过程,获得测试集的精确度:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Test Accuracy

get_1test_bag <- function(sample_indices) {
     df <- heart_test; 
     df$ID <- row.names(df); 
     if (length(unique(heart_train[sample_indices,]$ECG)) < 3) 
         df[df$ECG == 1,"ECG"] = NA; 
     return(df)
 }

test_bags <- apply(sample_vectors,2, get_1test_bag)
test_predictions<-mapply(glm_predictions, models, test_bags, 1 : M, SIMPLIFY = F)
test_pred_df <- Reduce(function(x, y) merge(x, y, by="ID",all=T), test_predictions)
test_pred_vote<-apply(test_pred_df[,-1],1,function(x) as.numeric(mean(x,na.rm=TRUE)>0.5))
(test_accuracy<-mean(test_pred_vote==heart_test[test_pred_df$ID,"OUTPUT"],na.rm=TRUE))

使用不同的M值把前面的代码运行多次,并把结果保存在一个数据框里

注意:只有遇到非线性模型时,采用装袋才有意义。因为装袋过程就是对产生的模型进行一次取平均值(线性运算)的处理,从而我们在线性回归里就不会看到任何改善,因为我们没有增加模型的表达力。

bagging的局限性

  1. bagging涉及对多个模型作出的预测结果取平均值,当目标函数比较平滑时,这样会减少偏误。但如果目标函数并不平滑,使用装袋实际上会引入更多偏误。
  2. 当某些类别特征包含了罕见的值时,它们可能在某些自助样本里根本不出现。如果发生了这种情况,用这些样本构建的模型在测试集里遇到这个新特征的时候就无法作出预测了。
  3. 构建的不同模型严格地说并不是真的相互独立的,因为它们使用的还是同一套输入特征。如果模型是相互独立的,取平均的过程会更有效果。
  4. 弱化模型的解释力。

增强

增强(Boosting)对如何把模型组合起来达到更大性能的问题给出了另一个思路。它尤其适用于弱学习器(weak learner)。弱学习器指能产生比随机猜测模型的精确度略好一点的模型。创建弱学习器的一种方式是采用一个复杂度可配置的模型。

预测大气中伽马射线的辐射

伽马射线会留下独特的椭圆形态,所以我们可以创建一组特征来描述这些模式。我们要用到的数据集是MAGIC Gamma Telescope data set,它存放在UCI机器学习知识库(UCI Machine Learning repository)里。数据包含了19020条观测数据.

首先,我们要把数据加载到一个名为magic的数据框里,并把CLASS输出变量,重新编码,对伽马射线和背景辐射分别采用1和-1类别:

1
2
3
4
5
6
magic <- read.csv("./data/magic04.dat", header=FALSE)
names(magic)=c("FLENGTH", "FWIDTH", "FSIZE", "FCONC", "FCONC1", "FASYM", "FM3LONG", "FM3TRANS", "FALPHA", "FDIST", "CLASS")
magic$CLASS = as.factor(ifelse(magic$CLASS=='g',1,-1))
renderTable({
 head(magic, 10)
}) 

采用典型的80%/20%拆分比例,把这个数据框拆分为一个训练数据框和一个测试数据框:

1
2
3
4
5
6
# 划分训练集测试集
magic_sampling_vector <- createDataPartition(magic$CLASS, p = 0.80, list = FALSE)
magic_train <- magic[magic_sampling_vector,1:10]
magic_train_output <- magic[magic_sampling_vector,11]
magic_test <- magic[-magic_sampling_vector,1:10]
magic_test_output <- magic[-magic_sampling_vector,11]

对于神经网络,我们在对输入进行归一化的情况下往往会得到良好的精确度,因此,在训练任何模型之前,我们会进行这个预处理步骤

1
2
3
4
5
# 数据预处理
magic_pp <- preProcess(magic_train, method = c("center", "scale"))
magic_train_pp <- predict(magic_pp, magic_train)
magic_train_df_pp <- cbind(magic_train_pp, CLASS = magic_train_output)
magic_test_pp <- predict(magic_pp, magic_test)

增强是设计用来和弱学习器搭配发挥最佳效果的,因此,我们会在隐藏层里使用很小数量的隐藏神经元。具体而言,我们会从能做出的最简单的只采用单个隐藏神经元的多层感知器入手。为了理解使用增强的效果,我们会训练单个神经网络并衡量其精确度,从而得到基准性能。做法如下:

1
2
3
4
library(nnet)
n_model <- nnet(CLASS~., data = magic_train_df_pp, size = 1)
n_test_predictions = predict(n_model, magic_test_pp, type = "class")
(n_test_accuracy <- mean(n_test_predictions == magic_test_output))

精度不错,看能不能通过adaboost提高。

编写自己的函数AdaBoostNN(),它接收的输入是一个数据框、输出变量名、我们要创建的单隐藏层神经网络模型的个数,以及这些神经网络里含有的隐藏神经元个数。然后,这个函数会返回一个模型及其权值的列表. 我们还需要一个作出预测的函数。这个函数会接受训练函数AdaBoostNN()产生的输出和一个测试数据集作为输入。我们把这个函数命名为AdaBoostNN.predict()。 在这个函数里,我们首先会从前面函数产生的列表中提取模型及模型权值。然后我们会创建一个预测结果的矩阵,其中的每列对应某个模型产生的预测结果向量。这样,在这个矩阵里的列数就会和用于增强的模型数量相等。 然后,把每个模型产生的预测结果乘以它们对应的模型权值。例如,第一个模型的每个预测结果都在预测矩阵第一列里,它们的值会乘以第一个模型的权值α1。最终,在最后一个步骤中,通过对每条观测数据求其加权预测结果之和并取结果的正负号,把加权观测数据的矩阵缩减为一条观测数据的向量。然后这个预测结果的向量就会被我们的函数返回。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
AdaBoostNN <- function(training_data, output_column, M, hidden_units) {
  require("nnet")
  models<-list()
  alphas<-list()
  n <- nrow(training_data)
  model_formula <- as.formula(paste(output_column,'~.',sep=''))
  w <- rep((1/n),n)
  for (m in 1:M) {
    model<-nnet(model_formula,data=training_data,size=hidden_units, weights=w)
    models[[m]]<-model
    predictions<-as.numeric(predict(model,training_data[,-which(names(training_data) == output_column)],type = "class"))
    errors<-predictions!=training_data[,output_column]
    error_rate<-sum(w*as.numeric(errors))/sum(w)
    alpha<-0.5*log((1-error_rate)/error_rate)
    alphas[[m]]<-alpha
    temp_w<-mapply(function(x,y) if (y) {x*exp(alpha)} else {x*exp(-alpha)},w,errors)
    w<-temp_w/sum(temp_w)
  }
  return(list(models=models, alphas=unlist(alphas)))
}

AdaBoostNN.predict<-function(ada_model,test_data) {
  models<-ada_model$models
  alphas<-ada_model$alphas
  prediction_matrix<-sapply(models,function (x) as.numeric(predict(x,test_data,type = "class")))
  weighted_predictions<-t(apply(prediction_matrix,1,function (x) mapply(function(y,z) y*z,x,alphas)))
  final_predictions<-apply(weighted_predictions,1,function(x) sign(sum(x)))
  return(final_predictions)
}

训练10个带有单个隐藏神经元的神经网络模型,观察增强是否会提高精确度

1
2
3
ada_model <- AdaBoostNN(magic_train_df_pp, 'CLASS', 10, 1)
predictions <- AdaBoostNN.predict(ada_model, magic_test_pp)
mean(predictions == magic_test_output)

增强10个模型看起来会给我们在精确度上带来微小的改进,但也许训练更多模型会带来更多不同。我们还对弱学习器之间复杂度的关系感兴趣,它是用隐藏神经元的数量以及在这个数据集上进行增强后能够预期的性能提高来衡量的。

使用不同参数来调整神经网络:

1
magic_boost <- expand.grid(M = c(5, 10, 20, 50, 100, 200), hidden_units = c(1, 3))

tp

对于只有1个隐藏神经元的神经网络,随着增强模型的数量增加,我们看到了精确度的提高,但在100个模型之后,精确度会逐渐减小,到200个模型时实际上精确度反而更小了一点。对于这些网络,性能比单个模型的基准有明显提高。当通过采用带有3个隐藏神经元的隐藏层来增加分类器的复杂度时,我们在性能上得到的提高就小得多。在模型数量为200个时,两种集成方法的性能都在相似水平上,这表明在这个点上,我们的精确度会受限于我们所训练的模型的类型。