数据科学学习指南（五）

原文：zh.annas-archive.org/md5/9950f82165bee426119f99afc3ab612d

译者：飞龙

协议：CC BY-NC-SA 4.0

第五部分：线性建模

第十五章：线性模型

在本书的这一部分，我们已经以不同的程度涵盖了数据科学生命周期的四个阶段。我们已经讨论了如何提出问题、获取和清理数据，并使用探索性数据分析来更好地理解数据。在本章中，我们将常数模型（在第四章中引入）扩展为线性模型。线性模型是生命周期的最后阶段中的一种流行工具：理解世界。

了解如何拟合线性模型为各种有用的数据分析打开了大门。我们可以使用这些模型进行预测——例如，环境科学家开发了一个线性模型，根据空气传感器测量和天气条件预测空气质量（参见第十二章）。在那个案例研究中，理解两个仪器测量值的变化帮助我们校准廉价传感器，并改善它们的空气质量读数。我们还可以使用这些模型来进行推断，例如，在第十八章中我们将看到兽医如何使用线性模型推断驴的体重与长度和胸围的系数：L e n g t h + 2 × G i r t h − 175。在那个案例研究中，该模型使得在现场工作的兽医能够为生病的驴子开具药物处方。模型还可以帮助描述关系并提供见解——例如，在本章中，我们探讨了与上升流动性相关的因素之间的关系，如通勤时间、收入不平等和 K-12 教育质量。我们进行了描述性分析，按照社会科学家用来塑造公众对话和制定政策建议的分析方法。

我们首先描述简单线性模型，它总结了两个特征之间的关系，并用一条直线表示。我们解释如何使用在第四章介绍的损失最小化方法来拟合这条直线到数据上。然后我们介绍多元线性模型，它使用多个其他特征来模拟一个特征。为了拟合这样的模型，我们使用线性代数，并揭示用平方误差损失拟合线性模型背后的几何学。最后，我们涵盖了特征工程技术，这些技术可以在构建模型时包括分类特征和转换特征。

简单线性模型

像常数模型一样，我们的目标是通过常数近似特征中的信号。现在我们有来自第二特征的额外信息来帮助我们。简而言之，我们希望利用第二特征的信息制定比常数模型更好的模型。例如，我们可能会通过房屋的大小描述其销售价格，或者根据驴的长度预测其重量。在这些例子中，我们有一个结果特征（销售价格，重量），我们希望用解释变量特征（房屋大小，长度）来解释、描述或预测。

注意

我们使用结果来指代我们试图建模的特征，解释变量来指代我们用来解释结果的特征。不同领域已经采用了描述这种关系的惯例。有些人称结果为因变量，解释变量为自变量。其他人使用响应和协变量；回归和回归器；被解释和解释性；内生和外生。在机器学习中，目标和特征或预测和预测因子是常见的。不幸的是，许多这些对都暗示了因果关系。解释或预测的概念并不一定意味着一个因素导致了另一个因素。特别令人困惑的是独立-依赖的用法，我们建议避免使用它。

我们可能会使用的一个可能模型是一条直线。数学上来说，这意味着我们有一个截距θ 0和一个斜率θ 1，并且我们使用解释特征x来通过直线上的一个点来近似结果y：

y ≈ θ 0 + θ 1 x

随着x的变化，对y的估计也会变化，但仍然落在直线上。通常情况下，估计并不完美，使用模型会有一些误差；这就是为什么我们使用符号≈来表示“大约”。

要找到一条能够很好地捕捉结果信号的线，我们使用了在第四章介绍的相同方法，并最小化平均平方损失。具体来说，我们按照以下步骤进行：

1. 找到误差：y i − ( θ 0 + θ 1 x i ) , i = 1 , … , n

2. 平方误差（即使用平方损失）：[ y i − ( θ 0 + θ 1 x i ) ] 2

3. 计算数据的平均损失：

   1 n ∑ i [ y i − ( θ 0 + θ 1 x i ) ] 2

为了拟合模型，我们寻找能够给出最小平均损失的斜率和截距；换句话说，我们最小化均方误差，简称 MSE。我们称截距和斜率的最小化值为θ ^ 0和θ ^ 1。

注意，在步骤 1 中计算的误差是沿垂直方向测量的，这意味着对于特定的x，误差是数据点( x , y )与线上点( x , θ 0 + θ 1 x )之间的垂直距离。图 15-1 展示了这一概念。左侧是点的散点图，带有用于从x估计y的直线。我们用方块标记了两个特定点，并用菱形表示它们在直线上的估计值。从实际点到直线的虚线段显示了误差。右侧的图是所有误差的散点图；作为参考，我们还标记了左图中两个方块点对应的误差在右图中也用方块标记了。

图 15-1. 左侧是( x i , y i )对的散点图和我们用于从x估计y的直线。具体的两个点用方块表示，它们的估计值用菱形表示。右侧是误差的散点图：y i − ( θ 0 + θ 1 x i )。

在本章的后面，我们推导出使均方误差最小化的值 θ ^ 0 和 θ ^ 1。我们表明这些值分别为：

θ ^ 0 = y ¯ − θ ^ 1 x ¯ θ ^ 1 = r ( x , y ) S D ( y ) S D ( x )

在这里，x 表示值 x 1 , … , x n，而 y 类似定义；r ( x , y ) 是 ( x i , y i ) 对的相关系数。

将两者结合起来，线性方程变为：

θ ^ 0 + θ ^ 1 x = y ¯ − θ ^ 1 x ¯ + θ ^ 1 x = y ¯ + r ( x , y ) S D ( y ) ( x − x ¯ ) S D ( x )

这个方程有一个很好的解释：对于给定的 x 值，我们找出其高出（或低于）平均值多少个标准偏差，然后预测（或解释，具体取决于环境） y 将是 r 倍的标准偏差高出（或低于）其平均值。

从最优线的表达式中我们看到 样本相关系数 发挥了重要作用。回想一下，r 衡量了线性关联的强度，定义如下：

r ( x , y ) = ∑ i ( x i − x ¯ ) S D ( x ) ( y i − y ¯ ) S D ( y )

以下是有助于我们拟合线性模型的相关性的几个重要特征：

• r 是无量纲的。注意，x，x ¯ 和 S D ( x ) 都有相同的单位，所以下面的比率是无单位的（涉及 y i 的项同理）：

  ( x i − x ¯ ) S D ( x )

• r 介于 − 1 和 + 1 之间。只有当所有点恰好位于一条直线上时，相关性才为 + 1 或 − 1 ，具体取决于线的斜率是正还是负。

• r 衡量线性关联的强度，而不是数据是否具有线性关联。图 15-2 中的四个散点图都有约为 0.8 的相同相关系数（以及相同的平均值和标准偏差），但只有一个图，即左上角的那个，具有我们认为的带有随机误差的线性关联。

图 15-2。这四组点，称为安斯康姆的四分位数，具有相同的相关性 0.8，以及相同的均值和标准差。左上角的图表展示线性关联；右上角显示完美的非线性关联；左下角除了一个点外，是完美的线性关联；右下角除了一个点外，没有关联。

再次强调，我们不希望数据点对精确地落在一条直线上，但我们期望点的分布能被直线合理描述，并且我们期望y i与估计值θ ^ 0 + θ ^ 1 x i之间的偏差大致对称分布在直线周围，并且没有明显的模式。

线性模型是在第十二章介绍的，我们在那里使用了由环境保护局操作的高质量空气监测器与邻近的廉价空气质量监测器之间的关系来校准廉价监测器，以进行更准确的预测。我们重新审视那个例子，以使简单的线性模型概念更加具体。

示例：空气质量的简单线性模型

请回想一下第十二章，我们的目标是利用美国政府操作的精确空气质量系统（AQS）传感器的空气质量测量来预测由 PurpleAir（PA）传感器进行的测量。数据值对来自同一天测量的邻近仪器，测量空气中直径小于 2.5mm 颗粒物的平均每日浓度。（测量单位是每立方升空气中的 24 小时内颗粒物的平均计数。）在本节中，我们关注乔治亚州一个位置的空气质量测量。这些是我们在第十二章案例研究中检验的数据子集。这些测量是 2019 年 8 月至 2019 年 11 月中旬的日均值：

	日期	id	区域	pm25aqs	pm25pa
5258	2019-08-02	GA1	东南	8.65	16.19
5259	2019-08-03	GA1	东南	7.70	13.59
5260	2019-08-04	GA1	东南	6.30	10.30
...	...	...	...	...	...
5439	2019-10-18	GA1	东南	6.30	12.94
5440	2019-10-21	GA1	东南	7.50	13.62
5441	2019-10-30	GA1	东南	5.20	14.55

184 rows × 5 columns

特征pm25aqs包含来自 AQS 传感器的测量值，pm25pa来自 PurpleAir 监测器。由于我们有兴趣研究 AQS 测量如何预测 PurpleAir 测量，我们的散点图将 PurpleAir 读数放在 y 轴上，AQS 读数放在 x 轴上。我们还添加了趋势线：

`px``.``scatter``(``GA``,` `x``=``"``pm25aqs``"``,` `y``=``"``pm25pa``"``,` `trendline``=``'``ols``'``,`
           `trendline_color_override``=``"``darkorange``"``,`
           `labels``=``{``'``pm25aqs``'``:``'``AQS PM2.5``'``,` `'``pm25pa``'``:``'``PurpleAir PM2.5``'``}``,`
           `width``=``350``,` `height``=``250``)`

这个散点图显示了这两种仪器测量值之间的线性关系。我们要拟合的模型具有以下形式：

P A ≈ θ 0 + θ 1 A Q

其中P A表示 PurpleAir 的平均日测量值，A Q表示其伙伴 AQS 的测量值。

由于pandas.Series对象具有计算标准偏差（SDs）和相关系数的内置方法，因此我们可以快速定义计算最佳拟合线的函数：

`def` `theta_1``(``x``,` `y``)``:`
    `r` `=` `x``.``corr``(``y``)`
    `return` `r` `*` `y``.``std``(``)` `/` `x``.``std``(``)`

`def` `theta_0``(``x``,` `y``)``:`
    `return` `y``.``mean``(``)` `-` `theta_1``(``x``,` `y``)` `*` `x``.``mean``(``)`

现在我们可以通过计算这些数据的θ ^ 0和θ ^ 1来拟合模型：

`t1` `=` `theta_1``(``GA``[``'``pm25aqs``'``]``,` `GA``[``'``pm25pa``'``]``)`
`t0` `=` `theta_0``(``GA``[``'``pm25aqs``'``]``,` `GA``[``'``pm25pa``'``]``)`

Model: -3.36 + 2.10AQ

这个模型与散点图中显示的趋势线相匹配。这并非偶然。在调用scatter()时，trendline的参数值为"ols"，表示普通最小二乘法，即通过最小化平方误差来拟合线性模型的另一个名称。

让我们检查一下误差。首先，我们找出了给定 AQS 测量值的 PA 测量值的预测值，然后计算误差—实际 PA 测量值与预测值之间的差异：

`prediction` `=` `t0` `+` `t1` `*` `GA``[``"``pm25aqs``"``]`
`error` `=` `GA``[``"``pm25pa``"``]` `-` `prediction`
`fit` `=` `pd``.``DataFrame``(``dict``(``prediction``=``prediction``,` `error``=``error``)``)`

让我们将这些误差绘制成预测值的图：

`fig` `=` `px``.``scatter``(``fit``,` `y``=``'``error``'``,` `x``=``'``prediction``'``,`
                 `labels``=``{``"``prediction``"``:` `"``Prediction``"``,`
                         `"``error``"``:` `"``Error``"``}``,`
                 `width``=``350``,` `height``=``250``)`

`fig``.``add_hline``(``0``,` `line_width``=``2``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`
`fig``.``update_yaxes``(``range``=``[``-``12``,` `12``]``)`

误差为 0 意味着实际测量值落在拟合线上；我们也称这条线为最小二乘线或回归线。正值意味着它在线上方，负值意味着它在线下方。你可能想知道这个模型有多好，以及它对我们的数据说了什么。我们接下来考虑这些话题。

解释线性模型

原始的成对测量散点图显示，PurpleAir 记录往往比更准确的 AQS 测量值要高得多。事实上，我们简单线性模型的方程的斜率约为 2.1。我们解释斜率意味着 AQS 监视器测得的 1 ppm 的变化平均对应于 PA 测量的 2 ppm 的变化。因此，如果一天 AQS 传感器测量 10 ppm，第二天它高出 5 ppm，即 15 ppm，那么我们对于下一天的 PA 测量的预测将增加2 × 5 = 10 ppm。

任何 PurpleAir 读数的变化都不是由 AQS 读数的变化引起的。相反，它们都反映了空气质量，而我们的模型捕捉了这两个设备之间的关系。通常情况下，术语预测被认为是因果关系，但在这里并非如此。相反，预测只是指我们对 PA 和 AQS 测量之间线性关联的使用。

至于模型中的截距，我们可能期望它为 0，因为当空气中没有颗粒物时，我们认为两个仪器都应该测量 0 ppm。但对于 AQS 为 0 的情况，模型预测 PurpleAir 为− 3.4 ppm，这是没有意义的。空气中不可能有负的颗粒物。这突显了在超出测量范围时使用模型的问题。我们观察到 AQS 记录在 3 到 18 ppm 之间，并且在这个范围内，模型拟合良好。虽然在理论上线应该有一个截距为 0，但在实际中这样的模型却不适用，预测往往会差得多。

著名统计学家 George Box 曾经说过：“所有模型都是错误的，但有些是有用的。” 在这里，尽管线的截距不通过 0，但简单线性模型在预测 PurpleAir 传感器的空气质量测量方面是有用的。事实上，我们两个特征之间的相关性非常高：

`GA``[``[``'``pm25aqs``'``,` `'``pm25pa``'``]``]``.``corr``(``)`

	pm25aqs	pm25pa
pm25aqs	1.00	0.92
pm25pa	0.92	1.00

除了查看相关系数之外，还有其他评估线性模型质量的方法。

评估拟合

早期的误差图针对拟合值给出了拟合质量的视觉评估。（这种图称为残差图，因为误差有时被称为残差。）一个良好的拟合应该显示一群点围绕着 0 的水平线，没有明显的模式。当出现模式时，我们通常可以得出简单线性模型并没有完全捕捉到信号的结论。我们之前看到残差图中没有明显的模式。

另一种有用的残差图类型是残差与不在模型中的特征的图。如果我们看到模式，那么我们可能希望在模型中加入这个特征，除了已经在模型中的特征之外。此外，当数据具有时间组成部分时，我们希望检查残差随时间的模式。对于这些特定的数据，由于测量是在四个月内的每日平均值，我们将错误绘制为测量记录日期：

看起来在八月底和九月底附近有几天数据远低于预期。回顾原始散点图（以及第一个残差图），我们可以看到两个小的水平点簇在主要点云下方。我们刚刚制作的图表表明，我们应该检查原始数据以及关于设备的任何可用信息，以确定这些天是否正常运行。

残差图还可以让我们大致了解模型在预测中的准确性。大多数误差在线路的 ± 6 ppm 之间。我们发现误差的标准偏差约为 2.8 ppm：

`error``.``std``(``)`

2.796095864304746

相比之下，PurpleAir 测量的标准偏差要大得多：

`GA``[``'``pm25pa``'``]``.``std``(``)`

6.947418231019876

如果我们发现监测器在八月底和九月份的某些日子不工作，并因此将其排除在数据集之外，可能会进一步减少模型误差。无论如何，在空气非常清洁的情况下，误差相对较大，但在绝对值上并不重要。我们通常更关心空气污染的情况，此时 2.8 ppm 的误差似乎是合理的。

让我们回到如何找到这条线的过程，即模型拟合的过程。在接下来的部分，我们通过最小化均方误差来推导截距和斜率。

拟合简单线性模型

我们在本章早些时候提到，当我们最小化数据的平均损失时：

1 n ∑ i [ y i − ( θ 0 + θ 1 x i ) ] 2

最佳拟合线具有截距和斜率：

θ ^ 0 = y ¯ − θ ^ 1 x ¯ θ ^ 1 = r ( x , y ) S D ( y ) S D ( x )

在本节中，我们使用微积分来推导这些结果。

对于简单线性模型，均方误差是两个模型参数的函数，即截距和斜率。这意味着如果我们使用微积分来找到最小化的参数值，我们需要找到均方误差对 θ 0 和 θ 1 的偏导数。我们也可以通过其他技术找到这些最小值：

梯度下降

当损失函数更复杂且找到近似解更快时，我们可以使用数值优化技术，如梯度下降（参见 第二十章）。

二次公式

由于平均损失是关于 θ 0 和 θ 1 的二次函数，我们可以使用二次公式（以及一些代数）来求解最小化参数值。

几何论证

在本章后面，我们使用最小二乘法的几何解释来拟合多元线性模型。这种方法与毕达哥拉斯定理相关，并具有几个直观的优点。

我们选择使用微积分来优化简单线性模型，因为这是快速且直接的方法。首先，我们对平方误差的偏导数进行计算（我们可以忽略 MSE 中的 e 1 / n ，因为它不影响最小值的位置）：

∂ ∂ θ 0 ∑ i [ y i − ( θ 0 + θ 1 x i ) ] 2 = ∑ i 2 ( y i − θ 0 − θ 1 x i ) ( − 1 ) ∂ ∂ θ 1 ∑ i [ y i − ( θ 0 + θ 1 x i ) ] 2 , = ∑ i 2 ( y i − θ 0 − θ 1 x i ) ( − x i )

然后我们将偏导数设为零，并通过将方程两边乘以 − 1 / 2 进行简化，得到：

0 = ∑ i ( y i − θ ^ 0 − θ ^ 1 x i ) 0 = ∑ i ( y i − θ ^ 0 − θ ^ 1 x i ) x i

这些方程式称为正规方程式。在第一个方程中，我们看到θ ^ 0可以表示为θ ^ 1的函数：

θ ^ 0 = y ¯ − θ ^ 1 x ¯

将这个值代入第二个方程中给出我们：

0 = ∑ i ( y i − y ¯ + θ ^ 1 x ¯ − θ ^ 1 x i ) x i = ∑ i [ ( y i − y ¯ ) − θ ^ 1 ( x i − x ¯ ) ] x i θ ^ 1 = ∑ i ( y i − y ¯ ) x i ∑ i ( x i − x ¯ ) x i

经过一些代数运算，我们可以用我们熟悉的量来表示θ ^ 1：

θ ^ 1 = r ( x , y ) S D ( y ) S D ( x )

正如本章前面所示，这个表示法表明拟合线上的点在x处可以写成如下形式：

θ ^ 0 + θ ^ 1 x = y ¯ + r ( x , y ) S D ( y ) ( x − x ¯ ) S D ( x )

我们已经推导出了在前一节中使用的最小二乘线方程。在那里，我们使用了pandas内置方法来计算S D ( x )，S D ( y )和r ( x , y )，以便轻松计算这条线的方程。然而，在实践中，我们建议使用scikit-learn提供的功能来进行模型拟合：

`from` `sklearn``.``linear_model` `import` `LinearRegression` 

`y` `=` `GA``[``'``pm25pa``'``]`
`x` `=` `GA``[``[``'``pm25aqs``'``]``]`
`reg` `=` `LinearRegression``(``)``.``fit``(``x``,` `y``)`

我们的拟合模型是：

Model: PA estimate = -3.36 + 2.10AQS

注意，我们将y作为数组和x作为数据框传递给LinearRegression。当我们在模型中引入多个解释特征时，很快就会看到原因。

LinearRegression方法提供了稳定的数值算法来通过最小二乘法拟合线性模型。当拟合多个变量时，这一点尤为重要，接下来我们将介绍。

多元线性模型

到目前为止，在本章中，我们使用单个输入变量预测结果变量。现在我们介绍使用多个特征的多元线性模型来预测（或描述或解释）结果。具有多个解释特征可以改善模型对数据的拟合并提高预测能力。

我们从一个简单的线性模型推广到包括第二个解释变量的模型，称为v。这个模型在x和v上都是线性的，这意味着对于x和v的一对数值，我们可以用线性组合来描述、解释或预测y：

y ≈ θ 0 + θ 1 x + θ 2 v

注意，对于特定的v值，比如v ⋆，我们可以将上述方程表示为：

y ≈ ( θ 0 + θ 2 v ⋆ ) + θ 1 x

换句话说，当我们将v固定在v ⋆时，x和y之间有一个简单的线性关系，斜率为θ 1，截距为θ 0 + θ 2 v ⋆。对于另一个v的值，比如v †，我们同样有x和y之间的简单线性关系。x的斜率保持不变，唯一的变化是截距，现在是θ 0 + θ 2 v †。

