python线性加权回归_python--线性回归、局部加权回归

1 #!/usr/bin/python

2 #-*- coding:utf-8 -*-

3

4 from numpy import *

5 importmatplotlib.pyplot as plt6 '''

7 解决python matplotlib画图无法显示中文的问题！8 '''

9 from pylab import *

10 mpl.rcParams['font.sans-serif']=['SimHei']11 mpl.rcParams['axes.unicode_minus']=False12 #############################################################

13

14 defloadDataSet(fileName):15 numFeat=len(open(fileName).readline().split('\t'))-1

16 dataMat=[];labelMat=[] #空列表

17 #为了减少内存消耗，一行行读入

18 #fr=open(fileName)

19 #for line in fr.readlines():

20 for line inopen(fileName):21 lineArr=[]22 curLine=line.strip().split('\t')23 for i inrange(numFeat):24 lineArr.append(float(curLine[i])) #list方法append()

25 dataMat.append(lineArr)26 labelMat.append(float(curLine[-1]))27 returndataMat,labelMat28 #线性回归(lr)主函数

29 defstandRegression(xArr,yArr):30 xMat=mat(xArr);yMat=mat(yArr).T #注意此处需要转置

31 xTx=xMat.T*xMat32 if linalg.det(xTx)==0.0: #若行列式为0，则不可逆

33 print 'This matrix is singular,cannot do inverse'

34 return

35 ws=xTx.I*(xMat.T*yMat)36 returnws37 #lr绘图

38 deflrPlot(xArr,yArr,ws):39 xMat=mat(xArr);yMat=mat(yArr)40 fig=plt.figure()41 ax=fig.add_subplot(111)42 #scatter()画散点图，yMat绘制原始图

43 ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0],label='yMat') #scatter(x,y),x和y必须转化为1-D

44 #yHat绘制拟合图

45 xCopy=xMat.copy()46 xCopy.sort(0) #按列排序，画直线图需要排序

47 yHat=xCopy*ws48 #plot()画直线图

49 ax.plot(xCopy[:,1],yHat,label='yHat')50

51 plt.legend(loc='down left') #指定方框位置

52 plt.show()53

54 #局部线性回归(lwlr)主函数

55 def lwlr(testPoint,xArr,yArr,k=1.0):56 xMat=mat(xArr);yMat=mat(yArr).T57 m=shape(xMat)[0]58 weights=mat(eye(m))59 for j inrange(m):60 diffMat=testPoint-xMat[j,:]61 weights[j,j]=exp(diffMat*diffMat.T/(-2*k**2))62 xTwx=xMat.T*(weights*xMat) #63 if linalg.det(xTwx)==0.0:64 print 'This matrix is singular,cannot do inverse'

65 return

66 ws=xTwx.I*(xMat.T*(weights*yMat))67 return testPoint*ws68 def lwlrTest(testArr,xArr,yArr,k=1.0):69 #testArr=mat(testArr) #转换为矩阵

70 m=shape(testArr)[0]71 yHat=zeros(m)72 for i inrange(m):73 yHat[i]=lwlr(testArr[i],xArr,yArr,k) #testArr[i]列表索引

74 returnyHat75 #lwlr绘图测试版

76 deflwlrPlot(xArr,yArr,yHat):77 xMat=mat(xArr)78 srtInd=xMat[:,1].argsort(0) #等价于argsort(xMat[:,1],0)

79 xSort=xMat[srtInd][:,0,:] #等价于xMat[srtInd.flatten().A[0]]

80

81 fig=plt.figure()82 ax=fig.add_subplot(111)83

84 #直线图plt.plot(),画plot前要排序

85 #ax.plot(xMat[:,1],yHat[:].T)

86 ax.plot(xSort[:,1],yHat[srtInd])87

88 #画散点图不需要排序

89 ax.scatter(xMat[:,1].flatten().A[0],mat(yHat).T.flatten().A[0],s=2,c='k')90 ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r') #散点图plt.scatter()

91

92 plt.show()93 #lwlr绘图比较3个yHat(k取1,0.01,0.003)

94 def lwlrPlot3(xArr,yArr): #输入：xArr是n×d矩阵/数组/列表；yArr是n×1

95

96 xMat=mat(xArr)97 srtInd=xMat[:,1].argsort(0) #等价于argsort(xMat[:,1],0)

98 xSort=xMat[srtInd][:,0,:] #等价于xMat[srtInd.flatten().A[0]]

99

100 yHat1=lwlrTest(xArr,xArr,yArr,1) #调用局部加权回归(lwlr)主函数

101 yHat2=lwlrTest(xArr,xArr,yArr,0.01)102 yHat3=lwlrTest(xArr,xArr,yArr,0.03)103

104 fig=plt.figure()105 ax1=fig.add_subplot(311)106 ax2=fig.add_subplot(312)107 ax3=fig.add_subplot(313)108

109 #画直线图需要排序

110 #直线图plt.plot(),plot前要排序

111 #ax1.plot(xMat[:,1],yHat[:].T)

112 ax1.plot(xSort[:,1],yHat1[srtInd])113 ax2.plot(xSort[:,1],yHat2[srtInd])114 ax3.plot(xSort[:,1],yHat3[srtInd])115

116 #画散点图不需要排序

117 ax1.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'欠拟合') #散点图plt.scatter()

118 ax2.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'最好')119 ax3.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'过拟合')120

121 ax1.legend(loc='upper left')122 ax2.legend(loc='upper left')123 ax3.legend(loc='upper left')124

125 plt.show()