线性回归

客官°小女子只卖身不卖艺 2022-05-31 23:24 333阅读 0赞

准备工作，需要了解矩阵的一些概念，不熟悉的需要自预习：

1矩阵转置

2矩阵内积

3矩阵的代数余子式

在n阶行列式中，把元素aₒₑi所在的第o行和第e列划去后，留下来的n-1阶行列式叫做元素aₒₑi的余子式，记作Mₒₑ，将余子式Mₒₑ再乘以\-1的o+e次幂记为Aₒₑ，Aₒₑ叫做元素aₒₑ的代数余子式。

4矩阵的行列式|A|

[https://zhidao.baidu.com/question/341289987.html?qbl=relate\_question\_0&word=%C7%F3%BE%D8%D5%F3%B5%C4%D0%D0%C1%D0%CA%BD][https_zhidao.baidu.com_question_341289987.html_qbl_relate_question_0_word_C7_F3_BE_D8_D5_F3_B5_C4_D0_D0_C1_D0_CA_BD]

记住二阶和三阶是怎么计算的，高阶的先降阶再处理。

5矩阵奇异性

首先奇异性是方阵特有的属性，行列式|A|是否等于0，若等于0，称矩阵A为奇异矩阵；若不等于0，称矩阵A为非奇异矩阵。只有奇异矩阵才可以求逆。

6伴随矩阵A\*是矩阵元素所对应的代数余子式，所构成的矩阵，转置后得到的新矩阵

7矩阵求逆A⁻¹

A⁻¹=A\*/|A|

8协方差

如果有X,Y两个变量，每个时刻的“X值与其均值之差”乘以“Y值与其均值之差”得到一个乘积，再对这每时刻的乘积求和并求出均值

协方差为正，说明X，Y同向变化，协方差越大说明同向程度越高；如果协方差为负，说明X，Y反向运动，协方差越小说明反向程度越高

![20180114223840707][]

9 X与Y的相关系数

![20180114223904032][]

公式中Cov(X,Y)为X，Y的协方差，D(X)、D(Y)分别为X、Y的方差

线性回归意味着可以将输入项分别乘以一些常量，再将结果加起来得到输出，求这些常量的过程就是线性回归要干的主要工作。假设这个结果是Y=WX，我们要求的就是这个W。预测的结果越贴近真实结果那么回归算法越好，于是变成了求最小平方误差的过程：

![20180114224025622][]

对w求导可得出：

![20180114224043834][]

普通最小二乘法ordinary least squares OLS，就是根据上面的推导简单的套公式求得w。

Python代码如下：

#加载训练集和结果集
    def loadDataSetByfile(fineName) :
        numFeat = len(open(fineName).readline().split('\t'))
        dataMat = []
        labelMat = []
        fr = open(fineName)
        for line in fr.readlines():
            lineArray = []
            curline = line.strip().split('\t')
            for i in range(numFeat-1) :
                lineArray.append(float(curline[i]))
            dataMat.append(lineArray)
            #默认最后一列是结果集
            labelMat.append(float(curline[-1]))
        return dataMat,labelMat
    
    #OLS，直接套公式
    def standRegres(xArr, yArr):
        xMat = mat(xArr)
        yMat = mat(yArr).T
        #求X转置和X的内积
        KX = xMat.T*xMat
        #判断行列式是否为0，如果为0求逆会出错
        if linalg.det(KX)==0 :
            print('该矩阵没有内积')
            return
        w = KX.I * (xMat.T * yMat)
        #下面代码跟上面是等价的
        #w = linalg.solve(KX,xMat.T * yMat)
        return w

代码中在对矩阵内积求逆之前有一个判断，要求行列式不能为零，因为行列式为零的矩阵是没有逆的。

现在需要回归的数据如下，已上传到git上（https://github.com/yejingtao/forblog/blob/master/MachineLearning/trainingSet/regression/ex0.txt），训练样本如下：

1.000000	0.067732	3.176513
    
    1.000000	0.427810	3.816464
    
    1.000000	0.995731	4.550095
    
    1.000000	0.738336	4.256571
    
    1.000000	0.981083	4.560815
    
    1.000000	0.526171	3.929515
    
    1.000000	0.378887	3.526170
    
    1.000000	0.033859	3.156393

利用上面的OLS算法求得ws结果是：\[\[ 3.00774324\]\[ 1.69532264\]\]

利用matplotlib画出我们的回归线型图如下：

dataMat,labelMat = loadDataSetByfile('C:\\2017\\提高\\机器学习\\训练样本\\regression\\ex0.txt')
    w = standRegres(dataMat,labelMat)
    print(w)
    xMat = mat(dataMat)
    yMat = mat(labelMat)
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(111)
    #flatten方法，把矩阵压平为数组
    ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0] )
    xCopy = xMat.copy()
    xCopy.sort(0)
    yPred=xCopy*w
    #绘制回归线
    ax.plot(xCopy[:,1],yPred)
    plt.show()

