• 《机器学习实战》8.预测数值型数据：回归

    目录

    预测数值型数据：回归

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

    2 局部加权线性回归

    3 示例：预测鲍鱼的年龄

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

    4.1 岭回归

    4.2 lasso

    4.3 前向逐步回归

    5 权衡偏差与方差

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

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

    6.2 训练算法：建立模型

    7 本章结束

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

    预测数值型数据：回归

    本章内容：

    ① 线性回归

    ② 局部加权线性回归

    ③ 岭回归和逐步线性回归

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

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

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

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

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

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

    回归的一般方法：

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

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

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

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

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

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

    采用最小平方误差

    1. # 标准回归函数和数据导入函数
    2. from numpy import *
    3. # 加载数据
    4. def loadDataSet(fileName):
    5.     # print(fileName)
    6.     numFeat=len(open(fileName).readline().split('\t'))-1
    7.     dataMat=[]
    8.     labelMat=[]
    9.     fr=open(fileName)
    10.     for line in fr.readlines():
    11.         lineArr=[]
    12.         curLine=line.strip().split('\t')
    13.         for i in range(numFeat):
    14.             lineArr.append(float(curLine[i]))
    15.             # print(lineArr)
    16.         # print(lineArr)
    17.         dataMat.append(lineArr)
    18.         labelMat.append(float(curLine[-1]))
    19.     # print(dataMat)
    20.     return dataMat,labelMat
    21.  
    22. # 线性回归函数
    23. def standRegres(xArr,yArr):
    24.     # 数组转为矩阵
    25.     xMat=mat(xArr)
    26.     yMat=mat(yArr).T
    27.     xTx=xMat.T*xMat
    28.     # linalg.det()函数可以直接用来计算行列式
    29.     if linalg.det(xTx)==0.0:
    30.         print("this matrix is singular,cannot do inverse")
    31.         return
    32.     # matrix.I获得与给定矩阵相同大小的乘法逆
    33.     ws=xTx.I*(xMat.T*yMat)
    34.     # 返回权重
    35.     return ws

    调用上述函数：

    1. xArr,yArr=loadDataSet('ex0.txt')
    2. # xArr
    3. # 回归系数
    4. ws=standRegres(xArr,yArr)
    5. ws

    输出得到的结果为：

     将得到的权重系数，画图表示为：

    1. # 绘图
    2. def plotData(xArr,yArr,ws):
    3.     import matplotlib.pyplot as plt
    4.     xMat = mat(xArr)
    5.     yMat = mat(yArr)
    6.     figure = plt.figure()
    7.     ax = figure.add_subplot(111)
    8.     # 取第二个特征绘图
    9.     # flatten()函数转化成一维数组
    10.     # matrix.A属性返回矩阵变成的数组，和getA()方法一样
    11.     # 绘制点
    12.     ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
    13.     # 返回给定数据的数组形式的拷贝
    14.     xCopy = xMat.copy()
    15.     # 升序排序
    16.     xCopy.sort(0)
    17.     print (ws.shape)
    18.     yHat = xCopy * ws # yHat 表示拟合直线的纵坐标，用回归系数求出
    19.     ax.plot(xCopy[:,1], yHat, c = 'green')
    20.     plt.show()
    21.  
    22.  
    23. plotData(xArr,yArr,ws)

    得到的输出图像为：

    可以看到该拟合，虽然在数据计算上结果较好，但可能单纯的线性拟合并不能够得到很好的结果，因为原始数据很明显有一个波浪形的分布，而不单单是线性分布。

     计算该拟合结果的相关系数为：

    1. # 计算相关系数
    2. def calcCorrel(xArr,yArr,ws):
    3.     xMat = mat(xArr)
    4.     yMat = mat(yArr)
    5.     yHat=xMat*ws
    6.     # 转置yHat向量，以保证都是行向量
    7.     # corrcoef函数计算皮尔逊相关系数
    8.     correl=corrcoef(yHat.T,yMat)
    9.     return correl
    10.  
    11. calcCorrel(xArr,yArr,ws)

     得到的结果为：

    2 局部加权线性回归

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

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

    1. def lwlr(testPoint,xArr,yArr,k=1.0):
    2.     xMat=mat(xArr)
    3.     yMat=mat(yArr).T
    4.     m=shape(xMat)[0]
    5.     # eye()函数返回一个二维数组，对角线为1，其余地方为0
    6.     weights=mat(eye((m)))
    7.     for j in range(m):
    8.         diffMat=testPoint-xMat[j,:]
    9.         weights[j,j]=exp(diffMat*diffMat.T/(-2.0*k**2))
    10.     xTx=xMat.T*(weights*xMat)
    11.     if linalg.det(xTx)==0.0:
    12.         print("this matrix is singular,cannot do inverse")
    13.         return
    14.     ws=xTx.I*(xMat.T*(weights*yMat))
    15.     
    16.     return testPoint*ws
    17.  
    18. def lwlrTest(testArr,xArr,yArr,k=1.0):
    19.     m=shape(testArr)[0]
    20.     yHat=zeros(m)
    21.     for i in range(m):
    22.         yHat[i]=lwlr(testArr[i],xArr,yArr,k)
    23.     return yHat

    调用上述函数：

    lwlr(xArr[0],xArr,yArr)

    得到的输出结果为：

     画出该拟合的结果曲线为：

    1. def plotDatalwlr(xArr,yArr,k=0.01):
    2.     yHat=lwlrTest(xArr,xArr,yArr,k)
    3.     xMat=mat(xArr)
    4.     # 对样本x升序排序并返回索引
    5.     srtInd=xMat[:,1].argsort(0)
    6.     xSort=xMat[srtInd][:,0,:]
    7.     import matplotlib.pyplot as plt
    8.     fig=plt.figure()
    9.     ax=fig.add_subplot(111)
    10.     ax.plot(xSort[:,1],yHat[srtInd])
    11.     ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='red')
    12.     plt.show()
    13. xArr,yArr=loadDataSet('ex0.txt')
    14. plotDatalwlr(xArr,yArr)

    得到的拟合图像为：

    局部加权线性回归也存在一个问题，即增加了计算量，因为他对每个点做预测时都必须使用整个数据集。如果避免这些计算将可以减少程序运行时间，从而减缓因计算量增加带来的问题。

    3 示例：预测鲍鱼的年龄

     

    1. # 计算误差
    2. def rssError(yArr,yHatArr):
    3.     # **表示幂运算，前面为底数，后面为指数
    4.     return ((yArr-yHatArr)**2).sum()
    5.  
    6. abX,abY=loadDataSet('abalone.txt')
    7.  
    8. plotDatalwlr(abX[0:99],abY[0:99],0.1)
    9. plotDatalwlr(abX[0:99],abY[0:99],1)
    10. plotDatalwlr(abX[0:99],abY[0:99],10)
    11.  
    12. # 训练集
    13. yHat01_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1)
    14. yHat1_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],1)
    15. yHat10_=lwlrTest(abX[0:99],abX[0:99],abY[0:99],10)
    16.  
    17. err01_=rssError(abY[0:99],yHat01_.T)
    18. err1_=rssError(abY[0:99],yHat1_.T)
    19. err10_=rssError(abY[0:99],yHat10_.T)
    20.  
    21. print(err01_,err1_,err10_)
    22.  
    23. # 测试集
    24. yHat01=lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
    25. yHat1=lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
    26. yHat10=lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)
    27.  
    28. err01=rssError(abY[100:199],yHat01.T)
    29. err1=rssError(abY[100:199],yHat1.T)
    30. err10=rssError(abY[100:199],yHat10.T)
    31.  
    32. print(err01,err1,err10)

    得到的输出结果为：

    可以看到，使用较小的核可以获得较低的误差，但是使用最小的核可能会造成过拟合，对新数据不一定能够达到最好的效果

    接下来再和简单的线性回归做个比较

    1. ws=standRegres(abX[0:99],abY[0:99])
    2. yHat=mat(abX[100:199])*ws
    3. rssError(abY[100:199],yHat.T.A)

     得到结果为：

    简单线性回归达到了与局部加权线性回归类似的效果，这也表明：必须在未知数据上比较效果才能够选取到最佳模型

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

    如果数据的特征比样本点还多，也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题，因此，引入了“岭回归”的概念，这就是第一种缩减方法，接着是lasso方法，效果很好但复杂，最后介绍了第二种缩减方法，称为前向逐步回归，可以得到与lasso差不多的效果，也更容易实现。

    4.1 岭回归

    岭回归就是在矩阵XTX上加入一个λI从而使矩阵非奇异，进而就能对矩阵求逆，其中矩阵I时一个单位矩阵，λ由用户自定义。

    缩减方法可以去掉不重要的参数，因此能够更好的理解数据，此外，与简单的线性回归相比，缩减法能够取得更好的效果。

    1. # 用于计算回归系数
    2. # 实现给定系数下的岭回归
    3. def ridgeRegres(xMat,yMat,lam=0.2):
    4.     xTx=xMat.T*xMat
    5.     denom=xTx+eye(shape(xMat)[1])*lam
    6.     if linalg.det(denom)==0.0:
    7.         print("this matrix is singlar,cannot do inverse")
    8.         return
    9.     ws=denom.I*(xMat.T*yMat)
    10.     return ws
    11.  
    12. # 在一组λ上测试结果
    13. def ridgeTest(xArr,yArr):
    14.     xMat=mat(xArr)
    15.     yMat=mat(yArr).T
    16.     yMean=mean(yMat,0)
    17.     yMat=yMat-yMean
    18.     xMeans=mean(xMat,0)
    19.     # var()函数求样本方差的无偏估计值，如果参数是1，就是有偏估计
    20.     xVar=var(xMat,0)
    21.     xMat=(xMat-xMeans)/xVar
    22.     numTestPts=30
    23.     wMat=zeros((numTestPts,shape(xMat)[1]))
    24.     for i in range(numTestPts):
    25.         ws=ridgeRegres(xMat,yMat,exp(i-10))
    26.         wMat[i,:]=ws.T
    27.     return wMat
    28.  
    29. abX,abY=loadDataSet('abalone.txt')
    30. ridgeWeights=ridgeTest(abX,abY)
    31.  
    32. # 绘制出回归系数与log（λ）的关系
    33. import matplotlib.pyplot as plt
    34. fig=plt.figure()
    35. ax=fig.add_subplot(111)
    36. ax.plot(ridgeWeights)
    37. plt.show()

     得到的结果为：

     

    4.2 lasso

    lasso对回归系数做出了限定，使用绝对值取代了平方和，在λ足够小的时候，一些系数会因此被迫缩减到0，这个特性将有助于理解数据

    4.3 前向逐步回归

    属于一种贪心的算法，每一步都尽可能的减小误差。一开始所有的权重都设置为1，然后每一步所作的决策是对某个权重增加或减少一个很小的值

    伪代码：

    数据标准化，使其分布满足0均值和单位方差

    在每轮迭代中：

        设置当前最小误差lowestError为正无穷

        对每个特征：

            增大或减小：

                改变一个系数得到一个新的w

                计算新w下的误差

                如果误差Error小于当前最小误差lowestError：

                    设置Wbest等于当前w

            将w设置为新的Wbest

    1. # 实现矩阵归一化
    2. def regularize(xMat):
    3.     inMat = xMat.copy()
    4.     # 得到平均数，压缩行，对每一列求平均值
    5.     inMeans = mean(inMat,0)  
    6.     
    7.     inVar = var(inMat,0)    
    8.     inMat = (inMat - inMeans)/inVar
    9.     return inMat
    10.  
    11. # 前向逐步线性回归
    12. def stageWise(xArr,yArr,eps=0.01,numIt=100):
    13.     xMat=mat(xArr)
    14.     yMat=mat(yArr).T
    15.     # mean()函数求平均值
    16.     yMean=mean(yMat,0)
    17.     yMat=yMat-yMean
    18.     xMat=regularize(xMat)
    19.     m,n=shape(xMat)
    20.     returnMat=zeros((numIt,n))
    21.     ws=zeros((n,1))
    22.     wsTest=ws.copy()
    23.     wsMax=ws.copy()
    24.     for i in range(numIt):
    25.         print(ws.T)
    26.         lowestError=inf
    27.         for j in range(n):
    28.             for sign in [-1,1]:
    29.                 wsTest=ws.copy()
    30.                 wsTest[j]+=eps*sign
    31.                 yTest=xMat*wsTest
    32.                 rssE=rssError(yMat.A,yTest.A)
    33.                 if rssE
    34.                     lowestError=rssE
    35.                     wsMax=wsTest
    36.         ws=wsMax.copy()
    37.         returnMat[i,:]=ws.T
    38.     return returnMat

     调用函数：

    1. xArr,yArr=loadDataSet('abalone.txt')
    2. stageWise(xArr,yArr,0.01,200)

    得到的输出结果为：

     与最小二乘法的结果进行对比得到

    1. # 最小二乘法的结果
    2. xMat=mat(xArr)
    3. yMat=mat(yArr).T
    4. xMat=regularize(xMat)
    5. yM=mean(yMat,0)
    6. yMat=yMat-yM
    7. weights=standRegres(xMat,yMat.T)
    8. weights.T

    结果为：

     

    5 权衡偏差与方差

    任何时候，一旦发现模型和测量值直接存在差异，就说明出现了误差。当考虑模型中的“噪声”或者说误差时，必须考虑其来源。可能会因为对复杂的模型进行简化，这将导致在模型和测量值之间出现“噪声”或者“误差”，若无法理解数据的真实生成过程，也会导致差异的发生。

    一般认为，误差一般由三个部分组成：偏差、测量误差和随机噪声

    方差是可以度量的，如果从总体数据中取一个随机样本集并用线性模型拟合，将会得到一组回归系数，同理，再取出另一组随机样本集并拟合，将会得到另一组回归系数。这些系数之间的差异大小也就是模型反差大小的反映。

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

    ①收集数据：用google Shopping的APi收集数据

    ②准备数据：从返回的json数据中抽取价格

    ③分析数据：可视化并观察数据

    ④训练算法：构建不同的模型，采用逐步线性回归和直线的线性回归模型

    ⑤测试算法：使用交叉验证来测试不同的模型，分析哪个效果更好

    ⑥使用算法：这次联系的目的就是生成数据模型

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

    1. from bs4 import BeautifulSoup
    2. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
    3.      # 打开并读取HTML文件
    4.     with open(inFile, encoding='utf-8') as f:
    5.         html = f.read()
    6.     # beautiful函数将复杂HTML文档转换成一个复杂的树形结构
    7.     soup = BeautifulSoup(html)
    8.     i = 1
    9.     # 根据HTML页面结构进行解析
    10.     currentRow = soup.find_all('table', r = "%d" % i)
    11.     while(len(currentRow) != 0):
    12.         currentRow = soup.find_all('table', r = "%d" % i)
    13.         title = currentRow[0].find_all('a')[1].text
    14.         lwrTitle = title.lower()
    15.         # 查找是否有全新标签
    16.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
    17.             newFlag = 1.0
    18.         else:
    19.             newFlag = 0.0
    20.         # 查找是否已经标志出售，我们只收集已出售的数据
    21.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
    22.         if len(soldUnicde) == 0:
    23.             print("商品 #%d 没有出售" % i)
    24.         else:
    25.             # 解析页面获取当前价格
    26.             soldPrice = currentRow[0].find_all('td')[4]
    27.             priceStr = soldPrice.text
    28.             priceStr = priceStr.replace('$','')
    29.             priceStr = priceStr.replace(',','')
    30.             if len(soldPrice) > 1:
    31.                 priceStr = priceStr.replace('Free shipping', '')
    32.             sellingPrice = float(priceStr)
    33.             # 去掉不完整的套装价格
    34.             if  sellingPrice > origPrc * 0.5:
    35.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
    36.                 retX.append([yr, numPce, newFlag, origPrc])
    37.                 retY.append(sellingPrice)
    38.         i += 1
    39.         currentRow = soup.find_all('table', r = "%d" % i)
    40.  
    41. def setDataCollect(retX, retY):
    42.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)  # 2006年的乐高8288,部件数目800,原价49.99
    43.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)  # 2002年的乐高10030,部件数目3096,原价269.99
    44.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)  # 2007年的乐高10179,部件数目5195,原价499.99
    45.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)  # 2007年的乐高10181,部件数目3428,原价199.99
    46.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)  # 2008年的乐高10189,部件数目5922,原价299.99
    47.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)  # 2009年的乐高10196,部件数目3263,原价249.99

    调用上面函数：

    1. lgX=[]
    2. lgY=[]
    3. setDataCollect(lgX,lgY)

    得到结果：

    6.2 训练算法：建立模型

    为上面收集到的数据建立模型。构建出来的模型可以对售价做出预测

    1. # 创建一个全为1 的矩阵
    2. lgX1=mat(ones((63,5)))
    3. print(lgX1[0])
    4. # 将原数据矩阵lgX复制到新数据矩阵lgX1的第一列到第五列
    5. lgX1[:,1:5]=mat(lgX)
    6. # 确认数据复制的正确性
    7. print(lgX[0])
    8. print(lgX1[0])
    9. ws=standRegres(lgX1,lgY)
    10. ws

    得到的输出结果为：

     交叉测试验证岭回归

    1. # 交叉测试验证岭回归
    2.                             #  交叉验证的次数
    3. def crossValidation(xArr,yArr,numVal=10):
    4.     m=len(yArr)
    5.     indexList=list(range(m))
    6.     errorMat=zeros((numVal,30))
    7.     for i in range(numVal):
    8.         trainX=[]
    9.         trainY=[]
    10.         testX=[]
    11.         testY=[]
    12.         random.shuffle(indexList)
    13.         # 将数据分为测试集和训练集
    14.         for j in range(m):
    15.             if j0.9:
    16.                 trainX.append(xArr[indexList[j]])
    17.                 trainY.append(yArr[indexList[j]])
    18.             else:
    19.                 testX.append(xArr[indexList[j]])
    20.                 testY.append(yArr[indexList[j]])
    21.             
    22.         wMat=ridgeTest(trainX,trainY)
    23.         for k in range(30):
    24.             matTestX=mat(testX)
    25.             matTrainX=mat(trainX)
    26.             meanTrain=mean(matTrainX)
    27.             varTrain=var(matTrainX,0)
    28.             matTeatX=(matTestX-meanTrain)/varTrain
    29.             yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
    30.             errorMat[i,k]=rssError(yEst.T.A,array(testY))
    31.     meanErrors=mean(errorMat,0)
    32.     minMean=float(min(meanErrors))
    33.     bestWeights=wMat[nonzero(meanErrors==minMean)]
    34.     xMat=mat(xArr)
    35.     yMat=mat(yArr).T
    36.     meanX=mean(xMat,0)
    37.     varX=var(xMat,0)
    38.     unReg=bestWeights/varX
    39.     # 岭回归的最佳模型是
    40.     print("the best model from Ridge Regression is :\n",unReg)
    41.     # 常数项
    42.     print("with constant term:",-1*sum(multiply(meanX,unReg))+mean(yMat))
    43.  
    44. crossValidation(lgX,lgY,10)

     得到的输出结果为：

    这些系数是经过不同程度的缩减得到的

    7 本章结束

  • 相关阅读:
    Oceanbase读写分离方案探索与优化
    [笔记] 线性筛质数（素数）
    golang 使用 make 创建 map 是否需要指定长度
    线程 Pthread API
    【构建ML驱动的应用程序】第 1 章 ：从产品目标到 ML 框架
    pytorh变量的定义和访问，detatch（）作用
    Elasticsearch Index Monitoring(索引监控)之Index Stats API详解
    【深度学习】实验4布置：脑部 MRI 图像分割
    Vue 2.0中引入的类型检查Flow
    对IP协议概念以及IP地址的概念进行简单整理
  • 原文地址：https://blog.csdn.net/weixin_46182244/article/details/127908310