回归---OLS、LWLR、缩减方法、岭回归、lasso、前向逐步回归【含python代码】

回归的目的就是预测数值型的目标值。最直接的办法就是写出一个目标值的计算公式，即所谓的回归方程，需要求方程中的回归系数。一旦有回归系数，就可以进行预测了，具体做法是用回归系数乘以输入值，再将结果全部加起来，就得到预测值了。

下面首先介绍找出最佳拟合直线的两种方法普通最小二乘法（OLS）和局部加权线性回归（LWLR），然后介绍缩减方法，如岭回归、lasso、前向逐步回归。

普通最小二乘法（OLS，Ordinary Least Squares）

核心思想：对于给定的数据集上，找出最佳拟合直线，使得真实值与预测值之间的误差的平方和最小。误差的平方和计算公式如下：
$$\sum _{i=1}^m(y_i-x_i^{T}w{)}^2=(y-Xw{)}^T(y-Xw)$$
如果对w求导，得到$x^T(y-Xw)$，令其等于零，解出w如下：
$$\hat{w}=(X^{T}X{)}^{-1}{X}^{T}y$$
$\hat{w}$表示估计出的w的最优解。因为需要对$(X^{T}X{)}^{-1}$求逆，所以这个方程必须在逆矩阵存在时才适用，但是，矩阵的逆可能也不存在，因此必须要在代码中对此作出判断。

python实现

from numpy import *
import matplotlib.pyplot as plt

def loadDataSet(fileName):
    """
    数据导入
    :param fileName:文件名
    :return: xArr,yArr
    """
    numFeat = len(open(fileName).readline().split('\t')) - 1
    xArr = []
    yArr = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = []
        currLine = line.strip().split('\t')
        for i in range(numFeat):
            lineArr.append(float(currLine[i]))  # 存入的数据需是浮点型的
        xArr.append(lineArr)
        yArr.append(float(currLine[-1])) # 存入的数据需是浮点型的
    return xArr,yArr

def standRegres(xArr, yArr):
    """
    标准回归函数，得到最佳拟合直线对应的回归系数
    :param xArr: 输入数据
    :param yArr: 目标数据
    :return:回归系数
    """
    xMat = mat(xArr)
    yMat = mat(yArr).T  # 转换为列矩阵，即m*1
    xTx = xMat.T * xMat
    # 判断xTx的行列式是否为0，如果是则无法求逆矩阵，如果行列式非0，则可以得到逆矩阵
    if linalg.det(xTx) == 0.0: # 避免后续计算ws时报错
        print("This matrix is singular, cannot do inverse")
        return

    ws = xTx.I*(xMat.T*yMat) # n*1
    # ws = linalg.solve(xTx,xMat.T*yMat) # 与上一句等价
    return ws

def plotFigure(xArr,yArr,yHat):
    """
    绘制出数据集散点图和最佳拟合直线图
    :param xArr: 输入值
    :param yArr: 真实值
    :param yHat: 预测值
    """
    xMat = mat(xArr)
    yMat = mat(yArr)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    # 绘制数据集散点图
    ax.scatter(xMat[:, 1].flatten().A[0], yMat.T[:, 0].flatten().A[0],s=2,c='red')
    # 绘制最佳拟合直线
    # 如果直线上的数据点次序混乱，绘图时将出现问题，首先需要将数据点按xArr升序排列
    srtInd = xMat[:,1].argsort(0) # 获取升序排列对应的索引
    xSort = xMat[srtInd][:,0,:]
    ax.plot(xSort[:,1],yHat[srtInd])
    plt.show()

输入如下代码：

xArr,yArr = loadDataSet('ex0.txt')
print("查看前两条数据：\n",xArr[:2])
ws = standRegres(xArr,yArr)
print("回归系数：\n", ws)
yHat = (mat(xArr) * ws).flatten().A[0] # 预测值
plotFigure(xArr, yArr, yHat)

运行结果如下：

查看前两条数据：
[[1.0, 0.067732], [1.0, 0.42781]]
回归系数：
[[3.00774324]
[1.69532264]]

局部加权线性回归

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

其中一个方法就是局部加权线性回归（Locally Weighted Linear Regression, LWLR）。该算法中，给待预测点附近的每个点赋予一定的权重，然后基于最小均方差来进行普通的回归。该算法解出的回归系数w如下：
$$\hat{w}=(X^{T}WX{)}^{-1}{X}^{T}Wy$$
其中，$W$是一个矩阵，用来给每个数据点赋予权重。

