Machine Learning in Action 读书笔记

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

文章目录

Machine Learning in Action 读书笔记
一、回归
- 1.回归的一般过程
- 2.线性回归
- 3.回归方程
- 4.回归系数
- 5.误差
- 6.普通最小二乘法
- 7.相关系数
二、回归函数
三、局部加权线性回归
四、缩减系数法
- 1.岭回归
- 2.lasso
- 3.前向逐步回归
五、全部代码

回归的目的是预测数值型的目标值。

一、回归

1.回归的一般过程

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

2.线性回归

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

3.回归方程

回归方程：是根据样本资料通过回归分析所得到的反映一个变量（因变量）对另一个或一组变量（自变量）的回归关系的数学表达式。回归直线方程用得比较多，可以用最小二乘法求回归直线方程中的a,b，从而得到回归直线方程。（百度百科）

4.回归系数

回归系数（regression coefficient）：在回归方程中表示自变量 x 对因变量 y 影响大小的参数。回归系数越大表示 x 对 y 影响越大，正回归系数表示 y 随 x 增大而增大，负回归系数表示y 随 x 增大而减小。例如回归方程式 Y=bX+a 中，斜率 b 称为回归系数，表示 X 每变动一单位，平均而言，Y 将变动 b 单位。（百度百科）

5.误差

误差：指预测 y 值和真实 y 值之间的差值，使用误差的简单累加将使得正差值和负差值相互抵消，所以该章采用平方误差。平方误差可以写做：
∑ i = 1 m ( y i − x i T w ) 2 \sum_{i=1}^m(y_i-x_i^Tw)^2 i=1∑m(yi−xiTw)2
用矩阵表示还可以写做：
( y − X w ) T ( y − X w ) (y-Xw)^T(y-Xw) (y−Xw)T(y−Xw)
对w求导，并令其等于0，解出：
w ^ = ( X T X ) − 1 X T y \hat{w}=(X^TX)^{-1}X^Ty w^=(XTX)−1XTy
w ^ \hat{w} w^表示这时当前估计出的w的最优解。

6.普通最小二乘法

普通最小二乘法(ordinary least squares，ols)：是回归分析（regression analysis）最根本的一个形式,对模型条件要求最少，也就是使散点图上的所有观测值到回归直线距离的平方和最小。（百度百科）

7.相关系数

计算预测值yHat和真实值y序列的匹配程度，就是计算这两个序列的相关系数。

二、回归函数

标准回归函数和数据导入函数

'''标准回归函数和数据导入函数'''
def loadDataSet(fileName):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 standRegres(xArr, yArr):  # 用于计算最佳拟合直线的回归系数矩阵xMat = mat(xArr); yMat = mat(yArr).TxTx = xMat.T*xMat  # 如果X的转置乘以X==0,那么X=0，无法求逆if linalg.det(xTx) == 0.0:  # linalg.det用于计算方阵的行列式值，行列式为0，方程不是满秩，不能计算逆print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T*yMat)  # xTx.I求矩阵的逆; wHat = （X^T*X）^(-1)*X^T*yreturn ws

使用标准回归函数计算出的回归系数矩阵，计算并绘出数据集散点图和最佳拟合直线图
其中ex0.txt部分数据集如下：

if __name__ == '__main__':xArr, yArr = loadDataSet('ex0.txt')# print(xArr[0:2])ws = standRegres(xArr, yArr)  # 计算回归系数# print(ws)xMat = mat(xArr)yMat = mat(yArr)yHat = xMat * wsfig = plt.figure()ax = fig.add_subplot(111)ax.scatter(xMat[:, 1].flatten().A[0], yMat.T[:, 0].flatten().A[0])xCopy = xMat.copy()xCopy.sort(0)  #axis=0，按列排列；axis=1，按行排列yHat = xCopy * wsax.plot(xCopy[:, 1], yHat)plt.show()

绘出的ex0.txt数据集散点图和它的最佳拟合直线如下所示：

三、局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。如果模型欠拟合，将不能取得最好的预测结果，所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。
局部加权线性回归（Locally Weighted Linear Regression, LWLR）：该算法为待预测点附近的每个点赋予一定的权重，该算法解出回归系数w的形式如下：
w ^ = ( X T W X ) X T W y \hat{w}=(X^TWX)X^TWy w^=(XTWX)XTWy
其中，W是一个矩阵，用来给每个数据点赋予权重。LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重，其中最常用的核就是高斯核，高斯核对应的权重如下：
w ( i , i ) = e x p ( ∣ x ( i ) − x ∣ − 2 k 2 ) w(i,i)=exp(\frac{|x^{(i)}-x|}{-2k^2}) w(i,i)=exp(−2k2∣x(i)−x∣)
其中，参数k决定了对附近的点赋予多大的权重。

