《机器学习实战》8.预测数值型数据：回归
目录
本节涉及的相关代码和数据
预测数值型数据：回归
本章内容：
① 线性回归
② 局部加权线性回归
③ 岭回归和逐步线性回归
④ 预测鲍鱼年龄和玩具售价
1 利用线性回归找到最佳拟合直线
线性回归的相关特点主要是
优点：结果易于理解，计算上不复杂
缺点：对非线性数据不友好
适用数据类型：数值型和标称型
回归的一般方法：
① 收集数据：采用任意方法收集数据
② 准备数据：回归需要数值型数据，标称型数据将会被转为二值型数据
③ 分析数据：绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合曲线绘在图上作为对比
④ 训练算法：找到回归系数
⑤ 测试算法：使用R2或者预测值和数据的拟合度，来分析模型的结果
⑥ 使用算法：使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续性数据而不仅仅是离散的数据标签
采用最小平方误差
# 标准回归函数和数据导入函数
from numpy import *
# 加载数据
def loadDataSet(fileName):
# print(fileName)
numFeat=len(open(fileName).readline().split('t'))-1
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]))
# print(lineArr)
# print(lineArr)
dataMat.append(lineArr)
labelMat.append(float(curLine[-1]))
# print(dataMat)
return dataMat,labelMat
# 线性回归函数
def standRegres(xArr,yArr):
# 数组转为矩阵
xMat=mat(xArr)
yMat=mat(yArr).T
xTx=xMat.T*xMat
# linalg.det()函数可以直接用来计算行列式
if linalg.det(xTx)==0.0:
print("this matrix is singular,cannot do inverse")
return
# matrix.I获得与给定矩阵相同大小的乘法逆
ws=xTx.I*(xMat.T*yMat)
# 返回权重
return ws调用上述函数：
xArr,yArr=loadDataSet('ex0.txt')
# xArr
# 回归系数
ws=standRegres(xArr,yArr)
ws
输出得到的结果为：
将得到的权重系数，画图表示为：
# 绘图
def plotData(xArr,yArr,ws):
import matplotlib.pyplot as plt
xMat = mat(xArr)
yMat = mat(yArr)
figure = plt.figure()
ax = figure.add_subplot(111)
# 取第二个特征绘图
# flatten()函数转化成一维数组
# matrix.A属性返回矩阵变成的数组，和getA()方法一样
# 绘制点
ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
# 返回给定数据的数组形式的拷贝
xCopy = xMat.copy()
# 升序排序
xCopy.sort(0)
print (ws.shape)
yHat = xCopy * ws # yHat 表示拟合直线的纵坐标，用回归系数求出
ax.plot(xCopy[:,1], yHat, c = 'green')
plt.show()
plotData(xArr,yArr,ws)得到的输出图像为：
可以看到该拟合，虽然在数据计算上结果较好，但可能单纯的线性拟合并不能够得到很好的结果，因为原始数据很明显有一个波浪形的分布，而不单单是线性分布。
计算该拟合结果的相关系数为：
# 计算相关系数
def calcCorrel(xArr,yArr,ws):
xMat = mat(xArr)
yMat = mat(yArr)
yHat=xMat*ws
# 转置yHat向量，以保证都是行向量
# corrcoef函数计算皮尔逊相关系数
correl=corrcoef(yHat.T,yMat)
return correl
calcCorrel(xArr,yArr,ws)得到的结果为：
2 局部加权线性回归
线性回归的一个问题就是有可能出现欠拟合的现象，因为他求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。
其中一个方法就是局部加权线性回归（LWLR）。我们给待预测点附近的每个点赋予一定的权重，然后在这个基础上基于最小均方差来进行普通的回归。与KNN一样，这种算法每次预测均需事先选取处对应的数据子集。LWLR使用‘核’来对附近的点赋予更高的权重。核的类型可以自由选择，常用的有高斯核。
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat=mat(xArr)
yMat=mat(yArr).T
m=shape(xMat)[0]
# eye()函数返回一个二维数组，对角线为1，其余地方为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))
return testPoint*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)
return yHat调用上述函数：
lwlr(xArr[0],xArr,yArr)得到的输出结果为：
画出该拟合的结果曲线为：
def plotDatalwlr(xArr,yArr,k=0.01):
yHat=lwlrTest(xArr,xArr,yArr,k)
xMat=mat(xArr)
# 对样本x升序排序并返回索引
srtInd=xMat[:,1].argsort(0)
xSort=xMat[srtInd][:,0,:]
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='red')
plt.show()
xArr,yArr=loadDataSet('ex0.txt')
plotDatalwlr(xArr,yArr)得到的拟合图像为：
局部加权线性回归也存在一个问题，即增加了计算量，因为他对每个点做预测时都必须使用整个数据集。如果避免这些计算将可以减少程序运行时间，从而减缓因计算量增加带来的问题。
3 示例：预测鲍鱼的年龄
# 计算误差
def rssError(yArr,yHatArr):
# **表示幂运算，前面为底数，后面为指数
return ((yArr-yHatArr)**2).sum()
abX,abY=loadDataSet('abalone.txt')
plotDatalwlr(abX[0:99],abY[0:99],0.1)
plotDatalwlr(abX[0:99],abY[0:99],1)
plotDatalwlr(abX[0:99],abY[0:99],10)
# 训练集
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)
err01_=rssError(abY[0:99],yHat01_.T)
err1_=rssError(abY[0:99],yHat1_.T)
err10_=rssError(abY[0:99],yHat10_.T)
print(err01_,err1_,err10_)
# 测试集
yHat01=lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
yHat1=lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
yHat10=lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)
err01=rssError(abY[100:199],yHat01.T)
err1=rssError(abY[100:199],yHat1.T)
err10=rssError(abY[100:199],yHat10.T)
print(err01,err1,err10)得到的输出结果为：
可以看到，使用较小的核可以获得较低的误差，但是使用最小的核可能会造成过拟合，对新数据不一定能够达到最好的效果
接下来再和简单的线性回归做个比较
ws=standRegres(abX[0:99],abY[0:99])
yHat=mat(abX[100:199])*ws
rssError(abY[100:199],yHat.T.A)得到结果为：
简单线性回归达到了与局部加权线性回归类似的效果，这也表明：必须在未知数据上比较效果才能够选取到最佳模型
4 缩减系数来“理解”数据
如果数据的特征比样本点还多，也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题，因此，引入了“岭回归”的概念，这就是第一种缩减方法，接着是lasso方法，效果很好但复杂，最后介绍了第二种缩减方法，称为前向逐步回归，可以得到与lasso差不多的效果，也更容易实现。
4.1 岭回归
岭回归就是在矩阵XTX上加入一个λI从而使矩阵非奇异，进而就能对矩阵求逆，其中矩阵I时一个单位矩阵，λ由用户自定义。
缩减方法可以去掉不重要的参数，因此能够更好的理解数据，此外，与简单的线性回归相比，缩减法能够取得更好的效果。
# 用于计算回归系数
# 实现给定系数下的岭回归
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 singlar,cannot do inverse")
return
ws=denom.I*(xMat.T*yMat)
return ws
# 在一组λ上测试结果
def ridgeTest(xArr,yArr):
xMat=mat(xArr)
yMat=mat(yArr).T
yMean=mean(yMat,0)
yMat=yMat-yMean
xMeans=mean(xMat,0)
# var()函数求样本方差的无偏估计值，如果参数是1，就是有偏估计
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
abX,abY=loadDataSet('abalone.txt')
ridgeWeights=ridgeTest(abX,abY)
# 绘制出回归系数与log（λ）的关系
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()得到的结果为：
4.2 lasso
lasso对回归系数做出了限定，使用绝对值取代了平方和，在λ足够小的时候，一些系数会因此被迫缩减到0，这个特性将有助于理解数据
4.3 前向逐步回归
属于一种贪心的算法，每一步都尽可能的减小误差。一开始所有的权重都设置为1，然后每一步所作的决策是对某个权重增加或减少一个很小的值
伪代码：
数据标准化，使其分布满足0均值和单位方差
在每轮迭代中：
设置当前最小误差lowestError为正无穷
对每个特征：
增大或减小：
改变一个系数得到一个新的w
计算新w下的误差
如果误差Error小于当前最小误差lowestError：
设置Wbest等于当前w
将w设置为新的Wbest
# 实现矩阵归一化
def regularize(xMat):
inMat = xMat.copy()
# 得到平均数，压缩行，对每一列求平均值
inMeans = mean(inMat,0)
inVar = var(inMat,0)
inMat = (inMat - inMeans)/inVar
return inMat
# 前向逐步线性回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):
xMat=mat(xArr)
yMat=mat(yArr).T
# mean()函数求平均值
yMean=mean(yMat,0)
yMat=yMat-yMean
xMat=regularize(xMat)
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)
lowestError=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<lowestError:
lowestError=rssE
wsMax=wsTest
ws=wsMax.copy()
returnMat[i,:]=ws.T
return returnMat调用函数：
xArr,yArr=loadDataSet('abalone.txt')
stageWise(xArr,yArr,0.01,200)得到的输出结果为：
与最小二乘法的结果进行对比得到
# 最小二乘法的结果
xMat=mat(xArr)
yMat=mat(yArr).T
xMat=regularize(xMat)
yM=mean(yMat,0)
yMat=yMat-yM
weights=standRegres(xMat,yMat.T)
weights.T结果为：
5 权衡偏差与方差
任何时候，一旦发现模型和测量值直接存在差异，就说明出现了误差。当考虑模型中的“噪声”或者说误差时，必须考虑其来源。可能会因为对复杂的模型进行简化，这将导致在模型和测量值之间出现“噪声”或者“误差”，若无法理解数据的真实生成过程，也会导致差异的发生。
一般认为，误差一般由三个部分组成：偏差、测量误差和随机噪声
方差是可以度量的，如果从总体数据中取一个随机样本集并用线性模型拟合，将会得到一组回归系数，同理，再取出另一组随机样本集并拟合，将会得到另一组回归系数。这些系数之间的差异大小也就是模型反差大小的反映。
6 示例：预测乐高玩具套装的价格
①收集数据：用google Shopping的APi收集数据
②准备数据：从返回的json数据中抽取价格
③分析数据：可视化并观察数据
④训练算法：构建不同的模型，采用逐步线性回归和直线的线性回归模型
⑤测试算法：使用交叉验证来测试不同的模型，分析哪个效果更好
⑥使用算法：这次联系的目的就是生成数据模型
6.1 收集数据：使用google购物的API
from bs4 import BeautifulSoup
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
# 打开并读取HTML文件
with open(inFile, encoding='utf-8') as f:
html = f.read()
# beautiful函数将复杂HTML文档转换成一个复杂的树形结构
soup = BeautifulSoup(html)
i = 1
# 根据HTML页面结构进行解析
currentRow = soup.find_all('table', r = "%d" % i)
while(len(currentRow) != 0):
currentRow = soup.find_all('table', r = "%d" % i)
title = currentRow[0].find_all('a')[1].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
newFlag = 0.0
# 查找是否已经标志出售，我们只收集已出售的数据
soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
if len(soldUnicde) == 0:
print("商品 #%d 没有出售" % i)
else:
# 解析页面获取当前价格
soldPrice = currentRow[0].find_all('td')[4]
priceStr = soldPrice.text
priceStr = priceStr.replace('$','')
priceStr = priceStr.replace(',','')
if len(soldPrice) > 1:
priceStr = priceStr.replace('Free shipping', '')
sellingPrice = float(priceStr)
# 去掉不完整的套装价格
if sellingPrice > origPrc * 0.5:
print("%dt%dt%dt%ft%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
retX.append([yr, numPce, newFlag, origPrc])
retY.append(sellingPrice)
i += 1
currentRow = soup.find_all('table', r = "%d" % i)
def setDataCollect(retX, retY):
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) # 2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) # 2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) # 2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) # 2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) # 2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) # 2009年的乐高10196,部件数目3263,原价249.99调用上面函数：
lgX=[]
lgY=[]
setDataCollect(lgX,lgY)得到结果：
6.2 训练算法：建立模型
为上面收集到的数据建立模型。构建出来的模型可以对售价做出预测
# 创建一个全为1 的矩阵
lgX1=mat(ones((63,5)))
print(lgX1[0])
# 将原数据矩阵lgX复制到新数据矩阵lgX1的第一列到第五列
lgX1[:,1:5]=mat(lgX)
# 确认数据复制的正确性
print(lgX[0])
print(lgX1[0])
ws=standRegres(lgX1,lgY)
ws得到的输出结果为：
交叉测试验证岭回归
# 交叉测试验证岭回归
# 交叉验证的次数
def crossValidation(xArr,yArr,numVal=10):
m=len(yArr)
indexList=list(range(m))
errorMat=zeros((numVal,30))
for i in range(numVal):
trainX=[]
trainY=[]
testX=[]
testY=[]
random.shuffle(indexList)
# 将数据分为测试集和训练集
for j in range(m):
if j<m*0.9:
trainX.append(xArr[indexList[j]])
trainY.append(yArr[indexList[j]])
else:
testX.append(xArr[indexList[j]])
testY.append(yArr[indexList[j]])
wMat=ridgeTest(trainX,trainY)
for k in range(30):
matTestX=mat(testX)
matTrainX=mat(trainX)
meanTrain=mean(matTrainX)
varTrain=var(matTrainX,0)
matTeatX=(matTestX-meanTrain)/varTrain
yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
errorMat[i,k]=rssError(yEst.T.A,array(testY))
meanErrors=mean(errorMat,0)
minMean=float(min(meanErrors))
bestWeights=wMat[nonzero(meanErrors==minMean)]
xMat=mat(xArr)
yMat=mat(yArr).T
meanX=mean(xMat,0)
varX=var(xMat,0)
unReg=bestWeights/varX
# 岭回归的最佳模型是
print("the best model from Ridge Regression is :n",unReg)
# 常数项
print("with constant term:",-1*sum(multiply(meanX,unReg))+mean(yMat))
crossValidation(lgX,lgY,10)得到的输出结果为：
这些系数是经过不同程度的缩减得到的