机器学习实战第8章预测数值型数据:回归2

Posted weiququ

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了机器学习实战第8章预测数值型数据:回归2相关的知识,希望对你有一定的参考价值。

1. Shrinkage(缩减) Methods

当特征比样本点还多时(n>m),输入的数据矩阵X不是满秩矩阵,在求解(XTX)-1时会出现错误。接下来主要介绍岭回归(ridge regression)和前向逐步回归(Foward Stagewise Regression)两种方法。

1.1 岭回归(ridge regression)

简单来说,岭回归就是在矩阵XTX上加上一个技术分享图片从而使得矩阵非奇异,进而能进行求逆。其中矩阵I是一个单位矩阵,技术分享图片是一个调节参数。

岭回归的回归系数计算公式为:

技术分享图片

岭回归最先用来处理特征数多于样本的情况,现在也用于在估计中加入偏差,来得到更好的估计。当技术分享图片=0时,岭回归和普通线性回归一样,当技术分享图片趋近于正无穷时,回归系数接近于0。可以通过交叉验证的方式得到使测试结果最好的技术分享图片值。

岭回归代码:

 1 # encoding: utf-8
 2 ‘‘‘
 3 @Author: shuhan Wei
 4 @File: ridgeRegression.py
 5 @Time: 2018/8/24 11:16
 6 ‘‘‘
 7 import numpy as np
 8 import matplotlib.pyplot as plt
 9 
10 def loadDataSet(fileName):
11     numFeat = len(open(fileName).readline().split(	)) - 1
12     dataMat = []; labelMat = []
13     fr = open(fileName)
14     for line in fr.readlines():
15         lineArr = []
16         curLine = line.strip().split(	)
17         for i in range(numFeat):
18             lineArr.append(float(curLine[i]))
19         dataMat.append(lineArr)
20         labelMat.append(float(curLine[-1]))
21     return dataMat, labelMat
22 
23 
24 """
25     函数说明:岭回归函数 w=(XT*X+lamI)-1 * XT*y 
26     Parameters:
27         xMat - x数据集
28         yMat - y数据集
29         lam - 用户自定义参数
30     Returns:
31         ws  - 回归系数
32 """
33 def ridgeRegres(xMat, yMat, lam=0.2):
34     xTx = xMat.T * xMat
35     denom = xTx + np.eye(np.shape(xMat)[1]) * lam
36     if np.linalg.det(denom) == 0.0:
37         print("The matrix is singular, cannot do inverse")
38         return
39     ws = denom.I * (xMat.T * yMat)
40     return ws
41 
42 
43 """
44     函数说明:测试岭回归
45     Parameter:
46         xArr - x数据集矩阵
47         yArr - y数据集矩阵
48     Returns:
49         wMat - 回归系数矩阵,30组系数
50 """
51 def ridgeTest(xArr, yArr):
52     xMat = np.mat(xArr); yMat = np.mat(yArr).T
53     #对数据进行标准化处理
54     yMean = np.mean(yMat, 0)
55     yVar = np.var(yMat,0)
56     yMat = yMat - yMean
57     xMeans = np.mean(xMat,0)    #行操作
58     xVar = np.var(xMat, 0)
59     xMat = (xMat - xMeans) / xVar
60     numTestPts = 30 # 在30个不同lambda值下
61     wMat = np.zeros((numTestPts, np.shape(xMat)[1]))  #系数矩阵
62     for i in range(numTestPts):
63         ws = ridgeRegres(xMat, yMat, np.exp(i-10))
64         wMat[i,:] = ws.T
65     return wMat
66 
67 """
68     函数说明:绘制回归系数与lambda关系曲线
69 """
70 def plotLambda(ridgeWeights):
71     plt.rcParams[font.sans-serif] = [SimHei]
72     plt.rcParams[axes.unicode_minus] = False
73     fig = plt.figure()
74     ax = fig.add_subplot(111)
75     ax.plot(ridgeWeights)
76     ax_title_text = ax.set_title(ulog(lambada)与回归系数的关系)
77     ax_xlabel_text = ax.set_xlabel(ulog(lambada))
78     ax_ylabel_text = ax.set_ylabel(u回归系数)
79     plt.setp(ax_title_text, size=20, weight=bold, color=red)
80     plt.setp(ax_xlabel_text, size=10, weight=bold, color=black)
81     plt.setp(ax_ylabel_text, size=10, weight=bold, color=black)
82     plt.show()
83 
84 
85 if __name__ == __main__:
86     abX, abY = loadDataSet(abalone.txt)
87     ridgeWeights = ridgeTest(abX, abY)
88     plotLambda(ridgeWeights)

