《机器学习实战笔记--第二部分利用回归预测数值型数据：回归1》_岭回归得出的最佳拟合直线(size = 0.9 + 0.8 x weight)怎么得出的-CSDN博客

本文链接：https://blog.csdn.net/qq_41635352/article/details/80736743

分类的目标变量是标称型数据，而本章将会对连续型的数据做出预测。

1. 利用线性回归找到最佳拟合直线

回归的目的是预测数值型的目标值。最直接的办法是依据输入写出一个目标值的计算公式。这个公式就是回归方程，其中的系数是回归系数，求回归系数的过程就是回归。

说到回归一般是指线性回归，所以本章里的回归与线性回归是一个意思。线性回意味着可以将输入项分别乘以一些常量，再将结果加起来得到输出。

输入数据存放在矩阵X中，回归系数存放在向量w中，对于给定的数据X1，预测结果会通过Y1 = X1.T * w给出。现在就是要找到w。一般我们找出使误差最下的w。这里误差是指预测值y和真实值t之间的差值，常用平方误差。

平方差可以写作：。

用矩阵形式还可以写作（y - Xw）.T * (y - Xw)，如果对w求导，得到X.T(Y-Xw)，令其等于0，解出w如下：

w上方的小标记表示，这是当前可以估计出的w的最优解。从现有的数据集上估计出的w可能不是数据中真是的w值，所以在这里加一个‘帽’来表示它仅仅是w的一个最佳估计。

值得注意的是，上述公式中有X的逆，因此这个方程仅仅只存在逆矩阵存在的时候使用。必须要在代码中对此作出判断。

上述求w的求解是统计学中最常见的问题。除了矩阵还有很多方法可以解决。通过调用numpy库里的库矩阵方法，我们可以仅使用几行代码就可以完成功能，该方法也称为‘普通最小二乘法’。

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

from numpy import *

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,labelMat

def standRegres(xArr,yArr):
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T*xMat
    if linalg.det(xTx) == 0.0:
        print ("This matrix is singular, cannot do inverse")
        return
    ws = xTx.I * (xMat.T*yMat)
    return ws

第一个函数就是打开用tab分隔的文本文件，文件每行的最后的一个值是目标值。

第二个函数用来计算最佳拟合直线。该函数首先保存x、y，然后计算X.T*X，然后判断行列式是否为0.如果行列式为0，那么计算就会出现错误。如果行列式非0，就计算返回w。

Numpy提供一个线性代数库linalg，可直接调用linalg.det()来计算行列式；还提供了一个函数来解未知矩阵，如果使用该函数，那么代码ws=(x.T*x).I * (xMat.T * yMat)应写成 ws = linalg.solve(xTx, xMat.T*yMatT)。

运行结果：

第一个值总是1.0，即X0。我们假定偏移量就是一个常数。第二个值X1，也就是我们图中的横坐标。

变量ws存放的就是回归系数。在用内积来预测y的时候，第一维将乘以前面的常数X0，第二维将乘以X1。因为前面假定了X0=1，所以最终会得到 y=ws[0]+ws[1]*X1。这里的y是预测值，为了和真实值y分开，我们将其记为yHat。

画出数据散点图和最佳拟合直线：

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
# flatten返回一个折叠成一维的数组。但是该函数只能适用于numpy对象，即array或者mat，普通的list列表是不行的。
ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
#plt.show()

xCopy = xMat.copy()
#print(xCopy)
xCopy.sort(0)
#print (xCopy)
yHat = xCopy*ws
ax.plot(xCopy[:,1], yHat)
plt.show()

几乎任何一个数据集都可以用上述方法建立模型，那么如何判断模型的好坏呢？比较下图的两个子图，如果在两个数据集上分别做线性回归，将得到一样的模型。

显然两个数据集是不一样的，那么模型分别在二者上的效果如何呢？应该如何比较二者的效果呢？有两种方法来计算预测值yhat和真实序列y之间的匹配程度，那就是计算两个序列的相关系数。

在python中Numpy库提供了相关函数的计算方法，可以通过命令corrcoef(yEstimate, yActual)来计算预测值和真实值的相关性。下面我们就在前面的数据集上做个试验。

计算相关系数（yHat转置，保证两个都是行向量）

yHat = xMat*ws
corrcoef(yHat.T, yMat)

可以看到对角线上都是1，因为yMat， yHat和自己匹配最完美，yHat和yMat的相关系数是0.98.

最佳拟合直线方法将数据视为直线进行建模。

2. 局部加权线性回归

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

其中的一个方法是 局部加权线性回归（LWLR）。在该算法中，我们给待预测点附近的每个点赋予一定的权重；然后与第一节类似，在这个子集上基于最小均方差来进行普通的回归。与KNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解出回归系数w的形式如下：

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

LWLR使用‘核’（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

这样就构建了一个只含对角元素的权重矩阵w，并且点x与x(i)越近，w(i,i)将会越大。k是用户指定参数，它决定了对附近的点赋予多大的权重，这也是使用LWLR时唯一要考虑的参数，在下图中可以看到参数k和权重的关系。

下面来看看模型效果：

局部加权线性回归函数：

# 给定x空间中的任意一点，计算对应的预测值yHat
def lwlr(testPoint,xArr,yArr,k=1.0):
    xMat = mat(xArr); yMat = mat(yArr).T
    m = shape(xMat)[0]
    # numpy.eye()生成对角权重矩阵  weights仅含有对角元素其他都为0，对应矩阵乘法，且w(i,i)值对应于Xi的权重
    '''
    第一个参数：输出方阵（行数=列数）的规模，即行数或列数
    第三个参数：默认情况下输出的是对角线全“1”，其余全“0”的方阵。
    如果k为正整数，则在右上方第k条对角线全“1”其余全“0”，k为负整数则在左下方第k条对角线全“1”其余全“0”
    '''
    weights = mat(eye((m))) # 只含对角元素的对角权重矩阵，方阵，阶数等于样本点数。为每个样本点初始化一个权重
    # 遍历数据集，计算每个样本点对应的权重值
    for j in range(m):                      # next 2 lines create weights matrix
        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

