本文寫於2018.10.28，原文地址，技術博客地址

如果覺得本文對你有幫助，點贊讓更多的同學看到哦~

說明：造輪子是為了深入理解演算法原理和參數，實際工作中不用造輪子

01 分類 v.s. 回歸

之前我們學習了很多分類方法，在機器學習中，還有一種任務叫回歸，回歸和分類其實挺像的，都是對樣本預測一個值，區別在於， - 分類：輸出為離散值 - 回歸：輸出為連續值

今天我們學習一波線性回歸的理論和演算法，不要小看線性回歸，其實很多商業模型都少不了線性回歸的功勞，把線性回歸用到極致你也是大神。

簡單來說，線性回歸就是在已知x,y的情況下，求解y=wx的回歸係數的過程。

02 標準線性回歸

標準線性回歸就是最簡單的線性回歸，其原理是最小化預測值與真實值的均方誤差，即最小二乘法。根據原理，得到係數公式：

標準線性回歸簡單粗暴，易於求解和解釋，但其缺點也顯而易見：對非線性樣本，欠擬合。

實現代碼如下：

#標準線性回歸函數，假設數據有m個樣本,n個特徵
def standRegres(xArr,yArr):
xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=1*m
xTx=xMat.T*xMat #shape=n*n
if np.linalg.det(xTx)==0.0: #矩陣行列式若為0，說明該矩陣不可逆
print (xTx為奇異矩陣，不可逆)
return
ws=xTx.I*(xMat.T*yMat) #shape=n*1,n為特徵數
return ws

測試一波：

dataMat,labelMat=loadDataSet(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh08ex0.txt)
ws=standRegres(dataMat,labelMat)

print (回歸係數w為：,round(float(ws[0]),4),round(float(ws[1]),4))
print (回歸方程：y={1}x+{0}.format(round(float(ws[0]),4),round(float(ws[1]),4)))

輸出如下，

均方誤差為1.355左右，

用於測試的數據集較簡單，可以作圖看看線性回歸的預測效果：

在上圖中我們可以看到，預測值非常「直」，雖然捕捉到了真實值的大趨勢，但沒有捕捉到真實數據的潛在特徵(非線性)，下面我們改進一下標準回歸演算法，嘗試捕捉潛在特徵。

03 局部加權線性回歸

局部加權線性回歸認為附近的點有相似的回歸特性，因此對待預測點附近的點加權。

之前我寫利用拉勾網的數據寫了一個人才價格計算器(相關文章在這裡)，其原理是將分類演算法KNN用於回歸，即對樣本分類後將同一類別的樣本標籤設置為該類別的均值，這個人才價格計算器的原理其實和今天這個局部加權線性回歸很相似。（當然，缺點也很相似，笑~）

利用這樣的想法，局部加權線性回歸的回歸係數公式如下：

其中，W為權重矩陣，

由權重矩陣公式可知， 1. W是對角矩陣 2. 樣本點x距離預測點xi越近，權重越大 3. 對同一個樣本點，k越大,權重矩陣W越大,賦予給每個數據點的權重越大,當k減小，W也減小，賦予給各點的權重減小 4. 權重最先減小到0的是距離預測點較遠的樣本點

下面我們實現此演算法

#局部加權線性回歸函數，假設數據有m個樣本,n個特徵
#針對某個測試點(testPoint),計算一個ws
"""
因為針對每個測試點，都有一個weight賦權矩陣，因此每個測試點都有一個ws回歸係數矩陣，
因此對於每個測試點，都需要利用lwlr()得到一個測試點的預測值testPoint*ws
因此lwlr()返回testPoint*ws
"""
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat=np.mat(xArr);yMat=np.mat(yArr) #shape(yMat)=m*1
m=np.shape(xMat)[0]
weight=np.mat(np.eye(m)) #構造一個m*m的對角矩陣
#針對每個樣本點j，計算針對測試點i的weight矩陣Wj,j的值
for j in range(m):
diffMat=testPoint-xMat[j,:]
weight[j,j]=np.exp(diffMat*diffMat.T/(-2.0*k**2))
xTx=xMat.T*weight*xMat
if np.linalg.det(xTx)==0.0: #矩陣行列式若為0，說明該矩陣不可逆
print (xTx為奇異矩陣，不可逆)
return
ws=xTx.I*(xMat.T*weight*yMat)
return testPoint*ws

def lwlrTest(testArr,xArr,yArr,k=1.0):
m=np.shape(testArr)[0]
yPred=np.zeros(m)
for i in range(m):
yPred[i]=lwlr(testArr[i],xArr,yArr,k)
return yPred

測試一下

dataMat,labelMat=loadDataSet(rD:DMpythondataMLiA_SourceCodemachinelearninginactionCh08ex0.txt)
yPred1=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,1.0)
yPred2=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,0.03)
yPred3=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,0.003)
rssError(labelMat,yPred1),rssError(labelMat,yPred2),rssError(labelMat,yPred3)

