一、前言

本篇文章讲解线性回归的缩减方法，岭回归以及逐步线性回归，同时熟悉sklearn的岭回归使用方法，对乐高玩具套件的二手价格做出预测。

二、岭回归

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

（可逆矩阵是满秩矩阵、不可逆矩阵（奇异矩阵）是降秩矩阵）

1、什么是岭回归？

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

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

2、编写代码

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

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

from numpy import *import matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesdef loadDataSet(fileName):"""加载数据集:param fileName: 文件名称:return: 返回特征数据和标签数据"""numFeat = len(open(fileName).readline().split("\t")) - 1 # 返回特征值dataMat = []labelMat = []fr = open(fileName)for line in fr.readlines():lineArr = []curLine = line.strip().split("\t")for i in range(numFeat):lineArr.append(float(curLine[i]))dataMat.append(lineArr)labelMat.append(float(curLine[-1]))return dataMat, labelMatdef ridgeRegres(xMat, yMat, lam = 0.2):"""函数说明:岭回归Parameters:xMat - x数据集yMat - y数据集lam - 缩减系数Returns:ws - 回归系数"""xTx = xMat.T * xMatdenom = xTx + eye(shape(xMat)[1]) * lamif linalg.det(denom) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = denom.I * (xMat.T * yMat)return wsdef ridgeTest(xArr, yArr):"""函数说明:岭回归测试Parameters:xMat - x数据集yMat - y数据集Returns:wMat - 回归系数矩阵"""xMat = mat(xArr) # 4177 x 8 矩阵yMat = mat(yArr).T # 4177 x 1 矩阵#数据标准化yMean = mean(yMat, axis = 0) #求均值 mean(axis=0) 其中 axis=0 求每列的均值 axis=1求没行的均值yMat = yMat - yMean #数据减去均值xMeans = mean(xMat, axis = 0) #求样本数据每个特征的均值xVar = var(xMat, axis = 0)#数据样本中每个特征，求方差xMat = (xMat - xMeans) / xVar #数据减去均值除以方差实现标准化numTestPts = 30 #30个不同的lambda测试wMat = zeros((numTestPts, shape(xMat)[1]))#初始回归系数矩阵 30 x 8 矩阵for i in range(numTestPts):#改变lambda计算回归系数ws = ridgeRegres(xMat, yMat, exp(i - 10)) #lambda以e的指数变化，最初是一个非常小的数，wMat[i, :] = ws.T #计算回归系数矩阵return wMatdef plotwMat():"""函数说明:绘制岭回归系数矩阵Website:/Modify:-11-20"""font = FontProperties(fname=r"./font/simsun.ttc", size=14)abX, abY = loadDataSet('./data/abalone.txt')redgeWeights = ridgeTest(abX, abY)fig = plt.figure()ax = fig.add_subplot(111)ax.plot(redgeWeights)ax_title_text = ax.set_title(u'log(lambada)与回归系数的关系', FontProperties = font)ax_xlabel_text = ax.set_xlabel(u'log(lambada)', FontProperties = font)ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties = font)plt.setp(ax_title_text, size = 20, weight = 'bold', color = 'red')plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black')plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black')plt.show()if __name__ == '__main__':xArr, yArr = loadDataSet("./data/abalone.txt")# 岭回归plotwMat()

运行结果如下：

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

三、前向逐步线性回归

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

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