使用多元线性回归时，我们需要记住在模型中的其他变量存在的情况下解释θ 1对x的系数。在保持模型中其他变量（在本例中仅为v）的值不变的情况下，x增加 1 个单位平均对应于y的θ 1的变化。一种可视化这种多元线性关系的方法是创建散点图的多面板，其中每个图中v的值大致相同。我们接下来为空气质量测量制作这样的散点图，并提供其他可视化和统计学示例以检验拟合多元线性模型时的情况。

研究空气质量监测仪的科学家们（参见第十二章）寻找一个包含天气因素的改进模型。他们检查的一种天气变量是相对湿度的每日测量值。让我们考虑一个双变量线性模型，以解释基于 AQS 传感器测量和相对湿度的 PurpleAir 测量。该模型具有以下形式：

P A ≈ θ 0 + θ 1 A Q + θ 2 R H

其中P A，A Q和R H分别指代变量：PurpleAir 平均每日测量、AQS 测量和相对湿度。

作为第一步，我们制作一个多面板图来比较固定湿度值下两种空气质量测量之间的关系。为此，我们将相对湿度转换为一个分类变量，使每个面板由湿度相似的观测组成。

`rh_cat` `=` `pd``.``cut``(``GA``[``'``rh``'``]``,` `bins``=``[``43``,``50``,``55``,``60``,``78``]``,` 
                `labels``=``[``'``<50``'``,``'``50-55``'``,``'``55-60``'``,``'``>60``'``]``)`

然后我们使用这个定性特征将数据划分为一个二乘二的散点图面板：

`fig` `=` `px``.``scatter``(``GA``,` `x``=``'``pm25aqs``'``,` `y``=``'``pm25pa``'``,` 
                 `facet_col``=``rh_cat``,` `facet_col_wrap``=``2``,`
                 `facet_row_spacing``=``0.15``,`
                 `labels``=``{``'``pm25aqs``'``:``'``AQS PM2.5``'``,` `'``pm25pa``'``:``'``PurpleAir PM2.5``'``}``,`
                 `width``=``550``,` `height``=``350``)`

`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`
`fig``.``show``(``)`

这四个图表显示了两种空气质量测量来源之间的线性关系。斜率看起来相似，这意味着多重线性模型可能非常合适。从这些图表中很难看出相对湿度是否对截距有显著影响。

我们还想检查三个特征之间的成对散点图。当两个解释性特征高度相关时，它们在模型中的系数可能不稳定。虽然三个或更多特征之间的线性关系在成对图中可能不明显，但检查这些图表仍然是一个好主意：

`fig` `=` `px``.``scatter_matrix``(`
    `GA``[``[``'``pm25pa``'``,` `'``pm25aqs``'``,` `'``rh``'``]``]``,`
    `labels``=``{``'``pm25aqs``'``:``'``AQS``'``,` `'``pm25pa``'``:``'``PurpleAir``'``,` `'``rh``'``:``'``Humidity``'``}``,`
    `width``=``550``,` `height``=``400``)`

`fig``.``update_traces``(``diagonal_visible``=``False``)`

湿度与空气质量之间的关系似乎并不特别强。我们应该检查的另一个成对测量是特征之间的相关性：

	pm25pa	pm25aqs	rh
pm25pa	1.00	0.95	-0.06
pm25aqs	0.95	1.00	-0.24
rh	-0.06	-0.24	1.00

一个小惊喜是，相对湿度与 AQS 测量的空气质量具有轻微的负相关。这表明湿度可能对模型有帮助。

在下一节中，我们将推导适合的方程。但现在，我们使用LinearRegression的功能来拟合模型。与之前不同的唯一变化是我们为解释变量提供了两列（这就是为什么x输入是一个数据框）：

`from` `sklearn``.``linear_model` `import` `LinearRegression`

`y` `=` `GA``[``'``pm25pa``'``]`
`X2` `=` `GA``[``[``'``pm25aqs``'``,` `'``rh``'``]``]`

`model2` `=` `LinearRegression``(``)``.``fit``(``X2``,` `y``)`

适合的多重线性模型，包括系数单位，是：

PA estimate = -15.8 ppm + 2.25 ppm/ppm x AQS +  0.21 ppm/percent x RH

模型中湿度的系数调整空气质量预测每百分点相对湿度 0.21 ppm。请注意，AQS 的系数与我们之前拟合的简单线性模型不同。这是因为系数反映了来自相对湿度的额外信息。

最后，为了检查拟合质量，我们制作了预测值和误差的残差图。这一次，我们使用LinearRegression来计算我们的预测：

`predicted_2var` `=` `model2``.``predict``(``X2``)`
`error_2var` `=` `y` `-` `predicted_2var`

`fig` `=` `px``.``scatter``(``y` `=` `error_2var``,` `x``=``predicted_2var``,`
                 `labels``=``{``"``y``"``:` `"``Error``"``,` `"``x``"``:` `"``Predicted PurpleAir measurement``"``}``,`
                 `width``=``350``,` `height``=``250``)`

`fig``.``update_yaxes``(``range``=``[``-``12``,` `12``]``)`
`fig``.``add_hline``(``0``,` `line_width``=``3``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`

`fig``.``show``(``)`

残差图表没有明显的模式，这表明模型拟合得相当好。还要注意，误差几乎都落在–4 和+4 ppm 之间，比简单线性模型的范围小。我们发现残差的标准偏差要小得多：

`error_2var``.``std``(``)`

1.8211427707294048

残差标准偏差从单变量模型的 2.8 ppm 降低到了 1.8 ppm，这是一个很好的尺寸缩减。

当我们有多个解释变量时，相关系数无法捕捉线性关联模型的强度。相反，我们调整 MSE 以了解模型的拟合程度。在下一节中，我们描述如何拟合多重线性模型并使用 MSE 来评估拟合。

拟合多重线性模型

在前一节中，我们考虑了两个解释变量的情况；其中一个我们称为x，另一个为v。现在我们希望将这种方法推广到p个解释变量。选择不同字母来表示变量的想法很快失效了。相反，我们使用一种更正式和通用的方法，将多个预测变量表示为一个矩阵，如图 15-3 所示。我们称X为设计矩阵。注意，X的形状为n × ( p + 1 )。X的每一列代表一个特征，每一行代表一个观察值。也就是说，x i , j是在观察值i上针对特征j的测量值。

图 15-3。在这个设计矩阵X中，每一行代表一个观察/记录，每一列代表一个特征/变量

注意

一个技术细节：设计矩阵被定义为数学矩阵，而不是数据框，因此您可能注意到矩阵不包括数据框具有的列或行标签。

也就是说，我们通常不必担心将数据框转换为矩阵，因为大多数用于建模的 Python 库将数字数据框视为矩阵。

对于给定的观察值，比如X中的第二行，我们通过线性组合近似得到结果y 2：

y 2 ≈ θ 0 + θ 1 x 2 , 1 + … + θ p x 2 , p

用矩阵表示线性近似更方便。为此，我们将模型参数写成一个p + 1列向量θ：

θ = [ θ 0 θ 1 ⋮ θ p ]

将这些符号定义放在一起，我们可以使用矩阵乘法为整个数据集编写预测向量：

X θ

如果我们检查X和θ的维度，我们可以确认Xθ是一个n维列向量。因此，使用这种线性预测的误差可以表示为向量：

e = y − X θ

其中结果变量也表示为列向量：

y = [ y 1 y 2 ⋮ y n ]

这种多重线性模型的矩阵表示可以帮助我们找到使均方误差最小化的模型。我们的目标是找到模型参数 ( θ 0 , θ 1 , … , θ p ) ，使均方误差最小化：

1 n ∑ i [ y i − ( θ 0 + θ 1 x i , 1 + ⋯ + θ p x i , p ) ] 2 = 1 n ‖ y − X θ ‖ 2

在这里，我们使用记号 ‖ v ‖ 2 表示向量 v 的长度的平方和的简写形式：‖ v ‖ 2 = ∑ i v i 2 。平方根 ‖ v ‖ 2 对应于向量 v 的长度，也称为向量 v 的 ℓ 2 范数。因此，最小化均方误差等同于找到最短的误差向量。

我们可以像简单线性模型那样使用微积分来拟合我们的模型。然而，这种方法变得笨重，我们改用更直观的几何论证，这更容易导致设计矩阵、误差和预测值的有用属性。

我们的目标是找到参数向量，我们称之为θ ^，使我们的平均平方损失最小化——我们希望使‖ y − X θ ‖ 2 在给定的X和y 下尽可能小。关键洞察力在于，我们可以以几何方式重新表述这个目标。由于模型预测和真实结果都是向量，我们可以将它们视为向量空间中的向量。当我们改变模型参数θ时，模型会进行不同的预测，但任何预测必须是X的列向量的线性组合；也就是说，预测必须在所谓的span ( X ) 中。这个概念在图 15-4 中有所体现，阴影区域代表可能的线性模型。请注意，y并没有完全包含在span ( X )中；这通常是情况。

图 15-4。在这个简化的图示中，所有可能的模型预测向量span ( X )被描绘为三维空间中的一个平面，而观测到的y作为一个向量。

尽管平方损失不能完全为零，因为y不在span ( X )中，我们可以找到一个尽可能接近y但仍在span ( X )中的向量。这个向量被称为y ^。

误差是向量 e = y − y ^ 。它的长度 ‖ e ‖ 表示真实结果与我们模型预测之间的距离。从视觉上看，当它与 span ( X ) 垂直 时，e 的大小最小，如 图 15-5 所示。关于此事实的证明被省略，我们依赖于图表来说服您。

图 15-5. 当预测值 y ^ 在 span ( X ) 垂直于 y 时，均方误差达到最小值。

最小误差 e 必须垂直于 y ^，这使我们能够推导出 θ ^ 的公式如下：

X θ ^ + e = y ( the definition of  y , y ^ , e ) X ⊤ X θ ^ + X ⊤ e = X ⊤ y ( left-multiply by  X ⊤ ) X ⊤ X θ ^ = X ⊤ y ( e ⊥ span ( X ) ) θ ^ = ( X ⊤ X ) − 1 X ⊤ y ( left-multiply by  ( X ⊤ X ) − 1 )

这种推导多元线性模型中θ ^的一般方法也给了我们简单线性模型中θ ^ 0和θ ^ 1。如果我们将X设置为包含截距列和一个特征列的两列矩阵，这个公式用于最小二乘拟合的简单线性模型的截距和斜率。实际上，如果X仅是1列的单列，那么我们可以使用这个公式表明θ ^只是y的均值。这与我们在第四章中介绍的常数模型很好地联系在一起。

注意

虽然我们可以编写一个简单的函数来根据公式推导θ ^

θ ^ = ( X ⊤ X ) − 1 X ⊤ y

我们建议使用优化调整方法来计算θ ^，这些方法由scikit-learn和statsmodels库提供。它们处理设计矩阵稀疏、高度共线性和不可逆的情况。

这个θ ^的解（以及图像）揭示了拟合系数和预测的一些有用性质：

• 残差e与预测值y ^正交。

• 如果模型有截距项，则残差的平均值为 0。

• 残差的方差就是均方误差。

这些属性解释了为什么我们要检查残差与预测值的图表。当我们拟合多元线性模型时，我们还会将残差与我们考虑添加到模型的变量绘制在一起。如果它们显示出线性模式，那么我们会考虑将它们添加到模型中。

除了检查错误的标准差之外，多元线性模型的均方误差与常数模型的均方误差比值可以衡量模型的拟合度。这被称为多元R 2，其定义如下：

R 2 = 1 − ‖ y − X θ ^ ‖ 2 ‖ y − y ¯ ‖ 2

随着模型越来越贴近数据，多个R 2接近于 1。这可能看起来是件好事，但这种方法可能存在问题，因为R 2即使在我们为模型添加无意义的特征时也会继续增长，只要这些特征扩展了span(X)。为了考虑模型的大小，我们通常通过模型中拟合系数的数量调整R 2的分子和分母。也就是说，我们通过1/[n−(p+1)]来标准化分子，并通过1/(n−1)来标准化分母。在选择模型的更好方法方面，详见第十六章。

接下来，我们考虑一个社会科学的例子，在这个例子中，我们有许多可用于建模的变量。

示例：什么是机会之地？

美国被称为“机会之地”，因为人们相信即使资源匮乏的人也可以在美国变得富有，经济学家称这种观念为“经济流动性”。在一项研究中，经济学家拉杰·切蒂及其同事对美国的经济流动性进行了大规模数据分析。他的基本问题是美国是否是一个机会之地。为了回答这个相对模糊的问题，切蒂需要一种衡量经济流动性的方法。

Chetty 可以访问 1980 年至 1982 年间出生于美国的每个人的 2011-2012 年联邦所得税记录，以及他们父母在他们出生年份的税务记录。他通过找到列出他们为家庭成员的父母的 1980-1982 年税务记录来将 30 岁的人与他们的父母配对。总共，他的数据集约有 1000 万人。为了衡量经济流动性，Chetty 将出生在特定地理区域、父母收入位于 1980-1982 年的第 25 个收入百分位的人群分组。然后，他找到了该组 2011 年的平均收入百分位数。Chetty 将这个平均值称为绝对向上流动（AUM）。如果一个地区的 AUM 为 25，那么出生在第 25 百分位的人通常会保持在第 25 百分位，即他们留在了父母出生时的位置。高 AUM 值意味着该地区具有更多的向上流动性。在这些地区出生在第 25 个收入百分位的人通常会进入比他们父母更高的收入阶层。作为参考，美国的平均 AUM 在撰写本文时约为 41。Chetty 计算了称为通勤区（CZs）的地区的 AUM，这些地区大致与县级相同的规模。

虽然原始数据的粒度是在个体级别，Chetty 分析的数据粒度是在通勤区域级别。由于隐私法规限制，收入记录不能公开，但通勤区域的 AUM 可以提供。然而，即使有了通勤区域的粒度，也并非所有通勤区域都包含在数据集中，因为在 40 个特征的数据中，可能会识别出小型通勤区域的个体。这一限制指向潜在的覆盖偏倚。测量偏倚是另一个潜在问题。例如，出生在第 25 收入百分位的儿童，如果成为极其富有的人，可能不会申报所得税。

我们还指出使用区域平均数据而不是个体测量数据的局限性。在聚合水平上，特征之间的关系通常比在个体水平上更高度相关。这种现象称为生态回归，需要谨慎解释从聚合数据中得出的发现。

Chetty 怀疑美国某些地方的经济流动性较高。他的分析证实了这一点。他发现一些城市，如加利福尼亚州圣何塞、华盛顿特区和西雅图，比其他城市如北卡罗来纳州夏洛特、密尔沃基和亚特兰大有更高的流动性。这意味着，例如，在圣何塞，人们从低收入阶层向高收入阶层的转移速度比在夏洛特要快。Chetty 使用线性模型发现社会和经济因素如隔离、收入不平等和当地学校系统与经济流动性相关。

在这个分析中，我们的结果变量是通勤区域的 AUM，因为我们有兴趣找出与 AUM 相关的特征。Chetty 的数据中可能有许多这样的特征，但我们首先调查了一个特别的特征：通勤区域内通勤时间在 15 分钟或更短的人口比例。

使用通勤时间解释向上流动性

我们开始通过将数据加载到名为cz_df的数据框中进行调查：

	aum	travel_lt15	gini	rel_tot	...	taxrate	worked_14	foreign	地区
0	38.39	0.33	0.47	0.51	...	0.02	3.75e-03	1.18e-02	南部
1	37.78	0.28	0.43	0.54	...	0.02	4.78e-03	2.31e-02	南部
2	39.05	0.36	0.44	0.67	...	0.01	2.89e-03	7.08e-03	南部
...	...	...	...	...	...	...	...	...	...
702	44.12	0.42	0.42	0.29	...	0.02	4.82e-03	9.85e-02	西部
703	41.41	0.49	0.41	0.26	...	0.01	4.39e-03	4.33e-02	西部
704	43.20	0.24	0.42	0.32	...	0.02	3.67e-03	1.13e-01	西部

705 rows × 9 columns

每一行代表一个通勤区。列 aum 是 1980–1982 年出生并且父母收入处于第 25 百分位数的人群的平均 AUM。数据框中有许多列，但现在我们专注于通勤区内通勤时间不超过 15 分钟的人群比例，即 travel_lt15。我们将 AUM 与这个比例进行绘图，以探索这两个变量之间的关系：

`px``.``scatter``(``cz_df``,` `x``=``'``travel_lt15``'``,` `y``=``'``aum``'``,` `width``=``350``,` `height``=``250``,`
           `labels``=``{``'``travel_lt15``'``:``'``Commute time under 15 min``'``,` 
                   `'``aum``'``:``'``Upward mobility``'``}``)`

散点图显示 AUM 与通勤时间之间存在大致的线性关系。事实上，我们发现它们之间的相关性非常强：

`cz_df``[``[``'``aum``'``,` `'``travel_lt15``'``]``]``.``corr``(``)`

	aum	travel_lt15
aum	1.00	0.68
travel_lt15	0.68	1.00

让我们用一个简单的线性模型来解释 AUM 与通勤时间的关系：

`from` `sklearn``.``linear_model` `import` `LinearRegression`

`y` `=` `cz_df``[``'``aum``'``]`
`X` `=` `cz_df``[``[``'``travel_lt15``'``]``]`

`model_ct` `=` `LinearRegression``(``)``.``fit``(``X``,` `y``)`

MSE 最小化得到的系数为：

Intercept: 31.3
    Slope: 28.7

有趣的是，通勤区的向上流动增加与通勤时间较短的人群比例增加相关联。

我们可以将 AUM 测量的标准差与残差的标准差进行比较。这种比较让我们了解模型在解释 AUM 方面的实用性：

`prediction` `=` `model_ct``.``predict``(``X``)`
`error` `=` `y` `-` `prediction`

`print``(``f``"``SD(errors):` `{``np``.``std``(``error``)``:``.2f``}``"``)`
`print``(``f``"` `SD(AUM):` `{``np``.``std``(``cz_df``[``'``aum``'``]``)``:``.2f``}``"``)`

SD(errors): 4.14
   SD(AUM): 5.61

围绕回归线的误差大小比常数模型减少了约 25%。

接下来，我们检查残差来判断拟合的不合适，因为在残差图中更容易看出拟合存在的潜在问题：

`fig` `=` `px``.``scatter``(``x``=``prediction``,` `y``=``error``,`
                 `labels``=``dict``(``x``=``'``Prediction for AUM``'``,` `y``=``'``Error``'``)``,`
                 `width``=``350``,` `height``=``250``)`

`fig``.``add_hline``(``0``,` `line_width``=``2``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`
`fig``.``update_yaxes``(``range``=``[``-``20``,` `15``]``)`

`fig``.``show``(``)`

看起来随着 AUM 的增加，误差也在增加。我们可以尝试对响应变量进行转换，或者拟合一个在通勤时间分数上是二次的模型。在下一节中我们将考虑变换和多项式。首先，我们看看包含额外变量是否能更准确地预测 AUM。

使用多个变量来关联向上流动

在他的原始分析中，Chetty 创建了几个与隔离、收入和 K–12 教育等因素相关的高级特征。我们考虑 Chetty 的七个预测因子，旨在构建一个更具信息性的模型来解释 AUM。这些在 Table 15-1 中描述。

Table 15-1. 解释 AUM 建模的潜在原因

列名	描述
travel_lt15	上班通勤时间不超过 15 分钟的人群比例。
gini	基尼系数，财富不平等的度量。取值介于 0 到 1 之间，数值较小表示财富分配较均匀，较大表示不平等程度更大。
rel_tot	自报宗教信仰的人群比例。
single_mom	单身母亲的子女比例。
taxrate	地方税率。
worked_14	14 到 16 岁工作的人群比例。
foreign	出生于美国以外的人群比例。

让我们首先检查 AUM 与解释性特征以及解释性特征之间的相关性：

	aum	travel_lt15	gini	rel_tot	single_mom	taxrate	worked_14	foreign
aum	1.00	0.68	-0.60	0.52	-0.77	0.35	0.65	-0.03
travel_lt15	0.68	1.00	-0.56	0.40	-0.42	0.34	0.60	-0.19
gini	-0.60	-0.56	1.00	-0.29	0.57	-0.15	-0.58	0.31
rel_tot	0.52	0.40	-0.29	1.00	-0.31	0.08	0.28	-0.11
single_mom	-0.77	-0.42	0.57	-0.31	1.00	-0.26	-0.60	-0.04
taxrate	0.35	0.34	-0.15	0.08	-0.26	1.00	0.35	0.26
worked_14	0.65	0.60	-0.58	0.28	-0.60	0.35	1.00	-0.15
foreign	-0.03	-0.19	0.31	-0.11	-0.04	0.26	-0.15	1.00

