用Python代碼來理解LinearRegression

回歸的目的是預測數值型特徵的值,根據已有的數據,找到一個最佳的預測模式,這個最佳的模式就是回歸方程；模型的訓練，就是求解回歸方程的回歸係數w的過程。

那麼怎麼找到那組最合適的回歸係數w呢?在這之前，首先得回答什麼樣的w才是最好，標準是什麼。標準就是預測的y值和真實的y值之間的平方誤差和，之所以採用平方誤差，而不是直接將每個點的誤差相加，是因為簡單累加將使得正誤差和負誤差值相互抵消。

loss =

將上面的式子可以寫成向量的內積的形式:

loss=(y-wx).Transpose*(y-wx)

對w求導,並令求導後的式子等於零,可以求得w:

w=(x.Transpose * X).Inverse * X.Transpose * y

基礎的公式就這麼多,後面的都用代碼來說話.

一、標準的線性回歸

所謂標準的線性的回歸，就是按照上面列出的求回歸係數ｗ的公式，然後乘以一個ｘ，就能得出預測值yHat了。先寫一個從文本文件裡面load數據的函數，本文所用的函數的詳細理解主要寫在代碼注釋裡面。

def loadDataSet(fileName): numFeat = len(open(fileName).readline().split( ))-1 # 特徵數 dataMat = []; labelMat = [] # 分別存儲特徵和標籤 fr = open(fileName) for line in fr.readlines(): linArr=[] currLine = line.strip().split( ) # str去除首位空格和縮進 currLine=map(float, currLine) # map將str格式轉換成float數值型 dataMat.append(currLine[:-1]) labelMat.append(currLine[-1]) return dataMat, labelMat

下面就是標準線性回歸的代碼：

def standRegression(xArr, yArr): xMat = np.matrix(xArr) # 數據轉換為numpy的matrix yMat = np.matrix(yArr).T # 記得將y的矩陣轉置 xTx = xMat.T*xMat if np.linalg.det(xTx)==0.0: print 訓練集為奇異矩陣，不能求逆！！ return None ws = xTx.I*xMat.T*yMat # 先求逆，再乘以x的轉置，再乘以y return ws

下面的代碼顯示標準線性回歸的擬合效果：

X, Y = loadDataSet(ex0.txt) X = np.array(X) ws = standRegression(X,Y) # 為了方便在二維圖像上面表達，數據總共2個特徵，且第一個特徵都為1 index = X.argsort(axis=0)[:,1] # 將數據按照第二個特徵從小到大排序，返回的對應位置的索引 yHat = X[index]*ws # 用排序後的X和回歸係數，得到預測的y print 回歸係數ws：
, ws plt.scatter(X[:,1],Y,s=8) # 散點圖，點的size為8 w0 = ws[0].tolist()[0][0] # 從二維矩陣的ws中提出回歸係數，方便在圖上顯示，不用matrix格式 w1 = ws[1].tolist()[0][0] plt.plot(X[index][:,1],yHat,c=r,label=yHat=%.3f+%.3f*x % (w0,w1)) plt.legend() #　顯示圖例 plt.xlabel(x[:,1]) # x坐標 plt.ylabel(y) # y坐標 plt.show()

回歸係數:

可以看到，結果為是一條直線，然後使得數據點均勻地分布在直線的兩側，regression的詞根為regress，就是後悔的意思,在一些離直線很遠的點之後,往往就會出現往直線靠攏的趨勢,和浪子回頭差不多一個意思.

從寫的線性回歸的函數中可以發現，用了一個if語句判斷(X.T * X)的行列式是否等於零，如果等於零，表示該矩陣不可逆，那麼計算就終止了，所以用這種方法求回歸係數，要求(X.T * X)必須為滿秩矩陣，條件比較苛刻。

二、局部加權線性回歸(Locally Weighted Linear Regression)

標準的線性分布在全局進行優化，可以認為預測值這個隨機變數的期望與期望值相等的，不存在系統偏差，因而求得的是具有最小方差的無偏估計，如果模型欠擬合，預測效果將不會太好。而局部加權線性回歸，則允許在估計中引入一些偏差，從而降低預測的均方誤差。類似於標準的線性回歸，局部加權線性回歸計算權重也是要考慮全局訓練的訓練樣本，但是在預測某個點的值時，會通過核函數為訓練集的樣本賦予了一個權重weigths。

ｗ＝(Ｘ.T *weights* X).Inverse * X.T * weights* y

高斯核函數如下,類比於正態分布，|xi-x|表示樣本點到預測點的絕對距離，ｋ可以理解為正態分布的sigma，值越小，則高斯核函數的曲線就越瘦高，那麼離預測點距離越近樣本點的權重就越大，計算回歸係數時就傾向於只考慮預測點附近的幾個樣本點的規律；如果ｋ的值越大，核函數的曲線就越矮胖，越傾向於為全局的樣本點分配權重，ｋ特別大的時候，其實就和標準的線性回歸一樣了：

weights = exp(|xi - x|/(-2k**2))

局部加權線性回歸的原理很簡單，看看下面的代碼和注釋就會理解得非常清楚。

