一、前言

本文将会对连续型的数据做出预测。主要讲解简单的线性回归和局部加权线性回归，并通过预测鲍鱼年龄和乐高玩具套件的二手价格的实例进行实战演练。

二、什么是回归？

回归的目的是预测数值型的目标值。最直接的办法是依据输入写出一个目标值的计算公式。假如你想预测小姐姐男友汽车的功率，可能会这么计算：
image.png
写成中文就是：
image.png
这就是所谓的回归方程（regression equation），其中的0.0015和-0.99称为回归系数（regression weights），求这些回归系数的过程就是回归。一旦有了这些回归系数，再给定输入，做预测就非常容易了。具体的做法是用回归系数乘以输入值，再将结果全部加在一起，就得到了预测值。

说到回归，一般都是指线性回归（linear regression），所以本文里的回归和线性回归代表同一个意思。线性回归意味着可以将输入项分别乘以一些常量，再将结果加起来得到输出。需要说明的是，存在另一种成为非线性回归的回归模型，该模型不认同上面的做法，比如认为输出可能是输入的乘积。这样，上面的功率计算公式也可以写做：
image.png
这就是一个非线性回归的例子，本文对此不做深入讨论。

三、揭开回归的神秘面纱

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

应该怎么从一大堆数据里求出回归方程呢？假定输入数据存放在矩阵X中，结果存放在向量y中：

1.png

而回归系数存放在向量w中：

2.png

那么对于给定的数据x1，即矩阵X的第一列数据，预测结果u1将会通过如下公式给出：

3.png

现在的问题是，手里有数据矩阵X和对应的标签向量y，怎么才能找到w呢？一个常用的方法就是找出使误差最小的w。这里的误差是指预测u值和真实y值之间的差值，使用该误差的简单累加将使得正差值和负差值相互抵消，所以我们采用平方误差。

平方误差和可以写做：

4.png

用矩阵表示还可以写做：

5.png

为啥能这么变化，记住一个前提：若x为向量，则默认x为列向量，x^T为行向量。将上述提到的数据矩阵X和标签向量y带进去，就知道为何这么变化了。

在继续推导之前，我们要先明确一个目的：找到w，使平方误差和最小。因为我们认为平方误差和越小，说明线性回归拟合效果越好。

现在，我们用矩阵表示的平方误差和对w进行求导：

6.png

令上述公式等于0，得到：

7.png

w上方的小标记表示，这是当前可以估计出的w的最优解。从现有数据上估计出的w可能并不是数据中的真实w值，所以这里使用了一个”帽”符号来表示它仅是w的一个最佳估计。

值得注意的是，上述公式中包含逆矩阵，也就是说，这个方程只在逆矩阵存在的时候使用，也即是这个矩阵是一个方阵，并且其行列式不为0。

述的最佳w求解是统计学中的常见问题，除了矩阵方法外还有很多其他方法可以解决。通过调用NumPy库里的矩阵方法，我们可以仅使用几行代码就完成所需功能。该方法也称作OLS， 意思是“普通小二乘法”（ordinary least squares）。

我们先看下数据集，数据下载地址：数据集下载

8.png

第一列都为1.0，即x0。第二列为x1，即x轴数据。第三列为x2，即y轴数据。首先绘制下数据，看下数据分布。编写代码如下：

  1. # -*- coding:utf-8 -*-
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  5. def loadDataSet(fileName):
  6.     """
  7.     函数说明:加载数据
  8.     Parameters:
  9.         fileName - 文件名
  10.     Returns:
  11.         xArr - x数据集
  12.         yArr - y数据集
  14.     """
  16.     numFeat = len(open(fileName).readline().split('\t')) - 1
  17.     xArr = []; yArr = []
  18.     fr = open(fileName)
  19.     for line in fr.readlines():
  20.         lineArr =[]
  21.         curLine = line.strip().split('\t')
  22.         for i in range(numFeat):
  23.             lineArr.append(float(curLine[i]))
  24.         xArr.append(lineArr)
  25.         yArr.append(float(curLine[-1]))
  26.     return xArr, yArr
  28. def plotDataSet():
  29.     """
  30.     函数说明:绘制数据集
  31.     Parameters:
  32.         无
  33.     Returns:
  34.         无
  36.     """
  37.     xArr, yArr = loadDataSet('ex0.txt')                                    #加载数据集
  38.     n = len(xArr)                                                        #数据个数
  39.     xcord = []; ycord = []                                                #样本点
  40.     for i in range(n):                                                   
  41.         xcord.append(xArr[i][1]); ycord.append(yArr[i])                    #样本点
  42.     fig = plt.figure()
  43.     ax = fig.add_subplot(111)                                            #添加subplot
  44.     ax.scatter(xcord, ycord, s = 20, c = 'blue',alpha = .5)                #绘制样本点
  45.     plt.title('DataSet')                                                #绘制title
  46.     plt.xlabel('X')
  47.     plt.show()
  49. if __name__ == '__main__':
  50.     plotDataSet()

运行代码如下：

9.png

通过可视化数据，我们可以看到数据的分布情况。接下来，让我们根据上文中推导的回归系数计算方法，求出回归系数向量，并根据回归系数向量绘制回归曲线，编写代码如下：

  1. # -*- coding:utf-8 -*-
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  5. def loadDataSet(fileName):
  6.     """
  7.     函数说明:加载数据
  8.     Parameters:
  9.         fileName - 文件名
  10.     Returns:
  11.         xArr - x数据集
  12.         yArr - y数据集
  14.     """
  15.     numFeat = len(open(fileName).readline().split('\t')) - 1
  16.     xArr = []; yArr = []
  17.     fr = open(fileName)
  18.     for line in fr.readlines():
  19.         lineArr =[]
  20.         curLine = line.strip().split('\t')
  21.         for i in range(numFeat):
  22.             lineArr.append(float(curLine[i]))
  23.         xArr.append(lineArr)
  24.         yArr.append(float(curLine[-1]))
  25.     return xArr, yArr
  27. def standRegres(xArr,yArr):
  28.     """
  29.     函数说明:计算回归系数w
  30.     Parameters:
  31.         xArr - x数据集
  32.         yArr - y数据集
  33.     Returns:
  34.         ws - 回归系数
  36.     """
  37.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  38.     xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
  39.     if np.linalg.det(xTx) == 0.0:
  40.         print("矩阵为奇异矩阵,不能求逆")
  41.         return
  42.     ws = xTx.I * (xMat.T*yMat)
  43.     return ws
  45. def plotRegression():
  46.     """
  47.     函数说明:绘制回归曲线和数据点
  48.     Parameters:
  49.         无
  50.     Returns:
  51.         无
  53.     """
  54.     xArr, yArr = loadDataSet('ex0.txt')                                    #加载数据集
  55.     ws = standRegres(xArr, yArr)                                        #计算回归系数
  56.     xMat = np.mat(xArr)                                                    #创建xMat矩阵
  57.     yMat = np.mat(yArr)                                                    #创建yMat矩阵
  58.     xCopy = xMat.copy()                                                    #深拷贝xMat矩阵
  59.     xCopy.sort(0)                                                        #排序
  60.     yHat = xCopy * ws                                                     #计算对应的y值
  61.     fig = plt.figure()
  62.     ax = fig.add_subplot(111)                                            #添加subplot
  63.     ax.plot(xCopy[:, 1], yHat, c = 'red')                                #绘制回归曲线
  64.     ax.scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue',alpha = .5)                #绘制样本点
  65.     plt.title('DataSet')                                                #绘制title
  66.     plt.xlabel('X')
  67.     plt.show()
  69. if __name__ == '__main__':
  70.     plotRegression()

运行结果如下：

10.png

可以看到，对角线上的数据是1.0，因为yMat和自己的匹配是完美的，而YHat和yMat的相关系数为0.98。

最佳拟合直线方法将数据视为直线进行建模，具有十分不错的表现。数据当中似乎还存在其他的潜在模式。那么如何才能利用这些模式呢？我们可以根据数据来局部调整预测，下面就会介绍这种方法。

2、局部加权线性回归

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

其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该方法中，我们给待预测点附近的每个点赋予一定的权重。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解除回归系数W的形式如下：