我们看到，在通勤区中单身母亲的比例与 AUM 有最强的相关性，这意味着它也是解释 AUM 的最佳特征。此外，我们看到几个解释变量彼此之间高度相关；基尼系数与工作的青少年比例、单身母亲比例以及 15 分钟以下通勤比例高度相关。由于这些高度相关的特征，我们在解释系数时需要谨慎，因为几种不同的模型可能同样能够用协变量来解释 AUM。

注意

我们在本章前面介绍的向量几何视角可以帮助我们理解这个问题。回顾一下，一个特征对应于 n 维空间中的一个列向量，如 x 。对于两个高度相关的特征 x 1 和 x 2 ，这些向量几乎是对齐的。因此，响应向量 y 在这些向量中的一个上的投影几乎与在另一个上的投影相同。当几个特征彼此相关时，情况变得更加混乱。

首先，我们可以考虑所有可能的两特征模型，看看哪一个具有最小的预测误差。Chetty 导出了 40 个潜在的变量作为预测变量，这将使我们检查 ( 40 × 39 ) / 2 = 780 个模型。拟合模型时，所有成对、三元组等变量很快就会失控。这可能导致找到伪相关性（见第十七章）。

在这里，我们保持事情稍微简单，只研究包含通勤时间和单身母亲特征的两变量模型。之后，我们查看包含数据框架中所有七个数值解释特征的模型：

`X2` `=` `cz_df``[``[``'``travel_lt15``'``,` `'``single_mom``'``]``]`
`y` `=` `cz_df``[``'``aum``'``]`

`model_ct_sm` `=` `LinearRegression``(``)``.``fit``(``X2``,` `y``)`

Intercept: 49.0
Fraction with under 15 minute commute coefficient: 18.10
Fraction of single moms coefficient: 18.10

请注意，旅行时间的系数与简单线性模型中这个变量的系数相比有很大不同。这是因为我们模型中的两个特征高度相关。

接下来我们比较这两个拟合的误差：

`prediction_ct_sm` `=` `model_ct_sm``.``predict``(``X2``)`
`error_ct_sm` `=` `y` `-` `prediction_ct_sm`

 SD(errors in model 1): 4.14
 SD(errors in model 2): 2.85

残差的标准偏差进一步减少了 30%。增加模型复杂性以添加第二个变量似乎是值得的。

让我们再次直观地检查残差。我们使用与前一个单变量模型相同的 y 轴刻度，以便与其残差图进行比较：

`fig` `=` `px``.``scatter``(``x``=``prediction_ct_sm``,` `y``=``error_ct_sm``,`
           `labels``=``dict``(``x``=``'``Two-variable prediction for AUM``'``,` `y``=``'``Error``'``)``,`
           `width``=``350``,` `height``=``250``)`

`fig``.``add_hline``(``0``,` `line_width``=``2``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`
`fig``.``update_yaxes``(``range``=``[``-``20``,` `15``]``)`

`fig``.``show``(``)`

对更高 AUM 的误差的较大变异性更为明显。这意味着估计值y ^不受影响，但其准确性取决于 AUM。可以通过加权回归来解决这个问题。

注意

再次强调，不同背景的数据科学家使用不同的术语来指代相同的概念。例如，将设计矩阵X中的每一行称为一个观察值，每一列称为一个变量的术语在统计背景的人群中更为常见。其他人则称设计矩阵的每一列代表一个特征，或者每一行代表一条记录。此外，我们称拟合和解释模型的整个过程为建模，而其他人则称其为机器学习。

现在让我们拟合一个多元线性模型，使用所有七个变量来解释上升的流动性。在拟合模型之后，我们再次使用与前两个残差图相同的 y 轴刻度绘制误差：

`X7` `=` `cz_df``[``predictors``]`
`model_7var` `=` `LinearRegression``(``)``.``fit``(``X7``,` `y``)`

`prediction_7var` `=` `model_7var``.``predict``(``X7``)`
`error_7var` `=` `y` `-` `prediction_7var`

fig = px.scatter(
    x=prediction_7var, y=error_7var,
    labels=dict(x='Seven-variable prediction for AUM', y='Error'),
    width=350, height=250)

`fig``.``add_hline``(``0``,` `line_width``=``2``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`
`fig``.``update_yaxes``(``range``=``[``-``20``,` `15``]``)`

`fig``.``show``(``)`

具有七个特征的模型似乎并没有比具有两个变量的模型好多少。事实上，残差的标准偏差仅减少了 8%：

`error_7var``.``std``(``)` 

2.588739233574256

我们可以比较这三个模型的多元R 2：

R² for 7-variable model: 0.79
R² for 2-variable model: 0.74
R² for 1-variable model: 0.46

对于我们来说，模型中特征数量的调整并没有太大差异，因为我们有超过 700 个观测值。现在我们已经确认了之前的发现，即使用两个变量大大改善了模型的解释能力，而七个变量模型几乎没有比两个变量模型有所改善。这种小的增益可能不值得模型的复杂性增加。

到目前为止，我们的模型仅使用了数值预测变量。但是类别数据在模型拟合中通常也很有用。此外，在第十章中，我们对变量进行了转换，并从变量的组合中创建了新的变量。接下来我们将讨论如何将这些变量纳入线性模型。

数值测量的特征工程

本章迄今为止我们拟合的所有模型都使用了最初在数据框中提供的数值特征。在本节中，我们将查看由数值特征变换创建的变量。将变量转换为建模使用的形式称为特征工程。

我们在第九章和 10 章中引入了特征工程。在那里，我们对特征进行了转换，使它们具有对称分布。变换可以捕捉数据中更多种类的模式，并导致更好和更准确的模型。

让我们回到我们在第十章中作为示例使用的数据集：旧金山湾区的房屋销售价格。我们将数据限制在 2006 年售出的房屋，当时房价相对稳定，因此我们不需要考虑价格趋势。

我们希望建立销售价格模型。回顾在第十章中的可视化结果，我们发现销售价格与多个特征相关，如房屋尺寸、地块尺寸、卧室数量和位置。我们对销售价格和房屋尺寸进行了对数变换以改善它们之间的关系，并且发现了关于卧室数量和城市的箱线图也显示了有趣的关系。在本节中，我们将在线性模型中包括转换后的数值特征。在下一节中，我们还将向模型添加序数特征（卧室数量）和名义特征（城市）。

首先，我们将在房屋尺寸上建立销售价格模型。相关矩阵告诉我们哪些数值解释变量（原始和转换后的）与销售价格最相关：

	price	br	lsqft	bsqft	log_price	log_bsqft	log_lsqft	ppsf	log_ppsf
price	1.00	0.45	0.59	0.79	0.94	0.74	0.62	0.49	0.47
br	0.45	1.00	0.29	0.67	0.47	0.71	0.38	-0.18	-0.21
lsqft	0.59	0.29	1.00	0.46	0.55	0.44	0.85	0.29	0.27
bsqft	0.79	0.67	0.46	1.00	0.76	0.96	0.52	-0.08	-0.10
log_price	0.94	0.47	0.55	0.76	1.00	0.78	0.62	0.51	0.52
log_bsqft	0.74	0.71	0.44	0.96	0.78	1.00	0.52	-0.11	-0.14
log_lsqft	0.62	0.38	0.85	0.52	0.62	0.52	1.00	0.29	0.27
ppsf	0.49	-0.18	0.29	-0.08	0.51	-0.11	0.29	1.00	0.96
log_ppsf	0.47	-0.21	0.27	-0.10	0.52	-0.14	0.27	0.96	1.00

销售价格与房屋尺寸最相关，称为bsqft（建筑面积）。我们制作了销售价格与房屋尺寸的散点图，以确认这种关联是线性的：

关系看起来大致是线性的，但非常大和昂贵的房屋远离分布中心，可能会对模型产生过度影响。如第十章所示，对数变换使得价格和尺寸的分布更对称（两者均为以对数 10 为底以便于将值转换为原始单位）：

理想情况下，使用变换的模型应当在数据的背景下有意义。如果我们基于对数（大小）拟合一个简单的线性模型，那么在检查系数时，我们可以考虑百分比增加。例如，x翻倍会使预测增加θlog ⁡ (2)，因为θlog ⁡ (2x)=θlog ⁡ (2)+θlog ⁡ (x)。

让我们从拟合一个通过房屋大小的对数变换解释的模型开始。但首先，我们注意到这个模型仍然被认为是一个线性模型。如果我们用y表示销售价格，x表示房屋大小，那么该模型是：

log ⁡ ( y ) = θ 0 + θ 1 log ⁡ ( x )

（请注意，在这个方程中，我们忽略了近似以使线性关系更加清晰。）这个方程可能看起来不是线性的，但如果我们将log ⁡ ( y )重命名为w，log ⁡ ( x )重命名为v，那么我们可以将这种“对数-对数”关系表达为w和v的线性模型：

w = θ 0 + θ 1 v

可以表达为转换特征的线性组合的其他模型示例是：

log ⁡ ( y ) = θ 0 + θ 1 x y = θ 0 + θ 1 x + θ 2 x 2 y = θ 0 + θ 1 x + θ 2 z + θ 3 x z

再次，如果我们将log ⁡ ( y )重命名为w，x 2重命名为u，xz重命名为t，那么我们可以将每个模型表示为这些重命名特征的线性组合。按顺序，前述模型现在是：

w = θ 0 + θ 1 x y = θ 0 + θ 1 x + θ 2 u y = θ 0 + θ 1 x + θ 2 z + θ 3 t

简言之，我们可以将包含特征的非线性变换和/或特征组合的模型视为其派生特征的线性。在实践中，当我们描述模型时，我们不会重命名转换后的特征；相反，我们使用原始特征的变换，因为在解释系数和检查残差图时保持追踪它们非常重要。

当我们提及这些模型时，我们包括对变换的提及。也就是说，当结果和解释变量都经过对数变换时，我们称之为对数-对数模型；当结果经过对数变换而解释变量没有时，我们称之为对数-线性；当解释变量包括二次幂变换时，我们描述模型具有二次多项式特征；当两个解释特征的乘积包含在模型中时，我们称之为交互项。

让我们拟合一个价格对大小的对数-对数模型：

`X1_log` `=` `sfh``[``[``'``log_bsqft``'``]``]`    
`y_log` `=` `sfh``[``'``log_price``'``]`

`model1_log_log` `=` `LinearRegression``(``)``.``fit``(``X1_log``,` `y_log``)`

此模型的系数和预测值不能直接与使用线性特征拟合的模型进行比较，因为其单位是美元和平方英尺的对数，而不是美元和平方英尺。

接下来，我们通过图表检查残差和预测值：

`prediction` `=` `model1_log_log``.``predict``(``X1_log``)`
`error` `=` `y_log` `-` `prediction` 

`fig` `=` `px``.``scatter``(``x``=``prediction``,` `y``=``error``,`
                 `labels``=``dict``(``x``=``'``Predicted sale price (log USD)``'``,` `y``=``'``Error``'``)``,`
                 `width``=``350``,` `height``=``250``)`

`fig``.``add_hline``(``0``,` `line_width``=``2``,` `line_dash``=``'``dash``'``,` `opacity``=``1``)`
`fig``.``show``(``)`

残差图看起来合理，但其中包含数千个点，这使得难以看到曲线。

为了查看是否有助于增加额外的变量，我们可以绘制拟合模型的残差图针对不在模型中的变量。如果看到模式，那就表明我们可能想要包括这个额外特征或其转换。之前，我们发现价格分布与房屋所在城市相关，因此让我们检查残差与城市之间的关系：

这个图表显示了错误的分布似乎受到城市的影响。理想情况下，每个城市的箱线图中位数应该与 y 轴上的 0 对齐。然而，皮德蒙特出售的房屋超过 75%存在正误差，这意味着实际销售价格高于预测值。而在另一个极端，里士满超过 75%的销售价格低于预测值。这些模式表明我们应该在模型中包括城市变量。从背景来看，地理位置影响销售价格是有道理的。在接下来的部分，我们展示了如何将名义变量纳入线性模型。

分类测量的特征工程

我们第一次拟合的模型是第四章中的常数模型。在那里，我们最小化平方损失以找到最适合的常数：

min c ∑ i ( y i − c ) 2

我们可以考虑以类似的方式在模型中包含名义特征。也就是说，我们找到每个数据子组中最适合的常数，对应于一个类别：

min c B ∑ i ∈ Berkeley ( y i − c B ) 2 min c L ∑ i ∈ Lamorinda ( y i − c L ) 2 min c P ∑ i ∈ Piedmont ( y i − c P ) 2 min c R ∑ i ∈ Richmond ( y i − c R ) 2

另一种描述这种模型的方式是独热编码。

独热编码将分类特征转换为多个只有 0 或 1 值的数值特征。为了对一个特征进行独热编码，我们创建新的特征，每个唯一的类别对应一个新的特征。在本例中，由于有四个城市——伯克利、Lamorinda、皮德蒙特和里士满——我们在一个设计矩阵中创建了四个新特征，称为 X c i t y 。 X c i t y 中的每一行包含一个值为 1，它出现在与城市对应的列中。 图 15-6 说明了这一概念。

图 15-6. 对一个分类特征进行独热编码（左）及其生成的设计矩阵（右）

现在我们可以简洁地表示模型如下：

θ B x i , B + θ L x i , L + θ P x i , P + θ R x i , R

在这里，我们用 B ，L ，P 和 R 对设计矩阵的列进行了索引，而不是 j ，以明确表示每一列代表一个只有 0 和 1 的列，例如，如果第 i 个房屋位于皮德蒙特，则 x i , P 为 1。

注意

独热编码创建的特征仅具有 0-1 值。这些特征也被称为虚拟变量或指示变量。在计量经济学中更常用“虚拟变量”这一术语，在统计学中更常用“指示变量”。

我们的目标是最小化关于 θ 的最小二乘损失：

‖ y − X θ ‖ 2 = ∑ i ( y i − θ B x i , B + θ L x i , L + θ P x i , P + θ R x i , R ) 2 = ∑ i ∈ B e r k e l e y ( y i − θ B x i , B ) 2 + ∑ i ∈ L a m o r i n d a ( y i − θ L x i , L ) 2 + ∑ i ∈ P i e d m o n t ( y i − θ P x i , P ) 2 + ∑ i ∈ R i c h m o n d ( y i − θ R x i , R ) 2

其中 θ 是列向量 [ θ B , θ L , θ P , θ R ] 。注意，这个最小化转化为四个最小化，每个城市对应一个。这正是我们在本节开始时提到的思路。

我们可以使用 OneHotEncoder 创建这个设计矩阵：

`from` `sklearn``.``preprocessing` `import` `OneHotEncoder`

`enc` `=` `OneHotEncoder``(`
    `# categories argument sets column order`
    `categories``=``[``[``"``Berkeley``"``,` `"``Lamorinda``"``,` `"``Piedmont``"``,` `"``Richmond``"``]``]``,`
    `sparse``=``False``,`
`)`

`X_city` `=` `enc``.``fit_transform``(``sfh``[``[``'``city``'``]``]``)`

`categories_city``=``[``"``Berkeley``"``,``"``Lamorinda``"``,` `"``Piedmont``"``,` `"``Richmond``"``]`
`X_city_df` `=` `pd``.``DataFrame``(``X_city``,` `columns``=``categories_city``)`

`X_city_df`

	伯克利	Lamorinda	皮德蒙特	里士满
0	1.0	0.0	0.0	0.0
1	1.0	0.0	0.0	0.0
2	1.0	0.0	0.0	0.0
...	...	...	...	...
2664	0.0	0.0	0.0	1.0
2665	0.0	0.0	0.0	1.0
2666	0.0	0.0	0.0	1.0

2667 rows × 4 columns

让我们使用这些独热编码特征拟合一个模型：

`y_log` `=` `sfh``[``'``log_price``'``]`

`model_city` `=` `LinearRegression``(``fit_intercept``=``False``)``.``fit``(``X_city_df``,` `y_log``)`

并且检查多重 R 2 :

R-square for city model: 0.57

如果我们只知道房屋所在的城市，该模型能够相当不错地估计其销售价格。以下是拟合的系数：

`model_city``.``coef_`

array([5.87, 6.03, 6.1 , 5.67])

正如盒图所示，估计的销售价格（以对数$表示）取决于城市。但是，如果我们知道房屋大小和城市，我们应该会有一个更好的模型。我们之前看到，简单的对数模型可以合理解释销售价格与房屋大小的关系，因此我们期望城市特征（作为独热编码变量）应该进一步改进模型。

这样的模型如下所示：

y i ≈ θ 1 x i + θ B x i , B + θ L x i , L + θ P x i , P + θ R x i , R

注意，这个模型描述了对数价格（表示为y）和对数大小（表示为x）之间的关系，对于每个城市的对数大小都具有相同的系数。但截距项取决于城市：

y i ≈ θ 1 x i + θ B for houses in Berkeley y i ≈ θ 1 x i + θ L for houses in Lamorinda y i ≈ θ 1 x i + θ P for houses in Piedmont y i ≈ θ 1 x i + θ R for houses in Richmond

接下来，我们制作一个散点图的多面板图，每个城市一个，看看这种关系大致成立：

`fig` `=` `px``.``scatter``(``sfh``,` `x``=``'``log_bsqft``'``,` `y``=``'``log_price``'``,` 
                 `facet_col``=``'``city``'``,` `facet_col_wrap``=``2``,`
                 `labels``=``{``'``log_bsqft``'``:``'``Building size (log ft²)``'``,`
                         `'``log_price``'``:``'``Sale price (log USD)``'``}``,`
                 `width``=``500``,` `height``=``400``)`

`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`
`fig`

在散点图中可以看出这种偏移。我们将两个设计矩阵连接在一起，以拟合包含大小和城市的模型：

`X_size` `=` `sfh``[``'``log_bsqft``'``]` 

`X_city_size` `=` `pd``.``concat``(``[``X_size``.``reset_index``(``drop``=``True``)``,` `X_city_df``]``,` `axis``=``1``)`
`X_city_size``.``drop``(``0``)`

	log_bsqft	伯克利	Lamorinda	皮德蒙特	里士满
1	3.14	1.0	0.0	0.0	0.0
2	3.31	1.0	0.0	0.0	0.0
3	2.96	1.0	0.0	0.0	0.0
...	...	...	...	...	...
2664	3.16	0.0	0.0	0.0	1.0
2665	3.47	0.0	0.0	0.0	1.0
2666	3.44	0.0	0.0	0.0	1.0

2666 rows × 5 columns

现在让我们拟合一个包含定量特征（房屋大小）和定性特征（城市）的模型：

`model_city_size` `=` `LinearRegression``(``fit_intercept``=``False``)``.``fit``(``X_city_size``,` `y_log``)`

这些截距反映了哪些城市的房屋更贵，即使考虑到房屋的大小：

`model_city_size``.``coef_`

array([0.62, 3.89, 3.98, 4.03, 3.75])

R-square for city and log(size):  0.79

这个拟合模型包括名义变量city和对数变换后的房屋大小，比仅包含房屋大小的简单对数模型以及为每个城市拟合常数的模型都要好。

注意，我们从模型中删除了截距项，以便每个子组都有自己的截距。然而，一个常见的做法是从设计矩阵中删除一个独热编码特征，并保留截距。例如，如果我们删除伯克利房屋的特征并添加截距，那么模型就是：

θ 0 + θ 1 x i + θ L x i , L + θ P x i , P + θ R x i , R

在这种表示中，虚拟变量的系数意义在这个表述中已经改变了。例如，考虑伯克利和皮德蒙特的房屋的方程式：

θ 0 + θ 1 x i for a house in Berkeley θ 0 + θ 1 x i + θ P for a house in Piedmont

在这个表示中，截距θ 0 是伯克利房屋的，而系数θ P 衡量皮德蒙特房屋与伯克利房屋之间的典型价格差异。在这种表述中，我们可以更容易地将θ P 与 0 比较，看看这两个城市的平均价格是否基本相同。

如果我们包括截距和所有城市变量，则设计矩阵的列是线性相关的，这意味着我们无法解出系数。我们的预测在任何情况下都将相同，但不会有唯一的最小化解。

当我们包含两个分类变量的独热编码时，我们还喜欢采用删除一个虚拟变量并包含截距项的模型表示。这种做法保持了系数解释的一致性。

我们演示如何使用statsmodels库构建一个具有两组虚拟变量的模型。该库使用公式语言描述要拟合的模型，因此我们无需自己创建设计矩阵。我们导入公式 API：

`import` `statsmodels``.``formula``.``api` `as` `smf`

首先，让我们重复使用名义变量city和房屋大小拟合模型，以展示如何使用公式语言并比较结果：

