R-统计与机器学习研讨会-四-

R 统计与机器学习研讨会（四）

原文：annas-archive.org/md5/e61f21209b8b5ad343780cefdc55a102

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：R 中的微积分

微积分是数学的一个分支，研究两个量之间的关系，这两个量通过一个函数连接，从微观或宏观的角度来看。从微观的角度来看，这种关系（通常表示为 y = f(x)）表现为输出 y 的一个非常小的变化，给定输入 x 的无穷小变化。当切换到宏观视角时，关系变为 y 随 x 变化的累积变化。微观视角对应于微分学，宏观视角对应于积分学，这两者都在本章中介绍。

在 R 中使用积分

本章的所有代码和数据均可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_9/working.R找到。

• 介绍微积分

• 在本章中，我们将涵盖以下主题：

• 在 R 中使用微积分

技术要求

与使用加法或减法等运算的基本数学相比，微积分应用函数和积分来研究变化率。在这里，变化率可以被视为速度，衡量 x 变化时 f(x) 的变化速度。这些变化也有一个方向，意味着当 x 增加时，f(x) 是增加还是减少。

要运行本章中的代码，您需要拥有最新版本的 mosaicCalc 包，写作时为 0.6.0。

介绍微积分

微积分是数学的一个分支，研究变化率，例如曲线在任何点的斜率。它是广泛应用于许多领域的一个基本学科，包括物理、经济学、金融、优化、人工智能（AI）等。微积分最初由两位先生在 17 世纪末开发：戈特弗里德·莱布尼茨和艾萨克·牛顿。牛顿首先发展微积分来分析物理系统，而莱布尼茨独立开发了我们今天使用的符号。

同样，我们可以测量当 x 变化时，f(x) 的变化率如何变化。这样的度量称为加速度，它评估当 x 增加时，f(x) 是否以增加或减少的速率增加或减少。

到本章结束时，您将掌握微分和积分微积分的基本概念。还将介绍 R 的自动微分和积分功能在实际应用中的实现。

微积分的两个主要分支是微分微积分和积分微积分。微分微积分研究特定量的变化率。它检查斜率和曲线的变化率，并研究结果变量 y 对输入变量 x 的微小变化的敏感性。另一方面，积分微积分研究体积或曲线下的****面积（AUC）。

如图 9**.1所示，我们绘制了一个样本函数 y = f(x)并在 x = x 0 处选择了一个任意点。然后我们向 x 0 添加了一个非常小的变化∆ x。这个变化如此之小，以至于在极限情况下会趋向于 0——也就是说，∆ x → 0。当 x 改变时，因变量 y 也会改变，产生一个结果变化∆ y。将这两个项相除给出了 y 在点 x 0 处相对于 x 的敏感性的指示。当∆ x → 0 时，我们称这种敏感性为当 x = x 0 时 y 相对于 x 的导数。用数学表达式表示，我们有 f ′ (x 0) = ∆ y / ∆ x，当∆ x → 0，或者简单地，f ′ (x 0) = lim(∆x→0) ∆ y / ∆ x：

图 9.1 – 展示研究任意点变化率的微分微积分的图形

在点 x 0 处的变化率，或函数 y = f(x)在点 x = x 0 处的敏感性，可以表示为 f ′ (x 0) = ∆ y / ∆ x，当∆ x 无限小的时候——也就是说，∆ x → 0。另一方面，积分微积分表明特定范围 x ∈ [a, b]的 AUC f(x)，如图图 9**.2所示：

图 9.2 – 展示积分微积分作为总 AUC 的图形

我们将积分的结果表示为 S = ∫ a b f(x)dx，在这种情况下报告一个正量。这意味着积分的结果也可以是负的，如果曲线位于x轴下方。

更多关于函数的内容

注意，函数是两个元素集合之间的映射机器。对于输入集中的每个元素，输出集中只有一个对应的元素。输入集中所有可能元素的集合称为定义域，所有对应元素的集合称为值域。图 9*.3*显示了四种不同的映射场景。前三个是根据函数的定义有效的映射，而最后一个由于一对一多映射关系而失败：

图 9.3 – 四种不同的映射场景

前三个映射是有效的函数，而最后一个由于最后一个元素中的一对多映射（3->6 和 3->2）而无效。

对于函数 y = f(x) = x 2，我们经常使用多个名称，如下所示：

• 变量 x 可以称为自变量、特征、协变量或输入。

• 变量 y 可以称为因变量、结果、目标、响应或输出。

• 映射函数 f(x) = x² 可以称为函数、映射、投影或假设。

图 9*.4*总结了这些术语，并提供了一组输入-输出样本及其相应的图形。注意，输入可以是任何数字或表示为一般变量：

图 9.4 – 示例映射函数 y = f(x) = x²

垂直线测试

评估一个映射或曲线是否为函数的一种技术是垂直线测试。特别是，平面上的任何垂直线最多只能与一个函数的图形相交一次。也就是说，对于任何函数 f : A → B，每个元素 x ∈ A 最多映射到 B 中的一个值 f(x) ∈ B。图 9*.5*说明了两个曲线，其中第一个曲线只与垂直线相交一次，因此是一个有效的函数，而第二个曲线与垂直线相交两次，因此不满足函数的定义：

图 9.5 – 使用垂直线测试评估曲线是否为函数

函数对称性

函数可能具有许多特殊性质，函数对称性是其中之一。例如，如果 f(x) = f(-x)，则函数被称为偶函数，例如 f(x) = |x|，f(x) = x²，f(x) = cos(x)。如果 f(-x) = -f(x)，则函数被称为奇函数，例如 f(x) = x³。

增减函数

函数也可以是增函数或减函数。如果对于任何 x_1 > x_2，f(x_1) > f(x_2)，则函数在区间[a, b]上是增函数；如果对于任何 x_1 > x_2，f(x_1) < f(x_2)，则函数在区间[a, b]上是减函数。

函数的斜率

一条直线的斜率 m 可以根据任意两点(x_1, y_1)和(x_2, y_2)使用以下定义来计算：

m = (y_2 - y_1) / (x_2 - x_1)

与图 9**.1中的导数定义相联系，导数本质上是在点 x = x⁰ 处的切线斜率。

函数复合

函数复合指的是一个函数由多个嵌套函数组成的情况，表示为 f(g(x)) = (f ∘ g)(x)。映射序列如下：输入变量 x 首先通过函数 g(x)进行转换，然后通过函数 f(x)进行转换。在 f(g(x))中，g(x)是内部函数，f(x)是外部函数。

例如，给定 f(x) = 2x - 1，g(x) = x³，我们可以得到 f(g(x))如下：

f(g(x)) = 2x³ - 1

注意，函数复合不是可交换的——也就是说，f(g(x)) ≠ g(f(x))。

常见函数

让我们来看看几种常见的函数类型。

幂函数

幂函数，f(x) = k x^p，其中 k 和 p 是常数，x 是变量。k 也称为系数。以下是一些例子：

• 立方函数：f(x) = x³

• 平方根函数：f(x) = √x

• 立方根函数：f(x) = 3√x

• 线性函数：f(x) = x

• 绝对值函数：f(x) = |x|

• 平方函数：f(x) = x²

多项式函数

多项式函数，f(x) = a_n x^n + a_{n−1} x^{n−1} + … + a_1 x + a_0，其中 a_n 是首项系数。n 是一个非负整数，称为多项式的次数，系数 a_0, … , a_n 是实数，且 a_n ≠ 0。

有理函数

有理函数，f(x) = n(x)/(d(x))，其中 n(x) 和 d(x) 都是多项式，且 d(x) ≠ 0。

指数函数

指数函数，f(x) = b^x，其中 b > 0 且 b ≠ 1。

对数函数

对数函数：f(x) = log_bx，其中 b > 0 且 b ≠ 1。注意，当 y = log_bx 时，我们有等价形式 x = b^y。

在快速了解了许多函数及其性质之后，让我们看看极限的概念，它与导数相关。

理解极限

函数 f(x) 的极限表示，当输入 x 接近但不等于一个数 c 时，f(x = c) 的值将接近一个实数 L。数学上，我们可以表示如下：

lim x→c f(x) = L

或者，它也可以表示为 f(x) → L 当 x → c。注意，极限 lim x→c f(x) 的值可能不一定等于函数在 x = c 处的值——即 f(c)。这种等式，lim x→c f(x) = f(c)，仅在 f(x) 是连续函数时发生。

无穷大极限

当极限 lim x→c f(x) 在 x → c 时不存在，我们说函数 f(x) 趋向于无穷大，因此当 x → c 时导致无穷大极限。我们也可以说 x = c 是函数 f(x) 的垂直渐近线。

一个例子是 f(x) = 1/(x - 1)。我们知道定义域是 x ≠ 1。当 x → 1 时，我们有 lim x→1 1/(x - 1) = ∞，这分解为从左侧接近 lim x→1− 1/(x - 1) = −∞ 和从右侧接近 lim x→1+ 1/(x - 1) = ∞。这两个结果共同导致 lim x→1 1/(x - 1) = ∞，因此 lim x→1 1/(x - 1) 不存在。垂直渐近线是 x = 1。

无穷大极限

无穷大极限关注的是当 x → ∞ 或 x → −∞ 时 f(x) 的值。这给出了函数的水平渐近线——即 lim x→∞ f(x) 和 lim x→−∞ f(x)。我们称 y = b 为水平渐近线，如果 lim x→∞ f(x) = b 或 lim x→−∞ f(x) = b。

让我们看看几个例子。对于平方函数 f(x) = x²，我们有 lim x→∞ x² = ∞ 和 lim x→−∞ x² = ∞。对于平方根函数 f(x) = √x，由于定义域是 x ≥ 0，我们有 lim x→∞ √x = ∞。对于自然对数函数 f(x) = lnx，由于定义域是 x > 0，我们有 lim x→∞ √x = ∞。对于一般结果，我们有 lim x→∞ x^p = ∞ 和 lim x→∞ 1/x^p = ∞。

下一个部分正式介绍了导数。

介绍导数

对于函数 y = f(x)，f 在 x 处的导数，记作 f′(x)，定义为以下内容：

f′(x) = lim(h→0) f(x + h) - f(x) / h，如果极限存在。

当 f′(x)在区间[a, b]内的每个 x 值上存在时，我们说函数 f 在 x ∈ [a, b]上是可导的。

导数有多种应用。除了将导数解释为切线斜率之外，最常见的应用是我们可以使用导数来确定函数在特定点上是上升/增加还是下降/减少。它还表示在给定点 x 的瞬时变化率，或速度。

一个图形说明将有助于我们的理解。图 9.6提供了一个 f(x)的样本曲线，其中有一条割线连接两个点(a, f(a))和(b, f(b))。在这里，割线是穿过函数曲线的两个或更多不同点的直线。我们通过在点(a, f(a))上加上一个小的量 h 得到第二个点(a + h, f(a + h))。连接这两个点给出通过这两个点穿过函数 f(x)的割线。当 h 趋近于 0 时，割线将逐渐接近在 x = a 处得到的切线。最终，割线将在极限中与切线重叠，割线的斜率，表示为 f(a + h) - f(a) / h，将变为它们在极限中重叠时的切线斜率：

图 9.6 – 展示导数推导过程的图形

在图 9.6中，第二个点(a + h, f(a + h))是通过在点(a, f(a))上加上一个小的量 h 得到的。连接这两个点给出割线。当 h 趋近于 0 时，割线将逐渐接近 x = a 处的切线，最终在极限中与切线重叠。

割线的斜率，也称为平均变化率，计算如下：

f(a + h) - f(a) / (a + h) - a = f(a + h) - f(a) / h

这是在 h ≠ 0 的前提下提供的。当 h 趋近于 0 时，割线无限接近于 x = a 处的切线，给出该点的瞬时变化率（如果极限存在的话）：

当 h 趋近于 0 时，函数 f 在 a + h 处的极限为 f(a + h)减去 f(a)除以 h

通常，如果 y = f(x)，我们可以将导数表示为 f′(x)，y′，dy/dx，或者 d/dx f(x)。

以下部分介绍一些常见的导数。

常见导数

在这里，我们提供了一些常见函数的导数列表：

• 常数函数：d/dx(c) = 0

• 幂函数：d/dx(x^n) = nx^(n-1)

• 指数函数：d/dx(e^x) = e^x, d/dx(b^x) = (lnb)b^x, 其中 b > 0

• 对数函数：d/dx(lnx) = 1/x, d/dx(log_bx) = 1/(lnb)x, 其中 b > 0, b ≠ 1

• 正弦和余弦函数：d/dx(sinx) = cosx, d/dx(cosx) = -sinx

计算导数通常涉及多个基本函数。下一节将介绍关于导数计算的一组性质和规则。

导数的常见性质和规则

我们引入两个常见的性质——常数倍性质和和差性质——然后是三个常见的规则：乘法法则、除法法则和链式法则。我们假设以下列表中的所有极限都存在：

• 常数倍性质：对于常数 c，如果 y = f(x) = ch(x)，那么 f′(x) = ch′(x)。用不同的方式表达，我们有 y′ = ch′，并且 dy/dx = c dh/dx。

• 和差性质：如果 y = f(x) = h(x) ± g(x)，那么 f′(x) = h′(x) ± g′(x)。用不同的方式表达，我们有 y′ = h′ + g′，并且 dy/dx = dh/dx + dg/dx。

• 乘法法则：如果 y = f(x) = h(x)g(x)，那么 f′(x) = h′(x)g(x) + h(x)g′(x)。用不同的方式表达，我们有 y′ = h′g + hg′，并且 dy/dx = g dh/dx + h dg/dx。

• 除法法则：如果 y = f(x) = h(x) / g(x)，那么 f′(x) = g(x)h′(x) - h(x)g′(x) / [g(x)]²。用不同的方式表达，我们有 y′ = gh′ - hg′ / g²，并且 dy/dx = g dh/dx - h dg/dx / g²。

• 链式法则：给定一个复合函数 y = f(x) = h(g(x))，我们有 f′(x) = h′(g(x))g′(x)。等价地，如果我们有 y = h(u)和 u = g(x)，那么 dy/dx = dy/du du/dx。

在下一节中，我们将转向积分微积分。

介绍积分微积分

积分是微分的逆运算。例如，对距离函数 S(t)关于时间的微分给出特定时间点的速度 v(t)。我们可以通过将时间周期分为 n 个等间距的区间来计算从时间 t = a 到 t = b 的累积距离 S = ∑i=0^(n-1) v(ti)Δt，在 n 个区间数趋于无穷大时——即 n → ∞，我们得到 S = lim(n→∞) ∑i=0^(n-1) v(ti)Δt = ∫a^b v(t)dt。

积分的结果称为不定积分，它允许我们从导数重建函数。这是微分运算的逆运算。形式上，如果函数 F(x)是 f(x)的不定积分，那么 F′(x) = f(x)。

图 9.7 说明了积分与 AUC 之间的等价性。我们可以使用积分来找到一般函数 f(x)的 AUC。这个面积可以通过添加趋近于零宽度的矩形切片来获得。每个矩形的宽度是 dx，根据定义，它指的是 x 的无穷小变化，即 x 的微分：

图 9.7 – 说明积分与 AUC 等价性的图表

注意，积分可以是定积分或不定积分。定积分有一个明确的积分起始点和结束点。例如，如果 x ∈ [a, b]，定积分变为∫a^b f(x)dx。另一方面，不定积分没有明确的边界，给出∫f(x)dx，它用来表示 f(x)的所有不定积分的集合。

下一节将更深入地探讨不定积分。

不定积分

通常，如果我们有 F′(x) = f(x)，那么我们可以将不定积分表示为 ∫ f(x)dx = F(x) + C，其中 C 是一个常数，在这里用来表示所有导数为 f(x) 的反导函数的集合。换句话说，我们有 d/dx[∫ f(x)dx] = f(x) 和 ∫ F′(x)dx = F(x) + C。

现在，让我们更深入地理解不定积分 ∫ f(x)dx 的表达式。∫ 是积分符号，f(x) 是被积函数，dx 是沿 x 的一个切片，表示反导数是在 x 上进行的。图 9.8 展示了命名约定：

图 9.8 – 总结不定积分的命名约定

我们现在将研究常见基本函数的不定积分。

基本函数的不定积分

以下列表提供了基本函数的不定积分。再次注意，常数 C 用来表示所有具有相同导数的反导函数的集合：

• ∫ x^n dx = (1/(n+1))x^(n+1) + C

• ∫ e xdx = e x + C

• ∫_b^xdx = (1/2)x² + C

• ∫ 1/x dx = ln|x| + C

• ∫ sinxdx = − cosx + C

• ∫ cosxdx = sinx + C

不定积分的性质

本节涵盖了与导数函数相对应的一组性质。以下两个性质是直接的线性运算：

• ∫ Kf(x)dx = K∫ f(x)dx

• ∫ [f(x) ± g(x)]dx = ∫ f(x)dx ± ∫ g(x)dx

以下列表包含了一些计算不定积分的一般公式：

∫ [f(x)]^n f'(x)dx = [(f(x))]^(n+1)/(n+1) + C, n ≠ 1

∫ e^(f(x)) f'(x)dx = e^(f(x)) + C

∫ 1/f(x)f'(x)dx = ln|f(x)| + C

接下来，我们介绍最广泛使用的技术之一：分部积分。

分部积分

根据之前引入的乘积法则，我们知道以下内容：

d/dx[f(x)g(x)] = f′(x)g(x) + f(x)g′(x)

我们可以对方程的两边进行积分，并重新排列项，得到以下结果：

∫ f(x)g'(x)dx = f(x)g(x) − ∫ g(x)f'(x)dx

这给出了分部积分公式。更简洁地说，我们可以定义 v = g(x) 和 u = f(x)。因此，dv = g′(x)dx 和 du = f′(x)dx。前面的方程变为：

∫ udv = uv − ∫ vdu

接下来，我们将更深入地探讨定积分。

定积分

如果 f(x) 是 x ∈ [a, b] 上的连续函数，且 F(x) 是 f(x) 的反导数，即 F(x) = ∫ f(x)dx + C。以下公式给出了微积分的基本定理：

∫ a bf(x)dx = F(b) − F(a)

这个表达式说明，在区间 [a, b] 内的 AUC f(x) 可以通过计算反导函数 F(x) 在两个端点的评估差来得到。

定积分与不定积分具有大部分相同的性质，以下是一个额外的性质：

∫_a^cf(x)dx = ∫_a^bf(x)dx + ∫_b^cf(x)dx

下一个部分将探讨使用 R 进行微积分的实现方面。

在 R 中使用微积分

在本节中，我们将使用mosaicCalc包来执行与微积分相关的操作。我们将使用的函数大多是分析性的（有显式表达式）且性质简单。以下代码片段检查此包是否已安装，如果条件评估为true，则会安装该包，然后将其导入当前会话：

if(!require("mosaicCalc")){
  install.packages("mosaicCalc")
}
library(mosaicCalc)

为了避免过多的想象，了解未知函数的一个好方法是绘制它们。让我们看看如何使用mosaicCalc包绘制函数。

绘制基本函数

在使用mosaicCalc包绘制函数时，我们需要指定一些输入参数。总的来说，我们需要指定函数的表达式、用于绘制函数的输入变量（s）、每个输入变量的域以及预先设置的参数值。

