机器学习实战 第八章 预测数值型数据：回归（Regression）

首先介绍线性回归，然后引入局部平滑技术，更好地拟合数据。接着探讨回归在“欠拟合”情况下的缩减（shrinkage），探讨偏差和方差的概念。

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

优点：结果已于理解，计算上并不复杂。

缺点：对非线性的数据拟合不好。

使用数据类型：数值型和标称型。

回归方程（regression equation），回归系数（regression weights），求回归系数的过程就是回归。说到回归，一般都是指线性回归（linear regression），还存在非线性回归模型。

假定输入数据存放在矩阵X中，而回归系数存放在向量w中，那么对于给定的数据X1，预测结果将会通过Y1=XT1w。现在的问题是，手里有一些x和对应的y，如何找到w呢？常用的方法是找到使误差（一般采用平方误差，否则正负误差会相互抵消）最小的w。平方误差可写作：

∑i=1m(yi−xTiw)2

用矩阵表示还可以写成(Y−Xw)T(Y−Xw)，如果对w求导，得到XT(Y−Xw)，令其等于零，解出w如下（这里的w^是w的一个最优解、最佳估计）：

w^=(XTX)−1XTY

公式中的(XTX)−1是对矩阵求逆，需要在代码中判断逆矩阵是否存在，可以调用NumPy中函数

linalg.det()

来计算行列式，如果行列式为零，那么计算逆矩阵的时候会出现错误。这种求解最佳w的方法称为OLS“普通最小二乘法”（ordinary least squares）。

求解完w^，可以计算Y^=Xw^，然后可以通过函数corrcoef(yHat, yActual)来计算预测值和真实值的相关性。

二、局部加权线性回归

线性回归会出现欠拟合现象，因为它求的是最小均方误差的无偏估计。可以在估计中引入一些偏差，从而降低预测的均方误差。其中一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR），该算法中给待测点附近的每个点赋予一定的权重，然后在这个子集上基于最小均方差来进行普通的回归。与kNN一样，此算法每次预测均需事先选取出对应的数据子集。该算法解出的回归系数的形式如下：

w^=(XTWX)−1XTWY

其中w是一个矩阵，用来给每个数据点赋予权重。LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的是高斯核，高斯核对应的权重如下：

w(i,i)=exp(|x(i)−x|−2k2)

这样就构建了一个只含对角元素的权重矩阵W，并且点x与x(i,i)越近，w(i,i)将会越大，公式中的k需要用户指定，它决定了对附近的点赋予多大的权重。

要注意区分权重W和回归系数w，与kNN一样，该加权模型认为样本点距离越近，越可能符合同一个线性模型。



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

上图中，上面的图中k=1.0，权重很大，将所有的数据视为等权重，得出的最佳拟合直线与标准回归一致，属于欠拟合。中间图k=0.01得到了非常好的效果，抓住了数据的潜在模式。下面图使用k=0.003纳入了太多的噪声点，拟合的直线与数据点过于贴近，属于过拟合。

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

三、缩减系数“理解”数据

3-1 岭回归

若数据的特征比样本点还多，在计算(XTX)−1的时候会出错，也就是输入数据的矩阵X不是满秩矩阵，非满秩矩阵在求逆是会出现问题。接下来介绍两种方法来解决这个问题：岭回归（ridge regression）与前向逐步回归（Forward stepwise regression），其中前向逐步回归与lasso法效果差不多。

岭回归就是在矩阵XTX上加上一个λI从而使得矩阵非奇异，进而能对XTX+λI求逆。其中矩阵I是一个m×m的单位矩阵，对角线元素全为1，其他元素全为0。λ是一个自定义的值。这种情况下，回归系数的计算公式将变成：

w^=(XTX+λI)−1XTY

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

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



图3　岭回归的回归系数变化图

