十六、交叉验证和正则化

Cross Validation and Regularization

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 认识到需要验证和测试集来预览模型在未知数据上的表现

• 应用交叉验证来选择模型超参数

• 了解 L1 和 L2 正则化的概念基础

在特征工程讲座结束时（第 14 讲），我们提出了调整模型复杂度的问题。我们发现一个过于复杂的模型会导致过拟合，而一个过于简单的模型会导致欠拟合。这带来了一个自然的问题：我们如何控制模型复杂度以避免欠拟合和过拟合？

为了回答这个问题，我们需要解决两件事：首先，我们需要通过评估模型在未知数据上的表现来了解何时我们的模型开始过拟合。我们可以通过交叉验证来实现这一点。其次，我们需要引入一种调整模型复杂度的技术 - 为此，我们将应用正则化。

16.1 训练、测试和验证集

从上一讲中，我们了解到增加模型复杂度减少了模型的训练误差，但增加了它的方差。这是很直观的：添加更多的特征使我们的模型更紧密地拟合了训练过程中遇到的数据，但对新数据的泛化能力更差。因此，低训练误差并不总是代表我们模型的基本性能 - 我们还需要评估它在未知数据上的表现，以确保它没有过拟合。

事实上，唯一知道我们的模型何时过拟合的方法是在未知数据上评估它。不幸的是，这意味着我们需要等待更多的数据。这可能非常昂贵和耗时。

我们应该如何进行？在本节中，我们将建立一个可行的解决方案来解决这个问题。

16.1.1 测试集

避免过拟合的最简单方法是将一些数据“保密”不让自己知道。我们可以将完整数据集的随机部分保留下来，仅用于测试目的。这个测试集中的数据点不会用于模型拟合过程。相反，我们将：

• 使用我们数据集的剩余部分 - 现在称为训练集 - 运行普通最小二乘法、梯度下降或其他一些技术来拟合模型参数

• 拿到拟合的模型并用它对测试集中的数据点进行预测。模型在测试集上的表现（以 MSE、RMSE 等表示）现在表明了它在未知数据上的预测能力有多好

重要的是，最佳模型参数是通过仅考虑训练集中的数据找到的。在模型拟合到训练数据之后，我们在进行测试集上的预测之前不改变任何参数。重要的是，我们在最终确定所有模型设计后，只对测试集进行一次预测。我们将测试集的表现视为模型表现的最终测试。

将我们的数据集分成训练集和测试集的过程被称为训练-测试拆分。通常，10%到 20%的数据被分配给测试集。

在sklearn中，model_selection模块的train_test_split函数允许我们自动生成训练-测试拆分。

在今天的工作中，我们将继续使用之前讲座中的vehicles数据集。与以往一样，我们将尝试从hp的转换中预测车辆的mpg。在下面的单元格中，我们将完整数据集的 20%分配给测试，剩下的 80%分配给训练。

代码

import pandas as pd
import numpy as np
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

# Load the dataset and construct the design matrix
vehicles = sns.load_dataset("mpg").rename(columns={"horsepower":"hp"}).dropna()
X = vehicles[["hp"]]
X["hp^2"] = vehicles["hp"]**2
X["hp^3"] = vehicles["hp"]**3
X["hp^4"] = vehicles["hp"]**4

Y = vehicles["mpg"]

from sklearn.model_selection import train_test_split

# `test_size` specifies the proportion of the full dataset that should be allocated to testing
# `random_state` makes our results reproducible for educational purposes
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=220)

print(f"Size of full dataset: {X.shape[0]} points")
print(f"Size of training set: {X_train.shape[0]} points")
print(f"Size of test set: {X_test.shape[0]} points")

Size of full dataset: 392 points
Size of training set: 313 points
Size of test set: 79 points

在进行训练-测试拆分后，我们将模型拟合到训练集，并评估其在测试集上的表现。

import sklearn.linear_model as lm
from sklearn.metrics import mean_squared_error

model = lm.LinearRegression()

# Fit to the training set
model.fit(X_train, Y_train)

# Make predictions on the test set
test_predictions = model.predict(X_test)

16.1.2 验证集

现在，如果我们对测试集的性能不满意怎么办？按照我们目前的框架，我们就会陷入困境。如前所述，评估模型在测试集上的性能是模型设计过程的最终阶段。我们不能回头根据新发现的过拟合来调整模型 - 如果这样做，我们就会考虑测试集的信息来设计我们的模型。测试误差将不再是模型在未见数据上性能的真实代表！

我们的解决方案是引入一个验证集。验证集是训练集的一个随机部分，用于在模型仍在开发中时评估模型性能。使用验证集的过程是：

• 进行训练-测试分割。将测试集放在一边；直到模型设计过程的最后才会使用。

• 设置一部分训练集用于验证。

• 将模型参数拟合到剩余训练集中包含的数据点。

• 评估模型在验证集上的性能。根据需要调整模型，重新拟合剩余部分的训练集，然后在验证集上重新评估。如有必要，重复此过程直到满意为止。

• 在所有模型开发完成后，评估模型在测试集上的性能。这是模型在未见数据上表现如何的最终测试。不应对模型进行进一步修改。

创建验证集的过程称为验证分割。

请注意，验证误差与之前探讨的训练误差行为有很大不同。回想一下，训练误差随着模型度数的增加而单调递减 - 随着模型变得更复杂，它在训练数据上做出了更好的预测。相反，验证误差在增加模型复杂度时先减少然后增加。这反映了从欠拟合到过拟合的转变。在低模型复杂度时，模型欠拟合，因为它不够复杂以捕捉数据的主要趋势。在高模型复杂度时，模型过拟合，因为它对训练数据进行了过于紧密的“记忆”。

我们可以更新我们对误差、复杂度和模型方差之间关系的理解：

我们的目标是训练一个复杂度接近橙色虚线的模型 - 这是我们的模型在验证集上达到最小误差的地方。请注意，这种关系是对现实世界的简化。但对于 Data 100 来说，这已经足够了。

16.2 K-Fold Cross-Validation

引入验证集给了我们一个“额外”的机会，评估模型在另一组未见数据上的性能。我们能够根据模型在这一组验证数据上的性能来微调模型设计。

但是，如果碰巧我们的验证集包含了很多异常值怎么办？可能我们设置的验证数据点实际上并不代表模型可能遇到的其他未见数据。理想情况下，我们希望在几个不同的未见数据集上验证模型的性能。这将让我们更加自信地了解模型在新数据上的行为。

让我们回顾一下我们的验证框架。之前，我们设置了训练数据的 x%（比如 20%）用于验证。

在上面的例子中，我们设置了前 20%的训练数据点作为验证集。这是一个任意的选择。我们可以将任何20%的训练数据用于验证。事实上，有 5 个不重叠的“块”训练数据点，我们可以指定为验证集。

其中一个这样的块的常见术语是折叠。在上面的示例中，我们有 5 个折叠，每个折叠包含 20%的训练数据。这给了我们一个新的视角：我们在我们的训练集中实际上有5个“隐藏”的验证集。

在交叉验证中，我们为训练集中的每个折叠执行验证拆分。对于具有$K$个折叠的数据集，我们：

• 选择一个折叠作为验证折叠

• 将模型拟合到除验证折叠之外的每个折叠的训练数据

• 计算验证折叠上的模型误差并记录它

• 对所有$K$个折叠重复

交叉验证误差然后是所有$K$个验证折叠的平均误差。

16.2.1 模型选择工作流程

在这个阶段，我们已经完善了我们的模型选择工作流程。我们首先执行训练-测试拆分，以设置一个测试集，用于最终评估模型性能。然后，我们在调整设计矩阵和计算交叉验证误差之间交替，以微调模型的设计。在下面的示例中，我们说明了使用 4 折交叉验证来帮助确定模型设计。

16.2.2 超参数

交叉验证的一个重要用途是进行超参数选择。超参数是模型中在模型拟合到任何数据之前选择的一些值。这意味着它与模型参数$\theta_i$不同，因为它的值是在训练过程开始之前选择的。我们不能使用我们通常的技术 - 微积分、普通最小二乘法或梯度下降 - 来选择它的值。相反，我们必须自己决定。

Data 100 中一些超参数的例子是：

• 我们的多项式模型的程度（回想一下，在创建设计矩阵和调用.fit之前我们选择了程度）

• 梯度下降中的学习率$\alpha$

• 正则化惩罚$\lambda$（稍后将介绍）

为了通过交叉验证选择超参数值，我们首先列出了几个关于最佳超参数可能是什么的“猜测”。对于每个猜测，我们然后运行交叉验证，计算模型在使用该超参数值时产生的交叉验证误差。然后我们选择导致最低交叉验证误差的超参数值。

例如，我们可能希望使用交叉验证来决定我们应该使用什么值作为$\alpha$，它控制每次梯度下降更新的步长。为此，我们列出了最佳$\alpha$的一些可能猜测：0.1、1 和 10。对于每个可能的值，我们执行交叉验证，看当我们使用该$\alpha$值来训练模型时，模型产生了什么错误。

16.3 正则化

我们现在已经解决了今天的两个目标中的第一个：创建一个框架来评估模型在未见数据上的性能。现在，我们将讨论我们的第二个目标：开发一种调整模型复杂性的技术。这将使我们能够直接解决欠拟合和过拟合的问题。

早些时候，我们通过调整超参数（多项式的程度）来调整多项式模型的复杂性。我们尝试了几个不同的多项式程度，计算了每个的验证误差，并选择了最小化验证误差的值。调整“复杂性”很简单；只需要调整多项式程度。

在大多数机器学习问题中，复杂性的定义与我们迄今为止所见的不同。今天，我们将探讨复杂性的两种不同定义：$\theta_i$系数的平方和绝对大小。

16.3.1 约束模型参数

回想一下我们使用梯度下降来下降损失曲面的工作。您可能会发现参考梯度下降笔记来提醒自己会很有帮助。我们的目标是找到导致模型损失最小的模型参数组合。我们通过在水平和垂直轴上绘制可能的参数值来使用等高线图来可视化这一点，这使我们可以从上方鸟瞰损失曲面。我们希望找到对应于损失曲面上最低点的模型参数。

让我们回顾一下我们当前的建模框架。

$\hat{\mathbb{Y}} = \theta_0 + \theta_1 \phi_1 + \theta_2 \phi_2 + \ldots + \theta_p \phi_p$

回想一下，我们用$\phi_i$表示我们的特征，以反映我们进行了特征工程的事实。

以前，我们通过限制模型中存在的特征的总数来限制模型的复杂性。我们一次只包括有限数量的多项式特征；所有其他多项式都被排除在模型之外。

如果我们不是完全删除特定特征，而是保留所有特征，并且每个特征只使用“一点点”，会怎么样？如果我们限制每个特征对预测的贡献量，我们仍然可以控制模型的复杂性，而无需手动确定应该删除多少个特征。

我们所说的“一点点”是什么意思？考虑某个参数$\theta_i$接近或等于 0 的情况。那么，特征$\phi_i$几乎不会影响预测 - 特征的权重值如此之小，以至于它的存在并不会显著改变\hat{\mathbb{Y}的值。如果我们限制每个参数$\theta_i$的大小，我们就限制了特征$\phi_i$对模型的贡献。这会减少模型的复杂性。

在正则化中，我们通过对模型参数$\theta_i$的大小设置限制来限制模型的复杂性。

这些限制看起来是什么样子？假设我们规定所有绝对参数值的总和不能大于某个数字$Q$。换句话说：

$\sum_{i=1}^p |\theta_i| \leq Q$

其中$p$是模型中参数的总数。您可以将这看作是我们为模型分配每个参数的大小的“预算”。如果模型为某些$\theta_i$分配了一个较大的值，它可能必须为其他一些$\theta_j$分配一个较小的值。这会增加特征$\phi_i$对预测的影响，同时减少特征$\phi_j$的影响。模型需要战略地分配参数权重 - 理想情况下，更“重要”的特征将获得更大的权重。

请注意，截距项$\theta_0$不受此约束的影响。我们通常不对截距项进行正则化。

现在，让我们回想一下梯度下降，并将损失曲面可视化为等高线图。损失曲面表示每个点代表模型对$\theta_1$、$\theta_2$的特定组合的损失。假设我们的目标是找到使我们获得最低损失的参数组合。

在没有约束的情况下，最优的$\hat{\theta}$位于中心。

应用这个约束限制了模型参数的有效组合。现在我们只能考虑总绝对和小于或等于我们的数字$Q$的参数组合。这意味着我们只能将我们的正则化参数向量$\hat{\theta}_{\text{Reg}}$分配到下面绿色菱形中的位置。

我们不能再选择真正最小化损失曲面的参数向量$\hat{\theta}_{\text{No Reg}}$，因为这组参数不在我们允许的区域内。相反，我们选择任何允许的组合，使我们尽可能接近真正的最小损失。

请注意，在正则化下，我们优化的$\theta_1$和$\theta_2$的值要比没有正则化时小得多（确实，$\theta_1$已经减少到 0）。模型的复杂度降低，因为我们限制了特征对模型的贡献。事实上，通过将其参数设置为 0，我们有效地从模型中完全删除了特征$\phi_1$的影响。

如果我们改变$Q$的值，我们就改变了允许的参数组合区域。模型仍然会选择产生最低损失的参数组合——最接近受约束区域中真正的最小化器$\hat{\theta}_{\text{No Reg}}$的点。

如果我们使$Q$更小：如果我们使$Q$更大：

• 当$Q$很小时，我们严重限制参数的大小。$\theta_i$的值很小，特征$\phi_i$对模型的贡献很小。模型参数的允许区域收缩，模型变得简单得多。

• 当$Q$很大时，我们并不严重限制参数的大小。$\theta_i$的值很大，特征$\phi_i$对模型的贡献更大。模型参数的允许区域扩大，模型变得更复杂。

考虑当$Q$极大时的极端情况。在这种情况下，我们的限制基本上没有效果，允许的区域包括 OLS 解！

现在如果$Q$非常小会怎么样？我们的参数将被设置为（基本上是 0）。如果模型没有截距项：$\hat{\mathbb{Y}} = (0)\phi_1 + (0)\phi_2 + \ldots = 0$。如果模型有一个截距项：$\hat{\mathbb{Y}} = (0)\phi_1 + (0)\phi_2 + \ldots = \theta_0$。请记住，截距项被排除在约束之外——这样我们就避免了总是预测 0 的情况。

让我们总结一下我们所看到的。

16.4 L1（LASSO）正则化

我们如何实际应用我们的约束$\sum_{i=1}^p |\theta_i| \leq Q$？我们将通过修改我们在拟合模型时寻求最小化的目标函数来实现。

回想一下我们的普通最小二乘目标函数：我们的目标是找到最小化模型均方误差的参数。

$\frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 \phi_{i, 1} + \theta_2 \phi_{i, 2} + \ldots + \theta_p \phi_{i, p}))^2$

