.. _ch_29:

########################################################
线性回归模型
########################################################

在统计分析中，经常需要分析变量之间的关系，
最常见的场景就是研究某个变量如何取决于其它一个或多个变量，
通常会把研究的目标变量称为响应变量、输出变量、被解释变量，
把影响目标变量的其它变量称为预测变量、输入变量、解释变量等。
习惯上，用符号 :math:`X` 表示输入变量，
用符号 :math:`Y`  表示输出变量，
通常二者间存在一定的关系
，统计学家研究内容就是找到二者之间合适的函数关系，
:math:`y=f(x)`
，然而二者之间最真实的关系是不得知的，
统计学家也只是尽可能的找到一个近似的关系。
线性回归模型是统计学中最基本的模型，
也是机器学习领域的入门模型，
得益于它的简单，应用十分广泛。
本书的主题，广义线性模型就是线性回归模型的扩展，
在讨论广义线性模型之前，先简单回顾一下线性回归模型。



最小二乘
########################################################

最小误差
====================================

通常情况下，机器学习就是要找到输出变量 :math:`Y` 和输入变量 :math:`X` 的函数关系，
变量之间最简单的函数关系，就是线性函数关系。
在线性回归模型中，假设 :math:`Y` 与 :math:`X` 之间是一个线性函数的关系。


.. math::
    :label: eq_linear_1

    y = \beta_0+ \beta_1 x_1 +\beta_2 x_2+\dots+ \beta_p x_p



由于输入特征数据 :math:`X` 通常是多维的，所以 :math:`f(x)` 是一个多元一次函数，
其中 :math:`x_j` 是第 :math:`j` 维输入数据，:math:`x_j` 可以是任意实数值。
:math:`\beta_j (j > 0)` 是 :math:`x_j` 的系数，:math:`\beta_0` 是这个线性方程的截距，
我们把 :math:`\beta_j` 看做这个函数的未知参数，其值暂时是未知的。
:eq:`eq_linear_1` 看上去不是很简洁，通常为了公式的简洁表达，我们令
:math:`\beta^T=[\beta_0,\beta_1,\dots,\beta_p],x^T=[1,x_1,x_2,\dots,x_p]` ，
然后把 :eq:`eq_linear_1` 看成是特征向量 :math:`x` 和参数向量 :math:`\beta` 的內积形式。
线性回归模型的简洁表示为：

.. math::
    y =  \beta_0 \times 1 + \beta_1x_1 +\beta_2 x_2+\dots+\beta_p x_p
    = x^T \beta

.. note::

    为了满足向量內积的形式，我们人为的增加一列特征数据 :math:`x_0=1`，所有输入样本的 :math:`x_0` 都为常量1。


:numref:`fg_linear_1` 是单一维度输入变量的线性回归模型的图形化展示，
为了方便图形化展示，我们假设输入特征变量 :math:`X` 只有一维。
图中蓝色的点代表一些 :math:`(x_i,y_i)` 的样本点，下标 :math:`i` 是样本点的编号。
线性回归模型的本质就是找到一条直线 :math:`y=x^T \beta` (比如图中的红色直线)
，并且这条直线和样本点的"走势"是一致的，这样我们就可以用这条直线去预测新的样。
比如根据图中蓝色样本点的分布，我们找到红色直线和样本的"走势"一致，
这样当输入一个新的 :math:`x_{new}` 时，就用直线上的点 :math:`y_{new}=x_{new}^T \beta`
作为预测值 :math:`\hat{y}=y_{new}=x_{new}^T \beta` 。


.. _fg_linear_1:

.. figure:: pictures/29_1.png
    :scale: 50 %
    :align: center


    单一输入特征的回归模型


然而空间上存在无数条直线，要如何确定和样本点"走势"相同的直线呢？
我们的最终目的是用这条之间预测新的样本点，那么理论上预测最准的直线是最优的直线。
对于一条样本数据 :math:`(x_i,y_i)` ，模型的预测值是 :math:`\hat{y}_i = x_i^T \beta`
,显然，最优的直线就是 **预测误差** 最小的直线。
我们把样本的真实值 :math:`y_i` 和预测值 :math:`\hat{y}_i` 之间的差值定义成残差(residual)，
我们的目标就是找到一条令所有观测样本残差最小的直线。
通常我们使用所有样本残差的平方和(residual sum of squares,RSS)做为整体的误差。


.. math::

    J(\beta)=\sum_{i=1}^N (\hat{y}_i-y_i)^2

.. _fg_linear_2:

.. figure:: pictures/29_2.png
    :scale: 50 %
    :align: center

    线性回归的残差



我们认为令 ``RSS`` 取得最小值的直线的最优的直线，
所以我们通过极小化 ``RSS`` 来确定这条最优的直线，
由于直线是由参数 :math:`\beta` 决定的，所以要确定这条直线就是等价于找到参数 :math:`\beta` 的值。
因此：


.. math::

    \hat{\beta} = \mathop{\arg \max}_{\beta} \sum_{i=1}^N (\hat{y}_i - y_i)^2
    =\mathop{\arg \max}_{\beta} \sum_{i=1}^N (x_i^T \beta-y_i)^2