上图绘出了回归系统与log(λ)的关系，在最左边λ最小时，得到所有系数的原始值（与线性回归一致）；而在右边全部系数缩减为0；中间部分的某值将可以取得更好的预测效果。为了定量地找到最佳参数值，需要进行交叉验证（cross validation）。

3-2 lasso

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

∑k=1nw2k≤λ

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

与岭回归类似，lasso缩减方法也对回归系数做了限定，其约束条件是：

∑k=1n|wk|≤λ

此约束条件使用绝对值取代了平方和。在λ足够小的时候，一些系数会因此缩减到0。两个约束条件看起来相差无几，但细微变化却极大地增加了计算复杂度，为了在这个新的约束条件下解出回归系数，需要使用二次规划算法。接下来的前向逐步回归算法效果与lasso差不多。

3-3 前向逐步回归

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

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

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

四、权衡偏差与误差

模型和测量值之间存在的差异，叫做误差。当考虑模型中的“噪声”或者说误差时，必须考虑其来源。

对复杂的过程简化，会导致模型和测量值之间出现“噪声”和误差。

无法理解数据的真实生成过程，也会导致差异的发生。

测量过程本身也可能产生“噪声”或问题。

训练误差和测试误差均有三部分组成：偏差、测量误差和随机噪声。局部加权线性回归通过引入越来越小的核来不断增大模型的方差。缩减法可以将一些系数缩减成很小的值或直接缩减为0，增大模型偏差。方差是模型之间的差异，而偏差指的是模型预测值与数据之间的差异。偏差与方差折中的概念在机器学习中十分流行。

程序代码

# coding=utf-8

from numpy import *

# 自适应数据加载函数，该函数能够自动检测出特征的数目，假定最后一个特征是类别标签
def loadDataSet(fileName) :
numFeat = len(open(fileName).readline().split('\t'))
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr = []
curLine = line.strip().split('\t')
for i in range(numFeat - 1) :
lineArr.append(float(curLine[i]))
dataMat.append(lineArr)
labelMat.append(float(curLine[-1]))
return dataMat, labelMat

# 计算最佳拟合直线
def standRegres(xArr, yArr) :
xMat = mat(xArr)
yMat = mat(yArr)
xTx = xMat.T * xMat
# 判断行列式是否为零，为零则计算逆矩阵出现错误
if(linalg.det(xTx) == 0) :
print "This matrix is singular, cannot do inverse."
return
ws = xTx.I * (xMat.T * yMat)
return ws

# 局部加权线性回归函数
def lwlr(testPoint, xArr, yArr, k=1.0) :
xMat = mat(xArr)
yMat = mat(yArr).T
m = 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))
xTx = xMat.T * (weights * xMat)
if linalg.det(xTx) == 0.0 :
print "This matrix is singular, cannot do inverse"
return
ws = xTx.I * (xMat.T * (weights * yMat))
return testPoint * ws

# 为数据集中的每个点调用lwlr()，有助于求解k的大小
def lwlrTest(testArr, xArr, yArr, k=1.0) :
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()

# 给定lambda下的岭回归
def ridgeRegres(xMat, yMat, lam=0.2) :
xTx = xMat.T * xMat
denom = xTx + eye(shape(xMat)[1]) * lam
if linalg.det(denom) == 0.0 :
print "This matrix is singular, cannot do inverse"
return
ws = denom.I * (xMat.T * yMat)
return ws

# 数据的标准化过程，具体过程就是所有特征都减去各自的均值并处理方差
def ridgeTest(xArr, yArr) :
xMat = mat(xArr)
yMat = mat(yArr).T
yMean = mean(yMat, 0)
yMat = yMat - yMean
xMeans = mean(xMat, 0)
# 计算xMat的样本方差
xVar = var(xMat, 0)
xMat = (xMat - xMeans) / xVar
numTestPts = 30
wMat = zeros((numTestPts, shape(xMat)[1]))
for i in range(numTestPts) :
ws = ridgeRegres(xMat, yMat, exp(i-10))
wMat[i,:] = ws.T
return wMat

