python实战之线性回归、局部加权回归

转载 http://www.cnblogs.com/cygalaxy/p/6810249.html

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代码实现

#!/usr/bin/python
#-*- coding:utf-8 -*-

from numpy import *
import matplotlib.pyplot as plt
'''
解决python matplotlib画图无法显示中文的问题！
'''
from pylab import *
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False
#############################################################

def loadDataSet(fileName):
numFeat=len(open(fileName).readline().split('\t'))-1
dataMat=[];labelMat=[] #空列表
#为了减少内存消耗，一行行读入
#fr=open(fileName)
#for line in fr.readlines():
for line in open(fileName):
lineArr=[]
curLine=line.strip().split('\t')
for i in range(numFeat):
lineArr.append(float(curLine[i])) #list方法append()
dataMat.append(lineArr)
labelMat.append(float(curLine[-1]))
return dataMat,labelMat
#线性回归(lr)主函数
def standRegression(xArr,yArr):
xMat=mat(xArr);yMat=mat(yArr).T #注意此处需要转置
xTx=xMat.T*xMat
if linalg.det(xTx)==0.0: #若行列式为0，则不可逆
print 'This matrix is singular,cannot do inverse'
return
ws=xTx.I*(xMat.T*yMat)
return ws
#lr绘图
def lrPlot(xArr,yArr,ws):
xMat=mat(xArr);yMat=mat(yArr)
fig=plt.figure()
ax=fig.add_subplot(111)
#scatter()画散点图，yMat绘制原始图
ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0],label='yMat') #scatter(x,y),x和y必须转化为1-D
#yHat绘制拟合图
xCopy=xMat.copy()
xCopy.sort(0) #按列排序，画直线图需要排序
yHat=xCopy*ws
#plot()画直线图
ax.plot(xCopy[:,1],yHat,label='yHat')

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

#局部线性回归(lwlr)主函数
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat=mat(xArr);yMat=mat(yArr).T
m=shape(xMat)[0]
weights=mat(eye(m))
for j in range(m):
diffMat=testPoint-xMat[j,:]
weights[j,j]=exp(diffMat*diffMat.T/(-2*k**2))
xTwx=xMat.T*(weights*xMat) #
if linalg.det(xTwx)==0.0:
print 'This matrix is singular,cannot do inverse'
return
ws=xTwx.I*(xMat.T*(weights*yMat))
return testPoint*ws
def lwlrTest(testArr,xArr,yArr,k=1.0):
#testArr=mat(testArr) #转换为矩阵
m=shape(testArr)[0]
yHat=zeros(m)
for i in range(m):
yHat[i]=lwlr(testArr[i],xArr,yArr,k) #testArr[i]列表索引
return yHat
#lwlr绘图测试版
def lwlrPlot(xArr,yArr,yHat):
xMat=mat(xArr)
srtInd=xMat[:,1].argsort(0) #等价于argsort(xMat[:,1],0)
xSort=xMat[srtInd][:,0,:] #等价于xMat[srtInd.flatten().A[0]]

fig=plt.figure()
ax=fig.add_subplot(111)

#直线图plt.plot(),画plot前要排序
#ax.plot(xMat[:,1],yHat[:].T)
ax.plot(xSort[:,1],yHat[srtInd])

#画散点图不需要排序
ax.scatter(xMat[:,1].flatten().A[0],mat(yHat).T.flatten().A[0],s=2,c='k')
ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r') #散点图plt.scatter()

plt.show()
#lwlr绘图比较3个yHat(k取1,0.01,0.003)
def lwlrPlot3(xArr,yArr): #输入：xArr是n×d矩阵/数组/列表；yArr是n×1

xMat=mat(xArr)
srtInd=xMat[:,1].argsort(0) #等价于argsort(xMat[:,1],0)
xSort=xMat[srtInd][:,0,:] #等价于xMat[srtInd.flatten().A[0]]

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

fig=plt.figure()
ax1=fig.add_subplot(311)
ax2=fig.add_subplot(312)
ax3=fig.add_subplot(313)

#画直线图需要排序
#直线图plt.plot(),plot前要排序
#ax1.plot(xMat[:,1],yHat[:].T)
ax1.plot(xSort[:,1],yHat1[srtInd])
ax2.plot(xSort[:,1],yHat2[srtInd])
ax3.plot(xSort[:,1],yHat3[srtInd])

#画散点图不需要排序
ax1.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'欠拟合') #散点图plt.scatter()
ax2.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'最好')
ax3.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='r',label=u'过拟合')

ax1.legend(loc='upper left')
ax2.legend(loc='upper left')
ax3.legend(loc='upper left')

plt.show()