机器学习之线性回归：OLS 无偏估计及相关性python分析

0 回顾

在最近的推送中，先后总结了最小二乘法的原理，两个求解方法：直接法和梯度下降，最后利用这两种思路进行了python实战。在用直接法求出权重参数时，有一个假设是某个矩阵不能为奇异矩阵。在实战中，我们发现如果它近似为奇异矩阵，然后再利用最小二乘法（OLS）去计算权重参数会出现bug。出现的是什么bug？在OLS算法的基础上应该怎么进行优化解决这个bug呢？

1 无偏估计

先看一个无偏估计的例子。工人师傅一天制造了1000个小零件，现在质检人员准备要检验这1000个件的合格数量和不合格数量，要求控制在10分钟内完成任务。如果质检员一个一个地检验这些零部件，想要在10分钟内完成对1000个零件的检验，基本是不可能的。那么，质检员想到了一个办法，他先是从这1000个零件中，随机地抽取了10个零件检验，结果9个合格，1个不合格，然后又随机选取了10个进行抽检，发现10个都是合格的，两遍过去了，用时正好10分钟，然后他说这1000个零件的合格数量为950个，不合格数量为50个。

在这个检验任务中，质检员随机地选取了两批共20个零件作为样本来检验，发现1个不合格，然后估计出总体中50个不合格，这个过程就叫做无偏估计。质检员如果想动一下手脚，然后在报表中填写40个不合格，那么这个过程就叫做有偏估计了。

造成无偏估计的误差来自于随机选择，如果运气好，样本中都是合格的，如果运气不好，可能还能多抽出几个不合格，但是整体上，质检员根据20个样本检验出1个不合格，进而估计整体的不合格数，是一个无偏估计的过程，与系统本身无关。

回想下最小二乘法，它是确确实实地无偏估计，忠实于原来的数据，根据样本估计出权重参数，进而推断新来的样本。

正因为太过忠实于数据，OLS遇到一类数据集，它们某几列具有很强的相关性，至于什么是共线性请参考本公众号的储备系列（1）。

下面来演示 OLS 如何惧怕共线性。

2 OLS 惧怕共线性

还得从直接法求解线性回归最小二乘项的权重参数公式入手，请看下面的公式，

如果X^TX 为奇异矩阵（不知道什么是奇异矩阵的请直接参考本公众号的知识储备系列机器学习：几个重要矩阵（5），那么显然X^TX 不能求逆了，自然就不可解了。

但是实战过程中，你会发现X^TX 很少正好是理论上的奇异矩阵，什么意思呢？它有可能是近似于，很接近于为奇异矩阵，比如长得像下面这个样子：

import numpy as np

import numpy.linalg as la

a = np.array([[0.9999999,1.0],[2,2]])

a

array([[ 0.9999999, 1. ],

[ 2. , 2. ]])

因此，当对上面这个矩阵求逆时，它的逆矩阵就会是一个元素值很大的矩阵，

ainv = la.inv(a)

ainv

array([[-10000000.00526356, 5000000.00263178],

[ 10000000.00526356, -4999999.50263178]])

那么得到这么大的一个矩阵，带来的影响是什么呢？

这种共线性矩阵最后导致的结果是得到一个元素取值很大的权重参数矩阵。因为线性回归的模型是 y = theta.transpose().dot(X)，因为X的系数很大，所以一个很小的样本X的扰动，会导致y的取值波动很大，这就是我们所说的方差会很大，取值不聚集，取值很散，会造成不小的误差值估计，这就是一个bug。

具体请看下面的测试。

3 python共线性测试

在Jupyter notebook中，我们快速实现测试最小二乘法直接求参数的公式的模拟，如下所示：

输入了一个矩阵 x，假定只有1个特征项，这样连上偏置项，x含有2列，再假定手上有2个样本，因此 x 是 2*2 的矩阵。它是线性相关的，x转置后得到 xt，xt的元素修改一个元素值，这样xt*x得到的矩阵为近似相关性矩阵，再求逆后看到得到一个如下数组，可以看到

这个数组的元素值非常大！

array([[ 2.00000001e+08, -9.99950004e+07],

[ -9.99950004e+07, 4.99950005e+07]])

测试代码如下所示：

import numpy as np

import numpy.linalg as la

'''

计算：(xt*x)的逆

xt是x的转置矩阵

* : 求内积

'''

'python的list转numpy的ndarray的包装函数'

x = np.array([[1.0, 2], [1, 2]])

'求转置'

xt = x.transpose()

xt[0, 0] = 0.9999

'求内积'

x2 = xt.dot(x)

'求逆矩阵'

x2inv = la.inv(x2)

'''

x2inv

array([[ 2.00000001e+08, -9.99950004e+07],

[ -9.99950004e+07, 4.99950005e+07]])

'''

xp = x2inv.dot(xt)

'标签值'

y = np.array([0.1, 0.2])

theta= xp.dot(y)

theta

'''

theta

array([ 1000.00000415, -499.90000207])

'''

可以看到得到的 theta 是一个元素取值很大的矩阵，偏置项大约为1000, 特征的权重参数大约为-500 。

因此，如下所示，当样本的取值原来为1.0，后来变为1.1后，y的变动从500到450，是500倍的变化关系。

theta = np.array([ 1000.00000415, -499.90000207])

x1 = np.array([1.0,1.0])

y1 = theta.transpose().dot(x1)

y1

500.10000207999997

x2 = np.array([1.0,1.1])

y2 = theta.transpose().dot(x2)

y2

450.11000187299987

4 总结和解决共线性

总结下上面的测试过程，当数据是一堆共线性数据时，再用OLS进行回归分析就会产生bug：属性值的扰动会与输出值间变得非常敏感，进而产生一个非常大的误差项的方差。

那么解决的办法呢？想办法调整线性相关列的权重参数，让原来线性相关的列变得不那么线性相关了，或者甚至直接过滤掉其中的某些列，都是可以解决这类多重线性相关问题的，这类算法早就被别人想到了，它们不就是 Ridge regression 和 Lasso regression 吗？

Ridge regression 是一种专用于共线性数据分析的有偏估计回归方法，它实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，获得回归系数更为符合实际、更可靠的回归方法。它的特征就是擅长回归共线性数据。

在明天的推送中，将阐述 Ridge regression 和 Lasso regression 具体怎么实施的这套delete 多重相关性的！ 然后，我们就开始一起讲讲机器学习应用非常广发的逻辑回归吧。