python实战之线性回归、局部加权回归
1.基本概念与思想
回归：求回归方程中回归系数的过程称为回归。
局部加权思想：给待预测点附近的每个点赋予一定的权重。
2.线性回归
回归方程的解： Θ=(XTX)-1XTY (1)
其中，Θ表示回归系数矩阵，X表示样本矩阵，Y表示样本类标矩阵。
3.局部加权回归
回归方程解： Θ=(XTWX)-1XTWY (2)
其中，（2）与（1）不同的是多了表示局部权重的矩阵W。另外，
w(i,i)=exp((x(i)-x)2/-2k2) (3)
4.python代码实现
1 #!/usr/bin/python 2 #-*- coding:utf-8 -*- 3 4 from numpy import * 5 import matplotlib.pyplot as plt 6 ''' 7 解决python matplotlib画图无法显示中文的问题！ 8 ''' 9 from pylab import * 10 mpl.rcParams['font.sans-serif']=['SimHei'] 11 mpl.rcParams['axes.unicode_minus']=False 12 ############################################################# 13 14 def loadDataSet(fileName): 15 numFeat=len(open(fileName).readline().split(' '))-1 16 dataMat=[];labelMat=[] #空列表 17 #为了减少内存消耗，一行行读入 18 #fr=open(fileName) 19 #for line in fr.readlines(): 20 for line in open(fileName): 21 lineArr=[] 22 curLine=line.strip().split(' ') 23 for i in range(numFeat): 24 lineArr.append(float(curLine[i])) #list方法append() 25 dataMat.append(lineArr) 26 labelMat.append(float(curLine[-1])) 27 return dataMat,labelMat 28 #线性回归(lr)主函数 29 def standRegression(xArr,yArr): 30 xMat=mat(xArr);yMat=mat(yArr).T #注意此处需要转置 31 xTx=xMat.T*xMat 32 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 return ws 37 #lr绘图 38 def lrPlot(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*ws 48 #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).T 57 m=shape(xMat)[0] 58 weights=mat(eye(m)) 59 for j in range(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*ws 68 def lwlrTest(testArr,xArr,yArr,k=1.0): 69 #testArr=mat(testArr) #转换为矩阵 70 m=shape(testArr)[0] 71 yHat=zeros(m) 72 for i in range(m): 73 yHat[i]=lwlr(testArr[i],xArr,yArr,k) #testArr[i]列表索引 74 return yHat 75 #lwlr绘图测试版 76 def lwlrPlot(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()
代码主要函数：loadDataSet(), standRegression(), lwlr(), lwlrTest()
转载请注明出处，谢谢！