精通数据科学(5)第4章 线性回归:模型之母(1)

第4章 线性回归:模型之母

道生一,一生二,二生三,三生万物。

——老子

线性回归是最重要的数学模型之一

线性回归模型看上去很简单,简单到让不少人觉得它并没有什么研究和使用的价值。其实并不是这样的,线性回归可以说是最重要的数学模型之一,其他很多模型都建立在它的基础之上。

1.一个有关数学家的笑话

一天,数学家觉得自己已受够了数学,于是他跑到消防队去宣布他想当消防员。

消防队长说:“您看上去不错,可是我得先给您一个测试。”消防队长带数学家到消防队后院小巷,巷子里有一个货栈、一只消防栓和一卷软管。

消防队长问:“假设货栈起火,您怎么办?”

数学家回答:“我把消防栓接到软管上,打开水龙,把火浇灭。”消防队长说:“完全正确!最后一个问题:假设您走进小巷,而货栈没有起火,您怎么办?”

数学家疑惑地思索了半天,终于答道:“我就把货栈点着。”消防队长大叫起来:“什么?太可怕了!您为什么要把货栈点着?”

数学家回答:“这样我就把问题化简为一个我已经解决过的问题了。”

建模型的思路和笑话中数学家的思路一样。当遇到一个新问题时,总是考虑通过某种数学上的变换,将未知问题转化为已知模型能解决的问题。

所以任何一个复杂模型,一层层拨开来,里面可能藏着好多个线性回归模型。因此,线性回归模型很有研究的必要。深入了解它的模型细节能帮助我们理解其他模型,进而指引我们根据实际场景搭建有效的模型。

2.线性回归的使用价值

读者可能会觉得,线性回归太过简单,应该没人会使用它来解决实际问题吧。但事实上,在经济学领域,有很多模型都是线性回归模型。

比如宏观经济学的基石之一——IS-LM模型 。这些模型的结果深刻地影响着政府的经济决策,进而影响着我们生活的方方面面。因此可见,线性回归并不缺乏实用例子。

IS-LM模型由两个线性模型组成:

  • 一个是IS模型,即投资-储蓄模型(Investment-Saving),它研究的是国民储蓄和市场利率之间的关系;
  • 另一个是LM模型,即流动性偏好-货币供给模型(Liquidity Preference-Money Supply),它研究的是国民收入与货币需求量之间的关系。

这两个模型都只涉及两个变量:

  • $Y$表示国民收入,可以简单地理解为日常生活中常听到的GDP;
  • $r$表示利率。

对于投资储蓄行为,假设在短期内,国民收入$Y$ 增加,这意味着人们可以用于储蓄的资金增加,市场上的利率会下降。因此在IS模型中, $r=-aY+b$,其中,$a>0$ 。

对于货币供给,假设在短期内,央行的货币供应量是不变的,而且国民收入$Y$ 增加。由于收入增加,人们对货币的需求量是增加的,而货币供应量并没有增加,这会抬高市场上的利率。因此在LM模型中,$r=cY+d$ ,其中$c>0$ 。

将这两个结合起来,就可以得到市场上的利率以及对应的国民收入。更进一步,根据想要达到的利率和国民收入可以倒推出央行的货币供应量。IS-LM模型充分说明:“模型不在复杂,有用则灵”。

而且搭建模型并不是要完全模拟现实世界,而是一个不断做近似的过程。这就是为什么数据科学家们常说“所有模型都是错的”。当然这句话还有后半句“但是,其中一些是有用的”。

一个“有用”的模型能过滤掉数据中那些不重要的细枝末节,抓住其中主要的内在关系,从而帮助我们更好地理解和解释数据。而在很多情况下,线性模型就是这么一个符合要求的“有用”模型,而且它还简洁、高效、容易理解。

既然线性模型已经够用了,那又何必再费劲去构造那些难以理解,同时又极易引入其他新问题的复杂模型呢 ?所以线性回归模型在现实世界里有很高的使用价值。