'''局部加权线性回归函数'''
def lwlr(testPoint, xArr, yArr, k=1.0):  # k用于控制权重衰减的速度xMat = mat(xArr); yMat = mat(yArr).Tm = shape(xMat)[0]weights = mat(eye((m)))  # 创建对角权重矩阵for j in range(m):diffMat = testPoint - xMat[j, :]     #计算样本点和待预测点的距离差值weights[j, j] = exp(diffMat*diffMat.T/(-2.0*k**2)) #利用高斯核公式对附近的点赋予权重# print(weights)xTx = xMat.T * (weights * xMat)if linalg.det(xTx) == 0.0:print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T * (weights * yMat))return testPoint * wsdef lwlrTest(testArr, xArr, yArr, k=1.0):  #将局部加权线性回归应用到每一个数据点上；给定x空间中任意一点，计算出对应的预测值m = shape(testArr)[0]yHat = zeros(m)for i in range(m):yHat[i] = lwlr(testArr[i], xArr, yArr, k)return yHat

使用三种不同的平滑值绘出的局部加权线性回归结果

    yHat = lwlrTest(xArr, xArr, yArr, 1)# print(yHat)# 对xArr进行排序srtInd = xMat[:, 1].argsort(0)xSort = xMat[srtInd][:, 0, :]fig = plt.figure()ax = fig.add_subplot(111)ax.plot(xSort[:, 1], yHat[srtInd])   # 横坐标为样本点的第二个特征（第一个特征都是1），纵坐标为预测的y值ax.scatter(xMat[:, 1].flatten().A[0], mat(yArr).T.flatten().A[0], s=2, c='red')#横坐标相同，纵轴为真实y值plt.show()

其中，k = 1：欠拟合

k = 0.01：拟合较好

k = 0.003：考虑了太多噪声，导致了过拟合。

局部加权线性回归类似于kNN，每次预测均需要事先取出对应的数据子集，因此计算量较大。

上面介绍了两种找出最佳拟合直线的方法：简单最小二乘法和局部加权线性回归。

四、缩减系数法

如果特征比样本点多这种情况，那么输入数据的矩阵x不是满秩矩阵，非满秩矩阵在求逆时会出现问题，为了解决这个问题，统计学家引入了岭回归（ridge regression）的概念。岭回归是一种缩减方法，该章提到的缩减方法有：

岭回归
lasso法
前向逐步回归

缩减方法可以去掉不重要的参数，因此可以更好地理解数据。

1.岭回归

岭回归就是在矩阵 X T X X^TX XTX上加一个 λ I λI λI，从而使得矩阵非奇异，进而能够求逆。回归系数计算公式：
w ^ = ( X T X + λ I ) − 1 X T y \hat{w}=(X^TX+λI)^{-1}X^Ty w^=(XTX+λI)−1XTy
其中，λ是一个用户定义的数值，I是单位矩阵。
岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。通过引入λ来限制所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也叫做缩减。

'''岭回归'''
def ridgeRegres(xMat, yMat, lam=0.2):  # 用于计算回归系数，实现给定λ下的岭回归求解xTx = xMat.T * xMatdenom = xTx + eye(shape(xMat)[1]) * lamif linalg.det(denom) == 0.0:print("This matrix is singular, cannot do inverse")returnws = denom.I * (xMat.T * yMat)return wsdef ridgeTest(xArr, yArr): # 用于在一组λ上测试结果xMat = mat(xArr)yMat = mat(yArr).TyMean = mean(yMat, 0)yMat = yMat - yMean# 标准化xMeans = mean(xMat, 0)  # 计算平均值，axis=0表示按列计算xVar = var(xMat, 0)  # 计算方差，同上，0 表示按列计算xMat = (xMat - xMeans) / xVar   # 数据标准化，所有特征都减去各自的均值并除以方差numTestPts = 30   # 定义λ的个数，使用30个不同的λwMat = zeros((numTestPts, shape(xMat)[1]))for i in range(numTestPts):ws = ridgeRegres(xMat, yMat, exp(i - 10))  # exp(i - 10)表示lambda以指数级变化，这样方便看出λ取非常小的值和取非常大的值时分别对结果造成的影响wMat[i, :] = ws.Treturn wMat  # 将所有的回归系数输出到一个矩阵上并返回

