目录

预测数值型数据：回归

1 利用线性回归找到最佳拟合直线

2 局部加权线性回归

3 示例：预测鲍鱼的年龄

4 缩减系数来“理解”数据

4.1 岭回归

4.2 lasso

4.3 前向逐步回归

5 权衡偏差与方差

6 示例：预测乐高玩具套装的价格

6.1 收集数据：使用google购物的API

6.2 训练算法：建立模型

7 本章结束

本节涉及的相关代码和数据

预测数值型数据：回归

本章内容：

① 线性回归

② 局部加权线性回归

③ 岭回归和逐步线性回归

④ 预测鲍鱼年龄和玩具售价

1 利用线性回归找到最佳拟合直线

线性回归的相关特点主要是

优点：结果易于理解，计算上不复杂

缺点：对非线性数据不友好

适用数据类型：数值型和标称型

回归的一般方法：

① 收集数据：采用任意方法收集数据

② 准备数据：回归需要数值型数据，标称型数据将会被转为二值型数据

③ 分析数据：绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合曲线绘在图上作为对比

④ 训练算法：找到回归系数

⑤ 测试算法：使用R2或者预测值和数据的拟合度，来分析模型的结果

⑥ 使用算法：使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续性数据而不仅仅是离散的数据标签

采用最小平方误差

# 标准回归函数和数据导入函数
from numpy import *
# 加载数据
def loadDataSet(fileName):# print(fileName)numFeat=len(open(fileName).readline().split('\t'))-1dataMat=[]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).TxTx=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 pltxMat = 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 correlcalcCorrel(xArr,yArr,ws)

 得到的结果为：

2 局部加权线性回归

线性回归的一个问题就是有可能出现欠拟合的现象，因为他求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。

其中一个方法就是局部加权线性回归（LWLR）。我们给待预测点附近的每个点赋予一定的权重，然后在这个基础上基于最小均方差来进行普通的回归。与KNN一样，这种算法每次预测均需事先选取处对应的数据子集。LWLR使用‘核’来对附近的点赋予更高的权重。核的类型可以自由选择，常用的有高斯核。

def lwlr(testPoint,xArr,yArr,k=1.0):xMat=mat(xArr)yMat=mat(yArr).Tm=shape(xMat)[0]# eye()函数返回一个二维数组，对角线为1，其余地方为0weights=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")returnws=xTx.I*(xMat.T*(weights*yMat))return testPoint*wsdef 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 pltfig=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*xMatdenom=xTx+eye(shape(xMat)[1])*lamif linalg.det(denom)==0.0:print("this matrix is singlar,cannot do inverse")returnws=denom.I*(xMat.T*yMat)return ws# 在一组λ上测试结果
def ridgeTest(xArr,yArr):xMat=mat(xArr)yMat=mat(yArr).TyMean=mean(yMat,0)yMat=yMat-yMeanxMeans=mean(xMat,0)# var()函数求样本方差的无偏估计值，如果参数是1，就是有偏估计xVar=var(xMat,0)xMat=(xMat-xMeans)/xVarnumTestPts=30wMat=zeros((numTestPts,shape(xMat)[1]))for i in range(numTestPts):ws=ridgeRegres(xMat,yMat,exp(i-10))wMat[i,:]=ws.Treturn wMatabX,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)/inVarreturn inMat# 前向逐步线性回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):xMat=mat(xArr)yMat=mat(yArr).T# mean()函数求平均值yMean=mean(yMat,0)yMat=yMat-yMeanxMat=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=inffor j in range(n):for sign in [-1,1]:wsTest=ws.copy()wsTest[j]+=eps*signyTest=xMat*wsTestrssE=rssError(yMat.A,yTest.A)if rssE

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

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].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else: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.textpriceStr = priceStr.replace('$','')priceStr = priceStr.replace(',','')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if  sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = 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.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)  # 2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)  # 2007年的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)  # 2007年的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)  # 2008年的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)  # 2009年的乐高10196,部件数目3263,原价249.99

lgX=[]
lgY=[]
setDataCollect(lgX,lgY)

# 创建一个全为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