第8章 预测数值型数据：回归

预测数值型数据回归首页

回归（Regression） 概述

我们前边提到的分类的目标变量是标称型数据，而回归则是对连续型的数据做出处理，回归的目的是预测数值型数据的目标值。

回归 场景

回归的目的是预测数值型的目标值。最直接的办法是依据输入写出一个目标值的计算公式。

假如你想要预测兰博基尼跑车的功率大小，可能会这样计算:

HorsePower = 0.0015 annualSalary - 0.99 hoursListeningToPublicRadio

这就是所谓的 回归方程(regression equation)，其中的 0.0015 和 -0.99 称作 回归系数（regression weights），求这些回归系数的过程就是回归。一旦有了这些回归系数，再给定输入，做预测就非常容易了。具体的做法是用回归系数乘以输入值，再将结果全部加在一起，就得到了预测值。我们这里所说的，回归系数是一个向量，输入也是向量，这些运算也就是求出二者的内积。

说到回归，一般都是指 线性回归(linear regression)。线性回归意味着可以将输入项分别乘以一些常量，再将结果加起来得到输出。

补充：
线性回归假设特征和结果满足线性关系。其实线性关系的表达能力非常强大，每个特征对结果的影响强弱可以由前面的参数体现，而且每个特征变量可以首先映射到一个函数，然后再参与线性计算。这样就可以表达特征与结果之间的非线性关系。

回归 原理

1、线性回归

我们应该怎样从一大堆数据里求出回归方程呢？ 假定输入数据存放在矩阵 x 中，而回归系数存放在向量 w 中。那么对于给定的数据 X1，预测结果将会通过 Y = X1^T w 给出。现在的问题是，手里有一些 X 和对应的 y，怎样才能找到 w 呢？一个常用的方法就是找出使误差最小的 w 。这里的误差是指预测 y 值和真实 y 值之间的差值，使用该误差的简单累加将使得正差值和负差值相互抵消，所以我们采用平方误差（实际上就是我们通常所说的最小二乘法）。

平方误差可以写做（其实我们是使用这个函数作为 loss function）:

平方误差

用矩阵表示还可以写做 平方误差_2 。如果对 w 求导，得到 平方误差_3 ，令其等于零，解出 w 如下（具体求导过程为: http://blog.csdn.net/nomadlx53/article/details/50849941 ）:

回归系数的最佳估计计算公式

1.1、线性回归 须知概念

1.1.1、矩阵求逆

因为我们在计算回归方程的回归系数时，用到的计算公式如下:

回归系数的最佳估计计算公式

需要对矩阵求逆，因此这个方程只在逆矩阵存在的时候适用，我们在程序代码中对此作出判断。
判断矩阵是否可逆的一个可选方案是:

判断矩阵的行列式是否为 0，若为 0 ，矩阵就不存在逆矩阵，不为 0 的话，矩阵才存在逆矩阵。

1.1.2、最小二乘法

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

1.2、线性回归 工作原理

  1. 读入数据，将数据特征x、特征标签y存储在矩阵x、y中
  2. 验证 x^Tx 矩阵是否可逆
  3. 使用最小二乘法求得 回归系数 w 的最佳估计

1.3、线性回归 开发流程

  1. 收集数据: 采用任意方法收集数据
  2. 准备数据: 回归需要数值型数据，标称型数据将被转换成二值型数据
  3. 分析数据: 绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合线绘在图上作为对比
  4. 训练算法: 找到回归系数
  5. 测试算法: 使用 R^2 或者预测值和数据的拟合度，来分析模型的效果
  6. 使用算法: 使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续型数据而不仅仅是离散的类别标签

1.4、线性回归 算法特点

  1. 优点：结果易于理解，计算上不复杂。
  2. 缺点：对非线性的数据拟合不好。
  3. 适用于数据类型：数值型和标称型数据。

1.5、线性回归 项目案例

1.5.1、线性回归 项目概述

根据下图中的点，找出该数据的最佳拟合直线。

线性回归数据示例图

数据格式为: 

  1. x0          x1          y 
  2. 1.000000    0.067732    3.176513
  3. 1.000000    0.427810    3.816464
  4. 1.000000    0.995731    4.550095
  5. 1.000000    0.738336    4.256571
