虽然在线性回归的求解过程中，通过借助最⼩⼆乘法能够迅速找到全域最优解，但最⼩⼆乘法本身的使⽤条件较为苛刻，必须要求当 为满秩矩阵或正定阵时才可进⾏逆矩阵或⼴义逆矩阵的求解， 在实际应⽤中经常会遇⻅矩阵不存在逆矩阵或⼴义逆矩阵的情况，并且当X的各列存在线性相关关系 （即多重共线性）的时候，最⼩⼆乘法的求解结果不为⼀，这⾥需要注意的是，在进⾏数据采集过程 中，数据集各列是对同⼀个客观事物进⾏描述，很难避免多重共线性的存在，因此存在共线性是很多 数据集的⼀般情况。当然更为极端的情况则是数据集的列⽐⾏多，此时最⼩⼆乘法⽆法对其进⾏求 解。 总的来说，解决共线性的问题的⽅法主要有以下几种： 其⼀是在建模之前对各特征进⾏相关性检验，若存在多重共线性，则可考虑进⼀步对数据集进⾏ SVD分解或PCA主成分分析，在SVD或PCA执⾏的过程中会对数据集进⾏正交变换，最终所得数 据集各列将不存在任何相关性。当然此举会对数据集的结构进⾏改变，且各列特征变得不可解 释。 其⼆则是采⽤逐步回归的⽅法，以此选取对因变量解释⼒度最强的⾃变量，同时对于存在相关性 的⾃变量加上⼀个惩罚因⼦，削弱其对因变量的解释⼒度，当然该⽅法不能完全避免多重共线性 的存在，但能够绕过最⼩⼆乘法对公线性较为敏感的缺陷，构建线性回归模型； 其三则是在原有的算法基础上进⾏修改，放弃对线性⽅程参数⽆偏估计的苛刻条件，使其能够容 忍特征列存在多重共线性的情况，并且能够顺利建模，且尽可能的保证RSS取得最⼩值。 通常来说，能够利⽤⼀个算法解决的问题尽量不⽤多个算的组合来解决，因此此处我们主要考虑后两 个解决⽅案，特别第三种方案我们用到岭回归算法和Lasso算法。 本文来自CDA数据分析学院（吴昊天）老师
所谓岭回归就是在原有⽅程系数计算公式中添加了⼀个扰动项aI,原先⽆法求⼴义逆的情况变成可以求出其⼴义逆，使得问题稳定并得以求解，
其中 a是⾃定义参数， I则是单位矩阵
In [ ]:
w=(x.T * x +aI).I * X.T * y
In [3]:
#对于对单位矩阵的⽣成，需要借助NumPy中的eye函数，该函数需要输⼊对⻆矩阵规模参数 import numpy as np import pandas as pd import matplotlib as mpl import matplotlib.pyplot as plt np.eye(3)
Out[3]:
In [4]:
def ridgeRegres(dataSet, lam=0.2): xMat = np.mat(dataSet.iloc[:, :-1].values) yMat = np.mat(dataSet.iloc[:, -1].values).T xTx = xMat.T * xMat denom = xTx + np.eye(xMat.shape[1])*lam ws = denom.I * (xMat.T * yMat) return ws
In [6]:
aba = pd.read_table('abalone.txt', header = None)#该数据集源于UCI，记录了鲍⻥的⽣物属性，⽬标字段是该⽣物的年龄 aba.head()
Out[6]:
In [7]:
#其中最后⼀列为标签列，同时，由于数据集第⼀列是分类变量，且没有⽤于承载截距系数的列，因此 #此处将第⼀列的值全部修改为1，然后再进⾏建模 aba.iloc[:, 0] = 1 aba.head()
Out[7]:
In [8]:
rws = ridgeRegres(aba) rws #返回各列的系数同样的第一行为截距
Out[8]:
In [9]:
#封装R**2 def rSquare(dataSet, regres):#设置参数为 数据集 与 回归方法 sse = sseCal(dataSet, regres) y = dataSet.iloc[:, -1].values sst = np.power(y - y.mean(), 2).sum() return 1 - sse / sst
In [10]:
#封装最小二乘 def standRegres(dataSet): xMat = np.mat(dataSet.iloc[:, :-1].values)#将DataFrame转换成array 再将array转换成matrix(矩阵) 因为只有矩阵可以进行计算 yMat = np.mat(dataSet.iloc[:, -1].values).T xTx = xMat.T*xMat if np.linalg.det(xTx) == 0:#判断xTx是否是满秩矩阵，若不满秩，则⽆法对其进⾏求逆矩阵的操作 print("This matrix is singular, cannot do inverse") return ws = xTx.I * (xMat.T*yMat) return ws #这⾥需要注意的是，当使⽤矩阵分解来求解多元线性回归时，必须添加⼀列全为1的列，⽤于表征线性⽅程截距b。
In [12]:
#将SSE做一个封装 def sseCal(dataSet, regres):#设置参数为 数据集 与 回归方法 n = dataSet.shape[0] y = dataSet.iloc[:, -1].values ws = regres(dataSet) yhat = dataSet.iloc[:, :-1].values * ws yhat = yhat.reshape([n,]) rss = np.power(yhat - y, 2).sum() return rss
In [15]:
rSquare(aba, ridgeRegres)
Out[15]:
In [14]:
rSquare(aba, standRegres)
Out[14]:
调⽤Scikit—Learn中岭回归算法验证⼿动建模有效性
In [16]:
from sklearn import linear_model ridge = linear_model.Ridge(alpha=.2) ridge.fit(aba.iloc[:, :-1], aba.iloc[:, -1])
Out[16]:
In [17]:
ridge.coef_#返回各行系数
Out[17]:
In [18]:
ridge.intercept_#返回截距
Out[18]:
绘制岭迹图
绘制岭迹图基本函数基本思路为从⼩到⼤依次取 a，然后查看各特征列系数的衰减速度，从中查看是否存在公线性，以及 的最佳取值
In [19]:
def ridgeTest(dataSet): xMat = np.mat(dataSet.iloc[:, :-1].values) yMat = np.mat(dataSet.iloc[:, -1].values).T yMean = np.mean(yMat, axis = 0) yMat = yMat - yMean xMeans = np.mean(xMat, axis = 0) xVar = np.var(xMat,axis = 0) xMat = (xMat - xMeans)/xVar numTestPts = 30 wMat = np.zeros((numTestPts,xMat.shape[1])) for i in range(numTestPts): ws = ridgeRegres(dataSet,np.exp(i-10)) wMat[i,:]=ws.T return wMat
In [20]:
aba = pd.read_table('abalone.txt', header = None) ridgeWeights = ridgeTest(aba) plt.plot(ridgeWeights) plt.xlabel('log(lambda)') plt.ylabel('weights')
Out[20]:
把所有回归系数的岭迹都绘制在一张图上，如果这些曲线比较稳定，如上图所示，利用最小二乘估计会有一定的把握。
在Scikit-Learn中执⾏Lasso算法
In [22]:
from sklearn import linear_model las = linear_model.Lasso(alpha = 0.01) las.fit(aba.iloc[:, :-1], aba.iloc[:, -1])
Out[22]:
In [23]:
ridge.coef_#返回各行系数
Out[23]:
In [24]:
ridge.intercept_#返回截距
Out[24]: