Python3《機器學習實戰》學習筆記（十一）：線性回歸基礎篇之預測鮑魚年齡

轉載聲明：

本文轉自我的個人博客：機器學習實戰教程（十一）：線性回歸基礎篇之預測鮑魚年齡

文章首發地址：Jack Cui | 關注人工智慧及互聯網的個人博客

一、前言

前面的文章介紹了很多分類演算法，分類的目標變數是標稱型數據，而本文將會對連續型的數據做出預測。主要講解簡單的線性回歸和局部加權線性回歸，並通過預測鮑魚年齡的實例進行實戰演練。

二、什麼是回歸？

回歸的目的是預測數值型的目標值。最直接的辦法是依據輸入寫出一個目標值的計算公式。假如你想預測小姐姐男友汽車的功率，可能會這麼計算：

HorsePower = 0.0015 * annualSalary - 0.99 * hoursListeningToPublicRadio

寫成中文就是：

小姐姐男友汽車的功率 = 0.0015 * 小姐姐男友年薪 - 0.99 * 收聽公共廣播的時間

這就是所謂的回歸方程（regression equation），其中的0.0015和-0.99稱為回歸係數（regression weights），求這些回歸係數的過程就是回歸。一旦有了這些回歸係數，再給定輸入，做預測就非常容易了。具體的做法是用回歸係數乘以輸入值，再將結果全部加在一起，就得到了預測值。

說到回歸，一般都是指線性回歸（linear regression），所以本文里的回歸和線性回歸代表同一個意思。線性回歸意味著可以將輸入項分別乘以一些常量，再將結果加起來得到輸出。需要說明的是，存在另一種成為非線性回歸的回歸模型，該模型不認同上面的做法，比如認為輸出可能是輸入的乘積。這樣，上面的功率計算公式也可以寫做：

HorsePower = 0.0015 * annualSalary / hoursListeningToPublicRadio

這就是一個非線性回歸的例子，本文對此不做深入討論。

三、揭開回歸的神秘面紗

1、用線性回歸找到最佳擬合直線

應該怎麼從一大堆數據里求出回歸方程呢？假定輸入數據存放在矩陣X中，結果存放在向量y中：

而回歸係數存放在向量w中：

那麼對於給定的數據x1，即矩陣X的第一列數據，預測結果u1將會通過如下公式給出：

現在的問題是，手裡有數據矩陣X和對應的標籤向量y，怎麼才能找到w呢？一個常用的方法就是找出使誤差最小的w。這裡的誤差是指預測u值和真實y值之間的差值，使用該誤差的簡單累加將使得正差值和負差值相互抵消，所以我們採用平方誤差。

平方誤差和可以寫做：

用矩陣表示還可以寫做：

為啥能這麼變化，記住一個前提：若x為向量，則默認x為列向量，x^T為行向量。將上述提到的數據矩陣X和標籤向量y帶進去，就知道為何這麼變化了。

在繼續推導之前，我們要先明確一個目的：找到w，使平方誤差和最小。因為我們認為平方誤差和越小，說明線性回歸擬合效果越好。

現在，我們用矩陣表示的平方誤差和對w進行求導：

如果對於矩陣求不熟悉的，可以移步這裡：點擊查看

令上述公式等於0，得到：

w上方的小標記表示，這是當前可以估計出的w的最優解。從現有數據上估計出的w可能並不是數據中的真實w值，所以這裡使用了一個"帽"符號來表示它僅是w的一個最佳估計。

值得注意的是，上述公式中包含逆矩陣，也就是說，這個方程只在逆矩陣存在的時候使用，也即是這個矩陣是一個方陣，並且其行列式不為0。

述的最佳w求解是統計學中的常見問題，除了矩陣方法外還有很多其他方法可以解決。通過調用NumPy庫里的矩陣方法，我們可以僅使用幾行代碼就完成所需功能。該方法也稱作OLS， 意思是「普通小二乘法」（ordinary least squares）。

我們先看下數據集，數據下載地址：數據集下載

第一列都為1.0，即x0。第二列為x1，即x軸數據。第三列為x2，即y軸數據。首先繪製下數據，看下數據分布。編寫代碼如下：