.. note::

    通常机器学习的过程，都是先根据场景或数据定义一个参数化的模型函数 :math:`y=f(x,\beta)`
    ，模型函数是对单条数据样本的建模。但它是含有未知参数的，
    我们需要找到一个最优的参数值使得这个模型函数尽可能好的拟合数据样本。
    因此就需要定义一个评价不同参数值模型好坏的标准或者说函数，通常我们称这个评价函数为目标函数(object function)，
    然后通过极大(小)化目标函数求得参数的最优解。
    比如似然函数就是目标函数的一种，除此之外，还可以定义某种损失(误差)函数(cost function| loss function | error function)
    作为目标函数，比如线性回归的平方损失、逻辑回归的交叉熵损失等等。







**以“残差平方和最小”确定直线位置的方法被称为最小二乘法。**
用最小二乘法除了计算比较方便外，得到的估计量还具有优良特性， **这种方法对异常值非常敏感**。




参数估计
==================


显然对于线性回归模型， ``RSS`` 是一个关于 :math:`\beta` 的二次函数，
我们知道二次函数一定是存在唯一的一个极值点的，所以参数 :math:`\beta` 一定存在唯一解，
并且在极值点函数的导数为 :math:`0`。
所以我们可以直接求出 ``RSS`` 的导数，并令导数为 :math:`0` 的方法求得 :math:`\hat{\beta}` 。


.. math::

    \frac{\partial J}{\partial \beta} = \sum_{i=1}^N 2 x_i (x_i^T \beta - y_i)


上述偏导结果中包含所有样本的求和符号 :math:`\sum_{i=1}^N` ，为了简单表达我们用矩阵和向量乘积的方式替换求和符号。
我们用符号 :math:`X` 表示训练样本集中所有的输入数据的矩阵，
用符号 :math:`y` 表示训练样本集中所有输出数据的向量。上述偏导用矩阵符号表示为：

.. math::

    \frac{\partial J}{\partial \beta} = 2X^T (X \beta - y)


然后我们令导数为0，可以得到：


.. math::
    :label: eq_linear_03

    X^TX\beta = X^T y


:eq:`eq_linear_03` 通常被称为正规方程组(normal equations)，
理论上，我们可以根据这个等式得到参数 :math:`\beta` 的解析解

.. math::

    \hat{\beta}=(X^TX)^{-1} X^Ty


然而，这个解析解需要求矩阵 :math:`X^TX` 的逆矩阵，
矩阵存在逆矩阵需要满足两个条件：(1)是方阵，(2)是满秩的。
虽然是 :math:`X^TX` 方阵，但是未必满秩，
那么也就不存在逆矩阵，当逆矩阵不存在时无法计算解析解。
当矩阵 :math:`X` 存在(行或列)共线性时，:math:`(X^TX)` 一定是不满秩的。

:math:`(X^TX)` 不存在逆矩阵不代表参数 :math:`\beta` 无解，上面已经讲过损失函数 :math:`J(\beta)`
是二次函数，一定存在唯一的极值点，所以参数 :math:`\beta` 一定有全局最优解的。
当无法求得解析解时，我们可以使用迭代法求解，比如基于一阶导数的梯度下降法和基于二阶导数的牛顿法。





线性回归的概率解释
########################################################

现在我们尝试在概率的框架下解释线性回归模型。
在概率的框架下，认为输入变量 :math:`X` 和输出变量 :math:`Y` 都是随机变量，
两个变量的联合概率为

.. math::

    P(X,Y)=P(X)P(Y|X)

在回归问题中，给定输入 :math:`x` ，模型输出 :math:`y` 的值，特征变量 :math:`X` 的值是已知的、确定的，
所以不需要边缘概率 :math:`P(X)` ，只需要得到条件概率 :math:`P(Y|X)` 即可。
换句话说，不需要对联合概率 :math:`P(X,Y)` 进行建模，只需要建模条件概率 :math:`P(Y|X)`
。在概率的框架下，用条件概率 :math:`P(Y|X)=f(x)` 去表达这种关系，
即当 :math:`X=x` 时，变量 :math:`Y=y` 的概率为 :math:`P(Y=y|X=x)` 。


高斯假设
============================================


.. _fg_linear_12:

.. figure:: pictures/29_12.jpg
    :scale: 30 %
    :align: center

    线性回归模型表示成条件均值函数加上一个高斯噪声。

观察 :numref:`fg_linear_12`
，变量 :math:`Y` 的值，可以看做是线性预测器 :math:`x^T \beta` 的值加上一个误差项（噪声）
:math:`\epsilon` 。

.. math::

    y=x^T \beta + \epsilon

对于噪声，最普遍常用的就是高斯噪声，
这里假设 :math:`\epsilon` 是一个均值为 :math:`0`
，方差为 :math:`\sigma^2` 的高斯噪声，
:math:`\epsilon \sim \mathcal{N}(0,\sigma^2)`
。线性部分 :math:`x^T \beta`
是一个数值量，不是随机量。
期望为 :math:`0` 的高斯项 :math:`\epsilon` 加上一个数值项 :math:`x^T \beta`
得到的就是一个期望值为 :math:`x^T \beta` 随机变量，
因此输出变量 :math:`Y` 是一个高斯变量。

