《机器学习实战》之八（二）——预测数值型数据：回归

一、前言

  前面我们介绍了线性回归和局部加权线性回归，并配有两示例——预测鲍鱼的年龄。但是在示例中数据集是多维的，有多个数据特征。如果数据的特征比样本点多应该怎么办？是否还可以使用线性回归和之前的方法来做预测呢？答案是否定的，这是因为我无无在计算的时候回出错。
  如果特征比样本点还多（n>m），也就是说输入数据的矩阵X不是忙秩矩阵。非满秩矩阵求逆时会出现问题。
  为了解决这个问题，统计学家引入了岭回归（ridge regression）的概念。这就是本文介绍的缩减系数的一种方法。另外，还有lasso法和前向逐步回归方法。

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

1、岭回归

  岭回归即我们所说的L2正则线性回归，在一般的线性回归最小化均方误差的基础上增加了一个参数w的L2范数的罚项，从而最小化罚项残差平方和：

  简单说来，岭回归就是在普通线性回归的基础上引入单位矩阵 I。回归系数的计算公式将变形如下：

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

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

  本文使用的标准化处理比较简单，具体的做法就是将所有特征都减去各自的均值并除以方差。

  新建文件regression2.py写入代码如下：

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import regression

"""
函数说明:岭回归
Parameters:
    xMat: x数据集
    yMat: y数据集
    lam: 缩减系数
Returns:
     ws: 回归系数
"""
def ridgeRegres(xMat, yMat, lam=0.2):
    xTx = xMat.T * xMat
    denom = xTx + np.eye(np.shape(xMat)[1])*lam
    if np.linalg.det(denom) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T * yMat)
    return ws



"""
函数说明:岭回归测试
Parameters:
    xMat: x数据集
    yMat: y数据集
Returns:
    wMat: 回归系数矩阵
""" 
def ridgeTest(xArr, yArr):
    xMat = np.mat(xArr)
    yMat = np.mat(yArr).T
    
    yMean = np.mean(yMat, axis = 0)  #行方向，求均值
    xMeans = np.mean(xMat, axis = 0)  #行方向，求均值
    xVar = np.var(xMat, axis = 0)  #行方向，求方差
    yMat = yMat - yMean
    xMat = (xMat -xMeans)/xVar  #标准化：数据减去均值除以方差实现标准化
    
    numTestPts = 30  #30个不同的lambda测试
    wMat = np.zeros((numTestPts, np.shape(xMat)[1])) #初始回归系数矩阵
    
    for i in range(numTestPts):
        ws = ridgeRegres(xMat, yMat, np.exp(i-10))  #lambda以e的指数变化，最初是一个非常小的数，

        wMat[i,:] = ws.T   #不同lamda值的回归系数输出到矩阵wMat
    return wMat
    
"""
函数说明：绘制岭回归系数矩阵
Parameters：
    无
Returns：
    无
"""
def plotMat():
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    abX, abY = regression.loadDataSet('abalone.txt')
    redgeweights = ridgeTest(abX, abY)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    
# 生成x轴上的数据:从-10到20，总共有30个点，因为上述λ是e的指数倍
    X = np.linspace(-10, 20, 30)
    ax.plot(X,redgeweights)
    
    ax_title = ax.set_title(u'log(λ)与回归系数的关系',FontProperties = font)
    ax_xlabel = ax.set_xlabel(u'log(λ)', FontProperties = font)
    ax_ylabel = ax.set_ylabel(u'回归系数',FontProperties = font)
    plt.setp(ax_title, size=20, weight='bold', color='red') #设置标题的属性
    plt.setp(ax_xlabel, size=15, weight='bold', color='black') #设置x坐标轴的属性
    plt.setp(ax_ylabel, size=15, weight='bold', color='black') #设置y坐标轴的属性
    plt.show()
    
    
if __name__ == '__main__':
    plotMat()
    

运行结果如下图所示：

  图绘制了回归系数与log(λ)的关系。在最左边，即λ最小时，可以得到所有系数的原始值（结果与线性回归一致）；而在最右边，系数全部缩减成0；在中间部分的某个位置，将会得到最好的预测结果。想要得到最佳的λ参数，可以使用交叉验证的方式获得，后面会继续介绍。

