集成学习算法应用比较广泛，只要包含Bagging，Boosting与Stacking三大类方法，这是本文的主要内容：
  1. 集成算法的主要思想；
  2. 集成算法的Bagging，Boosting与Stacking方法的数学模型与推导；
  3. 集成算法的实现与应用；
  4. XGBoost算法介绍；

---

  到目前为止机器学习入门告一个结束；后面继续写框架扩展部分并深入补充部分优化算法；并使用C++实现一个系列，包含OpenCV与Tensorflow；

集成学习算法说明

集成学习介绍

• 集成学习(ensemble learning)本身不是一个单独的机器学习算法，而是通过构建并结合多个机器学习器来完成学习任务。

• 集成学习可以用于：

  • 分类问题集成
  • 回归问题集成
  • 特征选取集成
  • 异常点检测集成

• 怎么通过多个机器学习的结果得到一个最终的学习结果，这是集成学习算法的关键，这需要组织不同机器学习的策略。一般从两个角度出发：

  • 数据集
  • 算法
• 集成学习思想示意图

• 名词说明：

  • 个体学习器，也称：基学习器
  • 在基学习器上形成的学习器称为元学习器。（下面会有详细的介绍）

• 集成学习的核心：

  • 个体机器学习器
  • 结合策略

其他集成学习的说明

• 集成学习算法的有效性：

  • 比如用线性函数去拟合曲线，效果当然不会很好；
  • 如果把数据集拆分开，使用多个线性回归集成，这些线性回归组合而成的折线其实也算集成学习。

• 通过集成学习，可以提高算法的泛化能力，所谓泛化能力就是通用性、广泛性、稳定性更强。

• 集成学习集成不同算法的条件：

  • 需要选择有差异性的机器学习算法； 比如没有差异的分类器，集成起来的分类效果也不会有什么变化。
  • 集成的学习算法的准确率需要50%以上；比如分类器的准确率小于0.5，随着集成规模的增加，分类准确率会不断下降；如果准确度大于0.5，那么最终分类准曲率可以趋近于1。

个体学习器构建

个体学习器构建方式

1. 算法不同：比如对相同数据集，使用贝叶斯分类，SVM分类算法，训练出多个基学习器。
2. 算法相同：比如对不同的数据集，使用决策树训练出多个决策树基学习器。
   • 训练集不同；
   • 参数不同

个体学习构建算法

• 常见的基学习器的构建算法有：
  1. boosting：用于减少偏差的boosting
     • AdaBoost算法
     • GradientBoosting算法
     • HistGradientBoosting算法
     • XGBoost
  2. Bagging：用于减少方差
     • 随机森林
  3. Stacking：用于提升预测结果

个体学习器结合策略

• 把多个基学习器（弱学习器）的按照不同的结合策略，形成强学习器，常见结合策略有两种：
  • Averaging 平均法
    • 一般在回归问题中采用
  • Voting 投票法
    • 一般在分类问题中采用

Averaging 平均法结合策略

• 假设基学习器数量是个，每个基学习器的预测输出是:

  • 

• 最终平均值输出是：

  • 

• 如果对基学习器考虑加权,每个基学习器的加权为，则加权输出为：

  • 

Voting 投票法结合策略

• 假设我们的预测类别是，对于任意一个预测样本，我们的个弱学习器的预测结果分别是:

  • 

• 投票法是相对多数投票法，也就是少数服从多数：

  • 也就是个弱学习器的对样本的预测结果中，数量最多的类别为最终的分类类别。
  • 如果不止一个类别获得最高票，则随机选择一个做最终类别。

• 投票方式：
  • 硬投票（hard）：基于分类标签投票
  • 软投票（soft）：基于分类标签概率投票

sklearn集成学习应用体验

import sklearn.datasets as ds
from sklearn.model_selection import train_test_split
# HistGradientBoostingClassifier 分类器不存在
from sklearn.ensemble import AdaBoostClassifier, GradientBoostingClassifier
from sklearn.ensemble import BaggingClassifier, ExtraTreesClassifier,RandomForestClassifier, IsolationForest
from sklearn.ensemble import VotingClassifier

# 加载数据集
data,target = ds.load_iris(return_X_y=True)  

# 切分数据集
data_train, data_test, target_train, target_test = train_test_split(
    data,    # 数据集
    target,    # 数据集标签
    test_size=0.2)  

Boosting 算法集成学习算法

AdaBoostClassifier

classifier = AdaBoostClassifier(
    base_estimator=None, 
    n_estimators=100, 
    learning_rate=0.01, 
    algorithm='SAMME.R',
    random_state=None)

classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：29

GradientBoostingClassifier

classifier = GradientBoostingClassifier(
    loss='deviance',   # 损失度量方式：默认偏差
    learning_rate=0.1,    # 学习率
    n_estimators=100,    # 基学习器
    subsample=1.0, 
    criterion='friedman_mse',   
    min_samples_split=2,   # 决策树的特性
    min_samples_leaf=1, 
    min_weight_fraction_leaf=0.0, 
    max_depth=3, 
    min_impurity_decrease=0.0, 
    min_impurity_split=None, 
    init=None, 
    random_state=None, 
    max_features=None, 
    verbose=0, 
    max_leaf_nodes=None, 
    warm_start=False, 
    # presort='auto', 
    # validation_fraction=0.1, 
    # n_iter_no_change=None, 
    # tol=0.0001    # 结束误差
)
classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：29

Bagging 与随机森林集成学习算法

BaggingClassifier

