常见线性回归|理论与算法实现

01 分类 v.s. 回归

之前我们学习了很多分类方法，在机器学习中，还有一种任务叫回归，回归和分类其实挺像的，都是对样本预测一个值，区别在于，

• 分类：输出为离散值

• 回归：输出为连续值

今天我们学习一波线性回归的理论和算法，不要小看线性回归，其实很多商业模型都少不了线性回归的功劳，把线性回归用到极致你也是大神。

简单来说，线性回归就是在已知x,y的情况下，求解y=wx的回归系数的过程。

02 标准线性回归

标准线性回归就是最简单的线性回归，其原理是最小化预测值与真实值的均方误差，即最小二乘法。根据原理，得到系数公式：

标准线性回归简单粗暴，易于求解和解释，但其缺点也显而易见：对非线性样本，欠拟合。

实现代码如下：

#标准线性回归函数，假设数据有m个样本,n个特征def standRegres(xArr,yArr):
    xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=1*m
    xTx=xMat.T*xMat #shape=n*n
    if np.linalg.det(xTx)==0.0: #矩阵行列式若为0，说明该矩阵不可逆
        print ('xTx为奇异矩阵，不可逆')        return
    ws=xTx.I*(xMat.T*yMat) #shape=n*1,n为特征数
    return ws

测试一波：

dataMat,labelMat=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch08\ex0.txt')
ws=standRegres(dataMat,labelMat)print ('回归系数w为：',round(float(ws[0]),4),round(float(ws[1]),4))print ('回归方程：y={1}x+{0}'.format(round(float(ws[0]),4),round(float(ws[1]),4)))

输出如下，

均方误差为1.355左右，

用于测试的数据集较简单，可以作图看看线性回归的预测效果：

在上图中我们可以看到，预测值非常“直”，虽然捕捉到了真实值的大趋势，但没有捕捉到真实数据的潜在特征(非线性)，下面我们改进一下标准回归算法，尝试捕捉潜在特征。

03 局部加权线性回归

局部加权线性回归认为附近的点有相似的回归特性，因此对待预测点附近的点加权。

之前我写利用拉勾网的数据写了一个人才价格计算器(相关文章在这里)，其原理是将分类算法KNN用于回归，即对样本分类后将同一类别的样本标签设置为该类别的均值，这个人才价格计算器的原理其实和今天这个局部加权线性回归很相似。（当然，缺点也很相似，笑~）

利用这样的想法，局部加权线性回归的回归系数公式如下：

其中，W为权重矩阵，

由权重矩阵公式可知，
1. W是对角矩阵
2. 样本点x距离预测点xi越近，权重越大
3. 对同一个样本点，k越大,权重矩阵W越大,赋予给每个数据点的权重越大,当k减小，W也减小，赋予给各点的权重减小
4. 权重最先减小到0的是距离预测点较远的样本点

下面我们实现此算法

