机器学习第八章笔记——预测数值型数据：回归

假设输入数据存放在矩阵x中，而回归系数存放在向量w中，对于给定的数据 $x_{1}$ 预测结果将会通过 $Y_{1} =X_{1}^{T}w$ 给出。对于如何找出w，一个常用的方法是找出使误差最小的w，误差指的是预测y值和真实y值之间的差值，使用该误差的简单累加将使正差值与负差值相互抵消，因此采用平方误差。

平方误差可以写成：

用矩阵表示还可以写作 $(y-Xw)^{T}(y-Xw)$ 。对w求导，得到 $X^{T}(Y-Xw)$ ，令其等于零，解出w如下：

w上方的小标记表示，这是当前可以估计出的w的最优解。上述公式中包含 $X^{T}X^{-1}$ ，因此需要对矩阵求逆。故需要在代码中对是否存在逆矩阵进行判断。

最佳w求解是统计学中的常见问题，可以通过调用Numpy库中的矩阵方法。这种方法被称作OLS，即普通最小二乘法。

对图1中的数据给出最佳拟合直线，代码为：

def loadDataSet(fileName):

    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 standRegres(xArr,yArr):
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T * xMat  # 根据文中推导的公示计算回归系数
    if linalg.det(xTx) == 0.0:  # 判断行列式是否为0，若为0则不存在逆矩阵
        print("矩阵为奇异矩阵,不能求逆")
        return
    ws = xTx.I * (xMat.T*yMat)   # 求出w，.I是对矩阵求逆矩阵
    return ws

第一个函数打开文本文件并提取其中数据，第二个函数用于计算最佳拟合函数，计算方法就为上面讨论中提及的w的计算公式。

接下来是绘制图像：

def plotRegression():
    xArr, yArr = loadDataSet('D:\\learning\\ex0.txt')                                    # 加载数据集
    ws = standRegres(xArr, yArr)                                        # 计算回归系数
    xMat = mat(xArr)                                                    # 创建xMat矩阵
    yMat = mat(yArr)                                                    # 创建yMat矩阵
    xCopy = xMat.copy()                                                    # 深拷贝xMat矩阵
    xCopy.sort(0)                                                        # 排序
    yHat = xCopy * ws                                                     # 计算对应的y值
    fig = plt.figure()
    ax = fig.add_subplot(111)                                            # 添加subplot
    ax.plot(xCopy[:, 1], yHat, c = 'red')                                # 绘制回归曲线
    ax.scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue',alpha = .5)                # 绘制样本点
    plt.title('DataSet')                                                # 绘制title
    plt.xlabel('X')
    plt.show()

调用函数，输出结果就能得到下图

二、局部加权线性回归

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

其中的一个方法是局部加权线性回归。算法中给待预测点附近的每个点赋予一定的权重，基于最小的均方差来进行普通的回归。同KNN一样，算法每次预测均需要事先选取出对应的数据子集，回归系数w的形式为：

w是一个矩阵，用来给每个数据点赋予权重。

LWLR使用“核”来对附近的点赋予更高的权重，核的类型可以自由选择最常用的核就是高斯核，对应权重如下：

这里构建的是一个只含对角元素的权重矩阵w，点x与x(i)越近，w(i,i)将会越大。上述公式包含一个需要用户指定的参数k，它决定了对附近的点赋予多大的权重，也是使用LWLR时唯一需要考虑的参数。

局部加权线性回归函数:

def lwlr(testPoint, xArr, yArr, k = 1.0):
    xMat = mat(xArr); yMat = mat(yArr).T
    m = shape(xMat)[0]
    weights = mat(eye((m)))                   # 创建权重对角矩阵，200行
    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("矩阵为奇异矩阵,不能求逆")
        return
    ws = xTx.I * (xMat.T * (weights * yMat))                            # 计算回归系数
    return testPoint * ws


def lwlrTest(testArr, xArr, yArr, k=1.0):
    m = shape(testArr)[0]                                            # 计算测试数据集大小,200行

    yHat = zeros(m)
    for i in range(m):                                                    # 对每个样本点进行预测
        yHat[i] = lwlr(testArr[i], xArr, yArr, k)  # 第一个参数testArr[i]表示训练集中的每一行，每次拿一个样本进行预测
    return yHat