在绘图方面，有三个常见的绘图函数：slice_plot()，用于绘制只有一个输入变量的函数，contour_plot()，用于绘制具有两个输入变量的函数，以及interactive_plot()，用于绘制交互式图形。这些函数允许我们将数学表达式转换为代码，然后绘制图形。让我们通过以下练习来探索每个函数。

练习 9.1 – 绘制基本函数

本练习将探索前面函数的绘图功能，从slice_plot()开始。请按照以下步骤进行：

1. 使用slice_plot()函数绘制函数 y = 2x + 1，其中 x ∈ [−5, 5]，如下所示：

   
   >>> slice_plot(2*x+1 ~ x, domain(x = range(-5, 5)))
   

   在这里，我们将函数表达式 2x + 1 放在波浪号（~）的左边，输入变量 x 放在右边。我们还通过传递边界给range()函数来指定域。运行此命令将生成如图 9**.9所示的输出：

图 9.9 – 使用 slice_plot()函数绘制直线

我们还可以通过设置直线 y = mx + b 的系数 m 和截距 b 的先验参数来使函数更通用。这使得代码更通用，因为我们只需要在以下代码片段中指定m和b来绘制任何直线：

m = 2
b = 1
slice_plot(m*x+b ~ x, domain(x = range(-5, 5)))

注意，如果没有设置输入参数的初始值，将会发生错误，如下面的代码片段所示：

>>> slice_plot(a*x+b ~ x, domain(x = range(-5, 5)))
Error in slice_plot(a * x + b ~ x, domain(x = range(-5, 5))) :
  Parameter <a> without specified numerical values.

此外，我们可以使用makeFun()函数给要绘制的函数命名，如下面的代码片段所示：

f = makeFun(2*x+1 ~ x)
slice_plot(f(x) ~ x, domain(x = range(-5, 5)))

一旦给函数命名，我们就可以传递任意输入值来评估函数。例如，设置 x = 2 将返回5，如下所示：

>>> f(x=2)
5

接下来，我们来看如何为具有两个输入变量的函数生成等高线图。

1. 使用contour_plot()函数绘制方程 z = 2x + 3y 的等高线图，如下所示：

   
   >>> contour_plot(2*x + 3*y ~ x & y, domain(x=-5:5, y=-5:5))
   

   在这里，我们使用 & 符号来表示多个输入变量，并为两个变量设置相应的范围。运行此命令将生成如图 9*.10* 所示的输出：

图 9.10 – 使用 contour_plot() 函数生成 z = 2x + 3y 的等高线图

接下来，我们看看如何使用 interactive_plot() 函数生成交互式绘图。

1. 使用 interactive_plot() 函数为相同的表达式生成交互式 3D 图，如下所示：

   
   >>> interactive_plot(2*x + 3*y ~ x & y, domain(x=-5:5, y=-5:5))
   

   运行此命令将生成如图 9*.11* 所示的输出。请注意，生成的图形是一个交互式 HTML 小部件，允许我们移动并当鼠标悬停时显示辅助信息：

图 9.11 – 使用 interactive_plot() 函数生成 z = 2x + 3y 的交互式 3D 图

下一节将介绍通过使用 D() 函数进行微分操作来处理导数。

处理导数

在 R 中，通过 D() 微分操作符完成微分操作，它输入单个表达式并输出导数函数。输入指定与之前绘图函数所需的相同表达式。例如，要指定函数 y = x² + 1，我们可以将 x²+1 ~ x 传递给 D() 函数，然后它将自动计算导数函数 y′ = 2x。然后我们可以将结果分配给另一个变量，该变量作为导数函数，并可用于在定义域中的任何点评估导数值。

以下代码片段说明了获取导数函数的过程。D() 操作符的返回值是导数函数，打印出来后为 y′ = 2x：

f_prime = D(x²+1 ~ x)
>>> f_prime
function (x)
2 * x
<bytecode: 0x144378300>

以下代码片段评估了两个输入值。结果显示，D() 函数能够正确计算导数函数并在任意输入位置进行评估：

>>> f_prime(1)
a
>>> f_prime(2)
4

注意，D() 函数可以在计算导数函数时执行之前提到的规则。例如，要获得 y = sin( x² − 5) 的导数，我们会调用链式法则并计算导数为 y' = 2x cos( x² − 5)。D() 函数为我们完成了这项工作，如下面的代码片段所示：

f_prime = D(sin(x²-5) ~ x)
>>> f_prime
function (x)
2 * x * cos(x² - 5)

让我们也验证商法则。在以下代码片段中，我们将函数 y = 2x² + x + 1 传递给 D()，理想情况下它应该返回导数函数 y' = 2(x + 1) − 2x(x + 1)²。结果显示这确实是正确的：

f_prime = D(2*x/(x+1) ~ x)
>>> f_prime
function (x)
{
    .e1 <- 1 + x
    (2 - 2 * (x/.e1))/.e1
}

下一节将介绍在函数中使用符号参数。

使用符号参数

符号参数在构建函数时提供了通用性。与之前一样，我们可以在将参数的先前值传递给函数表达式之前对其进行编码。以下练习说明了这一点。

练习 9.2 – 使用符号参数

在这个练习中，我们将研究一个一般函数 y = A x³ + Bx + 3，其中 A 和 B 是常数，x 是唯一的随机输入变量。导数函数将是 y′ = 3A x² + B。按照以下步骤进行：

1. 按如下方式计算函数 y = A x³ + Bx + 3 的导数：

   
   f_prime = D(A*x³+B*x+3 ~ x)
   >>> f_prime
   function (x, A, B)
   3 * A * x² + B
   

   我们看到导数函数计算正确。

2. 当 x = 2，A = 2，B = 3 时，按如下方式评估导数函数：

   
   >>> f_prime(x=2, A=2, B=3)
   27
   

   我们也可以通过多个点来评估函数，这相当于在特定范围内绘制导数函数。

3. 在 x ∈ [−5, 5] 的范围内，以 A = 2 和 B = 3 为条件绘制导数函数，如下所示：

   
   >>> slice_plot(f_prime(x, A=2, B=3) ~ x, domain(x=range(-5,5)))
   

   执行此命令将生成 图 9.12 中所示的输出：

图 9.12 – 可视化导数函数 y′ = 6x² + 3

下一个部分将介绍二阶导数。

处理二阶导数

二阶导数，表示为 f″(x)，就是原始函数 f(x) 的一阶导数的导数。在先前的例子中，当 y = A x³ + Bx + 3 时，我们有 y′ = 2A x² + B。再求一次导数得到 y″ = 4Ax。

让我们通过以下练习看看如何获得二阶导数。

练习 9.3 – 计算二阶导数

在这个练习中，我们将计算原始函数 f(x) 的二阶导数。二阶导数可以看作是对 f(x) 关于 x 进行两次求导，得到 f″(x) = y″ = d²/dx² f(x)。双重求导是通过在波浪号右侧有两个 x 来实现的。按照以下步骤进行：

1. 按如下方式计算二阶导数 f″(x)：

   
   f_pprime = D(A*x³+B*x+3 ~ x & x)
   >>> f_pprime
   function (x, A, B)
   6 * A * x
   

   结果显示，绘制函数的二阶导数计算正确。

2. 当 x = 2，A = 2，B = 3 时，按如下方式评估二阶导数函数：

   
   >>> f_pprime(x=2, A=2, B=3)
   24
   

   实际上，由于二阶导数函数与参数 B 无关，B 的任何值都会得到相同的结果。例如，以下代码即使在 B = 1 时也会返回相同的结果：

   >>> f_pprime(x=2, A=2, B=1)
   24
   

3. 在 x ∈ [−5, 5] 的范围内，以 A = 2 为条件绘制二阶导数函数，如下所示：

   
   >>> slice_plot(f_pprime(x, A=2) ~ x, domain(x=range(-5,5)))
   

   执行此命令将生成 图 9.13 中所示的输出。在这里，我们注意到，与一阶导数函数的钟形曲线相比，二阶导数函数是一条直线。这条直线表示一阶导数函数变化的速率：

图 9.13 – 可视化二阶导数函数 y″ = 12x

下一个部分将介绍更一般的偏导数函数。

处理偏导数

二次导数函数可以看作是偏导数函数的一种特殊情况。假设一个二维函数 z = f(x, y)。通过两次对 x 求导得到 x 的二次导数。我们也可以先对 x 求导，然后对 y 求导，得到 d²/dxdy z 或等价地，d²/dydx z。让我们通过以下练习看看这是如何工作的。

练习 9.4 – 计算偏导数

在这个练习中，我们将基于原始函数 z = A x² + Bxy + C y² 计算三个不同的偏导数：d²/dx² z，d²/dxdy z（或 d²/dydx z），以及 d²/dy² z。按照以下步骤进行：

1. 计算 z 关于 x 的二次导数。检查结果是否为 d²/dx² z = 2A：

   
   f_pprime = D(A*x² + B*x*y + C*y² ~ x & x)
   >>> f_pprime
   function (x, y, A, B, C)
   2 * A
   

   结果显示计算是正确的。

2. 计算 z 关于 x 的偏导数，然后是 y。检查结果是否为 d²/dxdy z = B：

   
   f_pprime = D(A*x² + B*x*y + C*y² ~ x & y)
   >>> f_pprime
   function (x, y, A, B, C)
   B
   

   结果显示计算是正确的。我们也可以先对 y 求导，然后对 x 求导。如下代码片段所示，这给出了相同的结果：

   f_pprime = D(A*x² + B*x*y + C*y² ~ y & x)
   >>> f_pprime
   B
   

3. 计算 z 关于 x 的偏导数，然后是 y。检查结果是否为 d²/dy² z = 2C：

   
   f_pprime = D(A*x² + B*x*y + C*y² ~ y & y)
   >>> f_pprime
   function (x, y, A, B, C)
   2 * C
   

   结果显示计算是正确的。

下一节将介绍如何使用 R 计算积分或不定积分。

在 R 中处理积分

回想一下，微分是通过 D() 函数进行的。假设 y = A x² + Bx + 3，我们有 y′ = f′(x) = 2Ax + B。我们可以绘制原始函数 f(x) 和其导数函数 f′(x)，以方便比较。

在以下代码片段中，我们使用 makeFun() 函数创建这个表达式，并将其命名为 f：

f = makeFun( A*x² + B*x + 3 ~ x)
>>> f
function (x, A, B)
A * x² + B * x + 3

我们可以进行如下简单评估：

>>> f(1, A=1, B=1)
5

现在，我们得到了导数函数并将其存储在 f_prime 中：

f_prime = D(f(x) ~ x)
>>> f_prime
function (x, A, B)
2 * A * x + B

这是一个从原始函数导出的新函数。我们对此函数也进行了简单的评估，如下所示：

>>> f_prime(x=1, A=1, B=1)
3

现在，让我们绘制原始函数 f(x) = x² + x + 3：

>>> slice_plot(f(x) ~ x, domain(x = -1:1)) %>%
  gf_labs(title = "Original function f(x)")

运行命令生成了如图 9*.14* 所示的输出。在这里，我们使用了 gf_labs() 函数来设置图表的标题：

图 9.14 – 可视化原始函数 f(x) = x² + x + 3

我们还可以如下绘制导数函数 f′(x) = 2x + 1：

>>> slice_plot(f_prime(x, A=1, B=1) ~ x, domain(x =-1:1), color = "red") %>%
  gf_labs(title = "Derivative function f'(x)")

运行命令生成了如图 9*.15* 所示的输出：

图 9.15 – 可视化导数函数 f′(x) = 2x + 1

在有了导数函数之后，我们可以使用 antiD() 函数来获取其不定积分，如下所示：

F_integral <- antiD(f_prime(x) ~ x)
>>> F_integral
function (x, A, B, C = 0)
A * x² + x * B + C

注意，除了作为不定积分函数的额外输入参数添加的附加常数 C（默认值为 0）外，原始函数的正确形式得到了恢复。

我们也可以评估不定积分函数，如下所示：

>>> F_integral(1, A=1, B=1)
2

将 A=1, B=1, 和 C=0 代入，不定积分函数 F(x) = x² + x 的可视化如下：

>>> slice_plot(F_integral(x, A=1, B=1) ~ x, domain(x=-1:1)) %>%
  gf_labs(title = "Antiderivative function F_integral(x)")

运行命令会生成如图 9.16 所示的输出：

图 9.16 – 可视化导数函数 f′(x) = 2x + 1

现在，如果我们再次对反导数函数 F(x) = A x² + Bx + C 求导，我们预计会得到相同的导数函数 f(x) = 2Ax + B。下面的代码片段验证了这一结果：

f_prime2 = D(F_integral(x) ~ x)
>>> f_prime2
function (x, A, B, C = 0)
2 * A * x + B

让我们暂停一下，看看导数函数和它们的原函数之间的关系。

更多关于反导数的内容

产生 F(x) = ∫ f'(x)dx + C 的反导数操作是生成 f'(x) 的导数的逆过程。它代表一个函数族，包括原始函数 f(x)。这些函数彼此之间非常相关。具体来说，导数函数 f'(x)，作为一个导出函数，告诉了原始函数 f(x) 在任意输入点 x 的变化率。它给出了 f(x) 在微观镜头中的局部性质，测量了 f(x) 在当前点的敏感性。

然而，我们并不总是在每种情况下都需要导出 f'(x)。有时，我们将从导数函数 f'(x) 开始工作，并会对导出原始函数——即反导数函数 F(x) 感兴趣。这给出了原始函数 f(x) 的全局性质，表示在特定范围内的累积值。

获得未知原始函数的过程称为反导数，或积分。结果是积分。根据积分是否在特定边界上计算，我们有一个不定积分和一个定积分。

积分操作生成一个反导数函数族。与前面的例子一样，反导数函数是 F(x) = A x² + Bx + C。将 A=1 和 B=1 设置，对于不同的 C 值，我们将在相同的输入位置获得不同的结果，如下面的代码片段所示：

>>> F_integral(x=1, A=1, B=1, C=0)
2
>>> F_integral(x=1, A=1, B=1, C=1)
3
>>> F_integral(x=1, A=1, B=1, C=2)
4

虽然这些是不同的原始函数，但它们都共享相同的导数函数，如下面的代码片段所示：

>>> D(F_integral(x, A=1, B=1, C=0) ~ x)
function (x)
2 * x + 1
>>> D(F_integral(x, A=1, B=1, C=1) ~ x)
function (x)
2 * x + 1
>>> D(F_integral(x, A=1, B=1, C=2) ~ x)
function (x)
2 * x + 1

以 F(x) = ∫ f′(x)dx + C 表示的反导数函数族，因此对应于无限多个函数，包括原始函数 f(x)。这些无限多个反导数函数实际上是 f(x) 的垂直平移，给出 F(x) = f(x) + C。

此外，请注意，导数函数与原始函数共享相同的输入参数集，如下面的代码片段所验证：

>>> f
function (x, A, B)
A * x² + B * x + 3
<bytecode: 0x128185eb0>
>>> f_prime
function (x, A, B)
2 * A * x + B
<bytecode: 0x11fbc6628>

然而，反导数函数需要一个额外的参数——常数 C，如下所示：

>>> F_integral
function (x, A, B, C = 0)
A * x² + x * B + C
<bytecode: 0x128303540>

下一节将探讨如何计算定积分。

计算定积分

计算定积分需要两次评估不定积分：一次在积分的起点，一次在积分的终点。

例如，假设不定积分函数为 F(x) = x² + x + C，对于 A=1 和 B=1。要计算不定积分 ∫₂³ f′(x)dx，我们需要在 x = 2 和 x = 3 处评估 F(x)，然后取它们的差，得到 F(3) − F(2)。以下代码片段显示了结果：

>>> F_integral(x=3, A=1, B=1) - F_integral(x=2, A=1, B=1)
6

这里需要注意的是，在两次评估中，常数 C 假设为 0。实际上，它假设的值并不重要，因为两次评估中的常数总会相互抵消。换句话说，我们总会得到一个固定的定积分，尽管在特定的输入点，不定积分对应着无限多个值。

摘要

在本章中，我们介绍了微积分的基础，包括微分学和积分学。在第一部分，我们介绍了这两个微积分分支的直观理解，并涵盖了常见函数及其性质的基础。我们开始介绍极限的概念及其与导数定义的联系，随后介绍了常见的导数规则和性质。我们还讨论了积分学，包括不定积分和定积分，以及它们的规则和性质。

第二部分和第三部分简要介绍了 R 语言中的实现。我们介绍了如何使用 D() 和 antiD() 函数进行常见的微分和积分操作，并通过几个示例说明了它们的用法以及导数函数与其反导数之间的转换。

在下一章中，我们将进入数学统计学的领域，从概率的基础开始。

第三部分：R 中的数学统计学基础

随着我们开始本书的第三部分，我们将从第二部分中探讨的数学基础过渡到高级统计方法和模型构建的领域。虽然前面的部分在基本统计和数学基础方面奠定了基础，但本部分旨在将你的能力提升到更高级的水平，同时利用 R 的强大功能进行实际应用。

到本部分结束时，你将掌握一套强大的高级统计和建模技能。本部分提供了你进入数据科学更高级主题所需的先进知识和实践经验。

本部分包含以下章节：

• 第十章, 概率基础

• 第十一章, 统计估计

• 第十二章, R 中的线性回归

• 第十三章, R 中的逻辑回归

• 第十四章, 贝叶斯统计

第十章：概率基础

概率分布是统计学和机器学习中的一个基本概念。它描述了控制实验或随机过程中潜在结果或事件生成的潜在分布。根据特定领域和数据的特点，存在不同类型的概率分布。一个合适的概率分布是理解和建模随机过程和事件行为的有用工具，在开发数据驱动的预测和优化模型时，提供了方便的工具进行决策和预测。

到本章结束时，你将理解常见的概率分布及其参数。你还将能够使用这些概率分布来执行常规任务，如 R 中的抽样和概率计算，以及常见的抽样分布和顺序统计量。

在本章中，我们将涵盖以下主题：

• 介绍概率分布

• 探索常见离散分布

• 发现常见连续分布

• 理解常见抽样分布

• 理解顺序统计量

技术要求

要运行本章中的代码，你需要拥有以下软件包的最新版本：

• ggplot2, 3.4.0

• dplyr, 1.0.10

请注意，前面提到的软件包版本是编写本章时的最新版本。

本章的代码和数据可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_10/working.R找到。

介绍概率分布

概率分布提供了一个理解和预测随机变量行为的框架。一旦我们知道了数据生成的潜在概率分布，我们就可以更明智地决定事物可能出现的可能性，无论是在预测还是优化环境中。换句话说，如果选定的概率分布可以很好地模拟观察到的数据，我们就有了一个强大的工具来预测潜在的未来值，以及这种发生的不确定性。

在这里，一个随机变量是一个值不固定且可能假设多个或无限多个可能值的变量，它代表了随机事件的结果（或实现）。概率分布使我们能够表示和分析这些结果的可能性，为各种场景下潜在的不确定性提供了一个全面的视角。概率分布将随机变量，表示为 x，转换为概率，P(x)，这是一个介于 0 和 1 之间的浮点数。概率分布可以是概率密度函数（PDF）或概率质量函数（PMF），它指定了观察到一个结果对于连续变量（或离散变量）的概率，或者是一个累积分布函数（CDF），它提供了随机变量小于或等于给定固定数量的总概率。

在以下示例中，我们使用 f(x) 来表示 x 的概率密度函数，并使用 F(x) 来表示 CDF。

概率分布主要有两大类：离散概率分布和连续概率分布。离散概率分布处理离散变量，这些是只能假设有限或可数个可能值的随机变量。例如，如果我们有一个概率分布，指定了一周内下雨天的概率，那么潜在的随机变量是星期几，只能取 1 到 7 之间的整数。假设离散随机变量有 C 个可能值。对于给定的可能值，x_i，其中 i ∈ {1, … , C}，相应的 PMF 是 f(x_i) ∈ [0,1]，并且所有概率的总和应为 1，即 ∑_i=1^C f(x_i) = 1。我们可以将 PMF f(x) 视为一个条形图，它指定了每个离散输入的概率输出。

我们将介绍一些常见的离散分布。例如，二项分布模型了在固定数量的具有相同成功概率的伯努利试验中成功的次数，泊松分布模型了在固定时间或空间间隔内的事件数量，几何分布模型了在一系列伯努利试验中第一次成功所需的试验次数。这些将在本章后面进行介绍。

相反，连续概率分布涉及连续变量，这些变量可以在指定的范围内取无限多个值，表示为 𝒳。随机变量 x ∈ 𝒳 现在是连续的，相应的 PDF f(x) ∈ [0,1] 满足 ∫ f(x)dx = 1，其中 x ∈ 𝒳，我们将求和符号改为积分，以考虑连续变量 x 的无限多个可能值。我们可以将 PDF f(x) 视为一个线形图，它指定了每个连续输入的概率输出。

我们将介绍一些广泛使用的连续分布，从正态分布（或高斯分布）开始，它描述了许多自然量的分布。其他连续分布的例子包括指数分布，它模拟泊松过程中独立事件之间的时间，以及均匀分布，它将等概率分配给指定范围内的所有结果。

图 10.1 总结了这两种类型的概率分布：

图 10.1 – 总结两种概率分布的分类。两种分布的总和为 1

注意，每个概率分布都有一个与之相关的闭式表达式和一组相应的参数。我们将在下一部分突出显示每个分布的表达式和参数，从离散分布开始。

探索常见的离散概率分布

离散概率分布以其对应的概率质量函数（PMF）为特征，为输入随机变量的每个可能结果分配一个概率。离散分布中所有可能结果的概率之和等于 1，导致 ∑ i=1 C f(x_i) = 1。这也意味着其中一个结果必须发生，即 f(x_i) > 0，∀ i = 1, … , C。

离散概率分布在各个领域都至关重要，例如金融。它们常用于统计分析，包括假设检验、参数估计和预测建模。我们可以使用离散概率分布来量化不确定性，进行预测，并深入了解观察到的结果背后的数据生成过程。

让我们从最基本的离散分布开始：伯努利分布。

伯努利分布

伯努利分布是一种基本的离散概率分布，它指定了在单个伯努利试验中二元随机变量的行为。在这里，伯努利试验是一个只有两种可能结果的单一实验，这些结果也可以标记为“成功”和“失败”。它是最简单的离散概率分布，并作为更复杂分布（如二项分布和几何分布）的基础。

伯努利分布广泛应用于具有二元结果的建模场景中，例如抛硬币、是/否调查问题，或在数据集中存在/不存在特定特征的场景。例如，在统计假设检验中，伯努利分布常用于比较两种治疗某种医疗状况的成功率。在金融领域，伯努利分布可以用来模拟二元结果，例如股票价格上升或下降。

按照惯例，伯努利分布中的两种结果通常编码为成功为 1，失败为 0。伯努利分布由一个参数 p 表示，它代表成功的概率。换句话说，我们有 f(x = 1) = p。同样，由于总概率加起来为 1，失败的概率由 1 − p 给出，即 f(x = 0) = 1 − p。

我们可以将伯努利分布的概率质量函数（PMF）表示如下：

