《机器学习实战》源码解析（六）：AdaBoost

本系列以书中源码为主，稍作修改并添加注释，均实际运行可行。为免后来者踩坑，特此公开！欢迎打赏！
转载请注明出处！

#coding:gbk
from numpy import *
import numpy as np
import matplotlib.pyplot as plt
#作用：载入数据
#输入：无
#输出：数据矩阵，标签向量
#stump:树桩
def loadSimpData():
    datMat = matrix([[1.,2.1],[2.,1.1],[1.3,1.],[1.,1.],[2.,1.]])
    classLabels = [1.0,1.0,-1.0,-1.0,1.0]
    return datMat,classLabels

#作用：通过阈值比较对数据进行分类
#输入：数据矩阵，维度，阈值，阈值不平等性（即关系符）
#输出：返回分类后的标签列表
def stumpClassify(dataMatrix,dimen,threshVal,threshIneq):
    retArray = ones((shape(dataMatrix)[0],1))#构建m*1矩阵，m为dataMatrix第二维的长度，即行数
#下面一步的意义在于，数据的类别有两类，
#通过在阈值左右两边分别采用这两种类别来判断哪一种分类使得错误率更低。详情见下一函数
    if threshIneq == 'lt':#小于等于情况（第一遍）
        retArray[dataMatrix[:,dimen] <= threshVal] = -1.0
    else:#大于情况（第二遍）
        retArray[dataMatrix[:,dimen] >threshVal] = -1.0
    return retArray

#遍历stumpClassify函数所有的可能输入值，并且对每一层特征值生成单层决策树，
#然后从中挑选出最佳的单层决策树
#输入：数据集，类标签向量，数据权值向量
#输出：最佳单层决策树的相关信息，最小错误率，最佳预测结果
def buildStump(dataArr,classLabels,D):
    dataMatrix = mat(dataArr)#可计算化
    labelMat = mat(classLabels).T#可计算化并转置
    m,n = shape(dataMatrix)#获取矩阵各维度的值，此处为行数和列数
    #numSteps:步数；bestStump保存最佳单层决策树的相关信息；
    #bestClasEst:最佳预测结果
    numSteps = 10.0;bestStump ={};bestClasEst = mat(zeros((m,1)))
    #minError取无穷大
    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 = mat(ones((m,1)))
                errArr[predictedVals == labelMat] = 0
                #计算加权错误率
                weightedError = D.T*errArr
                #print("split:dim %d,thresh %.2f,thresh inequal: \
                #%s,the weighted error is %.3f" %(i,threshVal,inequal,weightedError))
                #找到最小错误率的最佳单层决策树，并保存相关信息
                if weightedError < minError:
                    minError = weightedError
                    bestClasEst = predictedVals.copy()
                    bestStump['dim'] = i
                    bestStump['thresh'] = threshVal
                    bestStump['ineq'] = inequal
    return bestStump,minError,bestClasEst

#测试以上代码
datMat,classLabels = loadSimpData()
#D = mat(ones((5,1))/5)
#buildStump(datMat,classLabels,D)