11.png

其中W是一个矩阵，这个公式跟我们上面推导的公式的区别就在于W，它用来给每个点赋予权重。

LWLR使用”核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

12.png

这样我们就可以根据上述公式，编写局部加权线性回归，我们通过改变k的值，可以调节回归效果，编写代码如下：

  1. # -*- coding:utf-8 -*-
  2. from matplotlib.font_manager import FontProperties
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  5. def loadDataSet(fileName):
  6.     """
  7.     函数说明:加载数据
  8.     Parameters:
  9.         fileName - 文件名
  10.     Returns:
  11.         xArr - x数据集
  12.         yArr - y数据集
  14.     """
  15.     numFeat = len(open(fileName).readline().split('\t')) - 1
  16.     xArr = []; yArr = []
  17.     fr = open(fileName)
  18.     for line in fr.readlines():
  19.         lineArr =[]
  20.         curLine = line.strip().split('\t')
  21.         for i in range(numFeat):
  22.             lineArr.append(float(curLine[i]))
  23.         xArr.append(lineArr)
  24.         yArr.append(float(curLine[-1]))
  25.     return xArr, yArr
  27. def plotlwlrRegression():
  28.     """
  29.     函数说明:绘制多条局部加权回归曲线
  30.     Parameters:
  31.         无
  32.     Returns:
  33.         无
  35.     """
  36.     font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
  37.     xArr, yArr = loadDataSet('ex0.txt')                                    #加载数据集
  38.     yHat_1 = lwlrTest(xArr, xArr, yArr, 1.0)                            #根据局部加权线性回归计算yHat
  39.     yHat_2 = lwlrTest(xArr, xArr, yArr, 0.01)                            #根据局部加权线性回归计算yHat
  40.     yHat_3 = lwlrTest(xArr, xArr, yArr, 0.003)                            #根据局部加权线性回归计算yHat
  41.     xMat = np.mat(xArr)                                                    #创建xMat矩阵
  42.     yMat = np.mat(yArr)                                                    #创建yMat矩阵
  43.     srtInd = xMat[:, 1].argsort(0)                                        #排序，返回索引值
  44.     xSort = xMat[srtInd][:,0,:]
  45.     fig, axs = plt.subplots(nrows=3, ncols=1,sharex=False, sharey=False, figsize=(10,8))                                        
  46.     axs[0].plot(xSort[:, 1], yHat_1[srtInd], c = 'red')                        #绘制回归曲线
  47.     axs[1].plot(xSort[:, 1], yHat_2[srtInd], c = 'red')                        #绘制回归曲线
  48.     axs[2].plot(xSort[:, 1], yHat_3[srtInd], c = 'red')                        #绘制回归曲线
  49.     axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  50.     axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  51.     axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  52.     #设置标题,x轴label,y轴label
  53.     axs0_title_text = axs[0].set_title(u'局部加权回归曲线,k=1.0',FontProperties=font)
  54.     axs1_title_text = axs[1].set_title(u'局部加权回归曲线,k=0.01',FontProperties=font)
  55.     axs2_title_text = axs[2].set_title(u'局部加权回归曲线,k=0.003',FontProperties=font)
  56.     plt.setp(axs0_title_text, size=8, weight='bold', color='red')  
  57.     plt.setp(axs1_title_text, size=8, weight='bold', color='red')  
  58.     plt.setp(axs2_title_text, size=8, weight='bold', color='red')  
  59.     plt.xlabel('X')
  60.     plt.show()
  61. def lwlr(testPoint, xArr, yArr, k = 1.0):
  62.     """
  63.     函数说明:使用局部加权线性回归计算回归系数w
  64.     Parameters:
  65.         testPoint - 测试样本点
  66.         xArr - x数据集
  67.         yArr - y数据集
  68.         k - 高斯核的k,自定义参数
  69.     Returns:
  70.         ws - 回归系数
  72.     """
  73.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  74.     m = np.shape(xMat)[0]
  75.     weights = np.mat(np.eye((m)))                                        #创建权重对角矩阵
  76.     for j in range(m):                                                  #遍历数据集计算每个样本的权重
  77.         diffMat = testPoint - xMat[j, :]                                 
  78.         weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))
  79.     xTx = xMat.T * (weights * xMat)                                        
  80.     if np.linalg.det(xTx) == 0.0:
  81.         print("矩阵为奇异矩阵,不能求逆")
  82.         return
  83.     ws = xTx.I * (xMat.T * (weights * yMat))                            #计算回归系数
  84.     return testPoint * ws
  85. def lwlrTest(testArr, xArr, yArr, k=1.0):  
  86.     """
  87.     函数说明:局部加权线性回归测试
  88.     Parameters:
  89.         testArr - 测试数据集
  90.         xArr - x数据集
  91.         yArr - y数据集
  92.         k - 高斯核的k,自定义参数
  93.     Returns:
  94.         ws - 回归系数
  96.     """
  97.     m = np.shape(testArr)[0]                                            #计算测试数据集大小
  98.     yHat = np.zeros(m)    
  99.     for i in range(m):                                                    #对每个样本点进行预测
  100.         yHat[i] = lwlr(testArr[i],xArr,yArr,k)
  101.     return yHat
  102. if __name__ == '__main__':
  103.     plotlwlrRegression()

运行结果如下：

13.png

可以看到，当k越小，拟合效果越好。但是当k过小，会出现过拟合的情况，例如k等于0.003的时候。

四、预测鲍鱼的年龄

接下来，我们将回归用于真是数据。在abalone.txt文件中记录了鲍鱼（一种水生物）的年龄，这个数据来自UCI数据集合的数据。鲍鱼年龄可以从鲍鱼壳的层数推算得到。

数据集下载地址：数据集下载

数据集的数据如下：

14.png

可以看到，数据集是多维的，所以我们很难画出它的分布情况。每个维度数据的代表的含义没有给出，不过没有关系，我们只要知道最后一列的数据是y值就可以了，最后一列代表的是鲍鱼的真实年龄，前面几列的数据是一些鲍鱼的特征，例如鲍鱼壳的层数等。我们不做数据清理，直接用上所有特征，测试下我们的局部加权回归。