f(x) = {p, if x = 1 1 − p, if x = 0

或者，我们可以将 f(x) 以更紧凑的形式表示，如下所示。很容易验证这两种表示是等价的：

f(x = i) = p^i (1 − p)^(1−i) 对于 i ∈ {0,1}

注意到 p ∈ [0,1]。

对于伯努利分布的均值，μ（第一矩）和方差，σ²（第二矩），我们有以下结果：

μ = p

σ² = p(1 − p)

在下面的练习中，我们将使用 R 模拟和分析伯努利分布的随机变量。

练习 10.1 – 模拟和分析伯努利分布的随机变量

在这个练习中，我们将使用 rbinom() 函数模拟和分析伯努利分布的随机变量：

1. 模拟一个成功概率为 0.6 的单个伯努利试验：

   
   # The probability of success
   p = 0.6
   # Produce a random Bernoulli outcome
   outcome = rbinom(1, size = 1, prob = p)
   >>> print(outcome)
   0
   

   这里，结果将显示为 0 或 1。我们可以控制随机种子以确保结果的再现性：

   set.seed(8)
   >>> rbinom(1, size = 1, prob = p)
   1
   

2. 生成五个具有相同成功概率的随机伯努利结果：

   
   # Number of experiments
   n = 5
   # Generate corresponding outcomes
   outcomes = rbinom(n, size = 1, prob = p)
   >>> print(outcomes)
   1 0 0 1 0
   

3. 计算伯努利分布的均值和方差：

   
   # Get mean and variance
   mean_bernoulli = p
   var_bernoulli = p * (1 - p)
   >>> cat("Mean:", mean_bernoulli, "\nVariance:", var_bernoulli)
   Mean: 0.6
   Variance: 0.24
   

   在这里，我们使用 cat() 函数连接并打印结果。

4. 从观察/实证成功的概率来分析多次伯努利试验的结果：

   
   # Number of successes
   num_successes = sum(outcomes)
   # Empirical probability of success
   empirical_p = num_successes / n
   >>> cat("Number of successes:", num_successes, "\nEmpirical probability of success:", empirical_p)
   Number of successes: 2
   Empirical probability of success: 0.4
   

   由于我们只采样了 5 次，成功的结果的实证概率 (0.4) 与真实概率 (0.6) 相差甚远。我们可以增加随机试验的规模以获得更可靠的估计：

   n = 1000
   num_successes = sum(rbinom(n, size = 1, prob = p))
   empirical_p = num_successes / n
   >>> cat("Number of successes:", num_successes, "\nEmpirical probability of success:", empirical_p)
   Number of successes: 600
   Empirical probability of success: 0.6
   

   使用总共 1000 次试验，我们现在可以重现成功的确切概率。

下一个部分回顾了二项分布。

二项分布

1 或 0。具体来说，它是一个离散概率分布，指定了在给定数量的伯努利试验中的成功次数。这些试验是独立的，并且具有相同成功的概率。

二项分布 PMF 有两个参数：

• 试验次数，n

• 每次试验中成功的概率，p

注意

我们仍然假设只有两种可能的结果：成功 (1) 或失败 (0)。失败的概率也可以表示为 q，其中 q = 1 − p。

总共有 k 次成功的二项分布的 PMF 如下公式给出：

P(x = k) = C(n, k) p^k (1 − p)^(n−1)

这里，C(n, k) 表示从 n 次试验中选择 k 次成功的组合数，可以使用二项式系数公式计算：

C(n, k) = n! / (k! * (n − k)!)

前两个矩如下：

μ = np

σ² = np(1 − p)

使用二项分布的 PMF，我们可以计算在特定次数的试验中观察特定次数成功的概率。

二项分布与其他概率分布有一些重要关系。例如，当 n 趋向于无穷大且 p 保持不变时，二项分布收敛到正态分布，这是一种将在后面介绍的连续概率分布。

让我们通过一个练习来熟悉与二项分布相关的函数。

练习 10.2 – 模拟和分析二项随机变量

在这个练习中，我们将使用 dbinom() 和 pbinom() 函数模拟和分析二项分布的随机变量：

1. 使用 dbinom() 函数根据成功概率为 0.5 和总共 10 次试验的二项分布计算观察 0 到 10 次成功的概率：

   
   n = 10 # Number of trials
   p = 0.5 # Probability of success
   # Get binomial probabilities for different occurrences of successes
   binom_probs = dbinom(0:n, n, p)
   >>> binom_probs
   [1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000
    [5] 0.2050781250 0.2460937500 0.2050781250 0.1171875000
    [9] 0.0439453125 0.0097656250 0.0009765625
   

   在这里，我们使用 0:n 来创建一个从 0 到 10 的整数列表，每个整数随后将传递给 dbinom() 函数以评估相应的概率。

2. 使用 barplot() 函数创建二项概率的条形图：

   
   >>> barplot(binom_probs, names.arg = 0:n, xlab = "Number of Successes", ylab = "Probability", main = "Binomial Distribution (n = 10, p = 0.5)")
   

   运行此代码生成 图 10*.2*，显示中间发生 5 次的概率最高：

图 10.2 – 使用 n=10 和 p=0.5 可视化二项分布

1. 使用 pbinom() 函数计算累积二项概率：

   
   cum_binom_probs <- pbinom(0:n, n, p)
   >>> cum_binom_probs
   [1] 0.0009765625 0.0107421875 0.0546875000 0.1718750000
    [5] 0.3769531250 0.6230468750 0.8281250000 0.9453125000
    [9] 0.9892578125 0.9990234375 1.0000000000
   

   结果显示，累积二项概率是从 CDF 中计算出来的，它是 PMF 中前一个元素概率的累积和。

   我们还可以使用 CDF 来计算观察特定值或更高的概率。

2. 计算至少获得七次成功的概率：

   
   prob_at_least_7_successes = 1 - pbinom(6, n, p)
   >>> prob_at_least_7_successes
   0.171875
   

   在这里，通过观察 6 次或更少的成功来计算至少获得 7 次成功的概率，即我们有 P(X ≥ 7) = 1 − P(X ≤ 6)。

让我们通过另一个与实际应用相关的练习来将这些计算置于适当的背景中。

练习 10.3 – 计算获胜概率

在这个练习中，我们将计算体育队的获胜概率：

1. 假设体育队赢得比赛的概率为 80%。如果总共有五场比赛，那么赢得至少四场比赛的概率是多少？

   
   n = 5
   p = 0.8
   prob_at_least_4_wins = 1 - pbinom(3, n, p)
   >>> prob_at_least_4_wins
   0.73728
   

2. 计算最多赢得三场比赛的概率：

   
   prob_at_most_3_wins = pbinom(3, n, p)
   >>> prob_at_most_3_wins
   0.26272
   

   注意，这个概率是赢得至少四场比赛的前一个概率的补数。我们可以如下验证这个关系：

   >>> prob_at_most_3_wins == 1 – prob_at_least_4_wins
   TRUE
   

在下一节中，我们将暂停讨论二项分布的正态近似。这种关系在统计分析中得到广泛应用，其根源在于中心极限定理。

二项分布的正态近似

正态分布（或高斯分布）对二项分布的近似表明，二项分布可以被一个具有相同均值（μ = np）和方差值（σ² = np(1 - p)）的正态分布所近似。这种正态近似随着实验次数（n）的增加和成功概率（p）不接近 0 或 1 而变得更加准确。作为一个经验法则，我们通常在 np ≥ 10 和 nq ≥ 10 时使用正态近似。

要使用正态近似，我们需要通过以下公式将二项随机变量 x 标准化，将其转换为标准正态变量 z 的形式：

z = x - μ / σ

这里，μ = np，σ = √np(1 - p)。这也被称为z 分数。在尝试将不同数量比较在相同尺度上时，进行这种标准化是一种常见做法。我们可以利用标准正态分布（稍后介绍）来处理相应的 PDF 或 CDF，具体取决于具体任务。在这种情况下，我们可以使用标准正态分布（即 N(0,1)）来近似与二项分布相关的概率。

让我们来看一个具体的例子。假设我们抛掷一枚硬币 100 次（假设是公平的硬币，正面和反面的概率相等）并希望计算得到 40 到 60 个头（包括两端）的概率。令 x 表示头数。我们知道 x 服从参数为 n = 100 和 p = 0.5 的二项分布。

为了检查正态近似是否合适，我们计算 np = 100 * 0.5 = 50 > 10 和 np(1 - p) = 100 * 0.5 * 0.5 = 25 > 10。我们还可以使用 R 进行验证，如下面的代码片段所示：

n = 100
p = 0.5
# check conditions for normal approximation
>>> n*p > 10
TRUE
>>> n*p*(1-p) > 10
TRUE

两个条件都评估为TRUE。现在，我们可以将上限和下限（分别为60和40）标准化，将它们转换为标准化分数。为此，我们需要获得二项分布的参数，然后应用标准化公式，z = x - μ / σ，以获得 z 分数。以下代码片段完成了标准化：

# compute mean and std
mu = n*p
>>> mu
50
std = sqrt(n*p*(1-p))
>>> std
5
# compute P(lower_limit <= X <= upper_limit)
lower_limit = 40
upper_limit = 60
# Using z score
standard_lower_limit = (lower_limit – mu) / std
standard_upper_limit = (upper_limit – mu) / std
>>> standard_lower_limit
-2
>>> standard_upper_limit
2

使用标准化后的 z 分数，我们现在可以仅基于标准正态分布来计算原始概率。换句话说，我们有以下：

P(40 ≤ x ≤ 60) = P(40 - 50 / 5 ≤ x - 50 / 5 ≤ 60 - 50 / 5) = P(-2 ≤ z ≤ 2)

如以下代码片段所示，我们现在可以调用pnorm()函数来计算 z = -2 和 z = 2 处的累积分布函数（CDF），它们的差值给出了最终的概率：

# approximate using standard normal cdf
>>> pnorm(standard_upper_limit) - pnorm(standard_lower_limit)
0.9544997

让我们也计算使用二项分布对应的概率，看看正态近似有多接近。以下代码片段使用pbino()函数在两个极限处获得二项分布的 CDF，然后取差值给出在这个范围内的观察结果的总体概率：

# use binomial distribution
>>> pbinom(upper_limit, n, p) - pbinom(lower_limit, n, p)
0.9539559

结果显示，近似值精确到小数点后第二位。因此，归一化近似法为计算概率提供了一种替代方法，如果直接使用二项分布不方便的话。然而，尽管正态近似是一个强大的工具，我们仍然需要检查所需条件以确保良好的近似。

图 10.3 总结了计算观察特定值范围的总体概率的两种方法。我们可以通过计算范围两个边界之间的累积分布函数（CDF）的差来计算它。或者，我们可以依赖正态分布来近似二项分布，并在满足条件的情况下，根据这些参数获得标准化的 z 分数后计算总体概率：

图 10.3 – 总结在计算观察特定值范围的总体概率时对二项分布的正态近似

我们将在下一节回顾泊松分布。

泊松分布

另一个流行的离散概率分布也是泊松分布，它描述了在固定时间或空间区间内的事件数量。它有一个单一常数参数，指定了发生平均速率。具体来说，我们必须将λ表示为给定区间内事件发生的平均速率。我们可以将泊松分布的概率质量函数（PMF）表示如下：

P(x = k) = λ^k e^(-λ) / k!

在这里，P(x = k)表示在固定区间内经历总共 k 次发生的概率，e 是欧拉数（约等于 2.71），而 k!是 k 的阶乘（所有正整数乘积到 k）。

使用这个公式，我们可以计算在给定固定区间内观察到任何（整数）事件数量的概率。我们还可以进一步计算泊松分布的均值和方差如下：

μ = λ

σ² = λ

泊松分布通常用于模拟独立发生且平均发生速率恒定的事件。泊松过程的一些实际应用包括每小时前台收到的酒店预订数量和一天内收到的电子邮件数量。

让我们通过一个练习来熟悉与泊松分布相关的常见概率计算。

练习 10.4 – 模拟和分析泊松分布的随机变量

在这个练习中，我们将使用 R 来处理泊松分布，包括计算概率（使用概率质量函数 PMF）、绘制分布图和生成随机样本。具体来说，我们将计算平均发生率为每间隔 5 个事件的泊松概率（λ = 5），绘制概率质量函数 PMF，并计算累积概率。我们还将从这个泊松分布中生成 100 个观察值的随机样本：

1. 根据λ = 5 的泊松分布参数，使用 dpois() 函数计算每个间隔内观察到 0 到 15 个发生/事件的概率：

   
   lambda = 5 # distribution parameter
   # Calculate probabilities for each scenario
   pois_probs = dpois(0:15, lambda)
   >>> pois_probs
   [1] 0.0067379470 0.0336897350 0.0842243375 0.1403738958
    [5] 0.1754673698 0.1754673698 0.1462228081 0.1044448630
    [9] 0.0652780393 0.0362655774 0.0181327887 0.0082421767
   [13] 0.0034342403 0.0013208616 0.0004717363 0.0001572454
   

2. 创建泊松概率的条形图：

   
   >>> barplot(pois_probs, names.arg = 0:15, xlab = "Number of Events", ylab = "Probability", main = "Poisson Distribution (lambda = 5)")
   

   运行此代码生成 图 10*.4*。如预期，峰值概率出现在 5 次左右：

图 10.4 – 将泊松分布的概率质量函数 PMF 以条形图的形式可视化

1. 使用 ppois() 函数计算每个事件数（0 到 15）的累积泊松概率：

   
   cum_pois_probs = ppois(0:15, lambda)
   >>> cum_pois_probs
   [1] 0.006737947 0.040427682 0.124652019 0.265025915
    [5] 0.440493285 0.615960655 0.762183463 0.866628326
    [9] 0.931906365 0.968171943 0.986304731 0.994546908
   [13] 0.997981148 0.999302010 0.999773746 0.999930992
   >>> barplot(cum_pois_probs, names.arg = 0:15, xlab = "Number of Events", ylab = "Cumulative Probability", main = "CDF of Poisson Distribution (lambda = 5)")
   

   运行此代码生成 图 10*.5*。注意，CDF 曲线在平均发生 5 次附近急剧上升，并向图形的右侧逐渐饱和：

图 10.5 – 将泊松分布的累积分布函数 CDF 以条形图的形式可视化

1. 使用 rpois() 函数从这个泊松分布中生成 100 个随机样本：

   
   pois_samples = rpois(100, lambda)
   >>> pois_samples
     [1]  8  5  8  4  3  4  6  2  5  6  3  3  7  8  8  7
    [17]  5  9  6  1  4  2  7  7  5  5  5  2  7  4  6  5
    [33]  4  4  3  0  8  5  4  4  7  5 11  6  5  4  8  8
    [49]  2  5  6  2  3  4  6  4  6  2  5  3  6  0  5  8
    [65]  7  1  8  4  4  4  4  5  4  4  4  5  5  6  3  4
    [81]  3  0  8  9  2  3  4 13  2  6  8  9  6  4  7  7
    [97]  8  6  3  5
   

   如预期，大多数发生次数都在 5 次左右。

泊松分布对二项分布的近似

结果表明，我们还可以在特定条件下使用泊松分布来近似二项分布。例如，当二项分布中的试验次数（n）很大且成功概率（p）很小时，二项分布可以被近似为λ = np 的泊松分布。

让我们用一个例子来演示如何将泊松近似应用于二项分布。假设我们有一个 n = 1000 和 p = 0.01 的二项分布。我们想找到观察到恰好 15 次成功的概率：

1. 我们可以从二项概率开始，指定参数后使用 dbinom() 函数计算相应的概率：

# Binomial parameters
n = 1000
p = 0.01
# Probability of observing 15 successes
binom_prob = dbinom(15, n, p)
>>> binom_prob
0.03454173

接下来，我们必须计算近似的泊松参数，λ = np:

lambda_approx = n * p
>>> lambda_approx
10

现在，我们可以计算观察到 15 次成功的泊松概率：

pois_approx_prob <- dpois(15, lambda_approx)
>>> pois_approx_prob
0.03471807

结果表明，近似值在第三位小数点处相当准确。

另一个有趣的性质是，将几个独立的泊松分布随机变量相加也会产生泊松分布，其参数为相应单个λ值的总和。例如，如果 x₁和 x₂是分别具有λ₁和λ₂的独立泊松随机变量，它们的和（y = x₁ + x₂）也遵循λ = λ₁ + λ₂的泊松分布。这在处理多个泊松分布随机变量的和时提供了一个方便的性质。

在下一节中，我们将介绍另一个广泛使用的离散分布：几何分布。

几何分布

几何分布是一种离散概率分布，描述了在一系列独立伯努利试验中，每个试验具有相同的成功概率，第一次成功所需的试验次数。与二项分布类似，几何分布是多个独立伯努利试验的集合，尽管感兴趣的是试验序列中成功的第一次发生。在这里，成功的第一次发生意味着所有之前的试验都需要是非成功，而当前试验是在迄今为止进行的多个试验中第一次成功。

它通常用于模拟事件发生前的等待时间或达到期望结果所需的尝试次数。例如，包括通过驾驶考试所需的尝试次数直到成功、连续晴天观察的次数以及获得第一次正面所需的抛硬币次数。在模拟一系列独立伯努利试验中第一次成功所需的等待时间或尝试次数时，几何分布非常有用。

几何分布由一个参数 p 定义，该参数代表每次伯努利试验成功的概率。几何分布的概率质量函数（PMF）由以下公式给出：

P(x = k) = (1 − p)^(k−1) * p

此公式指定了在 k 次试验中观察到第一次成功的概率。因此，该概率是通过观察 k 个独立事件的联合概率来计算的，其中前 k-1 个事件是非成功，联合概率为(1 − p)^(k−1)，最后一个事件是成功，概率为 p。

几何分布的均值和方差参数可以表示如下：

μ = 1 / p

σ² = (1 − p) / p²

注意，几何分布是无记忆的，这意味着下一次试验成功的概率不依赖于过去的试验。换句话说，直到第一次成功所需的等待时间不会因已经进行的过去试验的数量而改变。

让我们通过一个练习来模拟和分析一个遵循几何分布的随机变量。

练习 10.5 – 模拟和分析几何分布的随机变量

在这个练习中，我们将使用 R 来处理几何分布，包括计算概率质量函数（PMF）和累积分布函数（CDF）概率，绘制分布，并生成随机样本：

1. 对于成功概率为 p = 0.25 的几何分布，使用 dgeom() 函数计算每个试验次数（从 1 到 10）的几何概率：

   
   # Parameters
   p = 0.25 # Probability of success
   # Get geometric probabilities
   geom_probs = dgeom(0:9, p)
   >>> geom_probs
   [1] 0.25000000 0.18750000 0.14062500 0.10546875 0.07910156 0.05932617
    [7] 0.04449463 0.03337097 0.02502823 0.01877117
   

   注意，0:9 参数表示 dgeom() 函数中的失败次数。

2. 为这些几何概率创建条形图：

   
   >>> barplot(geom_probs, names.arg = 1:10, xlab = "Number of Trials", ylab = "Probability", main = "Geometric Distribution (p = 0.25)")
   

   运行此代码生成图 10.6*. 如预期，随着试验次数的增加，获得更长的连续失败序列的概率降低：

图 10.6 – 将几何分布的概率质量函数（PMF）作为条形图可视化

1. 使用 pgeom() 函数计算之前试验的累积几何概率：

   
   cum_geom_probs = pgeom(0:9, p)
   >>> cum_geom_probs
   [1] 0.2500000 0.4375000 0.5781250 0.6835938 0.7626953
    [6] 0.8220215 0.8665161 0.8998871 0.9249153 0.9436865
   

2. 使用 rgeom() 函数从这个几何分布中生成 100 个随机样本：

   
   geom_samples = rgeom(100, p)
   >>> geom_samples
     [1]  0  0  0  2 10  1  1 10  0  0  1  1  5  3  1  0  2  0  0
    [20]  4  0  1  4  2  3  2  2  2  4  1  6 12  4  1  7  3  1  1
    [39]  0  2  1  2  3  0  8  0  0  2 10  3  2  8  0  3  1  2  3
    [58]  0  0  1  7  0  0  3  4 11  8  8  2  0  5  1  1  1  3  1
    [77]  3  1  3  3  6  0  0  7  1  0  0  1  0  1  0  0  0  3  0
    [96]  0  4 25  0  3
   

   结果显示，随着数字增大，频率逐渐降低，这与几何分布的概率质量函数（PMF）相匹配。

让我们通过另一个与概率相关的练习来查找计算机程序中的错误。

练习 10.6 – 模拟和分析几何分布的随机变量

在这个练习中，我们将访问一个涉及软件测试员尝试在计算机程序中查找错误的现实世界示例。在这个例子中，软件测试员每次尝试找到错误的概率为 0.1，且尝试是独立的。公司想知道在第一次五次尝试内找到第一个错误的概率，以及找到第一个错误所需的预期尝试次数：

1. 使用 pgeom() 函数计算在第一次五次尝试内找到第一个错误的概率：

   
   p = 0.1 # Probability of finding a bug on each attempt
   # Calculate the CDF for up to 5 attempts
   prob_within_5_attempts = pgeom(4, p)
   >>> prob_within_5_attempts
   0.40951
   

   这里，请注意我们使用 4，因为 dgeom() 函数使用零基索引。

   由于 pgeom() 函数返回特定输入的累积分布函数（CDF），我们可以通过将所有之前的概率加到当前输入上等效地计算这个概率，如下所示：

   >>> sum(dgeom(0:4, p))
   0.40951
   

2. 使用几何分布的均值（期望值）计算找到第一个错误所需的预期尝试次数：

   
   mean_attempts <- 1 / p
   >>> mean_attempts
   10
   

3. 在条形图中可视化在不同尝试次数（从 1 到 20）内找到第一个错误的概率：

   
   geom_probs <- dgeom(0:19, p)
   # Create a bar plot of probabilities
   barplot(geom_probs, names.arg = 1:20, xlab = "Number of Attempts", ylab = "Probability", main = "Geometric Distribution (p = 0.1)")
   

   运行此代码生成图 10.7*. 此图表明，随着我们向右移动，需要连续失败流的事件（那些具有较低概率的事件）的概率更低：

图 10.7 – 在不同尝试次数内可视化找到第一个错误的概率

如此图表所示，观察到第一个错误的概率随着尝试次数的增加而降低。

比较不同的离散概率分布

本章迄今为止介绍的离散概率分布是模拟结果变量取离散值场景的必要工具。每个离散分布都有特定的假设、属性和应用。图 10.8 通过比较二项分布、泊松分布和几何分布的主要特征提供了总结和解析：

图 10.8 – 总结和比较不同的离散分布

在手头有不同离散分布的情况下，了解每个分布的具体假设和要求，对于选择给定问题的合适分布非常重要。例如，二项分布适用于模拟在给定次数的实验中成功的次数，泊松分布适用于模拟在固定期间内发生的事件数，几何分布通常用于模拟达到第一次成功所需的试验次数。

在实践中，这些离散概率分布可以用来分析各种现实世界场景，进行预测和优化流程。了解每个分布的特征和应用，可以帮助你为特定问题选择合适的分布，并使用 R 进行相关分析。

下一个部分介绍了连续分布，包括正态分布、指数分布和均匀分布。

发现常见的连续概率分布

连续概率分布用于模拟随机变量在特定连续范围内取任何值的概率。换句话说，基础随机变量是连续的而不是离散的。这些分布描述了观察到的值落在连续区间内的概率，而不是在离散概率分布中等于单个离散结果的概率。具体来说，在连续概率分布中，随机变量等于任何特定值的概率通常为零，因为可能的结果是不可数的。相反，连续分布的概率是针对区间或值范围的。

我们可以使用 PDF 来描述连续分布。这对应于离散概率分布的 PMF。PDF 定义了在给定点周围无限小间隔内观察到一个值的概率。PDF 曲线在特定范围内的面积代表随机变量落在该范围内的概率——也就是说，通过在所需范围内对 PDF 进行积分来计算区间或值范围的概率。相比之下，PMF 为单个点分配概率，通过将它们的个别概率相加来计算一组离散值的概率。

此外，连续概率分布的可视化也有所不同。与用于离散概率分布 PMF 的条形图相比，连续分布的 PDF 以称为密度图的平滑曲线绘制。曲线在特定范围内的面积代表随机变量落在该范围内的概率。

图 10.9总结了离散和连续概率分布之间的主要差异：

图 10.9 – 总结离散和连续概率分布之间的差异

总结来说，离散和连续概率分布模型不同类型的随机变量。离散概率分布表示可数的输出结果，而连续分布表示在连续范围内不可数的可能性。这两种分布类型之间的差异决定了我们如何为特定问题选择合适的分布并执行相关分析。

在下一节中，我们将介绍最广泛使用的连续概率分布：正态概率分布。

正态分布

正态概率分布，也称为高斯分布，是一种连续概率分布，它模拟了连续结果围绕均值对称分布的场景。由于许多自然和社会现象由于中心极限定理的结果往往遵循正态分布，因此它是实践中最广泛使用的概率分布。

两个参数用于描述正态分布：均值（μ）和标准差（σ）。均值代表分布的中心趋势或平均值。分布中心的值获得最高的概率。标准差描述数据从均值的离散程度，作为分布变异性的度量。

正态分布的概率密度函数（PDF）如下所示：

f(x) = 1 / √(2πσ) e^(-(x−μ)² / (2σ²))

我们也可以写成 x ∼ N(μ, σ²)，这表示随机变量 x 遵循由 μ 和 σ² 参数化的正态分布。

从图形上看，正态分布看起来像钟形曲线，大多数值集中在均值附近，两端极端的值较少。一个经验法则，称为 68-95-99.7 规则，表示大约 68% 的值位于均值加减一个标准差范围内，95% 在加减两个标准差范围内，99.7% 在加减三个标准差范围内。

一种常用的特定正态分布是标准正态分布，表示为 z ∼ N(0, 1) – 即，标准正态分布具有 μ = 0 和 σ = 0。一个特殊性质是我们可以使用以下公式将任何正态分布的随机变量转换为标准正态变量：

z = x − μ / σ

然后，我们可以使用标准正态分布表（称为 Z 表）来找到感兴趣的概率和分位数。

让我们通过一个练习来练习与正态分布相关的计算。

练习 10.7 – 模拟和分析正态随机变量

在这个练习中，我们将模拟和分析正态分布的随机变量：

1. 使用 dnorm() 函数计算标准正态分布从 -4 到 4 的概率密度，步长为 0.1：

   
   # Parameters
   mu = 0      # Mean
   sigma = 1   # Standard deviation
   # Get the probability density for different x
   x = seq(-4, 4, by = 0.1)
   normal_density = dnorm(x, mu, sigma)
   

   在这里，我们使用 seq() 函数创建一个等间距的值向量，并使用 dnorm() 函数提取每个输入值的对应概率。

2. 使用 plot() 函数将正态分布绘制为连续曲线：

   
   # Plot the normal distribution
   >>> plot(x, normal_density, type = "l", xlab = "x", ylab = "Probability Density", main = "Normal Distribution (μ = 0, σ = 1)")
   

   运行此代码生成 图 10*.10*，它显示 PDF 以均值 0 为中心，具有标准差 1 的分布范围：

图 10.10 – 标准正态分布密度图的可视化

1. 使用 pnorm() 函数计算正态分布的累积概率：

   
   # Get cumulative probabilities for different x
   normal_cum_prob <- pnorm(x, mu, sigma)
   >>> plot(x, normal_cum_prob, type = "l", xlab = "x", ylab = "Cumulative Probability Density", main = "Cumulative Normal Distribution (μ = 0, σ = 1)")
   

   运行此代码生成 图 10*.11*：

图 10.11 – 标准正态分布累积密度函数的可视化

1. 使用 rnorm() 函数从正态分布中生成随机样本：

   
   # Generate 100 random samples from a normal distribution with μ = 0 and σ = 1
   normal_samples <- rnorm(100, mu, sigma)
   

2. 使用 qnorm() 函数找到一个给定概率的 90% 分位数（逆累积概率）：

   
   # Find the quantile corresponding to the 90th percentile
   quantile_90 <- qnorm(0.9, mu, sigma)
   >>> quantile_90
   1.281552
   

让我们看看另一个使用正态分布解决实际问题的练习。

练习 10.8 – 使用正态分布计算概率

假设一家公司生产的电池平均寿命为 100 小时，标准差为 10 小时。假设电池的寿命遵循正态分布：

1. 模拟一个包含 1,000 个电池的数据集：

   
   set.seed(8)
   mean_lifespan = 100
   sd_lifespan = 10
   n = 1000
   lifespans = rnorm(n, mean_lifespan, sd_lifespan)
   

   在这里，我们使用rnorm()函数从给定的正态分布中随机抽样。我们还指定了随机种子以确保可重复性。

2. 计算随机选择的一块电池使用时间超过 120 小时的概率：

   
   threshold = 120
   probability = 1 - pnorm(threshold, mean_lifespan, sd_lifespan)
   >>> probability
   0.02275013
   

   在这里，我们使用pnorm()函数计算小于120的总概率，然后取补数以获得大于120的概率。

   如预期，偏离平均值两个标准差外的概率相当小。

3. 绘制寿命概率密度函数，曲线下方的面积大于120小时的区域进行阴影处理：

   
   df <- data.frame(lifespan = lifespans)
   df_density <- density(lifespans)
   df_shaded <- data.frame(x = df_density$x, y = df_density$y)
   df_shaded <- df_shaded[df_shaded$x > threshold,]
   ggplot(df, aes(x=lifespan)) +
     geom_density(fill="lightblue") +
     geom_vline(xintercept = threshold, linetype="dashed", color="red") +
     geom_area(data = df_shaded, aes(x=x, y=y), fill="orange", alpha=0.5) +
     theme_minimal() +
     labs(title="Lifespan of batteries", x="Lifespan (hours)", y="Probability Density")
   

   在这里，我们构建一个 DataFrame，df，用于存储样本值和通过density()函数获得的相应密度。然后我们子集以获取要阴影处理的区域的相应 DataFrame，df_shaded。在ggplot中，我们使用geom_density()函数绘制密度曲线，geom_vline()添加表示阈值的垂直线，以及geom_area()对右侧区域进行阴影处理。

   运行此代码将生成图 10.12：

图 10.12 – 使用阴影区域可视化经验正态分布的概率密度函数

下一节将介绍另一种连续概率分布：指数分布。

指数分布

指数分布是一种连续概率分布，常用于模拟泊松过程中事件之间的时间或空间。如前所述，泊松过程模拟独立且以恒定平均速率发生的事件。指数分布常用于描述罕见事件之间的等待时间，例如呼叫中心的电话之间的时间。

指数分布的概率密度函数由以下公式给出：

f(x) = λ e^(-λx), x ≥ 0

在这里，λ是速率参数，表示每单位时间或空间发生的事件的平均数量，e 是自然对数的底数。

指数分布的一个定义特征是记忆性属性，它表明未来事件发生的概率与自上次事件以来已经过去的时间无关。这使得指数分布适合用于模拟独立且罕见事件之间的等待时间。

指数分布的均值和标准差如下：

μ = 1 / λ

σ² = 1 / λ²

指数分布包含一个参数 λ，并模拟事件之间等待时间每单位时间或空间发生的平均事件数。在模拟泊松过程的情况下，相同的参数 λ 指的是在固定间隔内（无论是时间还是空间）发生的平均事件数。这两个分布都用于模拟泊松过程的各个方面。

现在，让我们通过一个练习来复习与指数分布相关的概率计算。

练习 10.9 – 使用指数分布计算概率

在这个练习中，我们将探讨在遵循指数分布的同时生成随机样本，并计算和可视化超过某个阈值的总概率：

1. 使用 rexp() 函数从具有速率参数 0.01 的指数分布中生成 1,000 个数据点的随机样本：

   
   set.seed(8) # Set seed for reproducibility
   lambda = 0.01
   sample_size = 1000
   exponential_sample = rexp(sample_size, rate = lambda)
   

2. 使用 pexp() 函数计算事件之间的等待时间超过 150 个单位时的概率：

   
   threshold = 150
   probability_above_threshold = 1 - pexp(threshold, rate = lambda)
   >>> probability_above_threshold
   0.2231302
   

3. 绘制 PDF 并阴影表示超过阈值的等待时间下的曲线区域：

   
   # Create a data frame for the waiting times
   waiting_times = seq(0, max(exponential_sample), length.out = 1000)
   density_values = dexp(waiting_times, rate = lambda)
   df = data.frame(waiting_times, density_values)
   # Filter data for the shaded region
   df_shaded = df[df$waiting_times > threshold,]
   # Plot the PDF of the exponential distribution
   ggplot(df, aes(x = waiting_times, y = density_values)) +
     geom_line() +
     geom_area(data = df_shaded, aes(x = waiting_times, y = density_values), fill = "orange", alpha = 0.5) +
     geom_vline(xintercept = threshold, linetype = "dashed", color = "red") +
     theme_minimal() +
     labs(title = "Exponential Distribution ( = 0.01)", x = "Waiting Time", y = "Probability Density")
   

   在这里，我们创建一个用于存储在随机样本中生成的不同等待时间的 DataFrame，使用 dexp() 函数获取相应的密度，并构建 df 和 df_shaded DataFrame 以供 ggplot 使用。

   运行此代码生成 图 10*.13*：

图 10.13 – 以阴影表示阈值以上的区域来可视化指数分布的经验概率密度函数

下一节将介绍均匀分布。

均匀分布

如其名所示，均匀分布是一种连续概率分布，其中给定范围内的所有结果出现的可能性相同。PDF（如图所示为一条直线）由两个参数定义，即下限（a）和上限（b），它们定义了可能值的范围。此范围内的所有值具有相同的发生概率，而范围外的值具有零概率。

均匀分布的 PDF 由以下公式给出：

f(x) = { 1 _ b − a , if x ∈ [a, b]  0, otherwise

均匀分布的参数如下：

μ =  a + b _ 2

σ 2 =  (b − a) 2 _ 12

让我们通过一个类似的练习来分析均匀分布的随机变量。

练习 10.10 – 使用均匀分布计算概率

在这个练习中，我们将探讨生成随机样本、计算概率和绘制均匀分布的 PDF：

1. 使用 runif() 函数从下限 (a) 为 2 和上限 (b) 为 10 的均匀分布中生成 10,000 个数据点的随机样本：

   
   set.seed(8) # Set seed for reproducibility
   a = 2
   b = 10
   sample_size = 10000
   uniform_sample = runif(sample_size, min = a, max = b)
   

2. 使用punif()函数计算从均匀分布中选出的值大于7的概率：

   
   threshold = 7
   probability = 1 - punif(threshold, min = a, max = b)
   >>> probability
   0.375
   

   在这里，我们使用punif()函数计算a和threshold之间的 CDF。我们也可以使用经验样本来近似这个概率：

   probability2 = sum(uniform_sample > t
   hreshold) / length(uniform_sample)
   >>> probability2
   0.3771
   

   我们可以看到，近似相当接近，并且当样本量增大时，它将更加接近。

3. 使用ggplot()绘制其 PDF：

   
   library(ggplot2)
   # Create a data frame for the distribution
   x_values = seq(a, b, length.out = 1000)
   density_values = dunif(x_values, min = a, max = b)
   df = data.frame(x_values, density_values)
   # Plot the PDF of the uniform distribution
   ggplot(df, aes(x = x_values, y = density_values)) +
     geom_line() +
     theme_minimal() +
     labs(title = "Uniform Distribution (a = 2, b = 10)", x = "Value", y = "Probability Density")
   

   在这里，我们在x轴上构建一系列占位符。然后，我们获得相应的密度并将它们组合成一个 DataFrame 进行绘图。运行此代码生成图 10.14：

图 10.14 – 绘制均匀分布的 PDF

均匀分布可以用来生成正态分布的随机样本。让我们看看这是如何实现的。

生成正态分布的随机样本

到目前为止，我们已经学习了如何使用rnorm()函数从高斯分布中采样。这有助于了解这些随机样本是如何在幕后生成的。这种特定的技术被称为逆变换法，它依赖于目标分布（在这种情况下，正态分布）的逆 CDF 或分位数函数。

此过程涉及三个步骤。首先，我们将从均匀分布中生成随机样本，通常是 U(0,1)，其中0和1之间的所有值都是等可能的。接下来，对于每个均匀随机样本，我们将使用标准正态分布的逆 CDF（分位数函数）在目标分布中定位相应的值。最后，我们将应用缩放-定位变换将标准正态随机样本转换为目标正态分布的随机样本。

此方法依赖于连续随机变量的 CDF 是一个将样本空间（PDF 的定义域）映射到[0, 1]范围的函数。逆 CDF 是 CDF 的逆函数，将区间[0, 1]映射回样本空间。通过将均匀随机样本传递给逆 CDF，我们反转映射并得到遵循目标分布的随机样本，前提是进行额外的变换。

让我们更详细地看看这个过程。假设我们想要从正态分布 N(μ, σ²)中采样。我们如何从这个特定的分布中生成随机样本？一种方法是从标准正态分布 N(0,1)中生成一个随机样本 x，然后应用缩放-定位变换以获得最终的样本σx + μ。通过基于σ缩放样本，然后加上μ，得到的样本将遵循具有均值μ和方差σ²的正态分布。

从标准正态分布中采样的第一步是关键，而第二步是一个简单且确定的变换。我们之前介绍的方法是使用标准正态分布的逆 CDF 将均匀分布的变量进行变换。例如，如果 U 在[0,1]区间上均匀分布，那么Φ −1(U)遵循 N(0,1)，其中Φ −1 是标准正态分布累积函数的逆。

一个说明性的例子在图 10**.15中展示。首先，我们从[0,1]区间上的均匀分布中随机采样一个点。接下来，我们使用标准正态分布的逆累积分布函数（CDF）来获得 CDF 空间中的对应样本，考虑到 CDF 单调地将任意输入值映射到[0,1]区间的输出。从数学上讲，标准正态分布的随机样本由 x = Φ −1(U) 给出。考虑到标准正态分布的 PDF 和 CDF 之间的一一映射关系，我们也可以在 PDF 空间中获得相同的输入，x。我们还会期望大多数样本都集中在均值附近。最后，我们应用尺度-位置变换来转换为具有所需均值和方差的正态分布的随机样本：

图 10.15 – 从所需的单变量高斯分布中获得随机样本

让我们通过一个具体的例子来看如何使用逆变换方法生成正态分布的随机样本。在下面的代码片段中，我们首先设置种子以实现可重复性，并定义目标正态分布的参数。然后，我们在 0 和 1 之间遵循均匀分布生成 5 个样本。最后，我们使用正态分布的逆 CDF（分位数函数，qnorm()）计算均匀样本的对应分位数：

set.seed(8) # Set seed for reproducibility
# Define the target normal distribution parameters
mu = 5
sigma = 2
# Generate uniform random variables
n = 5
uniform_sample = runif(n)
# Calculate the corresponding quantiles for the uniform sample using the inverse CDF (quantile function) of the normal distribution
normal_sample = qnorm(uniform_sample, mean = mu, sd = sigma)
>>> normal_sample
[1] 4.830828 3.372006 6.680800 5.780755 4.073034

我们也可以应用尺度-位置变换来获得相同的随机样本，如下面的代码片段所示：

normal_sample2 = qnorm(uniform_sample, mean = 0, sd = 1)
>>> normal_sample2 * sigma + mu
[1] 4.830828 3.372006 6.680800 5.780755 4.073034

下一个部分将涵盖采样分布。

理解常见的采样分布

采样分布是基于从总体中抽取的许多样本的样本统计量的概率分布。换句话说，它是从同一总体的多个样本集中计算出的特定统计量（如均值、中位数或比例）的分布，其中每个集合的大小相同。这里有两点需要注意。首先，采样分布不是关于从 PDF 中抽取的随机样本。相反，它是由来自另一个分布的聚合统计量构成的分布，该分布是从 PDF 中抽取的。其次，我们需要在多个回合中从 PDF 中采样以创建采样分布，其中每个回合包括从 PDF 中抽取的多个样本。

让我们通过 R 中的练习来展示使用样本均值作为感兴趣统计量的抽样分布的概念。我们将从给定分布的总体中生成样本，并计算样本均值。然后，我们将创建样本均值的直方图来可视化样本均值的抽样分布。

练习 10.11 – 生成抽样分布

在这个练习中，我们首先从正态分布生成一个样本总体。然后，我们将从该总体中多次抽取样本，每次抽取包含多个样本。最后，我们将提取每轮样本的均值并将它们一起绘制在直方图中：

1. 从 N(50,10)生成 100,000 个样本。使用summary()检查样本的摘要：

   
   set.seed(8) # Set seed for reproducibility
   # Define the population parameters
   population_mean = 50
   population_sd = 10
   population_size = 100000
   # Generate the population using a normal distribution
   population <- rnorm(population_size, mean = population_mean, sd = population_sd)
   >>> summary(population)
      Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
     7.597  43.261  50.051  50.027  56.781  89.365
   

2. 定义一个函数，从之前的总体中抽取 50 个数字并返回样本的平均值：

   
   # Define the sample size in each round
   sample_size_per_round = 50
   # Function to draw a sample and calculate its mean
   get_sample_mean <- function(population, sample_size_per_round) {
     sample <- sample(population, size = sample_size_per_round, replace = FALSE)
     return(mean(sample))
   }
   

   在这里，我们使用sample()函数从之前创建的样本总体中抽取 50 个数字，不进行替换。我们可以测试这个函数，并观察到每次运行都会得到不同的返回值，这些值也接近于总体均值50：

   >>> get_sample_mean(population, sample_size_per_round)
   50.30953
   >>> get_sample_mean(population, sample_size_per_round)
   48.9098
   

3. 将此函数重复 1,000 次以获得相应的 1,000 个样本均值：

   
   # Generate multiple rounds of sample means
   num_rounds = 1000 # the number of rounds to sample
   sample_means = replicate(num_rounds, get_sample_mean(population, sample_size))
   >>> summary(sample_means)
      Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
     49.76   49.96   50.02   50.03   50.09   50.34
   

   在这里，我们使用replicate()函数重复应用get_sample_mean()函数指定次数。它常用于模拟、重抽样方法或任何需要多次执行相同操作并收集结果的情况。

   结果显示，样本均值从抽样分布中得到的值非常接近总体均值。

4. 使用直方图可视化样本均值的抽样分布：

   
   library(ggplot2)
   sampling_distribution_df = data.frame(sample_means)
   ggplot(sampling_distribution_df, aes(x = sample_means)) +
     geom_histogram(aes(y = after_stat(density)), bins = 30, color = "black", fill = "lightblue") +
     geom_density(color = "red", lwd = 1.2) +
     theme_minimal() +
     labs(title = "Sampling Distribution of the Sample Mean",
          x = "Sample Mean",
          y = "Density")
   

   运行此代码生成图 10.16：

图 10.16 – 使用直方图可视化样本均值的抽样分布

在下一节中，我们将介绍在统计估计和假设检验中常用的抽样分布。

常见的抽样分布

几种常见的抽样分布广泛应用于统计推断。这些分布源于从总体中抽取的随机样本计算出的不同统计量的特性。以下是一些最重要的抽样分布：

• 当n > 30时，样本均值的抽样分布趋近于正态分布，即使原始的总体分布不一定呈正态分布。这一特性允许我们使用正态分布对样本均值进行推断，例如假设检验和置信区间构建。

• 样本比例的抽样分布：在大量独立的伯努利试验（二元结果）中，样本比例（样本组中成功的比例）的分布遵循正态分布，前提是样本量足够大，成功概率不太接近 0 或 1。在这种情况下，均值保持不变，标准差是经过样本量调整的转换后的总体比例。特别是，在一定的条件下（np > 10 和 n(1 − p) > 10，基于我们之前关于通过正态分布近似二项分布的介绍），样本比例的抽样分布可以近似为具有相同均值和标准差σ = √p(1 − p) / n 的正态分布。

• t 分布：t 分布也称为学生 t 分布。当从正态分布的总体中抽取小样本（通常 n<30）估计总体均值，且未提供标准差时使用。t 分布的形状与正态分布相似，但尾部更厚，其确切形状由自由度（df）决定，通常等于 n − 1。t 分布可用于计算 t 分数，当未提供总体标准差时，t 分数用于假设检验和置信区间构建。

• 卡方分布：卡方分布用于假设检验和构建总体方差的置信区间。它还可以用于检验列联表的独立性。它是一族由自由度定义的右偏分布。在检验总体方差的假设检验背景下，检验统计量假设具有 n − 1 个自由度的卡方分布。在列联表的背景下，检验统计量遵循具有(r − 1)(c − 1)个自由度的卡方分布，其中 r 表示表中的行数，c 表示列数。

• F 分布：F 分布用于方差分析（ANOVA）中检验多个组均值相等的零假设。它是一族右偏分布，由两个参数定义，分子（df1）和分母（df2）的自由度。例如，在一元方差分析（ANOVA）的背景下，检验统计量遵循具有df1=k - 1和df2=N - k的 F 分布。

总体而言，这些抽样分布是各种推断过程的重要工具，例如，通过假设检验或构建置信区间，基于样本数据对总体分布的参数进行统计推断。四种类型的抽样分布也是统计推断策略中的更大列表的一个子集，每个都考察不同的假设、样本统计量和推断程序。

让我们看看一个关于使用 t 分布构建总体均值置信区间的练习。

练习 10.12 – 使用 t 分布估计总体均值

在这个练习中，我们将使用一个小样本来估计总体均值，并使用 R 中的 t 分布构建估计的置信区间：

1. 使用rnorm()函数从正态分布 N(50,10)中生成 10 个样本：

   
   sample_size = 10
   mu = 50
   sigma = 10
   samples = rnorm(sample_size, mean = mu, sd = sigma)
   >>> samples
   [1] 41.57424 39.61629 59.86689 58.94655 43.43934 28.41854 67.05759 50.36661 51.61680 37.71842
   

2. 计算样本均值和标准差：

   
   sample_mean = mean(samples)
   sample_sd = sd(samples)
   >>> sample_mean
   47.86213
   >>> sample_sd
   11.85024
   

   现在，我们将使用样本均值和标准差来估计总体均值。

3. 使用 t 分布计算 95%置信区间：

   
   alpha = 0.05
   t_critical = qt(1 - alpha/2, df = sample_size - 1)  # t-value for a two-tailed test with alpha = 0.05 and df = n - 1
   margin_of_error_t = t_critical * (sample_sd / sqrt(sample_size))
   ci_t = c(sample_mean - margin_of_error_t, sample_mean + margin_of_error_t)
   >>> ci_t
   39.38497 56.33928
   

   在这里，我们使用qt()函数找到对应于显著性水平为0.05和自由度 n - 1 的双尾检验的临界 t 值。然后，我们通过将临界 t 值乘以样本均值的标准误差（注意我们需要除以样本大小的平方根）来计算误差范围，并通过从样本均值中加减误差范围来构建置信区间。

我们将在下一章更详细地介绍围绕样本估计构建置信区间的过程。目前，只需了解如何推导和利用抽样分布来产生关于总体分布特征的估计即可。

下一节将介绍另一个有趣且重要的主题：顺序统计量。

理解顺序统计量

顺序统计量是在将样本按升序或降序排列后，样本集合中的值。这些有序样本提供了有关样本数据分布和特征的有用信息。通常，第 k 个顺序统计量是排序样本中的第 k 个最小值。

例如，对于大小为 n 的样本集合，顺序统计量表示为 X1, X2, … , Xn，其中 X1 是最小值（最小值），Xn 是最大值（最大值），而 Xk 代表排序样本中的第 k 个最小值。

让我们看看如何在 R 中提取顺序统计量。

提取顺序统计量

从样本集合中提取顺序统计量可能涉及两种类型的任务。我们可能对以有序方式收集样本感兴趣，这可以通过使用sort()函数实现。或者，我们可能对从有序样本集合中提取特定的顺序统计量感兴趣，例如找到集合中的第三个最大样本或计算集合的特定分位数。

让我们看看在 R 中如何进行此类提取的例子。

练习 10.13 – 提取顺序统计量

在这个练习中，我们将生成一组正态分布的随机样本，并查看对这些样本进行排序以及从中提取特定的顺序统计量：

1. 从均值为50和标准差为10的正态分布中生成 10 个随机样本：

   
   set.seed(8)
   samples = rnorm(10, mean = 50, sd = 10)
   >>> samples
   [1] 49.15414 58.40400 45.36517 44.49165 57.36040 48.92119 48.29711 39.11668 19.88948 44.06826
   

2. 使用sort()函数按升序排序：

   
   sorted_samples = sort(samples)
   >>> sorted_samples
   [1] 19.88948 39.11668 44.06826 44.49165 45.36517 48.29711 48.92119 49.15414 57.36040 58.40400
   

   我们可以看到，这些样本现在是按升序排列的。我们可以通过在sort()函数中设置decreasing = T来切换到降序：

   >>> sort(samples, decreasing = T)
   [1] 58.40400 57.36040 49.15414 48.92119 48.29711 45.36517 44.49165 44.06826 39.11668 19.88948
   

3. 找到最小值（第一阶统计量）：

   
   min_value = sorted_samples[1]
   >>> min_value
   19.88948
   

4. 找到最大值（最后一个顺序统计量）：

   
   max_value = sorted_samples[length(sorted_samples)]
   >>> max_value
   58.404
   

   这里，我们使用样本集合的长度来获取最后一个条目的索引。

5. 找到第三阶统计量（即k = 3）：

   
   k = 3
   kth_order_stat = sorted_samples[k]
   >>> kth_order_stat
   44.06826
   

6. 计算中位数（50 百分位数）：

   
   median_value = median(samples)
   >>> median_value
   46.83114
   

   注意，中位数（或任何其他顺序统计量）在原始样本或排序后的样本中保持不变：

   >>> median(sorted_samples)
   46.83114
   

7. 使用quantile()函数计算第 25 和第 75 百分位数（第一和第三四分位数）：

   
   quartiles = quantile(samples, probs = c(0.25, 0.75))
   >>> quartiles
       25%     75%
   44.1741 49.0959
   

   再次，将相同的函数应用于有序样本会得到相同的结果：

   >>> quantile(sorted_samples, probs = c(0.25, 0.75))
       25%     75%
   44.1741 49.0959
   

在下一节中，我们将介绍顺序统计量在金融中的一个非常重要的用途：风险价值。

计算风险价值

风险价值（VaR）是金融中广泛使用的风险管理指标，它估计在给定置信水平（如 95%或 99%）下，在指定期间内投资组合或投资的潜在损失。它用于量化下行风险并据此分配资本。在这里，置信水平表示潜在损失不会超过计算出的 VaR 的概率。更高的置信水平表示对风险的更保守估计。

计算 VaR 有不同的方法。我们将关注最简单的方法——即使用历史模拟。这种方法使用历史数据来模拟潜在损失。首先，我们必须按升序排序历史回报。然后，VaR 被计算为对应于指定置信水平的回报分位数。这是一个简单直观的方法，尽管它假设历史回报行为代表未来。

让我们通过一个练习来说明如何计算 VaR。

练习 10.14 – 计算 VaR

在这个练习中，我们将讨论如何计算 VaR，这是在给定置信水平下，特定时期内投资组合价值潜在损失的一个度量：

1. 从均值为0.08和标准差为0.05的正态分布中生成一年的每日回报率（252 个交易日）：

   
   # Set a seed for reproducibility
   set.seed(8)
   # Generate a random sample of daily returns from a normal distribution
   sample_size = 252  # Number of trading days in a year
   mu = 0.08          # Mean daily return
   sigma = 0.05       # Standard deviation of daily returns
   daily_returns = rnorm(sample_size, mean = mu, sd = sigma)
   

   我们可以如下检查每日回报率的摘要，它显示每日回报率可能高达 20%，也可能低至-7%。这意味着尽管基础资产平均来说是有利可图的（预期回报率为 8%），但在每日波动方面仍然存在显著风险。VaR 为我们提供了一个衡量在极端情况下风险大小的度量：

   >>> summary(daily_returns)
       Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
   -0.07073  0.04480  0.07926  0.07799  0.11424  0.20195
   

2. 计算总投资组合价值为 100 万美元时的 95%置信水平 VaR：

   
   confidence_level = 0.95
   portfolio_value = 1000000  # Portfolio value in USD
   sorted_returns = sort(daily_returns)
   VaR_index = ceiling(sample_size * (1 - confidence_level))
   VaR = sorted_returns[VaR_index]
   VaR_amount = portfolio_value * (1 - (1 + VaR))
   >>> VaR
   -0.006301223
   >>> VaR_amount
   6301.223
   

   在这里，我们将每日回报率按升序排序，并存储在sorted_returns中。然后，我们在VaR_index中获取底部 5%分位数的索引，然后使用它来检索相应的每日回报率 VaR。最后，我们将百分比回报率转换为投资组合价值的损失VaR_amount。请注意，尽管结果是负数，但我们通常将其报告为正数，以表示在极端情况下可能发生的潜在损失（甚至更多）。

3. 将每日回报率可视化为一个密度图，并阴影 VaR 区域：

   
   library(dplyr)
   daily_returns_df <- data.frame(DailyReturns = daily_returns)
   # Create the density plot
   density_plot <- ggplot(daily_returns_df, aes(x = DailyReturns)) +
     geom_density(fill = "blue", alpha = 0.5) +
     geom_vline(aes(xintercept = VaR), linetype = "dashed", color = "red") +
     labs(x = "Daily Returns", y = "Density", title = "Density Plot of Daily Returns with VaR") +
     theme_minimal()
   # Add shaded area below the VaR to the density plot
   density_data <- ggplot_build(density_plot)$data[[1]] %>%
     as.data.frame() %>%
     filter(x < VaR)
   density_plot +
     geom_ribbon(data = density_data, aes(x = x, ymin = 0, ymax = y), fill = "red", alpha = 0.5)
   

   在这里，我们首先将每日回报率转换为 DataFrame，该 DataFrame 由ggplot()使用geom_density()绘制密度曲线。我们还添加了一条表示 VaR 的垂直线，使用geom_vline()。为了阴影 VaR 下方的区域，我们使用ggplot_build()来过滤数据，并使用geom_ribbon()来着色满足过滤条件的区域。

   运行此代码将生成图 10.17：

图 10.17 – 可视化每日回报率和 VaR 区域

因此，我们可以根据观察到的每日回报率的经验分布来量化 VaR。

摘要

在本章中，我们介绍了常见的概率分布。我们首先介绍了离散概率分布，包括伯努利分布、二项分布、泊松分布和几何分布。然后，我们介绍了常见的连续概率分布，包括正态分布、指数分布和均匀分布。接下来，我们介绍了常见的抽样分布及其在统计推断中对总体统计的应用。最后，我们介绍了顺序统计及其在每日股票回报率 VaR 计算中的应用。

在下一章中，我们将介绍统计估计程序，包括点估计、中心极限定理和置信区间。

第十一章：统计估计

在本章中，我们将向您介绍一系列统计技术，这些技术使您能够使用数值和分类数据来进行推断和估计。我们将探讨关键概念和方法，如假设检验、置信区间和估计技术，这些技术使我们能够从给定的样本中对总体进行概括。

到本章结束时，您将掌握统计推断的核心概念，并能够在不同场景下进行假设检验。

本章将涵盖以下主要内容：

• 分类数据的统计推断

• 数值数据的统计推断

• 构建自助法置信区间

• 介绍 t 分布中使用的中心极限定理

• 使用 t 分布构建总体均值的置信区间

• 对两个均值进行假设检验

• 介绍方差分析（ANOVA）

要运行本章中的代码，您需要以下软件包的最新版本：

• dplyr，1.0.10

• ggplot2，3.4.0

• socviz，1.2

• infer，1.0.4

请注意，前面提到的版本是在我撰写本书时的最新版本。本章的所有代码和数据均可在github.com/PacktPublishing/The-Statistics-and-Machine-Learning-with-R-Workshop/blob/main/Chapter_11/working.R找到。

分类数据的统计推断

一个分类变量具有不同的类别或水平，而不是数值。分类数据在我们的日常生活中很常见，例如性别（男性或女性，尽管现代观点可能不同）、房产销售类型（新房或二手房）和行业。因此，对这些变量进行合理推断的能力对于在多种情境下得出有意义的结论和做出明智的决策至关重要。

作为分类变量，通常意味着我们无法将其传递给模型中的字符串值（如 "finance" 或 "technology"），一种常见的方法是将变量一热编码成多个列，每列对应一个特定行业，表示二进制值 0 或 1。

在本节中，我们将探讨专门用于处理分类数据的各种统计技术，使我们能够根据可用的样本得出有价值的见解并对总体进行推断。我们还将讨论重要概念，如比例、独立性和拟合优度，这些概念是理解和处理分类变量的基础，包括单参数和多参数的情况。

让我们从单个参数的推断开始讨论。

单个参数的统计推断

一个总体参数，我们感兴趣并要推断的主题，是一个描述总体特定统计属性的固定数量，包括均值、比例或标准差。这个数量通常对我们来说是隐藏的。例如，为了得到大学中最受欢迎的专业，我们需要计算整个大学每个专业注册的学生数量，然后返回计数最大的专业。

在单参数统计推断的背景下，我们的目标是估计这个未知参数或基于从样本收集的信息来测试其值的相关假设。换句话说，我们会使用统计推断工具，根据已知的样本来推断未知的人口参数。在前一个例子中，我们会通过特定学术年度注册的学生有限样本来推断整个大学最受欢迎的专业。

让我们先探索 一般社会调查（GSS）数据集。

介绍一般社会调查数据集

GSS 是一个综合数据集，被研究人员和政策制定者广泛用于理解美国的社会、文化和政治趋势。自 1972 年以来，GSS 由芝加哥大学的 国家民意研究中心（NORC）持续进行，目的是收集关于广泛主题的数据，包括态度、行为和对各种问题的观点。

让我们从 socviz 包中加载 GSS 数据集（请记住通过 install.packages("socviz") 安装此包）：

library(socviz)
data(gss_lon)

GSS 数据集现在存储在 gss_lon 变量中，包含总共 62,466 行和 25 列，如下所示：

>>> dim(gss_lon)
62466    25

GSS 数据集包含许多变量，涵盖了各种主题，如教育、收入、家庭结构、政治信仰和宗教归属。让我们使用 dplyr 包中的 glimpse() 函数来检查数据集的结构，该函数旨在帮助您快速探索和理解数据结构：

>>> glimpse(gss_lon)

图 11.1 显示了返回的前几个变量的截图。

图 11.1 – 展示运行 glimpse() 函数的结果的前几行

接下来，我们将查看基于分类变量的特定统计量的计算。

计算样本比例

数据集中的 siblings 列是一个分类变量，追踪家庭中的兄弟姐妹数量。在接下来的练习中，我们希望计算在最新的一年，即 2016 年，家庭有两个兄弟姐妹的受访者比例。

练习 11.1 – 计算兄弟姐妹的样本比例

在这个练习中，我们首先获取siblings列的摘要，并对数据集进行子集化以关注 2016 年，然后将其用于计算家庭中有特定数量兄弟姐妹的调查比例：

1. 使用summary()函数获取siblings列的摘要：

   
   >>> summary(gss_lon$siblings)
       0     1     2     3     4     5    6+  NA's
    3047 10152 11313  9561  7024  5066 14612  1691
   

   结果表明，大多数调查是在有六个或更多兄弟姐妹的家庭中进行的！

2. 对 2016 年的数据集进行子集化：

   
   gss2016 = gss_lon %>% filter(year == 2016)
   Plot the count of siblings in a bar chart using ggplot().
   ggplot(gss2016, aes(x = siblings)) +
     geom_bar() +
     labs(title = "Frequency count of siblings", x = "Number of siblings", y = "Count") +
     theme(text = element_text(size = 16))
   

   运行代码生成图11.2。

图 11.2 – 以柱状图可视化兄弟姐妹数量的频率计数

1. 计算有两个兄弟姐妹的调查比例：

   
   p_hat = gss2016 %>%
     summarize(prop_2sib = mean(siblings=="2", na.rm=TRUE)) %>%
     pull()
   >>> p_hat
   0.208246
   

   在这里，我们使用summarize()函数计算一系列二元值的平均值，这对应于有两个兄弟姐妹的调查比例。然后我们使用pull()函数从结果 DataFrame 中获取比例。

我们使用样本比例来估计总体统计量。换句话说，我们根据可用的样本数据计算有两个兄弟姐妹的家庭比例，以近似如果我们基于总体中的所有数据计算相同的比例时对应的比例。这样的估计伴随着一个置信区间，该区间量化了总体比例的可能值列表。

下一节将展示如何计算样本比例的置信区间。

计算置信区间

置信区间是基于样本数据对总体参数进行推断的重要工具。置信区间提供了一个估计范围，其中包含一个总体参数，如比例，在指定的置信水平（如 95%）下可能被找到。当处理样本比例时，计算置信区间使我们能够更好地理解总体中的真实比例，并评估与总体比例估计相关的不确定性。

我们可以使用以下步骤来计算置信区间：

1. 计算样本比例，ˆp（读作 p-hat）。这是我们基于 2016 年的样本数据计算出的值。在其他情况下，ˆp 是通过将成功的数量（对于感兴趣的属性）除以总样本量来计算的。

2. 确定所需的置信水平，通常表示为(1 − α) x 100%，其中α代表显著性水平。换句话说，这是在原假设为真时拒绝零假设的概率。最常用的置信水平是 90%，95%和 99%。

3. 计算样本比例的标准误差，其公式如下：

SE = √(ˆp(1 − ˆp) / n)

在这里，标准误差也对应于样本比例的标准差，假设它遵循成功概率为ˆp 的伯努利分布（回想一下前一章中伯努利分布的介绍）。这种计算依赖于两个假设：样本中的观测值是独立的，样本中有足够的观测值。检查第二个假设的一个常见经验法则是确保 nˆp > 10 和 n(1 − ˆp) > 10。

或者，我们不必假设伯努利分布，可以使用自举过程来估计标准误差，而不做任何分布假设。自举是一种非参数方法，涉及有放回地重新抽样数据以创建新的样本，对每个重新抽样的数据集计算感兴趣的统计量（在这种情况下，是比例），并从计算统计量的变异性中估计标准误差。

4. 找到对应于预设置信水平的临界值（z 分数）。这可以通过qnorm()函数完成，它给出了标准正态分布的分位数。

5. 计算误差范围（ME）为标准误差与临界值的乘积：

ME = SE * z _ score

6. 通过从样本比例中加减 ME 来计算置信区间，得到以下结果：

下限 = ˆp − ME

上限 = ˆp + ME

置信区间根据特定的置信水平提供了总体比例的可能值列表。参见图 11.3对计算过程的总结。

图 11.3 – 基于样本比例计算置信区间的过程总结

让我们再深入探讨一下自举过程。不假设任何特定分布，自举过程是一种灵活的方法，可以提供更准确的标准误差估计，特别是对于小样本量或数据表现不佳的情况。然而，它可能计算量很大，尤其是在大型数据集或生成许多自举复制时。

图 11.4 提供了自举过程的示意图。让我们回顾一下：

1. 首先，我们从整个数据集开始，并指定感兴趣的变量，在这个例子中是兄弟姐妹变量。这是通过specify()函数实现的。

2. 接下来，我们从变量中抽取样本，进行有放回抽样，其中新样本的大小与原始数据集相同。这种重抽样引入了随机性到结果数据集中。

3. 我们重复这个过程多次，使用generate()函数得到一系列通过自举得到的伪数据集。

4. 对于每个复制的数据集，我们将计算感兴趣的样本统计量，在这种情况下是具有两个兄弟姐妹的观察值的比例。这是通过calculate()函数完成的。

5. 这些通过重复抽样原始数据集得到的样本统计量将形成一个分布，称为自助法分布（通过ggplot()绘制），其标准差（通过summarize()函数提取）将是对标准误的良好近似。

图 11.4 – 使用自助法程序获取标准误的示意图

通过自助法得到的样本传达了样本统计量不同水平的不确定性，并共同形成多个人工样本统计量的密度分布。自助法分布的标准差随后给出样本统计量的标准误。请注意，如specify()、generate()和calculate()等函数均来自 R 中的infer包。在继续以下代码之前，请记住安装此包。

让我们通过以下练习来了解计算置信区间的自助法程序。

练习 11.2 – 通过自助法计算置信区间

在这个练习中，我们将探索计算样本比例的置信区间。置信区间包括一系列估计值，在这些估计值中，给定观察到的样本，真实总体比例可能取值。这是基于实际观察来量化估计总体比例不确定性的方法。除了使用自助法逐步计算过程外，我们还将比较使用假设的伯努利分布的替代方法的结果：

1. 使用前面描述的infer包中的指定-生成-计算程序构建一组自助法样本统计量。请记住构建一个二元变量来指示具有两个兄弟姐妹的观察值的二元条件：

   
   library(infer)
   gss2016 = gss2016 %>%
     mutate(siblings_two_ind = if_else(siblings=="2","Y","N")) %>%
     filter(!is.na(siblings_two_ind))
   bs = gss2016 %>%
     specify(response = siblings_two_ind,
             success = "Y") %>%
     generate(reps = 500,
              type = "bootstrap") %>%
     calculate(stat = "prop")
   

   在这里，我们首先使用if_else()函数创建一个二元指标变量，以表示当前调查中的家庭是否有两个兄弟姐妹。我们还删除了此列中具有NA值的行。接下来，我们使用specify()函数指出感兴趣的siblings_two_ind变量及其对应的成功级别。然后，我们使用generate()函数生成500个自助法样本，并使用calculate()函数通过设置stat = "prop"来获取每个自助法样本中相应的样本统计量（成功比例）。

   让我们观察自助法样本统计量中的内容：

   >>> bs
   Response: siblings_two_ind (factor)
   # A tibble: 500 × 2
      replicate  stat
          <int> <dbl>
    1         1 0.205
    2         2 0.209
    3         3 0.218
    4         4 0.189
    5         5 0.207
    6         6 0.205
    7         7 0.221
    8         8 0.214
    9         9 0.212
   10        10 0.212
   # … with 490 more rows
   # i Use `print(n = ...)` to see more rows
   

   结果显示，bs对象是一个包含 500 行（对应于 bootstrap 样本的总数）和 2 列的tibble DataFrame。第一列（replicate）表示 bootstrap 样本的数量，第二列（stat）表示 bootstrap 样本中成功的比例（即siblings_two_ind==2的行数除以总行数）。

2. 在密度图中绘制 bootstrap 样本统计量：

   
   >>> ggplot(bs, aes(x = stat)) +
     geom_density() +
     labs(title = "Density plot of the sample proportions", x = "Sample proportion", y = "Density") +
     theme(text = element_text(size = 16))
   

运行代码生成图 11.5的图表。这个分布的分散程度，与标准差相关，直接决定了标准误差的大小。此外，如果我们增加 bootstrap 样本的数量，我们预期密度曲线会更加平滑。

图 11.5 – 可视化所有 bootstrap 样本比例的密度图

1. 将标准误差计算为基于 bootstrap 样本比例的实证分布的标准差：

   
   SE = bs %>%
     summarise(sd(stat)) %>%
     pull()
   >>> SE
   0.007181953
   

   在这里，我们使用sd()函数计算bs中stat列的标准差，然后通过pull()函数返回该值。然后，标准误差将根据预定的 z 分数进行缩放，并从原始样本比例（p_hat）中减去和加上，以获得置信区间。

2. 使用 95%置信区间计算原始样本比例的置信区间：

   
   >>> c(p_hat - 2*SE, p_hat + 2*SE)
   0.1938821 0.2226099
   

   在这里，由于 95%置信水平对应于 z 分数 2，我们在从原始样本比例（p_hat）中减去和加上之前，将其与标准误差相乘以获得置信区间。

3. 假设成功的概率服从伯努利分布，使用结构信息计算置信区间：

   
   SE2 = sqrt(p_hat*(1-p_hat)/nrow(gss2016))
   >>> c(p_hat - 2*SE2, p_hat + 2*SE2)
   0.193079 0.223413
   

   在这里，我们使用伯努利分布的方差显式形式来计算标准误差。结果显示的置信区间与使用 bootstrap 方法获得的置信区间相当相似。

置信区间为我们使用观察到的样本比例估计未知总体比例的不确定性提供了一个不确定性度量。让我们看看如何在下一节中解释置信区间。

解释样本比例的置信区间

解释样本比例的置信区间涉及理解区间的含义和相关的置信水平。在我们之前的例子中，bootstrap 方法报告的置信区间为[0.1938821, 0.2226099]。对于这个置信区间有两个层面的解释。

首先，置信区间的范围表明，具有两个兄弟姐妹的家庭的真实比例很可能在 19.39%到 22.26%之间。这个范围基于样本数据，并估计了真实比例的不确定性。

其次，95%的置信区间意味着如果我们多次进行调查（无论是 2016 年还是其他年份），我们会生成不同大小的随机样本，基于这些样本我们可以计算每个样本的 95%置信区间。在这些人工样本中，我们将获得一系列区间，其中大约 95%的区间将包含真实的总体比例。

注意，置信区间仍然是一个估计值，真实的总体比例可能落在计算出的区间之外。然而，置信区间提供了一种有用的方法来量化估计的不确定性，并基于观察到的样本给出真实总体比例的可能值的列表。

下一节将介绍对样本比例的假设检验。

对样本比例进行假设检验

对样本比例的假设检验与前面章节中引入的置信区间密切相关，它捕捉了基于总体数据对未知比例估计的不确定性水平。自然地，观察值较少的样本会导致置信区间较宽。对样本比例的假设检验旨在确定样本中是否有足够的证据来支持或拒绝关于总体比例的声明。这个过程从零假设（H0）开始，它代表了关于总体比例的基本假设。相应地，存在一个备择假设（H1），它代表了我们对零假设进行测试的声明或陈述。然后，假设检验将观察到的样本比例与指定的零假设进行比较，以评估我们是否有足够的证据来拒绝零假设，并支持备择假设。

让我们概述一下进行假设检验所涉及的程序：

1. 制定假设。在这个步骤中，我们设定零假设（H0）和备择假设（H1）。零假设通常表示没有效果，情况保持现状，如 H0 中的等号所示。另一方面，备择假设表示存在效果或差异，如 H1 中的不等号所示。例如，我们可以为 H0 和 H1 设定以下假设：

   • H0: p = p 0（总体比例等于指定的值，p 0）

   • H1: p ≠ p 0（总体比例不等于 p 0）

2. 选择显著性水平（𝜶）。显著性水平是我们用来在原假设为真时拒绝原假设的概率阈值。广泛使用的显著性水平包括 0.05（5%）和 0.01（1%）。

3. 计算检验统计量。现在我们根据实际数据观察到一个样本比例，我们可以计算在零假设为真的情况下观察到这样一个样本比例的概率。这始于使用以下公式计算样本比例的检验统计量（z 分数）：

z = ˆp − p0 * √___________ * p0(1 − p0) / n

其中ˆp 是样本比例，p0 是在零假设下的总体比例，n 是样本大小。有两点需要注意。首先，分母类似于前面提到的基于样本比例的标准差。确实，我们假设一个伯努利分布，成功概率为 p0。在总共 n 个观察值中，样本比例变量的标准差是√___________ * p0(1 − p0) / n。其次，整个项对应于将一个数值转换为特定分布的 z 分数的过程，这是前一章讨论的主题。在这里，我们假设一个均值为 p0，标准差为√___________ * p0(1 − p0) / n 的正态分布。然后我们可以将观察到的样本比例ˆp 转换为相应的 z 分数，以便于后续计算。

注意，我们还可以使用自助法在零假设下计算经验 p 值。

4. 确定 p 值。z 分数是一个落在标准高斯分布上的度量。它是一个检验统计量，我们通常对观察到的检验统计量在此或更极端值上的概率感兴趣。这被称为 p 值，记为ˆp，当我们假设零假设为真时。换句话说，我们试图评估在零假设为真的情况下观察某些现象的可能性。如果观察到ˆp 或更极端数值的概率非常小，我们有信心零假设是错误的，我们可以拒绝 H0 并支持 H1。

注意，对于双尾检验，我们还可以使用标准正态分布并加倍单侧概率来计算 p 值：

p-value = 2P(Z > |z|)

5. 做出决定。最后，我们将 p 值与预设的显著性水平（α）进行比较，并使用以下规则做出决定。

如果 p 值 ≤ α，则拒绝零假设，支持备择假设。这样做表明有足够的证据表明总体比例与假设比例 p0 不同。

如果 p 值 > α，则不能拒绝零假设。这意味着没有足够的证据表明总体比例与 p0 不同。

进行假设检验的过程与上述类似。唯一的区别是使用hypothesise()函数（位于specify()之后），它作为零假设。然后我们执行相同的自举过程，以获得自举样本比例的密度图，随后计算至少与零假设中指出的比例一样极端的总体比例的概率。

让我们通过一个练习来回顾对样本比例进行假设检验的过程。

练习 11.3 – 对样本比例进行假设检验

在这个练习中，我们将设定一个零假设中的假设总体比例，并根据观察到的样本比例来检验这个假设的有效性：

1. 在条形图中绘制 2016 年有和无双胞胎的家庭的频率计数：

   
   gss2016 %>%
     ggplot(aes(x = siblings_two_ind)) +
     geom_bar() +
     labs(title = "Frequency count of families with two siblings", x = "Have two siblings", y = "Count") +
     theme(text = element_text(size = 16))
   

   运行代码生成图 11.6，该图显示有两个双胞胎的家庭占所有家庭的约 1/4。

图 11.6 – 可视化有两个双胞胎的家庭的频率计数

1. 计算有两个双胞胎的家庭的样本比例：

   
   p_hat = gss2016 %>%
     summarize(mean(siblings_two_ind=="Y")) %>%
     pull()
   >>> p_hat
   0.208246
   

   这里，我们首先使用siblings_two_ind=="Y"构建一系列二元结果。取该列的平均值给出TRUE值的比率，这在summarize()上下文中执行。然后我们使用pull()提取样本比例的值。

2. 使用指定-假设-生成-计算程序在零假设下生成一系列自举样本比例，该假设指定总体比例为0.19：

   
   null = gss2016 %>%
     specify(response = siblings_two_ind,
             success = "Y") %>%
     hypothesise(null = "point",
                 p = 0.19) %>%
     generate(reps = 500,
              type = "draw") %>%
     calculate(stat = "prop")
   >>> null
   Response: siblings_two_ind (factor)
   Null Hypothesis: point
   # A tibble: 500 × 2
      replicate  stat
      <fct>     <dbl>
    1 1         0.179
    2 2         0.193
    3 3         0.176
    4 4         0.181
    5 5         0.181
    6 6         0.198
    7 7         0.191
    8 8         0.189
    9 9         0.194
   10 10        0.189
   # … with 490 more rows
   # i Use `print(n = ...)` to see more rows
   

3. 通过垂直线生成自举样本比例的密度图以及零假设建议的比例：

   
   ggplot(null, aes(x = stat)) +
     geom_density() +
     geom_vline(xintercept = p_hat,
                color = "red") +
     labs(title = "Density plot using bootstrap", x = "Sample proportion", y = "Density") +
     theme(text = element_text(size = 16))
   

   运行代码生成图 11.7，根据零假设，红色线（根据零假设）所示值的观察概率是密度曲线右侧的总面积。然后我们将结果加倍，以考虑相反的方向。

图 11.7 – 可视化假设检验中自举样本比例的密度图

1. 计算 p 值：

   
   >>> null %>%
     summarise(mean(stat > p_hat)) %>%
     pull()* 2
   0.02
   

   由于这个结果小于预设的显著性水平 5%，我们有足够的证据支持备择假设并拒绝原假设。换句话说，假设的 19%与真实总体比例在 95%的置信水平下存在统计学上的差异。因此，我们可以得出结论，真实总体比例不是 19%。

下一个部分将探讨两个分类变量之间样本比例差异的推断。

样本比例差异的推断

现在的设置是我们有两个分类变量。以性别和学位为例。数据将报告女性和男性的学位持有者的比例。一个自然的问题是要问男性是否比女性更有可能获得学位。特定的数据集将报告这些比例的快照，这可能或可能不会表明学位持有者的百分比更高。假设检验的工具可以用来回答以下问题：如果数据集中男性是学位持有者的更高比例，这种差异在统计上是否显著？换句话说，男性是否比女性更有可能获得学位，反之亦然？本节试图回答这类问题。

对两个分类变量（例如，性别和学位）之间样本比例差异的推断涉及比较两个不同群体中每个级别的样本比例。这种分析通常用于实验或观察性研究，以确定两组之间比例差异的存在性。主要目标是估计群体比例之间的差异，并确定这种差异在统计上是否显著。

假设检验的程序与之前类似。我们首先提出零假设，即假设两个群体的比例之间没有差异，即 p1 = p2，或者 p1 - p2 = 0。备择假设随后表明它们的差异不是零；即 p1 - p2 ≠ 0。接下来，我们选择一个特定的显著性水平，计算样本统计量（样本比例的差异，包括两个分类变量之间的合并比例）和检验统计量（基于假设分布的闭式表达式或使用自助法）。最后，我们获得 p 值并决定在零假设下观察到的结果是否具有统计显著性。

让我们通过之前的例子进行一个具体的练习。

练习 11.4 – 对样本比例差异进行假设检验

在这个练习中，我们关注如何对性别和更高学位状态之间的样本比例差异进行假设检验。在这里，我们将更高学位定义为学士学位及以上。男性和女性群体中更高学位持有者的比例可能会不同，我们将测试在观察到的数据下这种差异是否显著：

1. 在之前的 DataFrame gss2016中添加一个名为higher_degree的二进制列，以指示更高学位的状态，包括学士及以上：

   
   gss2016 = gss2016 %>%
     mutate(higher_degree = if_else(degree %in% c("Bachelor","Graduate"), "Y", "N"))
   

2. 打印gender和higher_degree两个级别的比率：

   
   >>> table(gss2016$higher_degree)
      N    Y
   2008  854
   >>> table(gss2016$sex)
     Male Female
     1274   1588
   

3. 将这些计数在柱状图中绘制：

   
   ggplot(gss2016, aes(x = sex, fill=higher_degree)) +
     geom_bar() +
     labs(title = "Frequency count for gender and degree", x = "Gender", y = "Count") +
     theme(text = element_text(size = 16))
   

   运行代码生成图 11.8。

图 11.8 – 可视化性别和更高学位状态的频率计数

我们也可以通过在geom_bar()函数中指定position = "fill"来按百分比绘制它们：

ggplot(gss2016, aes(x = sex, fill=higher_degree)) +
  geom_bar(position = "fill") +
  labs(title = "Sample proportions for gender and degree", x = "Gender", y = "Ratio") +
  theme(text = element_text(size = 16))

运行代码生成图 11.9，该图表明男性和女性群体之间更高学位持有者比例没有明显差异。

图 11.9 – 可视化性别和更高学位状态的频率计数

1. 计算男性和女性之间更高学位持有者样本比例的差异：

   
   p_hats = gss2016 %>%
     group_by(sex) %>%
     summarise(mean(higher_degree=="Y", na.rm=TRUE)) %>%
     pull()
   d_hat = diff(p_hats)
   >>> d_hat
   0.007288771
   

   结果还显示，差异相当小，女性群体比男性群体高 0.7%（参考前一个图中女性群体略高的蓝色条形）。让我们看看这种差异是否具有统计学意义。

2. 在零假设下生成一个自助样本集，该假设表明男性和女性群体之间更高学位持有者比例没有差异：

   
   gss2016 %>%
     specify(
       response = higher_degree,
       explanatory = sex,
       success = "Y"
     ) %>%
     hypothesise(null = "independence") %>%
     generate(reps = 1, type = "permute")
   Response: higher_degree (factor)
   Explanatory: sex (factor)
   Null Hypothesis: independence
   # A tibble: 2,862 × 3
   # Groups:   replicate [1]
      higher_degree sex    replicate
      <fct>         <fct>      <int>
    1 N             Male           1
    2 N             Male           1
    3 Y             Male           1
    4 N             Female         1
    5 N             Female         1
    6 N             Female         1
    7 Y             Male           1
    8 N             Female         1
    9 N             Male           1
   10 N             Male           1
   # … with 2,852 more rows
   # i Use `print(n = ...)` to see more rows
   

   在这里，我们将higher_degree作为响应变量，将sex作为解释变量，在逻辑回归设置中使用（将在第十三章中介绍）。在零假设下，我们从原始数据集中随机抽样并创建一个具有相同形状的新的人工数据集。

3. 重复相同的自助抽样过程 500 次，并计算每对自助样本中女性和男性群体之间更高学位持有者样本比例的差异（注意这里的顺序）：

   
   null = gss2016 %>%
     specify(
       higher_degree ~ sex,
       success = "Y"
     ) %>%
     hypothesise(null = "independence") %>%
     generate(reps = 500, type = "permute") %>%
     calculate(stat = "diff in props", order = c("Female", "Male"))
   >>> null
   Response: higher_degree (factor)
   Explanatory: sex (factor)
   Null Hypothesis: independence
   # A tibble: 500 × 2
      replicate     stat
          <int>    <dbl>
    1         1  0.00870
    2         2  0.00587
    3         3 -0.00120
    4         4  0.0228
    5         5  0.00446
    6         6 -0.00827
    7         7 -0.0366
    8         8  0.0129
    9         9  0.0172
   10        10 -0.00261
   # … with 490 more rows
   # i Use `print(n = ...)` to see more rows
   

4. 将这些自助样本统计量绘制在密度曲线上，并将观察到的差异作为垂直红色线：

   
   ggplot(null, aes(x = stat)) +
     geom_density() +
     geom_vline(xintercept = d_hat, color = "red") +
     labs(x = "Difference in sample proportion (female - male)", y = "Count") +
     theme(text = element_text(size = 16))
   

   运行代码生成图 11.10，该图显示红色线并不位于经验分布的极端一侧。这表明接下来要计算的 p 值可能较高。

图 11.10 – 显示自助样本统计量和观察到的差异的密度图

1. 计算双尾 p 值：

   
   null %>%
     summarize(pval = 2 * mean(stat > d_hat)) %>%
     pull()
   0.608
   

   结果显示一个相当高的 p 值，这表明我们没有足够的证据来拒绝零假设。换句话说，没有足够的信息表明男性和女性之间更高学位持有者的比例不同。

假设检验依赖于一个预定义的显著性水平。这个显著性水平，用α表示，与该过程的统计误差有关。下一节介绍了在执行假设检验时常见的两种统计误差类型。

一类错误和二类错误

在进行假设检验并对零假设（H0）和备择假设（H1）做出决策时，存在两种类型的错误。它们被称为第一类错误和第二类错误。

第一类错误指的是假阳性。它发生在零假设为真但错误地被拒绝时。换句话说，我们在样本数据中找到了表明存在显著效应或差异的证据，我们倾向于备择假设，尽管它实际上在总体中并不存在。我们用α表示经历第一类错误的概率。它也被称为显著性水平，在先前的例子中设置为0.05。5%的显著性水平意味着当零假设为真时，有 5%的概率拒绝零假设。因此，显著性水平代表犯假阳性错误的概率。

第二类错误专注于假阴性情况。它发生在我们未能拒绝一个错误的零假设时。换句话说，即使在总体中存在，我们在样本数据中也没有找到拒绝零假设的证据。犯第二类错误的概率用β表示，也被称为测试的效力。效力的补数，表示为 1 − β，代表在假设错误时拒绝零假设的概率。

第一类错误涉及错误地拒绝零假设，而第二类错误涉及在假设错误时未能拒绝零假设。在这两种类型的错误都是假设检验中的重要考虑因素，因为它们可能导致错误的结论。为了最大限度地减少这些错误的风险，我们可以谨慎地选择显著性水平（α），并确保他们的研究具有足够的效力（1 − β）。测试的效力取决于样本大小、效应量（这是一个变量之间经验关系的量度）和选择的显著性水平。较大的样本大小和较大的效应量都会增加测试的效力，从而降低第二类错误的可能性。

图 11.11 提供了假设检验中不同类型结果的概述。请注意，假阳性和假阴性与决策的质量有关。根据错误决策的类型，我们将错误分类为第一类或第二类错误。

图 11.11 – 假设检验中不同类型结果概述

下一节介绍了卡方检验，该检验用于检验两个分类变量的独立性。

检验两个分类变量的独立性

要检查两个分类变量的独立性，涉及检查它们之间是否存在统计上显著的关系。一个常见的程序是独立性卡方检验。它通过比较列联表中的观察频率与独立性假设下的预期频率来工作。

让我们先回顾一下两个分类变量的列联表。

介绍列联表

列联表，也称为交叉表或 crosstab，是一种用于显示两个或更多分类变量频率分布的表格。它通过显示变量的类别如何在数据中相交或同时出现来总结变量之间的关系。它提供了分类变量之间关系的良好总结。

让我们继续以性别与学位之间的关系为例。这次，我们将查看所有类型的学位，如下面的代码所示：

>>> table(gss2016$degree)
Lt High School    High School Junior College       Bachelor        Graduate
           328           1459            215            536             318

为了表示其与性别的关联，我们可以像之前一样，将学位与性别一起绘制在堆积条形图中：

ggplot(gss2016, aes(x = sex, fill=degree)) +
  geom_bar() +
  labs(title = "Frequency count for gender and degree", x = "Gender", y = "Count") +
  theme(text = element_text(size = 16))

运行代码生成图图 11.12。

图 11.12 – 在条形图中可视化性别与学位之间的关系

然而，该图没有提供每个类别的确切计数信息。为了获得两个变量每个类别的确切频率，我们可以使用列联表：

tab = gss2016 %>%
  select(sex, degree) %>%
  table()
>>> tab
        degree
sex      Lt High School High School Junior College Bachelor Graduate
  Male              147         661             89      243      132
  Female            181         798            126      293      186

在这里，我们使用了table()函数在选择了sex和degree之后生成列联表。

下一节介绍了卡方检验，用于检验这两个分类变量之间的独立性。

应用两个分类变量之间的独立性卡方检验

卡方检验是一种统计检验，用于决定一组观察样本中两个分类变量之间可能存在显著关系（依赖性）。它可以用来检验独立性或拟合优度。在本章中，我们主要关注两个分类变量之间的独立性检验。该检验比较列联表中观察到的频率与预期的频率，假设变量是独立的。如果观察到的频率和预期的频率有显著差异，则检验表明变量不是独立的；换句话说，它们相互依赖。

按照之前的相同方法，我们可以生成一个人工的 bootstrap 数据集以获得一个样本统计量，称为卡方统计量。该数据集是在原数据集的形状下，假设在零假设下的独立性进行排列生成的。在下面的代码中，我们生成一个与原数据集形状相同的一个排列数据集，假设在零假设下的独立性：

perm_1 = gss2016 %>%
  # Specify the variables of interest
  specify(degree ~ sex) %>%
  # Set up the null hypothesis
  hypothesize(null = "independence") %>%
  # Generate a single permuted dataset
  generate(reps = 1, type = "permute")
>>> perm_1
Response: degree (factor)
Explanatory: sex (factor)
Null Hypothesis: independence
# A tibble: 2,856 × 3
# Groups:   replicate [1]
   degree         sex    replicate
   <fct>          <fct>      <int>
 1 Junior College Male           1
 2 Bachelor       Male           1
 3 High School    Male           1
 4 High School    Female         1
 5 High School    Female         1
 6 High School    Female         1
 7 Graduate       Male           1
 8 Bachelor       Female         1
 9 High School    Male           1
10 High School    Male           1
# … with 2,846 more rows
# i Use `print(n = ...)` to see more rows

接下来，我们创建 500 个排列数据集并提取相应的卡方统计量：

null_spac = gss2016 %>%
  specify(degree ~ sex) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 500, type = "permute") %>%
  calculate(stat = "Chisq")
>>> null_spac
Response: degree (factor)
Explanatory: sex (factor)
Null Hypothesis: independence
# A tibble: 500 × 2
   replicate  stat
       <int> <dbl>
 1         1  3.50
 2         2  1.11
 3         3 14.0
 4         4  4.62
 5         5  1.41
 6         6  1.41
 7         7  9.69
 8         8  4.17
 9         9  5.97
10        10  2.86
# … with 490 more rows
# i Use `print(n = ...)` to see more rows

为了运行测试，我们在假设分类变量之间相互独立的情况下，获取列联表中每个单元格的预期频率。单元格的预期频率计算为（行和 列和） / 总体和：

# calculate expected frequency table
row_totals = rowSums(tab)
col_totals = colSums(tab)
overall_total = sum(tab)
expected = outer(row_totals, col_totals) / overall_total
>>> expected
       Lt High School High School Junior College Bachelor Graduate
Male          146.084    649.8067        95.7563 238.7227 141.6303
Female        181.916    809.1933       119.2437 297.2773 176.3697

在这里，我们首先获取行和列的和以及总和。然后我们使用outer()函数获取这两个向量的外积，然后将其按总和缩放以获得每个单元格的预期频率计数。

现在，我们根据可用的样本计算观察到的卡方统计量：

# Compute chi-square statistic
observed_chi_square = sum((tab - expected)² / expected)
>>> observed_chi_square
2.536349

我们可以在之前通过自助法得到的样本统计量的密度曲线中绘制观察到的卡方统计量，以了解观察到的统计量所在的位置，据此我们可以计算出相应的 p 值：

ggplot(null_spac, aes(x = stat)) +
  geom_density() +
  geom_vline(xintercept = observed_chi_square, color = "red") +
  labs(title = "Density curve of bootstrapped chi-square statistic", x = "Chi-square statistic", y = "Density") +
  theme(text = element_text(size = 16))

运行代码生成了图 11.13，它显示了一个较高的 p 值。

图 11.13 – 可视化自助法卡方统计量的密度曲线和观察到的统计量

现在我们可以计算 p 值。如下面的代码所示，p 值为0.72确实相当高，因此没有足够的证据来拒绝零假设：

>>> null_spac %>%
  summarize(pval = 2 * mean(stat < observed_chi_square)) %>%
  pull()
0.72

在下一节中，我们将转向观察数值数据的统计推断。

数值数据的统计推断

在本节中，我们将转向使用数值数据进行统计推断。我们将介绍两种方法。第一种方法依赖于自助法程序，并排列原始数据集以创建额外的艺术数据集，然后可以使用这些数据集来推导置信区间。第二种方法使用对自助法样本分布的理论假设，并依赖于 t 分布来实现相同的结果。我们将学习如何进行 t 检验、推导置信区间以及进行方差分析（ANOVA）。

如前所述，自助法是一种非参数重采样方法，它允许我们估计特定统计量的抽样分布，例如均值、中位数或比例，就像上一节中那样。这是通过从原始数据中反复有放回地抽取随机样本来实现的。通过这样做，我们可以计算置信区间并执行假设检验，而不依赖于特定的分布假设。

此外，t 分布是一种概率分布，用于在样本量小且总体数据的标准差未知的情况下进行假设检验。它是一种更通用的方法，假设通过自助法得到的样本遵循特定的分布。然后我们将使用这个分布来估计置信区间并执行假设检验。

t 检验是一种广泛使用的统计检验方法，它允许我们比较两组数据的平均值，或者检验单个组数据的平均值是否等于一个特定的值。这次，我们的兴趣在于组数据的平均值，因为变量是数值型的。该检验依赖于 t 分布，并考虑样本大小、样本平均值和样本方差。

置信区间提供了一组可能的值，其中真实的总体统计量（如均值或比例）很可能位于，并且具有指定的置信水平（由显著性水平α指定）。

最后，ANOVA 在比较多于两组时扩展了 t 检验的应用。ANOVA 帮助我们通过将观察数据的总变异性分为两部分：组间变异性与组内变异性，来确定组均值之间可能存在的显著差异。它检验了所有组均值相等的零假设。如果零假设被拒绝，我们可以继续识别哪些特定的组均值彼此不同。

让我们从生成中位数自助分布开始。

生成中位数自助分布

如前所述，当构建单个统计量的自助分布时，我们首先通过有放回抽样生成一系列自助样本，然后记录每个分布的相关统计量（在本例中为中位数）。

让我们通过一个练习来构建自助样本集合。

练习 11.5 – 生成样本中位数的自助分布

在这个练习中，我们将使用infer包应用相同的指定-生成-计算工作流程，使用mtcars数据集生成样本中位数的自助分布。

加载mtcars数据集并查看其结构：

data(mtcars)
>>> str(mtcars)
'data.frame':  32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

结果显示，我们有一个包含 32 行和 11 列的数据集。在接下来的步骤中，我们将使用mpg变量并生成其中位数的自助分布：

1. 根据变量mpg生成 10,000 个自助样本并获取所有样本的中位数：

   
   bs <- mtcars %>%
     specify(response = mpg) %>%
     generate(reps = 10000, type = "bootstrap") %>%
     calculate(stat = "median")
   >>> bs
   Response: mpg (numeric)
   # A tibble: 10,000 × 2
      replicate  stat
          <int> <dbl>
    1         1  21.4
    2         2  22.2
    3         3  20.4
    4         4  17.8
    5         5  19.2
    6         6  19.2
    7         7  18.4
    8         8  20.4
    9         9  19.0
   10        10  21.4
   # … with 9,990 more rows
   # i Use `print(n = ...)` to see more rows
   

在这里，我们在calculate()函数中指定stat = "median"以从每个自助样本中提取中位数。

2. 将自助分布作为自助样本统计量的密度曲线进行绘图：

ggplot(bs, aes(x = stat)) +
  geom_density() +
  labs(title = "Density plot for bootstrapped median", x = "Median", y = "Probability") +
  theme(text = element_text(size = 16))

运行代码生成图 11.14.14 的图表。

图 11.14 – 可视化自助法样本中位数的密度曲线

下一节将探讨构建自举置信区间。

构建自举置信区间

我们已经探讨了如何使用标准误差方法构建自举置信区间。这涉及到从观察到的样本统计量中添加和减去缩放后的标准误差。结果发现，还有一种更简单的方法，它仅使用自举分布的百分位数来获得置信区间。

让我们继续之前的例子。假设我们想要计算之前自举分布的 95%置信区间。我们可以通过计算自举分布的上限和下限分位数（分别为 97.5%和 2.5%）来实现这一点。下面的代码实现了这一点：

>>> bs %>%
  summarize(
    l = quantile(stat, 0.025),
    u = quantile(stat, 0.975)
  )
# A tibble: 1 × 2
      l     u
  <dbl> <dbl>
1  16.6  21.4

让我们也使用标准误差方法计算自举置信区间，如下面的代码所示：

SE = bs %>%
  summarise(sd(stat)) %>%
  pull()
observed_median = median(mtcars$mpg)
>>> c(observed_median - 2*SE, observed_median + 2*SE)
16.64783 21.75217

如预期的那样，结果是使用百分位数方法获得的结果相近。然而，标准误差方法比百分位数方法更准确。

下一节将介绍在测试零假设后重新中心化自举分布。

重新中心化自举分布

之前章节中的自举分布是通过随机有放回抽样原始数据集生成的。每一组自举样本保持与原始样本集相同的大小。然而，我们无法直接使用这个自举分布进行假设检验。

在引入零假设后，我们在之前的假设检验部分为两个分类变量所做的是在零假设下重新生成一个新的自举分布。然后我们将观察到的样本统计量作为一条垂直的红线沿着自举分布放置，以计算 p 值，表示至少与观察到的样本统计量一样极端的现象发生的概率。唯一的额外步骤是生成零假设下的自举分布。

在零假设下生成自举样本时，主要思想是消除我们正在测试的效果并创建样本，假设零假设为真。换句话说，我们创建的样本应该是如果没有组间差异时预期的样本。例如，当比较两组之间的均值时，我们会在进行随机有放回抽样之前，从每个观察值中减去总体均值，以将数据中心在 0 周围。

有另一种实现方式。回想一下，原始的自助分布，按照设计，是围绕观察到的样本统计量中心化的。在引入零假设后，我们可以简单地移动原始自助分布，使其围绕零假设中的统计量中心化，即零值。这个移动后的自助分布代表了如果我们从原始数据集中移除效应并再次进行自助抽样，我们将得到的相同分布。然后我们可以将观察到的样本统计量放置在移动后的自助分布上，以计算相应的 p 值，这代表了生成至少与实际样本统计量一样有利的样本统计量的模拟比例。图 11.15 展示了这一过程。

图 11.15 – 将自助分布移至围绕零值中心化

让我们为前一个示例生成假设检验的自助分布。我们想要测试mpg变量的总体中位数为 16 的零假设。以下代码生成了自助样本统计量，我们在hypothesize()函数中通过med = 16指定零值，并通过null = "point"指定点估计：

bs = mtcars %>%
  specify(response = mpg) %>%
  hypothesize(null = "point", med = 16) %>%
  generate(reps = 10000, type = "bootstrap") %>%
  calculate(stat = "median")
>>> bs
Response: mpg (numeric)
Null Hypothesis: point
# A tibble: 10,000 × 2
   replicate  stat
       <int> <dbl>
 1         1  16
 2         2  16.2
 3         3  18
 4         4  16.2
 5         5  16
 6         6  16.2
 7         7  16
 8         8  16
 9         9  17.8
10        10  16
# … with 9,990 more rows
# i Use `print(n = ...)` to see more rows

现在，我们在密度图中绘制这些自助样本统计量，同时将观察到的样本统计量作为垂直红色线：

ggplot(bs, aes(x = stat)) +
  geom_density() +
  geom_vline(xintercept = median(mtcars$mpg), color = "red") +
  labs(title = "Density curve of bootstrapped median", x = "Sample median", y = "Density") +
  theme(text = element_text(size = 16))

运行代码生成了图 11.16中的图，它显示了一个较小的 p 值。

图 11.16 – 自助样本中位数和观察样本中位数的密度图（垂直红色线）

在下一节中，我们将介绍基于中心极限定理（CLT）的另一种基于分布的推断方法。

介绍在 t 分布中使用的中心极限定理

CLT（中心极限定理）指出，许多独立同分布的随机变量的和（或平均值）的分布将共同形成一个正态分布，无论这些个别变量的潜在分布如何。由于 CLT，正态分布常被用来近似各种统计量的抽样分布，例如样本均值和样本比例。

在统计推断的背景下，t 分布与中心极限定理（CLT）相关。当我们从样本中估计总体均值时，我们通常无法访问总体的真实标准差。相反，我们求助于样本标准差作为估计。在这种情况下，样本均值的抽样分布不遵循正态分布，而是 t 分布。换句话说，当我们从一组观察到的样本中提取样本均值，并且我们不确定总体标准差（在实际数据工作中这种情况很常见），样本均值可以建模为 t 分布的一个实现。

t 分布是一族连续概率分布，它们是对称的、钟形的，这显示了与正态分布的相似性。然而，t 分布的尾部更厚，这解释了由于从观察数据中估计总体标准差而产生的更大不确定性。也就是说，t 分布的观测值更有可能落入远离均值的尾部（例如，远离均值的两个标准差之外），而正态分布则不然。t 分布的形状依赖于自由度（df），这取决于样本大小并决定了尾部的厚度。随着收集到更多样本，df 增加，t 分布逐渐逼近正态分布。

在上一章中，我们简要介绍了用于在 t 分布下找到截断值的qt()函数。现在让我们通过一个练习来更熟悉与 t 分布相关的计算。

练习 11.6 – 理解 t 分布

在这个练习中，我们将使用pt()函数来找到 t 分布下的概率。对于给定的截断分位数值q和给定的df，pt(q, df)函数给出了t小于q的 t 分布下df的概率。换句话说，我们有 P(t df < T) = pt(q = T, df)。我们还可以使用qt()函数来找到 t 分布下特定概率的分位数。也就是说，如果 P(t df < T) = p，那么 T = qt(p, df)：

1. 找到在 t 分布下，自由度为 10 且低于T=3的概率：

   
   x = pt(3, df = 10)
   >>> x
   0.9933282
   

2. 找到在 t 分布下，自由度为 10 且T=3以上的概率：

   
   y = 1 - x
   >>> y
   0.006671828
   

   注意，我们首先计算在 t 分布下低于特定截断值的概率，然后取补数以找到高于阈值的概率。

3. 找到在 t 分布下，自由度为100且T=3以上的概率：

   
   z = 1 - pt(3, df = 100)
   >>> z
   0.001703958
   

   由于df=100比df=10更接近正态分布，因此得到的概率z因此比y小。

4. 找到具有 10 个自由度的 t 分布的 95 百分位数：

   
   d = qt(0.95, df = 10)
   >>> d
   1.812461
   

5. 找到界定中间 95 百分位数上限的 t 分布的截断值，该分布具有10个自由度：

   
   e = qt(0.975, df = 10)
   >>> e
   2.228139
   

   这里，中间 95 百分位数的上限指的是 97.5 百分位数。

6. 找到用100自由度（df）界定 t 分布中间 95%百分位上限的截断值：

   
   f = qt(0.975, df = 100)
   >>> f
   1.983972
   

下一节将讨论如何使用 t 分布来构建总体均值的置信区间。

构建使用 t 分布的总体均值置信区间

让我们回顾一下对总体均值进行统计推断的过程。我们从一个有限的样本开始，从这个样本中我们可以推导出样本均值。由于我们想要估计总体均值，我们希望基于观察到的样本均值进行统计推断，并量化总体统计量可能存在的范围。

例如，以下代码中显示的平均每加仑英里数，在mtcars数据集中大约是 20。

>>> mean(mtcars$mpg)
20.09062

给出这个结果，我们不会对遇到另一个具有平均mpg为 19 或 21 的类似数据集感到惊讶。然而，如果值是 5、50，甚至 100，我们就会感到惊讶。在评估新的样本集合时，我们需要一种方法来量化样本均值在多个样本中的变异性。我们已经学习了两种方法来做这件事：使用自助法来模拟人工样本或使用中心极限定理来近似这种变异性。在本节中，我们将重点关注中心极限定理方法。

根据中心极限定理，任何抽样分布的样本均值将大致呈正态分布，无论原始分布如何。换句话说，我们有以下情况：

_ x  ∼ N(mean = μ, SE =  σ _ √ _ n  )

注意，这是一个我们无法获得的理论分布。例如，总体标准差σ保持未知，我们只能获得观察到的样本。相反，我们将使用样本标准差 s 来估计标准误差，如下所示：

_ x  ∼ N(mean = μ, SE =  s _ √ _ n  )

然后，我们将使用 n - 1 个自由度的 t 分布来进行对总体均值的推断，因为它由于 s 引入的额外不确定性而具有更厚的尾部。

此外，请注意，使用中心极限定理（CLT）的近似依赖于几个假设。例如，样本之间需要相互独立。当样本是随机选择的时候，或者如果样本在无替换的情况下选取且占总人口的不到 10%，这通常可以得到满足。样本大小也需要更大，以便考虑样本中可能存在的偏度。

我们可以使用t.test()函数来构建 95%置信区间，如下所示：

# Construct 95% CI for avg mpg
>>> t.test(mtcars$mpg)
  One Sample t-test
data:  mtcars$mpg
t = 18.857, df = 31, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 17.91768 22.26357
sample estimates:
mean of x
 20.09062

在这里，我们正在进行单样本 t 检验，其中默认的零假设是总体均值为 0。结果显示一个非常小的 p 值，这表明我们可以拒绝零假设，支持备择假设；也就是说，总体均值不是 0。95%置信区间（在17.91768和22.26357之间）也是基于具有31个 df 和t 统计量为18.857的 t 分布构建的。

下一节回顾了使用自举模拟和 t 检验近似进行两个均值假设检验的过程。

对两个均值进行假设检验

在本节中，我们将探讨使用假设检验比较两个样本均值的过程。在比较两个样本均值时，我们想要确定两个不同总体或组之间是否存在显著的差异。

假设现在我们有两个样本组。这两组样本可以代表每个样本在处理前后的特定值。因此，我们的目标是比较这两组样本的样本统计量，例如样本均值，并确定治疗是否有影响。为此，我们可以通过自举模拟或 t 检验近似来执行假设检验，以比较两个独立分布的均值。

在假设检验中使用 t 检验比较两个独立样本的均值时，双样本 t 检验假设数据呈正态分布，并且两个总体的方差相等。然而，在这些假设可能不成立的情况下，可以采用替代的非参数检验或重采样方法，如自举，来对总体均值进行推断。

让我们通过一个练习来看看这两种假设检验方法是如何应用的。

练习 11.7 – 比较两个均值

在这个练习中，我们将探索两种方法（t 检验和自举）来比较两个样本均值并计算样本均值差异的置信区间：

1. 生成一个由两组样本组成的虚拟数据集：

   
   # Define two samples
   sample1 = c(10, 12, 14, 16, 18)
   sample2 = c(15, 17, 19, 21, 23)
   # Combine samples into a data frame
   data = tibble(
     value = c(sample1, sample2),
     group = factor(rep(c("Group 1", "Group 2"), each = length(sample1)))
   )
   >>> data
   # A tibble: 10 × 2
      value group
      <dbl> <fct>
    1    10 Group 1
    2    12 Group 1
    3    14 Group 1
    4    16 Group 1
    5    18 Group 1
    6    15 Group 2
    7    17 Group 2
    8    19 Group 2
    9    21 Group 2
   10    23 Group 2
   

   在这里，我们创建了一个tibble数据框，其中value列表示样本观测值，group列表示组号。我们希望评估这两组样本之间样本均值的差异。

2. 进行 1000 次自举抽样，并在以下假设下计算自举统计量：这两个组相互独立，它们的均值没有差异：

   
   bootstrap_results = data %>%
     specify(response = value, explanatory = group) %>%
     hypothesize(null = "independence") %>%
     generate(reps = 1000, type = "bootstrap") %>%
     calculate(stat = "diff in means", order = c("Group 1", "Group 2"))
   >>> bootstrap_results
   Response: value (numeric)
   Explanatory: group (factor)
   Null Hypothesis: independence
   # A tibble: 1,000 × 2
      replicate  stat
          <int> <dbl>
    1         1 -7.5
    2         2 -6.17
    3         3 -5
    4         4 -2.20
    5         5 -8.05
    6         6 -4.2
    7         7 -3.5
    8         8 -6.67
    9         9 -4.2
   10        10 -6.90
   # … with 990 more rows
   # i Use `print(n = ...)` to see more rows
   

   在这里，假设检验的管道首先通过指定响应变量（value）和解释变量（group），设置零假设，在零假设下生成自举样本，然后计算每个自举样本的检验统计量（在这种情况下，是均值差异）。零假设表明我们假设两组的样本均值值来自同一总体，并且任何观察到的差异仅仅是偶然的。

3. 根据自举统计量计算置信区间：

   
   ci = bootstrap_results %>%
     filter(!is.na(stat)) %>%
     get_confidence_interval(level = 0.95, type = "percentile")
   >>> ci
   # A tibble: 1 × 2
     lower_ci upper_ci
        <dbl>    <dbl>
   1       -9    -1.17
   

4. 使用t.test()函数进行双样本 t 检验：

   
   t_test_result = t.test(sample1, sample2)
   >>> t_test_result
      Welch Two Sample t-test
   data:  sample1 and sample2
   t = -2.5, df = 8, p-value = 0.03694
   alternative hypothesis: true difference in means is not equal to 0
   95 percent confidence interval:
    -9.6120083 -0.3879917
   sample estimates:
   mean of x mean of y
          14        19
   

   结果显示，基于 t 分布的 95%置信区间接近但仍然与通过自举抽样获得的结果不同。我们还可以通过传递模型形式来进行 t 检验：

   t_test_result2 = t.test(value ~ group, data = data)
   >>> t_test_result2
       Welch Two Sample t-test
   data:  value by group
   t = -2.5, df = 8, p-value = 0.03694
   alternative hypothesis: true difference in means between group Group 1 and group Group 2 is not equal to 0
   95 percent confidence interval:
    -9.6120083 -0.3879917
   sample estimates:
   mean in group Group 1 mean in group Group 2
                      14                    19
   

下一节介绍了方差分析（ANOVA）。

介绍方差分析（ANOVA）

ANOVA是一种用于比较两个以上组均值的统计假设检验方法，它扩展了上一节讨论的双样本 t 检验。ANOVA 的目标是在考虑每个组内的变异性（组内变异性）的同时，测试组均值之间是否存在潜在的显著差异（组间变异性）。

ANOVA 在假设检验中依赖于 F 统计量。F 统计量是两个方差估计值的比率：组间方差和组内方差。组间方差衡量组均值之间的差异，而组内方差表示每个组内的变异性。F 统计量可以根据这两个组方差计算得出。

在假设检验中，方差分析（ANOVA）的零假设指出所有组均值相等，任何观察到的差异都是由于偶然性造成的。另一方面，备择假设则表明至少有一个组的均值与其他组不同。如果 F 统计量足够大，组间方差显著大于组内方差，这为反对零假设提供了证据。

让我们来看一个具体的例子。我们首先加载PlantGrowth数据集，该数据集包含植物经过三种不同处理后重量：

data(PlantGrowth)
>>> str(PlantGrowth)
'data.frame':   30 obs. of  2 variables:
 $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
 $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

接下来，我们使用相同的指定-假设-生成-计算程序执行单因素 ANOVA 测试。具体来说，我们首先指定响应变量（重量）和解释变量（组别）。然后，我们使用hypothesize(null = "独立性")建立零假设，即各组均值没有差异。接下来，我们使用generate(reps = 1000, type = "permute")生成 1,000 个排列数据集。最后，我们使用calculate(stat = "F")计算每个排列数据集的 F 统计量：

anova_results = PlantGrowth %>%
  specify(response = weight, explanatory = group) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "F")
>>> anova_results
Response: weight (numeric)
Explanatory: group (factor)
Null Hypothesis: independence
# A tibble: 1,000 × 2
   replicate  stat
       <int> <dbl>
 1         1 0.162
 2         2 0.198
 3         3 1.18
 4         4 0.328
 5         5 1.21
 6         6 3.00
 7         7 1.93
 8         8 0.605
 9         9 0.446
10        10 1.10
# … with 990 more rows
# i Use `print(n = ...)` to see more rows

最后，我们可以使用观察到的 F 统计量和从排列数据集中获得的 F 统计量分布来计算 p 值。当 p 值小于预设的显著性水平（例如，0.05）时，我们可以拒绝零假设，并说各组均值之间存在显著差异：

p_value = anova_results %>%
  get_p_value(obs_stat = anova_results, direction = "right") %>%
  pull()
>>> p_value
0.376

结果表明，我们没有足够的信心拒绝零假设。

摘要

在本章中，我们介绍了用于假设检验的不同类型的统计推断，针对数值数据和分类数据。我们介绍了单变量、双变量和多变量的推断方法，使用比例（用于分类变量）或均值（用于数值变量）作为样本统计量。假设检验程序，包括基于模型近似的方法和基于自助法模拟的非参数方法，提供了如置信区间和 p 值等有价值的工具。这些工具使我们能够做出决定，是否可以拒绝原假设而接受备择假设。这样的决定也与第一类和第二类错误相关。

在下一章中，我们将介绍最广泛使用的统计和机器学习模型之一：线性回归。