UCB Data100：数据科学的原理和技巧：第十九章到第二十章

十九、偏差、方差和推断

学习成果

• 计算参数的估计器的偏差、方差和均方误差
  
• 介绍拟合模型的模型风险
  
• 将模型风险分解为偏差和方差项
  
• 构建假设检验的置信区间
  
• 了解我们所做的假设及其对回归推断的影响
  
• 比较回归和因果关系
  
• 实验设置、混杂变量、平均处理效应和协变量调整**上次，我们介绍了随机变量的概念及其对我们用来拟合模型的观察关系的影响。
  

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

19.1 偏差-方差权衡

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

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

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

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

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

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

19.1.1 估计器的性能

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

• 我们平均得到正确答案吗？(偏差)
  
• 答案有多大的变化？(方差)
  
• 我们离 Y YY 有多近？(风险/MSE)
  

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

19.1.2 模型风险

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

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

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

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

• 来源 1 (观测方差)：由于随机噪声 ϵ \epsilonϵ 导致新观测 Y YY 的随机性
  
• 来源 2 (模型方差)：在我们用来训练模型的样本中的随机性，因为样本 X 1 , X 2 , … , X n , Y X_1, X_2, \ldots, X_n, YX1,X2,…,Xn,Y 是随机的
  

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

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

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

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

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

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

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

Y ( x ) − Y ^ ( x ) = ϵ + ( g ( x ) − E [ Y ^ ( x ) ] ) + ( E [ Y ^ ( x ) ] − 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)Y(x)−Y^(x)=ϵ+(g(x)−E[Y^(x)])+(E[Y^(x)]−Y^(x))

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

E [ ( Y ( x ) − Y ^ ( x ) ) 2 ] = E [ ϵ 2 ] + ( g ( x ) − E [ Y ^ ( x ) ] ) 2 + E [ ( E [ Y ^ ( x ) ] − Y ^ ( x ) ) 2 ] 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]E[(Y(x)−Y^(x))2]=E[ϵ2]+(g(x)−E[Y^(x)])2+E[(E[Y^(x)]−Y^(x))2]

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

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

19.1.2.1 观察方差

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

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

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

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

19.1.2.2 模型方差

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

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

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

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

19.1.2.3 模型偏差

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

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

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

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

无偏估计

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

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

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

• 拟合不足：我们的模型对数据来说太简单。
  
• 缺乏领域知识：我们不了解哪些特征对响应变量有用
  

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

19.1.3 分解

总结一下：

• 模型风险，E [ ( Y ( x ) − Y ^ ( x ) ) 2 ] \mathbb{E}\left[(Y(x)-\hat{Y}(x))^2\right]E[(Y(x)−Y^(x))2]，是模型的平均预测误差的平方。
  
• 观测方差，σ 2 \sigma^2σ2，是观测中随机噪声的方差。它描述了每个观测中随机误差 ϵ \epsilonϵ 的变化程度。
  
• 模型偏差，E [ Y ^ ( x ) ] − g ( x ) \mathbb{E}\left[\hat{Y}(x)\right]-g(x)E[Y^(x)]−g(x)，是 Y ^ ( x ) \hat{Y}(x)Y^(x) 作为真实基本关系 g ( x ) g(x)g(x) 估计量的“偏离”程度。
  
• 模型方差，Var ( Y ^ ( x ) ) \text{Var}(\hat{Y}(x))Var(Y^(x))，描述了当我们在不同样本上拟合模型时，预测 Y ^ ( x ) \hat{Y}(x)Y^(x) 倾向于变化的程度。
  

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

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

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

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

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

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

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

19.2 解释回归系数

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

f θ ( x ) = g ( x ) + ϵ = θ 0 + θ 1 x 1 + … + θ p x p + ϵ f_{\theta}(x) = g(x) + \epsilon = \theta_0 + \theta_1 x_1 + \ldots + \theta_p x_p + \epsilonfθ(x)=g(x)+ϵ=θ0+θ1x1+…+θpxp+ϵ

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

f θ ^ ( x ) = θ ^ 0 + θ ^ 1 x 1 + … + θ ^ p x p f_{\hat{\theta}}(x) = \hat{\theta}_0 + \hat{\theta}_1 x_1 + \ldots + \hat{\theta}_p x_pfθ^(x)=θ^0+θ^1x1+…+θ^pxp

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

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

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

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

执行上述假设检验的一个等价方法是通过自举（可以通过对偶论证证明这种等价性，这超出了本课程的范围）。我们使用自举来计算每个θ i \theta_iθ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)

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