def lwlrTest(testArr,xArr,yArr,k=1.0):  #loops over all the data points and applies lwlr to each one
    m = shape(testArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

观察运行结果：

不同k对应不同的权重（每个点），预测出的结果也不同

    yHat = lwlrTest(xArry, xArry, yArry, 0.003)

xMat = mat(xArry)
srtInd = xMat[:,1].argsort(0)
xSort = xMat[srtInd][:,0,:]

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

对训练数据集中的所有样本都进行局部加权线性回归预测并画图观察：

k=0.003

k=1.0

k=0.01

使用3种不同平滑值绘制出的局部加权线性回归结果。上图中的平滑参数k分别取不同值。

可以看到k=1.0时与最小二乘法效果差不多，k=0.01时可以挖掘出数据的潜在规律，而k=0.003时则考虑了太多噪声，进而导致了过拟合。

局部加权回归也存在一个问题，即增加了计算量，因为它对每个点做也测时必须使用整个数据集。从上图中也可以看出，k=0.01时可以得到很好的估计，但同时看一下k=0.01的情况，发现大部分数据点的权重都接近0。如果可以避免这些计算可以减少程序的计算时间，从而缓解因计算量增加而带来的问题。

3. 示例：预测鲍鱼的年龄

接下来我们将应用于真实数据集。

为了分析预测误差的大小，可以使用函数rssError（）计算这一指标：

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

选择不同的核的运行结果：

可以看到使用较小的核将得到较低的误差。但是使用较小的核将会产生过拟合现象，对线数据的预测不一定能达到最好的效果。下面我们就来看一下它们在新数据集上的表现

上面三个参数中，核大小等于10时，测试误差最小，但是他在训练集上的误差确实最大的。接下来和简单的线性回归做一个比较：

可以看到简单的线性回归达到了与局部加权回归类似的效果，这也表明了，在未知的数据集上比较才能选取最好的模型。那么最佳的核大小是10吗？或许是吧，但是若想得到更好的效果，应该使用10个不同的样本集做十次测试来比较。

下面介绍另一种提高精度的方法，并分析它的优势所在。

4. 缩减系数来‘理解’数据

如果数据的特征比样本点还多怎么办呢？是否还可以使用之前的方法来做预测？答案是否定的，即不能使用之前的方法。这是因为在计算（X.T*X）.I时会出错。

特征比样本点多说明输入矩阵X不是满秩矩阵，非满秩矩阵在求逆时会出现问题。

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

4.1 岭回归

简单来说，岭回归就是在矩阵X.T*X上加上一个

从而使得矩阵非奇异，进而能对

求逆。其中I是m*m的单位阵，对角元素全为1，其他元素为0。而

是用户定义值。这个时候，回归系数的计算公式：

岭回归最先使用在处理特征多于样本数目的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过加入

来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术也称为缩减。

缩减方法可以去掉不重要的参数，因此能更好的理解数据，此外与简单的线性回归相比，缩减法能取得更好的预测效果。

与之前的训练其他参数所使用的方法类似，这里通过预测误差最小化得到：数据获取之后，首先抽一部分数据用于测试，剩余的数据来训练参数w。训练完毕之后再测试集上测试性能。通过选取不同的来重复上述的过程，最终得到一个使预测误差最小的。

岭回归:

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
    '''
     mean()函数功能：求取均值  （var()类似）
        经常操作的参数为axis，以m * n矩阵举例：
        axis 不设置值，对 m*n 个数求均值，返回一个实数
        axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
        axis = 1：压缩列，对各行求均值，返回 m *1 矩阵
    '''
    yMean = mean(yMat,0)
    yMat = yMat - yMean     #to eliminate X0 take mean off of Y
    #regularize X's
    xMeans = mean(xMat,0)   #calc mean then subtract it off
    xVar = var(xMat,0)      #calc variance of Xi then divide by it
    xMat = (xMat - xMeans)/xVar
    # 在30个不同的lam下调用ridgeRegres（）函数 
    # lam以指数级变化可以看到取值在非常小和非常大时分别对结果造成的影响
    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

运行一下程序：

得到30个不同的lam的值对应的回归系数，该图绘出了回归系数与的关系。在最左边，最小时，可以所有的系数的原始值（与线性回归一致）；右边系数全部缩减为0；在中间部分的某些值可以取得最好的预测效果。为了定量的找到最佳参数，还需要进行交叉验证。另外要判断哪些变量对预测结果最有影响力，在图中观察系数即可。

还有一些其他的缩减方法，如lasso，lar，pca回归以及子集选择。与岭回归一样，这样的方法不仅可以提高预测准确率，还可以解释回归系数。

4.2 lasso

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

上式限定了所有回归系数的平方和不能大于。与岭回归类似，lasso也对回归系数做了限定，对应约束条件如下：

两者不同点在于约束的形式变化了，结果大相径庭。在足够小的时候，一些系数会被缩减到0，这个特性可以帮助我们理解数据。

4.3 前向逐步线性回归

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

前向逐步线性回归

def stageWise(xArr,yArr,eps=0.01,numIt=100):
    '''
    输入数据，预测变量，迭代需要调整的步长，迭代次数
    '''
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0)
    yMat = yMat - yMean     #can also regularize ys but will get smaller coef
    xMat = regularize(xMat)
    m,n=shape(xMat)
    returnMat = zeros((numIt,n)) #testing code remove
    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()
    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

运行结果：