参考资料:
李航《统计学习方法》
《机器学习实战》

推荐阅读:
【机器学习】Boosting与AdaBoost分类与回归原理详解与公式推导
集成学习原理小结
集成学习之Adaboost算法原理小结

1. 基本算法思路

提升(boosting)算法实质上是基于“三个臭皮匠赛过诸葛亮”的思想进行的,也就是通过训练多个弱分类器的(有权)组合形成一个强分类器。

对于提升方法来说,有两个问题需要回答:

  • 一是在每一轮如何改变训练数据的权值或概率分布
  • 二是如何将弱分类器组合成一个强分类器

AdaBoost是adaptive boosting(自适应boosting)的缩写,在训练过程中,首先先训练一系列弱分类器,然后将这些分类器组合起来,就形成了强分类器。关于训练数据的权值或概率分布,AdaBoost的做法是提高前一轮弱分类器错误分类样本的权重,而降低被正确分类样本的权值。这样一来,没有得到正确分类的数据,由于其权值的加大而受到后一轮的弱分类器的更大关注。通俗来说,就是“爱哭的孩子有奶喝”。对于弱分类器的组合,AdaBoost采取加权多数表决的方法,加大分类误差概率小的弱分类器的权值,使其在表决中起较大的作用。

如下图所示:

左边是数据集,其中直方图的不同宽度表示每个样例上的不同权重。在经过一个分类器之后,加权的预测结果会通过三角形中的alpha值进行加权。每个三角形中输出的加权结果在圆形中求和,从而得到最终的输出结果。

假设数据集中有N个训练样本,每个样本都由一个权重\omega_{ki},表示第k次迭代时第i个样本的权重值,这些权重值组成一个向量D_k,表示第k次迭代时的权重向量。

D_k=(\omega_{k1},\omega_{k2},\omega_{k3},\cdots,\omega_{ki},\cdots,\omega_{kN})

关键步骤的计算如下:

(1)初始化样本权重D_1,一开始,这些权重都初始化成相等值。即D_1=(\omega_{11},\omega_{12},\omega_{13},\cdots,\omega_{1i},\cdots,\omega_{1N}),其中\omega_{1i}=\frac{1}{N}

(2)迭代训练。

m=1,2,\cdots,M,每个弱分类器,使用具有权值分布D_m的训练数据集学习,得到基本分类器

G_m(x)=\chi \to {-1,+1}

每次训练都会得到一个错误率\varepsilon

e_m=\sum_{i=1}^N P(G_m(x_i) \neq y_i)=\sum_{i=1}^N \omega_{mi}I(G_m(x_i) \neq y_i)

计算G_m(x)的系数

\alpha_m = \frac{1}{2} \log \frac{1-e_m}{e_m}

然后在同一数据集上再次训练弱分类器。在分类器的第二次训练当中,将会重新调整每个样本的权重,其中第一次分对的样本的权重将会降低,而第一次分错的样本的权重将会提高。训练数据集的权值分布更新:

D_{m+1}=(\omega_{m+1,1},\cdots,\omega_{m+1,i},\cdots,\omega_{m+1,N})

\omega_{m+1,i}=\frac{\omega_{mi}}{Z_m} \exp (-\alpha_m y_i G_m(x_i)), \text{其中} i=1,2,\cdots,N

这里的Z_m为规范化因子

Z_m=\sum_{i=1}^N \omega_{mi} \exp (-\alpha_m y_i G_m(x_i))

这个规范化因子使D_{m+1}成为一个概率分布。

(3)弱分类器的基本组合

设弱分类器的线性组合为

f(x)=\sum_{m=1}^M \alpha_m G_m(x)

则可以得到最终的分类器

G(x)=sign(f(x))=sign(\sum_{m=1}^M \alpha_m G_m(x))

2. AdaBoost算法流程

AdaBoost算法流程如下:

(1)假设训练数据集具有均匀的权值分布,即初始化每个训练样本在基本分类器的学习中的作用相同。

(2)AdaBoost反复学习基本分类器,在每一轮m=1,2,\cdots,M,依次进行如下操作。
(a)使用当前分布D_m加权的训练数据集,学习基本分类器G_m(x)
(b)计算基本分类器G_m(x)在加权训练数据集上的分类误差率:

e_m=\sum_{i=1}^N P(G_m(x_i) \neq y_i)=\sum_{G_m(x_i) \neq y_i} \omega_{}mi

(c)计算基本分类器G_m(x)的系数\alpha_m
(d)更新训练数据的权值,为下一轮作准备。

(3)得到基本分类器的线性组合f(x)

3. AdaBoost算例

假设训练数据为

序号 1 2 3 4 5 6 7 8 9 10
x 0 1 2 3 4 5 6 7 8 9
y 1 1 1 -1 -1 -1 1 1 1 -1

假设弱分类器由x>vx{<}v产生,阈值v使该分类器在训练数据集上分类误差率最低。

def basicClassifierTraining(x, y, D): ## 训练基本分类器
    v = 0.5
    error = 1000000

    best_v = None
    best_direct = None
    best_pred_y = None

    while v<11:
        pred_y_pos = np.array([1 if xi > v else -1 for xi in x])
        emPositive = sum([D[i] for i in range(len(x)) if pred_y_pos[i] != y[i]])

        pred_y_neg = np.array([-1 if xi > v else 1 for xi in x])
        emNegetive = sum([D[i] for i in range(len(x)) if pred_y_neg[i] != y[i]])

        if emPositive<emNegetive:
            direct = 'positive'
            em = emPositive
            pred_y = pred_y_pos
        else:
            direct = 'negetive'
            em = emNegetive
            pred_y = pred_y_neg

        if em<error:
            error = em
            best_v = v
            best_direct = direct
            best_pred_y = pred_y
        v = v + 1
    return best_v, error, best_direct, best_pred_y

在本算例中,

D_1=(\omega_{11},\omega_{12},\cdots,\omega_{110})

\omega_{1i}=0.1,\text{其中} i=1,2,\cdots,10

进行迭代操作如下。

m=1
(a)在权值分布为D_1的训练数据集上,阈值v=2.5时分类误差最低,所以基本分类器为

G_1(x)=
\begin{cases}
1, x{<}2.5 \\
-1, x>2.5 \\
\end{cases}

(b)G_1(x)在训练数据集上的误差率为e_1=P(G_1(x_i) \neq y_i)=0.3
(c)计算G_1(x)的系数:\alpha_1=\frac{1}{2} \log \frac{1-e_1}{e_1}=0.4236

def calalpha(error): ## 计算权重
    alpha = 0.5*np.log((1-error)/error)
    return alpha

(d)更新训练数据的权值分布:

D_2=(\omega_{21},\cdots,\omega_{2i},\cdots,\omega_{210})

\omega_{2i}=\frac{\omega_{1i}}{Z_1} \exp (-\alpha_1 y_i G_1(x_i)),\text{其中}i=1,2,\cdots,10

def updateD(D, alpha, v, x, y, best_pred_y): ## 更新D值
    newD = []
    Z = 0
    for i in range(len(D)):
        w = D[i]*np.exp(-alpha*y[i]*best_pred_y[i])
        newD.append(w)
        Z += w
    out = []
    for d in newD:
        out.append(d/Z)
    return out

D_2=(0.07143,0.07143,0.07143,0.07143,0.07143,0.07143,0.16667,0.16667,0.16667,0.07143)

f_1(x)=0.4236G_1(x)

训练器sign[f(x)]在训练数据集上有3个误分类点(点7,8,9)

def sign(x):
    if x>0:
        return 1
    if x<0:
        return -1

def calerrorpoint(gx, x, y):
    ## gx为组合分类器,是一个列表
    ## 列表的每个元素是一个二元组
    ## 二元组的第一个元素表示该分类器的权重alpha,第二个元素表示阈值v
    num = 0
    for i in range(len(x)):
        pred_y = 0
        for classifier in gx:
            alpha = classifier[0]
            v = classifier[1]
            direct = classifier[2]
            if direct == 'positive':
                if x[i]>v:
                    pred_y += 1
                else:
                    pred_y += -1
            else:
                if x[i]>v:
                    pred_y += -1
                else:
                    pred_y += 1
        ## 进行sign函数计算]
        pred_y = sign(pred_y)
        if pred_y != y[i]:
            num +=1
    return num

m=2

同样循环上述的操作(注意此时的训练数据权值为D_2),可以得到基本分类器为

G_1(x)=
\begin{cases}
1, x{<}8.5 \\
-1, x>8.5 \\
\end{cases}