`model_size_city` `=` `smf``.``ols``(``formula``=``'``log_price ~ log_bsqft + city``'``,`
                          `data``=``sfh``)``.``fit``(``)`

提供给formula参数的字符串描述了要拟合的模型。该模型以log_price作为结果，以log_bsqft和city的线性组合作为解释变量进行拟合。请注意，我们无需创建虚拟变量来拟合模型。方便地，smf.ols为我们执行了城市特征的独热编码。以下模型的拟合系数包括截距项，并且省略了伯克利指示变量：

`print``(``model_size_city``.``params``)`

Intercept            3.89
city[T.Lamorinda]    0.09
city[T.Piedmont]     0.14
city[T.Richmond]    -0.15
log_bsqft            0.62
dtype: float64

如果我们想要去除截距，我们可以在公式中添加 -1，这是一个指示从设计矩阵中去除 1 列的约定。在这个特定的例子中，所有独热编码特征所张成的空间等同于 1 向量和除了一个虚拟变量之外的所有虚拟变量所张成的空间，因此拟合是相同的。但是，系数是不同的，因为它们反映了设计矩阵的不同参数化：

`smf``.``ols``(``formula``=``'``log_price ~ log_bsqft + city - 1``'``,` `data``=``sfh``)``.``fit``(``)``.``params`

city[Berkeley]     3.89
city[Lamorinda]    3.98
city[Piedmont]     4.03
city[Richmond]     3.75
log_bsqft          0.62
dtype: float64

此外，我们可以在城市和大小变量之间添加交互项，以允许每个城市对大小具有不同的系数。我们在公式中指定此项，通过添加术语log_bsqft:city。我们在这里不详细说明。

现在让我们拟合一个具有两个分类变量的模型：卧室数量和城市。回想一下，我们之前重新分配了卧室数量大于 6 的卧室数量为 6，这实际上将 6、7、8、… 折叠到类别 6+ 中。我们可以在价格（对数 $）按卧室数量的箱线图中看到这种关系：

`px``.``box``(``sfh``,` `x``=``"``br``"``,` `y``=``"``log_price``"``,` `width``=``450``,` `height``=``250``,`
       `labels``=``{``'``br``'``:``'``Number of bedrooms``'``,``'``log_price``'``:``'``Sale price (log USD)``'``}``)`

关系看起来不是线性的：每增加一个卧室，销售价格并不会以相同的金额增加。鉴于卧室数量是离散的，我们可以将此特征视为分类的，这样每个卧室编码都可以为成本贡献不同的金额：

`model_size_city_br` `=` `smf``.``ols``(``formula``=``'``log_price ~ log_bsqft + city + C(br)``'``,`
                             `data``=``sfh``)``.``fit``(``)`

我们在公式中使用了术语C(br)来指示我们希望将卧室数量（数值型）视为分类变量对待。

让我们检查拟合的多重R 2：

`model_size_city_br``.``rsquared``.``round``(``2``)`

0.79

尽管我们添加了五个更多的独热编码特征，但多元 R 2 并没有增加。多元 R 2 根据模型参数数量进行调整，从这个度量来看，并不比之前仅包含城市和大小的模型更好。

在本节中，我们介绍了定性特征的特征工程。我们看到了一种独热编码技术，它让我们在线性模型中包含分类数据，并为模型参数提供了自然的解释。

总结

线性模型帮助我们描述特征之间的关系。我们讨论了简单线性模型，并将其扩展到多变量线性模型。在此过程中，我们应用了在建模中广泛有用的数学技术—用微积分来最小化简单线性模型的损失，用矩阵几何来处理多变量线性模型。

线性模型可能看起来很基础，但今天它们被用于各种任务。它们足够灵活，可以允许我们包括分类特征以及变量的非线性转换，如对数转换、多项式和比率。线性模型具有广泛的可解释性，非技术人员也能理解，同时足够复杂以捕捉数据中许多常见模式的特点。

诱人的做法是将所有可用的变量都放入模型中以获得“可能的最佳拟合”。但我们应该记住最小二乘法的几何特性来拟合模型。要记住，p 个解释变量可以被看作是 n 维空间中的 p 个向量，如果这些向量高度相关，那么在这个空间上的投影将类似于在由较少向量组成的较小空间上的投影。这意味着：

• 添加更多变量可能不会显著改善模型。

• 解释系数可能会很困难。

• 几种模型可以同样有效地预测/解释响应变量。

如果我们关心进行推断，希望解释/理解模型，那么我们应该倾向于使用简单的模型。另一方面，如果我们的主要关注是模型的预测能力，那么我们往往不会关心系数的数量及其解释。但是这种“黑箱”方法可能导致模型过度依赖数据中的异常值或在其他方面不足。因此，在预测可能对人们有害时，要小心使用这种方法。

在本章中，我们以描述性方式使用了线性模型。我们介绍了一些决定何时在模型中包含特征的概念，通过检查残差的模式、比较标准误差的大小和多重R 2 的变化来做出决策。通常情况下，我们选择了一个更简单、更容易解释的模型。在下一章中，我们将探讨其他更正式的工具，用于选择模型中要包含的特征。

第十六章：模型选择

到目前为止，当我们拟合模型时，我们已经采用了几种策略来决定包含哪些特征：

• 使用残差图评估模型拟合度。

• 将统计模型与物理模型连接起来。

• 保持模型简单。

• 比较在日益复杂的模型之间残差标准差和 MSE 的改进。

例如，在我们检查单变量模型中的上升流动性时，我们发现残差图中有曲率。增加第二个变量大大改善了平均损失（MSE 以及相关的多个R 2），但残差中仍然存在一些曲率。一个七变量模型在降低 MSE 方面几乎没有比两变量模型更多的改进，因此尽管两变量模型仍然显示出残差中的一些模式，我们选择了这个更简单的模型。

另一个例子是，在我们建立一只驴的体重模型时，在第十八章中，我们将从物理模型中获取指导。我们将忽略驴的附肢，并借鉴桶和驴身体的相似性开始拟合一个模型，通过其长度和围度（类似于桶的高度和周长）来解释体重。然后，我们将继续通过添加与驴的物理状况和年龄相关的分类特征来调整该模型，合并类别，并排除其他可能的特征，以保持模型简单。

我们在构建这些模型时所做的决策基于判断，而在本章中，我们将这些决策与更正式的标准结合起来。首先，我们提供一个示例，说明为什么在模型中包含太多特征通常不是一个好主意。这种现象称为过度拟合，经常导致模型过于贴近数据并捕捉到数据中的一些噪音。然后，当新的观察到来时，预测结果比较简单模型更糟糕。本章的其余部分提供了一些技术，如训练集-测试集分割、交叉验证和正则化，来限制过度拟合的影响。当模型可能包含大量潜在特征时，这些技术尤为有用。我们还提供一个合成示例，我们在其中了解到真实模型，以解释模型方差和偏差的概念及其与过拟合和欠拟合的关系。

过度拟合

当我们有许多可用于模型的特征时，选择包括或排除哪些特征会迅速变得复杂。在上升移动性的例子中，在第十五章中，我们选择了七个变量中的两个来拟合模型，但是对于一个双变量模型，我们可以检查并拟合 21 对特征。如果考虑所有一、二、...、七个变量模型，还有超过一百种模型可供选择。检查数百个残差图以决定何时简单到足够，以及确定一个模型是相当困难的。不幸的是，最小化均方误差的概念并非完全有用。当我们向模型添加一个变量时，均方误差通常会变小。回顾模型拟合的几何视角（第十五章），添加一个特征到模型中会添加一个 n -维向量到特征空间中，并且结果向量与其在由解释变量张成的空间内的投影之间的误差会减小。我们可能认为这是一件好事，因为我们的模型更接近数据，但过度拟合存在危险。

当模型过于紧密地跟随数据并捕捉到结果中的随机噪声变化时，就会发生过度拟合。当这种情况发生时，新的观察结果就无法很好地预测。举个例子可以帮助澄清这个概念。

示例：能源消耗

在这个例子中，我们研究一个可以下载的数据集，其中包含明尼苏达州一个私人住宅的公用事业账单信息。我们记录了一个家庭每月的气体使用量（立方英尺）和该月的平均温度（华氏度）。^(1) 我们首先读取数据：

`heat_df` `=` `pd``.``read_csv``(``"``data/utilities.csv``"``,` `usecols``=``[``"``temp``"``,` `"``ccf``"``]``)`
`heat_df`

	温度	立方英尺
0	29	166
1	31	179
2	15	224
...	...	...
96	76	11
97	55	32
98	39	91

99 rows × 2 columns

我们将从查看气体消耗作为温度函数的散点图开始：

这个关系显示了曲率（左图），但当我们尝试通过对数变换将其变成直线（右图）时，在低温区域出现了不同的曲率。此外，还有两个异常点。当我们查阅文档时，发现这些点代表记录错误，因此我们将它们移除。

看看二次曲线能否捕捉气体使用量和温度之间的关系。多项式仍然被认为是线性模型。它们在其多项式特征中是线性的。例如，我们可以将二次模型表示为：

θ 0 + θ 1 x + θ 2 x 2

这个模型对特征 x 和 x 2 是线性的，在矩阵表示中，我们可以将这个模型写成 X θ ，其中 X 是设计矩阵：

⌈ 1 x 1 x 1 2 1 x 2 x 2 2 ⋮ ⋮ ⋮ 1 x n x n 2 ⌉

我们可以使用scikit-learn中的Polynomial​Fea⁠tures工具创建设计矩阵的多项式特征：

`y` `=` `heat_df``[``'``ccf``'``]`
`X` `=` `heat_df``[``[``'``temp``'``]``]`

`from` `sklearn``.``preprocessing` `import` `PolynomialFeatures`

`poly` `=` `PolynomialFeatures``(``degree``=``2``,` `include_bias``=``False``)`
`poly_features` `=` `poly``.``fit_transform``(``X``)`
`poly_features`

array([[  29.,  841.],
       [  31.,  961.],
       [  15.,  225.],
       ...,
       [  76., 5776.],
       [  55., 3025.],
       [  39., 1521.]])

我们将参数include_bias设置为False，因为我们计划在scikit-learn中使用LinearRegression方法拟合多项式，默认情况下会在模型中包括常数项。我们用以下方法拟合多项式：

`from` `sklearn``.``linear_model` `import` `LinearRegression`

`model_deg2` `=` `LinearRegression``(``)``.``fit``(``poly_features``,` `y``)`

为了快速了解拟合的质量，让我们在散点图上叠加拟合的二次曲线，并查看残差：

二次多项式很好地捕捉了数据中的曲线，但残差显示出在 70°F 到 80°F 温度范围内略微上升的趋势，这表明有些拟合不足。此外，残差中也有些漏斗形状，在较冷的月份中，燃气消耗的变化性更大。我们可能会预期这种行为，因为我们只有月均温度。

为了比较，我们使用更高阶的多项式拟合了几个模型，并集体检查拟合曲线：

`poly12` `=` `PolynomialFeatures``(``degree``=``12``,` `include_bias``=``False``)`
`poly_features12` `=` `poly12``.``fit_transform``(``X``)`

`degrees` `=` `[``1``,` `2``,` `3``,` `6``,` `8``,` `12``]`

`mods` `=` `[``LinearRegression``(``)``.``fit``(``poly_features12``[``:``,` `:``deg``]``,` `y``)`
        `for` `deg` `in` `degrees``]`

警告

我们在本节中使用多项式特征来演示过拟合，但直接拟合x , x 2 , x 3 , …等多项式在实践中是不可取的。不幸的是，这些多项式特征往往高度相关。例如，能源数据中x和x 2之间的相关性为 0.98。高度相关的特征会导致不稳定的系数，即 x 值的微小变化可能会导致多项式系数的大幅变化。此外，当 x 值较大时，正规方程的条件很差，系数的解释和比较可能会很困难。

更好的做法是使用构造成彼此正交的多项式。这些多项式填充与原始多项式相同的空间，但它们彼此不相关，并提供更稳定的拟合。

让我们将所有的多项式拟合放在同一张图上，这样我们可以看到高阶多项式的弯曲越来越奇怪：

我们还可以在单独的面板中可视化不同的多项式拟合：

左上方面的一次曲线（直线）未能捕捉数据中的曲线模式。二次曲线开始捕捉，而三次曲线看起来有所改进，但请注意图表右侧的上升弯曲。六次、八次和十二次的多项式越来越贴近数据，因为它们变得越来越曲折。这些多项式似乎适应数据中的虚假波动。总体来看，这六条曲线说明了欠拟合和过拟合。左上角的拟合线欠拟合，完全错过了曲线。而右下角的十二次多项式明显过拟合，呈现了我们认为在这种情况下无意义的蜿蜒模式。

一般来说，随着特征的增加，模型变得更加复杂，均方误差（MSE）下降，但与此同时，拟合的模型变得越来越不稳定和对数据敏感。当我们过度拟合时，模型过于紧密地跟随数据，对新观测的预测效果很差。评估拟合模型的一种简单技术是在新数据上计算 MSE，这些数据未用于建模。由于通常情况下我们无法获取更多数据，因此我们会保留一些原始数据来评估拟合的模型。这个技术是下一节的主题。

训练-测试分离

虽然我们希望在建立模型时使用所有数据，但我们也想了解模型在新数据上的表现。通常情况下，我们无法收集额外的数据来评估模型，因此我们会将部分数据保留下来，称为测试集，以代表新数据。其余的数据称为训练集，我们使用这部分数据来建立模型。然后，在选择了模型之后，我们提取测试集，并查看模型（在训练集上拟合）在测试集中预测结果的表现。图 16-1 演示了这个概念。

图 16-1. 训练-测试分离将数据分为两部分：训练集用于建立模型，测试集评估该模型

通常，测试集包含数据的 10%到 25%。从图表中可能不清楚的是，这种分割通常是随机进行的，因此训练集和测试集彼此相似。

我们可以使用第 15 章 中介绍的概念来描述这个过程。设计矩阵，X ，和结果，y ，各自被分成两部分；标记为X T 的设计矩阵和相应的结果，标记为y T ，形成训练集。我们通过这些数据最小化θ 的平均平方损失：

min θ ‖ y T − X T θ ‖ 2

最小化训练误差的系数θ ^ T 用于预测测试集的结果，其中标记为X S 和 y S：

‖ y S − X S θ ^ T ‖ 2

由于X S和y S没有用于构建模型，它们可以合理估计我们可能对新观测到的损失。

我们使用上一节中的气耗多项式模型来演示训练-测试分离。为此，我们执行以下步骤：

1. 将数据随机分为两部分，训练集和测试集。

2. 对训练集拟合几个多项式模型并选择一个。

3. 计算在所选多项式（其系数由训练集拟合）上的测试集的 MSE。

对于第一步，我们使用scikit-learn中的train_test_split方法将数据随机分为两部分，并为模型评估设置了 22 个观测值：

`from` `sklearn``.``model_selection` `import` `train_test_split`

`test_size` `=` `22`

`X_train``,` `X_test``,` `y_train``,` `y_test` `=` `train_test_split``(`
    `X``,` `y``,` `test_size``=``test_size``,` `random_state``=``42``)`

`print``(``f``'``Training set size:` `{``len``(``X_train``)``}``'``)`
`print``(``f``'``Test set size:` `{``len``(``X_test``)``}``'``)`

Training set size: 75
Test set size: 22

与前一节类似，我们将气温与燃气消耗的模型拟合到各种多项式中。但这次，我们只使用训练数据：

`poly` `=` `PolynomialFeatures``(``degree``=``12``,` `include_bias``=``False``)`
`poly_train` `=` `poly``.``fit_transform``(``X_train``)`

`degree` `=` `np``.``arange``(``1``,``13``)`

`mods` `=` `[``LinearRegression``(``)``.``fit``(``poly_train``[``:``,` `:``j``]``,` `y_train``)`
        `for` `j` `in` `degree``]`

我们找出了每个模型的 MSE：

`from` `sklearn``.``metrics` `import` `mean_squared_error`

`error_train` `=` `[`
    `mean_squared_error``(``y_train``,` `mods``[``j``]``.``predict``(``poly_train``[``:``,` `:` `(``j` `+` `1``)``]``)``)`
    `for` `j` `in` `range``(``12``)`
`]`

为了可视化 MSE 的变化，我们将每个拟合的多项式的 MSE 绘制成其次数的图：

`px``.``line``(``x``=``degree``,` `y``=``error_train``,` `markers``=``True``,`
        `labels``=``dict``(``x``=``'``Degree of polynomial``'``,` `y``=``'``Train set MSE``'``)``,`
        `width``=``350``,` `height``=``250``)`

注意到随着模型复杂度的增加，训练误差逐渐减小。我们之前看到高阶多项式显示出了我们认为不反映数据中潜在结构的起伏行为。考虑到这一点，我们可能会选择一个更简单但 MSE 显著减小的模型。这可能是 3、4 或 5 次方。让我们选择 3 次方，因为这三个模型在 MSE 方面的差异非常小，而且它是最简单的。

现在我们已经选择了我们的模型，我们使用测试集提供了对其 MSE 的独立评估。我们为测试集准备设计矩阵，并使用在训练集上拟合的 3 次多项式来预测测试集中每一行的结果。最后，我们计算了测试集的 MSE：

`poly_test` `=` `poly``.``fit_transform``(``X_test``)`
`y_hat` `=` `mods``[``2``]``.``predict``(``poly_test``[``:``,` `:``3``]``)`

`mean_squared_error``(``y_test``,` `y_hat``)`

307.44460133992294

该模型的均方误差（MSE）比在训练数据上计算的 MSE 要大得多。这说明了在使用相同数据来拟合和评估模型时所存在的问题：MSE 并不充分反映出对新观测的 MSE。为了进一步说明过拟合问题，我们计算了这些模型的测试误差：

`error_test` `=` `[`
    `mean_squared_error``(``y_test``,` `mods``[``j``]``.``predict``(``poly_test``[``:``,` `:` `(``j` `+` `1``)``]``)``)`
    `for` `j` `in` `range``(``12``)`
`]`

在实践中，我们不会在承诺模型之前查看测试集。在训练集上拟合模型并在测试集上评估它之间交替可以导致过拟合。但出于演示目的，我们绘制了我们拟合的所有多项式模型在测试集上的 MSE：

注意，对于所有模型而言，测试集的均方误差（MSE）大于训练集的均方误差，而不仅仅是我们选择的模型。更重要的是，注意当模型从欠拟合到更好地拟合数据曲线时，测试集的均方误差最初是下降的。然后，随着模型复杂度的增加，测试集的均方误差增加。这些更复杂的模型对训练数据过拟合，导致预测测试集时出现较大的误差。这种现象的一个理想化示意图如图 16-2 所示。

图 16-2. 随着模型复杂度的增加，训练集的误差减少，而测试集的误差增加

测试数据提供新观察的预测误差评估。仅在我们已经选择了模型之后才使用测试集是至关重要的。否则，我们会陷入使用相同数据选择和评估模型的陷阱中。在选择模型时，我们回归到了简单性的论点，因为我们意识到越来越复杂的模型往往会过拟合。然而，我们也可以扩展训练-测试方法来帮助选择模型。这是下一节的主题。

交叉验证

我们可以使用训练-测试范式来帮助选择模型。其思想是进一步将训练集分成单独的部分，在其中一个部分上拟合模型，然后在另一个部分上评估模型。这种方法称为交叉验证。我们描述的是一种版本，称为k -折叠交叉验证。图 16-3 展示了这种数据划分背后的思想。

图 16-3. 一个例子展示了五折交叉验证，其中训练集被分为五部分，轮流用于验证在其余数据上构建的模型

交叉验证可以帮助选择模型的一般形式。这包括多项式的阶数，模型中的特征数量，或者正则化惩罚的截止（在下一节中介绍）。k -折叠交叉验证的基本步骤如下：

1. 将训练集分成k个大致相同大小的部分；每部分称为一个折叠。使用与创建训练集和测试集相同的技术来创建这些折叠。通常情况下，我们随机划分数据。

2. 将一个折叠保留作为测试集：

   • 在剩余的训练数据上拟合所有模型（训练数据减去特定折叠的数据）。

   • 使用您保留的折叠来评估所有这些模型。

3. 重复此过程共k次，每次将一个折叠保留出来，使用剩余的训练集来拟合模型，并在保留的折叠上评估模型。

4. 合并每个模型在折叠中的拟合误差，并选择具有最小误差的模型。