1.5.2、线性回归 编写代码
  1. def loadDataSet(fileName):                 
  2.     """ 加载数据
  3.         解析以tab键分隔的文件中的浮点数
  4.     Returns：
  5.         dataMat ：  feature 对应的数据集
  6.         labelMat ： feature 对应的分类标签，即类别标签
  8.     """
  9.     # 获取样本特征的总数，不算最后的目标变量 
  10.     numFeat = len(open(fileName).readline().split('\t')) - 1 
  11.     dataMat = []
  12.     labelMat = []
  13.     fr = open(fileName)
  14.     for line in fr.readlines():
  15.         # 读取每一行
  16.         lineArr =[]
  17.         # 删除一行中以tab分隔的数据前后的空白符号
  18.         curLine = line.strip().split('\t')
  19.         # i 从0到2，不包括2 
  20.         for i in range(numFeat):
  21.             # 将数据添加到lineArr List中，每一行数据测试数据组成一个行向量           
  22.             lineArr.append(float(curLine[i]))
  23.             # 将测试数据的输入数据部分存储到dataMat 的List中
  24.         dataMat.append(lineArr)
  25.         # 将每一行的最后一个数据，即类别，或者叫目标变量存储到labelMat List中
  26.         labelMat.append(float(curLine[-1]))
  27.     return dataMat,labelMat
  30. def standRegres(xArr,yArr):
  31.     '''
  32.     Description：
  33.         线性回归
  34.     Args:
  35.         xArr ：输入的样本数据，包含每个样本数据的 feature
  36.         yArr ：对应于输入数据的类别标签，也就是每个样本对应的目标变量
  37.     Returns:
  38.         ws：回归系数
  39.     '''
  41.     # mat()函数将xArr，yArr转换为矩阵 mat().T 代表的是对矩阵进行转置操作
  42.     xMat = mat(xArr)
  43.     yMat = mat(yArr).T
  44.     # 矩阵乘法的条件是左矩阵的列数等于右矩阵的行数
  45.     xTx = xMat.T*xMat
  46.     # 因为要用到xTx的逆矩阵，所以事先需要确定计算得到的xTx是否可逆，条件是矩阵的行列式不为0
  47.     # linalg.det() 函数是用来求得矩阵的行列式的，如果矩阵的行列式为0，则这个矩阵是不可逆的，就无法进行接下来的运算                   
  48.     if linalg.det(xTx) == 0.0:
  49.         print "This matrix is singular, cannot do inverse" 
  50.         return
  51.     # 最小二乘法
  52.     # http://www.apache.wiki/pages/viewpage.action?pageId=5505133
  53.     # 书中的公式，求得w的最优解
  54.     ws = xTx.I * (xMat.T*yMat)            
  55.     return ws
  58. def regression1():
  59.     xArr, yArr = loadDataSet("input/8.Regression/data.txt")
  60.     xMat = mat(xArr)
  61.     yMat = mat(yArr)
  62.     ws = standRegres(xArr, yArr)
  63.     fig = plt.figure()
  64.     ax = fig.add_subplot(111)               #add_subplot(349)函数的参数的意思是，将画布分成3行4列图像画在从左到右从上到下第9块
  65.     ax.scatter(xMat[:, 1].flatten(), yMat.T[:, 0].flatten().A[0]) #scatter 的x是xMat中的第二列，y是yMat的第一列
  66.     xCopy = xMat.copy() 
  67.     xCopy.sort(0)
  68.     yHat = xCopy * ws
  69.     ax.plot(xCopy[:, 1], yHat)
  70.     plt.show()

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

1.5.3、线性回归 拟合效果

线性回归数据效果图

2、局部加权线性回归

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

一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在这个算法中，我们给预测点附近的每个点赋予一定的权重，然后与 线性回归 类似，在这个子集上基于最小均方误差来进行普通的回归。我们需要最小化的目标函数大致为:

局部加权线性回归回归系数公式

与 kNN 一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解出回归系数 w 的形式如下:

局部加权线性回归回归系数公式

其中 w 是一个矩阵，用来给每个数据点赋予权重。

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

局部加权线性回归高斯核

这样就构建了一个只含对角元素的权重矩阵 w，并且点 x 与 x(i) 越近，w(i, i) 将会越大。上述公式中包含一个需要用户指定的参数 k ，它决定了对附近的点赋予多大的权重，这也是使用 LWLR 时唯一需要考虑的参数，下面的图给出了参数 k 与权重的关系。

参数k与权重的关系

上面的图是 每个点的权重图（假定我们正预测的点是 x = 0.5），最上面的图是原始数据集，第二个图显示了当 k = 0.5 时，大部分的数据都用于训练回归模型；而最下面的图显示当 k=0.01 时，仅有很少的局部点被用于训练回归模型。

2.1、局部加权线性回归 工作原理

  1. 读入数据，将数据特征x、特征标签y存储在矩阵x、y中
  2. 利用高斯核构造一个权重矩阵 W，对预测点附近的点施加权重
  3. 验证 X^TWX 矩阵是否可逆
  4. 使用最小二乘法求得 回归系数 w 的最佳估计