可以得到训练器f(x)=0.4236G_1(x)+0.6496G_2(x),在训练数据集上有3个误分类点。

m=3

可以得到基本分类器为

G_1(x)=
\begin{cases}
1, x{>}5.5 \\
-1, x{<}5.5 \\
\end{cases}

可以得到训练器f(x)=0.4236G_1(x)+0.6496G_2(x)+0.7514G_3(x),在训练数据集上有0个误分类点。

所以最终分类器为G(x)=sign[f(x)]=sign[0.4236G_1(x)+0.6496G_2(x)+0.7514G_3(x)]

主程序如下:

if __name__ == '__main__':
    x = [0,1,2,3,4,5,6,7,8,9]
    y = [1,1,1,-1,-1,-1,1,1,1,-1]

    D = [0.1]*10

    ## 初始化组合分类器, 每个元素是一个三元组
    ## 三元组的第一个元素表示该分类器的权重alpha,第二个元素表示阈值v,第三个元素表示比较的方向
    gx = [] 

    num = 10
    m = 1

    while num>0 and m<4:
        print(str(m) + '-'*10)
        v, error, direct, best_pred_y = basicClassifierTraining(x, y, D)
        alpha = calalpha(error)
        D = updateD(D, alpha, v, x, y, best_pred_y)
        gx.append((alpha, v, direct))
        print(v, error, direct, alpha)
        print(D)
        num = calerrorpoint(gx, x, y)
        print(num)
        m = m + 1

最后的输出:

1----------
2.5 0.30000000000000004 negetive 0.4236489301936017
[0.07142857142857142, 0.07142857142857142, 0.07142857142857142, 0.07142857142857142, 0.07142857142857142, 0.07142857142857142, 0.16666666666666663, 0.16666666666666663, 0.16666666666666663, 0.07142857142857142]
3
2----------
8.5 0.21428571428571427 negetive 0.6496414920651304
[0.04545454545454546, 0.04545454545454546, 0.04545454545454546, 0.16666666666666669, 0.16666666666666669, 0.16666666666666669, 0.10606060606060606, 0.10606060606060606, 0.10606060606060606, 0.04545454545454546]
6
3----------
5.5 0.18181818181818185 positive 0.752038698388137
[0.12499999999999996, 0.12499999999999996, 0.12499999999999996, 0.10185185185185185, 0.10185185185185185, 0.10185185185185185, 0.0648148148148148, 0.0648148148148148, 0.0648148148148148, 0.12499999999999996]
0

4. AdaBoost回归

以下算法流程摘自集成学习之Adaboost算法原理小结

Adaboost R2回归算法过程:

输入:样本集T=\{(x_1,y_1),(x_2,y_2),\cdots, (x_m,y_m)\},弱学习器算法G(x), 弱学习器迭代次数K

输出:最终的强学习器f(x)

(1)初始化样本权重为

D_1=(\omega_{11},\omega_{12},\cdots,\omega_{1m});\quad \text{其中 }\omega_{1i}=\frac{1}{m};i=1,2,\cdots,m

(2)对于k=1,2,\cdots,K

使用具有权重D_k的权重来训练样本,得到弱学习器G_k(x)

计算训练集上的最大误差E_k=\max |y_i-G_k(x_i)|,i=1,2,\cdots,m

计算每个样本的相对误差:

线性误差——e_{ki}=\frac{|y_i-G_k(x_i)|}{E_k}

平方误差——e_{ki}=\frac{(y_i-G_k(x_i))^2}{E_k^2}

指数误差——e_{ki}=1-\exp (\frac{-|y_i-G_k(x_i)|}{E_k})

之后,计算回归误差率e_k=\sum_{i=1}^m w_{ki} e_{ki}

然后,计算弱学习器的系数\alpha_k = \frac{e_k}{1-e_k}

接着,更新样本集的权重分布为w_{k+1,i}=\frac{w_{ki}}{Z_k} \alpha_k^{1-e_{ki}},其中Z_k = \sum_{i=1}^m w_{ki} \alpha_{k}^{1-e_{ki}}

(3)构建最后的强学习器为f(x)=G_{k^{\ast}}(x),其中G_{k^{\ast}}(x)是所有\ln \frac{1}{\alpha_k},k=1,2,\cdots,K的中位数值对应序号k^{\ast}对应的弱分类器。