这些拟合模型在不同的折叠中不会具有相同的系数。例如，当我们拟合一个三次多项式时，我们对k个折叠中的 MSE 取平均值，得到三次拟合多项式的平均 MSE。然后我们比较这些 MSE，并选择具有最低 MSE 的多项式次数。在三次多项式中，x，x2和x3项的实际系数在每个k个拟合中是不同的。一旦选择了多项式次数，我们使用所有训练数据重新拟合模型，并在测试集上评估它。（在选择模型的任何早期步骤中，我们没有使用测试集。）

通常，我们使用 5 或 10 个折叠。另一种流行的选择是将一个观察结果放入每个折叠中。这种特殊情况称为留一法交叉验证。其流行之处在于调整最小二乘拟合以减少一个观察结果的简单性。

通常，k折交叉验证需要一些计算时间，因为我们通常必须为每个折叠从头开始重新拟合每个模型。scikit-learn库提供了一个方便的sklearn.model_selection.KFold类来实现k折交叉验证。

为了让你了解 k 折交叉验证的工作原理，我们将在燃气消耗示例上演示这种技术。但是，这次我们将拟合不同类型的模型。在数据的原始散点图中，看起来点落在两条连接的线段上。在寒冷的温度下，燃气消耗与温度之间的关系看起来大致是负斜率，约为−4 立方英尺/度，而在温暖的月份，关系则似乎几乎是平坦的。因此，我们可以拟合一条弯曲的线条而不是拟合多项式。

让我们从 65 度处拟合一条有弯曲的线。为此，我们创建一个特征，使得温度高于 65°F 的点具有不同的斜率。该模型是：

y = θ 0 + θ 1 x + θ 2 ( x − 65 ) +

在这里，()+代表“正部分”，所以当x小于 65 时，它评估为 0，当x大于或等于 65 时，它就是x−65。我们创建这个新特征并将其添加到设计矩阵中：

`y` `=` `heat_df``[``"``ccf``"``]`
`X` `=` `heat_df``[``[``"``temp``"``]``]`
`X``[``"``temp65p``"``]` `=` `(``X``[``"``temp``"``]` `-` `65``)` `*` `(``X``[``"``temp``"``]` `>``=` `65``)`

然后，我们使用这两个特征拟合模型：

`bend_index` `=` `LinearRegression``(``)``.``fit``(``X``,` `y``)`

让我们将这条拟合的“曲线”叠加在散点图上，看看它如何捕捉数据的形状：

这个模型似乎比多项式更好地拟合了数据。但是可能有许多种弯线模型。线条可能在 55 度或 60 度处弯曲等。我们可以使用 k 折交叉验证来选择线条弯曲的温度值。让我们考虑在 40 , 41 , 42 , … , 68 , 69 度处弯曲的模型。对于这些模型的每一个，我们需要创建额外的特征来使线条在那里弯曲：

`bends` `=` `np``.``arange``(``40``,` `70``,` `1``)`

`for` `i` `in` `bends``:`
    `col` `=` `"``temp``"` `+` `i``.``astype``(``"``str``"``)` `+` `"``p``"`
    `heat_df``[``col``]` `=` `(``heat_df``[``"``temp``"``]` `-` `i``)` `*` `(``heat_df``[``"``temp``"``]` `>``=` `i``)`
`heat_df`

	temp	ccf	temp40p	temp41p	...	temp66p	temp67p	temp68p	temp69p
0	29	166	0	0	...	0	0	0	0
1	31	179	0	0	...	0	0	0	0
2	15	224	0	0	...	0	0	0	0
...	...	...	...	...	...	...	...	...	...
96	76	11	36	35	...	10	9	8	7
97	55	32	15	14	...	0	0	0	0
98	39	91	0	0	...	0	0	0	0

97 rows × 32 columns

交叉验证的第一步是创建我们的训练集和测试集。与之前一样，我们随机选择 22 个观测值放入测试集。这样剩下 75 个观测值作为训练集：

`y` `=` `heat_df``[``'``ccf``'``]`
`X` `=` `heat_df``.``drop``(``columns``=``[``'``ccf``'``]``)`

`test_size` `=` `22`

`X_train``,` `X_test``,` `y_train``,` `y_test` `=` `train_test_split``(`
    `X``,` `y``,` `test_size``=``test_size``,` `random_state``=``0``)`

现在我们可以将训练集分成折叠。我们使用三个折叠，以便每个折叠中有 25 个观测值。对于每个折叠，我们拟合 30 个模型，每个模型对应线条中的一个弯曲点。对于这一步骤，我们使用 scikit-learn 中的 KFold 方法来划分数据：

`from` `sklearn``.``model_selection` `import` `KFold`

`kf` `=` `KFold``(``n_splits``=``3``,` `shuffle``=``True``,` `random_state``=``42``)`

`validation_errors` `=` `np``.``zeros``(``(``3``,` `30``)``)`

`def` `validate_bend_model``(``X``,` `y``,` `X_valid``,` `y_valid``,` `bend_index``)``:`
    `model` `=` `LinearRegression``(``)``.``fit``(``X``.``iloc``[``:``,` `[``0``,` `bend_index``]``]``,` `y``)`
    `predictions` `=` `model``.``predict``(``X_valid``.``iloc``[``:``,` `[``0``,` `bend_index``]``]``)`
    `return` `mean_squared_error``(``y_valid``,` `predictions``)`

`for` `fold``,` `(``train_idx``,` `valid_idx``)` `in` `enumerate``(``kf``.``split``(``X_train``)``)``:`
    `cv_X_train``,` `cv_X_valid` `=` `(``X_train``.``iloc``[``train_idx``,` `:``]``,`
                              `X_train``.``iloc``[``valid_idx``,` `:``]``)`
    `cv_Y_train``,` `cv_Y_valid` `=` `(``y_train``.``iloc``[``train_idx``]``,`
                              `y_train``.``iloc``[``valid_idx``]``)`

    `error_bend` `=` `[`
        `validate_bend_model``(`
            `cv_X_train``,` `cv_Y_train``,` `cv_X_valid``,` `cv_Y_valid``,` `bend_index`
        `)`
        `for` `bend_index` `in` `range``(``1``,` `31``)`
    `]`

    `validation_errors``[``fold``]``[``:``]` `=` `error_bend`

然后我们找到三个折叠中的平均验证误差，并将它们绘制在弯曲位置的图表上：

`totals` `=` `validation_errors``.``mean``(``axis``=``0``)`

MSE 在 57 到 60 度之间看起来相当平缓。最小值出现在 58 度，因此我们选择那个模型。为了评估这个模型在测试集上的表现，我们首先在整个训练集上以 58 度拟合弯线模型：

`bent_final` `=` `LinearRegression``(``)``.``fit``(`
    `X_train``.``loc``[``:``,` `[``"``temp``"``,` `"``temp58p``"``]``]``,` `y_train`
`)`

然后我们使用拟合的模型来预测测试集的气体消耗：

`y_pred_test` `=` `bent_final``.``predict``(``X_test``.``loc``[``:``,` `[``"``temp``"``,` `"``temp58p``"``]``]``)`

`mean_squared_error``(``y_test``,` `y_pred_test``)`

71.40781435952441

让我们将弯线拟合叠加到散点图上，并检查残差，以了解拟合质量：

拟合曲线看起来合理，并且残差比多项式拟合的要小得多。

注意

在本节的教学目的上，我们使用 KFold 手动将训练数据分为三个折叠，然后使用循环找到模型验证误差。在实践中，我们建议使用 sklearn.model_selection.GridSearchCV 与 sklearn.pipeline.Pipeline 对象，它可以自动将数据分成训练集和验证集，并找到在折叠中平均验证误差最低的模型。

使用交叉验证来管理模型复杂度有几个关键的限制：通常需要使复杂度离散变化，并且可能没有自然的方式来对模型进行排序。与其改变一系列模型的维度，我们可以拟合一个大模型并对系数的大小施加约束。这个概念被称为正则化，将在下一节讨论。

正则化

我们刚刚看到交叉验证如何帮助找到适合的模型维度，从而平衡欠拟合和过拟合。与其选择模型的维度，我们可以构建一个包含所有特征的模型，但限制系数的大小。通过在均方误差上添加一个系数大小的惩罚项来防止过拟合。这个惩罚项称为正则化项，表达式为λ ∑ j = 1 p θ j 2 。我们通过最小化均方误差和这个惩罚项的组合来拟合模型：

1 n ∑ i = 1 n ( y i − x i θ ) 2 + λ ∑ j = 1 p θ j 2

当正则化参数λ很大时，会惩罚大的系数。（通常通过交叉验证来选择。）

对系数的平方进行惩罚称为L 2正则化，或称为岭回归。另一种流行的正则化方法惩罚系数的绝对大小：

1 n ∑ i = 1 n ( y i − x i θ ) 2 + λ ∑ j = 1 p | θ j |

这个L 1正则化的线性模型也被称为套索回归（lasso 代表最小绝对值收缩和选择运算符）。

为了了解正则化的工作原理，让我们考虑极端情况：当λ非常大或接近 0 时（λ从不为负）。当正则化参数很大时，系数会受到严重的惩罚，因此它们会收缩。另一方面，当λ很小时，系数不受限制。实际上，当λ为 0 时，我们回到了普通最小二乘法的世界。当我们考虑通过正则化控制系数大小时，会遇到几个问题：

• 我们不希望对截距项进行正则化。这样一来，一个大的惩罚就会拟合一个常数模型。

• 当特征具有非常不同的尺度时，惩罚可能会对它们产生不同的影响，具有较大值的特征会比其他特征受到更多的惩罚。为了避免这种情况，我们在拟合模型之前将所有特征标准化，使它们的均值为 0，方差为 1。

让我们看一个包含 35 个特征的例子。

模型的偏差和方差

在本节中，我们提供了一种不同的思考过拟合和欠拟合问题的方式。我们进行了模拟研究，从我们设计的模型中生成了合成数据。这样，我们知道真实模型，并可以看到在拟合数据时我们离真实情况有多近。

我们可以构建一个数据的通用模型如下：

y = g ( x ) + ϵ

这个表达式使得模型的两个组成部分很容易看出来：信号 g ( x ) 和噪声 ϵ 。在我们的模型中，我们假设噪声没有趋势或模式，方差恒定，并且每个观测值的噪声是独立的。

例如，让我们取 g ( x ) = sin ⁡ ( x ) + 0.3 x ，噪声来自均值为 0，标准差为 0.2 的正态曲线。我们可以从这个模型生成数据，使用以下函数：

`def` `g``(``x``)``:`
    `return` `np``.``sin``(``x``)` `+` `0.3` `*` `x`

`def` `gen_noise``(``n``)``:`
    `return` `np``.``random``.``normal``(``scale``=``0.2``,` `size``=``n``)`

`def` `draw``(``n``)``:`
    `points` `=` `np``.``random``.``choice``(``np``.``arange``(``0``,` `10``,` `0.05``)``,` `size``=``n``)`
    `return` `points``,` `g``(``points``)` `+` `gen_noise``(``n``)`

让我们生成 50 个数据点 ( x i , y i ) ，i = 1 , … , 50 ，从这个模型中：

`np``.``random``.``seed``(``42``)`

`xs``,` `ys` `=` `draw``(``50``)`
`noise` `=` `ys` `-` `g``(``xs``)`

我们可以绘制我们的数据，因为我们知道真实信号，我们可以找到错误并将它们绘制出来：

左边的图显示 g 作为虚线曲线。我们还可以看到 ( x , y ) 对形成了这条曲线的散点分布。右边的图显示了 50 个点的误差，y − g ( x ) 。请注意，它们没有形成模式。

当我们对数据进行模型拟合时，我们最小化均方误差。让我们用一般性写出这个最小化：

min f ∈ F 1 n ∑ i = 1 n [ y i − f ( x i ) ] 2

最小化是对函数集合 F 进行的。我们在本章中已经看到，这个函数集合可能是 12 阶多项式，或者简单的弯曲线。一个重要的点是真实模型 g 不必是集合中的一个函数。

让我们把 F 定为二次多项式的集合；换句话说，可以表示为 θ 0 + θ 1 x + θ 2 x 2 的函数。由于 g ( x ) = sin ⁡ ( x ) + 0.3 x ，它不属于我们正在优化的函数集合。

让我们对我们的 50 个数据点进行多项式拟合：

`poly` `=` `PolynomialFeatures``(``degree``=``2``,` `include_bias``=``False``)`
`poly_features` `=` `poly``.``fit_transform``(``xs``.``reshape``(``-``1``,` `1``)``)`

`model_deg2` `=` `LinearRegression``(``)``.``fit``(``poly_features``,` `ys``)`

Fitted Model: 0.98 + -0.19x + 0.05x²

再次，我们知道真实模型不是二次的（因为我们建立了它）。让我们绘制数据和拟合曲线：

二次函数不太适合数据，并且也不能很好地表示底层曲线，因为我们选择的模型集（二阶多项式）无法捕捉g中的曲率。

如果我们重复这个过程，并从真实模型中生成另外 50 个点，并将二次多项式拟合到这些数据上，那么二次多项式的拟合系数会因为依赖于新数据集而改变。我们可以多次重复这个过程，并对拟合曲线取平均。这个平均曲线将类似于从我们真实模型中取 50 个数据点拟合的典型二次多项式最佳拟合曲线。为了演示这个概念，让我们生成 25 组 50 个数据点，并对每个数据集拟合一个二次多项式：

`def` `fit``(``n``)``:`
    `xs_new` `=` `np``.``random``.``choice``(``np``.``arange``(``0``,` `10``,` `0.05``)``,` `size``=``n``)`
    `ys_new` `=` `g``(``xs_new``)` `+` `gen_noise``(``n``)`
    `X_new` `=` `xs_new``.``reshape``(``-``1``,` `1``)`
    `mod_new` `=` `LinearRegression``(``)``.``fit``(``poly``.``fit_transform``(``X_new``)``,` `ys_new``)`
    `return` `mod_new``.``predict``(``poly_features_x_full``)``.``flatten``(``)`

`fits` `=` `[``fit``(``50``)` `for` `j` `in` `range``(``25``)``]`

我们可以在图中显示所有 25 个拟合模型以及真实函数g和拟合曲线的平均值f ¯ 。为此，我们使用透明度来区分重叠的曲线：

我们可以看到 25 个拟合的二次多项式与数据有所不同。这个概念被称为模型变异。25 个二次多项式的平均值由实线表示。平均二次多项式与真实曲线之间的差异称为模型偏差。

当信号g不属于模型空间F时，我们有模型偏差。如果模型空间能很好地逼近g，那么偏差就很小。例如，一个 10 次多项式可以很接近我们示例中使用的g。另一方面，正如我们在本章前面看到的，高阶多项式可能会过度拟合数据并且在尝试接近数据时变化很大。模型空间越复杂，拟合模型的变化就越大。使用过于简单的模型会导致高模型偏差（g和f ¯之间的差异），而使用过于复杂的模型可能会导致高模型方差（f ^在f ¯周围的波动）。这个概念被称为偏差-方差权衡。模型选择旨在平衡这些竞争性的拟合不足来源。

总结

在本章中，我们看到当我们最小化均方误差来拟合模型并评估时会出现问题。训练集-测试集分割帮助我们避开这个问题，其中我们用训练集拟合模型，并在设置好的测试数据上评估我们的拟合模型。

“过度使用”测试集非常重要，因此我们保持其与模型分离直到我们决定了一个模型。为了帮助我们做出决定，我们可能使用交叉验证，模拟将数据分为测试和训练集。同样重要的是，只使用训练集进行交叉验证，并将原始测试集远离任何模型选择过程。

正则化采取了不同的方法，通过惩罚均方误差来防止模型过度拟合数据。在正则化中，我们利用所有可用数据来拟合模型，但是缩小系数的大小。

偏差-方差折衷使我们能够更准确地描述本章中看到的建模现象：拟合不足与模型偏差有关；过度拟合导致模型方差增大。在图 16-4 中，x 轴表示模型复杂度，y 轴表示模型不适配的这两个组成部分：模型偏差和模型方差。请注意，随着拟合模型的复杂性增加，模型偏差减少，模型方差增加。从测试误差的角度来看，我们看到这种误差首先减少，然后由于模型方差超过模型偏差的减少而增加。为了选择一个有用的模型，我们必须在模型偏差和模型方差之间取得平衡。

图 16-4. 偏差-方差折衷

如果模型能够完全拟合人口过程，收集更多观察数据会减少偏差。如果模型本质上无法建模人口（如我们的合成示例），即使是无限数据也无法消除模型偏差。就方差而言，收集更多数据也会减少方差。数据科学的一个最新趋势是选择具有低偏差和高内在方差（例如神经网络）的模型，但收集许多数据点以使模型方差足够低以进行准确预测。虽然在实践中有效，但为这些模型收集足够的数据通常需要大量的时间和金钱。

创建更多的特征，无论其是否有用，通常会增加模型的方差。具有许多参数的模型具有许多可能的参数组合，因此比具有少量参数的模型具有更高的方差。另一方面，在模型中添加一个有用的特征（例如，当基础过程是二次的时添加二次特征），可以减少偏差。但即使添加一个无用的特征，很少会增加偏差。

熟悉偏差-方差折衷可以帮助您更好地拟合模型。并且使用诸如训练-测试分割、交叉验证和正则化等技术可以改善这个问题。

建模的另一个部分考虑了拟合系数和曲线的变化。我们可能希望为系数提供置信区间或未来观测的预测带。这些区间和带子给出了拟合模型准确性的感觉。接下来我们将讨论这个概念。

^(1) 这些数据来自于丹尼尔·T·卡普兰（Daniel T. Kaplan）（CreateSpace Independent Publishing Platform, 2009）。

第十七章：推断和预测的理论

当您希望将您的研究结果推广到更大的背景中时，数据需要代表该更大的世界。例如，您可能希望基于传感器读数预测未来时间的空气质量（见第十二章），基于实验结果测试激励是否提高了贡献者的生产力（见第三章），或者构建一个时间间隔估计，用于估计等待公交车的时间（见第五章）。我们在前几章已经涉及了所有这些情境。在本章中，我们将正式化用于预测和推断的框架。

在这个框架的核心是分布的概念，无论是总体、经验（也称为样本）还是概率分布。理解这些分布之间的联系对于假设检验、置信区间、预测带和风险的基础至关重要。我们首先简要回顾了在第三章介绍的瓮模型，然后介绍了假设检验、置信区间和预测带的正式定义。我们在示例中使用模拟，包括引导样本作为特例。我们通过对期望、方差和标准误差的正式定义来结束本章——这些是测试、推断和预测理论中的基本概念。

分布：总体、经验、抽样

总体、抽样和经验分布是我们在对模型进行推断或对新观察进行预测时的重要概念。图 17-1 提供了一个图示，可以帮助区分它们。该图示使用了来自第二章的总体和访问框架以及来自第三章的瓮模型的概念。左侧是我们研究的总体，以一个瓮中的弹珠表示，每个单位一个。我们简化了情况，使得访问框架和总体相同；也就是说，我们可以访问总体中的每个单位。 （当这种情况不成立时所出现的问题在第二章和第三章中有所涉及。）从瓮到样本的箭头表示设计，意味着从框架中选择样本的协议。图示将此选择过程显示为一个机会机制，以从装有难以区分的弹珠的瓮中抽取来表示。图示的右侧，弹珠的集合构成了我们的样本（我们得到的数据）。

图 17-1. 数据生成过程的图示

我们通过考虑仅一个特征的测量来简化图表。在图表中的瓮下方是该特征的总体直方图。总体直方图表示整个总体中数值的分布。在图的最右边，经验直方图展示了我们实际样本的数值分布。请注意，这两个分布在形状上相似。当我们的抽样机制产生代表性样本时，就会出现这种情况。

我们通常对样本测量的摘要感兴趣，比如均值、中位数、简单线性模型的斜率等等。通常，这种摘要统计量是总体参数的估计，比如总体均值或中位数。总体参数在图的左侧显示为θ ∗；而在右侧，从样本计算出的摘要统计量显示为θ ^ 。

我们的样本生成机制很可能在我们再次进行调查时生成不同的数据集。但是如果协议设计良好，我们预期样本仍将类似于总体。换句话说，我们可以从样本计算出的摘要统计量推断总体参数。图中间的抽样分布是统计量的概率分布。它展示了不同样本可能取得的统计量数值及其概率。在第三章中，我们使用模拟来估计了几个示例中的抽样分布。在本章中，我们重访了这些以及前几章的其他示例，以形式化分析。

最后，关于这三个直方图的一个重要点：正如在第十章中介绍的，矩形提供了任何箱中观察结果的比例。在总体直方图的情况下，这是整个总体的比例；对于经验直方图，该区域表示样本中的比例；而对于抽样分布，则表示数据生成机制在该箱中产生样本统计量的机会。

最后，我们通常不了解总体分布或参数，并试图推断参数或预测未见单位的值。其他时候，可以使用样本测试关于总体的推测。测试将是下一节的主题。

