机器学习实战教程（十二）：线性回归提升篇之乐高玩具套件二手价预测
原文连接：cuijiahua.com/blog/2017/1…php
1、前言
本篇文章讲解线性回归的缩减方法，岭回归以及逐步线性回归，同时熟悉sklearn的岭回归使用方法，对乐高玩具套件的二手价格作出预测。html
2、岭回归
若是数据的特征比样本点还多应该怎么办？很显然，此时咱们不能再使用上文的方法进行计算了，由于矩阵X不是满秩矩阵，非满秩矩阵在求逆时会出现问题。为了解决这个问题，统计学家引入岭回归（ridge regression）的概念。git
一、什么是岭回归？
岭回归即咱们所说的L2正则线性回归，在通常的线性回归最小化均方偏差的基础上增长了一个参数w的L2范数的罚项，从而最小化罚项残差平方和：github
简单说来，岭回归就是在普通线性回归的基础上引入单位矩阵。回归系数的计算公式变形以下：算法
式中，矩阵I是一个mxm的单位矩阵，加上一个λI从而使得矩阵非奇异，进而能对矩阵求逆。windows
岭回归最早用来处理特征数多于样本数的状况，如今也用于在估计中加入误差，从而获得更好的估计。这里经过引入λ来限制了全部w之和，经过引入该惩罚项，可以减小不重要的参数，这个技术在统计学中也能够叫作缩减（shrinkage）。api
缩减方法能够去掉不重要的参数，所以能更好地裂解数据。此外，与简单的线性回归相比，缩减法可以取得更好的预测效果。数组
二、编写代码
为了使用岭回归和缩减技术，首先须要对特征作标准化处理。由于，咱们须要使每一个维度特征具备相同的重要性。本文使用的标准化处理比较简单，就是将全部特征都减去各自的均值并除以方差。bash
代码很简单，只须要稍作修改，其中，λ为模型的参数。咱们先绘制一个回归系数与log(λ)的曲线图，看下它们的规律，编写代码以下：网络
# -*-coding:utf-8 -*-
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
def loadDataSet(fileName):
""" 函数说明:加载数据 Parameters: fileName - 文件名 Returns: xArr - x数据集 yArr - y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
numFeat = len(open(fileName).readline().split('\t')) - 1
xArr = []; 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, yArr
def ridgeRegres(xMat, yMat, lam = 0.2):
""" 函数说明:岭回归 Parameters: xMat - x数据集 yMat - y数据集 lam - 缩减系数 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
xTx = xMat.T * xMat
denom = xTx + np.eye(np.shape(xMat)[1]) * lam
if np.linalg.det(denom) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = denom.I * (xMat.T * yMat)
return ws
def ridgeTest(xArr, yArr):
""" 函数说明:岭回归测试 Parameters: xMat - x数据集 yMat - y数据集 Returns: wMat - 回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
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 wMat
def plotwMat():
""" 函数说明:绘制岭回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
abX, abY = loadDataSet('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__':
plotwMat()
复制代码
运行结果以下：
图绘制了回归系数与log(λ)的关系。在最左边，即λ最小时，能够获得全部系数的原始值（与线性回归一致）；而在右边，系数所有缩减成0；在中间部分的某个位置，将会获得最好的预测结果。想要获得最佳的λ参数，可使用交叉验证的方式得到，文章的后面会继续讲解。
3、前向逐步线性回归
前向逐步线性回归算法属于一种贪心算法，即每一步都尽量减小偏差。咱们计算回归系数，再也不是经过公式计算，而是经过每次微调各个回归系数，而后计算预测偏差。那个使偏差最小的一组回归系数，就是咱们须要的最佳回归系数。
前向逐步线性回归实现也很简单。固然，仍是先进行数据标准化，编写代码以下：
# -*-coding:utf-8 -*-
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
def loadDataSet(fileName):
""" 函数说明:加载数据 Parameters: fileName - 文件名 Returns: xArr - x数据集 yArr - y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
numFeat = len(open(fileName).readline().split('\t')) - 1
xArr = []; 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, yArr
def regularize(xMat, yMat):
""" 函数说明:数据标准化 Parameters: xMat - x数据集 yMat - y数据集 Returns: inxMat - 标准化后的x数据集 inyMat - 标准化后的y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-11-23 """
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, inyMat
def rssError(yArr,yHatArr):
""" 函数说明:计算平方偏差 Parameters: yArr - 预测值 yHatArr - 真实值 Returns: Website: https://www.cuijiahua.com/ Modify: 2017-11-23 """
return ((yArr-yHatArr)**2).sum()
def stageWise(xArr, yArr, eps = 0.01, numIt = 100):
""" 函数说明:前向逐步线性回归 Parameters: xArr - x输入数据 yArr - y预测数据 eps - 每次迭代须要调整的步长 numIt - 迭代次数 Returns: returnMat - numIt次迭代的回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
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 = rssE
wsMax = wsTest
ws = wsMax.copy()
returnMat[i,:] = ws.T #记录numIt次迭代的回归系数矩阵
return returnMat
def plotstageWiseMat():
""" 函数说明:绘制岭回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
xArr, yArr = loadDataSet('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，同时也就减小了模型的复杂度。消除了多余的特征以后，模型更容易理解，同时也下降了预测偏差。可是当缩减过于严厉的时候，就会出现过拟合的现象，即用训练集预测结果很好，用测试集预测就糟糕不少。
4、预测乐高玩具套件的价格
乐高（LEGO）公司生产拼装类玩具，由不少大小不一样的塑料插块组成。它的样子以下图所示：
通常来讲，这些插块都是成套出售，它们能够拼装成不少不一样的东西，如船、城堡、一些著名建筑等。乐高公司每一个套装包含的部件数目从10件到5000件不等。
一种乐高套件基本上在几年后就会停产，但乐高的收藏者之间仍会在停产后彼此交易。本次实例，就是使用回归方法对收藏者之间的交易价格进行预测。
一、获取数据
书中使用的方法是经过Google提供的API进行获取数据，可是如今这个API已经关闭，咱们没法经过api获取数据了。不过幸运的是，我在网上找到了书上用到的那些html文件。
原始数据下载地址：数据下载
咱们经过解析html文件，来获取咱们须要的信息，若是学过个人《Python3网络爬虫》，那么我想这部分的内容会很是简单，解析代码以下：
# -*-coding:utf-8 -*-
from bs4 import BeautifulSoup
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
""" 函数说明:从页面读取数据，生成retX和retY列表 Parameters: retX - 数据X retY - 数据Y inFile - HTML文件 yr - 年份 numPce - 乐高部件数目 origPrc - 原价 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
# 打开并读取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].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
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.text
priceStr = 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 += 1
currentRow = soup.find_all('table', r = "%d" % i)
def setDataCollect(retX, retY):
""" 函数说明:依次读取六种乐高套装的数据，并生成数据矩阵 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99
if __name__ == '__main__':
lgX = []
lgY = []
setDataCollect(lgX, lgY)
复制代码
运行结果以下：
咱们对没有的商品作了处理。这些特征分别为：出品年份、部件数目、是否为全新、原价、售价（二手交易）。
html解析页面不会使用？那就学习一下爬虫知识吧~！若是对此不感兴趣，也能够跳过获取数据和解析数据，这个过程，看成已知数据，继续进行下一步。
二、创建模型
咱们已经处理好了数据集，接下来就是训练模型。首先咱们须要添加全为0的特征X0列。由于线性回归的第一列特征要求都是1.0。而后使用最简单的普通线性回归i，编写代码以下：
# -*-coding:utf-8 -*-
import numpy as np
from bs4 import BeautifulSoup
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
""" 函数说明:从页面读取数据，生成retX和retY列表 Parameters: retX - 数据X retY - 数据Y inFile - HTML文件 yr - 年份 numPce - 乐高部件数目 origPrc - 原价 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
# 打开并读取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].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
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.text
priceStr = 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 += 1
currentRow = soup.find_all('table', r = "%d" % i)
#
def setDataCollect(retX, retY):
""" 函数说明:依次读取六种乐高套装的数据，并生成数据矩阵 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99
def regularize(xMat, yMat):
""" 函数说明:数据标准化 Parameters: xMat - x数据集 yMat - y数据集 Returns: inxMat - 标准化后的x数据集 inyMat - 标准化后的y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
inxMat = xMat.copy() #数据拷贝
inyMat = yMat.copy()
yMean = np.mean(yMat, 0) #行与行操做，求均值
inyMat = yMat - yMean #数据减去均值
inMeans = np.mean(inxMat, 0) #行与行操做，求均值
inVar = np.var(inxMat, 0) #行与行操做，求方差
# print(inxMat)
print(inMeans)
# print(inVar)
inxMat = (inxMat - inMeans) / inVar #数据减去均值除以方差实现标准化
return inxMat, inyMat
def rssError(yArr,yHatArr):
""" 函数说明:计算平方偏差 Parameters: yArr - 预测值 yHatArr - 真实值 Returns: Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
return ((yArr-yHatArr)**2).sum()
def standRegres(xArr,yArr):
""" 函数说明:计算回归系数w Parameters: xArr - x数据集 yArr - y数据集 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-12 """
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T * xMat #根据文中推导的公示计算回归系数
if np.linalg.det(xTx) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = xTx.I * (xMat.T*yMat)
return ws
def useStandRegres():
""" 函数说明:使用简单的线性回归 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-11-12 """
lgX = []
lgY = []
setDataCollect(lgX, lgY)
data_num, features_num = np.shape(lgX)
lgX1 = np.mat(np.ones((data_num, features_num + 1)))
lgX1[:, 1:5] = np.mat(lgX)
ws = standRegres(lgX1, lgY)
print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0],ws[1],ws[2],ws[3],ws[4]))
if __name__ == '__main__':
useStandRegres()
复制代码
运行结果以下图所示：
能够看到，模型采用的公式如上图所示。虽然这个模型对于数据拟合得很好，可是看上不没有什么道理。套件里的部件数量越多，售价反而下降了，这是不合理的。
咱们使用岭回归，经过交叉验证，找到使偏差最小的λ对应的回归系数。编写代码以下：
# -*-coding:utf-8 -*-
import numpy as np
from bs4 import BeautifulSoup
import random
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
""" 函数说明:从页面读取数据，生成retX和retY列表 Parameters: retX - 数据X retY - 数据Y inFile - HTML文件 yr - 年份 numPce - 乐高部件数目 origPrc - 原价 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
# 打开并读取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].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
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.text
priceStr = 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 += 1
currentRow = soup.find_all('table', r = "%d" % i)
def ridgeRegres(xMat, yMat, lam = 0.2):
""" 函数说明:岭回归 Parameters: xMat - x数据集 yMat - y数据集 lam - 缩减系数 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
xTx = xMat.T * xMat
denom = xTx + np.eye(np.shape(xMat)[1]) * lam
if np.linalg.det(denom) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = denom.I * (xMat.T * yMat)
return ws
def setDataCollect(retX, retY):
""" 函数说明:依次读取六种乐高套装的数据，并生成数据矩阵 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99
def regularize(xMat, yMat):
""" 函数说明:数据标准化 Parameters: xMat - x数据集 yMat - y数据集 Returns: inxMat - 标准化后的x数据集 inyMat - 标准化后的y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
inxMat = xMat.copy() #数据拷贝
inyMat = yMat.copy()
yMean = np.mean(yMat, 0) #行与行操做，求均值
inyMat = yMat - yMean #数据减去均值
inMeans = np.mean(inxMat, 0) #行与行操做，求均值
inVar = np.var(inxMat, 0) #行与行操做，求方差
# print(inxMat)
print(inMeans)
# print(inVar)
inxMat = (inxMat - inMeans) / inVar #数据减去均值除以方差实现标准化
return inxMat, inyMat
def rssError(yArr,yHatArr):
""" 函数说明:计算平方偏差 Parameters: yArr - 预测值 yHatArr - 真实值 Returns: Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
return ((yArr-yHatArr)**2).sum()
def standRegres(xArr,yArr):
""" 函数说明:计算回归系数w Parameters: xArr - x数据集 yArr - y数据集 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-12 """
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T * xMat #根据文中推导的公示计算回归系数
if np.linalg.det(xTx) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = xTx.I * (xMat.T*yMat)
return ws
def crossValidation(xArr, yArr, numVal = 10):
""" 函数说明:交叉验证岭回归 Parameters: xArr - x数据集 yArr - y数据集 numVal - 交叉验证次数 Returns: wMat - 回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
m = len(yArr) #统计样本个数
indexList = list(range(m)) #生成索引值列表
errorMat = np.zeros((numVal,30)) #create error mat 30columns numVal rows
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个不一样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).T
meanX = 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]))
def ridgeTest(xArr, yArr):
""" 函数说明:岭回归测试 Parameters: xMat - x数据集 yMat - y数据集 Returns: wMat - 回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
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 wMat
if __name__ == '__main__':
lgX = []
lgY = []
setDataCollect(lgX, lgY)
crossValidation(lgX, lgY)
复制代码
运行结果以下图所示：
这里随机选取样本，由于其随机性，因此每次运行的结果可能略有不一样。不过总体如上图所示，能够看出，它与常规的最小二乘法，即普通的线性回归没有太大差别。咱们本指望找到一个更易于理解的模型，显然没有达到预期效果。
如今，咱们看一下在缩减过程当中回归系数是如何变化的。编写代码以下：
# -*-coding:utf-8 -*-
import numpy as np
from bs4 import BeautifulSoup
import random
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
""" 函数说明:从页面读取数据，生成retX和retY列表 Parameters: retX - 数据X retY - 数据Y inFile - HTML文件 yr - 年份 numPce - 乐高部件数目 origPrc - 原价 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
# 打开并读取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].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
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.text
priceStr = 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 += 1
currentRow = soup.find_all('table', r = "%d" % i)
def ridgeRegres(xMat, yMat, lam = 0.2):
""" 函数说明:岭回归 Parameters: xMat - x数据集 yMat - y数据集 lam - 缩减系数 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
xTx = xMat.T * xMat
denom = xTx + np.eye(np.shape(xMat)[1]) * lam
if np.linalg.det(denom) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = denom.I * (xMat.T * yMat)
return ws
def setDataCollect(retX, retY):
""" 函数说明:依次读取六种乐高套装的数据，并生成数据矩阵 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99
def regularize(xMat, yMat):
""" 函数说明:数据标准化 Parameters: xMat - x数据集 yMat - y数据集 Returns: inxMat - 标准化后的x数据集 inyMat - 标准化后的y数据集 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
inxMat = xMat.copy() #数据拷贝
inyMat = yMat.copy()
yMean = np.mean(yMat, 0) #行与行操做，求均值
inyMat = yMat - yMean #数据减去均值
inMeans = np.mean(inxMat, 0) #行与行操做，求均值
inVar = np.var(inxMat, 0) #行与行操做，求方差
# print(inxMat)
print(inMeans)
# print(inVar)
inxMat = (inxMat - inMeans) / inVar #数据减去均值除以方差实现标准化
return inxMat, inyMat
def rssError(yArr,yHatArr):
""" 函数说明:计算平方偏差 Parameters: yArr - 预测值 yHatArr - 真实值 Returns: Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
return ((yArr-yHatArr)**2).sum()
def standRegres(xArr,yArr):
""" 函数说明:计算回归系数w Parameters: xArr - x数据集 yArr - y数据集 Returns: ws - 回归系数 Website: https://www.cuijiahua.com/ Modify: 2017-11-12 """
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T * xMat #根据文中推导的公示计算回归系数
if np.linalg.det(xTx) == 0.0:
print("矩阵为奇异矩阵,不能转置")
return
ws = xTx.I * (xMat.T*yMat)
return ws
def ridgeTest(xArr, yArr):
""" 函数说明:岭回归测试 Parameters: xMat - x数据集 yMat - y数据集 Returns: wMat - 回归系数矩阵 Website: https://www.cuijiahua.com/ Modify: 2017-11-20 """
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 wMat
if __name__ == '__main__':
lgX = []
lgY = []
setDataCollect(lgX, lgY)
print(ridgeTest(lgX, lgY))
复制代码
运行结果以下图所示：
看运行结果的第一行，能够看到最大的是第4项，第二大的是第2项。
所以，若是只选择一个特征来作预测的话，咱们应该选择第4个特征，也就是原始加个。若是能够选择2个特征的话，应该选择第4个和第2个特征。
这种分析方法使得咱们能够挖掘大量数据的内在规律。在仅有4个特征时，该方法的效果也许并不明显；但若是有100个以上的特征，该方法就会变得十分有效：它能够指出哪一个特征是关键的，而哪些特征是不重要的。
5、使用Sklearn的linear_model
老规矩，最后让咱们用sklearn实现下岭回归吧。
官方英文文档地址：[点击查看]（cuijiahua.com/blog/tag/sk…
sklearn.linear_model提供了不少线性模型，包括岭回归、贝叶斯回归、Lasso等。本文主要讲解岭回归Ridge。
一、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 np
from bs4 import BeautifulSoup
import random
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
""" 函数说明:从页面读取数据，生成retX和retY列表 Parameters: retX - 数据X retY - 数据Y inFile - HTML文件 yr - 年份 numPce - 乐高部件数目 origPrc - 原价 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
# 打开并读取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].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
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.text
priceStr = 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 += 1
currentRow = soup.find_all('table', r = "%d" % i)
def setDataCollect(retX, retY):
""" 函数说明:依次读取六种乐高套装的数据，并生成数据矩阵 Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-03 """
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99
def usesklearn():
""" 函数说明:使用sklearn Parameters: 无 Returns: 无 Website: https://www.cuijiahua.com/ Modify: 2017-12-08 """
from sklearn import linear_model
reg = 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，其他参数使用默认便可。能够看到，得到的结果与上小结的结果相似。
6、总结
与分类同样，回归也是预测目标值的过程。回归与分类的不一样点在于，前者预测连续类型变量，然后者预测离散类型变量。 岭回归是缩减法的一种，至关于对回归系数的大小施加了限制。另外一种很好的缩减法是lasso。lasso难以求解，但可使用计算简便的逐步线性回归方法求的近似解。 缩减法还能够看作是对一个模型增长误差的同时减小方法。 下篇文章讲解树回归。
PS： 若是以为本篇本章对您有所帮助，欢迎关注、评论、赞！
本文出现的全部代码和数据集，都可在个人github上下载，欢迎Follow、Star：点击查看
参考文献：
《机器学习实战》的第五章内容。
相关文章和视频推荐
圆方圆学院聚集 Python + AI 名师，打造精品的 Python + AI 技术课程。 在各大平台都长期有优质免费公开课，欢迎报名收看。 公开课地址： ke.qq.com/course/3627…
相关文章
- 1. 机器学习实战_线性回归&逻辑回归（二）
- 2. 机器学习入门之房价预测（线性回归）
- 3. 机器学习之线性回归第二节线性回归
- 4. 机器学习之预测回归一：简单线性回归
- 5. 机器学习--线性回归算法预测房价
- 6. 机器学习做业-线性回归 南京房价预测
- 7. 机器学习全套教程（十六）-- 线性回归的改进-岭回归
- 8. python机器学习（二）线性回归
- 9. 《机器学习实战》-线性回归
- 10. Python3机器学习实践：线性回归【实例：波士顿房价预测】
- 更多相关文章...
- • Hibernate整合EHCache实现二级缓存 - Hibernate教程
- • C# 多线程 - C#教程
- • Kotlin学习（二）基本类型
- • Docker容器实战(七) - 容器眼光下的文件系统
相关标签/搜索
每日一句
-
每一个你不满意的现在，都有一个你没有努力的曾经。
欢迎关注本站公众号,获取更多信息