classifier = BaggingClassifier(
    base_estimator=None, 
    n_estimators=50, 
    max_samples=1.0, 
    max_features=1.0, 
    bootstrap=True, 
    bootstrap_features=False, 
    oob_score=False, 
    warm_start=False, 
    n_jobs=-1,    # 使用None 报异常
    random_state=None, 
    verbose=0)
classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：29

RandomForestClassifier

classifier =RandomForestClassifier(
    n_estimators= 100,   # 森林数量 
    criterion='gini', 
    max_depth=None, 
    min_samples_split=2, 
    min_samples_leaf=1, 
    min_weight_fraction_leaf=0.0, 
    max_features='auto', 
    max_leaf_nodes=None, 
    min_impurity_decrease=0.0, 
    min_impurity_split=None, 
    bootstrap=True, 
    oob_score=False, 
    n_jobs= -1,   # 指定数量，与默认值的区别 
    random_state=None, 
    verbose=0, 
    warm_start=False, 
    class_weight=None)

classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：29

ExtraTreesClassifier 极限树

classifier = ExtraTreesClassifier(
    n_estimators=100, 
    criterion='gini', 
    max_depth=None, 
    min_samples_split=2, 
    min_samples_leaf=1, 
    min_weight_fraction_leaf=0.0, 
    max_features='auto', 
    max_leaf_nodes=None, 
    min_impurity_decrease=0.0, 
    min_impurity_split=None, 
    bootstrap=False, 
    oob_score=False, 
    n_jobs=-1, 
    random_state=None, 
    verbose=0, 
    warm_start=False, 
    class_weight=None)
classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：29

IsolationForest 离群点检测

outlier=  IsolationForest(
     n_estimators=100, 
     max_samples='auto', 
     contamination=0.1, 
     max_features=1.0, 
     bootstrap=False, 
     # n_jobs=-1, 
     # behaviour='old', 
     # random_state=None, 
     # verbose=0, 
     # warm_start=False
)

outlier.fit(data_train)
y_pred_train = outlier.predict(data_train)
y_pred_test = outlier.predict(data_test)    # 返回 1（正常inlier ），-1（异常outlier）
y_pred_test

array([-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])

其他集成学习算法

VotingClassifier投票分类器

# import warnings
# warnings.filterwarnings(module='sklearn*', action='ignore', category=DeprecationWarning)

from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
import sklearn.datasets as ds
from sklearn.model_selection import train_test_split
from sklearn.ensemble import VotingClassifier


# 加载数据集
data,target = ds.load_iris(return_X_y=True)  
# data = data[0:100]
# target = target[0:100]
# 切分数据集
data_train, data_test, target_train, target_test = train_test_split(
    data,    # 数据集
    target,    # 数据集标签
    test_size=0.2)  
c1 =RandomForestClassifier(n_estimators= 50,  criterion='gini')
c2 = LogisticRegression()
c3 = SVC(C=1.0, kernel='rbf')

classifier = VotingClassifier(
    estimators=[        # 参数是元组列表
        ('c1', c1),
        ('c2', c2),
        ('c3', c3)], 
    voting='hard', 
    weights=None, 
    n_jobs=-1, 
    flatten_transform=True)

classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)
correct_num = (pre == target_test).sum()
print(F'测试样本总数： {len(target_test)}，分类正确数：{correct_num}' )

测试样本总数： 30，分类正确数：30

sklearn与bumpy版本不兼容的问题

• 因为sklearn默认调用当前python环境的numpy模块，如果sklearn与numpy版本不匹配，会导致出现一些警告。

• 一般版本不兼容，都是因为numpy版本更新速度超过sklearn的更新速度导致，一般官网提示都是最低版本，没有对最高版本进行说明，所以容易出现问题。下面是官网对版本的要求说明：

  • 我当前sklearn版本：0.19.2，匹配的版本是1.13.3（当前最新版本是：1.16.3）
  • 版本匹配

• 更新sklearn最新版本：

  • pip install -U scikit-learn

• 提示

  • 上面代码可能报警告如下：

    • numpy的不兼容问题

  • 上面的警告问题是numpy自身的问题，改问题已经被修正，但没有体现在最新的版本中。

• sklearn最新版本下载安装：

• 关闭警告

import warnings
warnings.filterwarnings(module='sklearn*', action='ignore', category=DeprecationWarning)

• 降低numpy的版本也可以解决这个问题，可能不同的sklearn版本匹配不同的numpy版本。
  • pip install numpy==1.13.3

Bagging集成

Bagging算法思想

• Bagging是Bootstrap Aggregating（自助结合）的简写

• Bagging是通过组合随机生成的训练集而改进分类的集成算法。

  1. Bagging把训练数据集采用放回抽样，产生个训练集（是原始数据集的子集），使用自定的分类器算法对个子集进行训练而得到个分类器（基分类器）；
  2. 对待分类样本进行分类是，调用个基学习器，得到个分类结果，然后采用投票结合策略，得到的分类结果作为最终的分类结果。

• Baggind算法的特点：

  • 数据集不同（对已知数据集随机抽样产生多个不同的数据集）
  • 算法相同

• Bagging算法的关键点

  • Bagging算法的关键在抽样的方式，一般采用的抽样方法是：自助采样法（Bootstrap Sampling）

Bootstrap采样

• Bootstrap方法最初由美国斯坦福大学统计学教授Efron在1977年提出。作为一种崭新的增广样本统计方法，Bootstrap方法为解决小规模子样试验评估问题提供了很好的思路。

