集成方法： 让多个训练好的模型作为一个整体协同工作，从而产生比其中任何单个模型都更为强大的模型。

同质个体学习器按照个体学习器之间是否存在依赖关系可以分为两类：
1）个体学习器之间存在强依赖关系，一系列个体学习器基本都需要串行生成，即：除了训练第一个之外，其他的学习器学习都需要依赖于前面生成的学习的结果。代表算法是boosting系列算法
2）个体学习器之间不存在强依赖关系，一系列个体学习器可以并行生成，代表算法是bagging和随机森林（Random Forest）系列算法。

装袋（bagging）

• 自助聚合的简写（bootstrap aggregation）
• 使用同一个数据集的不同样本来训练一个模型的多个版本。然后这些模型分类就投票，回归就取平均值

增强（boosting）

• 训练一系列模型，并给没有正确分类或远离其他预测值的观测数据分配权重，以便后续训练的模型把它们放到优先地位

1. 装袋（bagging）

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

1) 装袋（bagging）的方法
步骤1： 创建一个大小为𝑛的有放回的随机样本(bootstrapping)
意味着有些数据是重复的，有些根本没用过
步骤2： 利用这些抽样的数据集训练一个分类模型
步骤3： 对抽样数据集的每条观测数据，记录模型给它分配的类别
步骤4： 重复这个抽样过程M次，训练M个模型
步骤5： 对于原始训练数据集里的每条观测数据，通过不同模型的多数票计算预测出的类别
如果M=61，某条数据出现50次，37次是1，13次是-1，整体预测就是1
步骤6： 利用训练集提供的标签来计算模型的精确度

2) “每次进行有放回的抽样时，能得到多少不同的观测数据？”[𝑛条数据抽𝑛次]
每个做出抽样里最后有63%是不同的观测数据。
在抽样过程中，没有抽到某条观测数据𝑥1的概率就是伯努利实验失败𝑛次的结果

3)边缘和袋外观测数据

• 边缘：分类比例之间的差异
  比如，某特定的观测数据𝑥1 ，有85%预测正确，15%错误；某条观测数据𝑥2 ，有53%预测正确，47%错误。我们称之为对前一条观测数据具有更强的正确分类能力，我们都喜欢较大边缘(差异)的分类器。
  尽可能对所有观测数据有较大边缘的分类器

注意：
1.在给每个模型产生预测值集合时，使用的是和训练模型所采用的同一批数据。仔细思考步骤3，会发现分类用的是在步骤2训练模型使用的同一批抽样数据
2.即使最终依赖于一个取平均的过程来获得装袋分类器对未知数据的估计精确度，但我们没使用任何未知数据

• 袋外观测数据
  步骤1中，某次迭代中没有选中的观测数据称为袋外(out -of-bag, OOB)观测数据。
  实际上可以，用OOB观测数据来记录某个模型的精确度，最后对所有OOB精确率取平均值。

• 自助样本（bootstrapped sample）
  原始数据集有放回的抽样产生的样本
  和同一个分布中抽取多个样本相似

• 装袋可以降低方差
  例如，独立同分布的样本均值的方差是𝜎2/𝑛，观测数据量的增加，方差是递减的

• 包和函数
  ipred包，rpart()来构建的装袋预测工具
  bagging()函数，参数nbagg(默认25)，coob=T声明用袋外样本来计算精确度

用装袋预测复杂技能学习

# 获取数据
skillcraft <- read.csv("SkillCraft1_Dataset.csv")

#将一些特征因子化
skillcraft<-skillcraft[-1] #去掉ID列
skillcraft$TotalHours <- factor(skillcraft$TotalHours)
skillcraft$HoursPerWeek <- factor(skillcraft$HoursPerWeek)
skillcraft$Age <- factor(skillcraft$Age)
skillcraft$TotalHours=as.numeric(levels(skillcraft$TotalHours))[skillcraft$TotalHours]
skillcraft$HoursPerWeek=as.numeric(levels(skillcraft$HoursPerWeek))[skillcraft$HoursPerWeek]
skillcraft$Age=as.numeric(levels(skillcraft$Age))[skillcraft$Age]
#过滤所有列中有NA的行
skillcraft<-skillcraft[complete.cases(skillcraft),]