LWLR使用“核”来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：
$$w(i,i)=exp(\frac{|x^{(i)}-x|}{-2k^2})$$
这样就构建了一个只含对角元素的权重矩阵$W$，并且点$x$与$x^{(i)}$越近，$w(i,i)$将会越大。上述公式包含一个需要用户指定的参数k，它决定了对附近的点赋予多大的权重，这也是使用LWLR时唯一需要考虑的参数。

python实现

def lwlr(testPoint,xArr,yArr,k=1.0):
    """
    局部加权线性回归
    :param testPoint:测试点
    :param xArr:输入数据
    :param yArr:数值型目标值
    :param k:控制衰减的速度
    :return:单个预测值
    """
    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*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)) # n*1
    return testPoint * ws  # 单个预测值

def lwlrTest(testArr,xArr,yArr,k=1.0):
    """
    用局部加权线性回归对测试集进行预测
    :param testArr: 测试集输入数据
    :param xArr: 输入数据
    :param yArr: 数值型目标值
    :param k: 控制衰减的速度
    :return: 测试集的预测值
    """
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

输入如下代码：

xArr,yArr = loadDataSet('ex0.txt')
k=1.0  # 改变k值（如1.0 0.01 0.003），观察结果
yHat = lwlrTest(xArr, xArr, yArr, k)
plotFigure(xArr, yArr, yHat)

运行结果如下：

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

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

示例：预测鲍鱼的年龄

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

接下来用LWLR和线性回归两种方法来预测鲍鱼的年龄，首先增加一个计算误差的函数，具体代码如下：

def rssError(yArr,yHatArr):
    """
    计算平方误差
    :param yArr:真实值
    :param yHatArr: 预测值
    :return: 平方误差
    """
    return ((yArr-yHatArr)**2).sum()

输入如下代码：

abX,abY = loadDataSet('abalone.txt')
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.0)
yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
# 用rssError来分析预测误差的大小
print(rssError(abY[0:99],yHat01.T))
print(rssError(abY[0:99], yHat1.T))
print(rssError(abY[0:99], yHat10.T))

运行结果如下：

56.792296224920904
429.8905618702357
549.118170882701

可以看到，使用较小的核得到较低的误差。不能在所有数据集上都使用最小的核，因为会导致过拟合，对新数据不一定能达到最好的预测效果。

下面看下它们在新数据集上的表现，输入如下代码：

yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1.0)
yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
print(rssError(abY[100:199], yHat01.T))
print(rssError(abY[100:199], yHat1.T))
print(rssError(abY[100:199], yHat10.T))

运行结果如下：

39848.954743582595
573.5261441894779
517.5711905381465

可以发现，在上面三个参数中，核大小等于10时的测试误差最小，但它在训练集上的误差却是最大的。

接下来和简单的线性回归做个比较，输入如下代码

ws = standRegres(abX[0:99],abY[0:99])
yHat = mat(abX[100:199])*ws
print(rssError(abY[100:199],yHat.T.A))

运行结果如下：

518.6363153245036

可以看到，简单的线性回归达到了与局部加权线性回归类似的效果。这也表明，必须在未知数据集上比较效果才能选取到最佳的模型。建议用10个不同的样本集做10次测试来比较结果，得到更好的效果。

局部加权线性回归的问题在于，每次必须在整个数据集上运行。也就是说为了做出预测，必须保存所有的训练数据。下面介绍另一种提高预测精度的方法。

缩减系数来“理解”数据

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

为了解决这个问题，统计学家引入岭回归（ridge regression）的概念，这是介绍的第一种缩减方法。接着是lasso法，该方法效果很好但计算复杂。最后介绍第三种缩减方法，称为前向逐步回归，可以得到与lasso差不多的效果，且更容易实现。

注意： 在使用缩减法时，需要对特征作标准化处理，一般对于输入是 $x=\frac{x-xMean}{xVar}$，对输出是 $y-yMean$，使得每维特征具有相同的重要性。

岭回归