• Bootstrap采样的思想：

  • 假设原始样本数据集，样本容量为；
  • Bootstrap采样就是在原始数据的范围内作有放回的再抽样, 样本容量仍为，原始数据中每个样本每次被抽到的概率相等, 为 （这个在python编程中可以通过先获取之间的均匀随机数作为下标得到）, 所得样本称为Bootstrap样本。

• Bootstrap采样法解决了小样本时主成分分析的置信区间的计算方法：

  • 比如鸢尾花的四个特征，我们在做PCA主成分分析的时候，通过协方差矩阵得到特征值与特征向量，当样本数据集产生改变（数量与顺序）的时候，这个四个特征值也会产生改变，从而得到不同的结果，这个结果可能导致很大的学习误差。
  • 怎么评估这个特征值随着样本变化，有多大的可信度，Bootstrap方法就解决了这个问题：
    • 采用次boostrap抽样，可以得到个样本集（与原样本集容量一样：样本肯定有重复），从而得到个特征值
    • 采用特征值在某个分位区间出现的频次（记得matplotlib中密度图与直方图不？就是那个意思），可以得到特征值的置信区间的计算方式：值落在在某个范围内的概率。
    • 最终可以证明Bootstrap采样的有效性在于：只要采样的次数足够多，鸢尾花的主成分分析中最大特征值95%的概率落在一个比较窄的范围内，这种主成分就是在样本放生变化的时候就是可信的，特征的稳定性在未知的测试样本中也是可信的。
    • Bootstrap采样次数越多，其最后分类的结果置信度越高，方差越小。

• Bagging算法的每个基学习器之间没有联系，都是采用独立抽样的样本集训练出来的，我们称这种方式是并行方式。

• 随机森林是bagging的一个特例优化版本：

  • 算法都采用决策树。
  • 在样本随机采样基础上，又加上了特征的随机选择。

Bagging基础算法的实现

# 略

Boosting集成

Boosting算法思想

• Boosting（增长）算法是通过对数据集中样本采用不同的加权，形成不同数据集的方式改进分类效果：

  1. 首先从训练集用初始权重训练出一个基学习器1；
  2. 根据基学习器1的分类结果来更新训练样本的权重，分类错误的样本认为是分类困难样本，权重增加，反之权重降低；改变权重后得到新的样本；
  3. 对新的样本集，再训练一个新的基学习器2；
  4. 重复步骤2，得到个基学习器；
  5. 采用结合策略，得到最终的分类结果。

• Baggind算法的特点：

  • 数据集不同（对已知数据集采用不同加权产生多个不同的数据集）
  • 算法相同（可以采用不同的算法）

• Boosting算法的核心是权重的更新方式，因为权重的更新方式，产生不同的Boosting算法的

  • Ada Boosting算法
  • Gradient Boosting算法

AdaBoosting算法

AdaBoosting算法模型

分类与预测模型

• 假设，训练出个基学习器，则最终强学习器输出模型为：

  • 

• AdaBoosting算法就是要找出一组使得G(x)输出的损失最小。

损失函数模型

• AdaBoosting模型的损失最小的损失函数定义如下：

  • 

• 如果假设对所有样本：，损失函数最小模型为：

  • 
    • 其中：

• 最终完整表示损失函数为

  • 

损失函数模型简化

• 上面是复杂的全局优化问题，通常我们使用其简化版，即假设在第k次迭代中，前k-1次的系数与基学习器是已知与固定的，则第个基学习器是基于前面个与基学习器的。

  • 

• 损失最小模型，转换为如下问题求解：

  • 在已知与基学习器下，求使得下面公式最小：
    • 

备注：选择指数作为损失函数定义的理由

• 选择指数损失函数的理由是：可以通过数学推导，使得最后结果得到贝叶斯最优错误率：

  • 即对于每个样本都选择后验概率最大的类别。
  • 若指数损失最小化，则分类错误率也将最小化。

• 指数损失函数是分类任务损失函数的一致性替代函数。

  • 由于这个替代函数是单调连续可微函数，因此用它代替0-1损失函数作为优化目标。

• 下面算法过程都是由上面模型推导出来的，具体的推导过程，这里省略。

AdaBoosting算法实现过程

• AdaBoosting除了权重计算不同，最后的结合策略采用的是线性结合策略（加权表决），形成最终的强学习器。

• 假设训练集是：

  • ，
    • 其中为样本容量（样本总数）
    • 表示第个样本；
    • 标识第个样本的分类标签（假设是二分类，并且标签标准化为-1与1）；

1. 初始化样本的权重为：

   • 
     • 为样本容量。
     • ，目前是第个。

2. 使用权重得到新的数据集：

   • ，

3. 使用新的数据集，训练出基学习器：

   • 
     • ，目前是第个。

4. 使用基学习器对样本进行分类，并计算出样本分类的错误率：

   • 
     • ，目前是第个。
     • 如果，则错误率就是

5. 计算每个分类器在最终结合的线性系数：

   • 

     • 其中表示自然对数。

   • 【注意】上面函数有一个良好的性质：

     • 越小，越大，这表示误差率越小的基学习器在最后强学习器中的权重就越大。

6. 计算规范化因子

   • 

7. 新的权值计算

   • 

   • 【注意】上面公式具有一个特点：

     • 约大，就越大，这样下一个基学习器产生的时候，对应样本就强化（越难分类的样本，越强化）。

8. 基学习器的结合策略

   • 
     • 其中表示最终基学习器总数。