5. 提升树

提升树(boosting tree)是以分类树或者回归树为基本分类器的提升方法,被认为是统计学习中性能最好的方法之一。对分类问决策树是二叉分类树,对回归问题决策树是二叉回归树。

5.1 提升决策树算法

单层决策树(decision stump,也称决策树桩)是一种简单的决策树(关于决策树可以看我之前的文章——【机器学习系列】决策树)。我们前面提到的那个弱分类器——x>vx{<}v就可以看成是一个决策树桩。

提升树模型可以看成是决策树的加法模型(即基函数的线性组合):

f_M(x)=\sum_{m=1}^M T(x;\Theta_m)

其中,T(x;\Theta_m)表示决策树,\Theta_m表示决策树的参数,M为树的个数。

提升决策树算法采用前向分步算法。首先确定初识提升树f_0(x)=0,第m步的模型是

f_m(x)=f_{m-1}(x)+T(x;\Theta_m)

其中,f_{m-1}(x)为当前模型,经过经验风险极小化确定下一棵决策树的参数\Theta_m

\hat{\Theta}_m = \arg \min_{\Theta_m} \sum_{i=1}^N L(y_i,f_{m-1}(x_i)+T(x_i;\Theta_m))

5.2 基于单层决策树构建弱分类器

首先,我们先来看一下单层决策树这一弱分类器的生成和择优。

import numpy as np

def stumpClassify(dataMatrix, dimen, threshVal, threshIneq):
    retArray = np.ones((np.shape(dataMatrix)[0], 1))
    if threshIneq == 'lt':
        retArray[dataMatrix[:, dimen] <= threshVal] = -1.0
    else:
        retArray[dataMatrix[:, dimen] > threshVal] = -1.0
    return retArray

首先,stumpClassify()是通过阈值比较对数据进行分类的。所有在阈值一边的数据会分到类别-1,而在另外一边的数据分到类别+1。这个函数一共有四个输入参数,dataMatrix为输入的特征值,dimen表示分类是依据第几维特征进行的,threshVal为分类阈值,threshIneq为分类模式——‘lt’表示小于阈值的归为-1类,‘gt’表示大于阈值的归为-1类。

def buildStump(dataArr, classLabels, D):
    dataMatrix = np.mat(dataArr)
    labelMat = np.mat(classLabels).T
    m, n = np.shape(dataMatrix)

    numsteps = 10.0
    bestStump = {}
    bestClassEst = np.mat(np.zeros((m, 1)))
    minError = np.inf

    for i in range(n):
        rangeMin = dataMatrix[:,i].min()
        rangeMax = dataMatrix[:,i].max()
        stepSize = (rangeMax-rangeMin)/numsteps

        for j in range(-1, int(numsteps)+1):
            for inequal in ['lt', 'gt']:
                threshVal = (rangeMin + float(j) * stepSize)
                predictedVals = stumpClassify(dataMatrix, i, threshVal, inequal)
                errArr = np.mat(np.ones((m,1)))
                errArr[predictedVals == labelMat] = 0
                weightedError = D.T * errArr
                print("spit: dim %d, thresh %.2f, thresh ineqal: %s, the weighted error is %.3f"%(i, threshVal, inequal, weightedError))

                if weightedError < minError:
                    minError = weightedError
                    bestClassEst = predictedVals.copy()
                    bestStump['dim'] = i
                    bestStump['thresh'] = threshVal
                    bestStump['ineq'] = inequal
    return bestStump, minError, bestClassEst

buildStump()函数会遍历stumpClassify()函数所有的可能输入值,并找到数据集上最佳的单层决策树。这个函数一共有三层循环,第一层主要是遍历所有的特征维度,第二层以一定的步长进行该特征维度上的取值范围,第三层遍历阈值比较符号——‘gt’或者’lt’,从而找到最佳的单层决策树。

最佳决策树的判断标准为权重向量——D,表示该决策树中各个数据的权重。和前面提到的弱分类器一样,初始情况下,每个训练样本数据的权重是一样的:

m, n = np.shape(np.mat(dataArr))
D = np.mat(np.ones((m,1))/m)

在每一次迭代训练中,我们还是一样地更新\alpha值和权重向量D