2、lasson

  在增加如下约束时，普通的最小二乘法回归会得到与岭回归一样的公式。
  上式限定了所有的回归系数的平方和不能大于λ。使用普通的最小二乘法回归在当两个或更多的特征时，可能会得出一个很大的正系数和一个很大的负系数。正是因为上述限制条件的存在，使用岭回归可以避免这个问题。
  与岭回归类似，另一个缩减方法lasso也对回归系数做了限定，对于的约束条件如下：

  唯一不同点在于，这个约束条件使用绝对值取代了平方和。虽然约束形式只是稍作编号，结果且大相径庭：在λ足够小的时候，一些系数会因此被迫缩减到0，这个特征可以帮助我们更好地理解数据。这两个约束条件在公式上看起来相差无几，但细微的变化却极大地增加了计算复杂度。下面将介绍一个更为简单的方法来得到结果，该方法叫做前向逐步回归。

3、前向逐步回归

  前向逐步回归算法可以得到与lasso差不多的效果，但更加简单。它属于一种贪心算法，即每一步都尽可能减少误差。 一开始，所有权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。
前向逐步线性回归实现也很简单。当然，还是先进行数据标准化，编写代码如下：

"""
函数说明:数据标准化
Parameters:
    xMat - x数据集
    yMat - y数据集
Returns:
    inxMat - 标准化后的x数据集
    inyMat - 标准化后的y数据集
"""  
def regularize(xMat, yMat):
 
    inxMat = xMat.copy()                                                        #数据拷贝
    inyMat = yMat.copy()
    yMean = np.mean(yMat, 0)                                                    #行与行操作，求均值
    inyMat = yMat - yMean                                                        #数据减去均值
    inMeans = np.mean(inxMat, 0)                                                   #行与行操作，求均值
    inVar = np.var(inxMat, 0)                                                     #行与行操作，求方差
    inxMat = (inxMat - inMeans) / inVar                                            #数据减去均值除以方差实现标准化
    return inxMat, inyMat


"""
函数说明：前向逐步线性回归
Parameters：
    xArr: x输入数据
    yArr: y预测数据
    eps: 每次迭代需要调整的步长
    numIt: 迭代次数
Returns：
    returnMat: numIt次迭代的回归系数矩阵
"""
def stageWise(xArr, yArr, eps=0.01, numIt=100):
    xMat = np.mat(xArr)
    yMat = np.mat(yArr).T
    xMat,yMat= regularize(xMat, yMat) #标准化
    
    m,n = np.shape(xMat) #m个样本，每个样本有n列数据，也称为n个特征
    returnMat = np.zeros((numIt, n)) #初始化最终要返回的回归系数矩阵returnMat
    ws = np.zeros((n,1))
    wsTest = ws.copy()
    wsMax = ws.copy()
    for i in range(numIt): #循环迭代numIt次
        print(ws.T)
        lowestError = np.inf #设置当前最小的误差为正无穷
        for j in range(n): #对于每个特征
            for sign in [-1,1]: #对于增大或缩小
                wsTest = ws.copy()
                wsTest[j] += eps*sign #改变一个系数得到一个新的系数w
                yTest = xMat*wsTest   
                rssE = regression.rssError(yMat.A, yTest.A)  #计算新w下的误差
                if rssE < lowestError:  #选择最小误差
                    lowestError = rssE
                    wsMax = wsTest
            ws = wsMax.copy()
            returnMat[i,:] = ws.T   #存储新的ws
    return returnMat

"""
函数说明:绘制前向逐步线性回归系数矩阵
"""
def plotstageWiseMat():
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    xArr, yArr = regression.loadDataSet('abalone.txt')
    returnMat = stageWise(xArr, yArr, 0.005, 1000)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(returnMat)
    
    ax_title = ax.set_title(u'前向逐步回归：迭代次数与回归系数的关系',FontProperties = font)
    ax_xlabel = ax.set_xlabel(u'迭代次数', FontProperties = font)
    ax_ylabel = ax.set_ylabel(u'回归系数',FontProperties = font)
    plt.setp(ax_title, size=20, weight='bold', color='red') #设置标题的属性
    plt.setp(ax_xlabel, size=15, weight='bold', color='black') #设置x坐标轴的属性
    plt.setp(ax_ylabel, size=15, weight='bold', color='black') #设置y坐标轴的属性
    plt.show()

