Python线性回归的8种方法

Tirthajyoti Sarkar在Github上写了一篇文章，名叫Linear regression with various methods，概述了8种线性回归的Python实现方法和速度。文章很经典，在此对其进行翻译和注解。原文地址如下：

https://github.com/tirthajyoti/PythonMachineLearning/blob/master/LinearRegressionMethods.ipynb

0 说明

文章介绍了8种常用的线性回归方法，其原理是完全相同的，即最小二乘估计（OLS），因此结果也没有任何差异。开发人员会在意代码效率（即完成回归的时间），但我们往往不太在意；也就是说，其实这8种方法对我们来说是一样的。

与原文的写法不同，为了更加清晰地展示方法差异，我们在使用到某个库时才会载入这个库，不同的方法中会重复载入。

计算回归效率（时间）的方法很简单，即在回归的开始和结束处使用time库的time函数获取当前时间，并对两次获得的时间作差。为节约篇幅，我们省略这些代码，它们与回归本身毫无关系。

数据集

文章的数据集是自己造的，思路如下：

生成一个范围在-10到10之间的等差数列作为x，数列元素个数为$5*10^6$（用5e6表示）；

使用numpy库下的 函数（scipy库下也有），生成确定斜率和截距的y，这时y是一条直线；

使用 函数为y添加扰动，人为制造随机数据。

代码如下：

有必要解释一下 的用法。 是规律生成多项式的函数，其语法是

其中，参数x是多项式变量取值的序列，比如这里是一个从-10到10的等差数列；参数p是多项式变量前的系数列表，比如这里是[3.25, -6.5]。这两个参数确定后，我们就拥有了多项式3.25*x-6.5。不难想象，如果参数p是[1, 2, 3, 4]，多项式就成了x³+2x²+3x+4。这很酷对不对？

这时x和y形成了一条斜率为3.25、截距为-6.5的直线，准确地说是线段。我们等间距地选取50个点展示一下：

好了，准备工作就到这里。

1 方法一：摩尔-彭若斯伪逆

摩尔-彭若斯伪逆（Moore–Penrose Pseudoinverse）是对奇异矩阵（不满秩，行列式为0）或非方阵（行数不等于列数）求广义逆矩阵的方法。如果矩阵X是矩阵A的伪逆矩阵，则满足AXA=A和XAX=X。不难看出伪逆是兼容逆的，如果矩阵X和矩阵A都是非奇异矩阵，那么消去等式两边的相同元素，我们有AX=XA=I。

现在，系数θ与横坐标向量X的点乘等于纵坐标向量Y，即θ·X=Y，可得

那么求解思路很清晰了，先求出变量值矩阵的伪逆，再与y点乘。求伪逆用函数，点乘用方法。代码如下：

注意方程的常数项被视为x的零次方（即1）的系数，因此 的参数是2维的：[[x[0],1], [x[1],1], ..., [x[int(5e6)-1],1]]。我们可以使用列表推导式，但将[[i,1] for i in x]替换为纯Numpy方法np.vstack([x, np.ones(len(x))]).T后，运算用时将从4.921秒降到0.066秒。

2 方法二：投影矩阵

这个方法从空间几何的角度来理解最小二乘法。我们可以简单地思考如下：

假设这是一元线性回归，根据OLS的定义，如果变量x的系数为b，常数项为c，我们希望找到满足最小的b和c。

而在空间中，表示n-1维空间上任意一个点，表示n维空间的一个向量。我们假想这就发生在三维的空间里，那么x向量（v1）和常数向量（v0）确定了一个二维平面，我们希望找到y向量与该平面欧式距离（平方和）最小的点所确定的b和c。那很简单，我们只需要找到y向量在v0-v1平面上的投影p。

我们利用向量y-p（其实就是垂线）垂直于平面v0-v1就必然垂直于平面上任意向量（包含v0和v1）的性质，得到

令，我们有

可以解出系数C

这个结论同样很酷，就好像用二向箔估计线性回归系数。V就是方法一中的X，y就是因变量y，而必然是非奇异矩阵，我们可以用 函数求逆。代码如下：

方法三：statsmodels.api.OLS

这个方法来自statsmodels库，在之前的文章中有过介绍。模块的引入约定如下：

用statsmodels库估计参数要分成两步：

先使用 函数构建回归模型；