.. math::

    Y \sim \mathcal{N}(x^T \beta,\sigma^2)




因此条件概率 :math:`P(Y|X)` 是均值为 :math:`x^T \beta` 方差为 :math:`\sigma^2`
的高斯分布，
条件概率分布 :math:`P(Y|X)` 的概率密度函数为：


.. math::
    :label: eq_linear_10

    p(y|x;\beta) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \{ -\frac{1}{2\sigma^2}(y-\beta^Tx)^2 \}


.. note::

    实际上 :math:`Y` 的概率分布要根据数据的实际分布确定，并不是一定要高斯分布。
    只不过在线性回归这个模型中我们假设 :math:`Y` 是服从高斯分布的。
    如果你的数据不是(近似)高斯分布，那么就不应该使用线性回归模型。
    在后面的章节中我们会介绍当 :math:`Y` 是其它分布时，应该怎么处理，
    在广义线性模型的框架下，输出变量 :math:`Y` 可以扩展到指数族分布。


在传统线性回归模型中，方差项 :math:`\sigma^2` 被假设为常量 :math:`1` 。
:math:`\beta` 是模型的参数，是需要利用观测数据进行学习的。
高斯版的线性回归模型，是一个概率模型，可以使用最大似然估计参数 :math:`\beta`
，下一节详细阐述如何利用最大似然估计参数 :math:`\beta` 。

在得到参数 :math:`\beta` 的最大似然估计值 :math:`\hat{\beta}_{ML}`
后，可以用高斯变量 :math:`Y` 的期望值作为模型的预测值，
模型输入一个 :math:`x` 值，输出
:math:`\hat{y}=\mathbb{E}[P(Y|X=x;\hat{\beta}_{ML})]=x^T \hat{\beta}_{ML}`
。

实际上，线性回归模型比我们看上去的要广泛，线性组合部分 :math:`x^T \beta` 只要求对 :math:`\beta` 是线性的，
并不要求对 :math:`x` 是线性的，所以可以为 :math:`x` 加上一个非线性的函数 :math:`\phi(\cdot)`
使模型具有拟合非线性数据的能力。


.. math::

    y=\beta^T \phi(x) + \epsilon =\beta^T x' + \epsilon

:math:`\phi(x)` 可以看做是对特征数据的预处理 :math:`x \Rightarrow \phi(x)`
，转化之后的 :math:`x'=\phi(x)` 作为模型的输入特征，并不影响线性回归模型的定义和计算。


高斯假设的线性回归模型是建立在两个很强的假设之上的，(1)条件概率 :math:`P(Y|X)` 是高斯分布，
并且(2)不同的 :math:`x` 条件下方差是相同的。然而这种假设在很多时候是不满足的，尤其是第(2)点，
这也是线性回归的模型的局限性。

参数估计
==============================

假设我们有一个成对的观测数据集 :math:`\mathcal{D}=\{(x_i,y_i);i=1,2,\dots,N\}`
，其中 :math:`x_i` 是一条输入变量 :math:`X` 的观测样值，:math:`y_i` 是对应的输出变量 :math:`Y` 的观测值，
注意 :math:`x_i` 是一个 :math:`p` 维的向量(vector)，而 :math:`y_i` 是一个标量(scalar)。
样本集中的样本都是满足独立同分布(IID)的。
样本集的联合概率可以写成：

.. math::

    P(\mathcal{D}) &= \prod_{i=1}^N P(y_i|x_i;\beta)

    &=  {2\pi\sigma^2}^{-\frac{N}{2}} \exp \{ -\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i-x_i^T \beta)^2 \}


我们知道观测样本集的联合概率就是似然函数，我们可以通过最大似然估计法估计出模型的未知参数 :math:`\beta` ，
为了计算简单，通常我们采用极大化对数似然函数的方法估计参数。
线性回归模型的对数似然函数为：


.. math::


    \ell(\beta;x,y) = \underbrace{{-\frac{N}{2}} \ln (2\pi\sigma^2)}_{\text{常量}}
    -\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i-x_i^T\beta)^2


由于我们假设方差 :math:`\sigma^2` 是常量，所以上述公式的第一项是一个常量，
在极大化对数似然函数时不影响最终的求解，所以是可以去掉的。


.. math::

    \hat{\beta}_{ML} &={\arg \max}_{\beta}  \ell(\beta;x,y)

    &\triangleq  {\arg \max}_{\beta}
    \left \{ -\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i-x_i^T \beta)^2
    \right \}



我们发现这和最小二乘法的损失函数 :math:`J(\beta)` 是等价的，
对数似然函数中的 :math:`\sum_{i=1}^N (y_i-x_i^T \beta)^2` 就是残差的平方和(residual sum of squares,RSS)。
同理，我们可以使用迭代法求的最优解。