运行结果如下：

迭代次数与回归系数的关系曲线，可以看出，有些系数从始至终都是约为0的，这说明它们不对目标造成任何影响，也就是说这些特征很可能是不需要的。逐步线性回归算法的优点在于它可以帮助人们理解有的模型并做出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。

ps：注意这里书上面没有给出regularize函数，它是对输入数据集进行标准化的函数

总结
  

缩减方法（逐步线性回归或岭回归），就是将一些系数缩减成很小的值或者直接缩减为0。这样做，就增大了模型的偏差（减少了一些特征的权重），通过把一些特征的回归系数缩减到0，同时也就减少了模型的复杂度。消除了多余的特征之后，模型更容易理解，同时也降低了预测误差。但是当缩减过于严厉的时候，就会出现过拟合的现象，即用训练集预测结果很好，用测试集预测就糟糕很多。

三、权衡偏差与方差

这里需要再查资料。。。。。。

四、示例：预测乐高玩具套装的价格

  乐高（LEGO）公司生产拼装类玩具，由很多大小不同的塑料插块组成。一般来说，这些插块都是成套出售，它们可以拼装成很多不同的东西，如船、城堡、一些著名建筑等。乐高公司每个套装包含的部件数目从10件到5000件不等。

  一种乐高套件基本上在几年后就会停产，但乐高的收藏者之间仍会在停产后彼此交易。本次实例，就是使用回归方法对收藏者之间的交易价格进行预测。

1、获取数据

  书中使用的方法是通过Google提供的API进行获取数据，但是现在这个API已经关闭，我们无法通过api获取数据了。不过幸运的是，我在网上找到了书上用到的那些html文件。

  我们通过解析html文件，来获取我们需要的信息，如果学过python网络爬虫，那么这部分的内容会非常简单，新建legoData.py编写代码如下：

# -*- coding: utf-8 -*-

import numpy as np
from bs4 import BeautifulSoup
import regression
import regression2

"""
函数说明:从页面读取数据，生成retX和retY列表
Parameters:
    retX :数据X
    retY:数据Y
    inFile:HTML文件
    yr:年份
    numPce:乐高部件数目
    origPrc: 原价
Returns:
    无
"""
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
    # 打开并读取HTML文件
    with open(inFile, encoding='utf-8') as f:
        html = f.read()
    soup = BeautifulSoup(html)
    i = 1
    # 根据HTML页面结构进行解析
    currentRow = soup.find_all('table', r = "%d" % i)
    while(len(currentRow) != 0):
        currentRow = soup.find_all('table', r = "%d" % i)
        title = currentRow[0].find_all('a')[1].text
        lwrTitle = title.lower()
        # 查找是否有全新标签
        if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
            newFlag = 1.0
        else:
            newFlag = 0.0
        # 查找是否已经标志出售，我们只收集已出售的数据
        soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
        if len(soldUnicde) == 0:
            print("商品 #%d 没有出售" % i)
        else:
            # 解析页面获取当前价格
            soldPrice = currentRow[0].find_all('td')[4]
            priceStr = soldPrice.text
            priceStr = priceStr.replace('$','')
            priceStr = priceStr.replace(',','')
            if len(soldPrice) > 1:
                priceStr = priceStr.replace('Free shipping', '')
            sellingPrice = float(priceStr)
            # 去掉不完整的套装价格
            if  sellingPrice > origPrc * 0.5:
                print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
                retX.append([yr, numPce, newFlag, origPrc])
                retY.append(sellingPrice)
        i += 1
        currentRow = soup.find_all('table', r = "%d" % i)

"""
函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
Parameters:
    无
Returns:
    无
"""
def setDataCollect(retX, retY):
    scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)     #2006年的乐高8288,部件数目800,原价49.99
    scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
    scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
    scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)  #2007年的乐高10181,部件数目3428,原价199.99
    scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)  #2008年的乐高10189,部件数目5922,原价299.99
    scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)   #2009年的乐高10196,部件数目3263,原价249.99

if __name__ == '__main__':
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)

运行结果如下：

我们对没有的商品做了处理。这些特征分别为：出品年份、部件数目、是否为全新、原价、售价（二手交易）

