[ ML In Action ] Predicting numeric values : regression - Linear regression (3)

Preface : 
假如我們的 Data matrix X(m,n) 中的 n(feature size) < m(data size), 則因為該矩陣不是 full rank, 所以在求反矩陣可能會出現問題. 這時便需要考慮是否有些 feature 是沒有 contribute 到 weighting vector 所以可以拿掉以降低 X 的 rank. 

為了解決這個問題, 這邊要介紹 "Ridge regression" : 

Tikhonov regularization, named for Andrey Tikhonov, is the most commonly used method of regularization of ill-posed problems. In statistics, the method is known as ridge regression.

後面要介紹的還有相關的 lasso 與其較簡易的實作方法 "stagewise regression". 


Ridge regression : 
在 Ridge regression 導入了 λI 到原先的公式中; 而 λ 是一個 user-defined scalar value, 原先 weighting 的公式改寫如下 : 
 

底下是書上對於 Ridge regression 的說明 : 

Ridge regression was originally developed to deal with the problem of having more features than data points. But it can also be used to add bias into our estimations, giving us a better estimate. We can use the λ value to impose a maximum value on the sum of all our ws. By imposing this penalty, we can decrease unimportant parameters. This decreasing is known as shrinkage in statistics.

所謂的 "Shrinkage methods" 是要幫助我們更了解 training 的 data, 將沒有用的 features (weighting=0) 從 training data 移除已得到較好的 prediction. 

實務上 λ 的取的是透過 iteration 不斷的計算 ws, 並將有最小的 prediction error 對應的 λ 取出作為輸出 training mode 使用的 λ. 底下方法 ridgeRegres() 用於計算某個特定 λ 的 ws ; 而函數ridgeTest() 則從 λ 的範圍從 e^-10 到 e^20 中求出對應的 ws 並存在變數 wMat 後回傳.: 
- regression.py : 

view plaincopy to clipboardprint?

1. def ridgeRegres(xMat,yMat,lam=0.2):  
2.     xTx = xMat.T*xMat  
3.     denom = xTx + eye(shape(xMat)[1])*lam  
4.     if linalg.det(denom) == 0.0:  
5.         print "This matrix is singular, cannot do inverse"  
6.         return  
7.     ws = denom.I * (xMat.T*yMat)  
8.     return ws  
9.   
10. def ridgeTest(xArr,yArr):  
11.     xMat = mat(xArr); yMat=mat(yArr).T  
12.     yMean = mean(yMat,0)  
13.     yMat = yMat - yMean     #to eliminate X0 take mean off of Y  
14.     #regularize X's  
15.     xMeans = mean(xMat,0)   #calc mean then subtract it off  
16.     xVar = var(xMat,0)      #calc variance of Xi then divide by it  
17.     xMat = (xMat - xMeans)/xVar  
18.     numTestPts = 30  
19.     wMat = zeros((numTestPts,shape(xMat)[1]))  
20.     for i in range(numTestPts):  
21.         ws = ridgeRegres(xMat,yMat,exp(i-10))  
22.         wMat[i,:]=ws.T  
23.     return wMat  

下面代碼建立 wMat 矩陣, numTestPts 說明計算的 ws 數目 ; shape(xMat)[1] 則回得到 features 的數目 : 

wMat = zeros((numTestPts,shape(xMat)[1]))

有了這些理解, 接著來看看如何使用, 首先下面載入 training set "abalone.txt" 並建立不同 λ 的 ws 所組成的矩陣 (ridgeWeights): 

>>> reload(regression)

>>> abX, abY = regression.loadDataSet('abalone.txt')
>>> ridgeWeights=regression.ridgeTest(abX, abY)
>>> ridgeWeights[0] # 列印 λ=-10 的 ws
array([ 0.04304419, -0.02274163, 0.13214088, 0.02075182, 2.22403745,
-0.99895298, -0.11725424, 0.16622922])

接著我們將剛剛建立的矩陣 "ridgeWeights" 用圖形視覺化 : 

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(ridgeWeights)
>>> plt.show()

執行後可以得到下面圖形 : 
 
(x 軸為 λ 的 log 值; y 軸為 ws 的值. 不同顏色的線代表不同的 ws 欄位) 

可以發現 ws 隨著 λ 值得增大而逐漸減小. 這也是為什麼說可以透過 λ 讓不重要的 features fade out. 而讓更重要的 features 留下. 