9. 为了便于分类，使用符号函数处理下：

   • 因为上面假设我们的类别输出是，我们的强学习器的输出也表示成同样的格式，这个使用符号函数实现：
     • 

AdaBoosting的其他改进

离散与连续

• AdaBoost算法被称为"Discrete AdaBoost"，因为其各个基学习器Gm(x)的输出为{-1,1}
• 连续AdaBoost算法采用的是概率输出，sklearn中的连续与离散指定：SAMME.R，连续的情况，可以提高训练速度。

学习率与过拟合

• 防止学的太快产生过拟合，对每个基学习器乘以一个系数 (又称学习率)，

  • scikit-learn中由learning rate参数指定。

• 一般较小的产生更多的基学习器，可以通过设置基学习器数量来提早终止学习。

样本删除

• Weight trimming:
  • 在AdaBoost的每一轮基学习器训练过程中，只有小部分样本的权重较大，因而能产生较大的影响，而其他大部分权重小的样本则对训练影响甚微。
  • Weight trimming的思想是每一轮迭代中删除那些低权重的样本，只用高权重样本进行训练。具体是设定一个阈值 (比如90%或99%)，再将所有样本按权重排序，计算权重的累积和，累积和大于阈值的权重 (样本) 被舍弃，不会用于训练。注意每一轮训练完成后所有样本的权重依然会被重新计算，这意味着之前被舍弃的样本在之后的迭代中如果权重增加，可能会重新用于训练。

AdaBoosting算法实现

• AdaBoosting的算法实现，我们这里没有子集实现，是抄袭引用网络作者massquantity的文章与代码。

• 文章参考：http://www.cnblogs.com/massquantity/p/9063033.html

• 代码参考：https://github.com/massquantity/ML-algorithms-implemetation

import numpy as np
from sklearn import datasets
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.base import clone
from sklearn.ensemble import AdaBoostRegressor
from sklearn.metrics import zero_one_loss
import time


def clock(func):

    def wrapper(*args, **kwargs):
        start_time = time.time()
        func(*args, **kwargs)
        print("training time: ", time.time() - start_time)

    return wrapper


class AdaBoost(object):
    def __init__(self, M, clf, learning_rate=1.0, method="discrete", tol=None, weight_trimming=None):
        self.M = M
        self.clf = clf
        self.learning_rate = learning_rate
        self.method = method
        self.tol = tol
        self.weight_trimming = weight_trimming

    @clock
    def fit(self, X, y):
        if self.tol is not None:
            X, X_val, y, y_val = train_test_split(X, y, random_state=2)

        w = np.array([1 / len(X)] * len(X))
        self.clf_total = []
        self.alpha_total = []
        former_loss = 1
        count = 0
        tol_init = self.tol

        for m in range(self.M):
            classifier = clone(self.clf)
            if self.method == "discrete":
                if m >= 1 and self.weight_trimming is not None:
                    sort_w = np.sort(w)[::-1]
                    cum_sum = np.cumsum(sort_w)
                    percent_w = sort_w[np.where(cum_sum >= self.weight_trimming)][0]   # 0.999
                    w_fit, X_fit, y_fit = w[w >= percent_w], X[w >= percent_w], y[w >= percent_w]
                    y_pred = classifier.fit(X_fit, y_fit, sample_weight=w_fit).predict(X)
                #    if m % 100 == 0:
                #        print("round {}: {}".format(m, len(w_fit)))
                else:
                    y_pred = classifier.fit(X, y, sample_weight=w).predict(X)
                loss = np.zeros(len(X))
                loss[y_pred != y] = 1
                err = np.sum(w * loss)
                alpha = 0.5 * np.log((1 - err) / err) * self.learning_rate
                w = (w * np.exp(-y * alpha * y_pred)) / np.sum(w * np.exp(-y * alpha * y_pred))

                self.alpha_total.append(alpha)
                self.clf_total.append(classifier)

            elif self.method == "real":
                if m >= 1 and self.weight_trimming is not None:
                    sort_w = np.sort(w)[::-1]
                    cum_sum = np.cumsum(sort_w)
                    percent_w = sort_w[np.where(cum_sum >= self.weight_trimming)][0]
                    w_fit, X_fit, y_fit = w[w >= percent_w], X[w >= percent_w], y[w >= percent_w]
                    y_pred = classifier.fit(X_fit, y_fit, sample_weight=w_fit).predict_proba(X)[:, 1]
                else:
                    y_pred = classifier.fit(X, y, sample_weight=w).predict_proba(X)[:, 1]  
                y_pred = np.clip(y_pred, 1e-15, 1 - 1e-15)
                clf = 0.5 * np.log(y_pred / (1 - y_pred)) * self.learning_rate
                w = (w * np.exp(-y * clf)) / np.sum(w * np.exp(-y * clf))

                self.clf_total.append(classifier)

            '''early stopping'''
            if m % 10 == 0 and m > 300 and self.tol is not None:
                if self.method == "discrete":
                    p = np.array([self.alpha_total[m] * self.clf_total[m].predict(X_val) for m in range(m)])
                elif self.method == "real":
                    p = []
                    for m in range(m):
                        ppp = self.clf_total[m].predict_proba(X_val)[:, 1]
                        ppp = np.clip(ppp, 1e-15, 1 - 1e-15)
                        p.append(self.learning_rate * 0.5 * np.log(ppp / (1 - ppp)))
                    p = np.array(p)

                stage_pred = np.sign(p.sum(axis=0))
                # print("round {}".format(m), zero_one_loss(stage_pred,y_val))
                later_loss = zero_one_loss(stage_pred, y_val)
                if later_loss > (former_loss + self.tol):
                    count += 1
                    self.tol = self.tol / 2  
                #    print(self.tol)
                else:
                    count = 0
                    self.tol = tol_init
                if count == 2:
                    self.M = m - 20
                    print("early stopping in round {}, best round is {}, M = {}".format(m, m - 20, self.M))
                    break
                former_loss = later_loss

        return self

    def predict(self, X):
        if self.method == "discrete":
            pred = np.array([self.alpha_total[m] * self.clf_total[m].predict(X) for m in range(self.M)])

        elif self.method == "real":
            pred = []
            for m in range(self.M):
                p = self.clf_total[m].predict_proba(X)[:, 1]
                p = np.clip(p, 1e-15, 1 - 1e-15)
                pred.append(0.5 * np.log(p / (1 - p)))

        # pred = np.array([0.5 * np.log(self.clf_total[m].predict_proba(X)[:,1] / (1-self.clf_total[m].predict_proba(X)[:,1])) for m in range(self.M)])

        return np.sign(np.sum(pred, axis=0))

    def stage_predict(self, X):
        pred = None
        if self.method == "discrete":
            for alpha, clf in zip(self.alpha_total, self.clf_total):
                current_pred = alpha * clf.predict(X)

                if pred is None:
                    pred = current_pred
                else:
                    pred += current_pred

                yield np.sign(pred)

        elif self.method == "real":
            for clf in self.clf_total:
                p = clf.predict_proba(X)[:, 1]
                p = np.clip(p, 1e-15, 1 - 1e-15)
                current_pred = 0.5 * np.log(p / (1 - p))

                if pred is None:
                    pred = current_pred
                else:
                    pred += current_pred

                yield np.sign(pred)