#划分训练集和测试集
library(caret)
set.seed(133)
skillcraft_sampling_vector <- createDataPartition(skillcraft$LeagueIndex, p = 0.80, list = FALSE)
skillcraft_train <- skillcraft[skillcraft_sampling_vector,]
skillcraft_test <- skillcraft[-skillcraft_sampling_vector,]

#定义函数计算均方差SSE
compute_SSE <- function(correct,predictions) {
  return(sum((correct-predictions)^2))
}

# 装袋
library("ipred")
baggedtree <- bagging(LeagueIndex~., data=skillcraft_train, nbagg=100, coob=T)
baggedtree_predictions <- predict(baggedtree, skillcraft_test)
(baggedtree_SSE <- compute_SSE(baggedtree_predictions, skillcraft_test$LeagueIndex))

用装袋预测心脏病

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

#读取数据
heart <- read.table("heart.dat", quote="\"")
names(heart) <- c("AGE", "SEX", "CHESTPAIN", "RESTBP", "CHOL", "SUGAR", "ECG", "MAXHR", "ANGINA", "DEP", "EXERCISE", "FLUOR", "THAL", "OUTPUT")

#将一些特征因子化
heart$CHESTPAIN = factor(heart$CHESTPAIN)
heart$ECG = factor(heart$ECG)
heart$THAL = factor(heart$THAL)
heart$EXERCISE = factor(heart$EXERCISE)

#将输出范围限制在0和1之间（之前是1和2）
heart$OUTPUT = heart$OUTPUT-1

#划分测试集和训练集
library(caret)
set.seed(987954)
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]

#准备训练11个小型的训练集
M <- 11
seeds <- 70000 : (70000 + M - 1)
n <- nrow(heart_train) #心脏病数据的行数
#每次抽n行，有放回的抽样。抽M次。（抽出来的是行数）
sample_vectors<-sapply(seeds,function(x) { set.seed(x); return(sample(n,n,replace=T)) })
#根据序号，找到对应的数据，放到glm函数中进行逻辑回归
train_1glm <- function(sample_indices) { 
    data <- heart_train[sample_indices,]; 
    model <- glm(OUTPUT~., data=data, family=binomial("logit")); 
    return(model)
}
#生成11个逻辑回归的结果
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) 

#对每个数据，计算11个模型的预测结果
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])
}
#mapply跟sapply类似，只是是一个多变量版本。simplify=True表示返回一个矩阵，否则返回一个list
#由于glm_predictions需要传递多个参数，因此需要用mapply。
training_predictions <- mapply(glm_predictions,models,bags,1:M,SIMPLIFY=F)
training_predictions
#Reduce函数是将每次计算后的结果保留，并与下一个数字进行计算，这是和 apply 函数不同的
#这里是利用reduce函数将多个数据框按照同一列merge
train_pred_df <- Reduce(function(x, y) merge(x, y, by = "ID", all = T), training_predictions)
#进行结果投票，对于每一行，去掉ID列后，求这一行的均值。若均值大于0.5，则返回1，否则返回0
train_pred_vote<-apply(train_pred_df[,-1],1,function(x) as.numeric(mean(x,na.rm=TRUE)>0.5))
#计算精度（准确度），准确度从原来的0.86提高到0.88
(training_accuracy<-mean(train_pred_vote==heart_train$OUTPUT[as.numeric(train_pred_df$ID)]))

## 用袋外数据来进行预测（类似于测试集）
get_1oo_bag <- function(sample_indices) {
  #注意，这里的n是个全局变量，上面定义过，是心脏病数据的行数。
  #使用setdiff，将统计从1到n的数中与sample_indices不同的数，得到袋外数据的行号
    unique_sample <- setdiff(1:n,unique(sample_indices)); 
    df <- heart_train[unique_sample,];  #根据行号取出这些数据
    df$ID <- unique_sample; #赋予一个ID
    #如果ECG的数量小于3，则将ECG=1的数据的全部变为NA。
  #因为ECG=1的数量很少，很有可能训练集的时候没有。导致模型没有见过ECG的特征，所以将ECG判定为NA，忽略此特征
    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)