# 按照均值为0，方差为1，对列进行标准化处理
def regularize(xMat) :
inMat = xMat.copy()
inMeans = mean(inMat, 0)
# 求方差
inVar = var(inMat, 0)
inMat = (inMat - inMeans) / inVar
return inMat

# 前向逐步线性回归
# eps: 每次迭代需要调整的步长
# numIt: 迭代的次数
def stageWise(xArr, yArr, eps=0.01, numIt=100) :
xMat = mat(xArr)
yMat = mat(yArr).T
yMean = mean(yMat, 0)
yMat = yMat - yMean
# 按照均值为0，方差为1进行标准化处理
xMat = regularize(xMat)
m,n = shape(xMat)
returnMat = zeros((numIt, n))
ws = zeros((n, 1))
# 为实现贪心算法建立了ws的两份副本
wsTest = ws.copy()
wsMax = ws.copy()
# 贪心算法在所有特征上运行两次for循环，分别计算增加或减少该特征对误差的影响，这里使用的是平方误差
for i in range(numIt) :
# print ws.T
# 误差初始值设为正无穷
lowestError = 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
return returnMat

在命令行中运行：

>>> import ml.regression as regression
>>> from numpy import *
>>> xArr, yArr = regression.loadDataSet('c:\python27\ml\\ex0.txt')
# 第一个值是1.0，即X0，假定偏移量为一个常数；第二个值是X1，图中横坐标
>>> xArr[0:3]
[[1.0, 0.067732], [1.0, 0.42781], [1.0, 0.995731]]
>>> ws = regression.standRegres(xArr, yArr)
>>> ws
matrix([[ 3.00774324],
[ 1.69532264]])
# 使用ws值计算yHat(为了和y的真实值区分开来)
>>> xMat = mat(xArr)
>>> yMat = mat(yArr)
>>> yHat = xMat*ws
# 绘出数据集散点图
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
<matplotlib.collections.PathCollection object at 0x02E53AF0>
# 绘出最佳拟合直线，如页末图1
>>> xCopy = xMat.copy()
>>> xCopy.sort(0)
>>> yHat = xCopy*ws
>>> ax.plot(xCopy[:,1], yHat)
[<matplotlib.lines.Line2D object at 0x031128D0>]
>>> plt.show()
# 计算预测值和真实值的相关性
>>> yHat = xMat*ws
>>> corrcoef(yHat.T, yMat)
array([[ 1.        ,  0.98647356],
[ 0.98647356,  1.        ]])
# 局部加权线性回归，单点估计
>>> reload(regression)
<module 'ml.regression' from 'C:\Python27\ml\regression.pyc'>
>>> xArr, yArr = regression.loadDataSet('c:\python27\ml\\ex0.txt')
>>> yArr[0]
3.176513
>>> regression.lwlr(xArr[0], xArr, yArr, 1.0)
matrix([[ 3.12204471]])
>>> regression.lwlr(xArr[0], xArr, yArr, 0.001)
matrix([[ 3.20175729]])
# 数据集中，所有点的估计
>>> yHat = regression.lwlrTest(xArr, xArr, yArr, 0.003)
# 对xArr排序
>>> srtInd = xMat[:,1].argsort(0)
>>> xSort = xMat[srtInd][:,0,:]
# 使用Matplotlib绘图
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(xSort[:,1], yHat[srtInd])
[<matplotlib.lines.Line2D object at 0x03106270>]
>>> ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0], s=2, c='red')
<matplotlib.collections.PathCollection object at 0x02E53870>
>>> plt.show()
# 预测鲍鱼年龄
>>> reload(regression)
<module 'ml.regression' from 'C:\Python27\ml\regression.py'>
>>> abX, abY = regression.loadDataSet('c:\python27\ml\\abalone.txt')
# 使用较小的核将得到较低的误差，但易造成过拟合，对新数据预测效果不一定好
>>> yHat01 = regression.lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)
>>> yHat1 = regression.lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
>>> yHat10 = regression.lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
>>> regression.rssError(abY[0:99], yHat01.T)
56.818912090414067
>>> regression.rssError(abY[0:99], yHat1.T)
429.89056187030172
>>> regression.rssError(abY[0:99], yHat10.T)
549.11817088254838
# 使用三种核，在新数据上的预测效果
>>> yHat01 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
>>> regression.rssError(abY[100:199], yHat01.T)
29806.727112900877
>>> yHat1 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
>>> regression.rssError(abY[100:199], yHat1.T)
573.52614418958501
>>> yHat10 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
>>> regression.rssError(abY[100:199], yHat10.T)
517.57119053814722
# 和简单的线性回归作比较
>>> ws = regression.standRegres(abX[0:99], abY[0:99])
>>> yHat = mat(abX[100:199])*ws
>>> regression.rssError(abY[100:199], yHat.T.A)
518.63631532436455
# 鲍鱼数据集上的运行结果
>>> reload(regression)
<module 'ml.regression' from 'C:\Python27\ml\regression.py'>
>>> ridgeWeights = regression.ridgeTest(abX, abY)
# 绘制岭回归的回归系统变化图，见页末图3
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(ridgeWeights)
[<matplotlib.lines.Line2D object at 0x02C92670>, <matplotlib.lines.Line2D object
at 0x02C837F0>, <matplotlib.lines.Line2D object at 0x02C83F10>, <matplotlib.lin
es.Line2D object at 0x02C83E30>, <matplotlib.lines.Line2D object at 0x02E8CC10>,
<matplotlib.lines.Line2D object at 0x02E8CF50>, <matplotlib.lines.Line2D object
at 0x02E8C370>, <matplotlib.lines.Line2D object at 0x02E6C630>]
>>> plt.xlabel('log(lamdbda)')
>>> plt.show()
# 前向逐步回归
>>> reload(regression)
<module 'ml.regression' from 'C:\Python27\ml\regression.pyc'>
>>> xArr, yArr = regression.loadDataSet('c:\python27\ml\\abalone.txt')
>>> regression.stageWise(xArr, yArr, 0.01, 200)
array([[ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
[ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
[ 0.  ,  0.  ,  0.  , ...,  0.  ,  0.  ,  0.  ],
...,
[ 0.05,  0.  ,  0.09, ..., -0.64,  0.  ,  0.36],
[ 0.04,  0.  ,  0.09, ..., -0.64,  0.  ,  0.36],
[ 0.05,  0.  ,  0.09, ..., -0.64,  0.  ,  0.36]])
# 使用更小的步长和更多的步数
>>> regression.stageWise(xArr, yArr, 0.001, 5000)
array([[ 0.   ,  0.   ,  0.   , ...,  0.   ,  0.   ,  0.   ],
[ 0.   ,  0.   ,  0.   , ...,  0.   ,  0.   ,  0.   ],
[ 0.   ,  0.   ,  0.   , ...,  0.   ,  0.   ,  0.   ],
...,
[ 0.043, -0.011,  0.12 , ..., -0.963, -0.105,  0.187],
[ 0.044, -0.011,  0.12 , ..., -0.963, -0.105,  0.187],
[ 0.043, -0.011,  0.12 , ..., -0.963, -0.105,  0.187]])
# 最小二乘法结果，即标准regression方法
>>> xMat = mat(xArr)
>>> yMat = mat(yArr).T
>>> xMat = regression.regularize(xMat)
>>> yM = mean(yMat, 0)
>>> yMat = yMat - yM
>>> weights = regression.standRegres(xMat, yMat.T)
>>> weights.T
matrix([[ 0.0430442 , -0.02274163,  0.13214087,  0.02075182,  2.22403814,
-0.99895312, -0.11725427,  0.16622915]])

图1　ex0.txt数据集与它的最佳拟合直线



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



图3　岭回归的回归系数变化图