新建abalone.py文件，添加rssError函数，用于评价最后回归结果。编写代码如下：

  1. # -*- coding:utf-8 -*-
  2. from matplotlib.font_manager import FontProperties
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  6. def loadDataSet(fileName):
  7.     """
  8.     函数说明:加载数据
  9.     Parameters:
  10.         fileName - 文件名
  11.     Returns:
  12.         xArr - x数据集
  13.         yArr - y数据集
  15.     """
  16.     numFeat = len(open(fileName).readline().split('\t')) - 1
  17.     xArr = []; yArr = []
  18.     fr = open(fileName)
  19.     for line in fr.readlines():
  20.         lineArr =[]
  21.         curLine = line.strip().split('\t')
  22.         for i in range(numFeat):
  23.             lineArr.append(float(curLine[i]))
  24.         xArr.append(lineArr)
  25.         yArr.append(float(curLine[-1]))
  26.     return xArr, yArr
  28. def lwlr(testPoint, xArr, yArr, k = 1.0):
  29.     """
  30.     函数说明:使用局部加权线性回归计算回归系数w
  31.     Parameters:
  32.         testPoint - 测试样本点
  33.         xArr - x数据集
  34.         yArr - y数据集
  35.         k - 高斯核的k,自定义参数
  36.     Returns:
  37.         ws - 回归系数
  39.     """
  40.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  41.     m = np.shape(xMat)[0]
  42.     weights = np.mat(np.eye((m)))                                        #创建权重对角矩阵
  43.     for j in range(m):                                                  #遍历数据集计算每个样本的权重
  44.         diffMat = testPoint - xMat[j, :]
  45.         weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))
  46.     xTx = xMat.T * (weights * xMat)
  47.     if np.linalg.det(xTx) == 0.0:
  48.         print("矩阵为奇异矩阵,不能求逆")
  49.         return
  50.     ws = xTx.I * (xMat.T * (weights * yMat))                            #计算回归系数
  51.     return testPoint * ws
  53. def lwlrTest(testArr, xArr, yArr, k=1.0):
  54.     """
  55.     函数说明:局部加权线性回归测试
  56.     Parameters:
  57.         testArr - 测试数据集,测试集
  58.         xArr - x数据集,训练集
  59.         yArr - y数据集,训练集
  60.         k - 高斯核的k,自定义参数
  61.     Returns:
  62.         ws - 回归系数
  64.     """
  65.     m = np.shape(testArr)[0]                                            #计算测试数据集大小
  66.     yHat = np.zeros(m)
  67.     for i in range(m):                                                    #对每个样本点进行预测
  68.         yHat[i] = lwlr(testArr[i],xArr,yArr,k)
  69.     return yHat
  71. def standRegres(xArr,yArr):
  72.     """
  73.     函数说明:计算回归系数w
  74.     Parameters:
  75.         xArr - x数据集
  76.         yArr - y数据集
  77.     Returns:
  78.         ws - 回归系数
  80.     """
  81.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  82.     xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
  83.     if np.linalg.det(xTx) == 0.0:
  84.         print("矩阵为奇异矩阵,不能求逆")
  85.         return
  86.     ws = xTx.I * (xMat.T*yMat)
  87.     return ws
  89. def rssError(yArr, yHatArr):
  90.     """
  91.     误差大小评价函数
  92.     Parameters:
  93.         yArr - 真实数据
  94.         yHatArr - 预测数据
  95.     Returns:
  96.         误差大小
  97.     """
  98.     return ((yArr - yHatArr) **2).sum()
  100. if __name__ == '__main__':
  101.     abX, abY = loadDataSet('abalone.txt')
  102.     print('训练集与测试集相同:局部加权线性回归,核k的大小对预测的影响:')
  103.     yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)
  104.     yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
  105.     yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
  106.     print('k=0.1时,误差大小为:',rssError(abY[0:99], yHat01.T))
  107.     print('k=1  时,误差大小为:',rssError(abY[0:99], yHat1.T))
  108.     print('k=10 时,误差大小为:',rssError(abY[0:99], yHat10.T))
  110.     print('')
  112.     print('训练集与测试集不同:局部加权线性回归,核k的大小是越小越好吗？更换数据集,测试结果如下:')
  113.     yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
  114.     yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
  115.     yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
  116.     print('k=0.1时,误差大小为:',rssError(abY[100:199], yHat01.T))
  117.     print('k=1  时,误差大小为:',rssError(abY[100:199], yHat1.T))
  118.     print('k=10 时,误差大小为:',rssError(abY[100:199], yHat10.T))
  120.     print('')
  122.     print('训练集与测试集不同:简单的线性归回与k=1时的局部加权线性回归对比:')
  123.     print('k=1时,误差大小为:', rssError(abY[100:199], yHat1.T))
  124.     ws = standRegres(abX[0:99], abY[0:99])
  125.     yHat = np.mat(abX[100:199]) * ws
  126.     print('简单的线性回归误差大小:', rssError(abY[100:199], yHat.T.A))

运行结果如下：

15.png

可以看到，当k=0.1时，训练集误差小，但是应用于新的数据集之后，误差反而变大了。这就是经常说道的过拟合现象。我们训练的模型，我们要保证测试集准确率高，这样训练出的模型才可以应用于新的数据，也就是要加强模型的普适性。可以看到，当k=1时，局部加权线性回归和简单的线性回归得到的效果差不多。这也表明一点，必须在未知数据上比较效果才能选取到最佳模型。那么最佳的核大小是10吗？或许是，但如果想得到更好的效果，应该用10个不同的样本集做10次测试来比较结果。

本示例展示了如何使用局部加权线性回归来构建模型，可以得到比普通线性回归更好的效果。局部加权线性回归的问题在于，每次必须在整个数据集上运行。也就是说为了做出预测，必须保存所有的训练数据。

五、岭回归

如果数据的特征比样本点还多应该怎么办？很显然，此时我们不能再使用上文的方法进行计算了，因为矩阵X不是满秩矩阵，非满秩矩阵在求逆时会出现问题。为了解决这个问题，统计学家引入岭回归（ridge regression）的概念。

1、什么是岭回归？

岭回归即我们所说的L2正则线性回归，在一般的线性回归最小化均方误差的基础上增加了一个参数w的L2范数的罚项，从而最小化罚项残差平方和：

16.png

简单说来，岭回归就是在普通线性回归的基础上引入单位矩阵。回归系数的计算公式变形如下：

17.png

式中，矩阵I是一个mxm的单位矩阵，加上一个λI从而使得矩阵非奇异，进而能对矩阵求逆。

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入λ来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也可以叫做缩减（shrinkage）。

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

2、编写代码

为了使用岭回归和缩减技术，首先需要对特征做标准化处理。因为，我们需要使每个维度特征具有相同的重要性。本文使用的标准化处理比较简单，就是将所有特征都减去各自的均值并除以方差。

代码很简单，只需要稍做修改，其中，λ为模型的参数。我们先绘制一个回归系数与log(λ)的曲线图，看下它们的规律，编写代码如下：

  1. # -*-coding:utf-8 -*-
  2. from matplotlib.font_manager import FontProperties
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  5. def loadDataSet(fileName):
  6.     """
  7.     函数说明:加载数据
  8.     Parameters:
  9.         fileName - 文件名
  10.     Returns:
  11.         xArr - x数据集
  12.         yArr - y数据集
  14.     """
  15.     numFeat = len(open(fileName).readline().split('\t')) - 1
  16.     xArr = []; yArr = []
  17.     fr = open(fileName)
  18.     for line in fr.readlines():
  19.         lineArr =[]
  20.         curLine = line.strip().split('\t')
  21.         for i in range(numFeat):
  22.             lineArr.append(float(curLine[i]))
  23.         xArr.append(lineArr)
  24.         yArr.append(float(curLine[-1]))
  25.     return xArr, yArr
  26. def ridgeRegres(xMat, yMat, lam = 0.2):
  27.     """
  28.     函数说明:岭回归
  29.     Parameters:
  30.         xMat - x数据集
  31.         yMat - y数据集
  32.         lam - 缩减系数
  33.     Returns:
  34.         ws - 回归系数
  36.     """
  37.     xTx = xMat.T * xMat
  38.     denom = xTx + np.eye(np.shape(xMat)[1]) * lam
  39.     if np.linalg.det(denom) == 0.0:
  40.         print("矩阵为奇异矩阵,不能转置")
  41.         return
  42.     ws = denom.I * (xMat.T * yMat)
  43.     return ws
  44. def ridgeTest(xArr, yArr):
  45.     """
  46.     函数说明:岭回归测试
  47.     Parameters:
  48.         xMat - x数据集
  49.         yMat - y数据集
  50.     Returns:
  51.         wMat - 回归系数矩阵
  53.     """
  54.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  55.     #数据标准化
  56.     yMean = np.mean(yMat, axis = 0)                        #行与行操作，求均值
  57.     yMat = yMat - yMean                                    #数据减去均值
  58.     xMeans = np.mean(xMat, axis = 0)                    #行与行操作，求均值
  59.     xVar = np.var(xMat, axis = 0)                        #行与行操作，求方差
  60.     xMat = (xMat - xMeans) / xVar                        #数据减去均值除以方差实现标准化
  61.     numTestPts = 30                                        #30个不同的lambda测试
  62.     wMat = np.zeros((numTestPts, np.shape(xMat)[1]))    #初始回归系数矩阵
  63.     for i in range(numTestPts):                            #改变lambda计算回归系数
  64.         ws = ridgeRegres(xMat, yMat, np.exp(i - 10))    #lambda以e的指数变化，最初是一个非常小的数，
  65.         wMat[i, :] = ws.T                                 #计算回归系数矩阵
  66.     return wMat
  67. def plotwMat():
  68.     """
  69.     函数说明:绘制岭回归系数矩阵
  71.     """
  72.     font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
  73.     abX, abY = loadDataSet('abalone.txt')
  74.     redgeWeights = ridgeTest(abX, abY)
  75.     fig = plt.figure()
  76.     ax = fig.add_subplot(111)
  77.     ax.plot(redgeWeights)    
  78.     ax_title_text = ax.set_title(u'log(lambada)与回归系数的关系', FontProperties = font)
  79.     ax_xlabel_text = ax.set_xlabel(u'log(lambada)', FontProperties = font)
  80.     ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties = font)
  81.     plt.setp(ax_title_text, size = 20, weight = 'bold', color = 'red')
  82.     plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black')
  83.     plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black')
  84.     plt.show()
  85. if __name__ == '__main__':
  86.     plotwMat()