def adaBoostTrainDS(dataArr, classLabels, numIt = 40):
    weakClassArr = []
    m = np.shape(dataArr)[0]
    D = np.mat(np.ones((m, 1))/m)
    aggClassEst = np.mat(np.zeros((m,1)))

    for i in range(numIt):
        bestStump, error, classEst = buildStump(dataArr, classLabels, D)
        print("D: ", D.T)

        alpha = float(0.5*np.log((1.0-error)/max(error, 1e-16)))
        bestStump['alpha'] = alpha
        weakClassArr.append(bestStump)
        print("classEst: ", classEst.T)

        ## 更新D值
        expon = np.multiply(-1*alpha*np.mat(classLabels).T, classEst)
        D = np.multiply(D, np.exp(expon))
        D = D/D.sum()

        ## 错误率累加
        aggClassEst += alpha * classEst
        print("aggClassEst: ", aggClassEst.T)
        aggErrors = np.multiply(np.sign(aggClassEst) != np.mat(classLabels).T, np.ones((m,1)))
        errorRate = aggErrors.sum()/m
        print('total error: ', errorRate, "\n")
        if errorRate == 0.0:
            break
    return weakClassArr

我们可以自己准备一些数据集来测试一下:

dataArr = [[1],[2],[3],[4],[5],[6],[7],[8],[9],[10]]
classLabels = [1, 1, -1, 1, 1, -1, 1, 1, 1, -1]

weakClassArr = adaBoostTrainDS(dataArr, classLabels, numIt = 40)

需要注意的是,虽然都是二分类,但是上面的代码要求数据的标签应该是+1和-1,而不是0和1。

6. 回归问题——梯度提升回归(GBR, gradient boosting regression)

6.1 梯度提升算法

输入:训练数据集T=\{(x_1, y_1),(x_2,y_2),\cdots,(x_N,y_N)\}x_i \in \chi \subseteq R^ny_i \in \mathcal{Y} \subseteq R^n;损失函数L(y,f(x))

输出:回归树\tilde{f}(x)

(1)初始化

f_0(x)=\arg \min_c \sum_{i=1}^N L(y_i, c)

(2)对m=1,2,\cdots, M
i=1,2,\cdots,N,计算损失函数的负梯度在当前模型的值,将它作为残差的估计:

r_{mi} = -[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)}]_{f(x)=f_{m-1}(x)}

对于平方损失函数,它就是通常所说的残差;对于一般损失函数,它就是残差的近似值。

r_{mi}拟合一个回归树,得到第m棵树的叶结点区域R_{mj}j=1,2,\cdots,J,以拟合残差的近似值。
接下来,利用线性搜索估计叶结点区域的值,使损失函数极小化。即,对j=1,2,\cdots,J,计算

c_{mj}=\arg \min_{c} \sum_{x_i \in R_{mj}} L(y_i,f_{m-1}(x_i)+c)

最后,更新回归树f_m(x)=f_{m-1} (x)+\sum_{j=1}^J c_{mj} I(x \in R_{mj})

(3)得到回归树

\hat{f}(x)=f_M(x)=\sum_{m=1}^M \sum_{j=1}^J c_{mj} I(x \in R_{mj})

6.2 案例分析

案例来源于sklearn的官网教程——Gradient Boosting regression

参考资料:
sklearn的train_test_split()各函数参数含义解释(非常全)

数据集共有442个样本,每个输入有10个特征值。

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, ensemble
from sklearn.inspection import permutation_importance
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

diabetes = datasets.load_diabetes()
X, y = diabetes.data, diabetes.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, random_state=13)

首先,对数据集进行预处理,test_size=0.1表示训练集占数据集的0.1。test_size若为浮点时,表示测试集占总样本的百分比;若为整数时,表示测试样本样本数;若为None时,test size自动设置成0.25。random_state=13表示随机数的种子数为13;若为None时,每次生成的数据都是随机,可能不一样;若为整数时,每次生成的数据都相同。

在sklearn中使用GBR模型特别简单,只需要调用一下函数GradientBoostingRegressor(),设置一下相关参数就好了。关于数据的调参可以看Gradient Boosting Regressor机器学习超参数调整

params = {'n_estimators': 500,
          'max_depth': 4,
          'min_samples_split': 5,
          'learning_rate': 0.01,
          'loss': 'ls'}
reg = ensemble.GradientBoostingRegressor(**params)
reg.fit(X_train, y_train)