《统计学习方法》第八章，提升方法2

▶ 回归问题的提升树，采用分段常数函数进行回归，每次选择一个分界点及相应的两侧估值，最小化训练数据在两侧的残差平方和的和，然后将当前残差作为下一次训练的数据，最后把各次得到的分段函数加在一起，就是最终回归结果

● 代码

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

trainDataRatio = 0.3
dataSize = 1000
defaultEpsilon = 0.05                                                           # 需要进行分段的最小区间宽度
randomSeed = 103

def myColor(x):                                                                 # 颜色函数
    r = np.select([x < 1/2, x < 3/4, x <= 1, True],[0, 4 * x - 2, 1, 0])
    g = np.select([x < 1/4, x < 3/4, x <= 1, True],[4 * x, 1, 4 - 4 * x, 0])
    b = np.select([x < 1/4, x < 1/2, x <= 1, True],[1, 2 - 4 * x, 0, 0])
    return [r**2,g**2,b**2]

def dataSplit(dataX, dataY, part):                                              # 将数据集分割为训练集和测试集
    return dataX[:part],dataY[:part], dataX[part:], dataY[part:]

def judgeStrong(x, para):                                                       # 强分类判别函数，调用弱分类判别函数进行线性加和
    return para[1][targetIndex(x, para[0])]    

def targetIndex(x, xList):                                                      # 二分查找 xList 中大于 x 的最小索引
    lp = 0
    rp = len(xList) - 1    
    while lp < rp - 1:        
        mp = (lp + rp) >> 1
        if(xList[mp] >= x):
            rp = mp
        else:
            lp = mp
    return rp

def createData(count = dataSize):                                               # 创建数据
    np.random.seed(randomSeed)       
    X = np.random.rand(count)    
    Y = X * (32 / 3 * (X-1) * (X-1/2) + 1)
    return X, Y
    
def adaBoost(dataX, dataY, weakCount):                                          # 提升训练
    count = len(dataX)
    xSort, ySort = zip( *sorted(zip(dataX, dataY)) )
    xSort = np.array(xSort)
    ySort = np.array(ySort)    

    result = np.zeros([weakCount,3])
    for i in range(weakCount):
        table = np.zeros(count - 1)
        for j in range(count - 1):                                              # 找最佳分割点，损失函数为两侧残差平方和的和
            #table[j] = np.sum( (ySort[:j+1] - np.mean(ySort[:j+1]))**2 ) + np.sum( (ySort[j+1:] - np.mean(ySort[j+1:]))**2 )    
            table[j] = -np.sum(ySort[:j+1])**2 / (j+1) - np.sum(ySort[j+1:])**2 / (count - j - 1)    # 简化一点计算量        
        index = np.argmin(table)                    
        
        valueLeft = np.mean(ySort[:index+1])                                    # 两侧估值
        valueRight = np.mean(ySort[index+1:])                                
        result[i] = np.array([ xSort[index], valueLeft, valueRight ])           # 当前结果暂存，表示当前分类器以 xSort[index] 为分点，左右侧取值分别为 valueLeft 和 valueRight
        ySort[:index+1] -= valueLeft                                            # 两侧新残差
        ySort[index+1:] -= valueRight        
    
    result = np.array(sorted(result, key=lambda x:x[0]))
    #L = np.triu( np.tile(result[:,1], [weakCount+1,1]), 0 )                         # 分段结果加和为最终分类器，注意段数为 weakCount + 1 
    #R = np.tril( np.tile(result[:,2], [weakCount+1,1]), -1)                         # 这里第 i 行 等于 r[0] + r[1] + ... + r[i-1] + l[i] + l[i+1] + ... +l[n]，表示第 i 个分点
    #return np.concatenate(( result[:,0], np.array([np.inf]) )), np.sum(L + R, 1)    # 其中第 0 行是 np.sum(result[:,1])，第 weakCount 行是 np.sum(result[:,2])
    L = np.concatenate((                  np.cumsum(result[:,1].T[::-1])[::-1], np.array([0.0]) )) # 简化一点计算量和空间开销
    R = np.concatenate(( np.array([0.0]), np.cumsum(result[:,2].T) ))
    return np.concatenate(( result[:,0], np.array([np.inf]) )), L + R

def test(weakCount):                                                            # 单次测试
    allX, allY = createData()
    trainX, trainY, testX, testY = dataSplit(allX, allY, int(dataSize * trainDataRatio))
    
    para = adaBoost(trainX, trainY, weakCount)        

    testResult = [ judgeStrong(i, para) for i in testX ]
    errorRatio = np.sum( np.abs(np.array(testResult) -  testY) ) / (dataSize*(1-trainDataRatio))
    print( "weakCount = %d, errorRatio = %f"%(weakCount, round(errorRatio,4)) )         

    fig = plt.figure(figsize=(10, 8))                                           # 画图
    plt.xlim(-0.1, 1.1)
    plt.ylim(-0.1, 1.1)
    XX = [0.0] + para[0][:-1].flatten().tolist() + [1.0]
    YY = [ judgeStrong(i, para) for i in XX ]    
    plt.scatter(testX, testY, color = myColor(1), s = 8, label = "originalData")
    for i in range(weakCount+1):
        plt.plot(XX[i:i+2], [YY[i+1], YY[i+1]],color = myColor(0), label = "regressionData")        
    plt.text(0.1, 1.0, "weakCount = " + str(weakCount) + "
errorRatio = " + str(round(errorRatio,4)), size = 15, ha="center", va="center", bbox=dict(boxstyle="round", ec=(1., 0.5, 0.5), fc=(1., 1., 1.)))       
    R = [ Rectangle((0,0),0,0, color = myColor(1 - i)) for i in range(2) ]
    plt.legend(R, ["originalData", "regressionData"], loc=[0.79, 0.012], ncol=1, numpoints=1, framealpha = 1)

    fig.savefig("R:\weakCount" + str(weakCount) + ".png")
    plt.close()   

if __name__ == '__main__':
    test(1)
    test(2)
    test(3)
    test(4)           
    test(20)            
    test(100)            

● 输出结果，包含训练轮数和预测错误率

weakCount = 1, errorRatio = 0.184200
weakCount = 2, errorRatio = 0.142200
weakCount = 3, errorRatio = 0.115000
weakCount = 4, errorRatio = 0.102500
weakCount = 20, errorRatio = 0.035300
weakCount = 100, errorRatio = 0.016000

● 画图，总数据量 1000，训练集 300，测试集 700

● 画图，当总数据量减为 500，训练集 150 时，出现了过拟合？