技术分享图片

 1.2 前向逐步回归(Foward Stagewise Regression)

前向逐步回归算法属于一种贪心算法,即每一步都尽可能减少误差,一开始所有的权重都设置为1,然后每一步对权重增加或减少一个很小的值。

算法是伪代码如下:

 数据标准化,使其分布满足0均值和单位方差
  在每轮迭代过程中:
    设置当前最小误差lowestError为正无穷
    对每个特征:
      增大或缩小:
        改变一个系数得到一个新的W
        计算新W下的误差
        如果误差Error小于lowestError,设置wbest = 当前W
      将W设置为新的wbest

代码如下:

 1 # encoding: utf-8
 2 ‘‘‘
 3 @Author: shuhan Wei
 4 @File: stageWise.py
 5 @Time: 2018/8/24 13:10
 6 ‘‘‘
 7 import numpy as np
 8 import regression
 9 import matplotlib.pyplot as plt
10 
11 
12 def regularize(xMat):#regularize by columns
13     inMat = xMat.copy()
14     inMeans = np.mean(inMat,0)   #calc mean then subtract it off
15     inVar = np.var(inMat,0)      #calc variance of Xi then divide by it
16     inMat = (inMat - inMeans)/inVar
17     return inMat
18 
19 
20 def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays
21     return ((yArr-yHatArr)**2).sum()
22 
23 
24 """
25     函数说明:向前逐步回归算法
26     Parameters:
27         xArr - x数据集
28         yArr - y数据集
29         eps - 每次迭代步长
30         numIt - 迭代次数
31     Returns:
32         系数矩阵
33 """
34 def stageWise(xArr, yArr, eps=0.01, numIt = 100):
35     xMat = np.mat(xArr)
36     yMat = np.mat(yArr).T
37     yMean = np.mean(yMat, 0)
38     yMat = yMat - yMean
39     xMat = regularize(xMat)
40     m, n = np.shape(xMat)
41     returnMat = np.zeros((numIt, n))  # testing code remove
42     ws = np.zeros((n, 1));wsTest = ws.copy();wsMax = ws.copy()
43     for i in range(numIt):  # could change this to while loop
44         print(ws.T)
45         lowestError = np.inf;   #初始化最小误差为正无穷
46         for j in range(n):  #遍历每个特征
47             for sign in [-1, 1]:
48                 wsTest = ws.copy()
49                 wsTest[j] += eps * sign
50                 yTest = xMat * wsTest
51                 rssE = rssError(yMat.A, yTest.A)
52                 if rssE < lowestError:
53                     lowestError = rssE
54                     wsMax = wsTest
55         ws = wsMax.copy()
56         returnMat[i, :] = ws.T
57     return returnMat
58 
59 
60 def plot(weights):
61     plt.rcParams[font.sans-serif] = [SimHei]
62     plt.rcParams[axes.unicode_minus] = False
63     fig = plt.figure()
64     ax = fig.add_subplot(111)
65     ax.plot(weights)
66     ax_title_text = ax.set_title(u回归系数与迭代次数的关系)
67     ax_xlabel_text = ax.set_xlabel(u迭代次数)
68     ax_ylabel_text = ax.set_ylabel(u回归系数)
69     plt.setp(ax_title_text, size=20, weight=bold, color=red)
70     plt.setp(ax_xlabel_text, size=10, weight=bold, color=black)
71     plt.setp(ax_ylabel_text, size=10, weight=bold, color=black)
72     plt.show()
73 
74 
75 if __name__ == __main__:
76     abX, abY = regression.loadDataSet(abalone.txt)
77     weights = stageWise(abX, abY, 0.001, 5000)
78     plot(weights)