if __name__ == "__main__":
    X, y = datasets.make_hastie_10_2(n_samples=2000, random_state=1)
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

    model = AdaBoost(M=450, clf=DecisionTreeClassifier(max_depth=1, random_state=1), learning_rate=1.0, method="real",
                     tol=0.01, weight_trimming=0.999)
    model.fit(X_train, y_train)
    pred = model.predict(X_test)

    acc = np.zeros(pred.shape)
    acc[np.where(pred == y_test)] = 1
    accuracy = np.sum(acc) / len(pred)
    print('score: ', accuracy)  # 0.9656667

    base_model = DecisionTreeClassifier(max_depth=1,random_state=1)
    model = AdaBoostRegressor(n_estimators=450, learning_rate=1.0, base_estimator=base_model)
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    acc = np.zeros(pred.shape)
    acc[np.where(pred == y_test)] = 1
    accuracy = np.sum(acc) / len(pred)
    print('score_sklearn: ', accuracy)

training time:  1.3744990825653076
score:  0.918
score_sklearn:  0.734

GradientBoosting算法

• AdaBoost使用的是指数损失，这个损失函数的缺点是对于异常点非常敏感。因而通常在噪音比较多的数据集上表现不佳。

• GradientBoosting算法的特点就是在损失函数上面做了改进：

  • GradientBoosting可以使用任何损失函数 ，从而可以挑选适合数据集的损失函数来度量最终的分类与预测误差。

GradientBoosting算法模型

• GradientBoosting算法模型与AdaBoosting的算法模型是一样的。但有个很大的区别是GradientBoosting最终要得到最强学习器的方式不一样：

  • GradientBoosting算法把最终的最强学习器看成是一个最优的学习器，满足识别误差最小。
  • GradientBoosting算法的过程是函数梯度下降的算法，每次下降的梯度就是一个学习器，为了控制下降速度，提供一个学习率。

• 算法的示意图如下：

  • 梯度提升算法示意图

分类与预测模型

• 上述示意图最后的强学习器可以表示为：
  • 除第一个基学习器，其他基学习器都是使用上一次的残差作为标签作为训练得到。

损失函数模型

• 由于GradientBoosting算法可以使用任何损失函数，所以使用通用模式表示：
  • 
    • 其中表示训练样本的容量。
    • 其中表示最终的学习器

GradientBoosting算法思想

• GradientBoosting算法算法的基本思想是梯度下降，基于权重变量的梯度下降是：

  • 争对损失函数是：

• GradientBoosting算法中下降的是函数：

  • 争对损失函数是：

• 如果考虑第次基学习器的训练，公式表示为：

  • 

• 如果上面是第次使用样本与上次残差作为标签训练出来的基学习器,则第k次最终学习器是：

  • 

• 如果集成的基学习器是m个，则最终的最强学习器是：

  • 

• 总结：

  • 使用机器学习去拟合每次的识别误差，这就是GradientBoosting算法的特点。

GradientBoosting算法过程

• 假设训练集是：
  • ，
    • 其中为样本容量（样本总数）
    • 表示第个样本；
    • 标识第个样本的分类标签（假设是二分类，并且标签标准化为-1与1）；

1. 第一个基学习器

   • 使用指定的算法，与已知样本训练出一个基学习器：

2. 开始次基学习器的训练，使用表示第次

3. 计算第k次的残差：

   • 

   • 注意：

     • 如果损失函数是均方差函数，残差是：

4. 使用数据集训练得到基学习器

   • 这样得到的基学习器是使得最小。

5. 直到m步，完成集成学习过程，得到最强学习器：

   • 
   • 

GradientBoosting算法实现

• 备注:
  • 下面算法直接抄袭引用网络资源：
    • https://github.com/massquantity/ML-algorithms-implemetation/blob/master/Gradient%20Boosting/GradientBoosting_article.py