#对每个数据，计算11个模型的预测结果
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))
#查看投票后的精度（袋外数据的精度为0.8）
(oob_accuracy<-mean(oob_pred_vote==heart_train$OUTPUT[as.numeric(oob_pred_df$ID)],na.rm=TRUE))

#再来使用测试集进行验证
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))
#可以看到，测试集上的精度为0.925
(test_accuracy<-mean(test_pred_vote==heart_test[test_pred_df$ID,"OUTPUT"],na.rm=TRUE))

• 装袋的局限性
  1.装袋会引入偏误
  装袋本身取平均的过程会使得整体输出平滑化
  如果目标函数并不光滑，装袋会引入偏误
  2.高杠杆率点
  该数据对模型的输出函数的影响力比其他点高很多
  一个抽取的自助样本没有包括该高杠杆率点，模型会很不一样
  3.根源：模型本身并不相互独立
  4.装袋对低偏误和高方差的模型有效果

2. 增强（boosting）

增强(Boosting)提供了另一个思路，适用于弱学习器(weak learner)（能产生比随机猜测模型的精度略好一点的模型）。

比如，只有很少隐藏层神经元的多层感知器网络；树桩(stump)只有一个节点的决策树，只有一个分支。

• 增强（Boosting）
  通过在训练数据上构建一个模型，然后衡量在那些训练数据上的分类精确度。被模型错误分类的每条观测数据会得到一个比正确分类的观测数据更大的权值，然后这个模型会用这些新的权值重新训练。这个过程会重复很多次，每次根据迭代是否正确分类进行调整
  回归问题用预测值和标签值的距离衡量来调整权值

• 装袋和增强的区别
  1.装袋是取一组训练数据的随机自助样本，然后用这些不同的样本训练一个模型的多个版本
  2.增强没有随机成分，所有模型都会用所有的训练数据
  

• 自适应增强(adaptive boosting)
  用于二元分类的自适应增强
  输入：
  data：包含了输入特征和一个二元输出标签列的输入数据框
  M：表示我们需要训练的模型数量的一个整数
  输出：
  model：由M个训练好的模型组成的序列
  alpha：包含M个模型权值的向量
  方法：
  1.初始化一条观测数据权值的向量𝑤，它的长度为𝑛，初始元素为𝑤𝑖 = 1/𝑛。该向量会在每次迭代中更新
  2.使用观测数据权值的当前值和训练集里的所有数据训练成一个分类器模型𝐺𝑚
  3.计算权值的误差率，用所有错误分类观测数据乘以其权值然后求和，再除以权值向量的综合

  
  4.设置该模型的模型权值α𝑚为精确度和误差率之比的对数
  
  5.为下次迭代更新观测数据权值向量𝑤。错误分类的观测数据的权值会乘以𝑒^𝑎𝑚，从而增大它们在下一次迭代时的权值。正确分类的观测数据的权值会乘以𝑒^−𝑎𝑚，从而减少它们在下一次迭代时的权值
  6.对权值向量重新进行归一化，使权值的总和等于1
  7.把2-6步骤重复𝑀次，产生𝑀个模型
  8.定义集成分类器为所有增强模型输出的加权总和

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

望远镜照相机拍下的辐射中出现的模式，以此预测某个模式是来源于泄露到大气中的伽马射线，还是来自常规的背景辐射
伽马射线包含19020条数据，输出是2元的，伽马射线（g），背景强子辐射（b）

magic <- read.csv("magic04.data", header=FALSE)
names(magic)=c("FLENGTH", "FWIDTH", "FSIZE", "FCONC", "FCONC1", "FASYM", "FM3LONG", "FM3TRANS", "FALPHA", "FDIST", "CLASS")
#把输出变量进行因子化，g转化为1，b变为-1
magic$CLASS = as.factor(ifelse(magic$CLASS=='g',1,-1))

#划分训练集和测试集
set.seed(33711209)
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]