假设检验基础

根据我们的经验，假设检验是数据科学中较具挑战性的领域之一——学习起来困难，应用起来也具挑战性。这不一定是因为假设检验深奥技术性强；相反，假设检验可能令人感到反直觉，因为它利用了矛盾。顾名思义，我们通常从一个假设开始进行假设检验：关于我们想要验证的世界的一个陈述。

在理想的情况下，我们将直接证明我们的假设是正确的。不幸的是，我们通常无法获得确定真相所需的所有信息。例如，我们可能假设一种新的疫苗有效，但当代医学尚未完全理解支配疫苗效力的生物学细节。因此，我们转向概率、随机抽样和数据设计的工具。

假设检验之所以令人困惑的一个原因是，它很像“反证法”，我们假设与我们的假设相反是真的，并尝试证明我们观察到的数据与该假设不一致。我们采用这种方式来解决问题，因为通常，某事可能由许多原因导致，但我们只需要一个例子来反驳一个假设。我们称这个“相反假设”为 空假设，我们的原始假设为 备择假设。

使事情变得有点混乱的是，概率的工具并不直接证明或证伪事物。相反，它们告诉我们，在假设的情况下，我们观察到的某事有多么可能或不可能。这就是为什么设计数据收集如此重要的原因。

回顾 J&J 疫苗的随机临床试验（第三章），在这个试验中，有 43,738 人参加了试验，并被随机分成两组。治疗组接受了疫苗，对照组接受了一种假疫苗，称为安慰剂。这种随机分配创建了两个在除了疫苗之外的每个方面都相似的群体。

在这个试验中，治疗组中有 117 人生病，对照组中有 351 人生病。由于我们希望提供令人信服的证据证明疫苗有效，我们从一个空假设开始，即疫苗没有作用，这意味着随机分配导致治疗组中生病的人数如此之少纯属偶然。然后，我们可以使用概率来计算观察到如此少的治疗组生病的机会。概率计算基于一个有 43,738 个彩球的瓮，其中 468 个标记为 1 表示一个生病的人。然后我们发现，21,869 次有放回地从瓮中抽取最多 117 个彩球的概率几乎为零。我们把这个作为拒绝空假设的证据，支持备择假设即疫苗有效。因为 J&J 的实验设计良好，对空假设的拒绝使我们得出结论，即疫苗有效。换句话说，假设的真相取决于我们以及我们有多愿意可能是错误的。

在本节的其余部分，我们将介绍假设检验的四个基本步骤。然后我们提供两个示例，继续介绍来自 第三章 的两个示例，并深入探讨测试的形式。

假设检验有四个基本步骤：

第一步：建立

你已经有了你的数据，并且你想要测试一个特定的模型是否与数据相一致。 所以你指定一个统计量，θ ^，比如样本平均值、样本中零的比例或拟合的回归系数，目的是将你的数据的统计量与模型下可能产生的值进行比较。

第二步：模型

你要将要测试的模型以数据生成机制的形式呈现出来，以及关于总体的任何特定假设。 这个模型通常包括指定θ ∗，它可能是总体均值、零的比例或回归系数。 在这个模型下统计量的抽样分布称为零分布，而模型本身称为零假设。

第三步：计算

根据第二步中的零模型，你得到的数据（以及产生的统计量）至少与你在第一步中实际得到的数据有多相似？ 在形式推断中，这个概率被称为p -值。 为了近似p -值，我们经常使用计算机使用模型中的假设生成大量重复的随机试验，并找到给定统计量的值至少与我们观察到的值一样极端的样本的比例。 其他时候，我们可以使用数学理论来找到p -值。

第四步：解释

p -值用作惊奇的度量。 如果你在第二步中规定的模型是可信的，那么你得到的数据（和总结统计量）应该有多么惊讶？ 一个中等大小的p -值意味着观察到的统计量基本上是你预期的在零模型生成的数据中得到的。 一个微小的p -值会对零模型提出质疑。 换句话说，如果模型是正确的（或者近似正确的），那么从模型生成的数据中得到这样一个极端值是非常不寻常的。 在这种情况下，要么零模型是错误的，要么发生了一个非常不可能的结果。 统计逻辑表明要得出结论，模式是真实的，不仅仅是巧合。 然后就轮到你解释为什么数据生成过程会导致这样一个不寻常的值。 这是什么时候仔细考虑范围的重要性。

让我们用几个例子来演示测试过程中的这些步骤。

例子：用于比较维基百科贡献者生产力的等级测试

回想一下来自第二章的维基百科示例，在英语维基百科上过去 30 天活跃的前 1%的贡献者中，从未接受过奖励。这 200 名贡献者被随机分成两组各 100 人。一个组的贡献者，即实验组，每人获得一份非正式奖励，而另一个组则没有。所有 200 名贡献者被跟踪记录了 90 天，记录他们在维基百科上的活动。

有人推测非正式奖励对志愿者工作有一种强化作用，这个实验旨在正式研究这一推测。我们基于数据排名进行假设检验。

首先，我们将数据读入数据框架：

`wiki` `=` `pd``.``read_csv``(``"``data/Wikipedia.csv``"``)`
`wiki``.``shape`

(200, 2)

`wiki``.``describe``(``)``[``3``:``]`

	experiment	postproductivity
min	0.0	0.0
25%	0.0	57.5
50%	0.5	250.5
75%	1.0	608.0
max	1.0	2344.0

数据框架有 200 行，每行代表一个贡献者。特征experiment要么是 0 要么是 1，取决于贡献者是在对照组还是实验组，而postproductivity是在授予奖励后 90 天内贡献者编辑次数的计数。四分位数（下、中、上）之间的差距表明生产力的分布是倾斜的。我们制作直方图进行确认：

px.histogram(
	wiki, x='postproductivity', nbins=50,
	labels={'postproductivity': 'Number of actions in 90 days post award'},
	width=350, height=250)

实际上，授奖后的生产力直方图高度倾斜，在接近零的地方有一个峰值。这种偏斜性表明基于两个样本值排序的统计量。

要计算我们的统计量，我们将所有生产力值（来自两个组）从最小到最大排序。最小值的排名为 1，第二小的排名为 2，依此类推，直到最大值，其排名为 200。我们使用这些排名来计算我们的统计量，θ ^ ，它是实验组的平均排名。我们选择这个统计量是因为它对高度倾斜的分布不敏感。例如，无论最大值是 700 还是 700,000，它仍然获得相同的排名，即 200。如果非正式奖励激励贡献者，那么我们期望实验组的平均排名通常高于对照组。

零模型假设非正式奖励对生产力没有任何影响，观察到的实验组和对照组之间的任何差异都是由于将贡献者分配到组的偶然过程。零假设被设置为拒绝现状；也就是说，我们希望找到一个假设没有效果的惊喜。

零假设可以用从一个装有 200 个标有 1、2、3、...、200 的彩球的罐子中抽取 100 次的结果来表示。在这种情况下，平均等级将是 ( 1 + 200 ) / 2 = 100.5 。

我们使用 scipy.stats 中的 rankdata 方法对这 200 个值进行排名，并计算治疗组中等级的总和：

`from` `scipy``.``stats` `import` `rankdata`
`ranks` `=` `rankdata``(``wiki``[``'``postproductivity``'``]``,` `'``average``'``)`

让我们确认 200 个值的平均等级是 100.5：

`np``.``average``(``ranks``)`

100.5

并找出治疗组中 100 个生产力分数的平均等级：

`observed` `=` `np``.``average``(``ranks``[``100``:``]``)`
`observed`

113.68

治疗组的平均等级高于预期，但我们想要弄清楚这是否是一个异常高的值。我们可以使用模拟来找出这个统计量的抽样分布，以确定 113 是一个常规值还是一个令人惊讶的值。

为了进行这个模拟，我们将数据中的 ranks 数组设为罐子中的彩球。对数组中的 200 个值进行洗牌，并取前 100 个值代表一个随机抽取的治疗组。我们编写一个函数来洗牌排名数组并找出前 100 个的平均值。

`rng` `=` `np``.``random``.``default_rng``(``42``)`
`def` `rank_avg``(``ranks``,` `n``)``:`
    `rng``.``shuffle``(``ranks``)`
    `return` `np``.``average``(``ranks``[``n``:``]``)`      

我们的模拟混合了罐子中的彩球，抽取 100 次，计算 100 次抽取的平均等级，然后重复 100,000 次。

`rank_avg_simulation` `=` `[``rank_avg``(``ranks``,` `100``)` `for` `_` `in` `range``(``100_000``)``]` 

这里是模拟平均值的直方图：

正如我们预期的那样，平均等级的抽样分布以 100（实际上是 100.5）为中心，并呈钟形曲线。这个分布的中心反映了治疗没有影响的假设。我们观察到的统计量远远超出了模拟平均等级的典型范围，我们使用这个模拟的抽样分布来找到观察到的统计量的近似 p -值：

`np``.``mean``(``rank_avg_simulation` `>` `observed``)`

0.00058

这真是一个大惊喜。根据零假设，我们看到一个平均等级至少与我们的一样大的机会约为 10,000 分之 5。

这个测试对零模型提出了质疑。统计逻辑使我们得出结论，这种模式是真实存在的。我们如何解释这个结果？这个实验设计得很好。从顶尖 1% 中随机选择了 200 名贡献者，然后随机将他们分为两组。这些随机过程表明，我们可以依赖于这 200 个样本代表了顶尖贡献者，并且治疗组和对照组在除了治疗（奖励）的应用之外的其他方面上都是相似的。考虑到精心的设计，我们得出结论，非正式的奖励对顶尖贡献者的生产力有积极的影响。

早些时候，我们实施了一个模拟来找到我们观察到的统计量的 p -值。在实践中，排名测试通常被广泛使用，并在大多数统计软件中提供：

`from` `scipy``.``stats` `import` `ranksums`

`ranksums``(``x``=``wiki``.``loc``[``wiki``.``experiment` `==` `1``,` `'``postproductivity``'``]``,`
         `y``=``wiki``.``loc``[``wiki``.``experiment` `==` `0``,` `'``postproductivity``'``]``)`

RanksumsResult(statistic=3.220386553232206, pvalue=0.0012801785007519996)

这里的p值是我们计算的p值的两倍，因为我们只考虑大于观察值的值，而ranksums测试计算了分布两侧的p值。在我们的示例中，我们只对提高生产力感兴趣，因此使用单侧p值，这是报告值（0.0006）的一半，接近我们模拟的值。

这种使用排名而不是实际数据值的不太常见的检验统计量是在 1950 年代和 1960 年代开发的，在当今强大笔记本电脑时代之前。排名统计量的数学属性已经很好地发展，并且抽样分布表现良好（即使对于小数据集也是对称的，形状像钟形曲线）。排名测试在 A/B 测试中仍然很受欢迎，其中样本倾向于高度偏斜，并且通常会进行许多许多测试，可以从正态分布中快速计算出p值。

下一个示例重新讨论了来自第三章关于疫苗有效性的例子。在那里，我们遇到了一个假设检验，尽管并未明确称其为假设检验。

示例：疫苗有效性的比例检验

疫苗的批准需满足比我们之前执行的简单测试更严格的要求，其中我们比较了治疗组和对照组中的疾病计数。CDC 要求更强的成功证据，基于每组中患病者比例的比较。为了解释，我们将控制组和治疗组中患病人数的样本比例表示为p ^ C 和 p ^ T，并使用这些比例计算疫苗有效性：

θ ^ = p ^ C − p ^ T p ^ C = 1 − p ^ T p ^ C

在 J&J 试验中观察到的疫苗有效性数值为：

1 − 117 / 21869 351 / 21869 = 1 − 117 351 = 0.667

如果治疗不起作用，有效性将接近于 0。CDC 设定了疫苗有效性的标准为 50%，意味着有效性必须超过 50%才能获得批准用于分发。在这种情况下，零模型假设疫苗有效性为 50%（θ ∗ = 0.5），观察值与预期值的任何差异都归因于随机过程中将人们分配到不同组中。再次强调，我们设定零假设为当前状态，即疫苗不足以获得批准，我们希望发现一个意外并拒绝零假设。

简单代数运算后，零模型 0.5 = 1 − p T / p C 可以简化为 p T = 0.5 p C 。也就是说，零假设暗示接受治疗的人群中患病的比例最多是对照组的一半。请注意，零假设并不假设治疗无效，而是假设其有效性不超过 0.5。

在这种情况下，我们的瓮模型与我们在第三章中设定的有所不同。这个瓮仍然有 43,738 颗弹珠，对应于实验中的受试者。但现在每颗弹珠上有两个数字，为了简单起见，这些数字以一对的形式出现，比如( 0 , 1 ) 。左边的数字是如果人接受治疗后的反应，右边的数字对应于未接受治疗时的反应（对照组）。通常情况下，1 表示他们生病了，0 表示他们保持健康。

零模型假设一对中左边数字的比例是右边数字的一半。由于我们不知道这两个比例，我们可以使用数据来估计它们。瓮中有三种类型的弹珠：( 0 , 0 ) ，( 0 , 1 ) 和( 1 , 1 ) 。我们假设( 1 , 0 ) ，即治疗后患病但对照组未患病的情况是不可能发生的。我们观察到在对照组中有 351 人生病，在治疗组中有 117 人生病。在假设治疗组的患病率是对照组的一半的情况下，我们可以尝试瓮中弹珠的构成情况。例如，我们可以研究这样一种情况，即 117 人在治疗组没有生病，但如果他们在对照组，则所有 585 人（351 + 117 + 117 ）都会感染病毒，其中一半不会感染病毒。表格 17-1 展示了这些计数。

表格 17-1. 疫苗试验瓮模型

标签	计数
(0, 0)	43,152
(0, 1)	293
(1, 0)	0
(1, 1)	293
总数	43,738

我们可以利用这些计数来进行临床试验的模拟，并计算疫苗有效性。如第三章所示，多变量超几何函数模拟从一个装有不止两种颜色的弹珠的罐中抽取弹珠。我们建立这个罐和抽样过程：

`N` `=` `43738`
`n_samp` `=` `21869`
`N_groups` `=` `np``.``array``(``[``293``,` `293``,` `(``N` `-` `586``)``]``)`

`from` `scipy``.``stats` `import` `multivariate_hypergeom`

`def` `vacc_eff``(``N_groups``,` `n_samp``)``:`
    `treat` `=` `multivariate_hypergeom``.``rvs``(``N_groups``,` `n_samp``)`
    `ill_t` `=` `treat``[``1``]`
    `ill_c` `=` `N_groups``[``0``]` `-` `treat``[``0``]` `+` `N_groups``[``1``]` `-` `treat``[``1``]`
    `return` `(``ill_c` `-` `ill_t``)` `/` `ill_c`

现在我们可以对临床试验进行 10 万次模拟，并计算每次试验的疫苗有效性：

`np``.``random``.``seed``(``42``)`
`sim_vacc_eff` `=` `np``.``array``(``[``vacc_eff``(``N_groups``,` `n_samp``)` `for` `_` `in` `range``(``100_000``)``]``)`

`px``.``histogram``(``x``=``sim_vacc_eff``,` `nbins``=``50``,`
            `labels``=``dict``(``x``=``'``Simulated vaccine efficacy``'``)``,`
            `width``=``350``,` `height``=``250``)`

抽样分布的中心在 0.5 处，这与我们的模型假设一致。我们看到 0.667 远离了这个分布的尾部：

`np``.``mean``(``sim_vacc_eff` `>` `0.667``)`

1e-05

只有极少数的 10 万次模拟中，疫苗有效性达到了观察到的 0.667。这是一个罕见事件，这就是为什么 CDC 批准了强生公司的疫苗进行分发。

在这个假设检验的例子中，我们无法完全指定模型，并且必须根据我们观察到的p C和p T的近似值来提供近似值p ^ C和p ^ T 。有时，零模型并没有完全被指定，我们必须依靠数据来建立模型。下一节介绍了一种通用的方法，称为自举法（bootstrap），用于利用数据近似模型。

推断的自举

在许多假设检验中，零假设的假设导致对一个假设总体和数据设计进行完全规范化（参见图 17-1），我们利用这个规范化来模拟统计量的抽样分布。例如，维基百科实验的秩检验导致我们抽样整数 1, …, 200，这很容易模拟。不幸的是，我们并不能总是完全指定总体和模型。为了弥补这种情况，我们用数据代替总体。这种替代是自举概念的核心。图 17-2 更新了图 17-1 以反映这个想法；这里总体分布被经验分布替代，以创建所谓的自举总体。

第 17-2 图。引导数据生成过程的示意图

自举的理由如下：

• 你的样本看起来像是总体，因为它是一个代表性样本，所以我们用样本代替总体，并称之为自举总体。

• 使用产生原始样本的相同数据生成过程来获得一个新样本，这被称为自助样本，以反映人群的变化。以与之前相同的方式计算自助样本上的统计量，并称之为自助统计量。自助统计量的自助抽样分布在形状和传播上应与统计量的真实抽样分布类似。

• 模拟数据生成过程多次，使用自助人群，以获得自助样本及其自助统计量。模拟自助统计量的分布近似于自助统计量的自助抽样分布，后者本身近似于原始抽样分布。

仔细观察图 17-2 并将其与图 17-1 进行比较。基本上，自助模拟涉及两个近似值：原始样本近似于人群，而模拟则近似于抽样分布。到目前为止，我们在例子中一直在使用第二种近似值；样本通过样本来近似人群的近似是自助法的核心概念。注意，在图 17-2 中，自助人群的分布（左侧）看起来像原始样本的直方图；抽样分布（中间）仍然是基于与原始研究中相同的数据生成过程的概率分布，但现在使用的是自助人群；样本分布（右侧）是从自助人群中抽取的一个样本的直方图。

您可能想知道如何从自助人群中简单随机抽取样本，而不是每次都得到完全相同的样本。毕竟，如果您的样本中有 100 个单位，并且将其用作您的自助人群，那么从自助人群中抽取的 100 个单位（不重复）将获取所有单位，并且每次都给出相同的自助样本。解决此问题有两种方法：

• 当从自助人群中抽样时，使用带替换的方式从自助人群中抽取单位。基本上，如果原始人群非常大，则使用替换和不使用替换之间几乎没有区别。这是迄今为止更常见的方法。

• 将样本“扩展”到与原始人群相同的大小。即，计算样本中每个唯一值的比例，并向自助人群中添加单位，使其与原始人群大小相同，同时保持这些比例。例如，如果样本大小为 30，并且样本值的 1/3 是 0，则包含 750 个单位的自助人群应包括 250 个零。一旦有了这个自助人群，就使用原始数据生成过程来进行自助抽样。

疫苗有效性的例子使用了类似 Bootstrap 的过程，称为参数化 Bootstrap。我们的空模型指定了 0-1 瓮，但我们不知道要在瓮中放多少个 0 和 1。我们使用样本来确定 0 和 1 的比例；也就是说，样本指定了多元超几何分布的参数。接下来，我们使用校准空气质量监测器的例子来展示如何使用 Bootstrap 方法来测试假设。

警告

常见的错误是认为 Bootstrap 采样分布的中心与真实采样分布的中心相同。如果样本的均值不为 0，则 Bootstrap 总体的均值也不为 0。这就是为什么在假设检验中我们使用 Bootstrap 分布的传播范围而不是它的中心。下一个例子展示了我们如何使用 Bootstrap 方法来测试假设。

在校准空气质量监测仪器的案例研究中（见第十二章），我们拟合了一个模型来调整廉价监测器的测量值，使其更准确地反映真实的空气质量。这种调整包括与湿度相关的模型项。拟合系数约为0.2，这意味着在高湿度的日子里，测量值的调整幅度比低湿度的日子大。然而，这个系数接近于 0，我们可能会怀疑在模型中是否真的需要包含湿度。换句话说，我们想要检验线性模型中湿度系数是否为 0。不幸的是，我们无法完全指定模型，因为它是基于一组空气监测器（包括 PurpleAir 和 EPA 维护的监测器）在特定时间段内采集的测量数据。这就是 Bootstrap 方法可以帮助的地方。

我们的模型假设所采集的空气质量测量值类似于整体测量值的总体。请注意，天气条件、时间和监测器的位置使得这种说法有些模糊；我们的意思是在原始测量数据采集时，相同条件下的其他测量值类似于这些测量值。此外，由于我们可以想象有一个几乎无限的空气质量测量数据供应，我们认为生成测量数据的过程类似于从瓮中重复有放回地抽取。回顾一下，在第二章中，我们将瓮建模为从测量误差瓮中重复有放回地抽取的过程。这种情况有些不同，因为我们还包括了已经提到的其他因素（天气、季节、位置）。