import numpy as np
from sklearn.base import clone
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
from sklearn.metrics import zero_one_loss
from sklearn import datasets
from scipy import stats

def SquaredLoss_NegGradient(y_pred, y):
    return y - y_pred

def Huberloss_NegGradient(y_pred, y, alpha):
    diff = y - y_pred
    delta = stats.scoreatpercentile(np.abs(diff), alpha * 100)
    g = np.where(np.abs(diff) > delta, delta * np.sign(diff), diff)
    return g

def logistic(p):
    return 1 / (1 + np.exp(-2 * p))

def LogisticLoss_NegGradient(y_pred, y):
    g = 2 * y / (1 + np.exp(1 + 2 * y * y_pred))  # logistic_loss = log(1+exp(-2*y*y_pred))
    return g

def modified_huber(p):
    return (np.clip(p, -1, 1) + 1) / 2

def Modified_Huber_NegGradient(y_pred, y):
    margin = y * y_pred
    g = np.where(margin >= 1, 0, np.where(margin >= -1, y * 2 * (1-margin), 4 * y))
    # modified_huber_loss = np.where(margin >= -1, max(0, (1-margin)^2), -4 * margin)
    return g


class GradientBoosting(object):
    def __init__(self, M, base_learner, learning_rate=1.0, method="regression", tol=None, subsample=None,
                 loss="square", alpha=0.9):
        self.M = M
        self.base_learner = base_learner
        self.learning_rate = learning_rate
        self.method = method
        self.tol = tol
        self.subsample = subsample
        self.loss = loss
        self.alpha = alpha

    def fit(self, X, y):
        # tol为early_stopping的阈值，如果使用early_stopping，则从训练集中分出验证集
        if self.tol is not None:
            X, X_val, y, y_val = train_test_split(X, y, random_state=2)
            former_loss = float("inf")
            count = 0
            tol_init = self.tol

        init_learner = self.base_learner
        y_pred = init_learner.fit(X, y).predict(X)   # 初始值
        self.base_learner_total = [init_learner]
        for m in range(self.M):

            if self.subsample is not None:  # subsample
                sample = [np.random.choice(len(X), int(self.subsample * len(X)), replace=False)]
                X_s, y_s, y_pred_s = X[sample], y[sample], y_pred[sample]  
            else:
                X_s, y_s, y_pred_s = X, y, y_pred

            # 计算负梯度
            if self.method == "regression":
                if self.loss == "square":
                    response = SquaredLoss_NegGradient(y_pred_s, y_s)
                elif self.loss == "huber":
                    response = Huberloss_NegGradient(y_pred_s, y_s, self.alpha)
            elif self.method == "classification":
                if self.loss == "logistic":
                    response = LogisticLoss_NegGradient(y_pred_s, y_s)
                elif self.loss == "modified_huber":
                    response = Modified_Huber_NegGradient(y_pred_s, y_s)

            base_learner = clone(self.base_learner)
            y_pred += base_learner.fit(X_s, response).predict(X) * self.learning_rate
            self.base_learner_total.append(base_learner)

            '''early stopping'''
            if m % 10 == 0 and m > 300 and self.tol is not None:
                p = np.array([self.base_learner_total[m].predict(X_val) for m in range(1, m+1)])
                p = np.vstack((self.base_learner_total[0].predict(X_val), p))
                stage_pred = np.sum(p, axis=0)
                if self.method == "regression":
                    later_loss = np.sqrt(mean_squared_error(stage_pred, y_val))
                if self.method == "classification":
                    stage_pred = np.where(logistic(stage_pred) >= 0.5, 1, -1)  
                    later_loss = zero_one_loss(stage_pred, y_val)

                if later_loss > (former_loss + self.tol):
                    count += 1
                    self.tol = self.tol / 2  
                    print(self.tol)          
                else:
                    count = 0
                    self.tol = tol_init
                    
                if count == 2:
                    self.M = m - 20
                    print("early stopping in round {}, best round is {}, M = {}".format(m, m - 20, self.M))
                   # print("loss: ", later_loss)
                    break
                former_loss = later_loss

        return self

    def predict(self, X):
        pred = np.array([self.base_learner_total[m].predict(X) * self.learning_rate for m in range(1, self.M + 1)])
        pred = np.vstack((self.base_learner_total[0].predict(X), pred))    # 初始值 + 各基学习器
        if self.method == "regression":
            pred_final = np.sum(pred, axis=0)
        elif self.method == "classification":
            if self.loss == "modified_huber":
                p = np.sum(pred, axis=0)
                pred_final = np.where(modified_huber(p) >= 0.5, 1, -1)
            elif self.loss == "logistic":
                p = np.sum(pred, axis=0)
                pred_final = np.where(logistic(p) >= 0.5, 1, -1)
        return pred_final

    def stage_predict(self, X):
        pred = None
        for base_learner in self.base_learner_total:
            current_pred = base_learner.predict(X)

            if pred is None:
                pred = current_pred
            else:
                pred += current_pred * self.learning_rate

            if self.method == "regression":
                yield pred
            elif self.method == "classification":
                if self.loss == "logistic":
                    yield np.where(logistic(pred) >= 0.5, 1, -1)
                elif self.loss == "modified_huber":
                    yield np.where(modified_huber(pred) >= 0.5, 1, -1)