# -*- coding:utf-8 -*-nimport matplotlib.pyplot as pltnimport numpy as npnndef loadDataSet(fileName):n """n 函數說明:載入數據n Parameters:n fileName - 文件名n Returns:n xArr - x數據集n yArr - y數據集n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """nn numFeat = len(open(fileName).readline().split(t)) - 1n xArr = []; yArr = []n fr = open(fileName)n for line in fr.readlines():n lineArr =[]n curLine = line.strip().split(t)n for i in range(numFeat):n lineArr.append(float(curLine[i]))n xArr.append(lineArr)n yArr.append(float(curLine[-1]))n return xArr, yArrnndef plotDataSet():n """n 函數說明:繪製數據集n Parameters:n 無n Returns:n 無n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n xArr, yArr = loadDataSet(ex0.txt) #載入數據集n n = len(xArr) #數據個數n xcord = []; ycord = [] #樣本點n for i in range(n): n xcord.append(xArr[i][1]); ycord.append(yArr[i]) #樣本點n fig = plt.figure()n ax = fig.add_subplot(111) #添加subplotn ax.scatter(xcord, ycord, s = 20, c = blue,alpha = .5) #繪製樣本點n plt.title(DataSet) #繪製titlen plt.xlabel(X)n plt.show()nnif __name__ == __main__:n plotDataSet()n

運行代碼如下：

通過可視化數據，我們可以看到數據的分布情況。接下來，讓我們根據上文中推導的回歸係數計算方法，求出回歸係數向量，並根據回歸係數向量繪製回歸曲線，編寫代碼如下：

# -*- coding:utf-8 -*-nimport matplotlib.pyplot as pltnimport numpy as npnndef loadDataSet(fileName):n """n 函數說明:載入數據n Parameters:n fileName - 文件名n Returns:n xArr - x數據集n yArr - y數據集n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n numFeat = len(open(fileName).readline().split(t)) - 1n xArr = []; yArr = []n fr = open(fileName)n for line in fr.readlines():n lineArr =[]n curLine = line.strip().split(t)n for i in range(numFeat):n lineArr.append(float(curLine[i]))n xArr.append(lineArr)n yArr.append(float(curLine[-1]))n return xArr, yArrnndef standRegres(xArr,yArr):n """n 函數說明:計算回歸係數wn Parameters:n xArr - x數據集n yArr - y數據集n Returns:n ws - 回歸係數n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n xMat = np.mat(xArr); yMat = np.mat(yArr).Tn xTx = xMat.T * xMat #根據文中推導的公示計算回歸係數n if np.linalg.det(xTx) == 0.0:n print("矩陣為奇異矩陣,不能求逆")n returnn ws = xTx.I * (xMat.T*yMat)n return wsnndef plotRegression():n """n 函數說明:繪製回歸曲線和數據點n Parameters:n 無n Returns:n 無n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n xArr, yArr = loadDataSet(ex0.txt) #載入數據集n ws = standRegres(xArr, yArr) #計算回歸係數n xMat = np.mat(xArr) #創建xMat矩陣n yMat = np.mat(yArr) #創建yMat矩陣n xCopy = xMat.copy() #深拷貝xMat矩陣n xCopy.sort(0) #排序n yHat = xCopy * ws #計算對應的y值n fig = plt.figure()n ax = fig.add_subplot(111) #添加subplotn ax.plot(xCopy[:, 1], yHat, c = red) #繪製回歸曲線n ax.scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = blue,alpha = .5) #繪製樣本點n plt.title(DataSet) #繪製titlen plt.xlabel(X)n plt.show()nnif __name__ == __main__:n plotRegression()n

運行代碼如下：

如何判斷擬合曲線的擬合效果的如何呢？當然，我們可以根據自己的經驗進行觀察，除此之外，我們還可以使用corrcoef方法，來比較預測值和真實值的相關性。編寫代碼如下：

# -*- coding:utf-8 -*-nimport numpy as npnndef loadDataSet(fileName):n """n 函數說明:載入數據n Parameters:n fileName - 文件名n Returns:n xArr - x數據集n yArr - y數據集n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n numFeat = len(open(fileName).readline().split(t)) - 1n xArr = []; yArr = []n fr = open(fileName)n for line in fr.readlines():n lineArr =[]n curLine = line.strip().split(t)n for i in range(numFeat):n lineArr.append(float(curLine[i]))n xArr.append(lineArr)n yArr.append(float(curLine[-1]))n return xArr, yArrnndef standRegres(xArr,yArr):n """n 函數說明:計算回歸係數wn Parameters:n xArr - x數據集n yArr - y數據集n Returns:n ws - 回歸係數n Website:n http://www.cuijiahua.com/n Modify:n 2017-11-12n """n xMat = np.mat(xArr); yMat = np.mat(yArr).Tn xTx = xMat.T * xMat #根據文中推導的公示計算回歸係數n if np.linalg.det(xTx) == 0.0:n print("矩陣為奇異矩陣,不能求逆")n returnn ws = xTx.I * (xMat.T*yMat)n return wsnnif __name__ == __main__:n xArr, yArr = loadDataSet(ex0.txt) #載入數據集n ws = standRegres(xArr, yArr) #計算回歸係數n xMat = np.mat(xArr) #創建xMat矩陣n yMat = np.mat(yArr) #創建yMat矩陣n yHat = xMat * wsn print(np.corrcoef(yHat.T, yMat))n