2.2、局部加权线性回归 项目案例

2.2.1、局部加权线性回归 项目概述

我们仍然使用上面 线性回归 的数据集，对这些点进行一个 局部加权线性回归 的拟合。

局部加权线性回归数据示例图

数据格式为: 

  1. 1.000000    0.067732    3.176513
  2. 1.000000    0.427810    3.816464
  3. 1.000000    0.995731    4.550095
  4. 1.000000    0.738336    4.256571
2.2.2、局部加权线性回归 编写代码
  1.     # 局部加权线性回归
  2. def lwlr(testPoint,xArr,yArr,k=1.0):
  3.     '''
  4.         Description：
  5.             局部加权线性回归，在待预测点附近的每个点赋予一定的权重，在子集上基于最小均方差来进行普通的回归。
  6.         Args：
  7.             testPoint：样本点
  8.             xArr：样本的特征数据，即 feature
  9.             yArr：每个样本对应的类别标签，即目标变量
  10.             k:关于赋予权重矩阵的核的一个参数，与权重的衰减速率有关
  11.         Returns:
  12.             testPoint * ws：数据点与具有权重的系数相乘得到的预测点
  13.         Notes:
  14.             这其中会用到计算权重的公式，w = e^((x^((i))-x) / -2k^2)
  15.             理解：x为某个预测点，x^((i))为样本点，样本点距离预测点越近，贡献的误差越大（权值越大），越远则贡献的误差越小（权值越小）。
  16.             关于预测点的选取，在我的代码中取的是样本点。其中k是带宽参数，控制w（钟形函数）的宽窄程度，类似于高斯函数的标准差。
  17.             算法思路：假设预测点取样本点中的第i个样本点（共m个样本点），遍历1到m个样本点（含第i个），算出每一个样本点与预测点的距离，
  18.             也就可以计算出每个样本贡献误差的权值，可以看出w是一个有m个元素的向量（写成对角阵形式）。
  19.     '''
  20.     # mat() 函数是将array转换为矩阵的函数， mat().T 是转换为矩阵之后，再进行转置操作
  21.     xMat = mat(xArr)
  22.     yMat = mat(yArr).T
  23.     # 获得xMat矩阵的行数
  24.     m = shape(xMat)[0]
  25.     # eye()返回一个对角线元素为1，其他元素为0的二维数组，创建权重矩阵weights，该矩阵为每个样本点初始化了一个权重                   
  26.     weights = mat(eye((m)))
  27.     for j in range(m):
  28.         # testPoint 的形式是 一个行向量的形式
  29.         # 计算 testPoint 与输入样本点之间的距离，然后下面计算出每个样本贡献误差的权值
  30.         diffMat = testPoint - xMat[j,:]
  31.         # k控制衰减的速度
  32.         weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))
  33.     # 根据矩阵乘法计算 xTx ，其中的 weights 矩阵是样本点对应的权重矩阵
  34.     xTx = xMat.T * (weights * xMat)
  35.     if linalg.det(xTx) == 0.0:
  36.         print ("This matrix is singular, cannot do inverse")
  37.         return
  38.     # 计算出回归系数的一个估计
  39.     ws = xTx.I * (xMat.T * (weights * yMat))
  40.     return testPoint * ws
  42. def lwlrTest(testArr,xArr,yArr,k=1.0):
  43.     '''
  44.         Description：
  45.             测试局部加权线性回归，对数据集中每个点调用 lwlr() 函数
  46.         Args：
  47.             testArr：测试所用的所有样本点
  48.             xArr：样本的特征数据，即 feature
  49.             yArr：每个样本对应的类别标签，即目标变量
  50.             k：控制核函数的衰减速率
  51.         Returns：
  52.             yHat：预测点的估计值
  53.     '''
  54.     # 得到样本点的总数
  55.     m = shape(testArr)[0]
  56.     # 构建一个全部都是 0 的 1 * m 的矩阵
  57.     yHat = zeros(m)
  58.     # 循环所有的数据点，并将lwlr运用于所有的数据点 
  59.     for i in range(m):
  60.         yHat[i] = lwlr(testArr[i],xArr,yArr,k)
  61.     # 返回估计值
  62.     return yHat
  64. def lwlrTestPlot(xArr,yArr,k=1.0):  
  65.     '''
  66.         Description:
  67.             首先将 X 排序，其余的都与lwlrTest相同，这样更容易绘图
  68.         Args：
  69.             xArr：样本的特征数据，即 feature
  70.             yArr：每个样本对应的类别标签，即目标变量，实际值
  71.             k：控制核函数的衰减速率的有关参数，这里设定的是常量值 1
  72.         Return：
  73.             yHat：样本点的估计值
  74.             xCopy：xArr的复制
  75.     '''
  76.     # 生成一个与目标变量数目相同的 0 向量
  77.     yHat = zeros(shape(yArr))
  78.     # 将 xArr 转换为 矩阵形式
  79.     xCopy = mat(xArr)
  80.     # 排序
  81.     xCopy.sort(0)
  82.     # 开始循环，为每个样本点进行局部加权线性回归，得到最终的目标变量估计值
  83.     for i in range(shape(xArr)[0]):
  84.         yHat[i] = lwlr(xCopy[i],xArr,yArr,k)
  85.     return yHat,xCopy
  88. #test for LWLR
  89. def regression2():
  90.     xArr, yArr = loadDataSet("input/8.Regression/data.txt")
  91.     yHat = lwlrTest(xArr, xArr, yArr, 0.003)
  92.     xMat = mat(xArr)
  93.     srtInd = xMat[:,1].argsort(0)           #argsort()函数是将x中的元素从小到大排列，提取其对应的index(索引)，然后输出
  94.     xSort=xMat[srtInd][:,0,:]
  95.     fig = plt.figure()
  96.     ax = fig.add_subplot(111)
  97.     ax.plot(xSort[:,1], yHat[srtInd])
  98.     ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0] , s=2, c='red')
  99.     plt.show()

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