class GBRegression(GradientBoosting):
    def __init__(self, M, base_learner, learning_rate, method="regression", loss="square",tol=None, subsample=None, alpha=0.9):
        super(GBRegression, self).__init__(M=M, base_learner=base_learner, learning_rate=learning_rate, method=method, 
                                            loss=loss, tol=tol, subsample=subsample, alpha=alpha)

class GBClassification(GradientBoosting):
    def __init__(self, M, base_learner, learning_rate, method="classification", loss="logistic", tol=None, subsample=None):
        super(GBClassification, self).__init__(M=M, base_learner=base_learner, learning_rate=learning_rate, method=method,
                                                loss=loss, tol=tol, subsample=subsample)


if __name__ == "__main__":

    X, y = datasets.make_regression(n_samples=20000, n_features=10, n_informative=4, noise=1.1, random_state=1)
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
    model = GBRegression(M=1000, base_learner=DecisionTreeRegressor(max_depth=2, random_state=1), learning_rate=0.1,
                         loss="huber")
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    rmse = np.sqrt(mean_squared_error(y_test, pred))
    print('RMSE: ', rmse)

    X, y = datasets.make_classification(n_samples=20000, n_features=10, n_informative=4, flip_y=0.1, 
                                    n_clusters_per_class=1, n_classes=2, random_state=1)
    y[y==0] = -1
    X_train, X_test, y_train, y_test = train_test_split(X, y)
    model = GBClassification(M=1000, base_learner=DecisionTreeRegressor(max_depth=1, random_state=1), learning_rate=1.0,
                             method="classification", loss="logistic")
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    acc = np.zeros(pred.shape)
    acc[np.where(pred == y_test)] = 1
    accuracy = np.sum(acc) / len(pred)
    print('accuracy logistic score: ', accuracy)

    model = GBClassification(M=1000, base_learner=DecisionTreeRegressor(max_depth=1, random_state=1), learning_rate=1.0,
                             method="classification", loss="modified_huber")
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    acc = np.zeros(pred.shape)
    acc[np.where(pred == y_test)] = 1
    accuracy = np.sum(acc) / len(pred)
    print('accuracy modified_huber score: ', accuracy)

RMSE:  8.45962295476
accuracy logistic score:  0.9374
accuracy modified_huber score:  0.9356

关于GradientBoosting算法的说明

• Gradient Boosting算法理论上可以选择多种不同的学习算法作为基学习器，但实际使用地最多的无疑是决策树，这并非偶然。

  • 决策树有很多优良的特性：
    • 能灵活处理各种类型的数据，包括连续值和离散值；
    • 对缺失值不敏感；
    • 不需要做特征标准化/归一化；
    • 可解释性好；

  -决策树致命缺点是不稳定，导致容易过拟合，因而很多时候准确率不如其他算法。

  • 决策树是非参数模型，这并不意味着其没有参数，而是在训练之前参数数量是不确定的，因此完全生长的决策树有着较大的自由度，能最大化地拟合训练数据。然而单颗决策树是不稳定的，样本数相同的训练集产生微小变动就能导致最终模型的较大差异，即模型的方差大，泛化性能不好。

  • 集成学习的另一代表Bagging可以解决这个问题。而Bagging的拓展算法 —— 随机森林，通过在树内部结点的分裂过程中，随机选取固定数量的特征纳入分裂的候选项，这样就进一步降低了单模型之间的相关性，总体模型的方差也比Bagging更低。

  • 决策树和Gradient Boosting结合诞生了GBDT，GBDT继承了决策树的诸多优点，同时也改进了其缺点。

    • 由于GBDT采用的树都是复杂度低的树，所以方差很小，通过梯度提升的方法集成多个决策树，最终能够很好的解决过拟合的问题。
    • 然而Boosting共有的缺点为训练是按顺序的，难以并行，这样在大规模数据上可能导致速度过慢；
    • 近年来XGBoost和LightGBM的出现都极大缓解了这个问题。

Stacking集成

Stacking算法思路

• Stacking(有时候也称之为stacked generalization)是指训练一个模型用于组合(combine)其他各个模型。

  • 首先训练多个不同的模型；
  • 然后用训练好的多个模型的输出作为为输入来训练一个模型，以得到一个最终的输出。

• 我们称组合的训练模型为元（Meta）学习器；

• 注意：

  • 通常使用单层logistic回归作为组合模型。

Stacking算法过程

• 假设训练集是：
  • ，
    • 其中为样本容量（样本总数）
    • 表示第个样本；

1. 使用不同算法，训练多个基学习器

   • 
     • 其中表示基学习器的个数。

2. 使用基学习器得到新的训练集

   • 

3. 使用作为训练集训练一个元学习器

   • 

4. 就是最后的强学习器

Stacking算法的其他方式-Blend结合

• Blend结合策略：

  • 把训练集分成A，B两个人部分：
  • 使用A训练对个基学习器；
  • 使用多个基学习器多数据集B做预测；
  • 使用预测值作为元学习器的训练集；
  • 最终训练出元学习器。

• 上述算法还有很多变形处理，比如使用交叉验证的方式：

  • 使用折验证，把数据集分成个训练方案；

  • 使用训练集训练，使用测试集得到测试输出，并使用测试集输出作为元学习器的训练样本，元学习器的输入维数为基学习器个数，K折中每折预测值一共K，取K个的平均值作为输出。

  • 使用K折验证输出（K）构成的数据集训练的元学习器就是最终最强学习器。

Stacking算法的实现

import numpy as np
from sklearn.model_selection import KFold