要应用我们的约束，我们需要重新表述我们的最小化目标。

$\frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 \phi_{i, 1} + \theta_2 \phi_{i, 2} + \ldots + \theta_p \phi_{i, p}))^2\:\text{such that} \sum_{i=1}^p |\theta_i| \leq Q$

不幸的是，我们不能直接使用这个公式作为我们的目标函数——在约束上进行数学优化并不容易。相反，我们将应用拉格朗日对偶的魔力。这方面的细节超出了范围（如果你有兴趣了解更多，请参加 EECS 127 课程），但最终结果非常有用。事实证明，最小化以下增广目标函数等同于我们上面的最小化目标。

$\frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 \phi_{i, 1} + \theta_2 \phi_{i, 2} + \ldots + \theta_p \phi_{i, p}))^2 + \lambda \sum_{i=1}^p \vert \theta_i \vert = ||\mathbb{Y} - \mathbb{X}\theta||_2^2 + \lambda \sum_{i=1}^p |\theta_i|$

这两个表达式中的第二个包括使用向量表示的 MSE。

注意，我们已经用目标函数中的第二项替换了约束。我们现在正在最小化一个带有额外正则化项的函数，该项惩罚大的系数。为了最小化这个新的目标函数，我们最终会平衡两个组成部分：

• 保持模型在训练数据上的误差低，表示为术语$\frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 x_{i, 1} + \theta_2 x_{i, 2} + \ldots + \theta_p x_{i, p}))^2$

• 同时，保持模型参数的幅度低，表示为术语$\lambda \sum_{i=1}^p |\theta_i|$

$\lambda$因子控制正则化的程度。粗略地说，$\lambda$与之前的$Q$约束相关，规则为$\lambda \approx \frac{1}{Q}$。

为了理解原因，让我们考虑两个极端的例子：

• 假设$\lambda \rightarrow \infty$。那么，$\lambda \sum_{j=1}^{d} \vert \theta_j \vert$ 主导成本函数。为了最小化这个项，我们对所有$j \ge 1$设置$\theta_j = 0$。这是一个非常受限的模型，从数学上讲等同于常数模型。早些时候，我们解释了当 L2 范数球半径$Q \rightarrow 0$时，常数模型也会出现。

• 假设$\lambda \rightarrow 0$。那么，$\lambda \sum_{j=1}^{d} \vert \theta_j \vert$为 0。最小化成本函数等价于$\min_{\theta} \frac{1}{n} || Y - X\theta ||_2^2$，我们通常的 MSE 损失函数。最小化 MSE 损失的行为就是我们熟悉的 OLS，最优解是全局最小值$\hat{\theta} = \hat\theta_{No Reg.}$。我们表明当 L2 范数球半径$Q \rightarrow \infty$时，全局最优解被实现。

我们称$\lambda$为正则化惩罚超参数，并通过交叉验证选择其值。

找到最优$\hat{\theta}$以最小化我们的新目标函数的过程称为L1 正则化。它有时也被称为首字母缩写“LASSO”，代表“最小绝对收缩和选择算子”。

与普通最小二乘法不同，可以通过封闭形式解$\hat{\theta}_{OLS} = (\mathbb{X}^{\top}\mathbb{X})^{-1}\mathbb{X}^{\top}\mathbb{Y}$来解决，L1 正则化下的最优参数向量没有封闭形式解。相反，我们使用sklearn的Lasso模型类。

import sklearn.linear_model as lm

# The alpha parameter represents our lambda term
lasso_model = lm.Lasso(alpha=2)
lasso_model.fit(X_train, Y_train)

lasso_model.coef_

array([-2.54932056e-01, -9.48597165e-04,  8.91976284e-06, -1.22872290e-08])

注意所有模型系数的幅度都非常小。实际上，其中一些系数非常小，基本上为 0。L1 正则化的一个重要特征是许多模型参数被设置为 0。换句话说，LASSO 有效地只选择了一部分特征。这一原因可以追溯到我们先前的损失曲面和允许的“菱形”区域 - 我们通常可以在菱形的一个角附近更接近最低损失轮廓，而不是沿着边缘。

当模型参数设置为 0 或接近 0 时，其对应的特征基本上从模型中移除了。我们说 L1 正则化执行特征选择，因为通过将不重要特征的参数设置为 0，LASSO“选择”了哪些特征对建模更有用。

16.5 特征缩放用于正则化

我们刚刚执行的正则化过程有一个微妙的问题。为了看清楚，让我们来看看我们的lasso_model的设计矩阵。

X_train.head()

	hp	hp^2	hp^3	hp^4
259	85.0	7225.0	614125.0	52200625.0
129	67.0	4489.0	300763.0	20151121.0
207	102.0	10404.0	1061208.0	108243216.0
302	70.0	4900.0	343000.0	24010000.0
71	97.0	9409.0	912673.0	88529281.0

我们的特征——hp、hp^2、hp^3和hp^4——在数值尺度上有着截然不同的差异！hp^4中的值比hp中的值大几个数量级！这可能是一个问题，因为hp^4的值自然上会对每个预测的$\hat{y}$贡献更多，因为它比其他特征的值大得多。对于hp对每个预测产生影响，它必须被一个大的模型参数所缩放。

通过检查我们模型的拟合参数，我们发现这种情况确实存在——hp的参数的数量级远大于hp^4的参数。

pd.DataFrame({"Feature":X_train.columns, "Parameter":lasso_model.coef_})

	Feature	Parameter
0	hp	-2.549321e-01
1	hp^2	-9.485972e-04
2	hp^3	8.919763e-06
3	hp^4	-1.228723e-08

通过应用正则化，我们给我们的模型一个“预算”，来分配模型参数的值。为了让hp对每个预测产生影响，LASSO 被迫在hp的参数上“花费”更多的预算。

我们可以通过在正则化之前对数据进行缩放来避免这个问题。这是一个过程，我们将所有特征转换为相同的数值尺度。一个常见的缩放数据的方法是进行标准化，使得所有特征的均值为 0，标准差为 1；基本上，我们用 Z 分数替换所有内容。

$z_k = \frac{x_k - \mu_k}{\sigma_k}$

16.6 L2（岭）正则化

在我们上面的所有工作中，我们考虑了约束$\sum_{i=1}^p |\theta_i| \leq Q$来限制模型的复杂性。如果我们应用了不同的约束会怎样呢？

在 L2 正则化中，也被称为岭回归，我们约束模型，使得平方参数的总和必须小于某个数$Q$。这个约束的形式如下：

$\sum_{i=1}^p \theta_i^2 \leq Q$

与以前一样，我们通常不对截距项进行正则化。

对于给定的$Q$值，参数的允许区域现在呈球状。

如果我们像之前一样修改我们的目标函数，我们发现我们的新目标是最小化函数：$\frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 \phi_{i, 1} + \theta_2 \phi_{i, 2} + \ldots + \theta_p \phi_{i, p}))^2\:\text{such that} \sum_{i=1}^p \theta_i^2 \leq Q$

请注意，我们所做的只是改变了模型参数的约束。表达式中的第一项，均方误差，没有改变。

使用拉格朗日对偶性，我们可以重新表达我们的目标函数为：$\frac{1}{n} \sum_{i=1}^n (y_i - (\theta_0 + \theta_1 \phi_{i, 1} + \theta_2 \phi_{i, 2} + \ldots + \theta_p \phi_{i, p}))^2 + \lambda \sum_{i=1}^p \theta_i^2 = ||\mathbb{Y} - \mathbb{X}\theta||_2^2 + \lambda \sum_{i=1}^p \theta_i^2$

应用 L2 正则化时，我们的目标是最小化这个更新的目标函数。

与 L1 正则化不同，L2 正则化在应用正则化时确实有一个最佳参数向量的封闭形式解：

$\hat\theta_{\text{ridge}} = (\mathbb{X}^{\top}\mathbb{X} + n\lambda I)^{-1}\mathbb{X}^{\top}\mathbb{Y}$

即使$\mathbb{X}$不是完全列秩，这个解仍然存在。这是 L2 正则化经常被使用的一个主要原因——即使特征中存在共线性，它也可以产生一个解。我们将在未来的讲座中讨论共线性的概念。我们不会在 Data 100 中推导这个结果，因为它涉及相当多的矩阵微积分。

在sklearn中，我们使用Ridge类来执行 L2 正则化。请注意，在正则化之前我们会对数据进行缩放。

ridge_model = lm.Ridge(alpha=1) # alpha represents the hyperparameter lambda
ridge_model.fit(X_train, Y_train)

ridge_model.coef_

array([ 5.89130559e-02, -6.42445915e-03,  4.44468157e-05, -8.83981945e-08])

16.7 回归总结

我们的回归模型总结如下。请注意，目标函数是梯度下降优化器最小化的内容。

类型	模型	损失	正则化	目标函数	解决方案
普通最小二乘法（OLS）	$\hat{\mathbb{Y}} = \mathbb{X}\theta$	均方误差	无	$\frac{1}{n} \|\mathbb{Y}-\mathbb{X} \theta\|^2_2$	$\hat{\theta}_{OLS} = (\mathbb{X}^{\top}\mathbb{X})^{-1}\mathbb{X}^{\top}\mathbb{Y}$ 如果 $\mathbb{X}$ 是满秩的
岭回归	$\hat{\mathbb{Y}} = \mathbb{X} \theta$	均方误差	L2	$\frac{1}{n} \	\mathbb{Y}-\mathbb{X}\theta\	^2_2 + \lambda \sum_{i=1}^p \theta_i^2$	$\hat{\theta}_{ridge} = (\mathbb{X}^{\top}\mathbb{X} + n \lambda I)^{-1}\mathbb{X}^{\top}\mathbb{Y}$
LASSO	$\hat{\mathbb{Y}} = \mathbb{X} \theta$	均方误差	L1	$\frac{1}{n} \	\mathbb{Y}-\mathbb{X}\theta\	^2_2 + \lambda \sum_{i=1}^p \vert \theta_i \vert$	无闭式解

十七、随机变量

原文：Random Variables

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 以其分布的形式定义随机变量

• 计算随机变量的期望和方差

• 熟悉伯努利和二项式随机变量

在过去的几节课中，我们已经考虑了复杂性对模型性能的影响。我们已经考虑了模型复杂性在两个竞争因素之间的权衡：模型方差和训练误差。

到目前为止，我们的分析大部分是定性的。我们已经承认我们对模型复杂性的选择需要在模型方差和训练误差之间取得平衡，但我们还没有讨论为什么会存在这种权衡。

为了更好地理解这种权衡的起源，我们需要引入随机变量的语言。接下来的两节关于概率的讲座将是对我们在建模工作中的一个简短的离题，这样我们就可以建立起理解这所谓的偏差-方差权衡所需的概念。我们接下来几节的路线图将是：

1. 随机变量估计器：引入随机变量，考虑期望、方差和协方差的概念

2. 估计器、偏差和方差：用随机变量的术语重新表达模型方差和训练误差的概念，并利用这种新的视角来研究我们对模型复杂性的选择

Data 8 复习

1. 回顾 Data 8 中的以下概念：

2. 样本均值：你的随机样本的均值

3. 中心极限定理：如果你进行了一个带有替换的大样本随机抽样，那么无论总体分布如何，样本均值的概率分布

   1. 大致是正态的

   2. 以总体均值为中心

   3. 有一个$SD = \frac{\text{总体 SD}}{\sqrt{\text{样本大小}}}$

17.1 随机变量和分布

假设我们生成了一组随机数据，比如从某个总体中随机抽取的一个随机样本。随机变量是数据中随机性的数值函数。它是随机的，因为我们的样本是随机抽取的；它是变量的，因为它的确切值取决于这个随机样本的结果。因此，我们的随机变量的定义域或输入是一个样本空间中所有可能的（随机的）结果，它的值域或输出是数轴。我们通常用大写字母表示随机变量，如$X$或$Y$。

17.1.1 分布

对于任何随机变量$X$，我们需要能够指定两件事：

1. 可能的值：随机变量可以取得的值的集合。

2. 概率：描述 100%总概率如何分布在可能值上的概率集合。

如果$X$是离散的（有有限个可能的值），随机变量$X$取值$x$的概率由$P(X=x)$给出，概率必须总和为 1：$\sum_{\text{all} x} P(X=x) = 1$,

我们通常可以使用概率分布表来显示这一点，你将在下面的抛硬币示例中看到。

随机变量$X$的分布是对 100%总概率如何分布在$X$的所有可能值上的描述，它完全定义了一个随机变量。离散随机变量的分布也可以用直方图表示。如果一个变量是连续的 - 它可以取无限多个值 - 我们可以用密度曲线来说明它的分布。

概率是区域。对于离散随机变量，红色条形的面积表示离散随机变量$X$落在这些值范围内的概率。对于连续随机变量，曲线下的面积表示离散随机变量$Y$落在这些值范围内的概率。

如果我们将条形图/密度曲线下的总面积相加，应该得到 100%，或 1。

17.1.2 例子：抛硬币

举个具体的例子，让我们正式定义一个公平的硬币抛掷。一枚公平的硬币可以正面朝上($H$)或反面朝上($T$)，每种情况的概率都是 0.5。有了这些可能的结果，我们可以将随机变量$X$定义为：$X = \begin{cases} 1, \text{如果硬币正面朝上} \\ 0, \text{如果硬币反面朝上} \end{cases}$

$X$是一个具有域或输入$\{H, T\}$和值域或输出$\{1, 0\}$的函数。我们可以用函数符号表示为$\begin{cases} X(H) = 1 \\ X(T) = 0 \end{cases}$ $X$的概率分布表如下。

$x$	$P(X=x)$
0	$\frac{1}{2}$
1	$\frac{1}{2}$

假设我们从 Data 100 中所有注册学生中随机抽取一个大小为 3 的样本$s$。我们可以将$Y$定义为我们样本中数据科学学生的数量。它的域是大小为 3 的所有可能样本，其值域是$\{0, 1, 2, 3\}$。

我们可以在下表中显示$Y$的分布。左侧的表列出了所有可能的样本$s$及其出现次数($Y(s)$)。我们可以使用这个来计算右侧的表的值，即概率分布表。