# -*-coding:utf-8 -*-from matplotlib.font_manager import FontPropertiesimport matplotlib.pyplot as pltimport numpy as npdef loadDataSet(fileName):"""函数说明:加载数据Parameters:fileName - 文件名Returns:xArr - x数据集yArr - y数据集"""numFeat = len(open(fileName).readline().split('\t')) - 1xArr = []; yArr = []fr = open(fileName)for line in fr.readlines():lineArr =[]curLine = line.strip().split('\t')for i in range(numFeat):lineArr.append(float(curLine[i]))xArr.append(lineArr)yArr.append(float(curLine[-1]))return xArr, yArrdef regularize(xMat, yMat):"""函数说明:数据标准化Parameters:xMat - x数据集yMat - y数据集Returns:inxMat - 标准化后的x数据集inyMat - 标准化后的y数据集""" inxMat = xMat.copy()#数据拷贝inyMat = yMat.copy()yMean = np.mean(yMat, 0) #行与行操作，求均值inyMat = yMat - yMean#数据减去均值inMeans = np.mean(inxMat, 0) #行与行操作，求均值inVar = np.var(inxMat, 0) #行与行操作，求方差inxMat = (inxMat - inMeans) / inVar #数据减去均值除以方差实现标准化return inxMat, inyMatdef rssError(yArr,yHatArr):"""函数说明:计算平方误差Parameters:yArr - 预测值yHatArr - 真实值Returns:"""return ((yArr-yHatArr)**2).sum()def stageWise(xArr, yArr, eps = 0.01, numIt = 100):"""函数说明:前向逐步线性回归Parameters:xArr - x输入数据yArr - y预测数据eps - 每次迭代需要调整的步长numIt - 迭代次数Returns:returnMat - numIt次迭代的回归系数矩阵"""xMat = np.mat(xArr); yMat = np.mat(yArr).T #数据集xMat, yMat = regularize(xMat, yMat) #数据标准化m, n = np.shape(xMat)returnMat = np.zeros((numIt, n)) #初始化numIt次迭代的回归系数矩阵ws = np.zeros((n, 1))#初始化回归系数矩阵wsTest = ws.copy()wsMax = ws.copy()for i in range(numIt):#迭代numIt次# print(ws.T) #打印当前回归系数矩阵lowestError = float('inf'); #正无穷for j in range(n):#遍历每个特征的回归系数for sign in [-1, 1]:wsTest = ws.copy()wsTest[j] += eps * sign #微调回归系数yTest = xMat * wsTest #计算预测值rssE = rssError(yMat.A, yTest.A)#计算平方误差if rssE < lowestError: #如果误差更小，则更新当前的最佳回归系数lowestError = rssEwsMax = wsTestws = wsMax.copy()returnMat[i,:] = ws.T #记录numIt次迭代的回归系数矩阵return returnMatdef plotstageWiseMat():"""函数说明:绘制岭回归系数矩阵"""font = FontProperties(fname=r"./font/simsun.ttc", size=14)xArr, yArr = loadDataSet('./data/abalone.txt')returnMat = stageWise(xArr, yArr, 0.005, 1000)fig = plt.figure()ax = fig.add_subplot(111)ax.plot(returnMat) ax_title_text = ax.set_title(u'前向逐步回归:迭代次数与回归系数的关系', FontProperties = font)ax_xlabel_text = ax.set_xlabel(u'迭代次数', FontProperties = font)ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties = font)plt.setp(ax_title_text, size = 15, weight = 'bold', color = 'red')plt.setp(ax_xlabel_text, size = 10, weight = 'bold', color = 'black')plt.setp(ax_ylabel_text, size = 10, weight = 'bold', color = 'black')plt.show()if __name__ == '__main__':plotstageWiseMat()

运行结果如下：

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

总结一下：

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

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

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

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

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

1、获取数据

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

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

我们通过解析html文件，来获取我们需要的信息，如果学过我的《Python3网络爬虫》，那么我想这部分的内容会非常简单，解析代码如下：

# -*-coding:utf-8 -*-from bs4 import BeautifulSoupdef scrapePage(retX, retY, inFile, yr, numPce, origPrc):"""函数说明:从页面读取数据，生成retX和retY列表Parameters:retX - 数据XretY - 数据YinFile - HTML文件yr - 年份numPce - 乐高部件数目origPrc - 原价Returns:无"""# 打开并读取HTML文件with open(inFile, encoding='utf-8') as f:html = f.read()soup = BeautifulSoup(html)i = 1# 根据HTML页面结构进行解析currentRow = soup.find_all('table', r = "%d" % i)while(len(currentRow) != 0):currentRow = soup.find_all('table', r = "%d" % i)title = currentRow[0].find_all('a')[1].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else:newFlag = 0.0# 查找是否已经标志出售，我们只收集已出售的数据soldUnicde = currentRow[0].find_all('td')[3].find_all('span')if len(soldUnicde) == 0:print("商品 #%d 没有出售" % i)else:# 解析页面获取当前价格soldPrice = currentRow[0].find_all('td')[4]priceStr = soldPrice.textpriceStr = priceStr.replace('$','')priceStr = priceStr.replace(',','')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = soup.find_all('table', r = "%d" % i)def setDataCollect(retX, retY):"""函数说明:依次读取六种乐高套装的数据，并生成数据矩阵Parameters:无Returns:无"""scrapePage(retX, retY, './lego/lego8288.html', , 800, 49.99)#的乐高8288,部件数目800,原价49.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)#2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', , 5195, 499.99)#的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', , 3428, 199.99)#的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', , 5922, 299.99)#的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', , 3263, 249.99)#的乐高10196,部件数目3263,原价249.99if __name__ == '__main__':lgX = []lgY = []setDataCollect(lgX, lgY)