我们的模型专注于线性模型中湿度系数：

PA ≈ θ 0 + θ 1 AQ + θ 2 RH

在这里，PA指的是 PurpleAir PM2.5 测量，RH是相对湿度，AQ表示更精确的 PM2.5 测量，由更准确的 AQS 监测器进行。零假设是θ 2 = 0；也就是说，零模型是简单模型：

PA ≈ θ 0 + θ 1 AQ

为了估计θ 2，我们使用从第十五章的线性模型拟合过程中获取的结果。

我们的 bootstrap 总体由我们在第十五章中使用的来自乔治亚州的测量值组成。现在我们使用randint的机会机制从数据框中（相当于我们的乌尔恩）有放回地抽样行。这个函数从一组整数中有放回地随机抽取样本。我们使用随机索引样本创建数据框的 bootstrap 样本。然后我们拟合线性模型，并得到湿度系数（我们的 bootstrap 统计量）。以下的boot_stat函数执行这个模拟过程：

`from` `scipy``.``stats` `import` `randint`

`def` `boot_stat``(``X``,` `y``)``:`
    `n` `=` `len``(``X``)`
    `bootstrap_indexes` `=` `randint``.``rvs``(``low``=``0``,` `high``=``(``n` `-` `1``)``,` `size``=``n``)`
    `theta2` `=` `(`
        `LinearRegression``(``)`
        `.``fit``(``X``.``iloc``[``bootstrap_indexes``,` `:``]``,` `y``.``iloc``[``bootstrap_indexes``]``)`
        `.``coef_``[``1``]`
    `)`
    `return` `theta2`

我们设置设计矩阵和结果变量，并检查我们的boot_stat函数一次以测试它：

`X` `=` `GA``[``[``'``pm25aqs``'``,` `'``rh``'``]``]`
`y` `=` `GA``[``'``pm25pa``'``]`

`boot_stat``(``X``,` `y``)`

0.21572251745549495

当我们重复这个过程 10,000 次时，我们得到了对 bootstrap 统计量（拟合湿度系数）的 bootstrap 抽样分布的近似：

`np``.``random``.``seed``(``42``)`
`boot_theta_hat` `=` `np``.``array``(``[``boot_stat``(``X``,` `y``)` `for` `_` `in` `range``(``10_000``)``]``)`

我们对这个 bootstrap 抽样分布的形状和传播感兴趣（我们知道中心将接近原系数0.21）：

`px``.``histogram``(``x``=``boot_theta_hat``,` `nbins``=``50``,`
             `labels``=``dict``(``x``=``'``Bootstrapped humidity coefficient``'``)``,`
             `width``=``350``,` `height``=``250``)`

按设计，bootstrap 抽样分布的中心将接近θ ^，因为 bootstrap 总体由观测数据组成。因此，我们不是计算观察统计量至少大于某个值的机会，而是找到至少小于 0 的值的机会。假设值 0 远离抽样分布。

10,000 个模拟回归系数中没有一个像假设的系数那么小。统计逻辑引导我们拒绝零假设，即我们不需要调整湿度模型。

我们在这里执行的假设检验形式与早期的检验看起来不同，因为统计量的抽样分布不是集中在零假设上。这是因为我们使用 bootstrap 创建抽样分布。实际上，我们使用系数的置信区间来测试假设。在下一节中，我们更广泛地介绍区间估计，包括基于 bootstrap 的区间估计，并连接假设检验和置信区间的概念。

置信区间基础

我们已经看到，建模可以导致估计，例如公交车迟到的典型时间（见第 4 章），空气质量测量的湿度调整（见第 15 章），以及疫苗效力的估计（见第 2 章）。这些例子是未知值的点估计，称为参数：公交车迟到的中位数为 0.74 分钟；空气质量的湿度调整为每湿度百分点 0.21 PM2.5；而疫苗效力中 COVID 感染率的比率为 0.67。然而，不同的样本会产生不同的估计。简单提供点估计并不能反映估计的精度。相反，区间估计可以反映估计的准确性。这些区间通常采用以下两种形式：

• 一个基于 bootstrap 采样分布百分位数创建的自举置信区间。

• 使用抽样分布的标准误差（SE）和关于分布为正态曲线形状的额外假设构建的正态置信区间

我们描述了这两种类型的区间，然后给出了一个例子。回想一下抽样分布（见图 17-1）是反映不同θ ^值观察概率的概率分布。置信区间是根据θ ^的抽样分布扩展构建的，因此区间端点是随机的，因为它们基于θ ^。这些区间被设计成 95%的时间覆盖θ ∗。

正如其名称所示，基于百分位数的 bootstrap 置信区间是从 bootstrap 采样分布的百分位数创建的。具体来说，我们计算θ ^ B的抽样分布的分位数，其中θ ^ B是 bootstrap 统计量。对于 95th 百分位数区间，我们确定 2.5 和 97.5 分位数，分别称为q 2.5 , B和q 97.5 , B，其中 95%的时间，bootstrap 统计量在以下区间内：

q 2.5 , B ≤ θ ^ B ≤ q 97.5 , B

这个 bootstrap 百分位数置信区间被认为是一种快速而粗糙的区间估计方法。有许多替代方法可以调整偏差，考虑分布形状，并且更适合于小样本。

百分位置信区间不依赖于抽样分布具有特定形状或分布中心为θ ∗。相比之下，正态置信区间通常不需要引导抽样来计算，但它确实对θ ^的抽样分布形状做了额外的假设。

我们在抽样分布可以很好地近似于正态曲线时使用正常的置信区间。对于正态概率分布，以中心μ和扩展σ，有 95%的概率，从这个分布中随机取得的值在区间μ ± 1.96 σ内。由于抽样分布的中心通常是θ ∗，对于随机生成的θ ^，有 95%的概率是：

| θ ^ − θ ∗ | ≤ 1.96 S E ( θ ^ )

其中S E ( θ ^ )是θ ^的抽样分布的扩展。我们使用这个不等式来为θ ∗做一个 95%的置信区间：

[ θ ^ − 1.96 S E ( θ ^ ) , θ ^ + 1.96 S E ( θ ^ ) ]

可以使用不同倍数的S E ( θ ^ )基于正态曲线形成其他尺寸的置信区间。例如，99%置信区间为± 2.58 S E，单侧上限 95%置信区间为[ θ ^ − 1.64 S E ( θ ^ ) , ∞ ]。

注意

参数估计的标准差通常称为标准误差或 SE，以区别于样本、总体或来自容器的一次抽样的标准差。在本书中，我们不加区分地称之为标准差。

接下来我们提供每种类型的示例。

在本章的早些时候，我们测试了线性空气质量模型中湿度系数为 0 的假设。这些数据的拟合系数为0.21 。由于空模型未完全指定数据生成机制，我们转而使用自举法。也就是说，我们将数据视为总体，从自举总体中有放回地抽取了 11,226 条记录，并拟合模型以找到湿度的自举样本系数。我们的模拟重复了这一过程 10,000 次，以获得近似的自举抽样分布。

我们可以使用这个自举抽样分布的百分位数来创建θ ∗的 99% 置信区间。为此，我们找到自举抽样分布的分位数q 0.5 和q 99.5：

`q_995` `=` `np``.``percentile``(``boot_theta_hat``,` `99.5``,` `method``=``'``lower``'``)`
`q_005` `=` `np``.``percentile``(``boot_theta_hat``,` `0.05``,` `method``=``'``lower``'``)`

`print``(``f``"``Lower 0.05th percentile:` `{``q_005``:``.3f``}``"``)`
`print``(``f``"``Upper 99.5th percentile:` `{``q_995``:``.3f``}``"``)`

Lower 0.05th percentile: 0.099
Upper 99.5th percentile: 0.260

或者，由于抽样分布的直方图形状大致呈正态分布，我们可以基于正态分布创建一个 99% 置信区间。首先，我们找到θ ^的标准误差，这只是θ ^的抽样分布的标准偏差：

`standard_error` `=` `np``.``std``(``boot_theta_hat``)`
`standard_error`

0.02653498609330345

然后，θ ∗的 99% 置信区间是观察到的θ ^的2.58 个标准误差的距离，向任何方向：

`print``(``f``"``Lower 0.05th endpoint:` `{``theta2_hat` `-` `(``2.58` `*` `standard_error``)``:``.3f``}``"``)`
`print``(``f``"``Upper 99.5th endpoint:` `{``theta2_hat` `+` `(``2.58` `*` `standard_error``)``:``.3f``}``"``)`

Lower 0.05th endpoint: 0.138
Upper 99.5th endpoint: 0.275

这两个区间（自举百分位和正态）非常接近但显然不相同。考虑到自举抽样分布中的轻微不对称性，我们可能会预期到这一点。

还有其他基于正态分布的置信区间版本，反映了使用数据的抽样分布的标准误差的估计的变异性。还有一些针对统计量而不是平均数的百分位数的其他置信区间。（还要注意，对于置换测试，自举法的准确性不如正态近似。）

注意

置信区间很容易被误解为参数θ ∗ 在区间内的概率。然而，置信区间是从抽样分布的一个实现中创建的。抽样分布给出了一个不同的概率陈述；以这种方式构建的区间将在 95% 的时间内包含θ ∗ 。不幸的是，我们不知道这个特定时间是否是那些发生了 100 次中的 95 次之一。这就是为什么使用术语置信度而不是概率或机会，并且我们说我们有 95% 的置信度参数在我们的区间内。

置信区间和假设检验有以下关系。比如说，如果一个 95%的置信区间包含了假设值 θ ∗ ，那么这个检验的 p 值小于 5%。也就是说，我们可以反过来利用置信区间来创建假设检验。我们在前一节中使用了这种技术，当我们进行了对空气质量模型中湿度系数是否为 0 的检验。在本节中，我们基于自举百分位数创建了一个 99%的系数置信区间，由于 0 不在区间内，所以 p 值小于 1%，根据统计逻辑，我们可以得出结论系数不为 0。

另一种区间估计是预测区间。预测区间侧重于观测值的变化而不是估计量的变化。我们接下来来探索这些内容。

预测区间的基础

置信区间传达了估计量的准确性，但有时我们想要对未来观测的预测准确性。例如，有人可能会说：我的公交车三分之一的时间最多晚到三四分之一分钟，但它可能晚多久？另一个例子是，加利福尼亚州鱼类和野生动物部门将达摩尼斯蟹的最小捕捞尺寸设定为 146 毫米，一个休闲钓鱼公司可能会想知道他们顾客的捕捞品是否比 146 毫米更大。还有一个例子，兽医根据驴的长度和腰围估算其重量为 169 公斤，并使用这一估算来给驴注射药物。为了驴的安全，兽医希望知道驴的真实重量可能与这一估算相差多少。

这些例子的共同点是对未来观测的预测以及量化未来观测可能与预测之间的距离。就像置信区间一样，我们计算统计量（估计量），并在预测中使用它，但现在我们关心的是未来观测与预测之间的典型偏差。在接下来的几节中，我们通过基于分位数、标准偏差以及条件协变量的示例来工作，展示预测区间的例子。沿途我们还提供了关于观测值典型变化的额外信息。

示例：预测公交车晚点情况

第四章 模型化了西雅图公交车到达特定站点的晚点情况。我们观察到分布高度倾斜，并选择以中位数 0.74 分钟估算典型的晚点情况。我们在此重现了该章节的样本直方图：

`times` `=` `pd``.``read_csv``(``"``data/seattle_bus_times_NC.csv``"``)`
`fig` `=` `px``.``histogram``(``times``,` `x``=``"``minutes_late``"``,` `width``=``350``,` `height``=``250``)`
`fig``.``update_xaxes``(``range``=``[``-``12``,` `60``]``,` `title_text``=``"``Minutes late``"``)`
`fig`

预测问题解决了公交车可能晚点多久的问题。虽然中位数信息丰富，但它并不提供关于分布偏斜程度的信息。也就是说，我们不知道公交车可能会有多晚。第 75 百分位数，甚至第 95 百分位数，会增加有用的考虑信息。我们在这里计算这些百分位数：

median:          0.74 mins late
75th percentile: 3.78 mins late
95th percentile: 13.02 mins late

从这些统计数据中，我们了解到，超过一半的时间，公交车甚至不晚一分钟，四分之一的时间几乎晚了四分钟，有时几乎可以发生公交车晚了近 15 分钟。这三个值共同帮助我们制定计划。

示例：预测螃蟹大小

对于捕捞遠缅蟹的高度规定，包括限制螃蟹捕捞的贝壳宽度为 146 毫米。为了更好地理解加州渔猎和野生动物部与北加州和南俄勒冈州的商业螃蟹捕捞者合作捕捉、测量和释放螃蟹。这里是约 450 只捕获螃蟹的贝壳大小直方图：

分布略微向左倾斜，但平均值和标准差是分布的合理摘要统计数据：

`crabs``[``'``shell``'``]``.``describe``(``)``[``:``3``]`

count    452.00
mean     131.53
std       11.07
Name: shell, dtype: float64

平均值 132 毫米是典型螃蟹大小的良好预测。然而，它缺乏有关个体螃蟹可能与平均值相差多远的信息。标准差可以填补这一空白。

除了个别观察值围绕分布中心的可变性之外，我们还考虑了我们对平均贝壳大小估计的可变性。我们可以使用自助法来估计这种可变性，或者我们可以使用概率论（我们在下一节中会这样做）来展示估计值的标准差是S D ( p o p ) / n。我们还在下一节中展示，这两种变化来源结合如下：

S D ( p o p ) 2 + S D ( p o p ) 2 n = S D ( p o p ) 1 + 1 n

我们用S D ( s a m p l e )替换S D ( p o p )并将此公式应用于我们的螃蟹：

`np``.``std``(``crabs``[``'``shell``'``]``)` `*` `np``.``sqrt``(``1` `+` `1``/``len``(``crabs``)``)`

11.073329460297957

我们看到包括样本平均值的 SE 基本上不会改变预测误差，因为样本如此之大。我们得出结论，螃蟹通常与 132 毫米的典型大小相差 11 至 22 毫米。这些信息有助于制定围绕螃蟹捕捞的政策，以维护螃蟹种群的健康，并为娱乐捕鱼者设定期望。

示例：预测螃蟹的增长

加州渔业和野生动物部希望更好地理解蟹的生长，以便能够制定更好的捕捞限制，从而保护蟹的种群。在前述例子中提到的研究中捕获的蟹即将脱壳，除了它们的大小外，还记录了换壳前后壳尺寸的变化：

`crabs``.``corr``(``)`

	壳	增量
壳	1.0	-0.6
增量	-0.6	1.0

这两个测量值呈负相关，意味着蟹越大，它们换壳时的生长就越少。我们绘制生长增量与壳尺寸的关系图，以确定这些变量之间的关系是否大致为线性：

px.scatter(crabs, y='inc', x= 'shell', width=350, height=250,
	labels=dict(shell='Dungeness crab shell width (mm)',
				inc='Growth (mm)'))

关系看起来是线性的，我们可以拟合一个简单的线性模型来解释壳前换壳大小的生长增量。在本例中，我们使用statsmodels库，它提供了使用get_prediction生成预测区间的功能。我们首先设置设计矩阵和响应变量，然后使用最小二乘法拟合模型：

`import` `statsmodels``.``api` `as` `sm`

`X` `=` `sm``.``add_constant``(``crabs``[``[``'``shell``'``]``]``)`
`y` `=` `crabs``[``'``inc``'``]`

`inc_model` `=` `sm``.``OLS``(``y``,` `X``)``.``fit``(``)`

`print``(``f``"``Increment estimate =` `{``inc_model``.``params``[``0``]``:``0.2f``}` `+` `"``,` 
      `f``"``{``inc_model``.``params``[``1``]``:``0.2f``}` `x Shell Width``"``)`

Increment estimate = 29.80 +  -0.12 x Shell Width

建模时，我们为解释变量的给定值创建预测区间。例如，如果一只新捕获的蟹的壳宽度为 120 毫米，那么我们使用我们拟合的模型来预测其壳的生长。

如前例所示，我们对单个观测值的预测的可变性包括我们对蟹的生长估计的可变性以及蟹壳尺寸的蟹对蟹变异。我们可以再次使用自助法来估计这种变异，或者可以使用概率理论来展示这两种变异源如下组合：

S D ( e ) 1 + x 0 ( X ⊤ X ) − 1 x 0 ⊤

在这里，X是由原始数据组成的设计矩阵，e是来自回归的残差的n × 1列向量，而x 0是新观测值的1 × ( p + 1 )行特征向量（在本例中，这些是[ 1 , 120 ]）：

`new_data` `=` `dict``(``const``=``1``,` `shell``=``120``)`
`new_X` `=` `pd``.``DataFrame``(``new_data``,` `index``=``[``0``]``)`
`new_X`

	常数	壳
0	1	120

我们使用statsmodels中的get_prediction方法为壳宽度为 120 毫米的蟹找到 95%的预测区间：

`pred` `=` `inc_model``.``get_prediction``(``new_X``)`
`pred``.``summary_frame``(``alpha``=``0.05``)`

	均值	均值标准误	均值置信区间下限	均值置信区间上限	观测置信区间下限	观测置信区间上限
0	15.86	0.12	15.63	16.08	12.48	19.24

在这里，我们有一个对贝壳大小为 120 毫米的螃蟹平均生长增量的置信区间为[15.6, 16.1]，以及生长增量的预测区间为[12.5, 19.2]。预测区间要宽得多，因为它考虑了个体螃蟹的变异性。这种变异性体现在点围绕回归线的分布中，我们通过残差的标准偏差来近似。贝壳大小和生长增量之间的相关性意味着对于特定贝壳大小的生长增量预测的变异性要小于生长增量的整体标准偏差：

`print``(``f``"``Residual SD:` `{``np``.``std``(``inc_model``.``resid``)``:``0.2f``}``"``)`
`print``(``f``"``Crab growth SD:` `{``np``.``std``(``crabs``[``'``inc``'``]``)``:``0.2f``}``"``)`

Residual SD:    1.71
Crab growth SD: 2.14

get_prediction提供的区间依赖于增长增量分布的正态近似。这就是为什么 95% 的预测区间端点大约是预测值两倍的残差标准偏差之外。在下一节中，我们将深入探讨这些标准差、估计值和预测的计算。我们还讨论了在计算它们时所做的一些假设。

推断和预测的概率

假设检验、置信区间和预测区间依赖于从抽样分布和数据生成过程计算的概率计算。这些概率框架还使我们能够对假设调查、实验或其他机会过程进行模拟和自举研究，以研究其随机行为。例如，我们在维基百科实验中找到了排名平均值的抽样分布，假设该实验中的处理方式不有效。使用模拟，我们量化了期望结果的典型偏差和摘要统计量可能值的分布。图 17-1 中的三联画提供了一个图表，指导我们进行这个过程；它有助于区分人群、概率和样本之间的差异，并展示它们之间的联系。在本节中，我们为这些概念带来了更多的数学严谨性。

我们正式介绍了期望值、标准偏差和随机变量的概念，并将它们与本章中用于检验假设、制定置信区间和预测区间的概念联系起来。我们从维基百科实验的具体例子开始，然后进行概括。在此过程中，我们将这种形式主义与贯穿本章的三联画连接起来，作为我们的指导。

平均等级统计理论的形式化

回想一下在 Wikipedia 实验中，我们汇总了治疗组和对照组的奖后生产力值，并将它们转换为排名，1 , 2 , 3 , … , 200 ，因此总体仅由整数 1 到 200 组成。图 17-3 是代表这种特定情况的图表。注意，总体分布是平坦的，范围从 1 到 200（图 17-3 的左侧）。此外，我们使用的总体总结（称为 总体参数）是平均排名：

θ ∗ = Avg ( pop ) = 1 200 Σ k = 1 200 k = 100.5

另一个相关的总结是关于 θ ∗ 的分布，定义为总体标准偏差：

SD ( pop ) = 1 200 Σ k = 1 200 ( k − θ ∗ ) 2 = 1 200 Σ k = 1 200 ( k − 100.5 ) 2 ≈ 57.7

SD(pop) 表示一个排名与整体平均值的典型偏差。为了计算这个例子的 SD(pop)，需要进行一些数学手工操作：

图 17-3. Wikipedia 实验数据生成过程的示意图；这是一个我们知道总体的特殊情况

观察到的样本由治疗组的整数排名组成；我们将这些值称为 k 1 , k 2 , … , k 100. 样本分布显示在 图 17-3 的右侧（100 个整数每个出现一次）。

与总体平均值相对应的是样本平均值，这是我们感兴趣的统计量：

Avg ( sample ) = 1 100 Σ i = 1 100 k i = k ¯ = 113.7