bird_weight = θ 0 + θ 1 egg_weight + θ 2 egg_length + θ 3 egg_breadth + ϵ \text{bird\_weight} = \theta_0 + \theta_1 \text{egg\_weight} + \theta_2 \text{egg\_length} + \theta_3 \text{egg\_breadth} + \epsilonbird_weight=θ0+θ1egg_weight+θ2egg_length+θ3egg_breadth+ϵ

• 对于每个i ii，参数θ i \theta_iθi是一个固定的数字，但是不可观测的。我们只能估计它。
  
• 随机误差ϵ \epsilonϵ也是不可观测的，但假定其期望为 0，并且在蛋中是独立且同分布的。
  

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

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

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

接下来，我们使用我们的数据来拟合一个模型Y ^ = f θ ^ ( x ) \hat{Y} = f_{\hat{\theta}}(x)Y^=fθ^(x)，该模型近似上面的关系。这给我们了θ ^ 1 \hat{\theta}_1θ^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))

RMSE 0.04547085380275766

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

• 具有与收集到的数据相同的样本大小
  
• 用替换的方式抽取（这确保我们不会每次抽取完全相同的样本！）
  

我们抽取一个自举样本，使用这个样本来拟合一个模型，并记录在这个自举样本上的θ ^ 1 \hat{\theta}_1θ^1的结果。然后我们重复这个过程很多次，以生成θ ^ 1 \hat{\theta}_1θ^1的自举经验分布。这给我们一个估计，即真实分布θ ^ 1 \hat{\theta}_1θ^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.259, 1.103][−0.259,1.103]。立即可以看到 0 确实包含在这个区间内 - 这意味着我们无法断定θ 1 \theta_1θ1不为零！更正式地说，我们未能拒绝零假设（即θ 1 \theta_1θ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"])

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

我们如何解释这个结果？回想一下我们如何解释线性模型的参数。我们将每个θ i \theta_iθi都视为斜率，其中x i x_ixi的单位增加导致Y YY的θ i \theta_iθ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

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

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

• 推断：斜率不能用于推断任务。
  
• 模型方差：如果特征彼此强烈影响，即使在采样数据中进行微小的变化也可能导致估计斜率的大幅变化。
  
• 唯一解：如果一个特征是其他特征的线性组合，设计矩阵将不是满秩的，X ⊤ X \mathbb{X}^{\top}\mathbb{X}X⊤X就不可逆。这意味着最小二乘法没有唯一解。
  

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

19.4.1 一个更简单的模型

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

f θ ( x ) = θ 0 + θ 1 egg_weight + ϵ f_\theta(x) = \theta_0 + \theta_1 \text{egg\_weight} + \epsilonfθ(x)=θ0+θ1egg_weight+ϵ

代码

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"])

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

然而，真实参数θ 1 \theta_{1}θ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（奖励）相关性与因果

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

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

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

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

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

例子：

• 拥有花岗岩台面的房屋是否更值钱？
  
• 获得某项奖学金的学生的大学 GPA 是否更高？
  
• 母乳喂养的婴儿更不容易患哮喘吗？
  
• 接受某种激进治疗的癌症患者是否有更高的 5 年生存率？
  
• 吸烟的人更容易患癌症吗？
  

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

19.5.1 预测与因果

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

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

• 拥有花岗岩台面的房屋是否更值钱？
  
• 获得某项奖学金的学生的大学 GPA 是否更高？
  
• 母乳喂养的婴儿更不容易患哮喘吗？
  
• 接受某种激进治疗的癌症患者是否有更高的 5 年生存率？
  
• 吸烟的人更容易患癌症吗？
  

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

• 花岗岩台面会提高房屋的价值多少？
  
• 获得奖学金是否提高了学生的 GPA？
  
• 母乳喂养是否保护婴儿免受哮喘？
  
• 治疗是否提高了癌症的生存率？
  
• 吸烟是否导致癌症？
  

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

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

1. 因果效应： 母乳喂养使婴儿更健康
   
2. 逆因果关系： 更健康的婴儿更有可能成功母乳喂养
   
3. 共同原因： 更健康/更富有的父母有更健康的婴儿，并更有可能母乳喂养
   

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

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

19.5.2 混杂因素

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

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

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

19.5.3 术语

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

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

• 响应 (Y YY)：我们试图预测的内容
  
• 预测变量 (X XX)：我们预测的输入
  

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

• 响应 (Y YY)：我们感兴趣的结果
  
• 处理 (T TT)：我们可能进行干预的变量
  
• 协变量 (X XX)：我们测量的其他可能影响T TT和/或Y YY的变量
  

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

19.5.4 Neyman-Rubin 因果模型

因果问题涉及反事实：

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

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

• Y i ( 1 ) Y_{i}(1)Yi(1)：如果T i = 1 T_{i} = 1Ti=1的话，y i y_{i}yi的值（受治疗结果）
  