运行结果如下：

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

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

2、建立模型

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

# -*-coding:utf-8 -*-from bs4 import BeautifulSoupimport numpy as npdef scrapePage(retX, retY, inFile, yr, numPce, origPrc):"""函数说明:从页面读取数据，生成retX和retY列表Parameters:retX - 数据XretY - 数据YinFile - HTML文件yr - 年份numPce - 乐高部件数目origPrc - 原价Returns:无"""# 打开并读取HTML文件with open(inFile, encoding='utf-8') as f:html = f.read()soup = BeautifulSoup(html)i = 1# 根据HTML页面结构进行解析currentRow = soup.find_all('table', r="%d" % i)while (len(currentRow) != 0):currentRow = soup.find_all('table', r="%d" % i)title = currentRow[0].find_all('a')[1].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else:newFlag = 0.0# 查找是否已经标志出售，我们只收集已出售的数据soldUnicde = currentRow[0].find_all('td')[3].find_all('span')if len(soldUnicde) == 0:print("商品 #%d 没有出售" % i)else:# 解析页面获取当前价格soldPrice = currentRow[0].find_all('td')[4]priceStr = soldPrice.textpriceStr = priceStr.replace('$', '')priceStr = priceStr.replace(',', '')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = soup.find_all('table', r="%d" % i)def setDataCollect(retX, retY):"""函数说明:依次读取六种乐高套装的数据，并生成数据矩阵Parameters:无Returns:无"""scrapePage(retX, retY, './lego/lego8288.html', , 800, 49.99) # 的乐高8288,部件数目800,原价49.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) # 2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', , 5195, 499.99) # 的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', , 3428, 199.99) # 的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', , 5922, 299.99) # 的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', , 3263, 249.99) # 的乐高10196,部件数目3263,原价249.99def regularize(xMat, yMat):"""函数说明:数据标准化Parameters:xMat - x数据集yMat - y数据集Returns:inxMat - 标准化后的x数据集inyMat - 标准化后的y数据集"""inxMat = xMat.copy() # 数据拷贝yMean = np.mean(yMat, 0) # 行与行操作，求均值inyMat = yMat - yMean # 数据减去均值inMeans = np.mean(inxMat, 0) # 行与行操作，求均值inVar = np.var(inxMat, 0) # 行与行操作，求方差inxMat = (inxMat - inMeans) / inVar # 数据减去均值除以方差实现标准化return inxMat, inyMatdef rssError(yArr, yHatArr):"""函数说明:计算平方误差Parameters:yArr - 预测值yHatArr - 真实值Returns:"""return ((yArr - yHatArr) ** 2).sum()def standRegres(xArr, yArr):"""函数说明:计算回归系数wParameters:xArr - x数据集yArr - y数据集Returns:ws - 回归系数"""xMat = np.mat(xArr)yMat = np.mat(yArr).TxTx = xMat.T * xMat # 根据文中推导的公示计算回归系数if np.linalg.det(xTx) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = xTx.I * (xMat.T * yMat)return wsdef useStandRegres():"""函数说明:使用简单的线性回归Parameters:无Returns:无"""lgX = []lgY = []setDataCollect(lgX, lgY)data_num, features_num = np.shape(lgX) # （63 x 4）矩阵lgX1 = np.mat(np.ones((data_num, features_num + 1))) # 创建一个矩阵 63 x 5 矩阵 各元素值为1lgX1[:, 1:5] = np.mat(lgX) # 将 lgX 63 x 4 矩阵 变为 63 x 5 矩阵 矩阵第一列的值为1ws = standRegres(lgX1, lgY) # 求解w值print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0], ws[1], ws[2], ws[3], ws[4]))if __name__ == '__main__':lgX = []lgY = []useStandRegres()

运行结果如下图所示：

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

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