样本 ) 的观察值是 θ ^ 。同样，关于 样本 ) 的分布，被称为样本的标准偏差，表示一个样本中排名与样本平均值的典型偏差：

SD ( sample ) = 1 100 Σ i = 1 100 ( k i − k ¯ ) 2 = 553.

注意在平均值情况下样本统计量和总体参数的定义之间的平行。两个标准偏差之间的平行也值得注意。

接下来我们转向数据生成过程：从乌龙球中抽取 100 个弹珠（其值为 1 , 2 , … , 200 ），无放回地创建治疗排名。我们用大写字母 Z 1 表示从乌龙球中抽取第一个弹珠以及我们得到的整数。这个 Z 1 被称为 随机变量。它有一个由乌龙球模型确定的概率分布。也就是说，我们可以列出 Z 1 可能取到的所有值以及每个值对应的概率：

P ( Z 1 = k ) = 1 200  for  k = 1 , … , 200

在本例中，Z 1的概率分布由一个简单的公式确定，因为从罐子中抽取的所有整数都有相等的可能性。

我们经常通过期望值和标准差来总结随机变量的分布。就像在人群和样本中一样，这两个量给了我们预期结果和实际值可能偏离预期的程度的感觉。

对于我们的例子，Z 1的期望值简单地是：

E [ Z 1 ] = 1 P ( Z 1 = 1 ) + 2 P ( Z 1 = 2 ) + ⋯ + 200 P ( Z 1 = 200 ) = 1 × 1 200 + 2 × 1 200 + ⋯ + 200 × 1 200 = 100.5

注意，E [ Z 1 ] = θ ∗ ，即来自于罐子的总体平均值。在人群中的平均值以及代表从包含该人群的罐子中随机抽取的随机变量的期望值始终相同。通过将人群平均值表达为人群中唯一值的加权平均值，权重为具有该值的单位比例，这一点更容易理解。从人群罐子中随机抽取的随机变量的期望值使用完全相同的权重，因为它们匹配选择特定值的几率。

注意

术语期望值可能会有点令人困惑，因为它不一定是随机变量的可能值。例如，E [ Z 1 ] = 100.5，但Z 1只可能是整数。

接下来，Z 1的方差定义如下：

V ( Z 1 ) = E [ Z 1 − E ( Z 1 ) ] 2 = [ 1 − E ( Z 1 ) ] 2 P ( Z 1 = 1 ) + ⋯ + [ 200 − E ( Z 1 ) ] 2 P ( Z 1 = 200 ) = ( 1 − 100.5 ) 2 × 1 200 + ⋯ + ( 200 − 100.5 ) 2 × 1 200 = 3333.25

此外，我们定义Z 1的标准差如下所示：

SD ( Z 1 ) = V ( Z 1 ) = 57.7

再次指出Z 1的标准差与SD（pop）相匹配。

为了描述图 17-3 中的整个数据生成过程，我们还定义Z 2，Z 3，…，Z 100作为从罐子中进行的其余 99 次抽样的结果。基于对称性，这些随机变量应该都具有相同的概率分布。也就是说，对于任意的k = 1，…，200：

P ( Z 1 = k ) = P ( Z 2 = k ) = ⋯ = P ( Z 100 = k ) = 1 200

这意味着每个Z i的期望值和标准差均相同，分别为 100.5 和 57.7。然而，这些随机变量并不独立。例如，如果你知道Z 1 = 17，那么Z 2 = 17 是不可能的。

要完成图 17-3 的中间部分，涉及θ ^的抽样分布，我们将平均秩统计量表示如下：

θ ^ = 1 100 Σ i = 1 100 Z i

我们可以利用Z 1的期望值和标准差以及我们对数据生成过程的了解，找到θ ^的期望值和标准差。首先我们找到θ ^的期望值：

E ( θ ^ ) = E [ 1 100 Σ i = 1 100 Z i ] = 1 100 Σ i = 1 100 E [ Z i ] = 100.5 = θ ∗

换句话说，从总体中随机抽取的平均值的期望值等于总体平均值。这里我们提供了关于总体方差和平均值方差的公式，以及标准差：

V ( θ ^ ) = V [ 1 100 Σ i = 1 100 Z i ] = 200 − 100 100 − 1 × V ( Z i ) 100 = 16.75 SD ( θ ^ ) = 100 199 SD ( Z 1 ) 10 = 4.1

这些计算依赖于随机变量的期望值和方差以及随机变量和随机变量和的若干性质。接下来，我们提供用于推导刚才提出的公式的随机变量和随机变量平均数的性质。

随机变量的一般性质

一般来说，随机变量表示机会事件的数值结果。在本书中，我们使用像X、Y或Z这样的大写字母来表示随机变量。X的概率分布是指定的P ( X = x ) = p x，适用于随机变量可能取的所有值x。

因此，X的期望值定义如下：

E [ X ] = ∑ x x p x

方差X的定义如下：

V ( X ) = E [ ( X − E [ X ] ) 2 ] = ∑ x [ x − E ( X ) ] 2 p x

而SD ( X )是V ( X )的平方根。

注意

尽管随机变量可以表示数量，这些数量可以是离散的（例如从总体中随机抽取的家庭中孩子的数量）或连续的（例如由空气监测仪器测量的空气质量），但我们在本书中只讨论具有离散结果的随机变量。由于大多数测量都具有一定程度的精确性，这种简化并不会太大限制我们。

简单的公式提供了期望值、方差和标准差，当我们对随机变量进行比例和移位变化时，比如对于常数a和b，如a + b X：

E ( a + b X ) = a + b E ( X ) V ( a + b X ) = b 2 V ( X ) S D ( a + b X ) = | b | S D ( X )

要确信这些公式是合理的，请考虑如果向每个值添加常数a或将每个值按比例b缩放，分布如何变化。向每个值添加a只会简单地移动分布，从而移动期望值但不改变偏离期望值的大小。另一方面，例如将值按 2 倍缩放会扩展分布，并且基本上会使期望值和偏离期望值的大小加倍。

我们还对两个或更多随机变量的和的属性感兴趣。让我们考虑两个随机变量，X和Y。那么：

E ( a + b X + c Y ) = a + b E ( X ) + c E ( Y )

但是要找到a + b X + c Y的方差，我们需要知道X和Y如何一起变化，这称为X和Y的联合分布。X和Y的联合分布将概率分配给它们结果的组合：

P ( X = x , Y = y ) = p x , y

描述X和Y如何一起变化的总结，称为协方差，定义如下：

C o v ( X , Y ) = E [ ( X − E [ X ] ) ( Y − E [ Y ] ) ] = E [ ( X Y ) − E ( X ) E ( Y ) ] = Σ x , y [ ( x y ) − E ( X ) E ( Y ) ] p x , y

协方差进入到( a + b X + c Y )的计算中，如下所示：

V ( a + b X + c Y ) = b 2 V ( X ) + 2 b c C o v ( X , Y ) + c 2 V ( Y )

在X和Y独立的特殊情况下，它们的联合分布简化为p x , y = p x p y。在这种情况下，C o v ( X , Y ) = 0，因此：

V ( a + b X + c Y ) = b 2 V ( X ) + c 2 V ( Y )

这些性质可以用来表明，对于独立的随机变量X 1 , X 2 , … , X n，其期望值为μ，标准差为σ，平均值X ¯的期望值、方差和标准差如下：

E ( X ¯ ) = μ V ( X ¯ ) = σ 2 / n S D ( X ¯ ) = σ / n

这种情况发生在乌恩模型中，X 1 , … , X n是带替换随机抽取的结果。在这种情况下，μ代表乌恩的平均值，σ代表标准差。

然而，当我们从乌恩中无替换随机抽取时，X i不是独立的。在这种情况下，X ¯具有以下期望值和方差：

E ( X ¯ ) = μ V ( X ¯ ) = N − n N − 1 × σ 2 n

请注意，虽然期望值与无替换抽取时相同，但方差和标准差较小。这些量由称为有限总体修正因子的公式调整，公式为( N − n ) / ( N − 1 )。我们先前在维基百科示例中使用了这个公式来计算S D ( θ ^ )。

返回到图 17-3，我们看到在图表中央的X ¯的抽样分布具有与总体平均值匹配的期望；标准差像1 / n减小，但更快，因为我们是无替换抽取；并且分布形状类似正态曲线。我们在模拟研究中早些时候已经看到了这些特性。

现在我们已经概述了随机变量及其总和的一般性质，我们将这些想法与测试、置信区间和预测区间联系起来。

测试和区间背后的概率

正如本章开头提到的那样，概率是进行假设检验、为估计器提供置信区间和未来观测的预测区间的基础。

现在我们有了技术装备来解释这些概念，在本章中我们已经仔细定义了这些概念，没有使用正式的技术术语。这一次，我们用随机变量及其分布来展示结果。

请记住，假设检验依赖于一个提供统计量θ ^的概率分布的空模型。我们进行的测试本质上是计算（有时是近似地）以下概率。鉴于空分布的假设：

P ( θ ^ ≥ observed statistic )

通常，随机变量被标准化以使计算更容易和标准化：

P ( θ ^ − θ ∗ S D ( θ ^ ) ≥ observed stat − θ ∗ S D ( θ ^ ) )

当S D ( θ ^ )未知时，我们通过模拟来近似它，并且当我们有关于S D ( θ ^ )的公式，以S D ( p o p )为替代。这种标准化很流行，因为它简化了零分布。例如，如果θ ^近似服从正态分布，那么标准化版本将具有以 0 为中心和标准差 1 的标准正态分布。这些近似在进行大量假设检验时非常有用，比如 A/B 测试，因为不需要为每个统计量模拟，只需使用正态曲线的概率即可。

信心区间背后的概率陈述与测试中使用的概率计算非常相似。特别是，为了创建一个 95%的信心区间，其中估计量的抽样分布大致为正态分布，我们进行标准化并使用概率：

P ( | θ ^ − θ ∗ | S D ( θ ^ ) ≤ 1.96 ) = P ( θ ^ − 1.96 S D ( θ ^ ) ≤ θ ∗ ≤ θ ^ + 1.96 S D ( θ ^ ) ) ≈ 0.95

请注意，在前述概率陈述中θ ^是一个随机变量，而θ ∗被认为是一个固定的未知参数值。通过将观察到的统计量替换为θ ^来创建 95%的信心区间：

[ observed stat − 1.96 S D ( θ ^ ) , observed stat + 1.96 S D ( θ ^ ) ]

一旦观察到的统计量替换了随机变量，那么我们说我们有 95%的信心，我们创建的区间包含真值θ ∗。换句话说，在这样计算区间的 100 个案例中，我们期望 95 个覆盖我们正在估计的总体参数。

接下来，我们考虑预测区间。基本概念是提供一个区间，表示未来观察的期望变化约估器。在简单情况下，统计量为X¯，我们有一个假设的新观测值X0，其期望值相同，即μ，和标准差，即σ，与Xi的相同。然后我们找到平方损失的期望变化：

E [ ( X 0 − X ¯ ) 2 ] = E { [ ( X 0 − μ ) − ( X ¯ − μ ) ] 2 } = V ( X 0 ) + V ( X ¯ ) = σ 2 + σ 2 / n = σ 1 + 1 / n

注意，变化有两部分：一部分是由于X0的变化，另一部分是由于E (X0)由X¯逼近。

对于更复杂的模型，预测的变化也分解为两个组成部分：关于模型的数据固有变化以及由于模型估计的抽样分布的变化。假设模型大致正确，我们可以表达如下：

Y = X θ ∗ + ϵ

其中θ ∗是一个(p+1)×1列向量，X是一个n×(p+1)设计矩阵，ϵ由n个独立的随机变量组成，每个变量的期望值为 0，方差为σ2。在这个方程中，Y是一个随机变量向量，每个变量的期望值由设计矩阵确定，方差为σ2。也就是说，关于直线的变化是恒定的，即它不随x的变化而变化。

当我们在回归中创建预测区间时，它们会被赋予一个协变量的 1 × ( p + 1 ) 行向量，称为 x 0 。然后预测是 x 0 θ ^ ，其中 θ ^ 是基于原始 y 和设计矩阵 X 的估计参数向量。这种预测的预期平方误差是：

E [ ( Y 0 − x 0 θ ^ ) 2 ] = E { [ ( Y 0 − x 0 θ ∗ ) − ( x 0 θ ^ − x 0 θ ∗ ) ] 2 } = V ( ϵ 0 ) + V ( x 0 θ ^ ) = σ 2 [ 1 + x 0 ( X ⊤ X ) − 1 x 0 ⊤ ]

我们用最小二乘拟合的残差的方差来近似ϵ 的方差。

我们使用正态曲线创建的预测区间还依赖于另一个假设，即错误的分布近似正态。这是一个比我们为置信区间所做的假设更强的假设。对于置信区间，X i 的概率分布不需要看起来正态，X ¯ 的概率分布就会近似正态。同样，在线性模型中，ϵ 的概率分布不需要看起来正态，估计参数 θ ^ 就会近似正态。

当我们制作这些预测区间时，我们还假设线性模型近似正确。在第十六章中，我们考虑了拟合模型不匹配产生数据的模型的情况。我们现在有技术机器可以推导出该章节中介绍的模型偏差-方差权衡。它与带有一些小变化的预测区间推导非常相似。

模型选择背后的概率

在第十六章中，我们用均方误差（MSE）介绍了模型欠拟合和过拟合。我们描述了一般设置，其中数据可能表示为如下所示：

y = g ( x ) + ϵ

假设ϵ 表现为没有趋势或模式、方差恒定且彼此独立的随机误差。模型中的 信号 是函数 g ( ) 。数据是 ( x i , y i ) 对，我们通过最小化均方误差来拟合模型：

min f ∈ F 1 n ∑ i = 1 n ( y i − f ( x i ) 2

这里F是我们正在最小化的模型集合。这个集合可能是所有$m$次或更低次多项式，有一个在点$k$处弯曲的曲线，等等。注意$g$不一定要在我们用来拟合模型的函数集合中。

我们在模型选择中的目标是选择一个能够很好地预测新观测的模型。对于一个新的观测，我们希望期望损失很小：

E [ y 0 − f ( x 0 ) ] 2

这种期望是关于可能的分布( x 0 , y 0 ) 的，并被称为风险。由于我们不知道( x 0 , y 0 ) 的总体分布，因此我们无法计算风险，但我们可以通过收集到的数据上的平均损失来近似计算：

E [ y 0 − f ( x 0 ) ] 2 ≈ 1 n ∑ i = 1 n ( y i − f ( x i ) ) 2

这个近似被称为经验风险。但希望你能认出它是 MSE：

我们通过最小化所有可能模型上的经验风险（或 MSE）来拟合模型，F = { f } ：

min f ∈ F 1 n ∑ i = 1 n ( y i − f ( x i ) ) 2

拟合的模型称为f ^，它是线性模型X θ ^ 的稍微更一般的表示。这个技术被称为经验风险最小化。

在第十六章中，我们看到当我们使用经验风险来既拟合模型又评估新观测的风险时出现问题。理想情况下，我们希望估计风险（期望损失）：

E [ ( y 0 − f ^ ( x 0 ) ) 2 ]

期望值是关于新观测( x 0 , y 0 ) 和f ^（其中包括原始数据( x i , y i )，$i=1，\ldots，n$）。

为了理解问题，我们将这种风险分解为代表模型偏差、模型方差和来自$\epsilon$的不可约误差的三个部分：

E [ y 0 − f ^ ( x 0 ) ] 2 = E [ g ( x 0 ) + ϵ 0 − f ^ ( x 0 ) ] 2 definition of y 0 = E [ g ( x 0 ) + ϵ 0 − E [ f ^ ( x 0 ) ] + E [ f ^ ( x 0 ) ] − f ^ ( x 0 ) ] 2 adding ± E [ f ^ ( x 0 ) ] = E [ g ( x 0 ) − E [ f ^ ( x 0 ) ] − ( f ^ ( x 0 ) − E [ f ^ ( x 0 ) ] ) + ϵ 0 ] 2 rearranging terms = [ g ( x 0 ) − E [ f ^ ( x 0 ) ] ] 2 + E [ f ^ ( x 0 ) − E [ f ^ ( x 0 ) ] ] 2 + σ 2 expanding the square = model bias 2 + model variance + error

要推导标记为“展开平方”的等式，我们需要正式证明展开中的交叉项都是 0。这需要一些代数，我们在此不详细介绍。但主要思想是术语ϵ 0和( f ^ ( x 0 ) − E [ f ^ ( x 0 ] )是独立的，并且两者的期望值都为 0。最终方程中的其余三项——模型偏差、模型方差和不可减少的错误——如下所述：

模型偏差

最终方程中的第一个术语是模型偏差（平方）。当信号g不属于模型空间时，我们有模型偏差。如果模型空间能很好地逼近g，那么偏差就很小。请注意，我们的预测区间中不存在（或只有极少量的）模型偏差。

模型方差

第二项代表来自数据的拟合模型的变异性。在早期的例子中，我们看到高阶多项式可能会过拟合，因此在不同数据集之间变化很大。模型空间越复杂，拟合模型的变异性越大。

不可减少的错误

最后，最后一个术语是错误的可变性，即ϵ 0，也称为“不可减少的错误”。这种错误会保留下来，无论我们是使用简单模型（高偏差）还是复杂模型（高方差）进行欠拟合或过拟合。

这种预期损失的表示显示了拟合模型的偏差-方差分解。模型选择旨在平衡这两种竞争性误差来源。在第十六章介绍的训练-测试分割、交叉验证和正则化是模仿新观察预期损失或惩罚模型过拟合的技术。

虽然本章涵盖了大量的理论内容，我们尝试将其与乌尔姆模型的基础以及三种分布（总体、样本和抽样）联系起来。在进行假设检验和制作置信度或预测区间时，请牢记以下几点。

总结

本章的整个过程基于乌尔姆模型来发展推断和预测的理论。乌尔姆模型引入了估计量的概率分布，例如样本均值和最小二乘回归系数。我们在此结束本章，并提出一些关于这些统计程序的注意事项。

我们看到，估计量的标准差在分母中有样本大小的平方根因子。当样本很大时，标准差可能会很小，从而导致拒绝假设或非常窄的置信区间。在这种情况下，考虑以下几点是很有必要的：

• 您检测到的差异是否重要？也就是说，

  -值可能会很小，表明一个令人惊讶的结果，但观察到的实际效果可能不重要。统计显著性并不意味着实际显著性。

• 请记住，这些计算不包括偏差，如非响应偏差和测量偏差。这种偏差可能比由抽样分布中的偶然变化引起的任何差异都要大。

有时，我们知道样本不是来自于随机机制，但进行假设检验仍然有用。在这种情况下，零模型将检验样本（和估计量）是否像是随机的。当这个检验被拒绝时，我们确认某些非随机因素导致了观察到的数据。这可能是一个有用的结论：我们观察到的差异不能仅仅通过偶然性来解释。

在其他时候，样本包含整个人群。这种情况下，我们可能不需要进行置信区间或假设检验，因为我们已经观察到了人群中的所有值。也就是说，不需要进行推断。但是，我们可以对假设检验进行不同的解释：我们可以假设观察到的两个特征之间的任何关系是随机分布的，而彼此之间没有关联。

我们还看到了当我们没有足够的关于人群的信息时，如何利用自助法。自助法是一种强大的技术，但它也有局限性：

• 确保原始样本是大且随机的，以使样本类似于人群。

• 重复自助法过程多次。通常进行 10,000 次重复是一个合理的数量。

• 自助法在以下情况下往往会遇到困难：

  • 估计量受异常值的影响。

  • 参数基于分布的极端值。

  • 统计量的抽样分布远非钟形。

或者，我们依赖于抽样分布大致呈正态形状的情况。有时，抽样分布看起来大致正态，但尾部较厚。在这些情况下，使用-分布家族可能比正态分布更合适。

模型通常只是底层现实的一种近似，而 θ ∗ 精确等于 0 的说法与这种模型的概念相左。推断取决于我们模型的正确性。我们可以部分地检查模型的假设，但任何模型都会存在一定的怀疑。事实上，往往数据表明可能存在多个可能的模型，这些模型甚至可能是矛盾的。

最后，有时候假设检验或置信区间的数量可能会相当大，我们需要小心避免产生虚假结果。这个问题被称为 p -hacking，是科学中描述的可重复性危机的另一个例子，详见第十章。 P -hacking 基于这样的观点，即如果我们测试，比如说，100 个假设，所有这些假设都是真实的，那么我们预计会得到一些意外的结果并拒绝其中一些假设。这种现象在多元线性回归中可能会发生，特别是在模型中具有大量特征时，已经开发出技术来限制这些虚假发现的危险。

我们接下来通过一个案例研究来总结建模过程。