2.2.3、局部加权线性回归 拟合效果

局部加权线性回归数据效果图

上图使用了 3 种不同平滑值绘出的局部加权线性回归的结果。上图中的平滑系数 k =1.0，中图 k = 0.01，下图 k = 0.003 。可以看到，k = 1.0 时的模型效果与最小二乘法差不多，k=0.01时该模型可以挖出数据的潜在规律，而 k=0.003时则考虑了太多的噪声，进而导致了过拟合现象。

2.3、局部加权线性回归 注意事项

局部加权线性回归也存在一个问题，即增加了计算量，因为它对每个点做预测时都必须使用整个数据集。

3、线性回归 & 局部加权线性回归 项目案例

到此为止，我们已经介绍了找出最佳拟合直线的两种方法，下面我们用这些技术来预测鲍鱼的年龄。

3.1、项目概述

我们有一份来自 UCI 的数据集合的数据，记录了鲍鱼（一种介壳类水生动物）的年龄。鲍鱼年龄可以从鲍鱼壳的层数推算得到。

3.2、开发流程

  1. 收集数据: 采用任意方法收集数据
  2. 准备数据: 回归需要数值型数据，标称型数据将被转换成二值型数据
  3. 分析数据: 绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合线绘在图上作为对比
  4. 训练算法: 找到回归系数
  5. 测试算法: 使用 rssError()函数 计算预测误差的大小，来分析模型的效果
  6. 使用算法: 使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续型数据而不仅仅是离散的类别标签

收集数据: 采用任意方法收集数据

准备数据: 回归需要数值型数据，标称型数据将被转换成二值型数据

数据存储格式:

  1. 1    0.455    0.365    0.095    0.514    0.2245    0.101    0.15    15
  2. 1    0.35    0.265    0.09    0.2255    0.0995    0.0485    0.07    7
  3. -1    0.53    0.42    0.135    0.677    0.2565    0.1415    0.21    9
  4. 1    0.44    0.365    0.125    0.516    0.2155    0.114    0.155    10
  5. 0    0.33    0.255    0.08    0.205    0.0895    0.0395    0.055    7

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

训练算法: 找到回归系数

使用上面我们讲到的 局部加权线性回归 训练算法，求出回归系数

测试算法: 使用 rssError()函数 计算预测误差的大小，来分析模型的效果

  1. # test for abloneDataSet
  2. def abaloneTest():
  3.     '''
  4.     Desc:
  5.         预测鲍鱼的年龄
  6.     Args:
  7.         None
  8.     Returns:
  9.         None
  10.     '''
  11.     # 加载数据
  12.     abX, abY = loadDataSet("input/8.Regression/abalone.txt")
  13.     # 使用不同的核进行预测
  14.     oldyHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)
  15.     oldyHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
  16.     oldyHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)   
  17.     # 打印出不同的核预测值与训练数据集上的真实值之间的误差大小
  18.     print "old yHat01 error Size is :" , rssError(abY[0:99], oldyHat01.T)
  19.     print "old yHat1 error Size is :" , rssError(abY[0:99], oldyHat1.T)
  20.     print "old yHat10 error Size is :" , rssError(abY[0:99], oldyHat10.T)
  22.     # 打印出 不同的核预测值 与 新数据集（测试数据集）上的真实值之间的误差大小
  23.     newyHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
  24.     print "new yHat01 error Size is :" , rssError(abY[0:99], newyHat01.T)
  25.     newyHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
  26.     print "new yHat1 error Size is :" , rssError(abY[0:99], newyHat1.T)
  27.     newyHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
  28.     print "new yHat10 error Size is :" , rssError(abY[0:99], newyHat10.T)
  30.     # 使用简单的 线性回归 进行预测，与上面的计算进行比较
  31.     standWs = standRegres(abX[0:99], abY[0:99])
  32.     standyHat = mat(abX[100:199]) * standWs
  33.     print "standRegress error Size is:", rssError(abY[100:199], standyHat.T.A)

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