2、建立模型

  处理好数据集后，接下来就是建立模型。首先我们需要添加全为1的特征X₀列。因为线性回归的第一列特征要求都是1.0。然后使用线性回归，在legoData.py文件中编写代码如下：

"""
函数说明:使用简单的线性回归
Parameters:
    无
Returns:
    无
"""
def useStandRegres():
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)
    data_num, features_num = np.shape(lgX)
    lgX1 = np.mat(np.ones((data_num, features_num+1)))
    lgX1[:,1:5] = np.mat(lgX) #第0列仍为1，第[1,4)列设置为lgX
    ws = regression.standRegres(lgX1, lgY)
    print("%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原件" % (ws[0],ws[1],ws[3],ws[3],ws[4]))

if __name__ == '__main__':
    useStandRegres()

运行结果如下：

  从运行结果来看，这个模型对于数据拟合得很好，但是看上没有什么道理。套件里的部件数量越多，售价反而降低了，这是不合理的。

  接下来，我们使用岭回归，通过交叉验证，找到使误差最小的λ对应的回归系数。在legoData.py文件中继续编写代码如下：

"""
函数说明:交叉验证岭回归
Parameters:
    xArr: x数据集
    yArr: y数据集
    numVal: 交叉验证次数
Returns:
    wMat: 回归系数矩阵
"""
def crossValidation(xArr, yArr, numVal=10):
    m = len(yArr)                       #统计样本个数
    indexList = list(range(m))          #生成索引值列表
    errorMat = np.zeros((numVal, 30))   #create error mat 30columns numVal rows
    for i in range(numVal):             #交叉验证numVal次
        trainX = []; trainY = []        #训练集
        testX = []; testY = []          #测试集
        np.random.shuffle(indexList)    #混洗，打乱次序
        
        #划分数据集:90%训练集，10%测试集
        for j in range(m):            
            if j < m * 0.9:
                trainX.append(xArr[indexList[j]])
                trainY.append(yArr[indexList[j]])
            else:
                testX.append(xArr[indexList[j]])
                testY.append(yArr[indexList[j]])
        wMat = regression2.ridgeTest(trainX, trainY)  #获得30个不同lambda下的岭回归系数
       
        #遍历30个不同lambda下的岭回归系数
        for k in range(30):
            matTestX = np.mat(testX); matTrainX = np.mat(trainX)
            meanTrain = np.mean(matTrainX, axis=0)
            varTrain = np.var(matTrainX, axis=0)
            matTestX = (matTestX-meanTrain)/varTrain     #用训练时的参数将测试数据标准化
            yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY)     #根据ws预测y值
            errorMat[i,k] = regression.rssError(yEst.T.A, np.array(testY))  #存储误差
    meanErrors = np.mean(errorMat, axis=0)   #计算每次交叉验证的平均误差
    minMean = float(min(meanErrors))        #找到最小误差
    bestWeights = wMat[np.nonzero(meanErrors == minMean)]   #找到最佳回归系数
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    meanX = np.mean(xMat, axis=0)
    varX = np.var(xMat, axis=0)
    unReg = bestWeights/varX           #数据经过标准化，因此需要还原
    print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % ((-1 * np.sum(np.multiply(meanX,unReg)) + np.mean(yMat)), unReg[0,0], unReg[0,1], unReg[0,2], unReg[0,3]))
    
if __name__ == '__main__':
    #useStandRegres()
    
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)
    crossValidation(lgX, lgY)

运行结果如下：

  这里随机选取样本，因为其随机性，所以每次运行的结果可能略有不同。不过整体如上图所示，可以看出，它与常规的最小二乘法，即普通的线性回归没有太大差异。我们本期望找到一个更易于理解的模型，显然没有达到预期效果。

  现在，我们看一下在缩减过程中回归系数是如何变化的。在legoData.py文件中运行如下代码：

if __name__ == '__main__':
    #useStandRegres()
    
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)
    #crossValidation(lgX, lgY)
    print(regression2.ridgeTest(lgX, lgY))

运行结果如下：

  看运行结果的第一行，可以看到最大的是第4项，第二大的是第2项。

  因此，如果只选择一个特征来做预测的话，我们应该选择第4个特征，也就是原始价格。如果可以选择2个特征的话，应该选择第4个和第2个特征。