# -*-coding:utf-8 -*-from bs4 import BeautifulSoupimport numpy as npimport randomdef scrapePage(retX, retY, inFile, yr, numPce, origPrc):"""函数说明:从页面读取数据，生成retX和retY列表Parameters:retX - 数据XretY - 数据YinFile - HTML文件yr - 年份numPce - 乐高部件数目origPrc - 原价Returns:无"""# 打开并读取HTML文件with open(inFile, encoding='utf-8') as f:html = f.read()soup = BeautifulSoup(html)i = 1# 根据HTML页面结构进行解析currentRow = soup.find_all('table', r="%d" % i)while (len(currentRow) != 0):currentRow = soup.find_all('table', r="%d" % i)title = currentRow[0].find_all('a')[1].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else:newFlag = 0.0# 查找是否已经标志出售，我们只收集已出售的数据soldUnicde = currentRow[0].find_all('td')[3].find_all('span')if len(soldUnicde) == 0:print("商品 #%d 没有出售" % i)else:# 解析页面获取当前价格soldPrice = currentRow[0].find_all('td')[4]priceStr = soldPrice.textpriceStr = priceStr.replace('$', '')priceStr = priceStr.replace(',', '')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = soup.find_all('table', r="%d" % i)def setDataCollect(retX, retY):"""函数说明:依次读取六种乐高套装的数据，并生成数据矩阵Parameters:无Returns:无"""scrapePage(retX, retY, './lego/lego8288.html', , 800, 49.99) # 的乐高8288,部件数目800,原价49.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) # 2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', , 5195, 499.99) # 的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', , 3428, 199.99) # 的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', , 5922, 299.99) # 的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', , 3263, 249.99) # 的乐高10196,部件数目3263,原价249.99def regularize(xMat, yMat):"""函数说明:数据标准化Parameters:xMat - x数据集yMat - y数据集Returns:inxMat - 标准化后的x数据集inyMat - 标准化后的y数据集"""inxMat = xMat.copy() # 数据拷贝yMean = np.mean(yMat, 0) # 行与行操作，求均值inyMat = yMat - yMean # 数据减去均值inMeans = np.mean(inxMat, 0) # 行与行操作，求均值inVar = np.var(inxMat, 0) # 行与行操作，求方差inxMat = (inxMat - inMeans) / inVar # 数据减去均值除以方差实现标准化return inxMat, inyMatdef rssError(yArr, yHatArr):"""函数说明:计算平方误差Parameters:yArr - 预测值yHatArr - 真实值Returns:"""return ((yArr - yHatArr) ** 2).sum()def standRegres(xArr, yArr):"""函数说明:计算回归系数wParameters:xArr - x数据集yArr - y数据集Returns:ws - 回归系数"""xMat = np.mat(xArr)yMat = np.mat(yArr).TxTx = xMat.T * xMat # 根据文中推导的公示计算回归系数if np.linalg.det(xTx) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = xTx.I * (xMat.T * yMat)return wsdef useStandRegres():"""函数说明:使用简单的线性回归Parameters:无Returns:无"""lgX = []lgY = []setDataCollect(lgX, lgY)data_num, features_num = np.shape(lgX) # （63 x 4）矩阵lgX1 = np.mat(np.ones((data_num, features_num + 1))) # 创建一个矩阵 63 x 5 矩阵 各元素值为1lgX1[:, 1:5] = np.mat(lgX) # 将 lgX 63 x 4 矩阵 变为 63 x 5 矩阵 矩阵第一列的值为1ws = standRegres(lgX1, lgY) # 求解w值print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0], ws[1], ws[2], ws[3], ws[4]))def ridgeRegres(xMat, yMat, lam = 0.2):"""函数说明:岭回归系数Parameters:xMat - x数据集yMat - y数据集lam - 缩减系数Returns:ws - 回归系数"""xTx = xMat.T * xMatdenom = xTx + np.eye(np.shape(xMat)[1]) * lamif np.linalg.det(denom) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = denom.I * (xMat.T * yMat)return wsdef ridgeTest(xArr, yArr):"""函数说明:岭回归测试Parameters:xMat - x数据集yMat - y数据集Returns:wMat - 回归系数矩阵"""xMat = np.mat(xArr); yMat = np.mat(yArr).T#数据标准化yMean = np.mean(yMat, axis = 0) #行与行操作，求均值yMat = yMat - yMean#数据减去均值xMeans = np.mean(xMat, axis = 0)#行与行操作，求均值xVar = np.var(xMat, axis = 0) #行与行操作，求方差xMat = (xMat - xMeans) / xVar #数据减去均值除以方差实现标准化numTestPts = 30#30个不同的lambda测试wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #初始回归系数矩阵for i in range(numTestPts): #改变lambda计算回归系数ws = ridgeRegres(xMat, yMat, np.exp(i - 10)) #lambda以e的指数变化，最初是一个非常小的数，wMat[i, :] = ws.T #计算回归系数矩阵return wMatdef crossValidation(xArr, yArr, numVal = 10):"""函数说明:交叉验证岭回归Parameters:xArr - x数据集yArr - y数据集numVal - 交叉验证次数Returns:wMat - 回归系数矩阵"""m = len(yArr) #统计样本个数indexList = list(range(m)) #生成样本标签索引值列表errorMat = np.zeros((numVal,30))#create error mat 10 x 30 列for i in range(numVal): #交叉验证numVal次trainX = []; trainY = []#训练集testX = []; testY = []#测试集random.shuffle(indexList)#打乱次序for j in range(m): #划分数据集:90%训练集，10%测试集if j < m * 0.9:trainX.append(xArr[indexList[j]])trainY.append(yArr[indexList[j]])else:testX.append(xArr[indexList[j]])testY.append(yArr[indexList[j]])wMat = ridgeTest(trainX, trainY) # 30 x 4 矩阵系数#获得30个不同lambda下的岭回归系数for k in range(30): #遍历所有的岭回归系数matTestX = np.mat(testX); matTrainX = np.mat(trainX) #训练集meanTrain = np.mean(matTrainX,0) #训练集均值varTrain = np.var(matTrainX,0) #训练集方差matTestX = (matTestX - meanTrain) / varTrain #测试集标准化yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY) #根据ws预测y值errorMat[i, k] = rssError(yEst.T.A, np.array(testY)) #统计误差meanErrors = np.mean(errorMat,0) #计算每次交叉验证的平均误差minMean = float(min(meanErrors)) #找到最小误差bestWeights = wMat[np.nonzero(meanErrors == minMean)] #找到最佳回归系数xMat = np.mat(xArr); yMat = np.mat(yArr).TmeanX = np.mean(xMat,0); varX = np.var(xMat,0)unReg = bestWeights / varX#数据经过标准化，因此需要还原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]))if __name__ == '__main__':lgX = []lgY = []setDataCollect(lgX, lgY)crossValidation(lgX, lgY)