運行結果如下：

可以看到，對角線上的數據是1.0，因為yMat和自己的匹配是完美的，而YHat和yMat的相關係數為0.98。

最佳擬合直線方法將數據視為直線進行建模，具有十分不錯的表現。數據當中似乎還存在其他的潛在模式。那麼如何才能利用這些模式呢？我們可以根據數據來局部調整預測，下面就會介紹這種方法。

2、局部加權線性回歸

線性回歸的一個問題是有可能出現欠擬合現象，因為它求的是具有小均方誤差的無偏估 計。顯而易見，如果模型欠擬合將不能取得好的預測效果。所以有些方法允許在估計中引入一 些偏差，從而降低預測的均方誤差。

其中的一個方法是局部加權線性回歸（Locally Weighted Linear Regression，LWLR）。在該方法中，我們給待預測點附近的每個點賦予一定的權重。與kNN一樣，這種演算法每次預測均需要事先選取出對應的數據子集。該演算法解除回歸係數W的形式如下：

其中W是一個矩陣，這個公式跟我們上面推導的公式的區別就在於W，它用來給每個店賦予權重。

LWLR使用"核"（與支持向量機中的核類似）來對附近的點賦予更高的權重。核的類型可以自由選擇，最常用的核就是高斯核，高斯核對應的權重如下：

這樣我們就可以根據上述公式，編寫局部加權線性回歸，我們通過改變k的值，可以調節回歸效果，編寫代碼如下：

# -*- coding:utf-8 -*-nfrom matplotlib.font_manager import FontPropertiesnimport matplotlib.pyplot as pltnimport numpy as npndef loadDataSet(fileName):n """n 函數說明:載入數據n Parameters:n fileName - 文件名n Returns:n xArr - x數據集n yArr - y數據集n Website:http://www.cuijiahua.com/n Modify:n 2017-11-12n """n numFeat = len(open(fileName).readline().split(t)) - 1n xArr = []; yArr = []n fr = open(fileName)n for line in fr.readlines():n lineArr =[]n curLine = line.strip().split(t)n for i in range(numFeat):n lineArr.append(float(curLine[i]))n xArr.append(lineArr)n yArr.append(float(curLine[-1]))n return xArr, yArrnndef plotlwlrRegression():n """n 函數說明:繪製多條局部加權回歸曲線n Parameters:n 無n Returns:n 無n Website:http://www.cuijiahua.com/n Modify:n 2017-11-15n """n font = FontProperties(fname=r"c:windowsfontssimsun.ttc", size=14)n xArr, yArr = loadDataSet(ex0.txt) #載入數據集n yHat_1 = lwlrTest(xArr, xArr, yArr, 1.0) #根據局部加權線性回歸計算yHatn yHat_2 = lwlrTest(xArr, xArr, yArr, 0.01) #根據局部加權線性回歸計算yHatn yHat_3 = lwlrTest(xArr, xArr, yArr, 0.003) #根據局部加權線性回歸計算yHatn xMat = np.mat(xArr) #創建xMat矩陣n yMat = np.mat(yArr) #創建yMat矩陣n srtInd = xMat[:, 1].argsort(0) #排序，返回索引值n xSort = xMat[srtInd][:,0,:]n fig, axs = plt.subplots(nrows=3, ncols=1,sharex=False, sharey=False, figsize=(10,8)) n axs[0].plot(xSort[:, 1], yHat_1[srtInd], c = red) #繪製回歸曲線n axs[1].plot(xSort[:, 1], yHat_2[srtInd], c = red) #繪製回歸曲線n axs[2].plot(xSort[:, 1], yHat_3[srtInd], c = red) #繪製回歸曲線n axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = blue, alpha = .5) #繪製樣本點n axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = blue, alpha = .5) #繪製樣本點n axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = blue, alpha = .5) #繪製樣本點n #設置標題,x軸label,y軸labeln axs0_title_text = axs[0].set_title(u局部加權回歸曲線,k=1.0,FontProperties=font)n axs1_title_text = axs[1].set_title(u局部加權回歸曲線,k=0.01,FontProperties=font)n axs2_title_text = axs[2].set_title(u局部加權回歸曲線,k=0.003,FontProperties=font)n plt.setp(axs0_title_text, size=8, weight=bold, color=red) n plt.setp(axs1_title_text, size=8, weight=bold, color=red) n plt.setp(axs2_title_text, size=8, weight=bold, color=red) n plt.xlabel(X)n plt.show()ndef lwlr(testPoint, xArr, yArr, k = 1.0):n """n 函數說明:使用局部加權線性回歸計算回歸係數wn Parameters:n testPoint - 測試樣本點n xArr - x數據集n yArr - y數據集n k - 高斯核的k,自定義參數n Returns:n ws - 回歸係數n Website:http://www.cuijiahua.com/n Modify:n 2017-11-15n """n xMat = np.mat(xArr); yMat = np.mat(yArr).Tn m = np.shape(xMat)[0]n weights = np.mat(np.eye((m))) #創建權重對角矩陣n for j in range(m): #遍曆數據集計算每個樣本的權重n diffMat = testPoint - xMat[j, :] n weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))n xTx = xMat.T * (weights * xMat) n if np.linalg.det(xTx) == 0.0:n print("矩陣為奇異矩陣,不能求逆")n returnn ws = xTx.I * (xMat.T * (weights * yMat)) #計算回歸係數n return testPoint * wsndef lwlrTest(testArr, xArr, yArr, k=1.0): n """n 函數說明:局部加權線性回歸測試n Parameters:n testArr - 測試數據集n xArr - x數據集n yArr - y數據集n k - 高斯核的k,自定義參數n Returns:n ws - 回歸係數n Website:nhttp://www.cuijiahua.com/n Modify:n 2017-11-15n """n m = np.shape(testArr)[0] #計算測試數據集大小n yHat = np.zeros(m) n for i in range(m): #對每個樣本點進行預測n yHat[i] = lwlr(testArr[i],xArr,yArr,k)n return yHatnif __name__ == __main__:n plotlwlrRegression()n