运行结果如下：

18.png

图绘制了回归系数与log(λ)的关系。在最左边，即λ最小时，可以得到所有系数的原始值（与线性回归一致）；而在右边，系数全部缩减成0；在中间部分的某个位置，将会得到最好的预测结果。想要得到最佳的λ参数，可以使用交叉验证的方式获得，文章的后面会继续讲解。

六、前向逐步线性回归

前向逐步线性回归算法属于一种贪心算法，即每一步都尽可能减少误差。我们计算回归系数，不再是通过公式计算，而是通过每次微调各个回归系数，然后计算预测误差。那个使误差最小的一组回归系数，就是我们需要的最佳回归系数。

前向逐步线性回归实现也很简单。当然，还是先进行数据标准化，编写代码如下：

  1. # -*-coding:utf-8 -*-
  2. from matplotlib.font_manager import FontProperties
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  5. def loadDataSet(fileName):
  6.     """
  7.     函数说明:加载数据
  8.     Parameters:
  9.         fileName - 文件名
  10.     Returns:
  11.         xArr - x数据集
  12.         yArr - y数据集
  13.     Website:
  14.         http://www.cuijiahua.com/
  15.     Modify:
  16.         2017-11-20
  17.     """
  18.     numFeat = len(open(fileName).readline().split('\t')) - 1
  19.     xArr = []; yArr = []
  20.     fr = open(fileName)
  21.     for line in fr.readlines():
  22.         lineArr =[]
  23.         curLine = line.strip().split('\t')
  24.         for i in range(numFeat):
  25.             lineArr.append(float(curLine[i]))
  26.         xArr.append(lineArr)
  27.         yArr.append(float(curLine[-1]))
  28.     return xArr, yArr
  30. def regularize(xMat, yMat):
  31.     """
  32.     函数说明:数据标准化
  33.     Parameters:
  34.         xMat - x数据集
  35.         yMat - y数据集
  36.     Returns:
  37.         inxMat - 标准化后的x数据集
  38.         inyMat - 标准化后的y数据集
  39.     Website:
  40.         http://www.cuijiahua.com/
  41.     Modify:
  42.         2017-11-23
  43.     """    
  44.     inxMat = xMat.copy()                                                        #数据拷贝
  45.     inyMat = yMat.copy()
  46.     yMean = np.mean(yMat, 0)                                                    #行与行操作，求均值
  47.     inyMat = yMat - yMean                                                        #数据减去均值
  48.     inMeans = np.mean(inxMat, 0)                                                   #行与行操作，求均值
  49.     inVar = np.var(inxMat, 0)                                                     #行与行操作，求方差
  50.     inxMat = (inxMat - inMeans) / inVar                                            #数据减去均值除以方差实现标准化
  51.     return inxMat, inyMat
  53. def rssError(yArr,yHatArr):
  54.     """
  55.     函数说明:计算平方误差
  56.     Parameters:
  57.         yArr - 预测值
  58.         yHatArr - 真实值
  59.     Returns:
  61.     Website:
  62.         http://www.cuijiahua.com/
  63.     Modify:
  64.         2017-11-23
  65.     """
  66.     return ((yArr-yHatArr)**2).sum()
  67. def stageWise(xArr, yArr, eps = 0.01, numIt = 100):
  68.     """
  69.     函数说明:前向逐步线性回归
  70.     Parameters:
  71.         xArr - x输入数据
  72.         yArr - y预测数据
  73.         eps - 每次迭代需要调整的步长
  74.         numIt - 迭代次数
  75.     Returns:
  76.         returnMat - numIt次迭代的回归系数矩阵
  77.     Website:
  78.         http://www.cuijiahua.com/
  79.     Modify:
  80.         2017-12-03
  81.     """
  82.     xMat = np.mat(xArr); yMat = np.mat(yArr).T                                         #数据集
  83.     xMat, yMat = regularize(xMat, yMat)                                                #数据标准化
  84.     m, n = np.shape(xMat)
  85.     returnMat = np.zeros((numIt, n))                                                #初始化numIt次迭代的回归系数矩阵
  86.     ws = np.zeros((n, 1))                                                            #初始化回归系数矩阵
  87.     wsTest = ws.copy()
  88.     wsMax = ws.copy()
  89.     for i in range(numIt):                                                            #迭代numIt次
  90.         # print(ws.T)                                                                    #打印当前回归系数矩阵
  91.         lowestError = float('inf');                                                 #正无穷
  92.         for j in range(n):                                                            #遍历每个特征的回归系数
  93.             for sign in [-1, 1]:
  94.                 wsTest = ws.copy()
  95.                 wsTest[j] += eps * sign                                                #微调回归系数
  96.                 yTest = xMat * wsTest                                                #计算预测值
  97.                 rssE = rssError(yMat.A, yTest.A)                                    #计算平方误差
  98.                 if rssE < lowestError:                                                #如果误差更小，则更新当前的最佳回归系数
  99.                     lowestError = rssE
  100.                     wsMax = wsTest
  101.         ws = wsMax.copy()
  102.         returnMat[i,:] = ws.T                                                         #记录numIt次迭代的回归系数矩阵
  103.     return returnMat
  104. def plotstageWiseMat():
  105.     """
  106.     函数说明:绘制岭回归系数矩阵
  107.     Website:
  108.         http://www.cuijiahua.com/
  109.     Modify:
  110.         2017-11-20
  111.     """
  112.     font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
  113.     xArr, yArr = loadDataSet('abalone.txt')
  114.     returnMat = stageWise(xArr, yArr, 0.005, 1000)
  115.     fig = plt.figure()
  116.     ax = fig.add_subplot(111)
  117.     ax.plot(returnMat)    
  118.     ax_title_text = ax.set_title(u'前向逐步回归:迭代次数与回归系数的关系', FontProperties = font)
  119.     ax_xlabel_text = ax.set_xlabel(u'迭代次数', FontProperties = font)
  120.     ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties = font)
  121.     plt.setp(ax_title_text, size = 15, weight = 'bold', color = 'red')
  122.     plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black')
  123.     plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black')
  124.     plt.show()
  125. if __name__ == '__main__':
  126.     plotstageWiseMat()

运行结果如下：

19.png

还是，我们打印了迭代次数与回归系数的关系曲线。可以看到，有些系数从始至终都是约为0的，这说明它们不对目标造成任何影响，也就是说这些特征很可能是不需要的。逐步线性回归算法的优点在于它可以帮助人们理解有的模型并做出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。

总结一下：

缩减方法（逐步线性回归或岭回归），就是将一些系数缩减成很小的值或者直接缩减为0。这样做，就增大了模型的偏差（减少了一些特征的权重），通过把一些特征的回归系数缩减到0，同时也就减少了模型的复杂度。消除了多余的特征之后，模型更容易理解，同时也降低了预测误差。但是当缩减过于严厉的时候，就会出现过拟合的现象，即用训练集预测结果很好，用测试集预测就糟糕很多。

七、预测乐高玩具套件的价格

乐高（LEGO）公司生产拼装类玩具，由很多大小不同的塑料插块组成。它的样子如下图所示：

20.png

一般来说，这些插块都是成套出售，它们可以拼装成很多不同的东西，如船、城堡、一些著名建筑等。乐高公司每个套装包含的部件数目从10件到5000件不等。

一种乐高套件基本上在几年后就会停产，但乐高的收藏者之间仍会在停产后彼此交易。本次实例，就是使用回归方法对收藏者之间的交易价格进行预测。

1、获取数据

书中使用的方法是通过Google提供的API进行获取数据，但是现在这个API已经关闭，我们无法通过api获取数据了。不过幸运的是，在网上找到了书上用到的那些html文件。

原始数据下载地址：数据下载