运行结果如下图所示：

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

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

# -*-coding:utf-8 -*-from bs4 import BeautifulSoupimport numpy as npimport randomdef scrapePage(retX, retY, inFile, yr, numPce, origPrc):"""函数说明:从页面读取数据，生成retX和retY列表Parameters:retX - 数据XretY - 数据YinFile - HTML文件yr - 年份numPce - 乐高部件数目origPrc - 原价Returns:无"""# 打开并读取HTML文件with open(inFile, encoding='utf-8') as f:html = f.read()soup = BeautifulSoup(html)i = 1# 根据HTML页面结构进行解析currentRow = soup.find_all('table', r="%d" % i)while (len(currentRow) != 0):currentRow = soup.find_all('table', r="%d" % i)title = currentRow[0].find_all('a')[1].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else:newFlag = 0.0# 查找是否已经标志出售，我们只收集已出售的数据soldUnicde = currentRow[0].find_all('td')[3].find_all('span')if len(soldUnicde) == 0:print("商品 #%d 没有出售" % i)else:# 解析页面获取当前价格soldPrice = currentRow[0].find_all('td')[4]priceStr = soldPrice.textpriceStr = priceStr.replace('$', '')priceStr = priceStr.replace(',', '')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = soup.find_all('table', r="%d" % i)def setDataCollect(retX, retY):"""函数说明:依次读取六种乐高套装的数据，并生成数据矩阵Parameters:无Returns:无"""scrapePage(retX, retY, './lego/lego8288.html', , 800, 49.99) # 的乐高8288,部件数目800,原价49.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) # 2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', , 5195, 499.99) # 的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', , 3428, 199.99) # 的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', , 5922, 299.99) # 的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', , 3263, 249.99) # 的乐高10196,部件数目3263,原价249.99def regularize(xMat, yMat):"""函数说明:数据标准化Parameters:xMat - x数据集yMat - y数据集Returns:inxMat - 标准化后的x数据集inyMat - 标准化后的y数据集"""inxMat = xMat.copy() # 数据拷贝yMean = np.mean(yMat, 0) # 行与行操作，求均值inyMat = yMat - yMean # 数据减去均值inMeans = np.mean(inxMat, 0) # 行与行操作，求均值inVar = np.var(inxMat, 0) # 行与行操作，求方差inxMat = (inxMat - inMeans) / inVar # 数据减去均值除以方差实现标准化return inxMat, inyMatdef rssError(yArr, yHatArr):"""函数说明:计算平方误差Parameters:yArr - 预测值yHatArr - 真实值Returns:"""return ((yArr - yHatArr) ** 2).sum()def standRegres(xArr, yArr):"""函数说明:计算回归系数wParameters:xArr - x数据集yArr - y数据集Returns:ws - 回归系数"""xMat = np.mat(xArr)yMat = np.mat(yArr).TxTx = xMat.T * xMat # 根据文中推导的公示计算回归系数if np.linalg.det(xTx) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = xTx.I * (xMat.T * yMat)return wsdef useStandRegres():"""函数说明:使用简单的线性回归Parameters:无Returns:无"""lgX = []lgY = []setDataCollect(lgX, lgY)data_num, features_num = np.shape(lgX) # （63 x 4）矩阵lgX1 = np.mat(np.ones((data_num, features_num + 1))) # 创建一个矩阵 63 x 5 矩阵 各元素值为1lgX1[:, 1:5] = np.mat(lgX) # 将 lgX 63 x 4 矩阵 变为 63 x 5 矩阵 矩阵第一列的值为1ws = standRegres(lgX1, lgY) # 求解w值print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0], ws[1], ws[2], ws[3], ws[4]))def ridgeRegres(xMat, yMat, lam = 0.2):"""函数说明:岭回归系数Parameters:xMat - x数据集yMat - y数据集lam - 缩减系数Returns:ws - 回归系数"""xTx = xMat.T * xMatdenom = xTx + np.eye(np.shape(xMat)[1]) * lamif np.linalg.det(denom) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = denom.I * (xMat.T * yMat)return wsdef ridgeTest(xArr, yArr):"""函数说明:岭回归测试Parameters:xMat - x数据集yMat - y数据集Returns:wMat - 回归系数矩阵"""xMat = np.mat(xArr); yMat = np.mat(yArr).T#数据标准化yMean = np.mean(yMat, axis = 0) #行与行操作，求均值yMat = yMat - yMean#数据减去均值xMeans = np.mean(xMat, axis = 0)#行与行操作，求均值xVar = np.var(xMat, axis = 0) #行与行操作，求方差xMat = (xMat - xMeans) / xVar #数据减去均值除以方差实现标准化numTestPts = 30#30个不同的lambda测试wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #初始回归系数矩阵for i in range(numTestPts): #改变lambda计算回归系数ws = ridgeRegres(xMat, yMat, np.exp(i - 10)) #lambda以e的指数变化，最初是一个非常小的数，wMat[i, :] = ws.T #计算回归系数矩阵return wMatdef crossValidation(xArr, yArr, numVal = 10):"""函数说明:交叉验证岭回归Parameters:xArr - x数据集yArr - y数据集numVal - 交叉验证次数Returns:wMat - 回归系数矩阵"""m = len(yArr) #统计样本个数indexList = list(range(m)) #生成样本标签索引值列表errorMat = np.zeros((numVal,30))#create error mat 10 x 30 列for i in range(numVal): #交叉验证numVal次trainX = []; trainY = []#训练集testX = []; testY = []#测试集random.shuffle(indexList)#打乱次序for j in range(m): #划分数据集:90%训练集，10%测试集if j < m * 0.9:trainX.append(xArr[indexList[j]])trainY.append(yArr[indexList[j]])else:testX.append(xArr[indexList[j]])testY.append(yArr[indexList[j]])wMat = ridgeTest(trainX, trainY) # 30 x 4 矩阵系数#获得30个不同lambda下的岭回归系数for k in range(30): #遍历所有的岭回归系数matTestX = np.mat(testX); matTrainX = np.mat(trainX) #训练集meanTrain = np.mean(matTrainX,0) #训练集均值varTrain = np.var(matTrainX,0) #训练集方差matTestX = (matTestX - meanTrain) / varTrain #测试集标准化yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY) #根据ws预测y值errorMat[i, k] = rssError(yEst.T.A, np.array(testY)) #统计误差meanErrors = np.mean(errorMat,0) #计算每次交叉验证的平均误差minMean = float(min(meanErrors)) #找到最小误差bestWeights = wMat[np.nonzero(meanErrors == minMean)] #找到最佳回归系数xMat = np.mat(xArr); yMat = np.mat(yArr).TmeanX = np.mean(xMat,0); varX = np.var(xMat,0)unReg = bestWeights / varX#数据经过标准化，因此需要还原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]))if __name__ == '__main__':lgX = []lgY = []setDataCollect(lgX, lgY)crossValidation(lgX, lgY)print(ridgeTest(lgX, lgY))