给定x空间中的任意一点，计算出对应的预测值yHat。函数lwlr的开头用于读入数据并创建所需矩阵之后创建对角权重剧中weights。权重矩阵是一个方阵，阶数等于样本个数。即该矩阵为每个样本点初始化了一个权重，利用算法遍历数据集，计算每个样本点对应的权重值；随样本点与预测值点距离的递增，权重将以指数级衰减。k用于控制衰减的速度。计算完毕可以得到对回归系数w的一个估计lwlrtest则用于为每一个点调用lwlr()，方便求解k的大小。

xArr, yArr = regression.loadDataSet("D:\\learning\\ex0.txt")
result1 = regression.lwlr(xArr[0], xArr, yArr, 1.0)
result2 = regression.lwlr(xArr[0], xArr, yArr, 0.001)
print(result1)
print(result2)

得到的单点估计：

若想得到所有点的估计，可以调用lwlrtest()函数。

绘制出估计值与原始值，以观察拟合效果：

xArr, yArr = regression.loadDataSet("D:\\learning\\ex0.txt")
result1 = regression.lwlr(xArr[0], xArr, yArr, 1.0)
result2 = regression.lwlr(xArr[0], xArr, yArr, 0.001)
yHat = regression.lwlrTest(xArr, xArr, yArr, 0.003)

xMat = mat(xArr)
srtInd = xMat[:, 1].argsort(0)
xSort = xMat[srtInd][:, 0, :]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort[:, 1], yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0], s=2, c='red')
plt.show()

最终结果为：

当然局部加权线性回归也存在一个问题，即增加了计算量。

三、预测鲍鱼的年龄

新添一个rssError()函数用于分析预测误差的大小：

def rssError(yArr, yHatArr):
    return ((yArr-yHatArr)**2).sum()

输入：

abX,abY = regression.loadDataSet('D:\\learning\\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)

result1 = regression.rssError(abY[0:99],yHat01.T)
result2 = regression.rssError(abY[0:99], yHat1.T)
result3 = regression.rssError(abY[0:99], yHat10.T)
print(result1)
print(result2)
print(result3)

运行后：

从上面的数据可以看出：使用较小的核将得到较低的误差，但是在操作中并不能使用最小的核处理数据集，因为使用最小核处理将会造成过拟合，对新数据不一定能达到最好的预测效果。

新数据：

abX,abY = regression.loadDataSet('D:\\learning\\abalone.txt')
yHat01 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
yHat1 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
yHat10 = regression.lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)

result1 = regression.rssError(abY[100:199],yHat01.T)
result2 = regression.rssError(abY[100:199], yHat1.T)
result3 = regression.rssError(abY[100:199], yHat10.T)
print(result1)
print(result2)
print(result3)

从结果上看，核大小等于10时的测试误差最小，不过其在训练集上的误差却是最大的。

四、缩减系数来“理解”数据

当数据的特征比样本点还多时就不能在使用之前的方法进行预测了，否则计算 $(X^{T}X)^{-1}$ 时会出错。

当特征比样本点还多时，输入数据的矩阵x就不是满秩矩阵了，非满秩矩阵在求逆矩阵时会出现错误。为解决这个问题引入岭回归的概念。

4.1 岭回归

岭回归是在矩阵 $(X^{T}X)^{-1}$ 上加上一个 $\lambda x$ 从而使得矩阵非奇异，对 $(X^{T}X)^{-1}+\lambda x$ 求逆，x是一个m*m的单位矩阵，对角线上全是1，其他元素全是0。 $\lambda$ 是随便定义的一个数值。这时回归系数的计算公式就为：

岭回归通过引入 $\lambda$ 来限制了所有w之和，通过引入该惩罚项就能够减少不重要的参数，统计学中将该技术称为缩减。

岭回归使用了单位矩阵乘以常量 $\lambda$ 观察其中单位矩阵I，可以看到1贯穿整个对角线。其余元素均为0.

