---------------------------------------------------------------------------------------
本系列文章為《機器學習實戰》學習筆記，內容整理自書本，網絡以及自己的理解，如有錯誤歡迎指正。
源碼在Python3.5上測試均通過，代碼及數據 --> https://github.com/Wellat/MLaction
---------------------------------------------------------------------------------------
1、線性回歸
現有一數據集，其分布如下圖所示，
通過觀察發現可以通過一個線性方程去擬合這些數據點。可設直線方程為 y=wx. 其中w稱為回歸系數。那么現在的問題是，如何從一堆x和對應的y中確定w？一個常用的方法就是找出使誤差最小的w。這里的誤差是指預測y值和真實y值之間的差值，我們采用平方誤差，寫作：
用矩陣還可以寫作：
，如果對w求導，得到
，令其等於零，解出w為：
注意此處公式包含對矩陣求逆，所以求解時需要先對矩陣是否可逆做出判斷。以上求解w的過程也稱為“普通最小二乘法”。
Python實現代碼如下：
1 from numpy import * 2 3 def loadDataSet(fileName): 4 '''導入數據''' 5 numFeat = len(open(fileName).readline().split('\t')) - 1 6 dataMat = []; labelMat = [] 7 fr = open(fileName) 8 for line in fr.readlines(): 9 lineArr =[] 10 curLine = line.strip().split('\t') 11 for i in range(numFeat): 12 lineArr.append(float(curLine[i])) 13 dataMat.append(lineArr) 14 labelMat.append(float(curLine[-1])) 15 return dataMat,labelMat 16 17 def standRegres(xArr,yArr): 18 '''求回歸系數''' 19 xMat = mat(xArr); yMat = mat(yArr).T 20 xTx = xMat.T*xMat 21 if linalg.det(xTx) == 0.0:#判斷行列式是否為0 22 print("This matrix is singular, cannot do inverse") 23 return 24 ws = xTx.I * (xMat.T*yMat)#也可以用NumPy庫的函數求解：ws=linalg.solve(xTx,xMat.T*yMatT) 25 return ws 26 27 if __name__ == "__main__": 28 '''線性回歸''' 29 xArr,yArr=loadDataSet('ex0.txt') 30 ws=standRegres(xArr,yArr) 31 xMat=mat(xArr) 32 yMat=mat(yArr) 33 #預測值 34 yHat=xMat*ws 35 36 #計算預測值和真實值得相關性 37 corrcoef(yHat.T,yMat)#0.986 38 39 #繪制數據集散點圖和最佳擬合直線圖 40 #創建圖像並繪出原始的數據 41 import matplotlib.pyplot as plt 42 fig=plt.figure() 43 ax=fig.add_subplot(111) 44 ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0]) 45 #繪最佳擬合直線，需先要將點按照升序排列 46 xCopy=xMat.copy() 47 xCopy.sort(0) 48 yHat = xCopy*ws 49 ax.plot(xCopy[:,1],yHat) 50 plt.show()
幾乎任一數據集都可以用上述方法建立模型，只是需要判斷模型的好壞，計算預測值yHat和實際值yMat這兩個序列的相關系數，可以查看它們的匹配程度。
2、局部加權線性回歸
局部加權線性回歸給待預測點附近的每個點賦予一定的權重，用於解決線性回歸可能出現的欠擬合現象。與ｋＮＮ法類似，這種算法每次預測均需要事先選取出對應的數據子集，然后在這個子集上基於最小均分差來進行普通的回歸。該算法解出回歸系數的形式如下：
其中w是一個權重矩陣，通常采用核函數來對附近的點賦予權重，最常用的核函數是高斯核，如下：
這樣就構建了一個只含對角元素的權重矩陣W並且點x與x(i)越近，w(i,i)將會越大，k值控制衰減速度，且k值越小被選用於訓練回歸模型的數據集越小。
Python實現代碼：
1 def lwlr(testPoint,xArr,yArr,k=1.0): 2 '''局部加權線性回歸函數''' 3 xMat = mat(xArr); yMat = mat(yArr).T 4 m = shape(xMat)[0] 5 weights = mat(eye((m)))#創建對角矩陣 6 for j in range(m): 7 diffMat = testPoint - xMat[j,:] 8 #高斯核計算權重 9 weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2)) 10 xTx = xMat.T * (weights * xMat) 11 if linalg.det(xTx) == 0.0: 12 print("This matrix is singular, cannot do inverse") 13 return 14 ws = xTx.I * (xMat.T * (weights * yMat)) 15 return testPoint * ws 16 17 def lwlrTest(testArr,xArr,yArr,k=1.0): 18 '''為數據集中每個點調用lwlr()''' 19 m = shape(testArr)[0] 20 yHat = zeros(m) 21 for i in range(m): 22 yHat[i] = lwlr(testArr[i],xArr,yArr,k) 23 return yHat 24 25 if __name__ == "__main__": 26 '''局部加權線性回歸''' 27 xArr,yArr=loadDataSet('ex0.txt') 28 #擬合 29 yHat=lwlrTest(xArr,xArr,yArr,0.01) 30 #繪圖 31 xMat=mat(xArr) 32 yMat=mat(yArr) 33 srtInd = xMat[:,1].argsort(0) 34 xSort=xMat[srtInd][:,0,:] 35 import matplotlib.pyplot as plt 36 fig=plt.figure() 37 ax=fig.add_subplot(111) 38 ax.plot(xSort[:,1],yHat[srtInd]) 39 ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0],s=2,c='red') 40 plt.show()
k取0.01的結果
實際上，對k取不同值時有如下結果：
3、嶺回歸
如果數據的特征比樣本點多（n>m），也就是說輸入數據的矩陣x不是滿秩矩陣。而非滿秩矩陣在求逆時會出錯，所以此時不能使用之前的線性回歸方法。為解決這個問題，統計學家引入了嶺回歸的概念。
簡單來說，嶺回歸就是在矩陣xTx上加一個λI從而使得矩陣非奇異，進而能對 xTx+λI 求逆，其中I是一個mxm的單位矩陣。在這種情況下，回歸系數的計算公式將變成：
這里通過引入λ來限制了所有w之和，通過引入該懲罰項，能減少不重要的參數，這個技術在統計學中也叫縮減。
Python實現代碼：
1 def ridgeRegres(xMat,yMat,lam=0.2): 2 '''計算嶺回歸系數''' 3 xTx = xMat.T*xMat 4 denom = xTx + eye(shape(xMat)[1])*lam 5 if linalg.det(denom) == 0.0: 6 print("This matrix is singular, cannot do inverse") 7 return 8 ws = denom.I * (xMat.T*yMat) 9 return ws 10 11 def ridgeTest(xArr,yArr): 12 '''用於在一組lambda上測試結果''' 13 xMat = mat(xArr); yMat=mat(yArr).T 14 yMean = mean(yMat,0) 15 yMat = yMat - yMean #數據標准化 16 xMeans = mean(xMat,0) 17 xVar = var(xMat,0) 18 xMat = (xMat - xMeans)/xVar #所有特征減去各自的均值並除以方差 19 numTestPts = 30 #取30個不同的lambda調用函數 20 wMat = zeros((numTestPts,shape(xMat)[1])) 21 for i in range(numTestPts): 22 ws = ridgeRegres(xMat,yMat,exp(i-10)) 23 wMat[i,:]=ws.T 24 return wMat 25 26 if __name__ == "__main__": 27 '''嶺回歸''' 28 abX,abY=loadDataSet('abalone.txt') 29 ridgeWeights = ridgeTest(abX,abY)#得到30組回歸系數 30 #縮減效果圖 31 import matplotlib.pyplot as plt 32 fig=plt.figure() 33 ax=fig.add_subplot(111) 34 ax.plot(ridgeWeights) 35 plt.show()
運行之后得到下圖，橫軸表示第i組數據，縱軸表示該組數據對應的回歸系數值。從程序中可以看出lambda的取值為 exp(i-10) 其中i=0~29。所以結果圖的最左邊，即λ最小時，可以得到所有系數的原始值（與線性回歸一致）；而在右邊，系數全部縮減為0；在中間部分的某些值可以取得最好的預測效果。
4、前向逐步回歸
前向逐步回歸算法屬於一種貪心算法，即每一步盡可能減少誤差。一開始，所有的權重都設為1，然后每一步所做的決策是對某個權重增加或減少一個很小的值。
該算法偽代碼如下所示：
Python實現代碼：
1 def regularize(xMat): 2 '''數據標准化函數''' 3 inMat = xMat.copy() 4 inMeans = mean(inMat,0) 5 inVar = var(inMat,0) 6 inMat = (inMat - inMeans)/inVar 7 return inMat 8 9 def rssError(yArr,yHatArr): 10 '''計算均方誤差大小''' 11 return ((yArr-yHatArr)**2).sum() 12 13 def stageWise(xArr,yArr,eps=0.01,numIt=100): 14 ''' 15 逐步線性回歸算法 16 eps：表示每次迭代需要調整的步長 17 ''' 18 xMat = mat(xArr); yMat=mat(yArr).T 19 yMean = mean(yMat,0) 20 yMat = yMat - yMean 21 xMat = regularize(xMat) 22 m,n=shape(xMat) 23 returnMat = zeros((numIt,n)) #testing code remove 24 #為了實現貪心算法建立ws的兩份副本 25 ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy() 26 for i in range(numIt): 27 print(ws.T) 28 lowestError = inf; 29 for j in range(n):#對每個特征 30 for sign in [-1,1]:#分別計算增加或減少該特征對誤差的影響 31 wsTest = ws.copy() 32 wsTest[j] += eps*sign 33 yTest = xMat*wsTest 34 rssE = rssError(yMat.A,yTest.A) 35 #取最小誤差 36 if rssE < lowestError: 37 lowestError = rssE 38 wsMax = wsTest 39 ws = wsMax.copy() 40 returnMat[i,:]=ws.T 41 return returnMat 42 43 if __name__ == "__main__": 44 '''前向逐步線性回歸''' 45 abX,abY=loadDataSet('abalone.txt') 46 stageWise(abX,abY,0.01,200)
運行結果如下：
上述結果中值得注意的是w1和w6都是0，這表明它們不對目標值造成任何影響，也就是說這些特征很可能是不需要的。另外，第一個權重在0.04和0.05之間來回震盪，這是因為步長eps太大的緣故，一段時間后系數就已經飽和並在特定值之間來回震盪。
5、實例：預測樂高玩具套裝的價格
5.1 收集數據
原書介紹了從Google上在線獲取數據的方式，但是經測試該網址已經不可用，此處采用從離線網頁中爬取的方式收集數據。實現代碼如下：
1 def setDataCollect(retX, retY): 2 '''數據獲取方式一（不可用）''' 3 # searchForSet(retX, retY, 8288, 2006, 800, 49.99) 4 # searchForSet(retX, retY, 10030, 2002, 3096, 269.99) 5 # searchForSet(retX, retY, 10179, 2007, 5195, 499.99) 6 # searchForSet(retX, retY, 10181, 2007, 3428, 199.99) 7 # searchForSet(retX, retY, 10189, 2008, 5922, 299.99) 8 # searchForSet(retX, retY, 10196, 2009, 3263, 249.99) 9 '''數據獲取方式二''' 10 scrapePage("setHtml/lego8288.html","data/lego8288.txt",2006, 800, 49.99) 11 scrapePage("setHtml/lego10030.html","data/lego10030.txt", 2002, 3096, 269.99) 12 scrapePage("setHtml/lego10179.html","data/lego10179.txt", 2007, 5195, 499.99) 13 scrapePage("setHtml/lego10181.html","data/lego10181.txt", 2007, 3428, 199.99) 14 scrapePage("setHtml/lego10189.html","data/lego10189.txt", 2008, 5922, 299.99) 15 scrapePage("setHtml/lego10196.html","data/lego10196.txt", 2009, 3263, 249.99) 16 17 def scrapePage(inFile,outFile,yr,numPce,origPrc): 18 from bs4 import BeautifulSoup 19 fr = open(inFile,'r',encoding= 'utf8'); fw=open(outFile,'a') #a is append mode writing 20 soup = BeautifulSoup(fr.read()) 21 i=1 22 currentRow = soup.findAll('table', r="%d" % i) 23 while(len(currentRow)!=0): 24 title = currentRow[0].findAll('a')[1].text 25 lwrTitle = title.lower() 26 if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): 27 newFlag = 1.0 28 else: 29 newFlag = 0.0 30 soldUnicde = currentRow[0].findAll('td')[3].findAll('span') 31 if len(soldUnicde)==0: 32 print("item #%d did not sell" % i) 33 else: 34 soldPrice = currentRow[0].findAll('td')[4] 35 priceStr = soldPrice.text 36 priceStr = priceStr.replace('$','') #strips out $ 37 priceStr = priceStr.replace(',','') #strips out , 38 if len(soldPrice)>1: 39 priceStr = priceStr.replace('Free shipping', '') #strips out Free Shipping 40 print("%s\t%d\t%s" % (priceStr,newFlag,title)) 41 fw.write("%d\t%d\t%d\t%f\t%s\n" % (yr,numPce,newFlag,origPrc,priceStr)) 42 i += 1 43 currentRow = soup.findAll('table', r="%d" % i) 44 fw.close() 45 46 if __name__ == "__main__": 47 '''樂高玩具價格預測'''
48 #爬取數據
49 setDataCollect() 50 #讀取數據，這里已將以上方式獲取到的數據文本整合成為一個文件即legoAllData.txt 51 xmat,ymat = loadDataSet("data/legoAllData.txt")
5.2 訓練算法
首先我們用普通的線性回歸模型擬合數據看效果，擬合之前需要先添加對應常數項的特征X0
1 if __name__ == "__main__": 2 '''樂高玩具價格預測''' 3 #爬取數據 4 # setDataCollect() 5 #讀取數據，這里已將以上方式獲取到的數據文本整合成為一個文件即legoAllData.txt 6 # xMat,yMat = loadDataSet("data/legoAllData.txt") 7 #添加對應常數項的特征X0(X0=1) 8 lgX=mat(ones((76,5))) 9 lgX[:,1:5]=mat(xmat) 10 lgY=mat(ymat).T 11 12 #用標准回歸方程擬合 13 ws1=standRegres(lgX,mat(ymat)) #求標准回歸系數 14 yHat = lgX*ws1 #預測值 15 err1 = rssError(lgY.A,yHat.A) #計算平方誤差 16 cor1 = corrcoef(yHat.T,lgY.T) #計算預測值和真實值得相關性
測試結果為相關性cor1：0.7922，平方誤差和err1：3552526，顯然擬合效果還可以進一步提升。
接下來我們用交叉驗證測試嶺回歸：
1 def crossValidation(xArr,yArr,numVal=10): 2 ''' 3 交叉驗證測試嶺回歸 4 numVal:交叉驗證次數 5 ''' 6 m = len(yArr) 7 indexList = list(range(m)) 8 errorMat = zeros((numVal,30)) 9 for i in range(numVal): 10 trainX=[]; trainY=[] 11 testX = []; testY = [] 12 random.shuffle(indexList)#打亂順序 13 for j in range(m):#構建訓練和測試數據，10%用於測試 14 if j < m*0.9: 15 trainX.append(xArr[indexList[j]]) 16 trainY.append(yArr[indexList[j]]) 17 else: 18 testX.append(xArr[indexList[j]]) 19 testY.append(yArr[indexList[j]]) 20 wMat = ridgeTest(trainX,trainY) #30組不同參數下的回歸系數集 21 for k in range(30):#遍歷30個回歸系數集 22 matTestX = mat(testX); matTrainX=mat(trainX) 23 meanTrain = mean(matTrainX,0) 24 varTrain = var(matTrainX,0) 25 matTestX = (matTestX-meanTrain)/varTrain #用訓練參數標准化測試數據 26 yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)#預測值 27 errorMat[i,k]=rssError(yEst.T.A,array(testY))#計算預測平方誤差 28 # print(errorMat[i,k]) 29 #在完成所有交叉驗證后，errorMat保存了ridgeTest()每個lambda對應的多個誤差值 30 meanErrors = mean(errorMat,0)#計算每組平均誤差 31 minMean = float(min(meanErrors)) 32 bestWeights = wMat[nonzero(meanErrors==minMean)]#平均誤差最小的組的回歸系數即為所求最佳 33 #嶺回歸使用了數據標准化，而strandRegres()則沒有，因此為了將上述比較可視化還需將數據還原 34 xMat = mat(xArr); yMat=mat(yArr).T 35 meanX = mean(xMat,0); varX = var(xMat,0) 36 unReg = bestWeights/varX #還原后的回歸系數 37 constant = -1*sum(multiply(meanX,unReg)) + mean(yMat) #常數項 38 print("the best model from Ridge Regression is:\n",unReg) 39 print("with constant term: ",constant) 40 return unReg,constant 41 42 43 if __name__ == "__main__": 44 '''樂高玩具價格預測''' 45 #用交叉驗證測試嶺回歸 46 ws2,constant = crossValidation(xmat,ymat,10) 47 yHat2 = mat(xmat)*ws2.T + constant 48 err2 = rssError(lgY.A,yHat2.A) 49 cor2 = corrcoef(yHat2.T,lgY.T)
測試結果為相關性cor2：0.7874，平方誤差和err2：3827083，與最小二乘法比較好並沒有太大差異。其實這種分析方法使得我們可以挖掘大量數據的內在規律。在僅有4個特征時，該方法的效果也許並不明顯；但如果有100個以上的特征，該方法就會變得十分有效：它可以指出哪些特征是關鍵的，而哪些特征是不重要的。
THE END．