運行結果如下：

可以看到，當k越小，擬合效果越好。但是當k過小，會出現過擬合的情況，例如k等於0.003的時候。

四、預測鮑魚的年齡

接下來，我們將回歸用於真是數據。在abalone.txt文件中記錄了鮑魚（一種水生物→__→）的年齡，這個數據來自UCI數據集合的數據。鮑魚年齡可以從鮑魚殼的層數推算得到。

數據集下載地址：數據集下載

數據集的數據如下：

可以看到，數據集是多維的，所以我們很難畫出它的分布情況。每個維度數據的代表的含義沒有給出，不過沒有關係，我們只要知道最後一列的數據是y值就可以了，最後一列代表的是鮑魚的真實年齡，前面幾列的數據是一些鮑魚的特徵，例如鮑魚殼的層數等。我們不做數據清理，直接用上所有特徵，測試下我們的局部加權回歸。

新建abalone.py文件，添加rssError函數，用於評價最後回歸結果。編寫代碼如下：

# -*- coding:utf-8 -*-nfrom matplotlib.font_manager import FontPropertiesnimport matplotlib.pyplot as pltnimport numpy as npndef loadDataSet(fileName):n """n 函數說明:載入數據n Parameters:n fileName - 文件名n Returns:n xArr - x數據集n yArr - y數據集n Website:nhttp://www.cuijiahua.com/n Modify:n 2017-11-19n """n numFeat = len(open(fileName).readline().split(t)) - 1n xArr = []; yArr = []n fr = open(fileName)n for line in fr.readlines():n lineArr =[]n curLine = line.strip().split(t)n for i in range(numFeat):n lineArr.append(float(curLine[i]))n xArr.append(lineArr)n yArr.append(float(curLine[-1]))n return xArr, yArrndef lwlr(testPoint, xArr, yArr, k = 1.0):n """n 函數說明:使用局部加權線性回歸計算回歸係數wn Parameters:n testPoint - 測試樣本點n xArr - x數據集n yArr - y數據集n k - 高斯核的k,自定義參數n Returns:n ws - 回歸係數n Website:nhttp://www.cuijiahua.com/n Modify:n 2017-11-19n """n xMat = np.mat(xArr); yMat = np.mat(yArr).Tn m = np.shape(xMat)[0]n weights = np.mat(np.eye((m))) #創建權重對角矩陣n for j in range(m): #遍曆數據集計算每個樣本的權重n diffMat = testPoint - xMat[j, :] n weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))n xTx = xMat.T * (weights * xMat) n if np.linalg.det(xTx) == 0.0:n print("矩陣為奇異矩陣,不能求逆")n returnn ws = xTx.I * (xMat.T * (weights * yMat)) #計算回歸係數n return testPoint * wsndef lwlrTest(testArr, xArr, yArr, k=1.0): n """n 函數說明:局部加權線性回歸測試n Parameters:n testArr - 測試數據集,測試集n xArr - x數據集,訓練集n yArr - y數據集,訓練集n k - 高斯核的k,自定義參數n Returns:n ws - 回歸係數n Website:nhttp://www.cuijiahua.com/n Modify:n 2017-11-19n """n m = np.shape(testArr)[0] #計算測試數據集大小n yHat = np.zeros(m) n for i in range(m): #對每個樣本點進行預測n yHat[i] = lwlr(testArr[i],xArr,yArr,k)n return yHatndef standRegres(xArr,yArr):n """n 函數說明:計算回歸係數wn Parameters:n xArr - x數據集n yArr - y數據集n Returns:n ws - 回歸係數n Website:nhttp://www.cuijiahua.com/n Modify:n 2017-11-19n """n xMat = np.mat(xArr); yMat = np.mat(yArr).Tn xTx = xMat.T * xMat #根據文中推導的公示計算回歸係數n if np.linalg.det(xTx) == 0.0:n print("矩陣為奇異矩陣,不能求逆")n returnn ws = xTx.I * (xMat.T*yMat)n return wsndef rssError(yArr, yHatArr):n """n 誤差大小評價函數n Parameters:n yArr - 真實數據n yHatArr - 預測數據n Returns:n 誤差大小n """n return ((yArr - yHatArr) **2).sum()nif __name__ == __main__:n abX, abY = loadDataSet(abalone.txt)n print(訓練集與測試集相同:局部加權線性回歸,核k的大小對預測的影響:)n yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)n yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)n yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)n print(k=0.1時,誤差大小為:,rssError(abY[0:99], yHat01.T))n print(k=1 時,誤差大小為:,rssError(abY[0:99], yHat1.T))n print(k=10 時,誤差大小為:,rssError(abY[0:99], yHat10.T))n print()n print(訓練集與測試集不同:局部加權線性回歸,核k的大小是越小越好嗎？更換數據集,測試結果如下:)n yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)n yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)n yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)n print(k=0.1時,誤差大小為:,rssError(abY[100:199], yHat01.T))n print(k=1 時,誤差大小為:,rssError(abY[100:199], yHat1.T))n print(k=10 時,誤差大小為:,rssError(abY[100:199], yHat10.T))n print()n print(訓練集與測試集不同:簡單的線性歸回與k=1時的局部加權線性回歸對比:)n print(k=1時,誤差大小為:, rssError(abY[100:199], yHat1.T))n ws = standRegres(abX[0:99], abY[0:99])n yHat = np.mat(abX[100:199]) * wsn print(簡單的線性回歸誤差大小:, rssError(abY[100:199], yHat.T.A))n