局部加權線性回歸，通過高斯核函數，根據其他樣本點與當前計算樣本點的距離，賦予每個樣本點不同的權重，可以更精確擬合局部信息，但是sigma太小，只考慮當前計算樣本附近很少的點，容易納入雜訊數據，發生過擬合。sigma如果很大，就會為大部分樣本點賦予差不多的權重，都考慮進來，類似於普通的線性回歸，容易欠擬合。def lwlr(testPoint, xArr, yArr, k): xMat = np.matrix(xArr);yMat=np.matrix(yArr).T m = np.shape(xMat)[0] weights = np.matrix(np.eye((m))) for j in range(m): diffMat = testPoint-xMat[j,:] disTance = diffMat*diffMat.T weights[j,j] = np.exp(disTance/(-2.0*k**2)) # 高斯核函數計算出權重 if np.linalg.det(xMat.T*weights*xMat)==0: print 訓練集是奇異矩陣哦～,不可逆 return None xTX = xMat.T*weights*xMat # 局部回歸演算法每個計算的樣本點，都會根據核函數計算其他樣本點對應的權重矩陣，矩陣大部分都為零。 ws = xTX.I*xMat.T*(weights*yMat)# 不同K之下的局部加權線性回歸X = np.matrix(X)Y = np.matrix(Y)kList = [0.1, 0.05, 0.003]yHatList = [] # 不同k值下的預測y值xSortList = [] # 因為需要畫曲線圖，所以將x從小到大排序for k in kList: yHat = lwlrTest(X, X, Y, k=k) sortInd = X[:, 1].argsort(0) # 按從小到大排列的索引 xSort = X[sortInd][:][:, 1] xSortList.append(xSort) yHat = yHat[sortInd] yHatList.append(yHat)plt.plot(xSortList[0], yHatList[0], c=y, label=k=0.1)plt.plot(xSortList[1], yHatList[1], c="m", label=k=0.05)plt.plot(xSortList[2], yHatList[2], c=b, label=k=0.003)plt.xlabel(u我是x軸)plt.ylabel(u我是y軸)plt.scatter(X[:, 1], Y, c=r, s=8)plt.legend()plt.show()

運行的結果為：

當K=0.1時，高斯核函數比較『矮胖』，這時候高斯核函數為大部分的訓練集樣本點都分配了權重，且權重相差不大，這就類似於標準線性回歸，可以看到回歸曲線近似為一條直線；當k=0.003時，高斯核函數比較瘦高，為預測點附近的訓練樣本賦予很高的權重，而離預測點距離較遠的點，隨著樣本點和預測點的距離增大，權重將以指數衰減，所以大部分樣本的權重幾乎為零，但是容易也給雜訊數據也分配很大比列的權重，過分表達，過擬合。

三、Ridge回歸和Shrinkage

前面兩種回歸演算法的代碼中，都有用if語句判斷X.T * X是否為可逆矩陣，但是很多情況下，往往都是不可逆的，比如有時候樣本數比特徵數還要少的時候。這時帶Shrinkage功能演算法就派上用場了。

Ridge回歸在X.T * X上加一個λ*I，是的矩陣可逆，這樣回歸係數的計算公式就變成：

下面用代碼來說明

def ridgeRegression(xArr, yArr,lam=0.5): xMat = np.matrix(xArr) yMat = np.matrix(yArr).T featNum = xMat.shape[1] eye = np.eye(featNum) # 創建單位矩陣 ws = (xMat.T*xMat+lam*eye).I*xMat.T*yMat return ws# ridge測試，用一個矩陣存貯不同lambda下的係數shrink情況.def ridgeTest(xArr, yArr): xMat = np.mat(xArr);yMat = np.mat(yArr) xMean =np.mean(xMat, axis=0) xMat = (xMat-xMean)/np.var(xMat,axis=0) # 首先需要標準化各個特徵。減去均值除以方差，標準正態分布 yMean = np.mean(yMat) yMat = (yMat-yMean) numTest = 30 # 30個lambda值 wMat = np.zeros((numTest, xMat.shape[1])) for i in range(numTest): ws = ridgeRegression(xMat, yMat,lam=np.exp(i-10)) # lambda值指數級變化 wMat[i,:] = ws.T return wMat # 返回一個矩陣，每一行為一個對應lambda值下的回歸係數abX, abY = loadDataSet(abalone.txt)ridgeWeights = ridgeTest(abX,abY)print ridgeWeightsfor i in range(ridgeWeights.shape[1]): plt.plot(ridgeWeights[:,i], label=featureNum_%s % i) plt.legend()plt.xlabel(lamb=exp(x-10))plt.ylabel(weights of features)plt.show()

運行結果為隨著λ不斷增加，不同特徵的的回歸係數的縮減趨勢：

λ很小時，係數與普通的線性回歸一樣，當λ很大時，8個特徵的回歸係數都傾向於縮減到零，可以在中間某處找到一個最佳λ值。ridge回歸會對特徵做選擇，其實是通過引入了一個L2範式來約束所有的回歸係數的加和。λ就是調整約束的權重，就是懲罰的力度，因而在高維特徵的數據，ridge回歸具有的shrink特性非常有用，也可以用於特徵提取。

四、總結

標準的線性回歸，考慮全局使得誤差平方和最小，計算出回歸係數，具有很低的方差，但是預測的偏差較大；局部加權線性回歸，用高斯核函數根據樣本點與預測點的遠近，賦予不同樣本不同的權重，可以通過調整sigma的值，來平衡模型的過度表達和泛化能力之間的trade-off；ridge回歸通過引入正則項，放棄了無偏性，以損失部分信息為代價，得到更好的精度。