再使用fit方法获取估计结果。

有两点需要特别注意：

statsmodels.api提供了专门的添加常数项的函数 ，不用我们写列表推导式或堆积数组来将x转化为带常数项的X；

fit方法给出的结果包含很多内容，我们最关心的系数被放在params下。

当然，建议大家通过 命令完整地看一下回归结果。或者通过输入“results.”再点击按钮查看你想获取的指标值，比如results.rsquared_adj是调整后的R²。

方法四和方法五都来自scipy库，官方手册对它们的描述如下：

我们可以对比一下这两个方法：linregress是专门用来做线性回归的，能够汇报系数、标准误、t检验、R²；curve_fit专门用来做非线性回归，拿来做线性回归当然也可以，未免有些杀鸡用牛刀。

该命令的参数按照先x后y的顺序排列，与一般计量软件中的习惯略有不同；但它不需要为x添加常数项。返回值依次为①斜率（x系数）、②截距（常数项）、③R²、④双边t检验的P值、⑤标准误，可以用一个元组一次性解包：

非线性回归需要明确模型的形式，这就需要先定义一个表示”待估计模型“的函数。注意其中的参数顺序是先写变量、再写待估计参数：

返回值中的第一个元素是系数组成的数组，我们还是用元组来解包：

估计系数的顺序与自定义函数f中系数的顺序一致。

方法六：numpy.polyfit

方法六和方法七都来自numpy库。虽然已经用过了，但还是介绍一下引入约定为

numpy.polyfit

上文中，我们使用过 函数制造多项式，它的姊妹函数 则是专门用来估计多项式的系数。因此，该函数有一个很有意思的参数deg，来表示自变量x的最高次方。返回值为系数组成的数组，可以通过元组解包。代码如下

这里的1就是deg=1的位置参数写法，表示变量x只有零次项（常数项）和一次项。

函数不会自动为变量x添加零次项（常数项），因此需要人为添加。

返回值中的第一个元素是估计系数组成的数组，因此我们通过 就可以将其提取出来；第二个元素是残差和，因此 就是标准误。

方法八：sklearn.linear_model.LinearRegression

这个方法来自sklearn库，顾名思义这是个机器学习的方法库。 是普通最小二乘估计函数，引入约定如下

与之前的函数使用方法都不同，LinearRegression是用方法（而非函数）调用估计结果的。我们简单介绍一下思路：

首先，用 函数初始化估计方案，常用的参数包括fitintercept（是否估计常数项，默认True）、normalize（是否标准化数据，默认False）、copyX（是否备份x，默认True）等；

然后，用原地操作的fit方法（并非函数）为初始化后的估计方案添加自变量x和因变量y，常数项的问题已经在初始化中解决，因此不需要额外添加常数项；

最后，从fit过的方案中提取系数（coef_）和常数项（intercept_），coef_是一个数组，一元线性回归自然只有一个系数，索引为0。

一维数组的shape是个坑，它本质上是个n维向量而非n×1的矩阵，但是fit方法要求参数X和y是矩阵（而非向量），因此我们用列表推导式或修改shape的方式构造n×1的矩阵。显然，列表推导式慢到想哭（3.672秒），而修改shape的用时可以忽略不计。

对线性回归方法的补充

当然，线性回归的方法远远不止以上8种，比如下面的方法九。但它的效率不高，一元线性回归很少动用这样的大杀器：

用这两个函数估计回归模型的参数有些麻烦，我们需要自己写明残差的形式，并给出一个初始的估计值（乱写的）；但它的好处是对模型没有线性或多项式要求。我们用线性回归来举个例子：

首先，定义估计模型的函数形式：

参数的顺序值得我们再梳理一遍：alpha是初始的估计值，其顺序与正确估计得出的beta顺序一致，即b在前c在后，这也就决定了估计结果result中系数的结果也是如此。func_resid的参数必须按照①估计系数beta、②实际数据args的顺序排列。其他的顺序都可以改变。

大家可能会在一些2015年以前的书籍或代码中看到pandas库下的OLS方法，但该函数已经被转移至statsmodels库中，参考：https://github.com/pandas-dev/pandas/pull/11898。

9种线性回归方法的效率对比

我们将载入库的时间排除在外，可得9种方法的运算用时对比如下：

附完整代码：