機器學習之局部加權、嶺回歸和前向逐步回歸

　　都說萬事開頭難，可一旦開頭，就是全新的狀態，就有可能收獲自己未曾預料到的成果。記錄是為了更好的監督、理解和推進，學習過程中用到的數據集和代碼都將上傳到github

　　回歸是對一個或多個自變量和因變量之間的關系進行建模，求解的一種統計方法，之前的博客中總結了在線性回歸中使用最小二乘法推導最優參數的過程和logistic回歸，接下來將對最小二乘法、局部加權回歸、嶺回歸和前向逐步回歸算法進行依次說明和總結

1. 用線性回歸找到最佳擬合曲線

（1）回歸的一般方法

　　收集數據：采用任意方法收集數據

　　准備數據：回歸需要數值型數據，標稱型數據將被轉成數值型數據

　　分析數據：繪出數據的可視化二維圖將有助於對數據做出理解和分析，在采用縮減法求得新回歸系數之后，可以將新擬合線繪在圖上作為對比

　　訓練算法：找到回歸系數

　　測試算法：使用R²或者預測值和數據的擬合度，來分析模型的效果

　　使用算法：使用回歸，可以在給定輸入的時候預測出一個數值，這是對分類方法的提升，因為這樣可以預測連續性數據而不僅僅是離散的類別標簽

（2）需求

　　對dataset/ex0.txt文件中的數據集采用線性回歸找到最擬合直線，數據集：https://github.com/lizoo6zhi/MachineLearn.samples/tree/master/dataset

（3）繪制數據的可視化二維圖

　　

（4）通過最小二乘法得到最優參數

　　前面我們已經推導過使用最小二乘法求最優參數：

　　代碼實現：

def stand_regres(dataset,label):
    """
    標准線性回歸函數
    dataset:list
    label:list
    """
    X = np.mat(dataset)
    Y = np.mat(label).T

    temp = X.T * X
    if np.linalg.det(temp) == 0.0:  #計算行式的值
        print("X.T*X沒有逆矩陣")
        return 
    else:
        temp = temp.I * X.T * Y
    return temp

　　注意上面代碼中紅色標記處，通過計算行列式的值判斷X.T*X是否有逆矩陣

（5）通過計算出的最優參數對樣本進行回歸測試，並繪制擬合曲線

def plot_stand_regress(dataset,label,ws,ax):
    """
    繪制數據集和擬合曲線
    """
    dataset_matrix = np.mat(dataset)
    x = dataset_matrix[:,1]
    y = label
    cal_y = dataset * ws
    ax.scatter(x.T.getA(),y,s=10,c='red')
    ax.plot(x,cal_y)

（6）在python中，NumPy庫提供了相關系數的計算方法，可以通過corrcoef(yEstimate,yActual)來計算預測值與真實值的相關性

print(np.corrcoef(cal_y.T,np.mat(y)))

　　從上圖可以看出，預測值與真實值之間的相關系數為0.98647356

　　最佳擬合直線方法將數據視為直線進行建模，具有十分不錯的效果，但若是數據的離散程度很大呢，那該怎么辦？

2. 局部加權線性回歸

　　我們現實生活中的很多數據不一定都能用線性模型描述。比如房價問題，很明顯直線非但不能很好的擬合所有數據點，而且誤差非常大，但是一條類似二次函數的曲線卻能擬合地很好。為了解決非線性模型建立線性模型的問題，我們預測一個點的值時，對這個點附近的點設置權值。基於這個思想，便產生了局部加權線性回歸算法。在這個算法中，其他離一個點越近，權重越大，對回歸系數的貢獻就越多。損失函數如下：
　　
　　其中w(i)是權重，它根據要預測的點與數據集中的點的距離來為數據集中的點賦權值，當某點離要預測的點越遠，其權重越小，否則越大。一個比較好的權重函數如下：

　　
 
　　該函數稱為指數衰減函數，其中k為波長參數，它控制了權值隨距離下降的速率，該函數形式上類似高斯分布(正態分布)，但並沒有任何高斯分布的意義。該算法解出回歸系數如下：

　　

注：局部加權線性回歸中的上述內容摘抄自 https://blog.csdn.net/herosofearth/article/details/51969517，非常感謝作者分享！

（1）需求

　　使用局部加權線性回歸對dataset/ex0.txt文件中的數據集進行曲線擬合

（2）實現

def lwrw(test_point, xarray, yarray, k = 1.0):
    """
    局部加權回歸函數(locally weighted linear regression)
    dataset:list
    label:list
    """
    mat_x = np.mat(xarray)
    mat_y = np.mat(yarray).T
    row,col = np.shape(mat_x)    
    weights = np.mat(np.eye(row))
    for i in range(row):
        diffmat = test_point - mat_x[i,:]
        weights[i,i] = np.exp((diffmat*diffmat.T)/(-2.0*k*k))
    temp = mat_x.T * weights * mat_x
    if np.linalg.det(temp) == 0.0:
        print("mat_x.T * weights * X沒有逆矩陣")
        return 
    else:
        temp = temp.I * mat_x.T * weights * mat_y
    return test_point*temp

（3）測試

　　分別設置k為0.01和0.003，當k=1時擬合效果和最小二乘法擬合效果一樣

def plot_lwrw_regress(dataset,label,cal_y_list,ax):
    """
    繪制數據集和擬合曲線
    """
    dataset_matrix = np.mat(dataset)
    x = dataset_matrix[:,1]
    y = label
    ax.scatter(x.T.getA(),y,s=10,c='red')
    #對dataset進行排序
    sorted_index = x.argsort(0)
    sorted_x = dataset_matrix[sorted_index][:,0,:] #關鍵
    ax.plot(sorted_x[:,1],cal_y_list[sorted_index])

if __name__ == "__main__":
    dataset,label = load_dataset('dataset/ex0.txt')
    fig = plt.figure()
    
    #當k為0.01時
    for row in range(rows):
        cal_y = lwrw(dataset[row],dataset,label,0.01)
        x_y_map[row] = cal_y
    ax = fig.add_subplot(323)
    ax.text(0.05,4.2,'k=0.01')
    plot_lwrw_regress(dataset,label,x_y_map,ax)
    
    #當k為0.003
    for row in range(rows):
        cal_y = lwrw(dataset[row],dataset,label,0.003)
        x_y_map[row] = cal_y
    ax = fig.add_subplot(325)
    ax.text(0.05,4.2,'k=0.003')
    plot_lwrw_regress(dataset,label,x_y_map,ax)
    
    plt.show()

　　效果：

　　大家可以通過不同的k系數測試效果，局部加權線性回歸存在一個問題就是增加了計算量，因為它對每個點做預測時都必須使用整個數據集

3. 嶺回歸

　　思考下，如果數據的特征比樣本還多應該怎么辦？是否還可以使用線性回歸和之前的方法來做預測？

　　答案是否定的，即不能再使用前面介紹的方法，這是因為在計算（X.T*X)的逆矩陣的時候會出錯

　　有逆矩陣的矩陣一定是方正，但方正不一定有逆矩陣，只有方正的行列式不為0時才有逆矩陣

　　為了解決整個問題，統計學家引入了嶺回歸（ridge regression）的概念，簡單的說，嶺回歸就是在X.T * X上加一個λI使得矩陣非奇異，進而確保（X.T * X + λI）一定可逆，其中I為m*m的單位矩陣，λ是用戶定義的一個數值，我們可以通過測試使用不同的λ參數來確定最優λ

　　嶺回歸與多項式回歸唯一的不同在於代價函數上的差別。嶺回歸的代價函數如下：

　　

　　為了方便計算導數，通常也寫成下面的形式：

　　

　　上式中的w是長度為n的向量，不包括截距項的系數θ0；θθ是長度為n+1的向量，包括截距項的系數θ0；m為樣本數；n為特征數

　　嶺回歸的代價函數仍然是一個凸函數，因此可以利用梯度等於0的方式求得全局最優解（正規方程）：

　　

（1）需求

　　通過可視化圖形描述λ與回歸系數之間的關系

（2）實現

　　通過以e為底的指數函數設置λ的值，如：exp(i-10)，i為[0,30]，同公式計算不同λ下回歸系數的變化

（3）代碼

def ridge_regres(xarray,yarray,lam=0.2):
    """
    嶺回歸
    處理當特征比樣本點還多，也酒水說輸入數據的矩陣X不是滿秩矩陣，非滿秩矩陣在求逆時會出現問題
    w = (X.T*X + rI).I * X.T * Y
    """
    xmat = np.mat(xarray)
    ymat = np.mat(yarray).T
    xtx = xmat.T * xmat
    denom = xtx + np.eye(np.shape(xmat)[1])*lam
    if np.linalg.det(denom) == 0.0:
        print('denom 沒有逆矩陣')
    ws = denom.I * xmat.T * ymat
    return ws
    
def ridge_test(xarray,yarray,ax):
    xmat = np.mat(xarray)
    row,col = np.shape(xmat)
    ymat = np.mat(yarray)
    xmean,xvar = data_normalize(xarray) #對數據進行標准化處理
    np.seterr(divide='ignore',invalid='ignore')
    xmat = (xmat - xmean) / xvar
    num_test_pts = 30
    wmat = np.ones((num_test_pts,col))
    
    lam_list = []
    for i in range(num_test_pts):
        lam = np.exp(i-10)
        lam_list.append(i-10)
        ws = ridge_regres(xarray,yarray,lam)
        wmat[i,:] = ws.T
    
    #繪制不同lam情況下參數變化情況
    for j in range(np.shape(wmat)[1]):
        ax.plot(lam_list,wmat[:,j])

　　效果：