#基于单层决策树的AdaBoost训练过程
#输入参数：数据集，类别标签，迭代次数（默认40）
def adaBoostTrainDS(dataArr,classLabels,numIt=40):
    weakClassArr = []#弱分类，该数组作为返回值，包含三部词典，待会再写
    m = shape(dataArr)[0]#获取数据集第二维个数（行数）
    D = mat(ones((m,1))/m)#初始数据权值均匀化
    #创建一个M*1，值为1的向量，存储类类别估计累计值，初始全为1，最后用作sign函数的输入
    #亦即分类结果，只要进行sign二值化后即为最终分类结果
    aggClassEst = mat(zeros((m,1)))
    for i in range(numIt):
        bestStump,error,classEst = buildStump(dataArr,classLabels,D)
        #print("D:",D.T)
        #alpha公式:alpha=1/2*ln((1-error)/error)
        #max(error.1e-16)是为了确保在没有错误是不会发生溢出
        alpha = float(0.5*log((1.0-error)/max(error,1e-16)))
        bestStump['alpha'] = alpha#alpha加入字典
        weakClassArr.append(bestStump)#字典加入列表
        #print("classEst:",classEst.T)
        #更新D，D是一个概率分布向量，因此，其所有元素之和为1
        #D=(D*exp(-alpha*classLabels*classEst))/D.sum()
        #其中，classLabels是实际样本标签，classEst是分类器分类后得到的预测分类结果
        #multiply是点积，其中，该例中，classLabels:1*m;classEst:m*1
        expon = multiply(-1*alpha*mat(classLabels).T,classEst)
        D = multiply(D,exp(expon))
        D = D/D.sum()
        #更新累计类别估计值
        #aggClassEst记录每个数据点的类别估计累计值
        #最终分类函数H(x) = sign(aggClassEst),
        #而aggClassEst += alpha*classEst
        aggClassEst += alpha*classEst
        #print("aggClassEst:",aggClassEst.T)
        #计算错误率
        aggErrors = multiply(sign(aggClassEst) != mat(classLabels).T,
                            ones((m,1)))
        errorRate = aggErrors.sum()/m
        print("total error:",errorRate,"\n")
        #当总错误率为0时，直接终止循环
        if errorRate == 0.0:
            break
    return weakClassArr,aggClassEst

#AdaBoost分类函数
def adaClassify(datToClass,classifierArr):
    dataMatrix = mat(datToClass)#待分类样例 转换成numpy矩阵
    m = shape(dataMatrix)[0] #获取第2维长度,即数据个数
    aggClassEst =  mat(zeros((m,1)))#aggClassEst记录每个数据点的类别估计累计值
    #遍历弱分类器数组，获取分类结果的带权值并累加
    for i in range(len(classifierArr)):
        #注意此处源码报错，在分类器后加下标[0]解决，原因尚不明，可能和数组结构相关
        #即每个classifier后面需要添加[0],
        #因为这里是根据字符串来索引，而list形式并不支持key
        classEst = stumpClassify(dataMatrix,classifierArr[0][i]['dim']\
                                ,classifierArr[0][i]['thresh']\
                                ,classifierArr[0][i]['ineq'])
        aggClassEst += classifierArr[0][i]['alpha']*classEst
        print("aggClassEst:",aggClassEst)#输出分类结果
    return sign(aggClassEst)#符号函数二值化

#在难数据集上应用
#自适应数据加载函数
def loadDataSet(fileName):
    #获取数据集每一条数据的特征值个数
    numFeat = len(open(fileName).readline().split('\t'))
    dataMat = [];labelMat =[]#用于存储清洗后的数据，以及对应的标签
    fr = open(fileName)
    for line in fr.readlines():#清洗数据
        lineArr = []
        curLine = line.strip().split('\t')
        for i in range(numFeat-1):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat,labelMat

#ROC曲线的绘制以及AUC计算函数
def plotROC(predStrengths,classLabels):
    cur = (1.0,1.0)
    ySum =0.0
    #计算真实结果中的所有正例
    numPosClass = sum(np.array(classLabels) == 1.0)
    yStep = 1/float(numPosClass)
    xStep = 1/float(len(classLabels)-numPosClass)
    #获得排好序的索引
    sortedIndicies = np.argsort(predStrengths)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for index in sortedIndicies.tolist()[0]:
        if classLabels[index] == 1.0:
            delX = 0;delY = yStep;
        else:
            delX = xStep;delY = 0;
            ySum += cur[1]
        ax.plot([cur[0],cur[0]-delX],[cur[1],cur[1]-delY],c='b')
        cur = (cur[0]-delX,cur[1]-delY)
    ax.plot([0,1],[0,1],'b--')
    plt.xlabel("Flase Positive Rate")
    plt.ylabel("True Positive Rate")
    plt.title("ROC curve for AdaBoost Horse Colic Detection System")
    ax.axis([0,1,0,1])
    plt.show()
    print("the Area Under the Curve is:",ySum*xStep)

#测试
datArr,labelArr=loadDataSet(r'选择你自己的数据集存储路径\horseColicTraining2.txt')
classifierArray,aggClassEst = adaBoostTrainDS(datArr,labelArr,10)
plotROC(aggClassEst.T,labelArr)