![20180114224929216][]

如何判断这些模型的好坏呢？计算预测值yHat序列和真实值y序列的匹配程度，那就是计算这两个序列的相关系数。利用corrcoef函数，代码如下：

yPred = xMat*w
    cor = corrcoef(yPred.T,yMat)

YHat和yMat的相关系数结果为0.98，这已经是个很高的值，OLS虽然名字带个简单两字，但是如此看来效果其实也不算差。

线性回归的一个问题是有可能出现欠拟合现象，如果模型欠拟合将不能取得好的预测效果。，我们给待预测点附近的每个点赋予一定的权重，这就是局部加权线性回归（Locally Weighted Linear Regression，LWLR。

回归系数这里用W，或者用D也可以，OLS中的公式就会变成：

![20180114230839309][]

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

这里还是利用高斯核来决定权重大小：

![20180114230917183][]

python代码改造成：

def lwlr(testPoint, xArr,yArr,k=1.0):
        xMat = mat(xArr)
        yMat = mat(yArr).T
        m = shape(xMat)[0]
        #D为对角权重矩阵
        #eye函数：构建对角为1的方阵
        D = mat(eye((m)))
        for j in range(m) :
            diffMat = testPoint - xMat[j,:]
            #距离越远，权重呈指数及衰减
            D[j,j] = exp((diffMat*diffMat.T)/(-2.0*k**2))
        xTx = xMat.T*(D*xMat)
        if linalg.det(xTx)==0 :
            print('该矩阵没有内积')
            return
        #回归公式
        w = xTx.I * (xMat.T*(D*yMat))
        return testPoint*w
    
    #testArr测试集
    #xArr训练集，yArr训练结果集
    #k 高斯系数
    def lwlrTest(testArr, xArr, yArr, k=1.0):
        m = shape(testArr)[0]
        yHat = zeros(m)
        for i in range(m) :
            yHat[i] = lwlr(testArr[i],xArr,yArr,k)
        return yHat

预测结果如下：

K=1时，与OLS一样，相关系数是0.98647703

![20180114231729013][]

K=0.01时，拟合度比较理想，相关系数是0.9985249

![20180114231805710][]

K=0.005时，考虑了太多的噪音，有些过拟合了，相关系数是0.99897028

![20180114231837147][]

上例带了我们对过拟合和欠拟合的思考。

欠拟合比较好理解，就是回归的精度不够，导致无论是训练误差还是测试误差都很大，对于欠拟合的处理需要提高拟合的精度或者用更精度更高的算法继续拟合。

过拟合，过于注重训练集的数据，考虑了太多的噪音和偏差，导致训练误差很小但是测试误差很大，对于过拟合的处理比较麻烦，需要对训练集做某些“钝化”处理，降低训练精度，提高测试精度。

另一个预测鲍鱼寿命的案例引入了更多的特性，来展示过拟合的问题。