运行结果如下图所示：

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

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

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

五、使用Sklearn的linear_model

老规矩，最后让我们用sklearn实现下岭回归吧。

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

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

1、Ridge

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

参数说明如下：

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方法重新进行设定。

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

# -*-coding:utf-8 -*-import numpy as npfrom bs4 import BeautifulSoupimport randomdef scrapePage(retX, retY, inFile, yr, numPce, origPrc):"""函数说明:从页面读取数据，生成retX和retY列表Parameters:retX - 数据XretY - 数据YinFile - HTML文件yr - 年份numPce - 乐高部件数目origPrc - 原价Returns:无"""# 打开并读取HTML文件with open(inFile, encoding='utf-8') as f:html = f.read()soup = BeautifulSoup(html)i = 1# 根据HTML页面结构进行解析currentRow = soup.find_all('table', r = "%d" % i)while(len(currentRow) != 0):currentRow = soup.find_all('table', r = "%d" % i)title = currentRow[0].find_all('a')[1].textlwrTitle = title.lower()# 查找是否有全新标签if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):newFlag = 1.0else:newFlag = 0.0# 查找是否已经标志出售，我们只收集已出售的数据soldUnicde = currentRow[0].find_all('td')[3].find_all('span')if len(soldUnicde) == 0:print("商品 #%d 没有出售" % i)else:# 解析页面获取当前价格soldPrice = currentRow[0].find_all('td')[4]priceStr = soldPrice.textpriceStr = priceStr.replace('$','')priceStr = priceStr.replace(',','')if len(soldPrice) > 1:priceStr = priceStr.replace('Free shipping', '')sellingPrice = float(priceStr)# 去掉不完整的套装价格if sellingPrice > origPrc * 0.5:print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc])retY.append(sellingPrice)i += 1currentRow = soup.find_all('table', r = "%d" % i)def setDataCollect(retX, retY):"""函数说明:依次读取六种乐高套装的数据，并生成数据矩阵Parameters:无Returns:无"""scrapePage(retX, retY, './lego/lego8288.html', , 800, 49.99)#的乐高8288,部件数目800,原价49.99scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)#2002年的乐高10030,部件数目3096,原价269.99scrapePage(retX, retY, './lego/lego10179.html', , 5195, 499.99)#的乐高10179,部件数目5195,原价499.99scrapePage(retX, retY, './lego/lego10181.html', , 3428, 199.99)#的乐高10181,部件数目3428,原价199.99scrapePage(retX, retY, './lego/lego10189.html', , 5922, 299.99)#的乐高10189,部件数目5922,原价299.99scrapePage(retX, retY, './lego/lego10196.html', , 3263, 249.99)#的乐高10196,部件数目3263,原价249.99def usesklearn():"""函数说明:使用sklearnParameters:无Returns:无"""from sklearn import linear_modelreg = linear_model.Ridge(alpha = .5)lgX = []lgY = []setDataCollect(lgX, lgY)reg.fit(lgX, lgY)print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3])) if __name__ == '__main__':usesklearn()

运行结果如下图所示：

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

六、总结

与分类一样，回归也是预测目标值的过程。回归与分类的不同点在于，前者预测连续类型变量，而后者预测离散类型变量。岭回归是缩减法的一种，相当于对回归系数的大小施加了限制。另一种很好的缩减法是lasso。lasso难以求解，但可以使用计算简便的逐步线性回归方法求的近似解。缩减法还可以看做是对一个模型增加偏差的同时减少方法。下篇文章讲解树回归。

参考文献：

《机器学习实战》的第五章内容。

如果觉得《机器学习实战（七）02-线性回归提高篇之乐高玩具套件二手价预测》对你有帮助，请点赞、收藏，并留下你的观点哦！