上面代码中有数据标准化处理，具体做法是所有特征都减去各自的均值并除以方差。

#数据标准化，所有特征都减去各自的均值并除以方差(经常会用到这个操作，所以可以将这部分拉出来作为函数)
def regularize(xMat):inMat = xMat.copy()inMeans = mean(inMat,0)   #计算均值inVar = var(inMat,0)      #计算方差inMat = (inMat - inMeans)/inVarreturn inMat

下面是绘制的岭回归回归系数变化图，其中λ非常小时，系数与普通回归一样，而λ非常大时，所有回归系数缩减为0。可以在中间找到预测结果最好的λ值。

2.lasso

lasso法也对回归系数做出了限定：
∑ k = 1 n ∣ w k ∣ < = λ \sum_{k=1}^n|w_k|<=λ k=1∑n∣wk∣<=λ
为了在这个约束条件下解出回归系数，需要使用二次规划算法，计算比较复杂，使用前向逐步回归可以得到相同的结果，但是计算不复杂。

3.前向逐步回归

前向逐步回归算法属于一种贪心算法，即每一步都尽可能减少误差，一开始所有的权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。
该算法的伪代码如下所示：

数据标准化，使其分布满足0均值和单位方差
在每轮迭代过程中：设置当前最小误差lowestError为正无穷对每个特征：增大或缩小：改变一个系数得到一个新的w计算新w下的误差如果误差Error小于当前最小误差lowestError：设置Wbest等于当前的W将W设置为新的Wbest

'''逐步线性回归算法'''
def stageWise(xArr, yArr, eps=0.01, numIt=100): # 参数：输入数据，预测遍历，每次迭代需要调整的步长，迭代次数xMat = mat(xArr)  # 输入数据并存入矩阵yMat=mat(yArr).TyMean = mean(yMat,0)yMat = yMat - yMeanxMat = regularize(xMat) # 把特征按照均值为0方差为1进行标准化处理m, n=shape(xMat)# print(m, n)  # 4177 * 8returnMat = zeros((numIt, n)) #用于测试ws = zeros((n, 1))wsTest = ws.copy()  # 为了实现贪心算法，建立w的两份副本wsMax = ws.copy()for i in range(numIt): # 优化过程print(ws.T)lowestError = inffor j in range(n):for sign in [-1, 1]:wsTest = ws.copy()wsTest[j] += eps*signyTest = xMat*wsTestrssE = rssError(yMat.A, yTest.A)if rssE < lowestError:lowestError = rssEwsMax = wsTestws = wsMax.copy()returnMat[i, :]=ws.Treturn returnMat

逐步线性回归算法的主要优点在于它可以帮助人们理解现有模型并做出改进，该算法每100次迭代后就可以构建出一个模型，可以使用类似于10折交叉验证的方法比较这些模型，最终选择使误差最小的模型。
在使用缩减方法时，模型也就增加了偏差（bias），与此同时却减少了模型的方差。

五、全部代码

from numpy import *
import matplotlib.pyplot as plt'''标准回归函数和数据导入函数'''
def loadDataSet(fileName):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 standRegres(xArr, yArr):  # 用于计算最佳拟合直线xMat = mat(xArr); yMat = mat(yArr).TxTx = xMat.T*xMat  # 如果X的转置乘以X==0,那么X=0，无法求逆if linalg.det(xTx) == 0.0:  # linalg.det用于计算方阵的行列式值print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T*yMat)  # xTx.I求矩阵的逆; wHat = （X^T*X）^(-1)*X^T*yreturn ws'''局部加权线性回归函数'''
def lwlr(testPoint, xArr, yArr, k=1.0):  # k用于控制权重衰减的速度xMat = mat(xArr); yMat = mat(yArr).Tm = shape(xMat)[0]weights = mat(eye((m)))  # 创建对角权重矩阵for j in range(m):diffMat = testPoint - xMat[j, :]     #计算样本点和待预测点的距离差值weights[j, j] = exp(diffMat*diffMat.T/(-2.0*k**2)) #利用高斯核公式对附近的点赋予权重# print(weights)xTx = xMat.T * (weights * xMat)if linalg.det(xTx) == 0.0:print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T * (weights * yMat))return testPoint * wsdef lwlrTest(testArr, xArr, yArr, k=1.0):  #将局部加权线性回归应用到每一个数据点上；给定x空间中任意一点，计算出对应的预测值m = shape(testArr)[0]yHat = zeros(m)for i in range(m):yHat[i] = lwlr(testArr[i], xArr, yArr, k)return yHat'''示例：根据鲍鱼壳的层数预测鲍鱼的年龄'''
def rssError(yArr, yHatArr):return ((yArr - yHatArr)**2).sum()'''岭回归'''
def ridgeRegres(xMat, yMat, lam=0.2):  # 用于计算回归系数，实现给定λ下的岭回归求解xTx = xMat.T * xMatdenom = xTx + eye(shape(xMat)[1]) * lamif linalg.det(denom) == 0.0:print("This matrix is singular, cannot do inverse")returnws = denom.I * (xMat.T * yMat)return wsdef ridgeTest(xArr, yArr): # 用于在一组λ上测试结果xMat = mat(xArr)yMat = mat(yArr).TyMean = mean(yMat, 0)yMat = yMat - yMean# 标准化xMeans = mean(xMat, 0)  # 计算平均值，axis=0表示按列计算xVar = var(xMat, 0)  # 计算方差，同上，0 表示按列计算xMat = (xMat - xMeans) / xVar   # 数据标准化，所有特征都减去各自的均值并除以方差numTestPts = 30   # 定义λ的个数，使用30个不同的λwMat = zeros((numTestPts, shape(xMat)[1]))for i in range(numTestPts):ws = ridgeRegres(xMat, yMat, exp(i - 10))  # exp(i - 10)表示lambda以指数级变化，这样方便看出λ取非常小的值和取非常大的值时分别对结果造成的影响wMat[i, :] = ws.Treturn wMat  # 将所有的回归系数输出到一个矩阵上并返回#数据标准化，所有特征都减去各自的均值并除以方差(经常会用到这个操作，所以可以将这部分拉出来作为函数)
def regularize(xMat):inMat = xMat.copy()inMeans = mean(inMat,0)   #计算均值inVar = var(inMat,0)      #计算方差inMat = (inMat - inMeans)/inVarreturn inMat'''逐步线性回归算法'''
def stageWise(xArr, yArr, eps=0.01, numIt=100): # 参数：输入数据，预测遍历，每次迭代需要调整的步长，迭代次数xMat = mat(xArr)  # 输入数据并存入矩阵yMat=mat(yArr).TyMean = mean(yMat,0)yMat = yMat - yMeanxMat = regularize(xMat) # 把特征按照均值为0方差为1进行标准化处理m, n=shape(xMat)# print(m, n)  # 4177 * 8returnMat = zeros((numIt, n)) #用于测试ws = zeros((n, 1))wsTest = ws.copy()  # 为了实现贪心算法，建立w的两份副本wsMax = ws.copy()for i in range(numIt): # 优化过程print(ws.T)lowestError = inffor j in range(n):for sign in [-1, 1]:wsTest = ws.copy()wsTest[j] += eps*signyTest = xMat*wsTestrssE = rssError(yMat.A, yTest.A)if rssE < lowestError:lowestError = rssEwsMax = wsTestws = wsMax.copy()returnMat[i, :]=ws.Treturn returnMat'''示例：预测乐高玩具套装的价格'''
# 购物信息获取函数
from time import sleep
import json
import urllib
def searchForSet(retX, retY, setNum, yr, numPce, origPrc):#样本玩具特征矩阵，样本玩具的真实价格，获取样本的数量，样本玩具的年份，玩具套装的零件数，原始价格sleep(10)myAPIstr = 'get from code.google.com'searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum)# 字符串拼接print(searchURL)pg = urllib.request.urlopen(searchURL)   # 访问url地址retDict = json.loads(pg.read())  # 利用json打开和解析url获得数据，数据信息存入字典中for i in range(len(retDict['items'])): # 遍历数据的每一个条目try:currItem = retDict['items'][i]  # 获得当前条目if currItem['product']['condition'] == 'new': # 如果当前条目对应的产品是新产品newFlag = 1  # 标记为新else:newFlag = 0listOfInv = currItem['product']['inventories']  # 得到当前目录产品的库存列表for item in listOfInv:sellingPrice = item['price'] # 得到该条目玩具商品价格if sellingPrice > origPrc * 0.5:  # 价格低于原价50% 视为不完整套装print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))retX.append([yr, numPce, newFlag, origPrc]) # 将符合套装信息作为特征存入数据矩阵retY.append(sellingPrice)  # 将对应套装的出售价格存入矩阵except:print('problem with item %d' % i)# 多次调用信息获取函数，获得不同年份不同价格的数据
def setDataCollect(retX, retY):searchForSet(retX, retY, 8288, 2006, 800, 49.99)searchForSet(retX, retY, 10030, 2002, 3096, 269.99)searchForSet(retX, retY, 10179, 2007, 5195, 499.99)searchForSet(retX, retY, 10181, 2007, 3428, 199.99)searchForSet(retX, retY, 10189, 2008, 5922, 299.99)searchForSet(retX, retY, 10196, 2009, 3263, 249.99)if __name__ == '__main__':xArr, yArr = loadDataSet('ex0.txt')# print(xArr[0:2])ws = standRegres(xArr, yArr)  # 计算回归系数# print(ws)xMat = mat(xArr)yMat = mat(yArr)yHat = xMat * ws# fig = plt.figure()# ax = fig.add_subplot(111)# ax.scatter(xMat[:, 1].flatten().A[0], yMat.T[:, 0].flatten().A[0])## xCopy = xMat.copy()# xCopy.sort(0)  #axis=0，按列排列；axis=1，按行排列# yHat = xCopy * ws# ax.plot(xCopy[:, 1], yHat)# plt.show()# yHat = xMat * ws# corrcoef(yHat, yMat)  #corrcoef 用于计算相关系数yHat = lwlrTest(xArr, xArr, yArr, 1)# print(yHat)# 对xArr进行排序srtInd = xMat[:, 1].argsort(0)xSort = xMat[srtInd][:, 0, :]# fig = plt.figure()# ax = fig.add_subplot(111)# ax.plot(xSort[:, 1], yHat[srtInd])   # 横坐标为样本点的第二个特征（第一个特征都是1），纵坐标为预测的y值# ax.scatter(xMat[:, 1].flatten().A[0], mat(yArr).T.flatten().A[0], s=2, c='red')#横坐标相同，纵轴为真实y值# plt.show()'''示例：根据鲍鱼壳的层数预测鲍鱼的年龄'''abX, abY = loadDataSet('abalone.txt')# print(abX)# print("===============")# print(abY)yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)error01 = rssError(abY[0:99], yHat01.T)error1 = rssError(abY[0:99], yHat1.T)error10 = rssError(abY[0:99], yHat10.T)# print(error01) # 56.78868743050092# print(error1) # 429.89056187038# print(error10) # 549.1181708827924# 使用新数据进行测试yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)error01 = rssError(abY[100:199], yHat01.T)# print(error01) #57913.51550155911yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)error1 = rssError(abY[100:199], yHat1.T)# print(error1)  # 573.5261441895982yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)error10 = rssError(abY[100:199], yHat10.T)# print(error10)  # 517.5711905381903# 使用简单线性回归与上面的局部加权线性回归 作比较ws = standRegres(abX[0:99], abY[0:99])yHat = mat(abX[100:199])*wserrorStandReg = rssError(abY[100:199], yHat.T.A)# print(errorStandReg)  # 518.6363153245542'''在鲍鱼数据集上测试岭回归函数'''abX, abY = loadDataSet('abalone.txt')ridgeWeights = ridgeTest(abX, abY)  # 得到30个不同的lambda所对应的回归系数# fig = plt.figure()# ax = fig.add_subplot(111)# ax.plot(ridgeWeights)# plt.show()'''在鲍鱼数据集上测试前向回归函数'''xArr, yArr = loadDataSet('abalone.txt')# stageWise(xArr, yArr, 0.01, 200)# 将上面的结果与最小二乘法进行比较xMat = mat(xArr)yMat = mat(yArr).TxMat = regularize(xMat)yM = mean(yMat, 0)# print(yM)  # [[9.93368446]]yMat = yMat - yMweights = standRegres(xMat, yMat.T)# print(weights.T)'''预测乐高玩具套装的价格'''lgX = []lgY = []setDataCollect(lgX, lgY)