#对数据进行归一化、中心化处理
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个神经元的神经网络进行拟合
library(nnet)
n_model <- nnet(CLASS~., data = magic_train_df_pp, size = 1)
#在测试集上进行预测
n_test_predictions = predict(n_model, magic_test_pp, type = "class")
#计算测试集上的精度（精度为78.67%）
(n_test_accuracy <- mean(n_test_predictions == magic_test_output))

#构造AdaBoost
AdaBoostNN <- function(training_data, output_column, M, hidden_units) {
  #library和require都可以载入包，但二者存在区别。
  #在一个函数中，如果一个包不存在，执行到library将会停止执行，require则会继续执行。
  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) {
    #构造hidden_units个隐藏单元的神经网络
    model<-nnet(model_formula,data=training_data,size=hidden_units, weights=w)
    models[[m]]<-model #将构造后的模型存入模型列表中
    #先将training_data中，输出变量去除，然后使用模型进行预测，预测后将结果转为数值型
    predictions<-as.numeric(predict(model,training_data[,-which(names(training_data) == output_column)],type = "class"))
    errors<-predictions!=training_data[,output_column] #如果分类正确，则返回false，分类错误返回true
#as.numeric(errors)，会是的分类正确变为0，分类错误变为1。然后再乘以权重，相当于是把所有分类错误的权重求和，再除以全部的权重之和，得到错误率。
    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)
}

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

# 探索不同隐藏单元和模型数
magic_boost <- expand.grid(M = c(5, 10, 20, 50, 100, 200), hidden_units = c(1, 3))
magic_boost$accuracy <- mapply(function(x,y) { ada_model <- AdaBoostNN(magic_train_df_pp, 'CLASS', x, y); predictions <- AdaBoostNN.predict(ada_model, magic_test_pp); mean(predictions == magic_test_output) }, magic_boost$M, magic_boost$hidden_units )
magic_boost$hidden_units <- factor(magic_boost$hidden_units)

library(ggplot2)
p <- ggplot(data = magic_boost, aes(x = M, y = accuracy, shape = hidden_units, linetype = hidden_units))
p <- p + geom_line()
p <- p + geom_point(size=3)
p <- p + ggtitle("Accuracy of Boosted Neural Networks on the Magic Telescope Data")
p <- p + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2), legend.position="bottom")
p <- p + xlab("Number of Models")  
p <- p + ylab("Accuracy")
p <- p + scale_linetype_manual(name = "Hidden Units", values = c(1,2))
p <- p + scale_shape_manual(name = "Hidden Units", values = c(15,16))
p

利用增强算法预测复杂技能学习

随机梯度增强的增强技术

• 在增强的每次迭代中，计算在当前迭代中训练的模型产生的误差所在方向的梯度
• 使用gbm()函数

• 构建树的数量n.trees；算法的学习率shrinkage

  
  

3.随机森林

随机森林=集成方法+决策树

在skillcraft数据集中，我们预期会看到APM在大部分装袋树顶端被选中作为特征

• 会阻碍我们从装袋中获得方差缩减的优势
• 换言之，我们构建的树不够相互独立
• 具有很多共同的特征和分裂点，最后取均值过程在减少集成方差方面也不会那么有效

随机森林：在树的构建过程中引入随机元素
在输入特征的总数里面随机抽取一个𝑚𝑡𝑟𝑦的随机样本，例如：所有特征总数的平方根或者三分之一
这样的抽样步骤可以强制让装袋的树结构互不同。当输入特征具有相关性时，对输入特征进行抽样也是有用的

• 在一个相关的特征集合中只选取其中的一个特征，而忽略其他特征
• 我们在对输入特征进行抽样时，得到相互竞争的相关特征的可能性会更小
• 树木的数量是决策变量

在输入特征的总数中抽取一些随机样本，然后构建不同的树，避免树的相似性。抽样的随机天性是为了避免过拟合。尤其特征数量超过观测数据数量，随机森林是个好的选择。randomForest()函数

library("randomForest")
library("e1071")
#在输入特征的总数里面随机抽取一个mtry的随机样本
rf_ranges <- list(ntree = c(500, 1000, 1500, 2000), mtry = 3:8)
#使用tune进行自动调参
rf_tune <- tune(randomForest, LeagueIndex~., data = skillcraft_train, ranges = rf_ranges)