我们通过解析html文件，来获取我们需要的信息,解析代码如下：

  1. # -*-coding:utf-8 -*-
  2. from bs4 import BeautifulSoup
  4. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  5.     """
  6.     函数说明:从页面读取数据，生成retX和retY列表
  7.     Parameters:
  8.         retX - 数据X
  9.         retY - 数据Y
  10.         inFile - HTML文件
  11.         yr - 年份
  12.         numPce - 乐高部件数目
  13.         origPrc - 原价
  14.     Returns:
  15.         无
  17.     """
  18.     # 打开并读取HTML文件
  19.     with open(inFile, encoding='utf-8') as f:
  20.         html = f.read()
  21.     soup = BeautifulSoup(html)
  22.     i = 1
  23.     # 根据HTML页面结构进行解析
  24.     currentRow = soup.find_all('table', r = "%d" % i)
  25.     while(len(currentRow) != 0):
  26.         currentRow = soup.find_all('table', r = "%d" % i)
  27.         title = currentRow[0].find_all('a')[1].text
  28.         lwrTitle = title.lower()
  29.         # 查找是否有全新标签
  30.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  31.             newFlag = 1.0
  32.         else:
  33.             newFlag = 0.0
  34.         # 查找是否已经标志出售，我们只收集已出售的数据
  35.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
  36.         if len(soldUnicde) == 0:
  37.             print("商品 #%d 没有出售" % i)
  38.         else:
  39.             # 解析页面获取当前价格
  40.             soldPrice = currentRow[0].find_all('td')[4]
  41.             priceStr = soldPrice.text
  42.             priceStr = priceStr.replace('$','')
  43.             priceStr = priceStr.replace(',','')
  44.             if len(soldPrice) > 1:
  45.                 priceStr = priceStr.replace('Free shipping', '')
  46.             sellingPrice = float(priceStr)
  47.             # 去掉不完整的套装价格
  48.             if  sellingPrice > origPrc * 0.5:
  49.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
  50.                 retX.append([yr, numPce, newFlag, origPrc])
  51.                 retY.append(sellingPrice)
  52.         i += 1
  53.         currentRow = soup.find_all('table', r = "%d" % i)
  55. def setDataCollect(retX, retY):
  56.     """
  57.     函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
  58.     Parameters:
  59.         无
  60.     Returns:
  61.         无
  63.     """
  64.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)                #2006年的乐高8288,部件数目800,原价49.99
  65.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)                #2002年的乐高10030,部件数目3096,原价269.99
  66.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)                #2007年的乐高10179,部件数目5195,原价499.99
  67.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)                #2007年的乐高10181,部件数目3428,原价199.99
  68.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)                #2008年的乐高10189,部件数目5922,原价299.99
  69.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)                #2009年的乐高10196,部件数目3263,原价249.99
  71. if __name__ == '__main__':
  72.     lgX = []
  73.     lgY = []
  74.     setDataCollect(lgX, lgY)

运行结果如下：

21.png

我们对没有的商品做了处理。这些特征分别为：出品年份、部件数目、是否为全新、原价、售价（二手交易）。

html解析页面不会使用？那就学习一下爬虫知识吧~！如果对此不感兴趣，也可以跳过获取数据和解析数据，这个过程，当作已知数据，继续进行下一步。

2、建立模型

我们已经处理好了数据集，接下来就是训练模型。首先我们需要添加全为0的特征X0列。因为线性回归的第一列特征要求都是1.0。然后使用最简单的普通线性回归i，编写代码如下：

  1. # -*-coding:utf-8 -*-
  2. import numpy as np
  3. from bs4 import BeautifulSoup
  4. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  5.     """
  6.     函数说明:从页面读取数据，生成retX和retY列表
  7.     Parameters:
  8.         retX - 数据X
  9.         retY - 数据Y
  10.         inFile - HTML文件
  11.         yr - 年份
  12.         numPce - 乐高部件数目
  13.         origPrc - 原价
  14.     Returns:
  15.         无
  17.     """
  18.     # 打开并读取HTML文件
  19.     with open(inFile, encoding='utf-8') as f:
  20.         html = f.read()
  21.     soup = BeautifulSoup(html)
  22.     i = 1
  23.     # 根据HTML页面结构进行解析
  24.     currentRow = soup.find_all('table', r = "%d" % i)
  25.     while(len(currentRow) != 0):
  26.         currentRow = soup.find_all('table', r = "%d" % i)
  27.         title = currentRow[0].find_all('a')[1].text
  28.         lwrTitle = title.lower()
  29.         # 查找是否有全新标签
  30.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  31.             newFlag = 1.0
  32.         else:
  33.             newFlag = 0.0
  34.         # 查找是否已经标志出售，我们只收集已出售的数据
  35.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
  36.         if len(soldUnicde) == 0:
  37.             print("商品 #%d 没有出售" % i)
  38.         else:
  39.             # 解析页面获取当前价格
  40.             soldPrice = currentRow[0].find_all('td')[4]
  41.             priceStr = soldPrice.text
  42.             priceStr = priceStr.replace('$','')
  43.             priceStr = priceStr.replace(',','')
  44.             if len(soldPrice) > 1:
  45.                 priceStr = priceStr.replace('Free shipping', '')
  46.             sellingPrice = float(priceStr)
  47.             # 去掉不完整的套装价格
  48.             if  sellingPrice > origPrc * 0.5:
  49.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
  50.                 retX.append([yr, numPce, newFlag, origPrc])
  51.                 retY.append(sellingPrice)
  52.         i += 1
  53.         currentRow = soup.find_all('table', r = "%d" % i)
  54. #         
  55. def setDataCollect(retX, retY):
  56.     """
  57.     函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
  58.     Parameters:
  59.         无
  60.     Returns:
  61.         无
  63.     """
  64.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)                #2006年的乐高8288,部件数目800,原价49.99
  65.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)                #2002年的乐高10030,部件数目3096,原价269.99
  66.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)                #2007年的乐高10179,部件数目5195,原价499.99
  67.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)                #2007年的乐高10181,部件数目3428,原价199.99
  68.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)                #2008年的乐高10189,部件数目5922,原价299.99
  69.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)                #2009年的乐高10196,部件数目3263,原价249.99
  70. def regularize(xMat, yMat):
  71.     """
  72.     函数说明:数据标准化
  73.     Parameters:
  74.         xMat - x数据集
  75.         yMat - y数据集
  76.     Returns:
  77.         inxMat - 标准化后的x数据集
  78.         inyMat - 标准化后的y数据集
  80.     """    
  81.     inxMat = xMat.copy()                                                        #数据拷贝
  82.     inyMat = yMat.copy()
  83.     yMean = np.mean(yMat, 0)                                                    #行与行操作，求均值
  84.     inyMat = yMat - yMean                                                        #数据减去均值
  85.     inMeans = np.mean(inxMat, 0)                                                   #行与行操作，求均值
  86.     inVar = np.var(inxMat, 0)                                                     #行与行操作，求方差
  87.     # print(inxMat)
  88.     print(inMeans)
  89.     # print(inVar)
  90.     inxMat = (inxMat - inMeans) / inVar                                            #数据减去均值除以方差实现标准化
  91.     return inxMat, inyMat
  92. def rssError(yArr,yHatArr):
  93.     """
  94.     函数说明:计算平方误差
  95.     Parameters:
  96.         yArr - 预测值
  97.         yHatArr - 真实值
  98.     Returns:
  101.     """
  102.     return ((yArr-yHatArr)**2).sum()
  103. def standRegres(xArr,yArr):
  104.     """
  105.     函数说明:计算回归系数w
  106.     Parameters:
  107.         xArr - x数据集
  108.         yArr - y数据集
  109.     Returns:
  110.         ws - 回归系数
  112.     """
  113.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  114.     xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
  115.     if np.linalg.det(xTx) == 0.0:
  116.         print("矩阵为奇异矩阵,不能转置")
  117.         return
  118.     ws = xTx.I * (xMat.T*yMat)
  119.     return ws
  121. def useStandRegres():
  122.     """
  123.     函数说明:使用简单的线性回归
  124.     Parameters:
  125.         无
  126.     Returns:
  127.         无
  129.     """
  130.     lgX = []
  131.     lgY = []
  132.     setDataCollect(lgX, lgY)
  133.     data_num, features_num = np.shape(lgX)
  134.     lgX1 = np.mat(np.ones((data_num, features_num + 1)))
  135.     lgX1[:, 1:5] = np.mat(lgX)
  136.     ws = standRegres(lgX1, lgY)
  137.     print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0],ws[1],ws[2],ws[3],ws[4]))     
  139. if __name__ == '__main__':
  140.     useStandRegres()