这种分析方法使得我们可以挖掘大量数据的内在规律。在仅有4个特征时，该方法的效果也许并不明显；但如果有100个以上的特征，该方法就会变得十分有效：它可以指出哪个特征是关键的，而哪些特征是不重要的。

五、使用Sklearn的linear_model

  最后让我们用sklearn实现下岭回归吧。

  官方英文文档地址：点击查看

  sklearn.linear_model提供了很多线性模型，包括岭回归、贝叶斯回归、Lasso等。本文主要讲解岭回归Ridge，要看其他模型的可以看官方文档（https://scikit-learn.org/stable/modules/linear_model.html#ridge-regression）。

参数说明如下：

• alpha：正则化系数，float类型，默认为1.0。正则化改善了问题的条件并减少了估计的方差。较大的值指定较强的正则化。
• fit_intercept：是否需要截距，bool类型，默认为True。也就是是否求解b。
• normalize：是否先进行归一化，bool类型，默认为False。如果为真，则回归X将在回归之前被归一化。 当fit_intercept设置为False时，将忽略此参数。 当回归量归一化时，注意到这使得超参数学习更加鲁棒，并且几乎不依赖于样本的数量。 相同的属性对标准化数据无效。然而，如果你想标准化，请在调用normalize = False训练估计器之前，使用preprocessing.StandardScaler处理数据。
• copy_X：是否复制X数组，bool类型，默认为True，如果为True，将复制X数组; 否则，它覆盖原数组X。
• max_iter：最大的迭代次数，int类型，默认为None，最大的迭代次数，对于sparse_cg和lsqr而言，默认次数取决于scipy.sparse.linalg，对于sag而言，则默认为1000次。
• tol：精度，float类型，默认为0.001。就是解的精度。
• solver：求解方法，str类型，默认为auto。可选参数为：auto、svd、cholesky、lsqr、sparse_cg、sag。
  – auto根据数据类型自动选择求解器。 svd使用X的奇异值分解来计算Ridge系数。对于奇异矩阵比cholesky更稳定。
  – cholesky使用标准的scipy.linalg.solve函数来获得闭合形式的解。
  – sparse_cg使用在scipy.sparse.linalg.cg中找到的共轭梯度求解器。作为迭代算法，这个求解器比大规模数据（设置tol和max_iter的可能性）的cholesky更合适。
  – lsqr使用专用的正则化最小二乘常数scipy.sparse.linalg.lsqr。它是最快的，但可能在旧的scipy版本不可用。它是使用迭代过程。
  – sag使用随机平均梯度下降。它也使用迭代过程，并且当n_samples和n_feature都很大时，通常比其他求解器更快。注意，sag快速收敛仅在具有近似相同尺度的特征上被保证。您可以使用sklearn.preprocessing的缩放器预处理数据。
• random_state：sag的伪随机种子。
  以上就是所有的初始化参数，当然，初始化后还可以通过set_params方法重新进行设定。

现在就可以编写代码了，在legoData.py文件中写入如下代码：

"""
函数说明:使用sklearn库中的线性模型拟合数据
Parameters:
    无
Returns:
    无
"""
def usesklearn():
    from sklearn import linear_model
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)
    reg = linear_model.Ridge(alpha = .5)
    reg.fit(lgX, lgY)
    print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3])) 
    
    
if __name__ == '__main__':
    #useStandRegres()
    
#    lgX = []
#    lgY = []
#    setDataCollect(lgX, lgY)
#    #crossValidation(lgX, lgY)
#    print(regression2.ridgeTest(lgX, lgY))

    usesklearn()

运行结果如下图所示：

  我们不搞太复杂，正则化项系数设为0.5，其余参数使用默认即可。可以看到，获得的结果与上小结的结果类似。

总结

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

参考资料

• matplotlib坐标轴设置方法：https://www.jb51.net/article/129823.htm
• ax.plot(xMat)中Mat类型参数是m*n结构的图像：
  https://blog.csdn.net/u014540876/article/details/79331341
• 《机器学习实战》第八章 预测数值型数据：回归
• 博客：https://blog.csdn.net/c406495762/article/details/82967529
• sklearn官方文档：https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html