本文深入探讨了机器学习中的回归分析技术，包括线性回归、局部加权线性回归、岭回归和前向逐步回归等方法，详细介绍了每种方法的原理、计算过程及其实现代码。
一、用线性回归预测最佳拟合直线
回归的目的是用已知的回归方程预测数值型的目标值，求回归方程系数的过程叫做回归。具体的做法是用回归系数乘以输入值，再将结果全部加起来，就得到预测值。
求回归系数的一个常用方法是找出使得误差最小的www，这里的误差指的是预测yyy值和真实yyy值之间的差值，使用误差的简单累加正差值和负差值相互抵消，所以采用平方误差，而基于均方误差最小化来进行模型求解的方法称为“最小二乘法”。其中平方误差可以写为：
Ew=∑i=1m(yi−xiw)2E_{w}=\sum_{i=1}^{m}(y_{i}-x_{i}w)^{2}Ew=∑i=1m(yi−xiw)2
用矩阵表示还可以写作Ew=(y−Xw)T(y−Xw)E_{w}=(y-Xw)^{T}(y-Xw)Ew=(y−Xw)T(y−Xw),对www求导得：
Ew=yTy−wTXTy−yXw+wTXTXwE_{w}=y^{T}y-w^{T}X^{T}y-yXw+w^{T}X^{T}XwEw=yTy−wTXTy−yXw+wTXTXw
∂∂wEw=−XTy−XTy+XTXw+XTXw\frac{\partial }{\partial w}E_{w}=-X^{T}y-X^{T}y+X^{T}Xw+X^{T}Xw∂w∂Ew=−XTy−XTy+XTXw+XTXw
∂∂wEw=2XTXw−2XTy=2XT(Xw−y)\frac{\partial }{\partial w}E_{w}=2X^{T}Xw-2X^{T}y=2X^{T}(Xw-y)∂w∂Ew=2XTXw−2XTy=2XT(Xw−y)
这里用到对矩阵的转置求导，即:∂∂AAB=∂∂ABA=BT\frac{\partial }{\partial A}AB=\frac{\partial }{\partial A}BA=B^{T}∂A∂AB=∂A∂BA=BT,∂∂AATB=∂∂ABAT=B\frac{\partial }{\partial A}A^{T}B=\frac{\partial }{\partial A}BA^{T}=B∂A∂ATB=∂A∂BAT=B
令∂∂wEw=0\frac{\partial }{\partial w}E_{w}=0∂w∂Ew=0可得www的最优的闭式解，当XTXX^{T}XXTX为满秩矩阵或者正定矩阵时可得：
w=(XTX)−1XTyw=(X^{T}X)^{-1}X^{T}yw=(XTX)−1XTy
基于Python实现的代码如下所示：
def standRegres(xArr,yArr):
xMat=mat(xArr);yMat=mat(yArr)
xTx=xMat.T*xMat #计算X.T*X
if linalg.det(xTx)==0.0:
print("This matrix is singular,cannot do inverse")
return
ws=xTx.I*(xMat.T*yMat.T) #xTx.I为xTx的逆
return ws
画出数据集和最佳拟合直线Python代码如下所示：
def plotScatter(xArr,yArr,ws):
xMat=mat(xArr)
yMat=mat(yArr)
yHat=xMat*ws
coef=corrcoef(yHat.T,yMat) #计算相关系数
print("相关系数为：",coef)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0],c='r',s=10) #flatten函数将矩阵降成一维，.A表示将矩阵转换成数组
xCopy=xMat.copy()
xCopy.sort(0) #按列排序结果
yHat=xCopy*ws
ax.plot(xCopy[:,1],yHat)
plt.savefig('最佳拟合直线示意图.png')
plt.show()
得到的结果如下图所示：
二、局部加权线性回归
线性回归可能出现欠拟合现象，因为它要求的是具有最小均方误差的无偏估计。如果模型欠拟合则不能去的较好的预测效果。所以局部加权线性回归在估计中引入一些偏差，从而降低预测的均方误差。
在局部加权线性回归(LWLR)算法中,给带预测的点附近的每一个点赋予一定的权重，在这个子集上基于最小均方误差进行普通回归。与KNN一样，该算法每次预测均需要事先选取出对应的数据子集，该算法解出回归系数www的形式如下：
w=(XtWX)−1XTWyw=(X^{t}WX)^{-1}X^{T}Wyw=(XtWX)−1XTWy
其中WWW是一个矩阵，用来给每个数据点赋予权重。
局部加权线性回归使用"核"来对附近的点赋予更高的权重。最常用的就是高斯核，高斯核对应的权重如下所示：
w(i,i)=exp(∣xi−x∣−2k2)w(i,i)=exp\left ( \frac{\left | x^{i}-x \right |}{-2k^{2}} \right )w(i,i)=exp(−2k2∣xi−x∣)
这样就构建了一个只含对角线元素的权重矩阵WWW,并且点xxx与x(i)x(i)x(i)越近，w(i,i)w(i,i)w(i,i)将会越大。kkk为用户指定的参数，它决定对附近点赋予多大的权重。
基于Python实现局部加权线性回归以及画出拟合直线的代码如下所示：
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat=mat(xArr);yMat=mat(yArr)
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.0*k**2)) #随着距离的增加，权重值的大小以指数级衰减
xTx=xMat.T*(weights*xMat)
if linalg.det(xTx)==0.0:
print('This matrix is singular,cannot do inverse')
return
ws=xTx.I*(xMat.T*(weights*yMat.T))
return testPoint*ws #得到回归系数ws的一个估计
def lwlrTest(testArr,xArr,yArr,k=1.0):
m=shape(testArr)[0]
yHat=zeros(m)
for i in range(m):
yHat[i]=lwlr(testArr[i],xArr,yArr,k) #给定空间中的任意一点，计算出对应的预测值yHat
return yHat
def plotScatterlwlr(xArr,yArr,yHat):
xMat = mat(xArr)
yMat = mat(yArr)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xMat[:, 1].flatten().A[0], yMat.flatten().A[0], c='r', s=2) # flatten函数将矩阵降成一维
srtInd=xMat[:,1].argsort(0)
xSort=xMat[srtInd][:,0,:] #返回的是数组值从小到大的索引值
ax.plot(xSort[:,1], yHat[srtInd])
plt.savefig('k=0.01时局部加权线性回归示意图.png')
plt.show()
当k=0.01时k=0.01时k=0.01时，拟合的直线结果如下所示
当k=0.003k=0.003k=0.003考虑了太多的噪声，从而出现过拟合现象，拟合直线如下图所示：
当k=1k=1k=1时，拟合的效果与最小二乘法差不多
三、岭回归
当特征比样本点还多(n>m)(n>m)(n>m),也就是说矩阵XXX不是满秩的，非满秩矩阵在求逆时会出现问题。统计学家引入了岭回归的概念。简单来说岭回归就是在矩阵XTXX^{T}XXTX上加一个λI\lambda IλI,从而使得矩阵非奇异，进而能使得(XTX+λI)(X^{T}X+\lambda I)(XTX+λI)求逆。其中III是一个m∗mm*mm∗m的单位矩阵。在这种情况下回归系数的计算公式变为：
w=(XTX+λI)−1XTyw=(X^{T}X+\lambda I)^{-1}X^{T}yw=(XTX+λI)−1XTy
这里通过引入λ\lambdaλ来限制所有www之和，通过引入该惩罚项，能够减少不重要的参数。岭回归以及回归系数变化图实现代码如下所示：
#岭回归
def ridgeRegres(xMat,yMat,lam=0.2):
xTx=xMat.T*xMat
denom=xTx+eye(shape(xMat)[1])*lam
if linalg.det(denom)==0.0:
print('This matrix is singular,connot do inverse')
return
ws=denom.I*(xMat.T*yMat.T)
return ws
def ridgeTest(xArr,yArr):
xMat=mat(xArr);yMat=mat(yArr)
yMean=mean(yMat.T,0)
yMat=yMat-yMean
xMeans=mean(xMat,0)
xVar=var(xMat,0)
xMat=(xMat-xMeans)/xVar
numTestPts=30
wMat=zeros((numTestPts,shape(xMat)[1]))
for i in range(numTestPts):
ws=ridgeRegres(xMat,yMat,exp(i-10))
wMat[i,:]=ws.T
return wMat
def plotRidge(ridgeWeights):
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.savefig('岭回归的回归系数变化图.png')
plt.show()
岭回归的系数变化图如下所示：
四、前向逐步回归
前向逐步回归属于一种贪心算法，对每一步都尽可能的减少误差。其伪代码如下所示：
数据标准化，使其分布满足0均值和单位方差
在每轮迭代过程中：
设置当前最小误差lowestError为正无穷
对每个特征：
增大或者缩小：
改变一个系数得到一个新的w
计算新的w下的误差
如果误差Error小于最小当前误差lowestError:
设置wbest等于当前的w
将w设置为新的wbest
基于Python 实现的代码如下所示：
def stageWise(xArr,yArr,eps,numIt):
#数据标准化
xMat=mat(xArr);yMat=mat(yArr).T
yMean=mean(yMat,0)
yMat=yMat-yMean
xMeans=mean(xMat,0)
xVar=var(xMat,0)
xMat=(xMat-xMeans)/xVar
m,n=shape(xMat)
returnMat=zeros((numIt,n))
ws=zeros((n,1))
wsTest=ws.copy()
wsMax=ws.copy()
#每轮迭代过程
for i in range(numIt):
print(ws.T)
lowesError=inf
for j in range(n):
for sign in [-1,1]:
wsTest=ws.copy()
wsTest[j]+=eps*sign
yTest=xMat*wsTest
rssE=rssError(yMat.A,yTest.A)
if rssE<lowesError:
lowesError=rssE
wsMax=wsTest
ws=wsMax.copy()
returnMat[i,:]=ws.T
return returnMat
#分析预测误差大小
def rssError(yArr,yHatArr):
return((yArr-yHatArr)**2).sum()
def plotlasso(returnWeights):
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(returnWeights)
plt.savefig('前向逐步回归系数变化图.png')
plt.show()
前向逐步回归得到的系数与迭代次数的关系如下所示：
五、全部代码实现
工程打包地址：
"""
内容：
1.机器学习-回归分析
2.标准回归函数
3.局部加权线性回归
4.岭回归
5.前向逐步回归
姓名：pcb
日期：2019.1.6
"""
from numpy import *
import matplotlib.pyplot as plt
#加载数据
def loadDataSet(filename):
numFeat=len(open(filename).readline().split('\t'))-1
datMat=[];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]))
datMat.append(lineArr)
labelMat.append(float(curLine[-1]))
return datMat,labelMat
#1.---------------计算最佳拟合直线---------------------------------------------
def standRegres(xArr,yArr):
xMat=mat(xArr);yMat=mat(yArr)
xTx=xMat.T*xMat #计算X.T*X
if linalg.det(xTx)==0.0:
print("This matrix is singular,cannot do inverse")
return
ws=xTx.I*(xMat.T*yMat.T) #xTx.I为xTx的逆
return ws
#画出线性回归的散点图以及最佳拟合曲线
def plotScatter(xArr,yArr,ws):
xMat=mat(xArr)
yMat=mat(yArr)
yHat=xMat*ws
coef=corrcoef(yHat.T,yMat) #计算相关系数
print("相关系数为：",coef)
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(xMat[:,1].flatten().A[0],yMat.T[:,0].flatten().A[0],c='r',s=10) #flatten函数将矩阵降成一维，.A表示将矩阵转换成数组
xCopy=xMat.copy()
xCopy.sort(0) #按列排序结果
yHat=xCopy*ws
ax.plot(xCopy[:,1],yHat)
plt.savefig('最佳拟合直线示意图.png')
plt.show()
#-----------------------------------------------------------------------------
#2.----------------局部加权线性回归---------------------------------------------
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat=mat(xArr);yMat=mat(yArr)
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.0*k**2)) #随着距离的增加，权重值的大小以指数级衰减
xTx=xMat.T*(weights*xMat)
if linalg.det(xTx)==0.0:
print('This matrix is singular,cannot do inverse')
return
ws=xTx.I*(xMat.T*(weights*yMat.T))
return testPoint*ws #得到回归系数ws的一个估计
def lwlrTest(testArr,xArr,yArr,k=1.0):
m=shape(testArr)[0]
yHat=zeros(m)
for i in range(m):
yHat[i]=lwlr(testArr[i],xArr,yArr,k) #给定空间中的任意一点，计算出对应的预测值yHat
return yHat
def plotScatterlwlr(xArr,yArr,yHat):
xMat = mat(xArr)
yMat = mat(yArr)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xMat[:, 1].flatten().A[0], yMat.flatten().A[0], c='r', s=2) # flatten函数将矩阵降成一维
srtInd=xMat[:,1].argsort(0)
xSort=xMat[srtInd][:,0,:] #返回的是数组值从小到大的索引值
ax.plot(xSort[:,1], yHat[srtInd])
plt.savefig('k=0.01时局部加权线性回归示意图.png')
plt.show()
#-----------------------------------------------------------------------------
#3.----------使用岭回归---------------------------------------------------------
#岭回归
def ridgeRegres(xMat,yMat,lam=0.2):
xTx=xMat.T*xMat
denom=xTx+eye(shape(xMat)[1])*lam
if linalg.det(denom)==0.0:
print('This matrix is singular,connot do inverse')
return
ws=denom.I*(xMat.T*yMat.T)
return ws
def ridgeTest(xArr,yArr):
xMat=mat(xArr);yMat=mat(yArr)
yMean=mean(yMat.T,0)
yMat=yMat-yMean
xMeans=mean(xMat,0)
xVar=var(xMat,0)
xMat=(xMat-xMeans)/xVar
numTestPts=30
wMat=zeros((numTestPts,shape(xMat)[1]))
for i in range(numTestPts):
ws=ridgeRegres(xMat,yMat,exp(i-10))
wMat[i,:]=ws.T
return wMat
def plotRidge(ridgeWeights):
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.savefig('岭回归的回归系数变化图.png')
plt.show()
#-----------------------------------------------------------------------------
#-----------前向逐步回归-------------------------------------------------------
"""
伪代码：
数据标准化，使其分布满足0均值和单位方差
在每轮迭代过程中：
设置当前最小误差lowestError为正无穷
对每个特征：
增大或者缩小：
改变一个系数得到一个新的w
计算新的w下的误差
如果误差Error小于最小当前误差lowestError:
设置wbest等于当前的w
将w设置为新的wbest
"""
def stageWise(xArr,yArr,eps,numIt):
#数据标准化
xMat=mat(xArr);yMat=mat(yArr).T
yMean=mean(yMat,0)
yMat=yMat-yMean
xMeans=mean(xMat,0)
xVar=var(xMat,0)
xMat=(xMat-xMeans)/xVar
m,n=shape(xMat)
returnMat=zeros((numIt,n))
ws=zeros((n,1))
wsTest=ws.copy()
wsMax=ws.copy()
#每轮迭代过程
for i in range(numIt):
print(ws.T)
lowesError=inf
for j in range(n):
for sign in [-1,1]:
wsTest=ws.copy()
wsTest[j]+=eps*sign
yTest=xMat*wsTest
rssE=rssError(yMat.A,yTest.A)
if rssE<lowesError:
lowesError=rssE
wsMax=wsTest
ws=wsMax.copy()
returnMat[i,:]=ws.T
return returnMat
#分析预测误差大小
def rssError(yArr,yHatArr):
return((yArr-yHatArr)**2).sum()
def plotlasso(returnWeights):
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(returnWeights)
plt.savefig('前向逐步回归系数变化图.png')
plt.show()
#----------------------------------------------------------------------------
def main():
# #1.---------线性回归---------------------------
# xArr,yArr=loadDataSet('ex0.txt')
# ws=standRegres(xArr,yArr)
# plotScatter(xArr,yArr,ws)
# #---------------------------------------------
# #2.--------局部加权线性回归---------------------
# xArr,yArr=loadDataSet('ex0.txt')
# #ws=lwlr(xArr[0],xArr,yArr,0.001)
# #print(ws)
# yHat=lwlrTest(xArr,xArr,yArr,0.01)
# plotScatterlwlr(xArr,yArr,yHat)
# #---------------------------------------------
# #3.---------岭回归------------------------------
# abX,abY=loadDataSet('abalone.txt')
# ridgeWeights=ridgeTest(abX,abY)
# plotRidge(ridgeWeights)
#
# #----------------------------------------------
#4.---------前向逐步回归--------------------------
xArr,yArr=loadDataSet('abalone.txt')
returnMat=stageWise(xArr,yArr,0.005,1000)
plotlasso(returnMat)
print(returnMat)
#-----------------------------------------------
if __name__=='__main__':
main()
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