运行结果如下图所示：

22.png
可以看到，模型采用的公式如上图所示。虽然这个模型对于数据拟合得很好，但是看上不没有什么道理。套件里的部件数量越多，售价反而降低了，这是不合理的。

我们使用岭回归，通过交叉验证，找到使误差最小的λ对应的回归系数。编写代码如下：

  1. # -*-coding:utf-8 -*-
  2. import numpy as np
  3. from bs4 import BeautifulSoup
  4. import random
  6. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  7.     """
  8.     函数说明:从页面读取数据，生成retX和retY列表
  9.     Parameters:
  10.         retX - 数据X
  11.         retY - 数据Y
  12.         inFile - HTML文件
  13.         yr - 年份
  14.         numPce - 乐高部件数目
  15.         origPrc - 原价
  16.     Returns:
  17.         无
  19.     """
  20.     # 打开并读取HTML文件
  21.     with open(inFile, encoding='utf-8') as f:
  22.         html = f.read()
  23.     soup = BeautifulSoup(html)
  24.     i = 1
  25.     # 根据HTML页面结构进行解析
  26.     currentRow = soup.find_all('table', r = "%d" % i)
  27.     while(len(currentRow) != 0):
  28.         currentRow = soup.find_all('table', r = "%d" % i)
  29.         title = currentRow[0].find_all('a')[1].text
  30.         lwrTitle = title.lower()
  31.         # 查找是否有全新标签
  32.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  33.             newFlag = 1.0
  34.         else:
  35.             newFlag = 0.0
  36.         # 查找是否已经标志出售，我们只收集已出售的数据
  37.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
  38.         if len(soldUnicde) == 0:
  39.             print("商品 #%d 没有出售" % i)
  40.         else:
  41.             # 解析页面获取当前价格
  42.             soldPrice = currentRow[0].find_all('td')[4]
  43.             priceStr = soldPrice.text
  44.             priceStr = priceStr.replace('$','')
  45.             priceStr = priceStr.replace(',','')
  46.             if len(soldPrice) > 1:
  47.                 priceStr = priceStr.replace('Free shipping', '')
  48.             sellingPrice = float(priceStr)
  49.             # 去掉不完整的套装价格
  50.             if  sellingPrice > origPrc * 0.5:
  51.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
  52.                 retX.append([yr, numPce, newFlag, origPrc])
  53.                 retY.append(sellingPrice)
  54.         i += 1
  55.         currentRow = soup.find_all('table', r = "%d" % i)
  57. def ridgeRegres(xMat, yMat, lam = 0.2):
  58.     """
  59.     函数说明:岭回归
  60.     Parameters:
  61.         xMat - x数据集
  62.         yMat - y数据集
  63.         lam - 缩减系数
  64.     Returns:
  65.         ws - 回归系数
  67.     """
  68.     xTx = xMat.T * xMat
  69.     denom = xTx + np.eye(np.shape(xMat)[1]) * lam
  70.     if np.linalg.det(denom) == 0.0:
  71.         print("矩阵为奇异矩阵,不能转置")
  72.         return
  73.     ws = denom.I * (xMat.T * yMat)
  74.     return ws
  76. def setDataCollect(retX, retY):
  77.     """
  78.     函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
  79.     Parameters:
  80.         无
  81.     Returns:
  82.         无
  84.     """
  85.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)                #2006年的乐高8288,部件数目800,原价49.99
  86.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)                #2002年的乐高10030,部件数目3096,原价269.99
  87.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)                #2007年的乐高10179,部件数目5195,原价499.99
  88.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)                #2007年的乐高10181,部件数目3428,原价199.99
  89.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)                #2008年的乐高10189,部件数目5922,原价299.99
  90.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)                #2009年的乐高10196,部件数目3263,原价249.99
  92. def regularize(xMat, yMat):
  93.     """
  94.     函数说明:数据标准化
  95.     Parameters:
  96.         xMat - x数据集
  97.         yMat - y数据集
  98.     Returns:
  99.         inxMat - 标准化后的x数据集
  100.         inyMat - 标准化后的y数据集
  102.     """    
  103.     inxMat = xMat.copy()                                                        #数据拷贝
  104.     inyMat = yMat.copy()
  105.     yMean = np.mean(yMat, 0)                                                    #行与行操作，求均值
  106.     inyMat = yMat - yMean                                                        #数据减去均值
  107.     inMeans = np.mean(inxMat, 0)                                                   #行与行操作，求均值
  108.     inVar = np.var(inxMat, 0)                                                     #行与行操作，求方差
  109.     # print(inxMat)
  110.     print(inMeans)
  111.     # print(inVar)
  112.     inxMat = (inxMat - inMeans) / inVar                                            #数据减去均值除以方差实现标准化
  113.     return inxMat, inyMat
  115. def rssError(yArr,yHatArr):
  116.     """
  117.     函数说明:计算平方误差
  118.     Parameters:
  119.         yArr - 预测值
  120.         yHatArr - 真实值
  121.     Returns:
  124.     """
  125.     return ((yArr-yHatArr)**2).sum()
  126. def standRegres(xArr,yArr):
  127.     """
  128.     函数说明:计算回归系数w
  129.     Parameters:
  130.         xArr - x数据集
  131.         yArr - y数据集
  132.     Returns:
  133.         ws - 回归系数
  135.     """
  136.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  137.     xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
  138.     if np.linalg.det(xTx) == 0.0:
  139.         print("矩阵为奇异矩阵,不能转置")
  140.         return
  141.     ws = xTx.I * (xMat.T*yMat)
  142.     return ws
  144. def crossValidation(xArr, yArr, numVal = 10):
  145.     """
  146.     函数说明:交叉验证岭回归
  147.     Parameters:
  148.         xArr - x数据集
  149.         yArr - y数据集
  150.         numVal - 交叉验证次数
  151.     Returns:
  152.         wMat - 回归系数矩阵
  154.     """
  155.     m = len(yArr)                                                                        #统计样本个数                       
  156.     indexList = list(range(m))                                                            #生成索引值列表
  157.     errorMat = np.zeros((numVal,30))                                                    #create error mat 30columns numVal rows
  158.     for i in range(numVal):                                                                #交叉验证numVal次
  159.         trainX = []; trainY = []                                                        #训练集
  160.         testX = []; testY = []                                                            #测试集
  161.         random.shuffle(indexList)                                                        #打乱次序
  162.         for j in range(m):                                                                #划分数据集:90%训练集，10%测试集
  163.             if j < m * 0.9:
  164.                 trainX.append(xArr[indexList[j]])
  165.                 trainY.append(yArr[indexList[j]])
  166.             else:
  167.                 testX.append(xArr[indexList[j]])
  168.                 testY.append(yArr[indexList[j]])
  169.         wMat = ridgeTest(trainX, trainY)                                                #获得30个不同lambda下的岭回归系数
  170.         for k in range(30):                                                                #遍历所有的岭回归系数
  171.             matTestX = np.mat(testX); matTrainX = np.mat(trainX)                        #测试集
  172.             meanTrain = np.mean(matTrainX,0)                                            #测试集均值
  173.             varTrain = np.var(matTrainX,0)                                                #测试集方差
  174.             matTestX = (matTestX - meanTrain) / varTrain                                 #测试集标准化
  175.             yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY)                        #根据ws预测y值
  176.             errorMat[i, k] = rssError(yEst.T.A, np.array(testY))                            #统计误差
  177.     meanErrors = np.mean(errorMat,0)                                                    #计算每次交叉验证的平均误差
  178.     minMean = float(min(meanErrors))                                                    #找到最小误差
  179.     bestWeights = wMat[np.nonzero(meanErrors == minMean)]                                #找到最佳回归系数
  180.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  181.     meanX = np.mean(xMat,0); varX = np.var(xMat,0)
  182.     unReg = bestWeights / varX                                                            #数据经过标准化，因此需要还原
  183.     print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % ((-1 * np.sum(np.multiply(meanX,unReg)) + np.mean(yMat)), unReg[0,0], unReg[0,1], unReg[0,2], unReg[0,3]))    
  185. def ridgeTest(xArr, yArr):
  186.     """
  187.     函数说明:岭回归测试
  188.     Parameters:
  189.         xMat - x数据集
  190.         yMat - y数据集
  191.     Returns:
  192.         wMat - 回归系数矩阵
  194.     """
  195.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  196.     #数据标准化
  197.     yMean = np.mean(yMat, axis = 0)                        #行与行操作，求均值
  198.     yMat = yMat - yMean                                    #数据减去均值
  199.     xMeans = np.mean(xMat, axis = 0)                    #行与行操作，求均值
  200.     xVar = np.var(xMat, axis = 0)                        #行与行操作，求方差
  201.     xMat = (xMat - xMeans) / xVar                        #数据减去均值除以方差实现标准化
  202.     numTestPts = 30                                        #30个不同的lambda测试
  203.     wMat = np.zeros((numTestPts, np.shape(xMat)[1]))    #初始回归系数矩阵
  204.     for i in range(numTestPts):                            #改变lambda计算回归系数
  205.         ws = ridgeRegres(xMat, yMat, np.exp(i - 10))    #lambda以e的指数变化，最初是一个非常小的数，
  206.         wMat[i, :] = ws.T                                 #计算回归系数矩阵
  207.     return wMat
  209. if __name__ == '__main__':
  210.     lgX = []
  211.     lgY = []
  212.     setDataCollect(lgX, lgY)
  213.     crossValidation(lgX, lgY)