#查看调参后的最优参数组合
rf_tune$best.parameters
rf_best <- rf_tune$best.model

rf_best_predictions = predict(rf_best,skillcraft_test)
(rf_best_SSE <- compute_SSE(rf_best_predictions, skillcraft_test$LeagueIndex))

par(mai=c(1.02,2.5,0.82,0.42))
rf_imp <- importance(rf_best)
rf_imp <- as.data.frame(rf_imp)[order(-rf_imp),,drop=F]
n<-rownames(rf_imp)
barplot(t(rf_imp),horiz=T,las=1,main="Variable Importance in Random Forest", names=n, col = "gray")

4.用R来比较几种算法的优劣（决策树、bagging、随机森林、回归决策树）：

#计算SSE
compute_SSE <- function(correct,predictions) {
  return(sum((correct-predictions)^2))
}

skillcraft <- read.csv("SkillCraft1_Dataset.csv")
#数据预处理
skillcraft<-skillcraft[-1]
skillcraft$TotalHours <- factor(skillcraft$TotalHours)
skillcraft$HoursPerWeek <- factor(skillcraft$HoursPerWeek)
skillcraft$Age <- factor(skillcraft$Age)
skillcraft$TotalHours=as.numeric(levels(skillcraft$TotalHours))[skillcraft$TotalHours]
skillcraft$HoursPerWeek=as.numeric(levels(skillcraft$HoursPerWeek))[skillcraft$HoursPerWeek]
skillcraft$Age=as.numeric(levels(skillcraft$Age))[skillcraft$Age]
skillcraft<-skillcraft[complete.cases(skillcraft),] #去除空行

#划分测试集和训练集
library(caret)
set.seed(133)
skillcraft_sampling_vector <- createDataPartition(skillcraft$LeagueIndex, p = 0.80, list = FALSE)
skillcraft_train <- skillcraft[skillcraft_sampling_vector,]
skillcraft_test <- skillcraft[-skillcraft_sampling_vector,]

#使用rpart进行决策树分类
library(rpart)
regtree <- rpart(LeagueIndex~., data=skillcraft_train)
#在测试集上进行预测
regtree_predictions = predict(regtree, skillcraft_test)
(regtree_SSE <- compute_SSE(regtree_predictions, skillcraft_test$LeagueIndex))#评判测试集上的预测效果

#baggin函数默认使用rpart包的决策树方法来进行分类
#coob=T表示用袋外样本来评估误差
#nbagg表示用100个袋子，也就是生成100个小的模型
library("ipred")
baggedtree <- bagging(LeagueIndex~., data=skillcraft_train, nbagg=100, coob=T)
baggedtree_predictions <- predict(baggedtree, skillcraft_test)
(baggedtree_SSE <- compute_SSE(baggedtree_predictions, skillcraft_test$LeagueIndex))

#使用随机森林来进行预测
library("randomForest")
rf<-randomForest(LeagueIndex~., data=skillcraft_train)
rf_predictions = predict(rf,skillcraft_test)
(rf_SSE <- compute_SSE(rf_predictions, skillcraft_test$LeagueIndex))

#使用回归决策树来进行预测
library("gbm")
boostedtree <- gbm(LeagueIndex~., data=skillcraft_train, distribution="gaussian", n.trees=10000, shrinkage=0.1)
best.iter <- gbm.perf(boostedtree,method="OOB")
boostedtree_predictions = predict(boostedtree, skillcraft_test,best.iter)
(boostedtree_SSE <- compute_SSE(boostedtree_predictions, skillcraft_test$LeagueIndex))

可以看到，使用决策树，最后的SSE为799；使用bagging，SSE为677；使用随机森林，SSE为575；使用GBDT回归决策树，SSE为549。

5.小结

装袋

• 训练数据的自助样本构建同一模型的多个版本
• 用这些模型做出预测结果取平均值来构建
• 消减由于过拟合产生的误差，减少了方差

增强

• 使用所有数据来构建模型，训练一系列模型，不断修正参数
• 自适应增强和随机梯度增强