添加岭回归函数到.py文件中：

def ridgeRegres(xMat, yMat, lam=0.2):
    xTx = xMat.T * xMat
    denom = xTx + eye(shape(xMat)[1]) * lam
    if linalg.dot(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)
    xVar = var(xMat, 0)
    xMat = (xMat - xMeans) / xVar
    numTestPts = 30
    wMat = zeros((numTestPts, shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeTest(xMat, yMat, exp(i- 10))
        wMat[i, :] = ws.T
    return wMat

包含两个函数，ridgeRegres()函数以及ridgeTest()函数，前者用于计算回归系数，后者则用于在一组 $\lambda$ 上测试结果。

ridgeRegres()实现了给定lambda下的岭回归求解。当没指定lambda时，默认为0.2。函数首先构建了矩阵 $X^{T}X$ ，用lam乘以单位矩阵。在普通回归方法可能产生错误时，岭回归仍然可以正常工作。为使用岭回归和缩减技术，需要对特征做标准化处理。具体做法就是将所有特征都减去各自的均值并除以方差。

利用鲍鱼数据进行测试：

abX,abY = regression.loadDataSet('D:\\learning\\abalone.txt')
ridgeWeights = regression.ridgeTest(abX, abY)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()

得到数据：

4.2 lasso

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

上式限定了所有回归系数的平方和不能大于 $\lambda$ 。使用普通的最小二乘法回归在两个或更多的特征相关时，可能会得到一个很大的正系数和一个很大的负系数。与岭回归系数，另一个缩减方法lasso也对回归系数做了限定，对应的约束条件如下：

两个公式的不同点在于这个约束条件使用绝对值取代了平方和。尽管约束形式只是改了一点，但结果却大相径庭。当 $\lambda$ 足够小时，一些系数会因此被迫缩减到0。

4.3 向前逐步回归

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

算法伪代码：

数据标准化，使其分布满足0均值和单位方差

在每轮迭代中：

设置当前最小误差lowestError为正无穷

对每个特征

增大或减小

改变一个系数得到一个新的w

计算新w下的误差

如果误差Error小于当前最小误差lowerstError；设置Wbest等于当前的W

将w设置为新的Wbest

实际代码为：

def stageWise(xArr, yArr, eps=0.01, numIt=100):
    # 标准化处理
    xMat = mat(xArr); yMat = mat(yArr).T
    yMean = mean(yMat,0)
    yMat = yMat - yMean
    xMat = regularize(xMat)  # 归一化xMat
    m,n = shape(xMat)
    returnMat = zeros((numIt,n))
    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()  # 为了实现贪心算法建立两个ws的两份副本
    for i in range(numIt):
        # 每次迭代都打印w向量
        print(ws.T)
        loweastError = 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 < loweastError:
                    loweastError = rssE
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i,:] = ws.T
    return returnMat

# xMat的归一化矩阵
def regularize(xMat):
    inMat = xMat.copy()
    inMeans = mean(inMat, 0)
    inVar = var(inMat, 0)      
    inMat = (inMat - inMeans)/inVar
    return inMat

函数是一个逐步线性回归算法的实现，它与lasso做法相近但计算简单。xArr为输入数据，yArr为预测变量。还有两个参数，一个eps，表示每次迭代需要调整的步长，另一个numIt，表示迭代次数。

数据转换并存入矩阵中，然后把特征按照均值为0方差为1，进行标准化处理。创建向量ws保存w值，并且为了实现贪心算法建立了ws的两份副本。

贪心算法在所有特征上运行两次for循环，分别计算增加或减少该特征对误差的影响。这里使用的是平方误差，通过之前的函数rssError()得到。误差初值设为正无穷经过与所有的误差比较后取最小误差。

xArr, yArr = regression.loadDataSet('D:\\learning\\abalone.txt')
result = regression.stageWise(xArr, yArr, 0.01, 200)
print(result)

输出结果为：

从结果中看到w1和w6都是0，标明它们不对目标造成任何影响，这些特征很可能就是不需要的。

将结果与最小二乘法进行比较，结果如下：