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

一、前言

  前面我们介绍了线性回归和局部加权线性回归，并配有两示例——预测鲍鱼的年龄。但是在示例中数据集是多维的，有多个数据特征。如果数据的特征比样本点多应该怎么办？是否还可以使用线性回归和之前的方法来做预测呢？答案是否定的，这是因为我无无在计算的时候回出错。
  如果特征比样本点还多（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()
    
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

运行结果如下图所示：

  图绘制了回归系数与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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

运行结果如下：

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

总结
  

三、权衡偏差与方差

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

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

  乐高（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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

运行结果如下：

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

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

运行结果如下：

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

  接下来，我们使用岭回归，通过交叉验证，找到使误差最小的λ对应的回归系数。在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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

运行结果如下：

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

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

if __name__ == '__main__':
    #useStandRegres()
    
    lgX = []
    lgY = []
    setDataCollect(lgX, lgY)
    #crossValidation(lgX, lgY)
    print(regression2.ridgeTest(lgX, lgY))
1
2
3
4
5
6
7
8

运行结果如下：

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

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

五、使用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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

运行结果如下图所示：

  我们不搞太复杂，正则化项系数设为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