Forward stagewise regression : 
從上面的說明, 我們還是不知道怎麼求出最佳的 ws. 所謂最佳的 ws 就是它能有最小的 error (計算預測值於實際值的差, 通常會取 square root 來避免正負數抵銷). 這邊要介紹的 stagewise regression 演算法是一種 greedy algorithm, 在每一次計算 ws, 它都會嘗試求出當下最佳解. 而初始的 ws 的值會是 0. 該演算法的 Pseudo-code 如下 : 
 

在一開始我們先來定義計算 error 的函數 rssError() : 

view plaincopy to clipboardprint?

1. def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays  
2.     return ((yArr-yHatArr)**2).sum()  

接著來看看 stagewise regression 演算法的實作代碼 : 

view plaincopy to clipboardprint?

1. def stageWise(xArr,yArr,eps=0.01,numIt=100):  
2.     xMat = mat(xArr); yMat=mat(yArr).T  
3.     yMean = mean(yMat,0)  
4.     yMat = yMat - yMean     #can also regularize ys but will get smaller coef  
5.     xMat = regularize(xMat)  
6.     m,n=shape(xMat)  
7.     returnMat = zeros((numIt,n)) #testing code remove  
8.     ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()  
9.     for i in range(numIt): # For each iteration  
10.         print ws.T  
11.         lowestError = inf;  
12.         for j in range(n): # For each feature in ws  
13.             for sign in [-1,1]:  # For +eps or -eps  
14.                 wsTest = ws.copy()  
15.                 wsTest[j] += eps*sign  
16.                 yTest = xMat*wsTest  
17.                 rssE = rssError(yMat.A,yTest.A)  
18.                 if rssE < lowestError:  
19.                     lowestError = rssE  
20.                     wsMax = wsTest  
21.         ws = wsMax.copy()  
22.         returnMat[i,:]=ws.T  
23.     return returnMat  

函數 stageWise() 會返回矩陣 returnMat, 它包含每一次 iteration 計算出來當下最佳的 ws ; 你可以透過指定參數 numIt 決定 iteration 的次數 ; 而參數 eps 可以設定每次移動 feature 值的間距. 接著來看看如何使用 : 

>>> xArr, yArr = regression.loadDataSet('abalone.txt')
>>> regression.stageWise(xArr, yArr, 0.01, 200) # 計算 200 次 iteration, 使用間距為 0.01
array([[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
...,
[ 0.05, 0. , 0.09, ..., -0.64, 0. , 0.36],
[ 0.04, 0. , 0.09, ..., -0.64, 0. , 0.36],
[ 0.05, 0. , 0.09, ..., -0.64, 0. , 0.36]])

從結果你可以知道 w1 與 w6 為 0, 說明它們對 prediction 沒有幫助, 另外 w0 會在 0.04~0.05 間震盪, 可能是 eps 過大, 你可以試著調小 eps 試試看會不會收斂. 

接著使用下面命令計算的 stagewise matrix 

>>> stageWiseMat = regression.stageWise(xArr,yArr,0.001,5000)

將之視覺化如下圖, 可以發現 ws 的各個 feature 逐漸收斂 : 
 

這邊我們使用的 stagewise linear regression 或是 ridge regression 事實上都是加上了 bias 到原先的 model 中 ; 同時我們也在減少 model variance. 

The bias/variance tradeoff : 
在實際的應用中, 我們取得的 training set 可能會包含 error. 而這會造成我們 training 出來的 model 為了要有最佳(最小)的 training error 而提升了 model 的複雜度 (增加 features 或增加 bias). 如前面我們使用了 Locally weighted linear regression 來讓預測的線型更貼近 training 的 data 分布, 而這樣的最法有可能會造成 overfit ; 另一種作法則是如上面介紹的 shrinkage methods, 透過引入 bias λ 來減少 features (減少複雜度, 降低 variance), 但這可能會造成 underfit. 問題是我們並不知道實際上的 model 的長相! 所以在取捨該增加 features (variance, 原先的 feature 不足以描述實際的 model) 或是透過這邊介紹的 shrinkage methods 來減少不必要的 features (增加 bias, 降低 training model 的 variance). 下圖說明 bias/variance 取捨上的 tradeoff : 
 
(在 bottom 的曲線是 training error; 在 top 的曲線是 testing error) 

從上圖可以知道透過增加我們 training model 的複雜度 (High variance, low bias) 雖然可以減少 training error 但卻會遇到 overfit 的問題而提昇 testing error; 同樣的就算你降低 training model 的複雜度 (增加 bias/減少 variance) 但同樣會遇到 underfit 的問題而增加 testing error.

Supplement : 
* [ ML In Action ] Predicting numeric values : regression - Linear regression (1) 
* [ ML In Action ] Predicting numeric values : regression - Linear regression (2) 
* [ ML In Action ] Predicting numeric values : regression - Linear regression (3)