　　x軸為i的取值，y軸為回歸系數值，從圖可以看出當λ非常小時，系數和普通回歸一樣，當λ非常大時，所有回歸系數都縮減為0，因此，我們可以在中間某處找到是的預測的結果最好的λ值

　　以下為λ=0.02的效果圖

　　這樣我們可以結合局部加權和嶺回歸計算最優回歸系數，還有其他一些縮減方法，如lasso、LAR、PCA回歸以及子集選擇等，這些方法不僅可以提高預測精確率，而且還可以解釋回歸系數，將后續進行總結

　　再使用嶺回歸對數據集dataset/abalone.txt測試，看看回歸系數變化

4. 前向逐步回歸

　　前向逐步回歸算法可以得到與lasso差不多的效果，但更加簡單。前向回歸算法屬於一種貪心算法，即每一步都盡可能減少誤差，一開始，所有的權重都設為1，然后每一步所作的決策是對某個權重增加或減少一個很小的值

　　偽代碼如下：　

數據標准化，使其分布滿足0均值和單位方差
在每輪迭代過程中：
    設置當前最小誤差lowestError為正無窮
    對每個特征：
        增大或減少：
            改變一個系數得到一個新的w
            計算新w下的誤差
            如果誤差error小於當前最小誤差：設置wbest等於當前的w
        將當前w設置為新的wbest

（1）需求

　　使用前向逐步回歸算法，算出最優回歸參數，並可視化顯示回歸參數隨迭代次數的變化情況，最后標准回歸算法進行比較

（2）代碼實現

def stage_wise(dataset,label,step,numiter=100):
    """
    前向逐步回歸
    """
    xmat = np.mat(dataset)
    xmean,xvar = data_normalize(dataset) #平均值和標准差
    ymat = np.mat(label).T
    ymat = ymat - np.mean(ymat,0)
    xmat = (xmat - xmean) / xvar
    m,n = np.shape(xmat)
    return_mat = np.zeros((numiter,n))
    ws = np.zeros((n,1))
    ws_test = ws.copy()
    ws_max = ws.copy()
    
    for iter in range(numiter):  #迭代次數
        lowest_error = np.inf
        for j in range(n):  #對每個特征進行變量
            for sign in [-1,1]:
                ws_test = ws.copy()
                ws_test[j] += step*sign
                ytest = xmat * ws_test
                rssE = ((ymat.A-ytest.A)**2).sum()
                if rssE < lowest_error:
                    lowest_error = rssE
                    ws_max = ws_test
        ws = ws_max.copy()
        return_mat[iter,:] = ws.T
    return return_mat    

def stage_wise_test(ws_mat,iternum,ax):
    """
    前向逐步回歸算法中回歸系數隨迭代次數的變化情況
    """    
    x = range(iternum)
    for i in range(np.shape(ws_mat)[1]):
        ax.plot(x,ws_mat[:,i])

（3）測試，輸出可視化效果

#前向逐步回歸
ax = fig.add_subplot(326)
ws = stage_wise(dataset,label,0.01,500)
print("前向逐步回歸經過迭代后的最后回歸系數：",ws[-1,:])
stage_wise_test(ws,500,ax)

　　效果：

　　上圖表示了回歸系數隨迭代次數的變化情況，當迭代150次時，回歸系數穩定分別為：

（4）接下來使用標准回歸算法，也就是最小二乘法對dataset/abalone.txt數據集進行測試，計算回歸系數

    #使用標准回歸對"dataset/abalone.txt"數據集進行測試
    #對數據進行標准化處理
    xman,xvar = data_normalize(dataset)
    xmat = np.mat(dataset)
    xmat = (xmat - xman) / xvar
    ymat = np.mat(label)
    ax = stand_regres(xmat,ymat)
    print("標准回歸后得到的回歸系數：",ws[-1,:])

　　最終輸出結果和前向逐步回歸算出系數一樣：

5. 總結

　　與分類一樣，回歸也是預測目標值的過程，回歸與分類的不同點在於，前者預測連續性變量，而后者預測離散型變量，回歸是統計學中最有力的工具之一，在回歸方程里，求得特征對應的最佳回歸系數的方法是最小 化誤差的平方和，給定輸入的矩陣X，如果X.T * X的逆存在並可以求得的化，回歸法可以直接使用，數據集上計算出的回歸方程並不一定意味着它是最佳的，可以使用預測值yHat和原始值y的相關性來度量回歸方程的好壞。　

　　當數據的樣本數比特征數還少的時候，矩陣X.T * X的逆不能直接計算，即便當樣本數比特征數多時，X.T * X的逆仍可能無法直接計算，這是因為特征有可能高度相關，這時可以考慮使用嶺回歸，因為當X.T * X的逆不能計算時，它仍保證能求得回歸系數

　　嶺回歸是縮減法的一種，相當於對回歸系數的大小施加了限制，另一種很好的縮減法是lasso,lasso難以求解，但可以 使用計算簡便的逐步線性回歸方法來求得近似解

　　最后來張綜合圖：