从平面线性拟合谈起

线性回归是最简单的一种数据拟合。 我们举平面上的例子来看，平面上有若干个样本点，我们的目标就是去画一条直线去拟合这些样本点。

当然这个图只是为了简化场景，让我们更直观地回忆起线性回归的概念，实际上对于样本点的维数，是不会仅仅限制于 1 维的，而是可以为任意的 $p$ 维。

拟合数据扩展到 $p$ 维

这里我们还是正式引入一些符号标记：

我们有一组样本：

$x_i$ 是一个 $p$ 维向量，它可以表达第 $i$ 个样本被观察的 $p$ 个特征，而 $y_i$ 则表示第 $i$ 个样本的取值，是一个1维的数值。拟合出来的直线也就是样本的所有 $p$ 维特征值到最终样本取值 $y_i$ 的一种近似的线性映射关系。

因此 $N$ 个样本集合写作：

而每一个样本 $x_i$ 则包含 $p$ 维特征，同样写成向量：

因此我们的目标就是通过估计出一个系数向量$w=\begin{bmatrix}w_1\w_2\w_3\…\w_p\end{bmatrix}$ 和一个偏置$b$（也就是一个数值）来近似地建立起一种如下的线性映射关系：

为了进一步简化描述，我们把偏置 $b$ 也写作是 $w_0x_{i0}$ 的形式，其中 $x_{i0}$ 恒为 1，那么上述两个向量就更新为：

原有的映射关系 $w^Tx+b$ 就简化成了 $w^Tx$。

最小二乘参数估计

那么，如何去估计这种参数向量 $w$ 就成了我们将要面对的问题。在最小二乘估计法中，定义了如下的一个目标函数：

最小二乘估计的目标就是找到一个系数向量 $w$，使得 $L(w)$ 的取值最小，换句话说就是针对这 $N$ 个样本，使得线性拟合出来的估计值和真实值之间误差的平方和最小。

从高斯噪声的角度看最小二乘

我们思考一个问题，通过线性拟合，我们能否让直线精确无误的通过每一个样本点，使得拟合出来的直线的误差为 0？

这种理想的情况在现实当中显然是不可能的，因为样本数据本身是带有噪声的，带有随机性的，这个直观点说，我们看上面那幅图，所有的样本点都围绕着直线上下，在直线附近上下夹杂着这种随机的噪声。

那么我们可以这么理解：拟合出来的直线代表了样本数据的确定性，而拟合值和真实值之间的误差，也可以称之为是噪声，则代表了随机性。

这个随机性的噪声我们是不是就可以用符合某个分布的随机变量来描述，没错，这就是我们一开始讲的高斯噪声，我们令噪声 $\epsilon$ 服从均值为0，方差为 $\sigma^2$ 的一维高斯分布，记作：

那么，直线的拟合值加上上面定义的随机噪声值，就等于样本的真实值。即：

那么很显然，在 $w$ 和 $x$ 确定的情况下，$y$ 也服从一个正态分布：

这是一个条件概率，写成概率密度函数的形式就是：

那么此时，我们手上有一组 $N$ 个独立同分布的样本集$D={(x_1,y_1),(x_2,y_2),(x_3,y_3),…,(x_n,y_n)}$，如何去估计参数$w$，很显然，使用我们第一节中介绍的极大似然估计。

我们的对数似然函数是：

那么待估计的参数 $w$ 的极大似然估计值就是：

其中，去掉所有与 $w$ 无关的常数项，那么就能迅速地化简为：

到此处已经豁然开朗，最终我们不就到达了最小二乘法的定义式了吗？二者是不谋而合。

因此，从概率的角度来说，最小二乘估计其中隐含了样本噪声服从0均值高斯分布的假设。