岭回归是在矩阵$X^{T}X$上加入一个$\lambda I$从而使得矩阵非奇异，进而能对$X^{T}X+\lambda I$求逆。其中矩阵$I$是一个$n\ast n$的单位矩阵，对角线元素全是1，其他元素全为0.而$\lambda$是一个用户定义的数值。在这种情况下，回归系数的计算公式将变成下式：
$$\hat{w}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}y$$
通过引入$\lambda$来限制所有w之和，通过引入该惩罚项，能够减少不重要的参数，统计学中叫做缩减(shrinkage).

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计
岭回归中的“岭”，体现在单位矩阵I上，对角线元素全为1，形状像一个山岭

缩减技术可以去掉不重要的参数，因此能更好的理解数据。这里通过预测误差最小化得到$\lambda$：数据获取之后，首先抽一部分数据用于测试，剩余的作为训练集用于训练参数w。训练完毕后在测试集上测试预测性能。通过选取不同的$\lambda$来重复上述测试过程，最终得到一个使预测误差最小的$\lambda$。

python实现

def ridgeRegres(xMat, yMat, lam=0.2):
    """
    给定lambda下的岭回归求解,得到回归系数
    :param xMat:输入数据
    :param yMat: 数值型目标值
    :param lam:惩罚因子
    :return:回归系数
    """
    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 regularize(xMat):
    """
    数据标准化，使其分布满足0均值和单位方差
    :param xMat:输入数据
    :return:标准化后的数据
    """
    # 所有特征减去均值并除以方差
    retMat = xMat.copy()
    xMeans = mean(retMat, 0)
    xVar = var(retMat, 0)
    retMat = (retMat - xMeans) / xVar
    return retMat