• Y i ( 0 ) Y_{i}(0)Yi(0)：如果T i = 0 T_{i} = 0Ti=0的话，y i y_{i}yi的值（对照结果）
  

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

• 协变量x i x_{i}xi
  
• 处理T i T_{i}Ti
  
• 响应y i = Y i ( T i ) y_{i} = Y_{i}(T_{i})yi=Yi(Ti)
  

我们假设对于i = 1 , . . . , n i = 1,..., ni=1,...,n，(x i x_{i}xi, T i T_{i}Ti, y i = Y i ( T i ) y_{i} = Y_{i}(T_{i})yi=Yi(Ti))元组是独立同分布的

19.5.5 平均处理效应

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

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

A T E = E [ Y ( 1 ) − Y ( 0 ) ] = E [ Y ( 1 ) ] − E [ Y ( 0 ) ] ATE = \mathbb{E}[Y(1)-Y(0)] = \mathbb{E}[Y(1)] - \mathbb{E}[Y(0)]ATE=E[Y(1)−Y(0)]=E[Y(1)]−E[Y(0)]

我们能否只取样本均值？

A T E ^ = 1 n ∑ i = 1 n Y i ( 1 ) − Y i ( 0 ) \hat{ATE} = \frac{1}{n}\sum_{i=1}^{n}Y_{i}(1) - Y_{i}(0)ATE^=n1i=1∑nYi(1)−Yi(0)

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

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

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

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

A T E ^ = 1 n 1 ∑ i : T i = 1 Y i ( 1 ) − 1 n 0 ∑ i : T i = 0 Y i ( 0 ) = 1 n 1 ∑ i : T i = 1 y i − 1 n 0 ∑ i : T i = 0 y i \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^=n11i:Ti=1∑Yi(1)−n01i:Ti=0∑Yi(0)=n11i:Ti=1∑yi−n01i:Ti=0∑yi

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

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

这意味着，

E [ 1 n 1 ∑ i : T i = 1 y i ] = E [ Y i ( 1 ) ] \mathbb{E}[\frac{1}{n_{1}}\sum_{i: T_{i} = 1}{y_{i}}] = \mathbb{E}[Y_{i}(1)]E[n11i:Ti=1∑yi]=E[Yi(1)]

同样地，

E [ 1 n 0 ∑ i : T i = 0 y i ] = E [ Y i ( 0 ) ] \mathbb{E}[\frac{1}{n_{0}}\sum_{i: T_{i} = 0}{y_{i}}] = \mathbb{E}[Y_{i}(0)]E[n01i:Ti=0∑yi]=E[Yi(0)]

这使我们得出结论，A T E ^ \hat{ATE}ATE^是A T E ATEATE的无偏估计：

E [ A T E ^ ] = A T E \mathbb{E}[\hat{ATE}] = ATEE[ATE^]=ATE

19.5.6 随机实验

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

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

实验：

观测研究：

19.5.7 协变量调整

对于混杂因素该怎么办？

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

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

Y i ( t ) = θ 0 + θ 1 x 1 + . . . + θ p x p + τ t + ϵ Y_{i}(t) = \theta_{0} + \theta_{1}x_{1} + ... + \theta_{p}x_{p} + \tau{t} + \epsilonYi(t)=θ0+θ1x1+...+θpxp+τt+ϵ

**问题：**在这个模型中A T E ATEATE是多少？τ \tauτ

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

• 重要的协变量缺失或者对x xx的真实依赖是非线性的
  
• 有时被贬低地称为**“因果推断”**
  

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

对混杂因素怎么办？

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

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

Y i ( t ) = f θ ( x , t ) + ϵ Y_{i}(t) = f_{\theta}(x, t) + \epsilonYi(t)=fθ(x,t)+ϵ

然后：

A T E = 1 n ∑ i = 1 n f θ ( x i , 1 ) − f θ ( x i , 0 ) ATE = \frac{1}{n}\sum_{i=1}^{n}{f_{\theta}(x_i, 1) - f_{\theta}(x_i, 0)}ATE=n1i=1∑nfθ(xi,1)−fθ(xi,0)

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

• 如果x xx是高维的/其函数形式高度非线性，则非常困难
  
• 需要额外的假设：重叠
  

19.5.8 其他方法

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

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

• 将处理 T 建模为协变量 x 的函数
  
• 将结果 y 建模为 x，T 的函数
  

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

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

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

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

点击显示

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

E [ ( Y ( x ) − Y ^ ( x ) ) 2 ] = E [ ϵ 2 ] + ( g ( x ) − E [ Y ^ ( x ) ] ) 2 + E [ ( E [ Y ^ ( x ) ] − Y ^ ( x ) ) 2 ] . \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*}E[(Y(x)−Y^(x))2]=E[ϵ2]+(g(x)−E[Y^(x)])2+E[(E[Y^(x)]−Y^(x))2].

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

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

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

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