17.1.3 模拟

给定随机变量$X$的分布，我们如何生成/模拟一个总体？为此，我们可以根据其分布随机选择$X$的值，使用np.random.choice或df.sample。

17.2 期望和方差

描述随机变量的方法有几种。上面显示的方法 - 所有样本$s, X(s)$的表，分布表$P(X=x)$和直方图 - 都是完全描述随机变量的定义。通常，用一些数值摘要来描述随机变量比完全定义其分布更容易。这些数值摘要是表征随机变量某些属性的数字。因为它们给出了随机变量的行为倾向的“摘要”，它们不是随机的 - 将它们视为描述随机变量某个属性的静态数字。在 Data 100 中，我们将关注随机变量的期望和方差。

17.2.1 期望

随机变量$X$的期望是$X$的值的加权平均值，其中权重是每个值发生的概率。有两种等效的计算期望的方法：

1. 一次应用一个样本的权重：$\mathbb{E}[X] = \sum_{\text{所有可能的} s} X(s) P(s)$。

2. 一次应用权重一个可能的值：$\mathbb{E}[X] = \sum_{\text{所有可能的} x} x P(X=x)$

我们要强调的是，期望是一个数字，不是一个随机变量。期望是平均值的一种概括，它与随机变量具有相同的单位。它也是概率分布直方图的重心，这意味着如果我们多次模拟变量，它是随机变量的长期平均值。

17.2.1.1 示例 1：抛硬币

回到我们抛硬币的例子，我们将随机变量$X$定义为：$X = \begin{cases} 1, \text{如果硬币正面朝上} \\ 0, \text{如果硬币反面朝上} \end{cases}$ 我们可以使用第二种方法，一次应用权重一个可能的值来计算其期望$\mathbb{E}[X]$：\begin{align} \mathbb{E}[X] &= \sum_{x} x P(X=x) \\ &= 1 * 0.5 + 0 * 0.5 \\ &= 0.5 \end{align} 请注意，$\mathbb{E}[X] = 0.5$不是$X$的可能值；这是一个平均值。X 的期望值不需要是 X 的可能值。

17.2.1.2 示例 2

考虑随机变量$X$：

$x$	$P(X=x)$
3	0.1
4	0.2
6	0.4
8	0.3

计算期望，\begin{align} \mathbb{E}[X] &= \sum_{x} x P(X=x) \\ &= 3 * 0.1 + 4 * 0.2 + 6 * 0.4 + 8 * 0.3 \\ &= 0.3 + 0.8 + 2.4 + 2.4 \\ &= 5.9 \end{align} 再次注意，$\mathbb{E}[X] = 5.9$ 不是 $X$ 的可能值；这是一个平均值。X 的期望值不需要是 X 的可能值。

17.2.2 方差

随机变量的方差是其随机误差的度量。它被定义为$X$的期望值的平方偏差。更简单地说，方差问：$X$通常从其平均值变化多少，仅仅是由于偶然？$X$的分布是如何传播的？

$\text{Var}(X) = \mathbb{E}[(X-\mathbb{E}[X])^2]$

方差的单位是$X$的单位的平方。要将其恢复到正确的比例，使用$X$的标准差：$\text{SD}(X) = \sqrt{\text{Var}(X)}$

与期望一样，方差是一个数，不是一个随机变量！它的主要用途是量化偶然误差。

根据切比雪夫不等式，你在 Data 8 中看到的，无论$X$的分布形状如何，绝大多数的概率都在“期望值加上或减去几个标准差”的区间内。

如果我们展开平方并使用期望的属性，我们可以重新表达方差作为方差的计算公式。当手动计算变量的方差时，这种形式通常更方便使用，并且在均方误差计算中也很有用，因为如果$X$是居中的，那么$\mathbb{E}[X^2] = \text{Var}(X)$。

$\text{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$

证明

\begin{align} \text{Var}(X) &= \mathbb{E}[(X-\mathbb{E}[X])^2] \\ &= \mathbb{E}(X^2 - 2X\mathbb{E}(X) + (\mathbb{E}(X))^2) \\ &= \mathbb{E}(X^2) - 2 \mathbb{E}(X)\mathbb{E}(X) +( \mathbb{E}(X))^2\\ &= \mathbb{E}[X^2] - (\mathbb{E}[X])^2 \end{align}

我们如何计算$\mathbb{E}[X^2]$? 任何随机变量的函数也是随机变量 - 这意味着通过平方$X$，我们创建了一个新的随机变量。要计算$\mathbb{E}[X^2]$，我们可以简单地将我们的期望值定义应用于随机变量$X^2$。

$\mathbb{E}[X^2] = \sum_{x} x^2 P(X = x)$

17.2.3 例子：骰子

设$X$是单次公平掷骰子的结果。$X$是一个随机变量，定义为$X = \begin{cases} \frac{1}{6}, \text{if } x \in \{1,2,3,4,5,6\} \\ 0, \text{otherwise} \end{cases}$

期望值 $\mathbb{E}[X]?$

\begin{align} \mathbb{E}[X] &= 1(\frac{1}{6}) + 2(\frac{1}{6}) + 3(\frac{1}{6}) + 4(\frac{1}{6}) + 5(\frac{1}{6}) + 6(\frac{1}{6}) \\ &= (\frac{1}{6}) ( 1 + 2 + 3 + 4 + 5 + 6) \\ &= \frac{7}{2} \end{align}

方差 $\text{Var}(X)?$

使用方法 1: \begin{align} \text{Var}(X) &= (\frac{1}{6})((1 - \frac{7}{2})^2 + (2 - \frac{7}{2})^2 + (3 - \frac{7}{2})^2 + (4 - \frac{7}{2})^2 + (5 - \frac{7}{2})^2 + (6 - \frac{7}{2})^2) \\ &= \frac{35}{12} \end{align}

使用方法 2：$\mathbb{E}[X^2] = \sum_{x} x^2 P(X = x) = \frac{91}{6}$

$\text{Var}(X) = \frac{91}{6} - (\frac{7}{2})^2 = \frac{35}{12}$

17.3 随机变量的和

通常，我们会同时处理多个随机变量。随机变量的函数也是随机变量；如果你基于样本创建多个随机变量，那么这些随机变量的函数也是随机变量。

例如，如果$X_1, X_2, ..., X_n$是随机变量，那么这些也是随机变量：

• $X_n^2$

• $\#\{i : X_i > 10\}$

• $\text{max}(X_1, X_2, ..., X_n)$

• $\frac{1}{n} \sum_{i=1}^n (X_i - c)^2$

• $\frac{1}{n} \sum_{i=1}^n X_i$

17.3.1 相等 vs. 相同分布 vs. i.i.d.

假设我们有两个随机变量 $X$ 和 $Y$：

• $X$ 和 $Y$ 如果对于每个样本 $s$ 都有 $X(s) = Y(s)$ 则它们是相等的。无论抽取的确切样本是什么，$X$ 总是等于 $Y$。

• 如果$X$和$Y$的分布相等，则$X$和$Y$是相同分布的。我们说“X 和 Y 在分布上相等”。也就是说，$X$和$Y$取相同的可能值集，并且每个可能值都以相同的概率取到。在任何特定的样本$s$上，相同分布的变量不一定共享相同的值。如果 X = Y，则 X 和 Y 是相同分布的；然而，反之则不成立（例如：Y = 7-X，X 是一个骰子）

• 如果

  1. 这些变量是相同分布的。

  2. 知道一个变量的结果不会影响我们对另一个变量结果的信念。

例如，让$X_1$和$X_2$是两个公平骰子的点数。$X_1$和$X_2$是 i.i.d，所以$X_1$和$X_2$具有相同的分布。然而，和$Y = X_1 + X_1 = 2X_1$和$Z=X_1+X_2$具有不同的分布，但是相同的期望值。

然而，$Y = X_1$的方差更大

17.3.2 期望值的性质

我们经常直接计算期望值和方差，而不是模拟完整的分布。回顾期望值的定义：$\mathbb{E}[X] = \sum_{x} x P(X=x)$ 从中，我们可以推导出期望值的一些有用的性质：

1. 期望的线性性。常数$a$和$b$的线性变换$aX+b$的期望值是：

$\mathbb{E}[aX+b] = aE[\mathbb{X}] + b$

证明

\begin{align} \mathbb{E}[aX+b] &= \sum_{x} (ax + b) P(X=x) \\ &= \sum_{x} (ax P(X=x) + bP(X=x)) \\ &= a\sum_{x}P(X=x) + b\sum_{x}P(X=x)\\ &= a\mathbb{E}(X) + b * 1 \end{align}

2. 期望值在随机变量的总和中也是线性的。

$\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]$

证明

\begin{align} \mathbb{E}[X+Y] &= \sum_{s} (X+Y)(s) P(s) \\ &= \sum_{s} (X(s)P(s) + Y(s)P(s)) \\ &= \sum_{s} X(s)P(s) + \sum_{s} Y(s)P(s)\\ &= \mathbb{E}[X] + \mathbb{E}[Y] \end{align}

3. 如果$g$是一个非线性函数，那么一般来说，$\mathbb{E}[g(X)] \neq g(\mathbb{E}[X])$ 例如，如果$X$以相等的概率为-1 或 1，那么$\mathbb{E}[X] = 0$，但$\mathbb{E}[X^2] = 1 \neq 0$.

17.3.3 方差的性质

回顾方差的定义：$\text{Var}(X) = \mathbb{E}[(X-\mathbb{E}[X])^2]$ 结合期望值的性质，我们可以推导出方差的一些有用的性质：

1. 与期望值不同，方差是非线性的。线性变换$aX+b$的方差是：$\text{Var}(aX+b) = a^2 \text{Var}(X)$

• 随后，$\text{SD}(aX+b) = |a| \text{SD}(X)$

• 可以通过方差的定义找到这个事实的完整证明。作为一般直觉，考虑$aX+b$通过因子$a$缩放变量$X$，然后将$X$的分布移位$b$单位。

证明

我们知道$\mathbb{E}[aX+b] = aE[\mathbb{X}] + b$

为了计算$\text{Var}(aX+b)$，考虑到 b 单位的移位不会影响扩展，因此$\text{Var}(aX+b) = \text{Var}(aX)$。

因此，\begin{align} \text{Var}(aX+b) &= \text{Var}(aX) \\ &= E((aX)^2) - (E(aX))^2 \\ &= E(a^2 X^2) - (aE(X))^2\\ &= a^2 (E(X^2) - (E(X))^2) \\ &= a^2 \text{Var}(X) \end{align}

• 将分布移位 b 不会 影响分布的扩展。因此，$\text{Var}(aX+b) = \text{Var}(aX)$。

• 通过$a$缩放分布会 影响分布的扩展。

1. 随机变量的和的方差受到随机变量的（不）独立性的影响。$\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y) + 2\text{cov}(X,Y)$

$\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y) \qquad \text{if } X, Y \text{ independent}$

证明

两个随机变量相加的方差受到它们之间的依赖关系的影响。让我们展开$\text{Var}(X + Y)$的定义，看看发生了什么。

为了简化数学，让$\mu_x = \mathbb{E}[X]$和$\mu_y = \mathbb{E}[Y]$。

\begin{align} \text{Var}(X + Y) &= \mathbb{E}[(X+Y- \mathbb{E}(X+Y))^2] \\ &= \mathbb{E}[((X - \mu_x) + (Y - \mu_y))^2] \\ &= \mathbb{E}[(X - \mu_x)^2 + 2(X - \mu_x)(Y - \mu_y) + (Y - \mu_y)^2] \\ &= \mathbb{E}[(X - \mu_x)^2] + \mathbb{E}[(Y - \mu_y)^2] + \mathbb{E}[(X - \mu_x)(Y - \mu_y)] \\ &= \text{Var}(X) + \text{Var}(Y) + \mathbb{E}[(X - \mu_x)(Y - \mu_y)] \end{align}

17.3.4 Covariance and Correlation

我们将两个随机变量的协方差定义为期望的偏差乘积。更简单地说，协方差是方差对两个随机变量的泛化：$\text{Cov}(X, X) = \mathbb{E}[(X - \mathbb{E}[X])^2] = \text{Var}(X)$

$\text{Cov}(X, Y) = \mathbb{E}[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])]$

我们可以将协方差视为一种关联度量。还记得我们在建立简单线性回归时给出的相关性定义吗？

$r(X, Y) = \mathbb{E}\left[\left(\frac{X-\mathbb{E}[X]}{\text{SD}(X)}\right)\left(\frac{Y-\mathbb{E}[Y]}{\text{SD}(Y)}\right)\right] = \frac{\text{Cov}(X, Y)}{\text{SD}(X)\text{SD}(Y)}$

事实证明我们一直在悄悄使用协方差！如果$X$和$Y$是独立的，那么$\text{Cov}(X, Y) =0$和$r(X, Y) = 0$。然而，请注意，逆命题并不总是成立：$X$和$Y$可能有$\text{Cov}(X, Y) = r(X, Y) = 0$但并不是独立的。

17.3.5 Summary

• 设$X$是一个具有分布$P(X=x)$的随机变量。

  • $\mathbb{E}[X] = \sum_{x} x P(X=x)$

  • $\text{Var}(X) = \mathbb{E}[(X-\mathbb{E}[X])^2] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$

• 设$a$和$b$是标量值。

  • $\mathbb{E}[aX+b] = aE[\mathbb{X}] + b$

  • $\text{Var}(aX+b) = a^2 \text{Var}(X)$

• 设$Y$是另一个随机变量。

  • $\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]$

  • $\text{Var}(X + Y) = \text{Var}(X) + \text{Var}(Y) 2\text{cov}(X,Y)$

十八、估计器、偏差和方差

Estimators, Bias, and Variance

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 探索常见的随机变量，如伯努利和二项式分布

• 应用中心极限定理来近似总体参数

• 使用抽样数据对真实的潜在分布进行建模估计

• 使用自助法技术从样本中估计真实总体分布

上次，我们介绍了随机变量的概念：样本的数值函数。在上一讲中，我们的大部分工作是建立概率和统计学的背景。现在我们已经建立了一些关键的思想，我们可以将我们学到的知识应用到我们最初的目标上 - 理解样本的随机性如何影响模型设计过程。

在本讲座中，我们将更深入地探讨将模型拟合到样本的想法。我们将探讨如何用随机变量重新表达我们的建模过程，并利用这种新的理解来引导模型的复杂性。

18.1 常见随机变量