这里涉及一个很重要的哲学准则,它叫作奥卡姆剃刀(Occam's Razor),也称为简约之法则。这个准则在机器学习领域应用得很广,其主要含义可概括为 “如无必要,切勿假定繁多”。 也就是说,如果关于同一个问题有许多种理论,每一种都能做出同样准确的预言,那么应该挑选其中使用假定最少的。尽管越复杂的方法通常能做出越好的预言,但是在不考虑预言能力(即结果大致相同)的情况下,假设越少越好。在建模的过程中,我们大多遵循这一原则。

奥卡姆剃刀最成功的案例来自于天文学。在日心说刚刚诞生时,太阳、月亮和其他太阳系行星的运动既可以用地心说来解释,也可以用日心说来解释。在当时看来,两种假说都同样有效,然而,日心说只需要7个基本假设,地心说却需要很多的假设。那么基于奥卡姆剃刀原则,应该选择日心说。

一个简单的例子

为了更形象地描述,假设我们有一个制作手工玩偶的朋友小潘,小潘想分析一下他制作玩偶的数量和成本之间的关系。于是他一边生产玩偶,一边记录了如图所示的数据。

小潘仔细研究了一下他所记录的数据,似乎玩偶个数和成本是线性关系。这也很符合小潘的预期,为了更加直观地了解数据,验证他的猜测,他决定将数据可视化。他把上面的数据点表示在一个直角坐标系里,如图所示。

根据图像,小潘发现,生产成本和生产个数并不呈一条严格的直线,似乎是在沿着某条直线上下随机地波动。

其实上面展示的这些数据,是由“自然之力”按照下面这个数学公式来产生的。

$$y_i=x_i+\varepsilon_i$$

这个公式其实并不严谨。因为正态分布的值域是整个实数,所以根据公式,$y_i$ 的取值可能为负数。这与生产成本必然大于0相矛盾。数学上完全严格的写法应为 $y_i=\max(x_i+\varepsilon_i,0)$ 。但是在书中这个场景里, 的概率非常小,几乎可以忽略,所以模型可近似为上述公式。

$x_i$是某一天小潘制作的玩偶个数, $y_i$是对应那天的生产成本。这个数学式子告诉我们,小潘制作玩偶的平均成本是1。

其中$\varepsilon_i$ 是一个随机变量,它服从期望为0,方差为1的正态分布。它表示在小潘生产玩偶时,一些随机产生或随机节约的成本。比如制作失败的玩偶,这时$\varepsilon_i$ 为正;又比如制作过程中刚好发现有可用的旧布料,这时候 $\varepsilon_i$为负。

$\varepsilon$所代表的随机成本和制作玩偶的个数是相互独立的。 但是,小潘并不是控制公式的“自然之力”。他所看到的只是一堆数据:玩偶个数 和对应的生产成本 。而他看不到的公式正是他想要的结果。

小潘并不知道如何对数据建模分析,所以他将问题提给了他的两个好朋友:一位是从事机器学习工作的小陈,另一位是进行统计工作的小郭。我们来看看这两位分别打算怎么解决这个问题呢?

从机器学习的角度看这个问题

1.确定场景类型

(1)在现有的数据集里,我们有玩偶的生产个数(记为 $x_i$)和生产成本(记为 $y_i$)。其中 $i$表示第$i$ 天的数据,比如 $x_i$表示第 $i$天生产的玩偶数。而我们需要通过生产个数的信息去预测生产成本。在数据里面,已经有需要被模型预测的量:生产成本,所以这是一个监督式学习。

(2)需要被预测的成本 $y_i$是一个数量。它是一个连续变化的量,而并非表示类别的离散量,所以这是一个回归问题。

2.定义损失函数(loss function)

(1)搭建模型的目标是使得模型预测的生产成本和实际成本接近。

(2)我们把上面这句话翻译成数学语言:已知实际成本为 $y_i$,假设模型预测的成本为$\hat{y_i}$ 。那么预测值和真实值的差距就定义了一个损失函数 $LL$。而搭建模型的目标就是使得损失函数达到最小值

$$LL=\sum_i |y_i-\hat{y_i}| \tag {4-2}$$

(3)上述公式并不是一个处处可导的函数,数学上处理起来比较麻烦。因此,我们重新定义一个数学上容易处理的损失函数,如下公式其实就是真实值与预测值之间的欧式距离平方和。模型参数的估计就依赖于这个损失函数:

$$L=\sum_i (y_i-\hat{y_i})^2 \tag {4-3}$$

对于线性回归模型,公式(4-3)所定义的损失函数除了在数学上更易处理外,它还隐含地满足了数据的一个假设,即随机扰动项服从正态分布。

事实上,公式(4-2)和公式(4-3)对应着不同的线性回归模型。

  • 公式(4-2)对应着Least Absolute Deviations Regression,
  • 公式(4-3)对应着Least Squares Regression。

3.提取特征

(1)数据是小潘自己收集记录的。经过检查,里面并没有记错或者特别异常的数据。换句话说,数据是“干净”的,可以直接使用。

(2)现有的数据里,只有一个原始特征 $X={x_i}$,它表示个数。在这个变量上面的加减乘除是有明确含义的。也就是说变量本身的数学运算是有意义的,所以 $X$可以直接在模型里面使用。

(3)当然也可以对$X$ 做某种数学变换,得到一个新的特征。比如对它做平方运算,得到新的特征 $X^2$。这些新提取的特征也可以被应用到模型里。但对小潘的问题,我们先只用原始特征建模。如果效果不好,则再考虑提取新的特征。

4.确定模型形式并估计参数

(1)根据小潘的分析,$x_i$ 和 $y_i$之间是线性关系,因此可以直接使用线性模型,不需要考虑非线性问题到线性问题的转化了。

(2)模型的定义如公式(4-4),其中, $a$表示生产一个玩偶的变动成本,$b$ 表示生产的固定成本 [5] 。

$$\hat{y_i}=ax_i+b \tag {4-4}$$

如果用模型的语言来描述,公式(4-4)中的 $b$称为截距,它是线性回归模型中非常重要的参数。在加入这个参数后,模型才能在数据平移和数据单位转换时保持稳定。具体来讲:

  • 当发生数据平移时,比如小潘发现,他原来统计的成本中少算了1元,现在需要加上,即 $y_i$将变为$y_i+1$ 。参数$a$ (平均成本)的估计值不变。
  • 当发生数据单位转换时,比如小潘将成本的单位由元改为了百元,即$y_i$ 将变成$\frac{y_i}{100}$ 。参数 $a$将成比例的放大或缩小,比如上面的例子里, $a$将变成$\frac{a}{100}$ 。
  • 上面的两点性质对线性回归,特别是多元线性回归十分重要。因为当模型有多个自变量时,它们的单位常常是不一致的。如果没有加入 $b$这个参数,这两点性质便不再成立。所以在搭建线性回归模型时,通常会在模型中加入截距。除非给定的数据集里已经包含了一个常量自变量。

(3)如上面提到的那样,参数 $(a,b)$ 的估计值$(\hat{a},\hat{b})$ 将使得损失函数 达到最小值,如公式(4-5)所示。

$$(\hat{a},\hat{b})=\argmin_{a,b}\sum_i(y_i-ax_i-b)^2 \tag {4-5}$$

5.评估模型效果

(1)从预测的角度看,我们希望模型的预测成本越接近真实成本越好。所以如公式(4-6)所示,定义线性模型的均方差 。均方差越小,模型效果越好。

$$MSE=\frac{1}{n}\sum_{i=1}^n(y_i-\hat{y_i})^2=\frac{1}{n}L \tag{4-6}$$

(2)从解释数据的角度来看,我们希望模型能最大程度地解释成本变化的原因。换句话说,未被模型解释的成本$(y_i-\hat{y_i})$占成本变化$(y_i-\frac{1}{n}\sum y_i$的比例越小越好。因此如公式(4-7)所示,定义模型的决定系数 。决定系数越接近1,模型的效果越好。

$$\overline{y}=\frac{1}{n}\sum_{i=1}^ny_i$$

$$SS_{tot}=\sum_i(y_i-\overline{y})^2 $$

$$SS_{res}=\sum_i(y_i-\hat{y})^2 $$

$$R^2=1-\frac{SS_{res}}{SS_{tot}} \tag {4-7}$$

决定系数的英文名为 coefficient of determination,通常记为 $R^2$,如图所示。

通过数学运算,在最小化损失函数时,有 $\sum_iy_i=\sum_i\hat{y_i}$。由此可以得到:

$$\sum_i(\hat{y_i}-\overline{y})^2 + \sum_i(y_i-\hat{y_i})^2=\sum_i(y_i-\overline{y})^2$$

其中:

  • $SS_{tot}=\sum_i(y_i-\overline{y})^2$, 是因变量的方差,
  • $SS_{reg}=\sum_i(\hat{y_i}-\overline{y})^2$ 是已被模型(或者说是自变量)解释的方差。所以在这种情况下,有:

$$R^2=\frac{SS_{res}}{SS_{tot}} \tag {4-7}$$

用学术语言来解释,这个公式说明:决定系数等于因变量的方差中可由自变量解释的比例。因此,可以用这个指标来判断模型的解释力。

从统计学的角度看这个问题

1.假设条件概率

(1)根据小潘的描述,数据集里有两个变量,一个是自变量 [9] 玩偶个数(记为$x_i$ ),一个是因变量生产成本(记为$y_i$)。其中$i$表示第$i$ 天的数据,比如$x_i$表示第$i$天生产的玩偶数。而且根据小潘前面的分析,这两者之间似乎是线性关系,但又带着一些随机波动。因此可以假设$y_i$ 和$x_i$ 之间的关系如下:

$$y_i=ax_i+b+\varepsilon_i \tag {4-8}$$

细心的读者会注意到,对于数据里的玩偶个数$x_i$,小陈使用“特征”这个术语来描述它,而小郭使用的术语是“自变量”。这两个术语的含义并无差异,只是 “特征”这个术语常用于机器学习领域,而“自变量”这个术语常用于统计学领域。

(2)其中 $a,b$是模型的参数,分别表示生产一个玩偶的变动成本和固定成本;而 $\varepsilon_i$被称为噪声项,表示没被已有数据所捕捉到的随机成本。它服从期望为0,方差为 $\sigma^2$( $\sigma^2$也是模型的参数)的正态分布,记为 $\varepsilon_i \sim N(0,\sigma^2)$ 。这里假设${\varepsilon_i}$ 之间相互独立,而且 ${\varepsilon_i}$ 和${x_i}$ 之间也是相互独立的,这两点假设非常重要 。如果假设不成立,情况会怎样呢?本书将在第7章中详细讨论相关问题。

(3)从左到右看公式(4-8),如果给定一组参数$a,b$ 以及噪声项的方差 $\sigma^2$ 。由于 $x_i$表示玩偶个数,是一个确定的量(小潘能控制他生产的玩偶数)。那么 $y_i$ 就和 $\varepsilon_i$ 一样是一个随机变量,服从期望为 $ax_i+b$,方差为 $\sigma^2$正态分布,即 $y_i \sim N(ax_i+b,\sigma^2)$。换句话说,小潘提供的数据 $y_i$,只是$N(ax_i+b,\sigma^2)$ 这个正态分布的一个观测值,如图所示,而且 ${y_i}$之间也是相互独立的。

(4)把上面的第(3)点翻译成数学语言就是: $y_i$在已知 $a,b,\sigma$时的条件概率是 $N(ax_i+b,\sigma^2)$ ,如公式(4-9)所示。

$$P(y_i|a,b,x_i,\sigma^2\sim N(ax_i+b,\sigma^2) \tag{4-9}$$

2.估计参数

(1)根据上面的分析, ${y_i}$之间也是相互独立的。所以得到 ${y_i}$出现的联合概率如公式(4-10)所示。在学术上,这个概率被称为模型的似然函数(likelihood function),通常也被记为$L$ 。

$$P(Y|a,b,X,\sigma^2)=\prod P(y_i|a,b,x_i,\sigma^2)$$

$$\ln P(Y|a,b,X,\sigma^2)=-0.5n\ln(2\pi\sigma^2)-(\frac{1}{2\sigma^2})\sum_i(y_i-ax_i-b)^2 \tag{4-10}$$

在学术上,通常用 $\theta$ 表示模型中的所有参数,比如正文中的 $a,b,\sigma^2$;用 $D$表示数据。则参数的似然函数记为 $L(\theta|D)=P(D|\theta)$。值得注意的是,等式右边的$P(D|\theta)$ 表示条件概率,而等式的左边 $L(\theta|D)$,虽然形式与右边相似,但不表示条件概率。这一点通常会让初学者感到疑惑。

(2)对于不同的模型参数, ${y_i}$出现的概率(即参数的似然函数)并不相同。这个概率当然是越大越好,所以使这个概率最大的参数将是参数估计的最佳选择。此方法也被称为最大似然估计法(Maximum Likelihood Estimation,MLE)。根据公式(4-10),参数 $(a,b)$的估计值 $(\hat{a},\hat{b})$如下

$$(\hat{a},\hat{b})=\argmax_{a,b}P(Y|a,b,X,\sigma^2)=\argmin_{a,b}\sum_i(y_i-ax_i-b)^2 \tag{4-11}$$

(3)同理,可以得到参数$\sigma^2$ 的估计值$\hat{\sigma^2}$ ,如公式(4-12)所示。

$$\hat{\sigma^2}=\argmax_{a,b}P(Y|a,b,X,\sigma^2)=\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{n} $$

$$\hat{y_i}=\hat{a}x_i+\hat{b} \tag{4-12}$$

3.推导参数的分布

(1)其实上面得到的参数估计值 $(\hat{a},\hat{b},\hat{\sigma^2})$ 都是随机变量。具体的推导过程有些繁琐,限于篇幅,这里略去数学细节。仅以 $\hat{a}$为例,用一个不太严谨的数学推导来说明这个问题 [13] 。

假设数据集里只有两对数据 $(x_k,y_k),(x_l,y_l)$,且 。这时我们可以通过解公式(4-13)中的方程组来得到$\hat{a}$ 的表达式。公式(4-13)的右半部分表示 $\hat{a}$ 是一个随机变量,而且服从一个以参数真实值 $a$为期望的正态分布 。

$$\begin{cases} y_k=\hat{a}x_k+\hat{b}\ y_l=\hat{a}x_l+\hat{b} \end{cases}\ \Longrightarrow \ \hat{a}=\frac{y_k-y_l}{x_k-x_l}\ =\frac{a(x_k-x_l)+\varepsilon_k-\varepsilon_l}{x_k-x_l}\ =a+\frac{\varepsilon_k-\varepsilon_l}{x_k-x_l} \tag{4-13} $$

严谨的数学证明如下:用矩阵来表示线性回归模型。令: $$Y= \begin{pmatrix} y_1\ y_2\ \vdots\ y_n \end{pmatrix}, X= \begin{pmatrix} x_1\ x_2\ \vdots\ x_n \end{pmatrix},\beta= \begin{pmatrix} a\ b \end{pmatrix},\varepsilon= \begin{pmatrix} \varepsilon_1\ \varepsilon_2\ \vdots\ \varepsilon_n \end{pmatrix}$$ 则模型可表达为 $Y=X\beta+\varepsilon$。而参数 $\beta$ 的估计为 $\hat{\beta}=\argmin_{\beta}(Y-X\beta)^2$ 。进一步的数学运算可以得到: $$\hat{\beta}=(X^TX)^{-1}X^TY=(X^TX)^{-1}X^TX\beta + (X^TX)^{-1}X^T\varepsilon$$ 这个式子表示,模型参数的估计值是随机变量,而且它们服从一个期望为参数真实值的正态分布。

(2)通过更加细致的数学运算可以得到:

$$\hat{b}\sim N(b,\frac{\sigma^2}{n})$$

$$\hat{a}\sim N(a,\frac{\sigma^2}{\sum_i(x_i-\overline{x})^2})$$

$$\hat{\sigma^2} \sim \mathcal{X}_{n-2}^2\frac{\sigma^2}{n} \tag{4-14}$$

(3)既然参数估计值都是随机变量,那么我们更需要关心的是这些估计值所服从的概率分布,而不仅仅是根据公式(4-11)和公式(4-12)得到的数值。因为这些数值只是对应分布的一次观测值,它们并不总是等于真实参数,而是严重依赖于估计参数时所使用的数据。比如针对小潘提供的数据集,只使用前3天数据估计出来的参数和使用4~6天数据估计出来的参数就不一样,如图所示。

统计方法和机器学习方法一样严重依赖所使用的数据。数据即模型。

(4)公式(4-14)告诉我们,参数估计值的方差随着数据量的增大而减少。换句话说,数据量越大,模型估计的参数就越接近真实值。这也是大数据的价值之一:数据量越大,模型预测的效果就越好,如图所示。

4.假设检验与置信区间

(1)参数的概率分布可以透露许多很有用的信息。比如在95%的情况下,参数 $a$的真实值(生产一个玩偶的变动成本)会落在一个怎样的区间里?在学术上它被称为参数 $a$的95%置信区间。类似地,也可以计算模型预测结果的置信区间,即对于被预测对象,真实值的大致范围是怎样。这一点非常重要,因为模型几乎不可能准确地预测真实值。知道真实值的概率分布情况,能使我们更有信心地使用模型结果。

(2)又比如在1%犯错的概率下,我们能不能拒绝参数$b$ 的真实值其实等于0这个假设?在学术上它被称为参数 $b$的99%显著性假设检验。

对于这个假设检验,换个思路更通俗一点理解:参数$b$ 的真实值等于0的概率是否小于1%?这可以帮助我们更好地理解数据之间的关系,比如小潘生产玩偶时,固定成本(参数$b$ )是否真实存在?或者模型估计的固定成本 $\hat{b}$只是由于模型搭建得不准确而导致的“错误”结论?

小郭分析完问题之后,就着手利用Python和第三方库Statsmodels编写自己的解决方案。

© Licensed under CC BY-NC-SA 4.0

别向医生和律师提供错误的消息。—— 本杰明·富兰克林

发表我的评论
取消评论
表情

Hi,您需要填写昵称和邮箱!