def ridgeTest(xArr, yArr):
    """
    用于在一组lambda上的测试结果，计算得到不同的回归系数
    :param xArr:输入数据
    :param yArr:数值型目标值
    :return:返回所有的回归系数
    """
    # 为了使用岭回归和缩减技术，首先需要对特征做标准化处理，使每维特征具有相同重要性
    xMat = mat(xArr)
    yMat = mat(yArr).T
    # yMat归一化
    yMean = mean(yMat, 0)
    yMat = yMat - yMean
    # 数据标准化
    xMat = regularize(xMat)
    # 对30个不同的lambda下调用ridgeRegres()函数，并将所有的回归系数输出到一个矩阵中
    numTestPts = 30
    wMat = zeros((numTestPts, shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(xMat, yMat, exp(i - 10))  # lambda是以指数级变化
        wMat[i, :] = ws.T
    return wMat

输入如下代码：

ridgeWeights = ridgeTest(abX, abY)
# 绘制30个不同lambda所对应的回归系数变化图
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()

运行结果如下：

上图绘出了回归系数与$log(\lambda )$的关系。可以看到，$\lambda$非常小时，岭回归回归系数和普通回归一样，而当$\lambda$非常大时，所有回归系数缩减为0，可以在中间某处找到使得预测的结果最好的$\lambda$值。为了定量地找到最佳参数值，还需要进行交叉验证。

注：

• lambda值是以指数级变化，从而可以看出极大和极小值对结果的影响
• 样本数据标准化处理：使每维特征具有相同的重要性

lasso

在增加约束时，普通的最小二乘法回归会得到与岭回归一样的公式：
$$\sum _{k=1}^n{w}_k^2\le \lambda$$
另一个缩减方法lasso也对回归系数做了限定，对应的约束条件如下：
$$\sum _{k=1}^n|w_k|\le \lambda$$
唯一的不同是lasso这个约束条件使用绝对值取代了平方和。当λ足够小时，一些系数会被迫缩减到0，这样可以帮助我们更好的理解数据。

这两个约束条件在公式上看起来相差无几，但是细微的变化却极大地增加了计算的复杂度（为了在这个新的约束条件下解出回归系数，需要使用二次规划算法）。

前向逐步回归

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

迭代多次，每次迭代对每个特征的系数增加或减少一个很小是值，看效果，保存本轮迭代中对因变量变化最大（使得预计值与真实值越接近）的这个特征值的系数。 比如，第$i$个特征值的本轮$w_i$保存，w=[。。。。。，$w_i$，。。。。]，以这个w进行下一轮的迭代。

python实现

def stageWise(xArr, yArr, eps=0.01, numIt=100):
    """
    前向逐步线性回归算法实现
    :param xArr:输入数据
    :param yArr:数值型目标值
    :param eps:每轮迭代需要调整的步长
    :param numIt:迭代次数
    :return:返回迭代后的系数矩阵，每行为一次迭代，最后一行为本地迭代的最优
    """
    # 将数据转换为矩阵，并进行标准化处理
    xMat = mat(xArr)
    yMat = mat(yArr).T
    yMean = mean(yMat, 0)
    yMat = yMat - yMean
    # 特征标准化
    xMat = regularize(xMat)
    m, n = shape(xMat)
    # 每次迭代的权重
    returnMat = zeros((numIt, n))
    # 创建向量ws来保存w的值
    ws = zeros((n, 1))
    wsBest = ws.copy()
    for i in range(numIt):  # 遍历每轮迭代
        print(ws.T)  # 打印w向量，用于分析执行的过程和效果
        # 设置当前最小误差
        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
                    wsBest = wsTest
        ws = wsBest.copy()
        returnMat[i, :] = ws.T
    return returnMat

输入如下代码：

abX, abY = loadDataSet('abalone.txt')
stageWise(abX,abY, 0.01, 200)

运行结果如下：

[[0. 0. 0. 0. 0. 0. 0. 0.]]
[[0. 0. 0. 0.01 0. 0. 0. 0. ]]
[[0. 0. 0. 0.02 0. 0. 0. 0. ]]

…

[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.05 0. 0.09 0.03 0.31 -0.64 0. 0.36]]
[[ 0.04 0. 0.09 0.03 0.31 -0.64 0. 0.36]]

可以看到，每次迭代只改变一个$w_i$的值。从最后一行（本地迭代的最优回归系数）可以发现，w1和w6都是0，说明它们不对目标值造成任何影响，也就是说这些特征很可能是不需要的。

另外，在eps=0.01下，一段时间后，系数已经饱和并在特定值之间来回震荡（如：第一个权重再0.04和0.05之间来回震荡），这是因为步长太大原因。

下面是这用更小的步长和更多的迭代次数，输入如下代码：

stageWise(abX, abY, 0.001, 5000)

运行结果如下：

[[0. 0. 0. 0. 0. 0. 0. 0.]]
[[0. 0. 0. 0.01 0. 0. 0. 0. ]]
[[0. 0. 0. 0.02 0. 0. 0. 0. ]]

…

[[ 0.044 -0.011 0.12 0.022 2.023 -0.963 -0.105 0.187]]
[[ 0.043 -0.011 0.12 0.022 2.023 -0.963 -0.105 0.187]]
[[ 0.044 -0.011 0.12 0.022 2.023 -0.963 -0.105 0.187]]

把这些结果与最小二乘法进行比较，后者的结果可以通过如下代码获得：

xMat = mat(abX)
xMat = regularize(xMat)
yMat = mat(abY).T
yM = mean(yMat,0)
yMat = yMat-yM
weights = standRegres(xMat,yMat.T)
print(weights.T)

运行结果如下：
[[ 0.0430442  -0.02274163  0.13214087  0.02075182  2.22403814 -0.99895312
  -0.11725427  0.16622915]]

可以看到5000次迭代以后，逐步线性回归算法与常规的最小二乘法效果类似。

输入如下代码，可以得到当eps=0.005，且迭代次数为1000时，在鲍鱼数据集上执行逐步线性回归法得到系数与迭代次数间的关系：

stageWiseWeights = stageWise(abX, abY, 0.005, 1000)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(stageWiseWeights)
plt.show()

运行结果如下图所示：

逐步线性回归算法的优点：它能帮助人们理解现有模型并作出改进。构建一个模型后，利用本算法找出重要的特征，及时停止对那些不重要特征的收集。

总结

回归与分类不同点在于，回归预测连续型变量，而分类预测离散型变量。在回归方程里，求得特征对应的最佳回归系数的方法是最小化误差的平方和。常用的求解最佳拟合直线的方法有普通最小二乘法和局部加权线性回归。

岭回归是缩减法的一种，相当于对回归系数的大小施加了限制。另一种很好的缩减法是lasso,lasso很难求解，但是可以使用计算简便的前向逐步线性回归方法求得近似结果。缩减法还可以看做是对一个模型增加偏差的同时减少方差。

参考链接

书籍：《机器学习实战》、周志华的西瓜书《机器学习》