运行结果如下图所示：

23.png

这里随机选取样本，因为其随机性，所以每次运行的结果可能略有不同。不过整体如上图所示，可以看出，它与常规的最小二乘法，即普通的线性回归没有太大差异。我们本期望找到一个更易于理解的模型，显然没有达到预期效果。

现在，我们看一下在缩减过程中回归系数是如何变化的。编写代码如下：

  1. # -*-coding:utf-8 -*-
  2. import numpy as np
  3. from bs4 import BeautifulSoup
  4. import random
  6. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  7.     """
  8.     函数说明:从页面读取数据，生成retX和retY列表
  9.     Parameters:
  10.         retX - 数据X
  11.         retY - 数据Y
  12.         inFile - HTML文件
  13.         yr - 年份
  14.         numPce - 乐高部件数目
  15.         origPrc - 原价
  16.     Returns:
  17.         无
  19.     """
  20.     # 打开并读取HTML文件
  21.     with open(inFile, encoding='utf-8') as f:
  22.         html = f.read()
  23.     soup = BeautifulSoup(html)
  24.     i = 1
  25.     # 根据HTML页面结构进行解析
  26.     currentRow = soup.find_all('table', r = "%d" % i)
  27.     while(len(currentRow) != 0):
  28.         currentRow = soup.find_all('table', r = "%d" % i)
  29.         title = currentRow[0].find_all('a')[1].text
  30.         lwrTitle = title.lower()
  31.         # 查找是否有全新标签
  32.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  33.             newFlag = 1.0
  34.         else:
  35.             newFlag = 0.0
  36.         # 查找是否已经标志出售，我们只收集已出售的数据
  37.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
  38.         if len(soldUnicde) == 0:
  39.             print("商品 #%d 没有出售" % i)
  40.         else:
  41.             # 解析页面获取当前价格
  42.             soldPrice = currentRow[0].find_all('td')[4]
  43.             priceStr = soldPrice.text
  44.             priceStr = priceStr.replace('$','')
  45.             priceStr = priceStr.replace(',','')
  46.             if len(soldPrice) > 1:
  47.                 priceStr = priceStr.replace('Free shipping', '')
  48.             sellingPrice = float(priceStr)
  49.             # 去掉不完整的套装价格
  50.             if  sellingPrice > origPrc * 0.5:
  51.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
  52.                 retX.append([yr, numPce, newFlag, origPrc])
  53.                 retY.append(sellingPrice)
  54.         i += 1
  55.         currentRow = soup.find_all('table', r = "%d" % i)
  57. def ridgeRegres(xMat, yMat, lam = 0.2):
  58.     """
  59.     函数说明:岭回归
  60.     Parameters:
  61.         xMat - x数据集
  62.         yMat - y数据集
  63.         lam - 缩减系数
  64.     Returns:
  65.         ws - 回归系数
  67.     """
  68.     xTx = xMat.T * xMat
  69.     denom = xTx + np.eye(np.shape(xMat)[1]) * lam
  70.     if np.linalg.det(denom) == 0.0:
  71.         print("矩阵为奇异矩阵,不能转置")
  72.         return
  73.     ws = denom.I * (xMat.T * yMat)
  74.     return ws
  75. def setDataCollect(retX, retY):
  76.     """
  77.     函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
  78.     Parameters:
  79.         无
  80.     Returns:
  81.         无
  83.     """
  84.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)                #2006年的乐高8288,部件数目800,原价49.99
  85.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)                #2002年的乐高10030,部件数目3096,原价269.99
  86.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)                #2007年的乐高10179,部件数目5195,原价499.99
  87.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)                #2007年的乐高10181,部件数目3428,原价199.99
  88.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)                #2008年的乐高10189,部件数目5922,原价299.99
  89.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)                #2009年的乐高10196,部件数目3263,原价249.99
  91. def regularize(xMat, yMat):
  92.     """
  93.     函数说明:数据标准化
  94.     Parameters:
  95.         xMat - x数据集
  96.         yMat - y数据集
  97.     Returns:
  98.         inxMat - 标准化后的x数据集
  99.         inyMat - 标准化后的y数据集
  101.     """    
  102.     inxMat = xMat.copy()                                                        #数据拷贝
  103.     inyMat = yMat.copy()
  104.     yMean = np.mean(yMat, 0)                                                    #行与行操作，求均值
  105.     inyMat = yMat - yMean                                                        #数据减去均值
  106.     inMeans = np.mean(inxMat, 0)                                                   #行与行操作，求均值
  107.     inVar = np.var(inxMat, 0)                                                     #行与行操作，求方差
  108.     # print(inxMat)
  109.     print(inMeans)
  110.     # print(inVar)
  111.     inxMat = (inxMat - inMeans) / inVar                                            #数据减去均值除以方差实现标准化
  112.     return inxMat, inyMat
  114. def rssError(yArr,yHatArr):
  115.     """
  116.     函数说明:计算平方误差
  117.     Parameters:
  118.         yArr - 预测值
  119.         yHatArr - 真实值
  120.     Returns:
  122.     """
  123.     return ((yArr-yHatArr)**2).sum()
  124. def standRegres(xArr,yArr):
  125.     """
  126.     函数说明:计算回归系数w
  127.     Parameters:
  128.         xArr - x数据集
  129.         yArr - y数据集
  130.     Returns:
  131.         ws - 回归系数
  133.     """
  134.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  135.     xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
  136.     if np.linalg.det(xTx) == 0.0:
  137.         print("矩阵为奇异矩阵,不能转置")
  138.         return
  139.     ws = xTx.I * (xMat.T*yMat)
  140.     return ws
  142. def ridgeTest(xArr, yArr):
  143.     """
  144.     函数说明:岭回归测试
  145.     Parameters:
  146.         xMat - x数据集
  147.         yMat - y数据集
  148.     Returns:
  149.         wMat - 回归系数矩阵
  151.     """
  152.     xMat = np.mat(xArr); yMat = np.mat(yArr).T
  153.     #数据标准化
  154.     yMean = np.mean(yMat, axis = 0)                        #行与行操作，求均值
  155.     yMat = yMat - yMean                                    #数据减去均值
  156.     xMeans = np.mean(xMat, axis = 0)                    #行与行操作，求均值
  157.     xVar = np.var(xMat, axis = 0)                        #行与行操作，求方差
  158.     xMat = (xMat - xMeans) / xVar                        #数据减去均值除以方差实现标准化
  159.     numTestPts = 30                                        #30个不同的lambda测试
  160.     wMat = np.zeros((numTestPts, np.shape(xMat)[1]))    #初始回归系数矩阵
  161.     for i in range(numTestPts):                            #改变lambda计算回归系数
  162.         ws = ridgeRegres(xMat, yMat, np.exp(i - 10))    #lambda以e的指数变化，最初是一个非常小的数，
  163.         wMat[i, :] = ws.T                                 #计算回归系数矩阵
  164.     return wMat
  166. if __name__ == '__main__':
  167.     lgX = []
  168.     lgY = []
  169.     setDataCollect(lgX, lgY)
  170.     print(ridgeTest(lgX, lgY))