#局部加权线性回归函数，假设数据有m个样本,n个特征#针对某个测试点(testPoint),计算一个ws"""
因为针对每个测试点，都有一个weight赋权矩阵，因此每个测试点都有一个ws回归系数矩阵，
因此对于每个测试点，都需要利用lwlr()得到一个测试点的预测值testPoint*ws
因此lwlr()返回testPoint*ws
"""def lwlr(testPoint,xArr,yArr,k=1.0):
    xMat=np.mat(xArr);yMat=np.mat(yArr) #shape(yMat)=m*1
    m=np.shape(xMat)[0]
    weight=np.mat(np.eye(m)) #构造一个m*m的对角矩阵
    #针对每个样本点j，计算针对测试点i的weight矩阵Wj,j的值
    for j in range(m):
        diffMat=testPoint-xMat[j,:]
        weight[j,j]=np.exp(diffMat*diffMat.T/(-2.0*k**2))
    xTx=xMat.T*weight*xMat    if np.linalg.det(xTx)==0.0: #矩阵行列式若为0，说明该矩阵不可逆
        print ('xTx为奇异矩阵，不可逆')        return
    ws=xTx.I*(xMat.T*weight*yMat)    return testPoint*wsdef lwlrTest(testArr,xArr,yArr,k=1.0):
    m=np.shape(testArr)[0]
    yPred=np.zeros(m)    for i in range(m):
        yPred[i]=lwlr(testArr[i],xArr,yArr,k)    return yPred

测试一下

dataMat,labelMat=loadDataSet(r'D:\DM\python\data\MLiA_SourceCode\machinelearninginaction\Ch08\ex0.txt')
yPred1=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,1.0)
yPred2=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,0.03)
yPred3=lwlrTest(dataMat,dataMat,np.mat(labelMat).T,0.003)
rssError(labelMat,yPred1),rssError(labelMat,yPred2),rssError(labelMat,yPred3)

可以看到，训练集的均方误差随着k的减小而减小，但k过小可能会过拟合(相当于KNN的k取很小的情况，预测点至于其附近几个点相关，就会过拟合)

缺点：类比KNN。对每个点做预测时，都必须用到整个数据集，计算量较大。

04 岭回归ridge

若特征数比样本数还多(n>m)，即特征矩阵X非满秩，那么xTx不可逆，在用上述推导的公式求解w时会出问题，此时怎么办呢？

我们让xTx加上一个值：xTx+lambda*I，通过加上lambda*I使矩阵非奇异，从而求逆矩阵(xTx+lambda*I)，这就是缩减的原理，通过引入lambda来限制回归系数w之和，通过引入lambda惩罚项，减少不重要的参数。

当约束条件如下时，回归就是岭回归：

岭回归的回归系数公式如下：

其中，
1. lambda是一个常数，可以通过k折交叉验证，寻找最小的均方误差，从而找到最佳的lambda
2. I是一个n*n的单位对角矩阵，n为特征数

实现一下，

#求岭回归ws函数def ridgeRegres(xMat,yMat,lam=0.2):
    xTx=xMat.T*xMat
    denom=xTx+lam*np.eye(np.shape(xMat)[1]) #xTx+lambda*I
    if np.linalg.det(denom)==0.0:        print ('xTx为奇异矩阵，不可逆')        return
    ws=denom.I*xMat.T*yMat    return ws#z-score标准化函数def regularize(xMat):
    reMat=xMat.copy()
    reMeans=np.mean(reMat,0) #按列求均值
    reVar=np.var(reMat,0) #按列求标准差
    reMat=(reMat-reMeans)/reVar    return reMat#遍历lambda，求各lambda对应的回归系数ws函数def ridgeTest(xArr,yArr):
    xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=m*1
    yMean=np.mean(yMat,0)
    yMat=yMat-yMean
    xMat=regularize(xMat) #对输入数据标准化
    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
        wMat[i,:]=ws.T #n*1>>1*n
    return wMat

测试，我们用一组鲍鱼的数据来预测鲍鱼年龄，图中8条曲线分别表示w0~w7在30个不同的lambda下的值。

从结果看，回归系数有8个w0-w7，在不同lambda下会得到不同的系数，随着lambda增大，到足够大时(缩减足够大)，各参数减小至接近0，趋于相等，这就是缩减，缩减特征。

05 套索回归lasso+前向逐步回归

同上节，当约束条件如下时，回归就是套索回归：

可以看到，岭回归对平方函数求极值，可使用最小二乘法，但lasso回归对绝对值求极值，求解较复杂，怎么办呢？

这时我们可以使用前向分步算法逼近求解

前向分布算法是一种贪心算法，每一步都对某个权重增加或减少一个很小的值，从而在每一步尽可能减小误差。这样做虽然循环会变多，但每一次循环的计算都变得十分简单。

实现一下

def stageWise(xArr,yArr,step=0.01,numIt=100):
    xMat=np.mat(xArr);yMat=np.mat(yArr).T #shape(yMat)=m*1
    yMean=np.mean(yMat,0)
    yMat=yMat-yMean
    xMat=regularize(xMat) #对输入数据标准化
    m,n=np.shape(xMat)
    returnMat=np.zeros((numIt,n))
    ws=np.zeros((n,1));wsBest=ws.copy    for i in range(numIt): #对每次迭代
        minErr=np.inf #将最小误差初始化为正无穷
        for j in range(n): #对每个特征
            for sign in [-1,1]: #对每个加/减
                wsTest=ws.copy()
                wsTest[j]+=step*sign #对第j个系数进行加/减
                yTest=xMat*wsTest
                rss=rssError(yMat.A,yTest.A)                if rss<minErr:
                    minErr=rss
                    wsBest=wsTest
        ws=wsBest.copy() #找到本次迭代最小误差对应的系数矩阵
        returnMat[i,:]=ws.T    return returnMat

测试，

可以看到，w0,w1,w4,w5,w6这四个系数在每次迭代中都保持在0附近，可以认为这几个系数不重要，可以剔除这几个维度，达到降维的作用(lasso回归的特性之一)

06 总结

本文总结了几种线性回归的原理，并用python3实现了这些回归。

常见的线性回归包括标准线性回归、局部加权线性回归、岭回归Ridge、套索回归Lasso、前向逐步回归(lasso回归的近似求解)。除了这些之外，还有弹性网络回归ElasticNet Regression，它综合了Ridge和Lasso的优点，惩罚项是L1范数和L2范数的融合。

下期我们再来看看树回归，这样回归板块的大部分知识就掌握了，敬请期待~

07 参考

• 《机器学习实战》 Peter Harrington Chapter8

作者：邓莎
链接：https://www.jianshu.com/p/31677806497f