可以看到，訓練集的均方誤差隨著k的減小而減小，但k過小可能會過擬合(相當於KNN的k取很小的情況，預測點至於其附近幾個點相關，就會過擬合)

缺點：類比KNN。對每個點做預測時，都必須用到整個數據集，計算量較大。

04 嶺回歸ridge

若特徵數比樣本數還多(n>m)，即特徵矩陣X非滿秩，那麼xTx不可逆，在用上述推導的公式求解w時會出問題，此時怎麼辦呢？

我們讓xTx加上一個值：xTx+lambda*I，通過加上lambda*I使矩陣非奇異，從而求逆矩陣(xTx+lambda*I)，這就是縮減的原理，通過引入lambda來限制回歸係數w之和，通過引入lambda懲罰項，減少不重要的參數。

當約束條件如下時，回歸就是嶺回歸：

嶺回歸的回歸係數公式如下：

其中， 1. lambda是一個常數，可以通過k折交叉驗證，尋找最小的均方誤差，從而找到最佳的lambda 2. I是一個n*n的單位對角矩陣，n為特徵數

實現一下，

#求嶺回歸ws函數
def ridgeRegres(xMat,yMat,lam=0.2):
xTx=xMat.T*xMat
denom=xTx+lam*np.eye(np.shape(xMat)[1]) #xTx+lambda*I
if np.linalg.det(denom)==0.0:
print (xTx為奇異矩陣，不可逆)
return
ws=denom.I*xMat.T*yMat
return ws

#z-score標準化函數
def regularize(xMat):
reMat=xMat.copy()
reMeans=np.mean(reMat,0) #按列求均值
reVar=np.var(reMat,0) #按列求標準差
reMat=(reMat-reMeans)/reVar
return reMat

#遍歷lambda，求各lambda對應的回歸係數ws函數
def ridgeTest(xArr,yArr):
xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=m*1
yMean=np.mean(yMat,0)
yMat=yMat-yMean
xMat=regularize(xMat) #對輸入數據標準化
numTestPts=30 #設定遍歷30個lambda
wMat=np.zeros((numTestPts,np.shape(xMat)[1]))
for i in range(numTestPts):
ws=ridgeRegres(xMat,yMat,np.exp(i-10)) #指數級別遍歷lambda
wMat[i,:]=ws.T #n*1>>1*n
return wMat

測試，我們用一組鮑魚的數據來預測鮑魚年齡，圖中8條曲線分別表示w0~w7在30個不同的lambda下的值。

從結果看，回歸係數有8個w0-w7，在不同lambda下會得到不同的係數，隨著lambda增大，到足夠大時(縮減足夠大)，各參數減小至接近0，趨於相等，這就是縮減，縮減特徵。

05 套索回歸lasso+前向逐步回歸

同上節，當約束條件如下時，回歸就是套索回歸：

可以看到，嶺回歸對平方函數求極值，可使用最小二乘法，但lasso回歸對絕對值求極值，求解較複雜，怎麼辦呢？

這時我們可以使用前向分步演算法逼近求解

前向分佈演算法是一種貪心演算法，每一步都對某個權重增加或減少一個很小的值，從而在每一步儘可能減小誤差。這樣做雖然循環會變多，但每一次循環的計算都變得十分簡單。

實現一下

def stageWise(xArr,yArr,step=0.01,numIt=100):
xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=m*1
yMean=np.mean(yMat,0)
yMat=yMat-yMean
xMat=regularize(xMat) #對輸入數據標準化
m,n=np.shape(xMat)
returnMat=np.zeros((numIt,n))
ws=np.zeros((n,1));wsBest=ws.copy
for i in range(numIt): #對每次迭代
minErr=np.inf #將最小誤差初始化為正無窮
for j in range(n): #對每個特徵
for sign in [-1,1]: #對每個加/減
wsTest=ws.copy()
wsTest[j]+=step*sign #對第j個係數進行加/減
yTest=xMat*wsTest
rss=rssError(yMat.A,yTest.A)
if rss<minErr:
minErr=rss
wsBest=wsTest
ws=wsBest.copy() #找到本次迭代最小誤差對應的係數矩陣
returnMat[i,:]=ws.T
return returnMat

測試，

可以看到，w0,w1,w4,w5,w6這四個係數在每次迭代中都保持在0附近，可以認為這幾個係數不重要，可以剔除這幾個維度，達到降維的作用(lasso回歸的特性之一)

06 總結

本文總結了幾種線性回歸的原理，並用python3實現了這些回歸。

常見的線性回歸包括標準線性回歸、局部加權線性回歸、嶺回歸Ridge、套索回歸Lasso、前向逐步回歸(lasso回歸的近似求解)。除了這些之外，還有彈性網路回歸ElasticNet Regression，它綜合了Ridge和Lasso的優點，懲罰項是L1範數和L2範數的融合。

下期我們再來看看樹回歸，這樣回歸板塊的大部分知識就掌握了，敬請期待~

07 參考

  • 《機器學習實戰》 Peter Harrington Chapter8

如果覺得本文對你有幫助，點贊讓更多的同學看到吧

推薦閱讀：

相關文章