运行结果如下图所示：

24.png

看运行结果的第一行，可以看到最大的是第4项，第二大的是第2项。

因此，如果只选择一个特征来做预测的话，我们应该选择第4个特征，也就是原始加个。如果可以选择2个特征的话，应该选择第4个和第2个特征。

这种分析方法使得我们可以挖掘大量数据的内在规律。在仅有4个特征时，该方法的效果也许并不明显；但如果有100个以上的特征，该方法就会变得十分有效：它可以指出哪个特征是关键的，而哪些特征是不重要的。

八、使用Sklearn的linear_model

最后让用sklearn实现下岭回归吧。

官方英文文档地址：点击查看

sklearn.linear_model提供了很多线性模型，包括岭回归、贝叶斯回归、Lasso等。本文主要讲解岭回归Ridge。

25.png

1、Ridge
让我们先看下Ridge这个函数，一共有8个参数：

26.png

参数说明如下：

  • alpha：正则化系数，float类型，默认为1.0。正则化改善了问题的条件并减少了估计的方差。较大的值指定较强的正则化。
  • fit_intercept：是否需要截距，bool类型，默认为True。也就是是否求解b。
  • normalize：是否先进行归一化，bool类型，默认为False。如果为真，则回归X将在回归之前被归一化。 当fit_intercept设置为False时，将忽略此参数。 当回归量归一化时，注意到这使得超参数学习更加鲁棒，并且几乎不依赖于样本的数量。 相同的属性对标准化数据无效。然而，如果你想标准化，请在调用normalize = False训练估计器之前，使用preprocessing.StandardScaler处理数据。
  • copy_X：是否复制X数组，bool类型，默认为True，如果为True，将复制X数组; 否则，它覆盖原数组X。
  • max_iter：最大的迭代次数，int类型，默认为None，最大的迭代次数，对于sparse_cg和lsqr而言，默认次数取决于scipy.sparse.linalg，对于sag而言，则默认为1000次。
  • tol：精度，float类型，默认为0.001。就是解的精度。
  • solver：求解方法，str类型，默认为auto。可选参数为：auto、svd、cholesky、lsqr、sparse_cg、sag。
    • auto根据数据类型自动选择求解器。
    • svd使用X的奇异值分解来计算Ridge系数。对于奇异矩阵比cholesky更稳定。
    • cholesky使用标准的scipy.linalg.solve函数来获得闭合形式的解。
    • sparse_cg使用在scipy.sparse.linalg.cg中找到的共轭梯度求解器。作为迭代算法，这个求解器比大规模数据（设置tol和max_iter的可能性）的cholesky更合适。
    • lsqr使用专用的正则化最小二乘常数scipy.sparse.linalg.lsqr。它是最快的，但可能在旧的scipy版本不可用。它是使用迭代过程。
    • sag使用随机平均梯度下降。它也使用迭代过程，并且当n_samples和n_feature都很大时，通常比其他求解器更快。注意，sag快速收敛仅在具有近似相同尺度的特征上被保证。您可以使用sklearn.preprocessing的缩放器预处理数据。
  • random_state：sag的伪随机种子。
    以上就是所有的初始化参数，当然，初始化后还可以通过set_params方法重新进行设定。

知道了这些，接下来就可以编写代码了：

  1. # -*-coding:utf-8 -*-
  2. import numpy as np
  3. from bs4 import BeautifulSoup
  4. import random
  6. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  7.     """
  8.     函数说明:从页面读取数据，生成retX和retY列表
  9.     Parameters:
  10.         retX - 数据X
  11.         retY - 数据Y
  12.         inFile - HTML文件
  13.         yr - 年份
  14.         numPce - 乐高部件数目
  15.         origPrc - 原价
  16.     Returns:
  17.         无
  19.     """
  20.     # 打开并读取HTML文件
  21.     with open(inFile, encoding='utf-8') as f:
  22.         html = f.read()
  23.     soup = BeautifulSoup(html)
  24.     i = 1
  25.     # 根据HTML页面结构进行解析
  26.     currentRow = soup.find_all('table', r = "%d" % i)
  27.     while(len(currentRow) != 0):
  28.         currentRow = soup.find_all('table', r = "%d" % i)
  29.         title = currentRow[0].find_all('a')[1].text
  30.         lwrTitle = title.lower()
  31.         # 查找是否有全新标签
  32.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  33.             newFlag = 1.0
  34.         else:
  35.             newFlag = 0.0
  36.         # 查找是否已经标志出售，我们只收集已出售的数据
  37.         soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
  38.         if len(soldUnicde) == 0:
  39.             print("商品 #%d 没有出售" % i)
  40.         else:
  41.             # 解析页面获取当前价格
  42.             soldPrice = currentRow[0].find_all('td')[4]
  43.             priceStr = soldPrice.text
  44.             priceStr = priceStr.replace('$','')
  45.             priceStr = priceStr.replace(',','')
  46.             if len(soldPrice) > 1:
  47.                 priceStr = priceStr.replace('Free shipping', '')
  48.             sellingPrice = float(priceStr)
  49.             # 去掉不完整的套装价格
  50.             if  sellingPrice > origPrc * 0.5:
  51.                 print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
  52.                 retX.append([yr, numPce, newFlag, origPrc])
  53.                 retY.append(sellingPrice)
  54.         i += 1
  55.         currentRow = soup.find_all('table', r = "%d" % i)
  57. def setDataCollect(retX, retY):
  58.     """
  59.     函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
  60.     Parameters:
  61.         无
  62.     Returns:
  63.         无
  65.     """
  66.     scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)                #2006年的乐高8288,部件数目800,原价49.99
  67.     scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)                #2002年的乐高10030,部件数目3096,原价269.99
  68.     scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)                #2007年的乐高10179,部件数目5195,原价499.99
  69.     scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)                #2007年的乐高10181,部件数目3428,原价199.99
  70.     scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)                #2008年的乐高10189,部件数目5922,原价299.99
  71.     scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)                #2009年的乐高10196,部件数目3263,原价249.99
  73. def usesklearn():
  74.     """
  75.     函数说明:使用sklearn
  76.     Parameters:
  77.         无
  78.     Returns:
  79.         无
  81.     """
  82.     from sklearn import linear_model
  83.     reg = linear_model.Ridge(alpha = .5)
  84.     lgX = []
  85.     lgY = []
  86.     setDataCollect(lgX, lgY)
  87.     reg.fit(lgX, lgY)
  88.     print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3]))
  90. if __name__ == '__main__':
  91.     usesklearn()

运行结果如下图所示：

27.png

我们不搞太复杂，正则化项系数设为0.5，其余参数使用默认即可。可以看到，获得的结果与上小结的结果类似。

九、总结

  • 训练的模型要在测试集比较它们的效果，而不是在训练集上。
  • 在局部加权线性回归中，过小的核可能导致过拟合现象，即训练集表现良好，测试集表现很差。
  • 与分类一样，回归也是预测目标值的过程。回归与分类的不同点在于，前者预测连续类型变量，而后者预测离散类型变量。
  • 岭回归是缩减法的一种，相当于对回归系数的大小施加了限制。另一种很好的缩减法是lasso。lasso难以求解，但可以使用计算简便的逐步线性回归方法求的近似解。
  • 缩减法还可以看做是对一个模型增加偏差的同时减少方法。

参考资料

  • 本文中提到的使用预测鲍鱼年龄和玩具售价，来自于《机器学习实战》的第八章线性回归。
  • 本文的理论部分，参考《机器学习实战》的第八章回归。
  • 本文出现的所有代码，均可在github上下载，欢迎Follow、Star：Github地址