E [ V W ]   =   ∑ v ∑ w v w P ( V = v  and  W = w ) =   ∑ v ∑ w v w P ( V = v ) P ( W = w )     根据独立性 =   ∑ v v P ( V = v ) ∑ w w P ( W = w ) =   E [ V ] E [ W ] \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*}E[VW]=v∑w∑vwP(V=v and W=w)=v∑w∑vwP(V=v)P(W=w)    根据独立性=v∑vP(V=v)w∑wP(W=w)=E[V]E[W]

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

19.6.1 目标

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

19.6.2 步骤 1

模型风险  =   E [ ( Y − Y ^ ( x ) ) 2 ] =   E [ ( g ( x ) + ϵ − Y ^ ( x ) ) 2 ] =   E [ ( ϵ + ( g ( x ) − Y ^ ( x ) ) ) 2 ] =   E [ ϵ 2 ] + 2 E [ ϵ ( g ( x ) − Y ^ ( x ) ) ] + E [ ( g ( x ) − Y ^ ( x ) ) 2 ] \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*}模型风险=E[(Y−Y^(x))2]=E[(g(x)+ϵ−Y^(x))2]=E[(ϵ+(g(x)−Y^(x)))2]=E[ϵ2]+2E[ϵ(g(x)−Y^(x))]+E[(g(x)−Y^(x))2]

右边：

• 第一项是观测方差σ 2 \sigma^2σ2。
  
• 交叉乘积项为 0，因为ϵ \epsilonϵ与g ( x ) − Y ^ ( x ) g(x) - \hat{Y}(x)g(x)−Y^(x)独立，且E ( ϵ ) = 0 E(\epsilon) = 0E(ϵ)=0
  
• 最后一项是我们预测值与x xx处真实函数值之间的均方差差
  

19.6.3 步骤 2

到这个阶段我们有

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

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

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

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

g ( x ) − Y ^ ( x )   =   ( g ( x ) − E [ Y ^ ( x ) ] ) + ( E [ Y ^ ( x ) ] − Y ^ ( x ) ) 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)g(x)−Y^(x)=(g(x)−E[Y^(x)])+(E[Y^(x)]−Y^(x))

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

g ( x ) − Y ^ ( x )   =  模型偏差 − D Y ^ ( x ) g(x) - \hat{Y}(x) ~ = ~ \text{模型偏差} - D_{\hat{Y}(x)}g(x)−Y^(x)=模型偏差−DY^(x)

19.6.4 步骤 3

记住，在x xx处的模型偏差是一个常数，不是一个随机变量。把它看作你最喜欢的数字，比如 10。那么E [ ( g ( x ) − Y ^ ( x ) ) 2 ]   =   模型偏差 2 − 2 ( 模型偏差 ) E [ D Y ^ ( x ) ] + E [ D Y ^ ( x ) 2 ] =   模型偏差 2 − 0 + 模型方差 =   模型偏差 2 + 模型方差 \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*}E[(g(x)−Y^(x))2]=模型偏差2−2(模型偏差)E[DY^(x)]+E[DY^(x)2]=模型偏差2−0+模型方差=模型偏差2+模型方差

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

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

在第 2 步中我们有

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

第 3 步显示

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

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

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

也就是说，

E [ ( Y ( x ) − Y ^ ( x ) ) 2 ] = σ 2 + ( E [ Y ^ ( x ) ] − g ( x ) ) 2 + E [ ( Y ^ ( x ) − E [ Y ^ ( x ) ] ) 2 ] 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]E[(Y(x)−Y^(x))2]=σ2+(E[Y^(x)]−g(x))2+E[(Y^(x)−E[Y^(x)])2]

二十、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.

上面的摘要显示了有关数据库的信息。数据库包含四个表，名称分别为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.

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

20.3 从表中选择

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

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

%%sql
SELECT *
FROM Dragon

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

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

%%sql
SELECT cute, year
FROM Dragon

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

就像这样，我们已经编写了两个 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.

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

%%sql
SELECT DISTINCT year
FROM Dragon

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

每个 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.

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

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

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

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

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

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

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

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

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

20.5 排序和限制输出

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

%%sql
SELECT *
FROM Dragon
ORDER BY cute

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

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

%%sql
SELECT *
FROM Dragon
ORDER BY cute DESC

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

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

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

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

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

%%sql
SELECT *
FROM Dragon
LIMIT 2

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

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

%%sql
SELECT *
FROM Dragon
LIMIT 2
OFFSET 1

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

让我们总结一下到目前为止我们学到的东西。我们知道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.

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

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

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

这里发生了什么？语句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.

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

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

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

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

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

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

有了这个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等）。