那剩下来去寻找 $w_{mle}$ 的过程，思路上非常简单，就是让 $L(w)=\sum_{i=1}^N(y_i-w^Tx_i)^2$ 对 $w$求导即可，但是这其中的运算涉及到向量的求导和一些矩阵向量的运算，包含一些小的技巧，我们一起来一步一步地看一下：

我们先处理向量：

我们进行一些简单处理：

而后面的 $\begin{bmatrix}w^Tx_1-y_1\w^Tx_2-y_2\…\w^Tx_N-y_N\end{bmatrix}$是其转置，因此：

最终，我们得到了：

后面的就是机械的求取 $L(w)$ 对向量 $w$ 的偏导。

首先我们来化简 $L(w)$：

一看这里是不是有点懵，不过仔细观察一下，由于 $L(w)$ 最终的结果是一个标量实数，那么式子中的每一项也必然是一个实数。

我们发现式子的第二项 $w^TX^TY$ 和第三项 $Y^TXw$，其实是互为转置的关系，加之二者都是实数，那么我们可以得出 $w^TX^TY=Y^TXw$的相等关系，因此式子得以化简：

由于 $Y^TY$ 是一个与 $w$ 无关的常量，那么最终向量 $w$ 的导数就如下进行计算，这里面涉及了一些向量求导的基本运算。

\(\frac{\partial}{\partial w}L(w)=\frac{\partial}{\partial w}(w^TX^TXw-2w^TX^TY)\)\(=2X^TXw-2X^TY=0\)

费尽周折，最终我们我们得到了用于估计向量 $w$ 的表达式：

这其中，$(X^TX)^{-1}X^T$ 又被称之为是矩阵 $X$ 的伪逆。

线性回归代码实践

最终，我们来快速看一下 Python 当中做线性回归的实现方法，我们借助 sklearn 工具包，实现起来是相当容易的，这里做一个简单示例：

这里我们举一个大家耳熟能详的例子，就是 sklearn 中自带的波士顿房价问题的，波士顿房价数据集中，x 拥有 13 维特征，对应一个数值型的房价数据 y，我们直接利用 sklearn 模块就能加载波士顿的房价数据，我们简单地看一下代码。

代码片段：

    from sklearn import datasets
    from sklearn import linear_model
    from sklearn.metrics import r2_score
    from sklearn.datasets import load_boston
    boston = load_boston()
    from sklearn.model_selection import train_test_split
    
    #获取波士顿房价数据集（506*13）
    boston = datasets.load_boston()
    
    #选取一半的数据做训练数据，一半做测试数据
    X_train,X_test,y_train,y_test = \
        train_test_split(boston.data,boston.target,test_size=0.5,random_state=33)
    
    #使用最小二乘线性回归进行拟合，导入相应的模块
    lr = linear_model.LinearRegression()
    lr.fit(X_train, y_train)   #进行线性拟合
    y_pred = lr.predict(X_test)#得到预测值集合y_pred
    
    w = lr.coef_               #获得该回该方程的回归系数与截距
    w0 = lr.intercept_
    print("估计系数w={}：".format(w))
    print("估计截距w0={}：".format(w0))
    #利用评价指标 r2_score，来评价模型的拟合优度，越接近 1 拟合情况越好
    print("拟合优度r2_score={}".format(r2_score(y_test, y_pred)))

运行结果：

    估计系数 w=[-1.13886467e-01  5.57751443e-02  3.19967368e-02  4.83656552e+00
     -1.73802992e+01  3.81631385e+00  8.77697010e-03 -1.48553776e+00
      3.16012600e-01 -1.14998706e-02 -8.54089583e-01  8.31673898e-03
     -6.00625497e-01]：
    估计截距w0=34.525556152532474：
    拟合优度r2_score=0.7206754178645804

结合注释，我们会发现利用 sklearn 的工具做最小二乘法的线性拟合还是相当容易的，最终如我们所愿，样本数据有 13 维特征，因此我们得到了一个 13维的系数矩阵 $w$ 和一个截距数值 $w_0$。