根据我们上边的测试，可以看出:

简单线性回归达到了与局部加权现行回归类似的效果。这也说明了一点，必须在未知数据上比较效果才能选取到最佳模型。那么最佳的核大小是 10 吗？或许是，但如果想得到更好的效果，应该用 10 个不同的样本集做 10 次测试来比较结果。

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

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

如果数据的特征比样本点还多应该怎么办？是否还可以使用线性回归和之前的方法来做预测？答案是否定的，即我们不能再使用前面介绍的方法。这是因为在计算 矩阵求逆 的时候会出错。

如果特征比样本点还多(n > m)，也就是说输入数据的矩阵 x 不是满秩矩阵。非满秩矩阵求逆时会出现问题。

为了解决这个问题，我们引入了 岭回归（ridge regression） 这种缩减方法。接着是 lasso法，最后介绍 前向逐步回归。

4.1、岭回归

简单来说，岭回归就是在矩阵 矩阵_1 上加一个 λI 从而使得矩阵非奇异，进而能对 矩阵_2 求逆。其中矩阵I是一个 n * n （等于列数） 的单位矩阵，
对角线上元素全为1，其他元素全为0。而λ是一个用户定义的数值，后面会做介绍。在这种情况下，回归系数的计算公式将变成：

岭回归的回归系数计算

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

岭回归

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

这里通过预测误差最小化得到 λ: 数据获取之后，首先抽一部分数据用于测试，剩余的作为训练集用于训练参数 w。训练完毕后在测试集上测试预测性能。通过选取不同的 λ 来重复上述测试过程，最终得到一个使预测误差最小的 λ 。

4.1.1、岭回归 原始代码
  1. def ridgeRegres(xMat,yMat,lam=0.2):
  2.     '''
  3.         Desc：
  4.             这个函数实现了给定 lambda 下的岭回归求解。
  5.             如果数据的特征比样本点还多，就不能再使用上面介绍的的线性回归和局部线性回归了，因为计算 (xTx)^(-1)会出现错误。
  6.             如果特征比样本点还多（n > m），也就是说，输入数据的矩阵x不是满秩矩阵。非满秩矩阵在求逆时会出现问题。
  7.             为了解决这个问题，我们下边讲一下：岭回归，这是我们要讲的第一种缩减方法。
  8.         Args：
  9.             xMat：样本的特征数据，即 feature
  10.             yMat：每个样本对应的类别标签，即目标变量，实际值
  11.             lam：引入的一个λ值，使得矩阵非奇异
  12.         Returns：
  13.             经过岭回归公式计算得到的回归系数
  14.     '''
  16.     xTx = xMat.T*xMat
  17.     # 岭回归就是在矩阵 xTx 上加一个 λI 从而使得矩阵非奇异，进而能对 xTx + λI 求逆
  18.     denom = xTx + eye(shape(xMat)[1])*lam
  19.     # 检查行列式是否为零，即矩阵是否可逆，行列式为0的话就不可逆，不为0的话就是可逆。
  20.     if linalg.det(denom) == 0.0:
  21.         print ("This matrix is singular, cannot do inverse")
  22.         return
  23.     ws = denom.I * (xMat.T*yMat)
  24.     return ws
  27. def ridgeTest(xArr,yArr):
  28.     '''
  29.         Desc：
  30.             函数 ridgeTest() 用于在一组 λ 上测试结果
  31.         Args：
  32.             xArr：样本数据的特征，即 feature
  33.             yArr：样本数据的类别标签，即真实数据
  34.         Returns：
  35.             wMat：将所有的回归系数输出到一个矩阵并返回
  36.     '''
  38.     xMat = mat(xArr)
  39.     yMat=mat(yArr).T
  40.     # 计算Y的均值
  41.     yMean = mean(yMat,0)
  42.     # Y的所有的特征减去均值
  43.     yMat = yMat - yMean
  44.     # 标准化 x，计算 xMat 平均值
  45.     xMeans = mean(xMat,0)
  46.     # 然后计算 X的方差
  47.     xVar = var(xMat,0)
  48.     # 所有特征都减去各自的均值并除以方差
  49.     xMat = (xMat - xMeans)/xVar
  50.     # 可以在 30 个不同的 lambda 下调用 ridgeRegres() 函数。
  51.     numTestPts = 30
  52.     # 创建30 * m 的全部数据为0 的矩阵
  53.     wMat = zeros((numTestPts,shape(xMat)[1]))
  54.     for i in range(numTestPts):
  55.         # exp() 返回 e^x 
  56.         ws = ridgeRegres(xMat,yMat,exp(i-10))
  57.         wMat[i,:]=ws.T
  58.     return wMat
  61. #test for ridgeRegression
  62. def regression3():
  63.     abX,abY = loadDataSet("input/8.Regression/abalone.txt")
  64.     ridgeWeights = ridgeTest(abX, abY)
  65.     fig = plt.figure()
  66.     ax = fig.add_subplot(111)
  67.     ax.plot(ridgeWeights)
  68.     plt.show()

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

