AdaBoost python代碼實現

本文參考自：（1）李航《統計學習與方法》  （2）https://github.com/apachecn/AiLearning/blob/master/src/py2.x/ml/7.AdaBoost/adaboost.py 

• 提升方法(boosting)是一種常用的統計學習方法，在分類問題中，他通過改變訓練樣本的權重，學習多個分類器，並將這些分類器進行線性組合，提高分類器的性能

• 具體來說，對於提升方法來說，有兩個問題需要回答：一是在每一輪如何改變訓練數據的權值或概率分布；二是如何將弱分類器組合成一個強分類器。

• 關於第一個問題，AdaBoost的方法事提高那些被前一輪弱分類器錯誤分類樣本的權值，而降低那些被正確分類的樣本的權值。這樣一來，那些沒有得到正確分類的數據，由於其權值的加大而受到后一輪的弱分類器更大的關注

• 對於第2個問題，即弱分類器的組合，AdaBoost采取加權多數表決的方式，具體的，加大分類誤差率小的弱分類器的權值，使其在表決中起較大的作用；減小分類誤差率大的弱分類器的權值，使其在表決中起較小的作用，AdaBoost的巧妙之處在於將這些想法自然且有效的實現在一種算法里
• 下圖是AdaBoost的工作原理：
•  

項目概述：

預測患有疝氣病的馬的存活問題，這里的數據包括368個樣本和28個特征，疝氣病是描述馬胃腸痛的術語，然而，這種病並不一定源自馬的胃腸問題，其他問題也可能引發疝氣病，該數據集中包含了醫院檢測馬疝氣病的一些指標，有的指標比較主觀，有的指標難以測量，例如馬的疼痛級別。另外，除了部分指標主觀和難以測量之外，該數據還存在一個問題，數據集中有30%的值是缺失的。

1.准備數據

def loadDataSet(filename):
    dim = len(open(filename).readline().split('\t'))  #獲取每個樣本的維度(包括標簽)
    data = []
    label = []
    fr = open(filename)
    for line in fr.readlines():  #一行行的讀取
        LineArr = []
        curline = line.strip().split('\t')  #以tab鍵划分，去除掉每個的空格
        for i in range(dim-1):
            LineArr.append(float(curline[i]))
        data.append(LineArr)
        label.append(float(curline[1]))
    return data,label

2. 構建弱分類器：單層決策樹

def buildStump(dataArr, labelArr, D):
    """buildStump(得到決策樹的模型)
    Args:
        dataArr   特征標簽集合
        labelArr  分類標簽集合
        D         最初的樣本的所有特征權重集合
    Returns:
        bestStump    最優的分類器模型
        minError     錯誤率
        bestClasEst  訓練后的結果集
    """
    # 轉換數據
    dataMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).T
    # m行 n列
    m, n = np.shape(dataMat)

    # 初始化數據
    numSteps = 10.0
    bestStump = {}
    bestClasEst = np.mat(np.zeros((m, 1)))
    # 初始化的最小誤差為無窮大
    minError = np.inf

    # 循環所有的feature列，將列切分成 若干份，每一段以最左邊的點作為分類節點
    for i in range(n):
        rangeMin = dataMat[:, i].min()
        rangeMax = dataMat[:, i].max()
        # print 'rangeMin=%s, rangeMax=%s' % (rangeMin, rangeMax)
        # 計算每一份的元素個數
        stepSize = (rangeMax-rangeMin)/numSteps
        # 例如： 4=(10-1)/2   那么  1-4(-1次)   1(0次)  1+1*4(1次)   1+2*4(2次)
        # 所以： 循環 -1/0/1/2
        for j in range(-1, int(numSteps)+1):
            # go over less than and greater than
            for inequal in ['lt', 'gt']:
                # 如果是-1，那么得到rangeMin-stepSize; 如果是numSteps，那么得到rangeMax
                threshVal = (rangeMin + float(j) * stepSize)
                # 對單層決策樹進行簡單分類，得到預測的分類值
                predictedVals = stumpClassify(dataMat, i, threshVal, inequal)
                # print predictedVals
                errArr = np.mat(np.ones((m, 1)))
                # 正確為0，錯誤為1
                errArr[predictedVals == labelMat] = 0
                # 計算 平均每個特征的概率0.2*錯誤概率的總和為多少，就知道錯誤率多高
                # 例如： 一個都沒錯，那么錯誤率= 0.2*0=0 ， 5個都錯，那么錯誤率= 0.2*5=1， 只錯3個，那么錯誤率= 0.2*3=0.6
                weightedError = D.T*errArr
                '''
                dim            表示 feature列
                threshVal      表示樹的分界值
                inequal        表示計算樹左右顛倒的錯誤率的情況
                weightedError  表示整體結果的錯誤率
                bestClasEst    預測的最優結果
                '''
                # print "split: dim %d, thresh %.2f, thresh ineqal: %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

    # bestStump 表示分類器的結果，在第幾個列上，用大於／小於比較，閾值是多少
    return bestStump, minError, bestClasEst

def stumpClassify(dataMat, dimen, threshVal, threshIneq):
    """stumpClassify(將數據集，按照feature列的value進行 二分法切分比較來賦值分類)
    Args:
        dataMat    Matrix數據集
        dimen      特征列
        threshVal  特征列要比較的值
    Returns:
        retArray 結果集
    """
    # 默認都是1
    retArray = np.ones((np.shape(dataMat)[0], 1))
    # dataMat[:, dimen] 表示數據集中第dimen列的所有值
    # threshIneq == 'lt'表示修改左邊的值，gt表示修改右邊的值
    # print '-----', threshIneq, dataMat[:, dimen], threshVal
    if threshIneq == 'lt':
        retArray[dataMat[:, dimen] <= threshVal] = -1.0
    else:
        retArray[dataMat[:, dimen] > threshVal] = -1.0
    return retArray 

3. AdaBoost算法實現

def adaBoostTrainDS(dataArr, labelArr, numIt=40):
    """
    Args:
        dataArr   特征標簽集合
        labelArr  分類標簽集合
        numIt     實例數
    Returns:
        weakClassArr  弱分類器的集合
        aggClassEst   預測的分類結果值
    """
    weakClassArr = []
    m = np.shape(dataArr)[0]  #樣本的個數 
    # 初始化 D，設置每個樣本的權重值，平均分為m份
    W = np.mat(np.ones((m, 1))/m)
    aggClassEst = np.mat(np.zeros((m, 1)))
    for i in range(numIt):
        # 得到決策樹的模型
        bestStump, error, classEst = buildStump(dataArr, labelArr, W)

        # alpha目的主要是計算每一個分類器實例的權重(組合就是分類結果)
        # 計算每個分類器的alpha權重值
        alpha = float(0.5*np.log((1.0-error)/max(error, 1e-16)))
        bestStump['alpha'] = alpha
        # store Stump Params in Array
        weakClassArr.append(bestStump)

        print("alpha=%s, classEst=%s, bestStump=%s, error=%s "% (alpha, classEst.T, bestStump, error))
        # 分類正確：乘積為1，不會影響結果，-1主要是下面求e的-alpha次方
        # 分類錯誤：乘積為 -1，結果會受影響，所以也乘以 -1
        expon = np.multiply(-1*alpha*np.mat(labelArr).T, classEst)
        print('(-1取反)預測值expon=', expon.T)
        # 計算e的expon次方，然后計算得到一個綜合的概率的值
        # 結果發現： 判斷錯誤的樣本，D中相對應的樣本權重值會變大。
        W = np.multiply(W, np.exp(expon))
        W = W/W.sum()

        # 預測的分類結果值，在上一輪結果的基礎上，進行加和操作
        print('當前的分類結果：', alpha*classEst.T)
        aggClassEst += alpha*classEst
        print("疊加后的分類結果aggClassEst: ", aggClassEst.T)
        # sign 判斷正為1， 0為0， 負為-1，通過最終加和的權重值，判斷符號。
        # 結果為：錯誤的樣本標簽集合，因為是 !=,那么結果就是0 正, 1 負
        aggErrors = np.multiply(np.sign(aggClassEst) != np.mat(labelArr).T, np.ones((m, 1)))
        errorRate = aggErrors.sum()/m
        # print "total error=%s " % (errorRate)
        if errorRate == 0.0:
            break
    return weakClassArr, aggClassEst

4.  分類准確率的計算

def adaClassify(datToClass, classifierArr):
    # do stuff similar to last aggClassEst in adaBoostTrainDS
    dataMat = np.mat(datToClass)
    m = np.shape(dataMat)[0]
    aggClassEst = np.mat(np.zeros((m, 1)))
    
    # 循環 多個分類器
    for i in range(len(classifierArr)):
        # 前提： 我們已經知道了最佳的分類器的實例
        # 通過分類器來核算每一次的分類結果，然后通過alpha*每一次的結果 得到最后的權重加和的值。
        classEst = stumpClassify(dataMat, classifierArr[i]['dim'], classifierArr[i]['thresh'], classifierArr[i]['ineq'])
        aggClassEst += classifierArr[i]['alpha']*classEst
        # print aggClassEst
    return np.sign(aggClassEst)

5.  繪制ROC曲線　　

def plotROC(predStrengths, classLabels):
    """plotROC(打印ROC曲線，並計算AUC的面積大小)
    Args:
        predStrengths  最終預測結果的權重值
        classLabels    原始數據的分類結果集
    """
    print('predStrengths=', predStrengths)
    print('classLabels=', classLabels)

    # variable to calculate AUC
    ySum = 0.0
    # 對正樣本的進行求和
    numPosClas = sum(np.array(classLabels)==1.0)
    # 正樣本的概率
    yStep = 1/float(numPosClas)
    # 負樣本的概率
    xStep = 1/float(len(classLabels)-numPosClas)
    # argsort函數返回的是數組值從小到大的索引值
    # get sorted index, it's reverse
    sortedIndicies = predStrengths.argsort()
    # 測試結果是否是從小到大排列
    print('sortedIndicies=', sortedIndicies, predStrengths[0, 176], predStrengths.min(), predStrengths[0, 293], predStrengths.max())

    # 開始創建模版對象
    fig = plt.figure()
    fig.clf() 
    ax = plt.subplot(111)
    # cursor光標值
    cur = (1.0, 1.0)
    # loop through all the values, drawing a line segment at each point
    for index in sortedIndicies.tolist()[0]:
        if classLabels[index] == 1.0:
            delX = 0
            delY = yStep
        else:
            delX = xStep
            delY = 0
            ySum += cur[1]
        # draw line from cur to (cur[0]-delX, cur[1]-delY)
        # 畫點連線 (x1, x2, y1, y2)
        print(cur[0], cur[0]-delX, cur[1], cur[1]-delY)
        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('False positive rate')
    plt.ylabel('True positive rate')
    plt.title('ROC curve for AdaBoost horse colic detection system')
    # 設置畫圖的范圍區間 (x1, x2, y1, y2)
    ax.axis([0, 1, 0, 1])
    plt.show()
    '''
    參考說明：http://blog.csdn.net/wenyusuran/article/details/39056013
    為了計算 AUC ，我們需要對多個小矩形的面積進行累加。
    這些小矩形的寬度是xStep，因此可以先對所有矩形的高度進行累加，最后再乘以xStep得到其總面積。
    所有高度的和(ySum)隨着x軸的每次移動而漸次增加。
    '''
    print("the Area Under the Curve is: ", ySum*xStep)

6. 主函數　　

if __name__ == "__main__":    
    # 馬疝病數據集
    # 訓練集合
    dataArr, labelArr = loadDataSet("./data/horseColicTraining.txt")
    weakClassArr, aggClassEst = adaBoostTrainDS(dataArr, labelArr, 40)
    print(weakClassArr, '\n-----\n', aggClassEst.T)
    # 計算ROC下面的AUC的面積大小
    plotROC(aggClassEst.T, labelArr)
    # 測試集合
    dataArrTest, labelArrTest = loadDataSet("./data/horseColicTest.txt")
    m = np.shape(dataArrTest)[0]
    predicting10 = adaClassify(dataArrTest, weakClassArr)
    errArr = np.mat(np.ones((m, 1)))
    # 測試：計算總樣本數，錯誤樣本數，錯誤率
    print(m, errArr[predicting10 != np.mat(labelArrTest).T].sum(), errArr[predicting10 != np.mat(labelArrTest).T].sum()/m)