def get_stacking(clf, x_train, y_train, x_test, n_folds=10):
    """
    这个函数是stacking的核心，使用交叉验证的方法得到次级训练集
    x_train, y_train, x_test 的值应该为numpy里面的数组类型 numpy.ndarray .
    如果输入为pandas的DataFrame类型则会把报错"""
    train_num, test_num = x_train.shape[0], x_test.shape[0]
    second_level_train_set = np.zeros((train_num,))
    second_level_test_set = np.zeros((test_num,))
    test_nfolds_sets = np.zeros((test_num, n_folds))
    kf = KFold(n_splits=n_folds)

    for i,(train_index, test_index) in enumerate(kf.split(x_train)):

        x_tra, y_tra = x_train[train_index], y_train[train_index]
        x_tst, y_tst =  x_train[test_index], y_train[test_index]
        clf.fit(x_tra, y_tra)

        second_level_train_set[test_index] = clf.predict(x_tst)

        test_nfolds_sets[:,i] = clf.predict(x_test)

    second_level_test_set[:] = test_nfolds_sets.mean(axis=1)

    return second_level_train_set, second_level_test_set

 

#我们这里使用5个分类算法，为了体现stacking的思想，就不加参数了
from sklearn.ensemble import (RandomForestClassifier, AdaBoostClassifier,
                              GradientBoostingClassifier, ExtraTreesClassifier)
from sklearn.svm import SVC


rf_model = RandomForestClassifier()
adb_model = AdaBoostClassifier()
gdbc_model = GradientBoostingClassifier()
et_model = ExtraTreesClassifier()
svc_model = SVC()


#在这里我们使用train_test_split来人为的制造一些数据
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split


iris = load_iris()
train_x, test_x, train_y, test_y = train_test_split(iris.data, iris.target, test_size=0.2)
train_sets = []
test_sets = []

for clf in [rf_model, adb_model, gdbc_model, et_model, svc_model]:
    train_set, test_set = get_stacking(clf, train_x, train_y, test_x)
    train_sets.append(train_set)
    test_sets.append(test_set)

meta_train = np.concatenate([result_set.reshape(-1,1) for result_set in train_sets], axis=1)
meta_test = np.concatenate([y_test_set.reshape(-1,1) for y_test_set in test_sets], axis=1)

#使用决策树作为我们的次级分类器
from sklearn.tree import DecisionTreeClassifier
dt_model = DecisionTreeClassifier()
dt_model.fit(meta_train, train_y)
df_predict = dt_model.predict(meta_test)
print((df_predict==test_y).sum())

29

XGBoost算法

XGBoost算法介绍

• XGBoost全名叫（eXtreme Gradient Boosting）极限梯度提升

  • 经常被用在一些比赛中，其效果显著。
  • 它是大规模并行boosted tree的工具，它是目前最快最好的开源boosted tree工具包。
  • XGBoost 所应用的算法就是 GBDT（gradient boosting decision tree）的改进，改进Boosting算法串行的局限性；
  • 既可以用于分类也可以用于回归问题中。

• XGBoost是以分类回归树(CART树)进行组合

XGBoost算法模型

预测与分类模型

• 与其他Boosting算法模型一样

• 

  • 是基学习器的个数；
  • 表示第个基学习器；

损失函数模型

• 

  • 其中表示树的复杂度的函数，越小复杂度越低，泛化能力越强；
  • 实际是GradientBoosting梯度下降算法的改进，后面增加了一个树的复杂度正则项（惩罚项）。

• 通常正则项（惩罚项）可以定义如下：

  • 
    • 其中是所有叶子节点输出构成的向量；
    • 其中表示叶子节点的个数；
    • 其中表示叶子切分的难度；
    • 其中是正则化系数；

• 如果使用标识每个节点的输出值，则损失函数可以表示如下：

  • 

• 由损失函数推导的算法以及算法过程，这类不深入推导与解释，我们直接使用第三方的实现模块；

XGBoost算法

官方网站介绍

• 官方地址：

  • https://xgboost.readthedocs.io/en/latest/

• 项目代码地址：

  • https://github.com/dmlc/xgboost

安装

• 所有安装的帮助在官网都有说明：

  • 安装说明

• 安装指令

  • pip3 install xgboost

• 安装过程：

  • 安装过程截图

• 安装后的版本：

  • 安装后的版本显示截图

• 安装注意：
  • 官方安装提示需要编译安装xgboost共享库与动态库；
  • 动态库安装提示

XGBoost的模块结构

• 可以使用多种方式的帮助获得模块结构
  • 交互式帮助

XGBoost使用例子

from numpy import loadtxt
from xgboost import XGBClassifier
from sklearn import datasets as ds
from sklearn.model_selection import train_test_split

data,target = ds.load_iris(return_X_y=True)

data_train, data_test, target_train, target_test = train_test_split(data, target, test_size=0.2)

classifier = XGBClassifier()
classifier.fit(data_train, target_train)
pre = classifier.predict(data_test)

(pre==target_test).sum()

28

• 注意：
  • 其中XGBClassifier类的构造器参数比较多，关注其中与算法有关的参数：

% matplotlib inline
from xgboost import plot_importance
plot_importance(classifier)

<matplotlib.axes._subplots.AxesSubplot at 0x113c134e0>

特征选择可视化

% matplotlib inline
from xgboost import plot_tree,to_graphviz
g = to_graphviz(classifier,num_trees=10)
#存储决策树到图像  
import codecs  
f=codecs.open('xgb_tree.svg', mode='wb')
f.write(g.pipe('svg'));  
f.close()  
# plot_tree 缺少png图像转换环境
g

导出XGBoost树的结果