4.1.2、岭回归在鲍鱼数据集上的运行效果

岭回归的运行效果

上图绘制出了回归系数与 log(λ) 的关系。在最左边，即 λ 最小时，可以得到所有系数的原始值（与线性回归一致）；而在右边，系数全部缩减为0；在中间部分的某值将可以取得最好的预测效果。为了定量地找到最佳参数值，还需要进行交叉验证。另外，要判断哪些变量对结果预测最具有影响力，在上图中观察它们对应的系数大小就可以了。

4.2、套索方法(Lasso，The Least Absolute Shrinkage and Selection Operator)

在增加如下约束时，普通的最小二乘法回归会得到与岭回归一样的公式:

lasso_1

上式限定了所有回归系数的平方和不能大于 λ 。使用普通的最小二乘法回归在当两个或更多的特征相关时，可能会得到一个很大的正系数和一个很大的负系数。正是因为上述限制条件的存在，使用岭回归可以避免这个问题。

与岭回归类似，另一个缩减方法lasso也对回归系数做了限定，对应的约束条件如下:

lasso_2

唯一的不同点在于，这个约束条件使用绝对值取代了平方和。虽然约束形式只是稍作变化，结果却大相径庭: 在 λ 足够小的时候，一些系数会因此被迫缩减到 0.这个特性可以帮助我们更好地理解数据。

4.3、前向逐步回归

前向逐步回归算法可以得到与 lasso 差不多的效果，但更加简单。它属于一种贪心算法，即每一步都尽可能减少误差。一开始，所有权重都设置为 1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。

伪代码如下:

  1. 数据标准化，使其分布满足 0 均值 和单位方差
  2. 在每轮迭代过程中: 
  3.     设置当前最小误差 lowestError 为正无穷
  4.     对每个特征:
  5.         增大或缩小:
  6.             改变一个系数得到一个新的 w
  7.             计算新 w 下的误差
  8.             如果误差 Error 小于当前最小误差 lowestError: 设置 Wbest 等于当前的 W
  9.         将 W 设置为新的 Wbest
4.3.1、前向逐步回归 原始代码
  1. def stageWise(xArr,yArr,eps=0.01,numIt=100):
  2.     xMat = mat(xArr); yMat=mat(yArr).T
  3.     yMean = mean(yMat,0)
  4.     yMat = yMat - yMean     # 也可以规则化ys但会得到更小的coef
  5.     xMat = regularize(xMat)
  6.     m,n=shape(xMat)
  7.     #returnMat = zeros((numIt,n)) # 测试代码删除
  8.     ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()
  9.     for i in range(numIt):
  10.         print (ws.T)
  11.         lowestError = inf; 
  12.         for j in range(n):
  13.             for sign in [-1,1]:
  14.                 wsTest = ws.copy()
  15.                 wsTest[j] += eps*sign
  16.                 yTest = xMat*wsTest
  17.                 rssE = rssError(yMat.A,yTest.A)
  18.                 if rssE < lowestError:
  19.                     lowestError = rssE
  20.                     wsMax = wsTest
  21.         ws = wsMax.copy()
  22.         returnMat[i,:]=ws.T
  23.     return returnMat
  26. #test for stageWise
  27. def regression4():
  28.     xArr,yArr=loadDataSet("input/8.Regression/abalone.txt")
  29.     stageWise(xArr,yArr,0.01,200)
  30.     xMat = mat(xArr)
  31.     yMat = mat(yArr).T
  32.     xMat = regularize(xMat)
  33.     yM = mean(yMat,0)
  34.     yMat = yMat - yM
  35.     weights = standRegres(xMat, yMat.T)
  36.     print (weights.T)

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