结果:

        技术分享图片

        中间省略

        技术分享图片

技术分享图片

当迭代次数很大时,接近于普通最小二乘回归

小结:当使用缩减方法时,把一些系数的回归系数缩减到0,减少了模型的复杂度,增加了模型的偏差,同时减小了方差

权衡偏差和方差

      技术分享图片

 

1.3 实例:预测乐高玩具价格

在交叉验证中使用90%数据集作为训练数据,10%做为测试数据,通过10次迭代之后,求30个不同lambda下的误差均值,再选择最小的误差,从而获取最优回归系数。

使用岭回归时需要对x数据集进行标准化,对y数据集进行中心化。最后需要对计算的回归系数进行还原,同时截距等于技术分享图片

# encoding: utf-8
‘‘‘
@Author: shuhan Wei
@File: legao.py
@Time: 2018/8/24 16:16
‘‘‘
#BeautifulSoup是一个可以从html或XML文件中提取数据的Python库
import random

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

def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays
    return ((yArr-yHatArr)**2).sum()

def scrapePage(inFile, outFile, yr, numPce, origPrc):
    # 打开并读取HTML文件
    with open(inFile, encoding=utf-8) as fr:
        html = fr.read()
    fw = open(outFile, a)  # a是追加模式
    soup = BeautifulSoup(html)
    i = 1
    currentRow = soup.findAll(table, r="%d" % i)
    while (len(currentRow) != 0):
        currentRow = soup.findAll(table, r="%d" % i)
        title = currentRow[0].findAll(a)[1].text
        lwrTitle = title.lower()
        print(lwrTitle)
        if (lwrTitle.find(new) > -1) or (lwrTitle.find(nisb) > -1):
            newFlag = 1.0   #全新
        else:
            newFlag = 0.0   #不是全新
        soldUnicde = currentRow[0].findAll(td)[3].findAll(span)
        if len(soldUnicde) == 0:
            print
            "item #%d did not sell" % i
        else:
            soldPrice = currentRow[0].findAll(td)[4]
            priceStr = soldPrice.text
            priceStr = priceStr.replace($, ‘‘)  # strips out $
            priceStr = priceStr.replace(,, ‘‘)  # strips out ,
            if len(soldPrice) > 1:
                priceStr = priceStr.replace(Free shipping, ‘‘)  # strips out Free Shipping
            print
            "%s	%d	%s" % (priceStr, newFlag, title)
            fw.write("%d	%d	%d	%f	%s
" % (yr, numPce, newFlag, origPrc, priceStr))
        i += 1
        currentRow = soup.findAll(table, r="%d" % i)
    fw.close()


def setDataCollect():
    scrapePage(setHtml/lego8288.html, out.txt, 2006, 800, 49.99)    #2006年出品,部件数目800,原价49.99
    scrapePage(setHtml/lego10030.html, out.txt, 2002, 3096, 269.99)
    scrapePage(setHtml/lego10179.html, out.txt, 2007, 5195, 499.99)
    scrapePage(setHtml/lego10181.html, out.txt, 2007, 3428, 199.99)
    scrapePage(setHtml/lego10189.html, out.txt, 2008, 5922, 299.99)
    scrapePage(setHtml/lego10196.html, out.txt, 2009, 3263, 249.99)


"""
    函数说明:交叉验证测试岭回归
    Parameters:
        xArr - x数据集
        yArr - y数据集
        numVal - 交叉验证的次数
"""
def crossValidation(xArr,yArr,numVal=10):
    m = len(yArr)
    indexList = list(range(m))  #生成索引列表
    errorMat = np.zeros((numVal,30))#存放误差值
    for i in range(numVal):
        trainX=[]; trainY=[]
        testX = []; testY = []
        random.shuffle(indexList)    #打乱indexList,为了下面随机选取数据
        for j in range(m):  #使用90%的数据做训练集,10%做测试集
            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 = ridgeRegression.ridgeTest(trainX,trainY)    #使用岭回归获取30组回归系数
        print(wMat)
        for k in range(30):#loop over all of the ridge estimates
            matTestX = np.mat(testX); matTrainX = np.mat(trainX)
            meanTrain = np.mean(matTrainX,0)
            varTrain = np.var(matTrainX,0)
            matTestX = (matTestX-meanTrain)/varTrain #正规化测试数据
            yEst = matTestX * np.mat(wMat[k,:]).T + np.mean(trainY)#计算预测y值
            errorMat[i,k] = rssError(yEst.T.A,np.array(testY))
            #print errorMat[i,k]
    meanErrors = np.mean(errorMat,0)#calc avg performance of the different ridge weight vectors
    minMean = float(min(meanErrors))    #获取最小的误差值
    bestWeights = wMat[np.nonzero(meanErrors==minMean)]
    #can unregularize to get model
    #when we regularized we wrote Xreg = (x-meanX)/var(x)
    #we can now write in terms of x not Xreg:  x*w/var(x) - meanX/var(x) +meanY
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    meanX = np.mean(xMat,0); varX = np.var(xMat,0)
    unReg = bestWeights / varX
    print("the best model from Ridge Regression is:
",unReg)
    print("with constant term: ",-1*np.sum(np.multiply(meanX,unReg),1) + np.mean(yMat))  #常数项


if __name__ == __main__:
    # setDataCollect()
    lgX, lgY = regression.loadDataSet(out.txt)
    #加上X0
    lgXMat = np.mat(np.ones((np.shape(lgX)[0],np.shape(lgX)[1]+1)))
    lgXMat[:,1:5] = np.mat(lgX)
    print(lgXMat[0])
    ws = regression.standRegres(lgXMat,lgY)
    print(ws)
    print(lgXMat[0]*ws,lgY[0])
    crossValidation(lgX,lgY)

  使用最小二乘回归预测结果为:

   技术分享图片

  它对数据拟合得很好,但看上去却没什么道理,从公式看零部件数量越多,反而销售价更低。

  使用岭回归:

  技术分享图片

  技术分享图片

  该结果和最小二乘结果没有太大差异,我们本期望找到一个更易于理解的模型,显然没有达到预期效果。为了达到这一点,我们看一下在缩减回归中系数是如何变化的

  技术分享图片

  看运行结果的第一行,可以看到最大的是第4项,第二大的是第2项。因此,如果只选择一个特征来做预测的话,我们应该选择第4个特征,也就是原始加个。如果可以选择2个特征的话,应该选择第4个和第2个特征。

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

使用sklearn中的岭回归

def usesklearn():
    from sklearn import linear_model
    reg = linear_model.Ridge(alpha=.5)
    lgX, lgY = regression.loadDataSet(out.txt)
    setDataCollect(lgX, lgY)
    reg.fit(lgX, lgY)
    print("$%f%f*Year%f*NumPieces+%f*NewOrUsed+%f*original price" % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3],reg.coef[4]))

 

 

            

 

 

 

以上是关于机器学习实战第8章预测数值型数据:回归2的主要内容,如果未能解决你的问题,请参考以下文章

机器学习实战第8章预测数值型数据:回归

机器学习实战笔记-预测数值型数据:回归

预测数值型数据-回归(Regression)

机器学习——预测数值型数据:回归

机器学习实战精读--------回归

机器学习实战 - ApacheCN