案例数据在git上（[https://github.com/yejingtao/forblog/blob/master/MachineLearning/trainingSet/regression/abalone.txt），数据如下，有8][https_github.com_yejingtao_forblog_blob_master_MachineLearning_trainingSet_regression_abalone.txt_8]个特性几千条记录，样本如下：

1	0.455	0.365	0.095	0.514	0.2245	0.101	0.15	15
    1	0.35	0.265	0.09	0.2255	0.0995	0.0485	0.07	7
    -1	0.53	0.42	0.135	0.677	0.2565	0.1415	0.21	9
    1	0.44	0.365	0.125	0.516	0.2155	0.114	0.155	10
    0	0.33	0.255	0.08	0.205	0.0895	0.0395	0.055	7
    0	0.425	0.3	0.095	0.3515	0.141	0.0775	0.12	8
    -1	0.53	0.415	0.15	0.7775	0.237	0.1415	0.33	20
    -1	0.545	0.425	0.125	0.768	0.294	0.1495	0.26	16

为了更精细的计算误差，我们增加了一个函数：

#预测误差
    def rssError(yArr,yHatArr):
        return ((yArr-yHatArr)**2).sum()

我们将庞大的数据集拆分成每100个一份，测试集和训练集不相同，然后用例不同系数的高斯核来验证下过拟合的问题：

X,Y = loadDataSetByfile('C:\\2017\\提高\\机器学习\\训练样本\\regression\\abalone.txt')
    #需要限定X和Y的范围，否则计算量太大了。。。
    yHat1 = lwlrTest(X[0:99],X[0:99],Y[0:99],1)
    yHat10 = lwlrTest(X[0:99],X[0:99],Y[0:99],10)
    yHat01 = lwlrTest(X[0:99],X[0:99],Y[0:99],0.1)
    print('训练集：')
    print(rssError(Y[0:99],yHat01.T))#56.7842091184
    print(rssError(Y[0:99],yHat1.T))#429.89056187
    print(rssError(Y[0:99],yHat10.T))#549.118170883
    print('测试集：')
    yHat1 = lwlrTest(X[100:199],X[0:99],Y[0:99],1)
    yHat10 = lwlrTest(X[100:199],X[0:99],Y[0:99],10)
    yHat01 = lwlrTest(X[100:199],X[0:99],Y[0:99],0.1)
    print(rssError(Y[100:199],yHat01.T))#25119.4591112
    print(rssError(Y[100:199],yHat1.T))#573.52614419
    print(rssError(Y[100:199],yHat10.T))#517.571190538
    print('简单线性回归：')
    w = standRegres(X[0:99],Y[0:99])
    yPred = mat(X[100:199])*w
    print(rssError(Y[100:199],yPred.T.A))#518.636315325

执行结果如下：

训练集：
    56.7842091184
    429.89056187
    549.118170883
    测试集：
    25119.4591112
    573.52614419
    517.571190538
    简单线性回归：
    518.636315325

可以看得出k=0.1时训练误差最小但是测试误差最大，这就是过拟合的问题。同时看得出k=10时与OLS的回归效果差不多，但是要知道局部加权线性回归的计算量确实简单线性回归的N倍，所以至少在这里回归效果相同的话，局部加权法显然是不划算的。

上面算法看似很美好，但是存在一个缺点，（X转置与X的内积）求逆必须是满秩矩阵才有解，如果数据的特征比样本点还多应该怎么办？于是岭回归的概念出线了。

岭回归就是在矩阵XTX上加一个λI从而使得矩阵非奇异，进而能对XTX + λI求逆。其中矩阵I是一个m×m的单位矩阵，对角线上元素全为1，其他元素全为0。而λ是一个用户定义的数值

回归系数的计算公式变成了：

![20180115211204923][]

主要的优点在于它可以帮助人们理解现有的模型并做出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。

再回到rssError这个函数，会发现一个问题，每个特性的量级不一样造成每个特性的误差求和的结果是不平等的，所以我们这次需要先对数据进行标准化处理（(X-均值)/方差）

Python代码：

#岭回归算法
    def ridgeRegres(xMat, yMat, lam=0.2):
        xTx = xMat.T*xMat
        #denom 8*8矩阵
        denom = xTx+eye(shape(xMat)[1])*lam
        if linalg.det(denom)==0.0:
            return
        #岭回归公式
        #m*n矩阵 乘以
        ws = denom.I * (xMat.T*yMat)
        #ws 8*1
        return ws
    
    #所有特征都减去各自的均值并除以方差
    def ridgeTest(xArr, yArr):
        xMat = mat(xArr)
        yMat = mat(yArr).T
        #axis = 0：压缩行，对各列求均值，返回 1* n 矩阵
        #y本来就是个单列，压缩完后是1*1，这里就返回y的平均值
        yMean = mean(yMat,0)
        #yMat现在变成了与平均值的差
        yMat = yMat-yMean
        #压缩完后是1*n
        xMeans = mean(xMat,0)
        #var = mean(abs(x - x.mean())**2)
        xVar = var(xMat,0)
        xMat = (xMat-xMeans)/xVar
        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

λ不同时8个特性对应的线性图：

![20180115211620029][]

最左边时λ为0，相当于没有引入岭回归，相当于普通回归；当λ增大到一定系数后，相当于人为增加的噪音λI远远超过了样本中的有效数据，所有特性的回归系数缩减为0.该图可以帮我们分析出每个特性对线性的影响和变化。

前向逐步回归：

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

代码如下：

#eps步长
    #numIt迭代次数
    def stageWise(xArr, yArr, eps=0.01, numIt=100):
        xMat = mat(xArr)
        yMat = mat(yArr).T
        yMean = mean(yMat,0)
        xMeans = mean(xMat, 0)
        # var = mean(abs(x - x.mean())**2)
        xVar = var(xMat, 0)
        xMat = (xMat - xMeans) / xVar
        #这行代码在原文中是错的
        #xMat=regularize(xMat)
        m,n = shape(xMat)
        returnMat = zeros((numIt,n))
        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()
                    #更新步长后的目前的w
                    wsTest[j] +=eps*sign
                    #预测结果
                    yTest = xMat*wsTest
                    #求预测的误差
                    rssE = rssError(yMat.A,yTest.A)
                    if rssE<lowestError:
                        #保留最小误差的情况，后续做判断
                        lowestError = rssE
                        wsMax = wsTest
            #wsMax是调整过一轮后的最优的w
            ws = wsMax.copy()
            returnMat[i,:] = ws.T
        return returnMat

X,Y = loadDataSetByfile('C:\\2017\\提高\\机器学习\\训练样本\\regression\\abalone.txt')
    returnMat=stageWise(X,Y,0.001,5000)

测试结果：

…….
    [[ 0.043 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]
    [[ 0.044 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]
    [[ 0.043 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]
    [[ 0.044 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]

8个特性都是从0开始，到最后会发现有几个特性来回微幅振荡，这个在前面有介绍过是由于重复的训练集引起的震荡，可以通过随机训练集的方式去掉这个震荡。从结果中也可以看出哪些特性对线性的影响较小，在建模时可以选择取舍，避免不必要的浪费。

偏差Bias、方差Variance与误差Error：

简单说Error = Bias + Variance  
Error反映的是整个模型的准确度，

Bias反映的是模型在样本上的输出与真实值之间的误差，即模型本身的精准度，

Variance反映的是模型每一次输出结果与模型输出期望之间的误差，即模型的稳定性。

举一个例子（来自网友[orangeprince][]），一次打靶实验，目标是为了打到10环，但是实际上只打到了7环，那么这里面的Error就是3。具体分析打到7环的原因，可能有两方面：一是瞄准出了问题，比如实际上射击瞄准的是9环而不是10环；二是枪本身的稳定性有问题，虽然瞄准的是9环，但是只打到了7环。那么在上面一次射击实验中，Bias就是1,反应的是模型期望与真实目标的差距，而在这次试验中，由于Variance所带来的误差就是2，即虽然瞄准的是9环，但由于本身模型缺乏稳定性，造成了实际结果与模型期望之间的差距。

**![20180115212055330][]**

**![20180115212105815][]**

2018-02-28补充：

岭回归之类的缩减法，本质上是增加偏差减少权重的来尝试让拟合度更高。增加的岭虽然在算法上是“噪声”，单只要能提高相似度适当的增加一些干扰也许会有不错的效果。当然了，“噪声的”λ太大，导致偏置太大，就会让原始数据失去了意义。

[https_zhidao.baidu.com_question_341289987.html_qbl_relate_question_0_word_C7_F3_BE_D8_D5_F3_B5_C4_D0_D0_C1_D0_CA_BD]: https://zhidao.baidu.com/question/341289987.html?qbl=relate_question_0&word=%C7%F3%BE%D8%D5%F3%B5%C4%D0%D0%C1%D0%CA%BD
[20180114223840707]: /images/20220601/64e486f4fe0d47ef940771d8e5d52d94.png
[20180114223904032]: /images/20220601/a2b2a7bb4f50440cb0f0f3fb932a269e.png
[20180114224025622]: /images/20220601/8f8dcd86d0634750b4c57dd632ccd4c7.png
[20180114224043834]: /images/20220601/252fa5266e8245ec90e89e18af38d82c.png
[20180114224929216]: /images/20220601/632d5abee3924b75b5d5e8aae997f605.png
[20180114230839309]: /images/20220601/cf4b4bb467644e8cad68f17db963bb3f.png
[20180114230917183]: /images/20220601/da98ae04430e4a33ba1cc4e03fea56a6.png
[20180114231729013]: /images/20220601/a2fa11ee1fbb43c28eb8abd83e1d0692.png
[20180114231805710]: /images/20220601/8e7c635229b04a5998106e053f694f0e.png
[20180114231837147]: /images/20220601/f43d0aec9886412faa90034174b3fd7b.png
[https_github.com_yejingtao_forblog_blob_master_MachineLearning_trainingSet_regression_abalone.txt_8]: https://github.com/yejingtao/forblog/blob/master/MachineLearning/trainingSet/regression/abalone.txt%EF%BC%89%EF%BC%8C%E6%95%B0%E6%8D%AE%E5%A6%82%E4%B8%8B%EF%BC%8C%E6%9C%898
[20180115211204923]: /images/20220601/46bc0e1587894bbeadf2a84aa5c6d1ac.png
[20180115211620029]: /images/20220601/09538aff779c495a87082d911da56a7e.png
[orangeprince]: https://www.zhihu.com/people/orangeprince
[20180115212055330]: /images/20220601/75242256812147b2aed5ed47c37551b7.png
[20180115212105815]: /images/20220601/be67ff9a47884bb495ee832b4b3c59b9.png