有几种经常出现并且具有有用特性的随机变量情况。以下是我们将在本课程中进一步探讨的情况。括号中的数字是随机变量的参数，这些参数是常数。参数定义了随机变量的形状（即分布）和其值。在本讲座中，我们将更加重点关注加粗的随机变量及其特殊性质，但你应该熟悉下面列出的所有随机变量：

• 伯努利(p)

  • 以概率 p 取值 1，以概率 1-p 取值 0。

  • 又称“指示”随机变量。

  • 设 X 是一个伯努利(p)随机变量

    • $\mathbb{E}[X] = 1 * p + 0 * (1-p) = p$

      • $\mathbb{E}[X^2] = 1^2 * p + 0 * (1-p) = p$

    • $\text{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2 = p - p^2 = p(1-p)$

• 二项式(n, p)

  • $n$ 独立伯努利(p)试验中的 1 的数量。

  • 设$Y$是一个二项式(n, p)随机变量。

    • $Y$的分布由二项式公式给出，我们可以写成$Y = \sum_{i=1}^n X_i$，其中：

      • $X_i$是第$i$次试验成功的指示。如果第$i$次试验成功，则$X_i = 1$，否则为 0。

      • 所有的$X_i$都是独立同分布的伯努利(p)。

    • $\mathbb{E}[Y] = \sum_{i=1}^n \mathbb{E}[X_i] = np$

    • $\text{Var}(X) = \sum_{i=1}^n \text{Var}(X_i) = np(1-p)$

      • $X_i$是独立的，所以对于所有的 i, j，$\text{Cov}(X_i, X_j) = 0$。

• 有限值集上均匀分布

  • 每个值的概率是 1 / (可能的值的数量)。

  • 例如，一个标准/公平的骰子。

• 单位区间(0, 1)上均匀分布

  • 密度在(0, 1)上为 1，在其他地方为 0。

• 正态($\mu, \sigma^2$)

  • $f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left( -\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{\!2}\,\right)$

18.1.1 例子

假设你根据 20 次抛硬币中得到的正面数量赢得现金。如果第$i$次抛硬币得到正面，则令$X_i = 1$，否则为 0。你会选择哪种支付策略？

A. $Y_A = 10 * X_1 + 10 * X_2$

B. $Y_B = \sum_{i=1}^{20} X_i$

C. $Y_C = 20 * X_1$

解决方案

设$X_1, X_2, ... X_{20}$是 20 个独立同分布的伯努利(0.5)随机变量。由于$X_i$是独立的，对于所有的$i, j$对，$\text{Cov}(X_i, X_j) = 0$。另外，由于$X_i$是伯努利(0.5)，我们知道$\mathbb{E}[X] = p = 0.5$和$\text{Var}(X) = p(1-p) = 0.25$。我们可以计算每种情况的如下内容：

	A. $Y_A = 10 * X_1 + 10 * X_2$	B. $Y_B = \sum_{i=1}^{20} X_i$	C. $Y_C = 20 * X_1$
期望	$\mathbb{E}[Y_A] = 10 (0.5) + 10(0.5) = 10$	$\mathbb{E}[Y_B] = 0.5 + ... + 0.5 = 10$	$\mathbb{E}[Y_C] = 20(0.5) = 10$
方差	$\text{Var}(Y_A) = 10^2 (0.25) + 10^2 (0.25) = 50$	$\text{Var}(Y_B) = 0.25 + ... + 0.25 = 5$	$\text{Var}(Y_C) = 20^2 (0.25) = 100$
标准差	$\text{SD}(Y_A) \approx 7.07$	$\text{SD}(Y_B) \approx 2.24$	$\text{SD}(Y_C) = 10$

正如我们所看到的，所有的情景都有相同的期望值，但方差不同。方差越大，风险和不确定性就越大，因此“正确”的策略取决于个人的偏好。你会选择“最安全”的选项 B，最“冒险”的选项 C，还是介于两者之间的选项 A？

18.2 样本统计

今天，我们已经广泛讨论了总体；如果我们知道随机变量的分布，我们可以可靠地计算期望、方差、随机变量的函数等。请注意：

• 总体的分布描述了随机变量在所有感兴趣的个体中的行为。

• 样本的分布描述了随机变量在来自总体的特定样本中的行为。

然而，在数据科学中，我们经常无法接触到整个总体，因此我们不知道它的分布。因此，我们需要收集一个样本，并使用它的分布来估计或推断总体的属性。在这种情况下，我们可以从总体中取几个大小为$n$的样本（一个简单的方法是使用df.sample(n, replace=True)），并计算每个样本的均值。在抽样时，我们做出（很大的）假设，即我们从总体中均匀随机地进行有放回抽样；我们样本中的每个观察都是从我们的总体分布中独立同分布地随机抽取的随机变量。请记住，我们的样本均值是一个随机变量，因为它取决于我们随机抽取的样本！另一方面，我们的总体均值只是一个数字（一个固定的值）。

18.2.1 样本均值

考虑一个从具有均值𝜇和标准差𝜎的总体中抽取的 i.i.d.样本$X_1, X_2, ..., X_n$。我们定义样本均值为$\bar{X}_n = \frac{1}{n} \sum_{i=1}^n X_i$

样本均值的期望值由以下公式给出：\begin{align} \mathbb{E}[\bar{X}_n] &= \frac{1}{n} \sum_{i=1}^n \mathbb{E}[X_i] \\ &= \frac{1}{n} (n \mu) \\ &= \mu \end{align}

方差由以下公式给出：\begin{align} \text{Var}(\bar{X}_n) &= \frac{1}{n^2} \text{Var}( \sum_{i=1}^n X_i) \\ &= \frac{1}{n^2} \left( \sum_{i=1}^n \text{Var}(X_i) \right) \\ &= \frac{1}{n^2} (n \sigma^2) = \frac{\sigma^2}{n} \end{align}

$\bar{X}_n$根据中心极限定理（CLT）呈正态分布。

18.2.2 中心极限定理

在Data 8和之前的讲座中，你遇到了中心极限定理（CLT）。这是一个强大的定理，用于从一系列较小的样本中估计具有均值$\mu$和标准差$\sigma$的总体的分布。中心极限定理告诉我们，如果一个大小为$n$的 i.i.d 样本很大，那么样本均值的概率分布大致正态，均值为$\mu$，标准差为$\frac{\sigma}{\sqrt{n}}$。更一般地，任何提供统计量粗略分布并且不需要总体分布的定理对于数据科学家来说都是有价值的！这是因为我们很少对总体了解很多。

重要的是，中心极限定理假设我们样本中的每个观察都是从总体的分布中抽取的 i.i.d。此外，中心极限定理仅在$n$“大”时才准确，但什么样的“大”样本量取决于特定的分布。如果一个总体高度对称和单峰，我们可能只需要$n=20$；如果一个总体非常倾斜，我们需要更大的$n$。如果有疑问，可以对样本均值进行自举，并查看自举分布是否呈钟形。像 Data 140 这样的课程会对这个想法进行详细的探讨。

要了解更详细的演示，请查看onlinestatbook。

18.2.3 使用样本均值估计总体均值

现在假设我们想使用样本均值来估计总体均值，例如，加州大学本科生的平均身高。通常我们可以收集一个单一样本，其中只有一个平均值。但是，如果我们碰巧以随机方式抽取了一个具有不同均值或扩展性的样本，会怎么样呢？我们可能会对总体行为有一个偏斜的看法（考虑极端情况，我们碰巧抽取了相同的值 $n$ 次！）。

例如，注意这两个分布之间的变化差异，这两个分布在样本大小上是不同的。样本量更大的分布（$n=800$）比样本量较小的分布（$n=200$）更紧密地围绕均值。尝试将这些值代入正态分布的标准偏差方程中，以理解这一点！

应用中心极限定理使我们能够理解所有这些并解决这个问题。通过抽取许多样本，我们可以考虑样本分布在数据的多个子集中的变化。这使我们能够近似总体的属性，而无需调查每个成员。

鉴于这种潜在的差异，我们还要考虑所有可能的样本均值的平均值和扩展性，以及这对 $n$ 应该有多大的影响。对于每个样本量，样本均值的期望值是总体均值：$\mathbb{E}[\bar{X}_n] = \mu$。我们称样本均值是总体均值的无偏估计量，并将在下一讲中更多地探讨这个想法。

Data 8 Recap: 平方根定律

平方根定律（Data 8）指出，如果将样本量增加一个因子，标准偏差将减少该因子的平方根。 $\text{SD}(\bar{X_n}) = \frac{\sigma}{\sqrt{n}}$ 如果我们有更大的样本量，样本均值更有可能接近总体均值。

18.3 预测和推断

在课程的这一阶段，我们花了大量时间研究模型。几周前我们首次介绍了建模的概念时，是在预测的背景下：使用模型对未知数据进行准确预测。我们构建模型的另一个原因是更好地理解我们周围复杂的现象。推断是使用模型推断特征和响应变量之间真实的基本关系的任务。例如，如果我们正在处理一组房屋数据，预测可能会问：根据房屋的属性，它值多少钱？推断可能会问：当地公园对房屋价值有多大影响？

推断的一个主要目标是仅凭随机样本对完整数据总体进行推断。为此，我们旨在估计参数的值，这是总体的数值函数（例如，总体均值 $\mu$）。我们使用收集的样本来构建统计量，这是随机样本的数值函数（例如，样本均值 $\bar{X}_n$）。将“p”视为“参数”和“总体”，将“s”视为“样本”和“统计量”是有帮助的。

由于样本代表总体的随机子集，我们生成的任何统计量可能会偏离真实的总体参数，并且可能会有所不同。我们说样本统计量是真实总体参数的估计量。在符号上，总体参数通常称为 $\theta$，而其估计量用 $\hat{\theta}$ 表示。

为了回答我们的推断问题，我们旨在构建能够紧密估计总体参数值的估计量。我们通过回答三个问题来评估估计量的“好坏”：

• 我们平均得到参数的正确答案吗？

• 答案有多大的变化？

• 我们的答案与参数有多接近？

18.3.1 建模作为估计

现在我们已经建立了估计量的概念，让我们看看如何将这种学习应用到建模过程中。为此，我们将花一点时间用随机变量的语言来形式化我们的数据收集和模型。

假设我们正在处理一个输入变量$x$和一个响应变量$Y$。我们假设$Y$和$x$通过某种关系$g$相关联；换句话说，$Y = g(x)$。$g$代表一些定义$x$和$Y$之间基础关系的“普遍真理”或“自然法则”。在下面的图像中，$g$由红线表示。

然而，作为数据科学家，我们无法直接“看到”基础关系$g$。我们能做的最好的事情就是收集在现实世界中观察到的数据，以尝试理解这种关系。不幸的是，数据收集过程总会存在一些固有的误差（想象一下在科学实验中进行测量时可能遇到的随机性）。我们说每个观察都伴随着一些随机误差或噪声项$\epsilon$。假定这个误差是一个随机变量，期望为$\mathbb{E}(\epsilon)=0$，方差为$\text{Var}(\epsilon) = \sigma^2$，并且在每个观察中都是独立同分布的。这种随机噪声的存在意味着我们的观察$Y(x)$是随机变量。

我们只能观察到我们的随机数据样本，用蓝色点表示。从这个样本中，我们想要估计真实关系$g$。我们通过构建模型$\hat{Y}(x)$来估计$g$。

$\text{真实关系: } g(x)$

$\text{观察到的关系: }Y = g(x) + \epsilon$

$\text{预测: }\hat{Y}(x)$

18.3.1.1 估计线性关系

如果我们假设真实关系$g$是线性的，我们可以将响应表示为$Y = f_{\theta}(x)$，其中我们的真实关系由$Y = g(x) + \epsilon$

$f_{\theta}(x) = Y = \theta_0 + \sum_{j=1}^p \theta_j x_j + \epsilon$来建模。

哪些表达式是随机的？

在我们上面的两个方程中，真实关系$g(x) = \theta_0 + \sum_{j=1}^p \theta_j x_j$不是随机的，但$\epsilon$是随机的。因此，$Y = f_{\theta}(x)$也是随机的。

*这个真实关系有真实的、不可观测的参数$\theta$，并且它有随机噪声$\epsilon$，所以我们永远无法观察到真实关系。相反，我们能做的下一个最好的事情就是获得一个样本$\Bbb{X}$，$\Bbb{Y}$的$n$个观察关系$(x, Y)$，并用它来训练一个模型并获得$\hat{\theta}$的估计$\hat{Y}(x) = f_{\hat{\theta}}(x) = \hat{\theta_0} + \sum_{j=1}^p \hat{\theta_j} x_j$

哪些表达式是随机的？

在我们上面的估计方程中，我们的样本$\Bbb{X}$，$\Bbb{Y}$是随机的。因此，我们从样本中计算的估计$\hat{\theta}$也是随机的，所以我们的预测$\hat{Y}(x)$也是随机的。

现在看一下我们的原始方程，我们可以看到它们都有不同的随机来源。对于我们观察到的关系，$Y = g(x) + \epsilon$，$\epsilon$代表测量误差并反映未来的随机性。对于估计模型，我们拥有的数据是从总体中收集的随机样本，因此是过去的随机性。

18.4 自助法重采样（复习）

确定估计量的抽样分布的属性，比如方差，我们需要访问总体，以便我们可以考虑所有可能的样本并计算每个样本的估计。

然而，我们无法访问总体；我们只有来自总体的一个随机样本。如果我们只有一个样本，我们如何考虑所有可能的样本呢？

自助法的想法是将我们的随机样本视为“总体”，并从中进行有放回的重新采样。直观地说，随机样本类似于总体，因此随机重新采样也重新对随机样本进行重新采样。

Bootstrap 重采样是一种估计估计量抽样分布的技术。要执行它，我们可以按照下面的伪代码进行：

collect a random sample of size n (called the bootstrap population)

initiate list of estimates

repeat 10,000 times:
    resample with replacement n times from bootstrap population

apply estimator f to resample

store in list

list of estimates is the bootstrapped sampling distribution of f

为什么我们必须进行有放回的重新采样？

给定大小为$n$的原始样本，我们希望得到与原始样本相同大小$n$的重新采样。不进行替换的抽样将给我们洗牌后的原始样本。因此，当我们计算像平均值这样的摘要统计时，我们不进行替换的样本将始终具有与原始样本相同的平均值，从而破坏了自助法的目的。

自助法实际上如何代表我们的总体？估计量的自助法抽样分布并不完全匹配该估计量的抽样分布，但通常是接近的。同样，自助法分布的方差通常接近于估计量的真实方差。下面的示例显示了使用样本大小$n=50$从已知总体进行不同自助法的结果。**

在现实世界中，我们不知道总体分布。自助法分布的中心是应用于我们原始样本的估计量，因此我们无法恢复估计量的真实期望值。我们的自助法分布的质量取决于我们原始样本的质量；如果我们的原始样本不代表总体，自助法几乎没有用处。

需要注意的一点是，自助法通常对某些统计量（如中位数或其他基于分位数的统计量）效果不佳，这些统计量严重依赖于较大样本中的少数观察结果。自助法无法克服小样本作为推断依据的弱点。事实上，对于非常小的样本，最好是做出额外的假设，比如参数族。

十九、偏差、方差和推断

原文：Bias, Variance, and Inference

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 计算参数的估计器的偏差、方差和均方误差

• 介绍拟合模型的模型风险

• 将模型风险分解为偏差和方差项

• 构建假设检验的置信区间

• 了解我们所做的假设及其对回归推断的影响

• 比较回归和因果关系

• 实验设置、混杂变量、平均处理效应和协变量调整**上次，我们介绍了随机变量的概念及其对我们用来拟合模型的观察关系的影响。

在本讲座中，我们将探讨从拟合模型中分解模型风险，通过假设检验进行回归推断，并考虑我们所做的假设以及理论和实践中理解因果关系的环境。

19.1 偏差-方差权衡

回顾上一节中建立的模型和我们从该模型中生成的数据：

$\text{真实关系：} g(x)$

$\text{观察到的关系：}Y = g(x) + \epsilon$

$\text{预测：} \hat{Y}(x)$

通过这个重新制定的建模目标，我们现在可以重新审视两次讲座前的偏差-方差权衡（如下所示）：

在今天的讲座中，我们将通过引入模型风险、观测方差、模型偏差和模型方差这些术语，探讨上面所见图表的更数学化版本。最终，我们将更新偏差-方差权衡图表，如下所示

19.1.1 估计器的性能

假设我们想要使用估计器 $\hat{Y}(x)$ 估计目标 $Y$。对于我们训练的每个估计器，我们可以通过以下问题来确定模型的好坏：

• 我们平均得到正确答案吗？(偏差)

• 答案有多大的变化？(方差)

• 我们离 $Y$ 有多近？(风险/MSE)

理想情况下，我们希望我们的估计器偏差和方差都很低，但我们如何在数学上量化呢？为此，让我们引入一些术语。

19.1.2 模型风险

模型风险 被定义为随机变量 $\hat{Y}$ 的均方预测误差。它是对我们拟合模型时可能得到的 所有 样本的期望，我们可以将其表示为随机变量 $X_1, X_2, \ldots, X_n, Y$。模型风险考虑了模型在理论上可能的任何样本上的表现，而不是我们收集到的具体数据。

$\text{模型风险} = E\left[(Y-\hat{Y(x)})^2\right]$

模型风险所编码的错误的起源是什么？请注意，有两种类型的错误：

• 偶然误差：仅由随机性引起

  • 来源 1 (观测方差)：由于随机噪声 $\epsilon$ 导致新观测 $Y$ 的随机性

  • 来源 2 (模型方差)：在我们用来训练模型的样本中的随机性，因为样本 $X_1, X_2, \ldots, X_n, Y$ 是随机的

• (模型偏差)：由于我们的模型与真实的基本函数 $g$ 不同而产生的非随机误差

回顾我们之前建立的数据生成过程。存在一个真实的基本关系 $g$，观察到的数据（带有随机噪声）$Y$，以及模型 $\hat{Y}$。

为了更好地理解模型风险，我们将放大上图中的一个数据点。

记住$\hat{Y}(x)$是一个随机变量 - 它是在用于训练的特定样本上拟合后对$x$的预测。如果我们使用不同的样本进行训练，可能会对这个值的预测进行不同的预测。为了捕捉这一点，上面的图考虑了对特定随机训练样本进行的预测$\hat{Y}(x)$，以及在所有可能的训练样本上的预期预测$E[\hat{Y}(x)]$。

我们可以使用这个简化的图表来将预测误差分解为更小的组件。首先，从单个预测的误差$Y(x)-\hat{Y}(x)$开始。

我们可以确定这个错误的三个组成部分。

也就是说，错误可以写成：

$Y(x)-\hat{Y}(x) = \epsilon + \left(g(x)-E\left[\hat{Y}(x)\right]\right) + \left(E\left[\hat{Y}(x)\right] - \hat{Y}(x)\right)$

模型风险是上述表达式的平方的期望值，$E\left[(Y(x)-\hat{Y}(x))^2\right]$。如果我们两边平方，然后取期望值，我们将得到模型风险的以下分解：

$E\left[(Y(x)-\hat{Y}(x))^2\right] = E[\epsilon^2] + \left(g(x)-E\left[\hat{Y}(x)\right]\right)^2 + E\left[\left(E\left[\hat{Y}(x)\right] - \hat{Y}(x)\right)^2\right]$

看起来当我们平方右边时，我们缺少一些交叉乘积项，但事实证明所有这些交叉乘积项都为零。这门课程的详细推导超出了范围，但在本笔记的末尾包括了一个证明供您参考。

这个表达式乍一看可能很复杂，但实际上我们在本讲座中已经定义了每个术语！让我们逐个术语地来看。

19.1.2.1 观察方差

上述分解中的第一项是$E[\epsilon^2]$。记住$\epsilon$是观察$Y$时的随机噪声，期望为$\mathbb{E}(\epsilon)=0$，方差为$\text{Var}(\epsilon) = \sigma^2$。我们可以证明$E[\epsilon^2]$是$\epsilon$的方差：\begin{align*} \text{Var}(\epsilon) &= E[\epsilon^2] + \left(E[\epsilon]\right)^2\\ &= E[\epsilon^2] + 0^2\\ &= \sigma^2. \end{align*}

这个术语描述了每个观察中随机误差$\epsilon$（和$Y$）的变量。这被称为观察方差。它存在于我们的观察$Y$的随机性中。这是我们在抽样讲座中谈到的偶然误差的一种形式。

$\text{观察方差} = \text{Var}(\epsilon) = \sigma^2.$

观察方差是由观察数据时的测量误差或行为像噪声一样的缺失信息引起的。要减少这种观察方差，我们可以尝试获得更精确的测量，但这通常超出了数据科学家的控制范围。因此，观察方差$\sigma^2$有时被称为“不可减少的误差”。

19.1.2.2 模型方差

然后我们来看最后一项：$E\left[\left(E\left[\hat{Y}(x)\right] - \hat{Y}(x)\right)^2\right]$。如果你回忆一下上一讲的方差的定义，这正是$\text{Var}(\hat{Y}(x))$。我们称之为模型方差。

它描述了当我们在不同样本上拟合模型时，预测$\hat{Y}(x)$往往变化多少。记住我们收集的样本可能会有很大的不同，因此预测$\hat{Y}(x)$也会有所不同。模型方差描述了由于我们抽样过程的随机性而产生的这种变异性。与观察方差一样，它也是一种偶然误差 - 即使随机性的来源是不同的。

$模型方差 = Var(\hat{Y}(x)) = E\left[\left(\hat{Y}(x) - E\left[\hat{Y}(x)\right]\right)^2\right]$

模型方差较大的主要原因是过拟合：我们过于关注样本中的细节，导致随机样本中的微小差异导致拟合模型中的大差异。为了解决这个问题，我们尝试减少模型复杂性（例如去掉一些特征和限制估计模型系数的大小），并且不要在噪声上拟合我们的模型。

19.1.2.3 模型偏差

最后，第二项是 $\left(g(x)-E\left[\hat{Y}(x)\right]\right)^2$。这是什么？术语 $E\left[\hat{Y}(x)\right] - g(x)$ 被称为模型偏差。

记住 $g(x)$ 是固定的基本真相，$\hat{Y}(x)$ 是我们拟合的模型，是随机的。因此，模型偏差衡量了 $g(x)$ 和 $\hat{Y}(x)$ 在所有可能样本上的平均偏差。

$\text{模型偏差} = E\left[\hat{Y}(x) - g(x)\right] = E\left[\hat{Y}(x)\right] - g(x)$

模型偏差不是随机的；它是特定个体 $x$ 的平均度量。如果偏差是正的，我们的模型倾向于高估 $g(x)$；如果是负的，我们的模型倾向于低估 $g(x)$。如果是 0，我们可以说我们的模型是无偏的。

无偏估计

一个无偏模型具有 $\text{模型偏差 } = 0$。换句话说，我们的模型平均预测 $g(x)$。

类似地，我们可以为估计量定义偏差，比如均值。样本均值是总体均值的无偏估计，因为根据中心极限定理，$\mathbb{E}[\bar{X}_n] = \mu$。因此，$\text{估计器偏差 } = \mathbb{E}[\bar{X}_n] - \mu = 0$.

模型偏差较大的两个主要原因是：

• 拟合不足：我们的模型对数据来说太简单。

• 缺乏领域知识：我们不了解哪些特征对响应变量有用

为了解决这个问题，我们增加模型复杂性（但我们不想过拟合！）或请教领域专家，看看哪些模型是合理的。你可以开始看到这里的权衡：如果我们增加模型复杂性，我们会减少模型偏差，但我们也会增加模型方差。

19.1.3 分解

总结一下：

• 模型风险，$\mathbb{E}\left[(Y(x)-\hat{Y}(x))^2\right]$，是模型的平均预测误差的平方。

• 观测方差，$\sigma^2$，是观测中随机噪声的方差。它描述了每个观测中随机误差 $\epsilon$ 的变化程度。

• 模型偏差，$\mathbb{E}\left[\hat{Y}(x)\right]-g(x)$，是 $\hat{Y}(x)$ 作为真实基本关系 $g(x)$ 估计量的“偏离”程度。

• 模型方差，$\text{Var}(\hat{Y}(x))$，描述了当我们在不同样本上拟合模型时，预测 $\hat{Y}(x)$ 倾向于变化的程度。

上述定义使我们能够在之前简化模型风险的分解：

$E[(Y(x) - \hat{Y}(x))^2] = \sigma^2 + (E[\hat{Y}(x)] - g(x))^2 + \text{Var}(\hat{Y}(x))$

$\text{模型风险 } = \text{观测方差} + (\text{模型偏差})^2 \text{+ 模型方差}$

这被称为偏差-方差权衡。这是什么意思？记住，模型风险是模型性能的一个度量。我们建模的目标是保持模型风险低；这意味着我们希望确保模型风险的每个组成部分都保持在一个小的值。

观测方差是数据收集过程中固有的随机部分。我们无法减少观测方差，所以我们将把注意力集中在模型偏差和模型方差上。

在特征工程讲座中，我们考虑了过拟合的问题。我们发现，随着模型复杂度的增加，模型的误差或偏差往往会减少 - 如果我们设计一个非常复杂的模型，它往往会倾向于做出更接近真实关系$g$的预测。与此同时，模型方差往往会增加随着模型复杂度的增加；复杂模型可能会对训练数据过拟合，这意味着用于训练的随机样本的微小差异会导致拟合模型的巨大差异。我们有一个问题。为了减少模型偏差，我们可以增加模型的复杂度，这将导致过拟合和模型方差的增加。或者，我们可以通过减少模型的复杂度来减少模型方差，但这会增加由于欠拟合而产生的模型偏差。

我们需要取得平衡。我们在模型创建中的目标是使用足够高的复杂度水平来保持偏差低，但不要太高以至于模型方差很大。

19.2 解释回归系数

回想一下我们在本讲座中建立的框架。如果我们假设观察值和输入特征之间的潜在关系是线性的，我们可以用未知的真实模型参数$\theta$来表达这种关系。

$f_{\theta}(x) = g(x) + \epsilon = \theta_0 + \theta_1 x_1 + \ldots + \theta_p x_p + \epsilon$

我们的模型试图使用从设计矩阵$\Bbb{X}$和响应向量$\Bbb{Y}$计算出的估计值$\hat{\theta}_i$来估计每个真实参数$\theta_i$。

$f_{\hat{\theta}}(x) = \hat{\theta}_0 + \hat{\theta}_1 x_1 + \ldots + \hat{\theta}_p x_p$

让我们暂停一下。在这一点上，我们非常习惯于使用模型参数的概念。但是每个系数$\theta_i$实际上意味着什么呢？我们可以将每个$\theta_i$看作线性模型的斜率 - 如果所有其他变量保持不变，$x_i$的单位变化将导致$f_{\theta}(x)$中的$\theta_i$变化。广义上讲，$\theta_i$的值越大，意味着特征$x_i$对响应的影响越大；相反，$\theta_i$的值越小，意味着$x_i$对响应的影响越小。在极端情况下，如果真实参数$\theta_i$为 0，则特征$x_i$对$Y(x)$没有影响。

如果某个特定特征的真实参数$\theta_i$为 0，这告诉我们一些非常重要的事情：$x_i$和$Y(x)$之间没有潜在关系！那么，我们如何测试参数是否为 0 呢？作为基线，我们按照通常的流程抽取样本，使用这些数据拟合模型，并计算估计值$\hat{\theta}_i$。然而，我们还需要考虑这样一个事实，即如果我们的随机样本结果不同，我们可能会得到不同的$\hat{\theta}_i$结果。为了推断真实参数$\theta_i$是否为 0，我们希望从我们可能在所有其他随机样本中抽取的$\hat{\theta}_i$估计的分布中得出结论。这就是假设检验派上用场的地方！

测试真实参数$\theta_i$是否为 0，我们构建一个假设检验，其中零假设表明真实参数$\theta_i$为 0，备择假设表明真实参数$\theta_i$ 不是 0。如果我们的 p 值小于我们的截断值（通常 p=0.05），我们拒绝零假设。

19.3 通过 Bootstrap 进行假设检验：PurpleAir 演示

执行上述假设检验的一个等价方法是通过自举（可以通过对偶论证证明这种等价性，这超出了本课程的范围）。我们使用自举来计算每个$\theta_i$的近似 95%置信区间。如果区间不包含 0，我们在 5%的水平上拒绝零假设。否则，数据与零假设一致，因为真实参数可能为 0。

为了展示这个假设检验过程的一个例子，我们将在本节中使用雪鸻数据集。这些数据是关于雪鸻的蛋和新孵出的雏鸟。这些数据是由伯克利的一位前学生在雷耶斯角国家海岸收集的。这是一个父母鸟和一些蛋。

请注意，蛋长和蛋宽（最宽直径）以毫米为单位测量，蛋重和鸟重以克为单位测量；作为比较，一个标准的回形针重约一克。

代码

# import numpy as np
# import pandas as pd
# import matplotlib
# import matplotlib.pyplot as plt
# import seaborn as sns
# import sklearn.linear_model as lm
# from sklearn.linear_model import LinearRegression

# # big font helper
# def adjust_fontsize(size=None):
#     SMALL_SIZE = 8
#     MEDIUM_SIZE = 10
#     BIGGER_SIZE = 12
#     if size != None:
#         SMALL_SIZE = MEDIUM_SIZE = BIGGER_SIZE = size

#     plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
#     plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
#     plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
#     plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
#     plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
#     plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
#     plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

# plt.style.use('fivethirtyeight')
# sns.set_context("talk")
# sns.set_theme()
# #plt.style.use('default') # revert style to default mpl
# adjust_fontsize(size=20)
# %matplotlib inline
# csv_file = 'data/Full24hrdataset.csv'
# usecols = ['Date', 'ID', 'region', 'PM25FM', 'PM25cf1', 'TempC', 'RH', 'Dewpoint']
# full_df = (pd.read_csv(csv_file, usecols=usecols, parse_dates=['Date'])
#         .dropna())
# full_df.columns = ['date', 'id', 'region', 'pm25aqs', 'pm25pa', 'temp', 'rh', 'dew']
# full_df = full_df.loc[(full_df['pm25aqs'] < 50)]

# bad_dates = ['2019-08-21', '2019-08-22', '2019-09-24']
# GA = full_df.loc[(full_df['id'] == 'GA1') & (~full_df['date'].isin(bad_dates)) , :]
# AQS, PA = GA[['pm25aqs']], GA['pm25pa']
# AQS.head()
# pd.DataFrame(PA).head()

import pandas as pd
eggs = pd.read_csv("data/snowy_plover.csv")
eggs.head(5)

	egg_weight	egg_length	egg_breadth	bird_weight
0	7.4	28.80	21.84	5.2
1	7.7	29.04	22.45	5.4
2	7.9	29.36	22.48	5.6
3	7.5	30.10	21.71	5.3
4	8.3	30.17	22.75	5.9

我们的目标是预测新生雪鸻雏鸟的重量，我们假设其遵循下面的真实关系$Y = f_{\theta}(x)$。

$\text{bird\_weight} = \theta_0 + \theta_1 \text{egg\_weight} + \theta_2 \text{egg\_length} + \theta_3 \text{egg\_breadth} + \epsilon$

• 对于每个$i$，参数$\theta_i$是一个固定的数字，但是不可观测的。我们只能估计它。

• 随机误差$\epsilon$也是不可观测的，但假定其期望为 0，并且在蛋中是独立且同分布的。

假设我们希望确定蛋重是否影响雏鸟的鸟重-我们想推断$\theta_1$是否等于 0。

首先，我们定义我们的假设：

• 零假设：真实参数$\theta_1$为 0；任何变化都是由随机机会引起的。

• 备择假设：真实参数$\theta_1$不为 0。

接下来，我们使用我们的数据来拟合一个模型$\hat{Y} = f_{\hat{\theta}}(x)$，该模型近似上面的关系。这给我们了$\hat{\theta}_1$的观察值，从我们的数据中找到。

from sklearn.linear_model import LinearRegression
import numpy as np

X = eggs[["egg_weight", "egg_length", "egg_breadth"]]
Y = eggs["bird_weight"]

model = LinearRegression()
model.fit(X, Y)

# This gives an array containing the fitted model parameter estimates
thetas = model.coef_

# Put the parameter estimates in a nice table for viewing
display(pd.DataFrame([model.intercept_] + list(model.coef_),
 columns=['theta_hat'],
 index=['intercept', 'egg_weight', 'egg_length', 'egg_breadth']))

print("RMSE", np.mean((Y - model.predict(X)) ** 2))

	$\hat{\theta}$
intercept	-4.605670
egg_weight	0.431229
egg_length	0.066570
egg_breadth	0.215914

RMSE 0.04547085380275766

现在我们有了$\hat{\theta}_1$的值，考虑到我们拥有的单个数据样本。为了了解如果我们抽取不同的随机样本，这个估计可能会如何变化，我们将使用**自举**。为了构建一个自举样本，我们将从收集到的数据中抽取一个重采样：

• 具有与收集到的数据相同的样本大小

• 用替换的方式抽取（这确保我们不会每次抽取完全相同的样本！）

我们抽取一个自举样本，使用这个样本来拟合一个模型，并记录在这个自举样本上的$\hat{\theta}_1$的结果。然后我们重复这个过程很多次，以生成$\hat{\theta}_1$的自举经验分布。这给我们一个估计，即真实分布$\hat{\theta}_1$在所有可能的样本中可能是什么样子。

# Set a random seed so you generate the same random sample as staff
# In the "real world", we wouldn't do this
import numpy as np
np.random.seed(1337)

# Set the sample size of each bootstrap sample
n = len(eggs)

# Create a list to store all the bootstrapped estimates
estimates = []

# Generate a bootstrap resample from `eggs` and find an estimate for theta_1 using this sample. 
# Repeat 10000 times.
for i in range(10000):
 bootstrap_resample = eggs.sample(n, replace=True)
 X_bootstrap = bootstrap_resample[["egg_weight", "egg_length", "egg_breadth"]]
 Y_bootstrap = bootstrap_resample["bird_weight"]

 bootstrap_model = LinearRegression()
 bootstrap_model.fit(X_bootstrap, Y_bootstrap)
 bootstrap_thetas = bootstrap_model.coef_

 estimates.append(bootstrap_thetas[0])

# calculate the 95% confidence interval 
lower = np.percentile(estimates, 2.5, axis=0)
upper = np.percentile(estimates, 97.5, axis=0)
conf_interval = (lower, upper)
conf_interval

(-0.25864811956848754, 1.1034243854204049)

我们发现，我们的自举近似 95%置信区间为$[-0.259, 1.103]$。立即可以看到 0 确实包含在这个区间内 - 这意味着我们无法断定$\theta_1$不为零！更正式地说，我们未能拒绝零假设（即$\theta_1$为 0）在 5%的 p 值截断下。

19.4 共线性

我们可以重复这个过程，为模型的其他参数构建 95%置信区间。

代码

np.random.seed(1337)

theta_0_estimates = []
theta_1_estimates = []
theta_2_estimates = []
theta_3_estimates = []

for i in range(10000):
 bootstrap_resample = eggs.sample(n, replace=True)
 X_bootstrap = bootstrap_resample[["egg_weight", "egg_length", "egg_breadth"]]
 Y_bootstrap = bootstrap_resample["bird_weight"]

 bootstrap_model = LinearRegression()
 bootstrap_model.fit(X_bootstrap, Y_bootstrap)
 bootstrap_theta_0 = bootstrap_model.intercept_
 bootstrap_theta_1, bootstrap_theta_2, bootstrap_theta_3 = bootstrap_model.coef_

 theta_0_estimates.append(bootstrap_theta_0)
 theta_1_estimates.append(bootstrap_theta_1)
 theta_2_estimates.append(bootstrap_theta_2)
 theta_3_estimates.append(bootstrap_theta_3)

theta_0_lower, theta_0_upper = np.percentile(theta_0_estimates, 2.5), np.percentile(theta_0_estimates, 97.5)
theta_1_lower, theta_1_upper = np.percentile(theta_1_estimates, 2.5), np.percentile(theta_1_estimates, 97.5)
theta_2_lower, theta_2_upper = np.percentile(theta_2_estimates, 2.5), np.percentile(theta_2_estimates, 97.5)
theta_3_lower, theta_3_upper = np.percentile(theta_3_estimates, 2.5), np.percentile(theta_3_estimates, 97.5)

# Make a nice table to view results
pd.DataFrame({"lower":[theta_0_lower, theta_1_lower, theta_2_lower, theta_3_lower], "upper":[theta_0_upper, \
 theta_1_upper, theta_2_upper, theta_3_upper]}, index=["theta_0", "theta_1", "theta_2", "theta_3"])

	lower	upper
theta_0	-15.278542	5.161473
theta_1	-0.258648	1.103424
theta_2	-0.099138	0.208557
theta_3	-0.257141	0.758155

这里有些不对劲。注意到在模型的每个参数的 95%置信区间中都包含了 0。根据我们上面概述的解释，这意味着我们无法确定任何输入变量对响应变量的影响！这似乎表明我们的模型无法进行任何预测 - 然而，我们在上面的自助法实验中拟合的每个模型都可以非常好地预测$Y$。

我们如何解释这个结果？回想一下我们如何解释线性模型的参数。我们将每个$\theta_i$都视为斜率，其中$x_i$的单位增加导致$Y$的$\theta_i$增加，如果所有其他变量保持不变。事实证明，最后这个假设非常重要。如果我们模型中的变量某种程度上相关，那么在保持其他变量不变的情况下可能不可能改变其中一个变量。这意味着我们的解释框架不再有效！在我们上面拟合的模型中，我们将egg_length、egg_breadth和egg_weight作为输入变量。这些变量很可能彼此相关 - 一个具有较大egg_length和egg_breadth的蛋很可能在egg_weight上很重。这意味着模型参数不能被有意义地解释为斜率。

为了支持这个结论，我们可以可视化我们的特征变量之间的关系。注意特征之间的强正相关。

代码

import seaborn as sns
sns.pairplot(eggs[["egg_length", "egg_breadth", "egg_weight", 'bird_weight']]);

/Users/Ishani/micromamba/lib/python3.9/site-packages/seaborn/axisgrid.py:118: UserWarning:

The figure layout has changed to tight 

这个问题被称为共线性，有时也被称为多重共线性。当一个特征与其他特征高度相关时，就会发生共线性，这意味着一个特征可以被其他特征的线性组合相当准确地预测。

为什么共线性是一个问题？它的后果涵盖了建模过程的几个方面：

• 推断：斜率不能用于推断任务。

• 模型方差：如果特征彼此强烈影响，即使在采样数据中进行微小的变化也可能导致估计斜率的大幅变化。

• 唯一解：如果一个特征是其他特征的线性组合，设计矩阵将不是满秩的，$\mathbb{X}^{\top}\mathbb{X}$就不可逆。这意味着最小二乘法没有唯一解。

重点是，我们需要小心选择建模的特征。如果两个特征很可能编码相似的信息，通常最好只选择其中一个作为输入变量。

19.4.1 一个更简单的模型

让我们现在考虑一个更易解释的模型：我们假设真实关系只使用蛋重：

$f_\theta(x) = \theta_0 + \theta_1 \text{egg\_weight} + \epsilon$

代码

from sklearn.linear_model import LinearRegression
X_int = eggs[["egg_weight"]]
Y_int = eggs["bird_weight"]

model_int = LinearRegression()

model_int.fit(X_int, Y_int)

# This gives an array containing the fitted model parameter estimates
thetas_int = model_int.coef_

# Put the parameter estimates in a nice table for viewing
pd.DataFrame({"theta_hat":[model_int.intercept_, thetas_int[0]]}, index=["theta_0", "theta_1"])

	theta_hat
theta_0	-0.058272
theta_1	0.718515

import matplotlib.pyplot as plt

# Set a random seed so you generate the same random sample as staff
# In the "real world", we wouldn't do this
np.random.seed(1337)

# Set the sample size of each bootstrap sample
n = len(eggs)

# Create a list to store all the bootstrapped estimates
estimates_int = []

# Generate a bootstrap resample from `eggs` and find an estimate for theta_1 using this sample. 
# Repeat 10000 times.
for i in range(10000):
 bootstrap_resample_int = eggs.sample(n, replace=True)
 X_bootstrap_int = bootstrap_resample_int[["egg_weight"]]
 Y_bootstrap_int = bootstrap_resample_int["bird_weight"]

 bootstrap_model_int = LinearRegression()
 bootstrap_model_int.fit(X_bootstrap_int, Y_bootstrap_int)
 bootstrap_thetas_int = bootstrap_model_int.coef_

 estimates_int.append(bootstrap_thetas_int[0])

plt.figure(dpi=120)
sns.histplot(estimates_int, stat="density")
plt.xlabel(r"$\hat{\theta}_1$")
plt.title(r"Bootstrapped estimates $\hat{\theta}_1$ Under the Interpretable Model");

注意可解释模型的表现几乎与我们的其他模型一样好：

代码

from sklearn.metrics import mean_squared_error

rmse = mean_squared_error(Y, model.predict(X))
rmse_int = mean_squared_error(Y_int, model_int.predict(X_int))
print(f'RMSE of Original Model: {rmse}')
print(f'RMSE of Interpretable Model: {rmse_int}')

RMSE of Original Model: 0.04547085380275766
RMSE of Interpretable Model: 0.046493941375556846

然而，真实参数$\theta_{1}$的置信区间不包含零。

代码

lower_int = np.percentile(estimates_int, 2.5)
upper_int = np.percentile(estimates_int, 97.5)

conf_interval_int = (lower_int, upper_int)
conf_interval_int

(0.6029335250209633, 0.8208401738546206)

回顾来看，母鸡的重量最能预测新生小鸡的重量，这并不奇怪。

高度相关变量的模型会阻止我们解释变量与预测之间的关系。

19.4.2 提醒：假设很重要

请记住：所有推断都假设回归模型成立。

• 如果模型不成立，推断可能无效。

• 如果自助法的假设不成立…

  • 样本量 n 很大

  • 样本代表人口分布（随机抽取，无偏）…那么自助法的结果可能无效。

19.5（奖励）相关性与因果

让我们考虑一个任意回归问题中的一些问题。

在我们的回归中，$\theta_{j}$代表什么？

• 在保持其他变量不变的情况下，我们的预测应该随着$X_{j}$的变化而变化多少？

对于简单线性回归，这归结为相关系数

• 拥有更多的$x$是否预测更多的$y$（以及预测多少）？

例子：

• 拥有花岗岩台面的房屋是否更值钱？

• 获得某项奖学金的学生的大学 GPA 是否更高？

• 母乳喂养的婴儿更不容易患哮喘吗？

• 接受某种激进治疗的癌症患者是否有更高的 5 年生存率？

• 吸烟的人更容易患癌症吗？

这些听起来像是因果问题，但实际上并不是！

19.5.1 预测与因果

相关性/预测与因果之间的区别最好通过例子来说明。

一些关于相关性/预测的问题包括：

• 拥有花岗岩台面的房屋是否更值钱？

• 获得某项奖学金的学生的大学 GPA 是否更高？

• 母乳喂养的婴儿更不容易患哮喘吗？

• 接受某种激进治疗的癌症患者是否有更高的 5 年生存率？

• 吸烟的人更容易患癌症吗？

一些关于因果关系的问题包括：

• 花岗岩台面会提高房屋的价值多少？

• 获得奖学金是否提高了学生的 GPA？

• 母乳喂养是否保护婴儿免受哮喘？

• 治疗是否提高了癌症的生存率？

• 吸烟是否导致癌症？

因果问题涉及干预的效果（不仅仅是被动观察）。但需要注意的是，回归系数有时被称为“效果”，这可能是误导性的！

仅使用数据时，可以回答预测性问题（即母乳喂养的婴儿更健康吗？），但无法回答因果问题（即母乳喂养是否改善了婴儿的健康？）。原因在于我们的预测问题有许多可能的原因。例如，母乳喂养的婴儿平均更健康的可能解释包括：

1. 因果效应： 母乳喂养使婴儿更健康

2. 逆因果关系： 更健康的婴儿更有可能成功母乳喂养

3. 共同原因： 更健康/更富有的父母有更健康的婴儿，并更有可能母乳喂养

我们不能仅通过观察（$x$,$y$）对来判断哪些解释是真实的（或在多大程度上是真实的）。

此外，因果问题隐含地涉及反事实，即未发生的事件。例如，我们可以问，如果母乳喂养的婴儿没有被母乳喂养，他们是否会更健康？上面的解释 1 意味着他们会更健康，但解释 2 和 3 则不是。

19.5.2 混杂因素

让 T 代表一种治疗（例如，饮酒），Y 代表一个结果（例如，肺癌）。

混杂变量是影响 T 和 Y 的变量，扭曲它们之间的相关性。使用上面的例子。混杂因素可以是一个已测量的协变量或者是我们不知道的未测量变量，它们通常会引起问题，因为 T 和 Y 之间的关系实际上受到我们看不到的数据的影响。

**常见假设：**所有混杂因素都是被观察到的（可忽略性）

19.5.3 术语

让我们定义一些术语，这些术语将帮助我们理解因果效应。

在预测中，我们有两种变量：

• 响应 ($Y$)：我们试图预测的内容

• 预测变量 ($X$)：我们预测的输入

因果推断中的其他变量包括：

• 响应 ($Y$)：我们感兴趣的结果

• 处理 ($T$)：我们可能进行干预的变量

• 协变量 ($X$)：我们测量的其他可能影响$T$和/或$Y$的变量

对于本讲座，$T$是一个**二元（0/1）**变量：

19.5.4 Neyman-Rubin 因果模型

因果问题涉及反事实：

• 如果 T 不同会发生什么？

• 如果我们将来设定 T 会发生什么？

我们假设每个个体都有两个潜在结果：

• $Y_{i}(1)$：如果$T_{i} = 1$的话，$y_{i}$的值（受治疗结果）

• $Y_{i}(0)$：如果$T_{i} = 0$的话，$y_{i}$的值（对照结果）

对于数据集中的每个个体，我们观察到：

• 协变量$x_{i}$

• 处理$T_{i}$

• 响应$y_{i} = Y_{i}(T_{i})$

我们假设对于$i = 1,..., n$，($x_{i}$, $T_{i}$, $y_{i} = Y_{i}(T_{i})$)元组是独立同分布的

19.5.5 平均处理效应

对于每个个体，处理效应是$Y_{i}(1)-Y_{i}(0)$

最常见的估计是平均处理效应（ATE）

$ATE = \mathbb{E}[Y(1)-Y(0)] = \mathbb{E}[Y(1)] - \mathbb{E}[Y(0)]$

我们能否只取样本均值？

$\hat{ATE} = \frac{1}{n}\sum_{i=1}^{n}Y_{i}(1) - Y_{i}(0)$

我们不能。为什么？我们只观察到$Y_{i}(1)$，$Y_{i}(0)$中的一个。

**因果推断的基本问题：**我们只能观察到一个潜在结果

要得出因果结论，我们需要一些关于观察到的和未观察到的单位之间的因果假设

与其$\frac{1}{n}\sum_{i=1}^{n}Y_{i}(1) - Y_{i}(0)$，不如我们取每组的样本均值之间的差异？

$\hat{ATE} = \frac{1}{n_{1}}\sum_{i: T_{i} = 1}{Y_{i}(1)} - \frac{1}{n_{0}}\sum_{i: T_{i} = 0}{Y_{i}(0)} = \frac{1}{n_{1}}\sum_{i: T_{i} = 1}{y_{i}} - \frac{1}{n_{0}}\sum_{i: T_{i} = 0}{y_{i}}$

这个$ATE$的估计是否无偏？因此，这个提出的$\hat{ATE}$不适合我们的目的。

如果处理分配来自随机抛硬币，那么受治疗单位是来自$Y_{i}(1)$总体的大小为$n_{1}$的独立同分布随机样本。

这意味着，

$\mathbb{E}[\frac{1}{n_{1}}\sum_{i: T_{i} = 1}{y_{i}}] = \mathbb{E}[Y_{i}(1)]$

同样地，

$\mathbb{E}[\frac{1}{n_{0}}\sum_{i: T_{i} = 0}{y_{i}}] = \mathbb{E}[Y_{i}(0)]$

这使我们得出结论，$\hat{ATE}$是$ATE$的无偏估计：

$\mathbb{E}[\hat{ATE}] = ATE$

19.5.6 随机实验

然而，通常随机分配处理是不切实际或不道德的。例如，分配香烟治疗可能是不切实际和不道德的。

绕过这个问题的另一种方法是利用观测研究。

实验：

观测研究：

19.5.7 协变量调整

对于混杂因素该怎么办？

• **可忽略性假设：**所有重要的混杂因素都在数据集中！

**一个想法：**提出一个包括它们的模型，比如：

$Y_{i}(t) = \theta_{0} + \theta_{1}x_{1} + ... + \theta_{p}x_{p} + \tau{t} + \epsilon$

**问题：**在这个模型中$ATE$是多少？$\tau$

这种方法可能有效，但是脆弱。如果：

• 重要的协变量缺失或者对$x$的真实依赖是非线性的

• 有时被贬低地称为**“因果推断”**

19.5.7.1 不需要参数假设的协变量调整

对混杂因素怎么办？

• **可忽略性假设：**数据集中包含所有可能的混杂因素！

**一个想法：**提出一个包括它们的模型，比如：

$Y_{i}(t) = f_{\theta}(x, t) + \epsilon$

然后：

$ATE = \frac{1}{n}\sum_{i=1}^{n}{f_{\theta}(x_i, 1) - f_{\theta}(x_i, 0)}$

有了足够的数据，我们可能能够非常准确地学习$f_{\theta}$

• 如果$x$是高维的/其函数形式高度非线性，则非常困难

• 需要额外的假设：重叠

19.5.8 其他方法

因果推断很难，协变量调整通常不是最佳方法

许多其他方法是一些组合：

• 将处理 T 建模为协变量 x 的函数

• 将结果 y 建模为 x，T 的函数

如果我们不相信可忽略性呢？其他方法寻找一个

• 最喜欢的例子：回归不连续

19.6（奖励）偏差-方差分解的证明

本节详细推导了偏差-方差分解在前面笔记中的偏差-方差权衡部分。

点击显示

我们想证明模型风险可以分解为

\begin{align*} E\left[(Y(x)-\hat{Y}(x))^2\right] &= E[\epsilon^2] + \left(g(x)-E\left[\hat{Y}(x)\right]\right)^2 + E\left[\left(E\left[\hat{Y}(x)\right] - \hat{Y}(x)\right)^2\right]. \end{align*}

为了证明这一点，我们首先需要以下引理：

如果$V$和$W$是独立的随机变量，则$E[VW] = E[V]E[W]$。

我们将在离散有限的情况下证明这一点。相信它在更广泛的情况下也是成立的。

工作是计算$VW$值的加权平均值，其中权重是这些值的概率。开始吧。

\begin{align*} E[VW] ~ &= ~ \sum_v\sum_w vwP(V=v \text{ and } W=w) \\ &= ~ \sum_v\sum_w vwP(V=v)P(W=w) ~~~~ \text{根据独立性} \\ &= ~ \sum_v vP(V=v)\sum_w wP(W=w) \\ &= ~ E[V]E[W] \end{align*}

现在我们进入实际的证明：

19.6.1 目标

将模型风险分解为可识别的组成部分。

19.6.2 步骤 1

\begin{align*} \text{模型风险} ~ &= ~ E\left[\left(Y - \hat{Y}(x)\right)^2 \right] \\ &= ~ E\left[\left(g(x) + \epsilon - \hat{Y}(x)\right)^2 \right] \\ &= ~ E\left[\left(\epsilon + \left(g(x)- \hat{Y}(x)\right)\right)^2 \right] \\ &= ~ E\left[\epsilon^2\right] + 2E\left[\epsilon \left(g(x)- \hat{Y}(x)\right)\right] + E\left[\left(g(x) - \hat{Y}(x)\right)^2\right]\\ \end{align*}

右边：

• 第一项是观测方差$\sigma^2$。

• 交叉乘积项为 0，因为$\epsilon$与$g(x) - \hat{Y}(x)$独立，且$E(\epsilon) = 0$

• 最后一项是我们预测值与$x$处真实函数值之间的均方差差

19.6.3 步骤 2

到这个阶段我们有

$\text{模型风险} ~ = ~ E\left[\epsilon^2\right] + E\left[\left(g(x) - \hat{Y}(x)\right)^2\right]$

我们还不太了解$g(x) - \hat{Y}(x)$。但我们了解偏差$D_{\hat{Y}(x)} = \hat{Y}(x) - E\left[\hat{Y}(x)\right]$。我们知道

• $E\left[D_{\hat{Y}(x)}\right] ~ = ~ 0$

• $E\left[D_{\hat{Y}(x)}^2\right] ~ = ~ \text{模型方差}$

因此，让我们添加并减去$E\left[\hat{Y}(x)\right]，看看是否有帮助。

$g(x) - \hat{Y}(x) ~ = ~ \left(g(x) - E\left[\hat{Y}(x)\right] \right) + \left(E\left[\hat{Y}(x)\right] - \hat{Y}(x)\right)$

右边的第一项是$x$处的模型偏差。第二项是$-D_{\hat{Y}(x)}$。所以

$g(x) - \hat{Y}(x) ~ = ~ \text{模型偏差} - D_{\hat{Y}(x)}$

19.6.4 步骤 3

记住，在$x$处的模型偏差是一个常数，不是一个随机变量。把它看作你最喜欢的数字，比如 10。那么 \begin{align*} E\left[ \left(g(x) - \hat{Y}(x)\right)^2 \right] ~ &= ~ \text{模型偏差}^2 - 2(\text{模型偏差})E\left[D_{\hat{Y}(x)}\right] + E\left[D_{\hat{Y}(x)}^2\right] \\ &= ~ \text{模型偏差}^2 - 0 + \text{模型方差} \\ &= ~ \text{模型偏差}^2 + \text{模型方差} \end{align*}

同样，交叉乘积项为 0，因为$E\left[D_{\hat{Y}(x)}\right] ~ = ~ 0$。

19.6.5 第 4 步：偏差-方差分解

在第 2 步中我们有

$\text{模型风险} ~ = ~ \text{观测方差} + E\left[\left(g(x) - \hat{Y}(x)\right)^2\right]$

第 3 步显示

$E\left[ \left(g(x) - \hat{Y}(x)\right)^2 \right] ~ = ~ \text{模型偏差}^2 + \text{模型方差}$

因此，我们已经展示了偏差-方差分解：

$\text{模型风险} = \text{观测方差} + \text{模型偏差}^2 + \text{模型方差}。$

也就是说，

$E\left[(Y(x)-\hat{Y}(x))^2\right] = \sigma^2 + \left(E\left[\hat{Y}(x)\right] - g(x)\right)^2 + E\left[\left(\hat{Y}(x)-E\left[\hat{Y}(x)\right]\right)^2\right]$

二十、SQL I

原文：SQL I

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 确定数据库可能优于 CSV 文件的情况

• 使用SELECT，FROM，WHERE，ORDER BY，LIMIT和OFFSET编写基本的 SQL 查询

• 使用GROUP BY执行聚合操作

到目前为止，在本课程中，我们已经完成了整个数据科学生命周期：我们学会了如何加载和探索数据集，制定问题，并使用预测和推断工具得出答案。在本学期的剩余几周中，我们将再次经历整个生命周期，这次使用不同的工具、思想和抽象。

20.1 数据库

有了这个目标，让我们回到生命周期的最开始。我们首先通过查看pandas库开始了数据分析工作，该库为我们提供了强大的工具，用于操作（主要是）CSV 文件中存储的表格数据。在研究和工业领域，数据科学家经常需要访问无法轻松存储在 CSV 格式中的大量数据。与他人合作处理 CSV 数据也可能会很棘手-真实世界的数据科学家可能会在多个用户尝试进行修改时遇到问题，或者更糟的是，数据应该由谁访问和谁不应该访问的安全问题。

数据库是一个大型的、有组织的数据集合。数据库由**数据库管理系统（DBMS）**管理，这是一种存储、管理和促进访问一个或多个数据库的软件系统。数据库有助于减轻使用 CSV 进行数据存储时出现的许多问题：它们提供可靠的存储，可以在系统崩溃或磁盘故障时幸存，被优化用于计算无法适应内存的数据，并包含特殊的数据结构以提高性能。使用数据库而不是 CSV 还可以从数据管理的角度获得进一步的好处。DBMS 可以应用设置来配置数据的组织方式，阻止某些数据异常（例如，强制执行非负权重或年龄），并确定谁有权访问数据。它还可以确保安全的并发操作，其中多个用户读取和写入数据库不会导致致命错误。

正如您可能已经猜到的那样，我们无法使用通常的pandas方法来处理数据库中的数据。相反，我们将转向结构化查询语言。

20.2 结构化查询语言和数据库模式

结构化查询语言，或SQL（通常发音为“sequel”，尽管这是激烈辩论的话题），是一种专门设计用于与数据库通信的编程语言。您可能在 CS 61A 或 Data C88C 等课程中遇到过它。与 Python 不同，它是一种声明性编程语言 - 这意味着与编写完成任务所需的确切逻辑不同，SQL 代码的一部分“声明”了所需的最终输出应该是什么，并且让程序确定应该实现什么逻辑。

重申一点，SQL 是一种与 Python 完全不同的语言。但是，Python 确实有特殊的引擎，允许我们在 Jupyter 笔记本中运行 SQL 代码。虽然这通常不是 SQL 在教育环境之外的使用方式，但我们将使用这种工作流程来说明如何使用我们本学期已经使用过的工具构建 SQL 查询。您将在实验 10 中了解如何在 Jupyter 中运行 SQL 查询。

下面的语法对您来说可能会很陌生；现在，只需专注于理解显示的输出。我们将稍后澄清 SQL 代码。

首先，我们将查看一个名为basic_examples.db的数据库。

# Load the SQL Alchemy Python library
import sqlalchemy
import pandas as pd

# load %%sql cell magic
%load_ext sql

连接到 SQLite 数据库basic_examples.db。

%%sql
sqlite:///data/basic_examples.db 

%%sql
SELECT * 
FROM sqlite_master
WHERE type="table"

 * sqlite:///data/basic_examples.db
Done.

type	Name	tbl_name	rootpage	sql
table	sqlite_sequence	sqlite_sequence	7	CREATE TABLE sqlite_sequence(name,seq)
table	Dragon	Dragon	2	CREATE TABLE Dragon (
				name TEXT PRIMARY KEY,
				year INTEGER CHECK (year >= 2000),
				cute INTEGER
				)
table	Dish	Dish	4	CREATE TABLE Dish (
				name TEXT PRIMARY KEY,
				type TEXT,
				cost INTEGER CHECK (cost >= 0)
				)
table	Scene	Scene	6	CREATE TABLE Scene (
				id INTEGER PRIMARY KEY AUTOINCREMENT,
				biome TEXT NOT NULL,
				city TEXT NOT NULL,
				visitors INTEGER CHECK (visitors >= 0),
				created_at DATETIME DEFAULT (DATETIME('now'))
				)

上面的摘要显示了有关数据库的信息。数据库包含四个表，名称分别为sqlite_sequence、Dragon、Dish和Scene。上面最右边的列列出了用于构造每个表的命令。

让我们更仔细地看一下用于创建Dragon表的命令（上面的第二个条目）。

CREATE TABLE Dragon (name TEXT PRIMARY KEY,
                     year INTEGER CHECK (year >= 2000),
                     cute INTEGER)

语句“CREATE TABLE”用于指定表的模式-表的组织逻辑的描述。模式遵循一组格式：

• “ColName”：列的名称

• “DataType”：要存储在列中的数据类型。一些最常见的 SQL 数据类型是INT（整数）、FLOAT（浮点数）、TEXT（字符串）、BLOB（任意数据，如音频/视频文件）和DATETIME（日期和时间）。

• “约束”：对要存储在列中的数据的一些限制。常见的约束有“CHECK”（数据必须遵守某个条件）、“PRIMARY KEY”（指定列为表的主键）、“NOT NULL”（数据不能为空）和“DEFAULT”（如果没有给定特定条目，则为默认填充值）。

我们看到Dragon包含三列。其中第一列“名称”包含文本数据。它被指定为表的主键；也就是说，“名称”中包含的数据唯一标识表中的每个条目。因为“名称”是表的主键，表中的两个条目不能具有相同的名称-“名称”的给定值对于每个龙是唯一的。列“年份”包含整数数据，约束条件是年份值必须大于或等于 2000。最后一列“可爱”包含整数数据，没有对允许的值施加限制。

我们可以通过查看Dragon本身来验证这一点。

%%sql
SELECT *
FROM Dragon

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
hiccup	2010	10
drogon	2011	-100
drogon 2	2019	0

数据库表（也称为关系）的结构与pandas中的DataFrame非常相似。每一行，有时称为元组，代表数据集中的单个记录。每一列，有时称为属性或字段，描述记录的某些特征。

20.3 从表中选择

为了提取和操作存储在 SQL 表中的数据，我们需要熟悉编写 SQL 代码片段的语法，我们称之为查询。

SQL 查询的基本单元是“SELECT”语句。“SELECT”指定我们想要从给定表中提取哪些列。我们使用“FROM”告诉 SQL 我们想要从哪个表中“SELECT”我们的数据。

%%sql
SELECT *
FROM Dragon

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
hiccup	2010	10
drogon	2011	-100
drogon 2	2019	0

在 SQL 中，*表示“所有”。上面的查询抓取Dragon中的所有列，并在输出的表中显示它们。我们也可以指定要SELECT的特定列的子集。请注意，输出的列按照它们被SELECT的顺序出现。

%%sql
SELECT cute, year
FROM Dragon

 * sqlite:///data/basic_examples.db
Done.

Cute	Year
10	2010
-100	2011
0	2019

就像这样，我们已经编写了两个 SQL 查询。上面的查询中有一些要注意的地方。首先，请注意每个“动词”都是大写写的。按照惯例，SQL 操作应以大写字母编写，但即使您选择保持小写，您的代码也会正常运行。其次，上面的查询将每个语句与新行分隔开。SQL 查询不受查询内部的空格影响；这意味着 SQL 代码通常是在每个语句后写上新行，以使其更易读。分号（;）表示查询的结束。在某些 SQL“风味”中，如果没有分号，查询将无法运行；但是，在 Data 100 中，我们将使用的 SQL 版本无论是否有结束分号都可以正常运行。这些笔记中的查询将以分号结束，以养成良好的习惯。

AS关键字允许我们在SELECT后为列指定一个新名称（称为别名）。一般语法是：

SELECT column_name_in_database_table AS new_name_in_output_table

%%sql
SELECT cute AS cuteness, year AS birth
FROM Dragon

 * sqlite:///data/basic_examples.db
Done.

cuteness	birth
10	2010
-100	2011
0	2019

要仅SELECT列中的唯一值，我们使用DISTINCT关键字。这将导致列中的任何重复条目被删除。如果我们只想找到Dragon中唯一的年份，而没有任何重复，我们将写：

%%sql
SELECT DISTINCT year
FROM Dragon

 * sqlite:///data/basic_examples.db
Done.

Year
2010
2011
2019

每个 SQL 查询必须包括SELECT和FROM语句。直观地说，这是有道理的 - 我们知道我们将要从表中提取一些信息；为了这样做，我们还需要指示我们要考虑哪个表。

重要的是要注意，SQL 强制执行严格的“操作顺序” - SQL 子句必须始终遵循相同的顺序。例如，SELECT语句必须始终在FROM之前。这意味着任何 SQL 查询都将遵循相同的结构。

SELECT <column list>
FROM <table>
[additional clauses]

我们使用的附加子句取决于要实现的具体任务。我们可以通过细化查询来过滤特定条件，聚合特定列或将多个表连接在一起。我们将在本讲座的其余时间中概述一些有用的子句，以增进我们对操作顺序的理解。

20.4 应用WHERE条件

WHERE关键字用于仅选择表的某些行，这些行基于给定的布尔条件进行过滤。

%%sql
SELECT name, year
FROM Dragon
WHERE cute > 0

 * sqlite:///data/basic_examples.db
Done.

Name	Year
hiccup	2010

我们可以使用AND，OR和NOT关键字给WHERE条件添加复杂性，就像在 Python 中一样。

%%sql
SELECT name, year
FROM Dragon
WHERE cute > 0 OR year > 2013

 * sqlite:///data/basic_examples.db
Done.

Name	Year
hiccup	2010
drogon 2	2019

为了避免需要通过组合多个条件来编写复杂的逻辑表达式，我们还可以过滤IN指定值列表中的条目。这类似于 Python 中的in或.isin的用法。

%%sql
SELECT name, year
FROM Dragon
WHERE name IN ("hiccup", "puff")

 * sqlite:///data/basic_examples.db
Done.

Name	Year
hiccup	2010

您可能已经注意到我们的表实际上缺少一个值。在 SQL 中，缺失的数据被赋予特殊值NULL。NULL的行为方式与其他数据类型根本不同。我们不能对NULL值使用典型的运算符（=，>和<）（实际上，NULL == NULL返回False！）；相反，我们检查值是否为NULL或NULL。

%%sql
SELECT *
FROM Dragon
WHERE cute IS NOT NULL

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
hiccup	2010	10
drogon	2011	-100
drogon 2	2019	0

20.5 排序和限制输出

如果我们希望输出表按特定顺序显示，ORDER BY关键字的行为类似于pandas中的.sort_values()。

%%sql
SELECT *
FROM Dragon
ORDER BY cute

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
drogon	2011	-100
drogon 2	2019	0
hiccup	2010	10

默认情况下，ORDER BY将按升序显示结果（最小值在表的顶部）。要按降序排序，我们在指定用于排序的列后使用DESC关键字。

%%sql
SELECT *
FROM Dragon
ORDER BY cute DESC

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
hiccup	2010	10
drogon 2	2019	0
drogon	2011	-100

我们还可以告诉 SQL 同时按两列ORDER BY。这将按照列出的第一列对表进行排序，然后使用第二列中的值来打破任何平局。

%%sql
SELECT *
FROM Dragon
ORDER BY name, cute

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
drogon 2	2019	0
drogon	2011	-100
hiccup	2010	10

在许多情况下，我们只关心输出表中的某些行（例如，想要找到表中的前两条龙）。LIMIT关键字将输出限制为指定数量的行。它的功能类似于pandas中的.head()。

%%sql
SELECT *
FROM Dragon
LIMIT 2

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
hiccup	2010	10
drogon	2011	-100

OFFSET关键字表示LIMIT应该开始的索引。换句话说，我们可以使用OFFSET来将LIMIT的开始位置向后移动指定数量的行。例如，我们可能关心表中位置为#2 和#3 的龙。

%%sql
SELECT *
FROM Dragon
LIMIT 2
OFFSET 1

 * sqlite:///data/basic_examples.db
Done.

Name	Year	Cute
drogon	2011	-100
drogon 2	2019	0

让我们总结一下到目前为止我们学到的东西。我们知道SELECT和FROM是任何 SQL 查询的基本构建模块。我们可以使用其他子句来完善输出表中的数据。

我们包含的任何子句必须在查询中遵循严格的顺序：

SELECT <column list>
FROM <table>
[WHERE <predicate>]
[ORDER BY <column list>]
[LIMIT <number of rows>]
[OFFSET <number of rows>]

在这里，方括号[ ]中包含的任何子句都是可选的 - 只有在与我们要执行的表操作相关时才需要使用关键字。还要注意，按照惯例，我们在 SQL 语句中使用大写字母表示关键字，并使用换行符使代码更易读。

20.6 使用GROUP BY进行聚合

到目前为止，我们已经看到 SQL 提供了与pandas给我们的许多相同功能。我们可以从表中提取数据，对其进行过滤和重新排序以满足我们的需求。

在pandas中，我们的许多分析工作在很大程度上依赖于能够使用.groupby()来对数据集的行进行聚合。SQL 对这一任务的回答是（非常方便地命名为）GROUP BY子句。虽然GROUP BY的输出与.groupby()的输出类似 - 在这两种情况下，我们都获得了一个输出表，其中某些列已被用于分组 - 但在 SQL 中用于分组数据的语法和逻辑与pandas的实现有相当大的不同。

为了说明GROUP BY，我们将考虑basic_examples.db数据库中的Dish表。

%%sql
SELECT * 
FROM Dish

 * sqlite:///data/basic_examples.db
Done.

Name	type	cost
ravioli	entree	10
ramen	entree	13
taco	entree	7
edamame	appetizer	4
fries	appetizer	4
potsticker	appetizer	4
ice cream	dessert	5

假设我们想要找到某种“类型”的菜肴的总成本。为了实现这一目标，我们将编写以下代码。

%%sql
SELECT type, SUM(cost)
FROM Dish
GROUP BY type

 * sqlite:///data/basic_examples.db
Done.

type	SUM(cost)
appetizer	12
dessert	5
entree	30

这里发生了什么？语句GROUP BY type告诉 SQL 根据type列中包含的值（记录是开胃菜、主菜还是甜点）对数据进行分组。SUM(cost)对每个type中的菜肴成本进行求和，并在输出表中显示结果。

你可能会想：为什么SUM(cost)在GROUP BY type命令之前？我们不需要在计算每个分组中的条目数量之前先形成分组吗？

请记住，SQL 是一种声明性编程语言 - SQL 程序员只需陈述他们想要看到的最终结果，然后将如何获得这个结果的任务留给 SQL 本身。这意味着 SQL 查询有时不遵循读者认为“逻辑”的思维顺序。相反，SQL 要求我们在构建查询时遵循其设定的操作顺序。只要我们遵循这个顺序，SQL 就会处理底层逻辑。

在实际操作中：我们这个查询的目标是输出每个“类型”的总成本。为了向 SQL 传达这一点，我们说我们要SELECT每个type组的SUMmedcost值。

有许多聚合函数可以用来聚合每个组中包含的数据。一些常见的例子是：

• COUNT: 计算与每个组相关联的行数

• MIN: 找到每个组的最小值

• MAX: 找到每个组的最大值

• SUM: 对每个组中的所有记录求和

• AVG: 找到每个组的平均值

我们可以轻松地一次性计算多个聚合（这在pandas中非常棘手）。

%%sql
SELECT type, SUM(cost), MIN(cost), MAX(name)
FROM Dish
GROUP BY type

 * sqlite:///data/basic_examples.db
Done.

type	SUM(cost)	MIN(cost)	MAX(name)
appetizer	12	4	potsticker
dessert	5	5	ice cream
entree	30	7	taco

要计算与每个组相关联的行数，我们使用COUNT关键字。调用COUNT(*)将计算每个组中的总行数，包括具有空值的行。它在pandas中的等价物是.groupby().size()。

%%sql
SELECT type, COUNT(*)
FROM Dish
GROUP BY type

 * sqlite:///data/basic_examples.db
Done.

type	COUNT(*)
appetizer	3
dessert	1
entree	3

在计算每个组中的行数时排除NULL值，我们在表中明确调用COUNT。这类似于在pandas中调用.groupby().count()。

%%sql
SELECT year, COUNT(cute)
FROM Dragon
GROUP BY year

 * sqlite:///data/basic_examples.db
Done.

Year	COUNT(cute)
2010	1
2011	1
2019	1

有了这个GROUP BY的定义，让我们更新一下我们的 SQL 操作顺序。记住：每个SQL 查询必须按照这个顺序列出子句。

SELECT <column expression list>
FROM <table>
[WHERE <predicate>]
[GROUP BY <column list>]
[ORDER BY <column list>]
[LIMIT <number of rows>]
[OFFSET <number of rows>];

请注意，我们可以在选择过程中使用AS关键字重命名列，并且列表达式可以包括聚合函数（MAX，MIN等）。