预测数值型数据:回归源码分析(1)
回歸模型比較簡(jiǎn)單,這里先簡(jiǎn)單介紹下,后面遇到難點(diǎn)再具體分析。
回歸的一般方法:
(1)收集數(shù)據(jù):采用任意方法收集數(shù)據(jù)。
(2)準(zhǔn)備數(shù)據(jù):回歸需要數(shù)值型數(shù)據(jù),標(biāo)稱型數(shù)據(jù)將被轉(zhuǎn)成二值型數(shù)據(jù)。
(3)分析數(shù)據(jù):繪出數(shù)據(jù)的可視化二維圖將有助于對(duì)數(shù)據(jù)做出理解和分析,在采用縮減法
求得新回歸系數(shù)之后,可以將新擬合線繪在圖上作為對(duì)比。
(4)訓(xùn)練算法:找到回歸系數(shù)。
(5)測(cè)試算法:使用幻或者預(yù)測(cè)值和數(shù)據(jù)的擬合度,來分析模型的效果。
(6)使用算法:使用回歸,可以在給定輸入的時(shí)候預(yù)測(cè)出一個(gè)數(shù)值,這是對(duì)分類方法的提
升,因?yàn)檫@樣可以預(yù)測(cè)連續(xù)型數(shù)據(jù)而不僅僅是離散的類別標(biāo)簽。
1. 用簡(jiǎn)單線性回歸找到最佳擬合直線
# -*- coding: utf-8 -*- """ Created on Wed Oct 25 16:49:50 2017 """ from numpy import * import matplotlib.pyplot as plt# 數(shù)據(jù)導(dǎo)入函數(shù) def loadDataSet(fileName): numFeat = len(open(fileName).readline().split('\t')) - 1 # 得到特征數(shù) dataMat = []; labelMat = []fr = open(fileName)for line in fr.readlines():lineArr =[]curLine = line.strip().split('\t')for i in range(numFeat):lineArr.append(float(curLine[i]))dataMat.append(lineArr)labelMat.append(float(curLine[-1])) # 得到最后一列目標(biāo)值return dataMat,labelMat# 標(biāo)準(zhǔn)回歸函數(shù) def standRegres(xArr,yArr):xMat = mat(xArr); yMat = mat(yArr).TxTx = xMat.T*xMat # 計(jì)算X'Xif linalg.det(xTx) == 0.0: # 判斷xTX行列式是否為0print "This matrix is singular, cannot do inverse"returnws = xTx.I * (xMat.T*yMat) # 求得當(dāng)前w的最優(yōu)估計(jì)值return ws# 局部加權(quán)線性回歸函數(shù) def lwlr(testPoint,xArr,yArr,k=1.0):xMat = mat(xArr); yMat = mat(yArr).Tm = shape(xMat)[0]weights = mat(eye((m))) # 產(chǎn)生對(duì)角單位權(quán)重矩陣for j in range(m): # 遍歷每一個(gè)樣本點(diǎn)diffMat = testPoint - xMat[j,:] # 所有的數(shù)據(jù)離測(cè)試點(diǎn)的距離weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2)) # 高斯核得到的權(quán)重,并且是對(duì)角陣xTx = xMat.T * (weights * xMat)if linalg.det(xTx) == 0.0:print "This matrix is singular, cannot do inverse"returnws = xTx.I * (xMat.T * (weights * yMat))return testPoint * ws# 為數(shù)據(jù)集每個(gè)點(diǎn)調(diào)用lwlr() def lwlrTest(testArr,xArr,yArr,k=1.0): # testArr,xArr是同一數(shù)據(jù)列m = shape(testArr)[0] # 得到數(shù)據(jù)樣本的個(gè)數(shù)yHat = zeros(m) # 0矩陣用于存儲(chǔ)預(yù)測(cè)值for i in range(m): # yHat[i] = lwlr(testArr[i],xArr,yArr,k)return yHat# 主函數(shù)xArr,yArr=loadDataSet('ex0.txt') ws=standRegres(xArr,yArr) # ws是回歸系數(shù) xMat=mat(xArr) # 是前兩維的數(shù)據(jù) yMat=mat(yArr) # 第三維的數(shù)據(jù),就是真實(shí)值 yHat=xMat*ws # 得到預(yù)測(cè)輸出的值yHat fig=plt.figure() ax=fig.add_subplot(111) # xMat[:,1].flatten().A[0]得到xMat第二維的數(shù)據(jù),yMat.T[:,0].flatten().A[0]是真實(shí)值 ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0]) # 創(chuàng)建出原始的數(shù)據(jù)圖像 xCopy=xMat.copy() # 拷貝 xCopy.sort(0) # 升排序 yHat=xCopy*ws # 得到預(yù)測(cè)值 ax.plot(xCopy[:,1],yHat,'r') # 畫出擬合直線 plt.show() # 計(jì)算相關(guān)系數(shù) yHat=xMat*ws print 'Correlation coefficient:',corrcoef(yHat.T,yMat) #numpy庫(kù)提供了相關(guān)系數(shù)的計(jì)算方法注意:
這里的回歸系數(shù)是用最小化平方誤差得到的,即
w^=(XTX)?1XTy
和代碼中的ws = xTx.I * (xMat.T*yMat)相對(duì)應(yīng)。通過計(jì)算兩個(gè)序列的相關(guān)系數(shù), 可以計(jì)算預(yù)測(cè)值和真實(shí)值之間的匹配程度。Numpy庫(kù)提供了相關(guān)系數(shù)的計(jì)算方法:通過命令corrcoef(yEstimate,yActual)來計(jì)算預(yù)測(cè)和真實(shí)值的相關(guān)性。
線性回歸是得到最佳的回歸系數(shù)后,再對(duì)原始輸入樣本進(jìn)行預(yù)測(cè),此處的回歸系數(shù)對(duì)于所有的樣本都是一樣的(yHat=xMat*ws),這點(diǎn)不同于局部線性回歸,然后畫出原始數(shù)據(jù)的散點(diǎn)圖和回歸線。
運(yùn)行結(jié)果:
相關(guān)系數(shù):
Correlation coefficient: [[ 1. 0.98647356][ 0.98647356 1. ]]2. 局部加權(quán)線性回歸
線性回歸的一個(gè)問題是有可能出現(xiàn)欠擬合現(xiàn)象,因?yàn)樗蟮氖蔷哂凶钚【秸`差的無偏估計(jì)。顯而易見,如果模型欠擬合將不能取得最好的預(yù)測(cè)效果。所以有些方法允許在估計(jì)中引人一些偏差,從而降低預(yù)測(cè)的均方誤差。
其中的一個(gè)方法是局部加權(quán)線性回歸(LWLR)。在該算法中,我們給待預(yù)測(cè)點(diǎn)附近的每個(gè)點(diǎn)賦予一定的權(quán)重;然后在這個(gè)子集上基于最小均方差來進(jìn)行普通的回歸。與KNN一樣,這種算法每次預(yù)測(cè)均需要事先選取出對(duì)應(yīng)的數(shù)據(jù)子集。該算法解出回歸系數(shù)w的形式如下:
w^=(XTWX)?1XTWy
其中W是一個(gè)矩陣,用來給每個(gè)數(shù)據(jù)點(diǎn)賦予權(quán)重。
LWLR使用“核”(與支持向量機(jī)中的核類似)來對(duì)附近的點(diǎn)賦予更高的權(quán)重。最常用的核就是高斯核,高斯核對(duì)應(yīng)的權(quán)重如下:
w(i,j)=exp(|xi?x|?2k2)
這樣就可以得到要預(yù)測(cè)的點(diǎn)x(i)附近的點(diǎn)的權(quán)重值。
把上述的代碼中的主函數(shù)改為如下:
# 主函數(shù) xArr,yArr=loadDataSet('ex0.txt') #print 'Weighted prediction:',lwlr(xArr[0],xArr,yArr,1.0) yHat_1=lwlrTest(xArr,xArr,yArr,1.0) yHat_2=lwlrTest(xArr,xArr,yArr,0.01) yHat_3=lwlrTest(xArr,xArr,yArr,0.003) # 對(duì)xArr排序 xMat=mat(xArr) # 轉(zhuǎn)化為矩陣 srtInd=xMat[:,1].argsort(0) # 按升序排列,并返回對(duì)應(yīng)的索引值 xSort=xMat[srtInd][:,0,:] # 得到二維的矩陣,重新升序后得到的矩陣 #xx=xMat[srtInd] # 三維的矩陣 fig=plt.figure() ax=fig.add_subplot(111) ax.plot(xSort[:,1],yHat_1[srtInd]) # 預(yù)測(cè)值也按相應(yīng)的坐標(biāo) ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0],s=2,c='r') plt.show()fig=plt.figure() ax=fig.add_subplot(111) ax.plot(xSort[:,1],yHat_2[srtInd]) # 預(yù)測(cè)值也按相應(yīng)的坐標(biāo) ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0],s=2,c='r') plt.show()fig=plt.figure() ax=fig.add_subplot(111) ax.plot(xSort[:,1],yHat_3[srtInd]) # 預(yù)測(cè)值也按相應(yīng)的坐標(biāo) ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0],s=2,c='r') plt.show()運(yùn)行結(jié)果:
可以看到從左到右的系數(shù)k分別為:1.0,0.01,0.003 的擬合效果。
參數(shù)k決定了對(duì)待測(cè)點(diǎn)附近的點(diǎn)賦予多大的權(quán)重,控制著權(quán)重衰減的速度。
分析:
其中待預(yù)測(cè)點(diǎn)附近點(diǎn)的權(quán)重的代碼:
diffMat*diffMat.T:是每個(gè)待測(cè)點(diǎn)與樣本的每個(gè)點(diǎn)的距離。
weights:是基于每個(gè)待測(cè)點(diǎn)與全部樣本點(diǎn)進(jìn)行計(jì)算得到的權(quán)重對(duì)角陣。
return testPoint * ws :返回待測(cè)點(diǎn)的預(yù)測(cè)值,而每個(gè)待測(cè)點(diǎn)的權(quán)重都要基于全部樣本進(jìn)行計(jì)算。
ax.plot(xSort[:,1],yHat_1[srtInd]) :畫出回歸線
ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0],s=2,c='r'):畫出原始數(shù)據(jù)的散點(diǎn)圖。
要注意的語法:
srtInd=xMat[:,1].argsort(0) # 按升序排列,并返回對(duì)應(yīng)的索引值 xx=xMat[srtInd] # 三維的矩陣 xSort=xMat[srtInd][:,0,:] # 得到二維的矩陣,重新升序后得到的矩陣xMat[索引值][:,0,;]其得到的是:二維的矩陣,而xMat[索引值]得到的是三維的形式,效果如下:
In [21]: xSort Out[21]: matrix([[ 1. , 0.014855],[ 1. , 0.015371],[ 1. , 0.033859],..., [ 1. , 0.990947],[ 1. , 0.993888],[ 1. , 0.995731]])In [22]: xx Out[22]: matrix([[[ 1. , 0.014855]],[[ 1. , 0.015371]],[[ 1. , 0.033859]],..., [[ 1. , 0.990947]],[[ 1. , 0.993888]],[[ 1. , 0.995731]]])In [23]: shape(xx) Out[23]: (200L, 1L, 2L)In [24]: shape(xSort) Out[24]: (200L, 2L)In [25]:3. 預(yù)測(cè)鮑魚的年齡
添加計(jì)算誤差大小的函數(shù),并改變主函數(shù)如下:
# 計(jì)算預(yù)測(cè)誤差的大小 def rssError(yArr,yHatArr): # yArr和yHatArr都需要是數(shù)組return ((yArr-yHatArr)**2).sum()# 主函數(shù)abX,abY=loadDataSet('abalone.txt') yHat01=lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1) yHat1=lwlrTest(abX[0:99],abX[0:99],abY[0:99],1) yHat10=lwlrTest(abX[0:99],abX[0:99],abY[0:99],10) print 'k=0.1 :',rssError(abY[0:99],yHat01.T) # 不同的核的預(yù)測(cè)預(yù)測(cè)誤差結(jié)果 print 'k=1 :',rssError(abY[0:99],yHat1.T) print 'k=10 :',rssError(abY[0:99],yHat10.T) # 測(cè)試新數(shù)據(jù)的表現(xiàn) print 'Forecast new data....' yHat01=lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1) # 此處的權(quán)重計(jì)算是基于前100個(gè)樣本的,但預(yù)測(cè)的是另100個(gè) yHat1=lwlrTest(abX[100:199],abX[0:99],abY[0:99],1) yHat10=lwlrTest(abX[100:199],abX[0:99],abY[0:99],10) print 'k=0.1 :',rssError(abY[100:199],yHat01.T) print 'k=1 :',rssError(abY[100:199],yHat1.T) print 'k=10 :',rssError(abY[100:199],yHat10.T) # 和簡(jiǎn)單的線性回歸對(duì)比 print 'compare to Simple regression:....' ws=standRegres(abX[0:99],abY[0:99]) yHat=mat(abX[100:199])*ws print 'Simple regression:',rssError(abY[100:199],yHat.T.A)運(yùn)行結(jié)果:
k=0.1 : 56.7836322443 k=1 : 429.89056187 k=10 : 549.118170883 Forecast new data.... k=0.1 : 14651.309989 k=1 : 573.52614419 k=10 : 517.571190538 compare to Simple regression.... Simple regression: 518.636315325分析:
- 使用較小的核將得到較低的誤差。那么,為什么不在所有數(shù)據(jù)集上都使用最小的核呢?這是因?yàn)槭褂米钚〉暮藢⒃斐蛇^擬合,對(duì)新數(shù)據(jù)不一定能達(dá)到最好的預(yù)測(cè)效果
- 簡(jiǎn)單線性回歸達(dá)到了與局部加權(quán)線性回歸類似的效果。這也表明一點(diǎn),必須在未知數(shù)據(jù)上比較效果才能選取到最佳模型。那么最佳的核大小是10嗎?或許是,但如果想得到更好的效果,應(yīng)該用10個(gè)不同的樣本集做10次測(cè)試來比較結(jié)果。
- 預(yù)測(cè)鮑魚的年齡的例子展示了如何使用局部加權(quán)線性回歸來構(gòu)建模型,可以得到比普通線性回歸更好的效果。局部加權(quán)線性回歸的問題在于,每次必須在整個(gè)數(shù)據(jù)集上運(yùn)行。也就是說為了做出預(yù)測(cè),必須保存所有的訓(xùn)練數(shù)據(jù)。這是其劣勢(shì)所在。
總結(jié)
以上是生活随笔為你收集整理的预测数值型数据:回归源码分析(1)的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 玩具手枪哪里最容易坏
- 下一篇: 预测数值型数据:回归 源码分析(2)