運行結果如下：

可以看到，當k=0.1時，訓練集誤差小，但是應用於新的數據集之後，誤差反而變大了。這就是經常說道的過擬合現象。我們訓練的模型，我們要保證測試集準確率高，這樣訓練出的模型才可以應用於新的數據，也就是要加強模型的普適性。可以看到，當k=1時，局部加權線性回歸和簡單的線性回歸得到的效果差不多。這也表明一點，必須在未知數據上比較效果才能選取到最佳模型。那麼最佳的核大小是10嗎？或許是，但如果想得到更好的效果，應該用10個不同的樣本集做10次測試來比較結果。

本示例展示了如何使用局部加權線性回歸來構建模型，可以得到比普通線性回歸更好的效果。局部加權線性回歸的問題在於，每次必須在整個數據集上運行。也就是說為了做出預測，必須保存所有的訓練數據。

五、總結

• 本文主要介紹了簡單的線性回歸和局部加權線性回歸。
• 訓練的模型要在測試集比較它們的效果，而不是在訓練集上。
• 在局部加權線性回歸中，過小的核可能導致過擬合現象，即測試集表現良好，訓練集表現就渣渣了。
• 下篇文章將繼續講解回歸，會介紹另一種提高預測精度的方法。
• 如有問題，請留言。如有錯誤，還望指正，謝謝！

PS： 如果覺得本篇本章對您有所幫助，歡迎關注、評論、贊！

本文出現的所有代碼和數據集，均可在我的github上下載，歡迎Follow、Star：點擊查看

參考資料：

• [1] 機器學習實戰第八章內容

推薦閱讀：