4.3.2、逐步线性回归在鲍鱼数据集上的运行效果

逐步线性回归运行效果

逐步线性回归算法的主要优点在于它可以帮助人们理解现有的模型并作出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。最后，如果用于测试，该算法每100次迭代后就可以构建出一个模型，可以使用类似于10折交叉验证的方法比较这些模型，最终选择使误差最小的模型。

4.4、小结

当应用缩减方法（如逐步线性回归或岭回归）时，模型也就增加了偏差（bias），与此同时却减小了模型的方差。

5、回归 项目案例

项目案例1: 预测乐高玩具套装的价格

项目概述

Dangler 喜欢为乐高套装估价，我们用回归技术来帮助他建立一个预测模型。

开发流程
  1. (1) 收集数据：用 Google Shopping 的API收集数据。
  2. (2) 准备数据：从返回的JSON数据中抽取价格。
  3. (3) 分析数据：可视化并观察数据。
  4. (4) 训练算法：构建不同的模型，采用逐步线性回归和直接的线性回归模型。
  5. (5) 测试算法：使用交叉验证来测试不同的模型，分析哪个效果最好。
  6. (6) 使用算法：这次练习的目标就是生成数据模型。

收集数据: 使用 Google 购物的 API

由于 Google 提供的 api 失效，我们只能自己下载咯，将数据存储在了 input 文件夹下的 setHtml 文件夹下

准备数据: 从返回的 JSON 数据中抽取价格

因为我们这里不是在线的，就不再是 JSON 了，我们直接解析线下的网页，得到我们想要的数据。

分析数据: 可视化并观察数据

这里我们将解析得到的数据打印出来，然后观察数据。

训练算法: 构建不同的模型

  1. from numpy import *
  2. from bs4 import BeautifulSoup
  4. # 从页面读取数据，生成retX和retY列表
  5. def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
  7.     # 打开并读取HTML文件
  8.     fr = open(inFile)
  9.     soup = BeautifulSoup(fr.read())
  10.     i=1
  12.     # 根据HTML页面结构进行解析
  13.     currentRow = soup.findAll('table', r="%d" % i)
  14.     while(len(currentRow)!=0):
  15.         currentRow = soup.findAll('table', r="%d" % i)
  16.         title = currentRow[0].findAll('a')[1].text
  17.         lwrTitle = title.lower()
  19.         # 查找是否有全新标签
  20.         if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
  21.             newFlag = 1.0
  22.         else:
  23.             newFlag = 0.0
  25.         # 查找是否已经标志出售，我们只收集已出售的数据
  26.         soldUnicde = currentRow[0].findAll('td')[3].findAll('span')
  27.         if len(soldUnicde)==0:
  28.             print "item #%d did not sell" % i
  29.         else:
  30.             # 解析页面获取当前价格
  31.             soldPrice = currentRow[0].findAll('td')[4]
  32.             priceStr = soldPrice.text
  33.             priceStr = priceStr.replace('$','') #strips out $
  34.             priceStr = priceStr.replace(',','') #strips out ,
  35.             if len(soldPrice)>1:
  36.                 priceStr = priceStr.replace('Free shipping', '')
  37.             sellingPrice = float(priceStr)
  39.             # 去掉不完整的套装价格
  40.             if  sellingPrice > origPrc * 0.5:
  41.                     print "%d\t%d\t%d\t%f\t%f" % (yr,numPce,newFlag,origPrc, sellingPrice)
  42.                     retX.append([yr, numPce, newFlag, origPrc])
  43.                     retY.append(sellingPrice)
  44.         i += 1
  45.         currentRow = soup.findAll('table', r="%d" % i)
  47. # 依次读取六种乐高套装的数据，并生成数据矩阵        
  48. def setDataCollect(retX, retY):
  49.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego8288.html', 2006, 800, 49.99)
  50.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego10030.html', 2002, 3096, 269.99)
  51.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego10179.html', 2007, 5195, 499.99)
  52.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego10181.html', 2007, 3428, 199.99)
  53.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego10189.html', 2008, 5922, 299.99)
  54.     scrapePage(retX, retY, 'input/8.Regression/setHtml/lego10196.html', 2009, 3263, 249.99)

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

  1. # 交叉验证测试岭回归
  2. def crossValidation(xArr,yArr,numVal=10):
  3.     # 获得数据点个数，xArr和yArr具有相同长度
  4.     m = len(yArr)
  5.     indexList = range(m)
  6.     errorMat = zeros((numVal,30))
  8.     # 主循环 交叉验证循环
  9.     for i in range(numVal):
  10.         # 随机拆分数据，将数据分为训练集（90%）和测试集（10%）
  11.         trainX=[]; trainY=[]
  12.         testX = []; testY = []
  14.         # 对数据进行混洗操作
  15.         random.shuffle(indexList)
  17.         # 切分训练集和测试集
  18.         for j in range(m):
  19.             if j < m*0.9: 
  20.                 trainX.append(xArr[indexList[j]])
  21.                 trainY.append(yArr[indexList[j]])
  22.             else:
  23.                 testX.append(xArr[indexList[j]])
  24.                 testY.append(yArr[indexList[j]])
  26.         # 获得回归系数矩阵
  27.         wMat = ridgeTest(trainX,trainY)
  29.         # 循环遍历矩阵中的30组回归系数
  30.         for k in range(30):
  31.             # 读取训练集和数据集
  32.             matTestX = mat(testX); matTrainX=mat(trainX)
  33.             # 对数据进行标准化
  34.             meanTrain = mean(matTrainX,0)
  35.             varTrain = var(matTrainX,0)
  36.             matTestX = (matTestX-meanTrain)/varTrain
  38.             # 测试回归效果并存储
  39.             yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)
  41.             # 计算误差
  42.             errorMat[i,k] = ((yEst.T.A-array(testY))**2).sum()
  44.     # 计算误差估计值的均值
  45.     meanErrors = mean(errorMat,0)
  46.     minMean = float(min(meanErrors))
  47.     bestWeights = wMat[nonzero(meanErrors==minMean)]
  49.     # 不要使用标准化的数据，需要对数据进行还原来得到输出结果
  50.     xMat = mat(xArr); yMat=mat(yArr).T
  51.     meanX = mean(xMat,0); varX = var(xMat,0)
  52.     unReg = bestWeights/varX
  54.     # 输出构建的模型
  55.     print "the best model from Ridge Regression is:\n",unReg
  56.     print "with constant term: ",-1*sum(multiply(meanX,unReg)) + mean(yMat)
  59. # predict for lego's price
  60. def regression5():
  61.     lgX = []
  62.     lgY = []
  64.     setDataCollect(lgX, lgY)
  65.     crossValidation(lgX, lgY, 10)

使用算法：这次练习的目标就是生成数据模型

完整代码地址: https://github.com/apachecn/MachineLearning/blob/master/src/python/8.PredictiveNumericalDataRegression/regression.py

6、附加 权衡偏差和方差

任何时候，一旦发现模型和测量值之间存在差异，就说出现了误差。当考虑模型中的 “噪声” 或者说误差时，必须考虑其来源。你可能会对复杂的过程进行简化，这将导致在模型和测量值之间出现 “噪声” 或误差，若无法理解数据的真实生成过程，也会导致差异的产生。另外，测量过程本身也可能产生 “噪声” 或者问题。下面我们举一个例子，我们使用 线性回归 和 局部加权线性回归 处理过一个从文件导入的二维数据。

生成公式

其中的 N(0, 1) 是一个均值为 0、方差为 1 的正态分布。我们尝试过禁用一条直线来拟合上述数据。不难想到，直线所能得到的最佳拟合应该是 3.0+1.7x 这一部分。这样的话，误差部分就是 0.1sin(30x)+0.06N(0, 1) 。在上面，我们使用了局部加权线性回归来试图捕捉数据背后的结构。该结构拟合起来有一定的难度，因此我们测试了多组不同的局部权重来找到具有最小测试误差的解。

下图给出了训练误差和测试误差的曲线图，上面的曲面就是测试误差，下面的曲线是训练误差。我们根据 预测鲍鱼年龄 的实验知道: 如果降低核的大小，那么训练误差将变小。从下图开看，从左到右就表示了核逐渐减小的过程。

偏差方差图

一般认为，上述两种误差由三个部分组成: 偏差、测量误差和随机噪声。局部加权线性回归 和 预测鲍鱼年龄 中，我们通过引入了三个越来越小的核来不断增大模型的方差。

在缩减系数来“理解”数据这一节中，我们介绍了缩减法，可以将一些系数缩减成很小的值或直接缩减为 0 ，这是一个增大模型偏差的例子。通过把一些特征的回归系数缩减到 0 ，同时也就减小了模型的复杂度。例子中有 8 个特征，消除其中两个后不仅使模型更易理解，同时还降低了预测误差。对照上图，左侧是参数缩减过于严厉的结果，而右侧是无缩减的效果。

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