再抽样和刀切法：计算资源有限条件下大规模数据集分析的一种方法
伍书缘，上海财经大学统计与管理学院新聘研究员。主要研究方向为再抽样方法、统计优化算法、大规模数据统计建模等。目前已在Journal of Business and Economic Statistics, Statistica Sinica， Journal of the Royal Statistical Society. Series C等国际权威期刊上发表多篇论文。
今天要跟大家分享Shuyuan Wu, Xuening Zhu Hansheng Wang于2023年发表在Statistica Sinica上的一篇文章：“Subsampling and Jackknifing: A Practically Convenient Solution for Large Data Analysis with Limited Computational Resources”。全文信息详见文章链接：https://www3.stat.sinica.edu.tw/statistica/J33N3/J33N312/J33N312.html或arxiv链接：https://arxiv.org/abs/2304.06231。
1引言
现代统计分析常常需要处理规模庞大的数据集。与此同时，许多研究者只拥有非常有限的计算资源，他们无法构建庞大的分布式计算系统，例如Hadoop或Spark分布式计算系统，只能依靠有限的计算资源（例如个人笔记本电脑）进行数据分析。因此，如何在这种情况下进行可行的大规模数据分析就是一个至关重要的问题。
为了解决这个问题，多种再抽样方法被提出 (Drineas et al., 2011; Ma et al., 2015, 2020; Mahoney, 2011; Wang, 2019; Wang et al., 2018; Yu et al., 2022)。已有文献的核心思想是设计一种新颖的抽样策略，使得估计量在较小的样本量下能够达到较优的统计效率。例如，Ma et al. (2015)提出根据杠杆得分来选择最优样本的方法；而Wang et al. (2018) 针对逻辑回归问题，提出一种-最优准则来挑选样本；Yu et al. (2022)提出一种最优的泊松抽样方法；Ma et al. (2020) 针针对线性回归，研究了再抽样估计量的渐近分布。
虽然这些再抽样方法都很有效，但这些方法存在两个主要缺点。首先，这些方法不具有普适性，针对不同的分析目标，需要精心设计对应的特定抽样策略。其次，这些方法的计算复杂度较高，在具体的实施过程中会消耗大量时间。对于绝大多数方法，抽样成本至少为，这里代表全样本量。
因此，本文提出一种新颖的再抽样方法以解决已有文献存在的问题。我们所提出的方法有如下的独特特点：首先，方法非常简单，因此能够在几乎所有计算设备上轻易实现。其次，本文使用了刀切法，和其他方法相比，估计量得到了明显的修偏。因此，只要抽样次数足够，我们提出的方法能够在更小的样本量下，达到和其他方法相同的统计效率。除此以外，我们所提出的方法支持全自动且一致的统计推断。实际上，对于大多数实际应用而言，有效的统计推断（例如，置信区间）是不可避免的。然而，当估计量的渐近分布较为复杂时，其解析形式将难以获得。在这里，全自动代表各式各样的统计的标准误差可以被自动计算，而无需参考其渐近分布的解析形式，而一致代表方法可以很容易地应用于许多不同的统计量。
具体而言，本文提出一种基于刀切法的再抽样方法，首先实施多次简单有放回随机抽样，随后，对于每一个再抽样子集，针对目标参数，计算一个刀切法修偏估计量，并将所有子集的刀切法估计量取平均，这形成了最终的估计。我们从理论上证明，估计量具有相合性和渐近正态性。在非常温和的假设下，估计量的统计效率能够渐近与全样本估计量相同。这个性质在子样本量很小的时候仍然成立。这样良好的性质主要归功于刀切法。在得到刀切法修偏估计的同时，还可以获得该估计的标准误差的刀切估计量，由此实现了自动统计推断。为了实际应用的可行性，我们进一步提出了一种基于图形处理器（graphics processing unit，GPU）的算法。经实验表明，该方法非常高效。大量的数值研究被实施以验证方法在有限样本下的表现。
需要指出的是，本文所提出的方法也有明显的局限性：与采用分布式方法计算的全样本均值估计 (Suresh et al., 2017)相比，其计算效率较低。然而，所提出的方法仍然具有独特的价值，因为在以下两种重要情况下，它可能是一种更方便可行的选择。第一种情况是全样本量非常大时，在这种情况下，需要计算所有样本的方法（例如，基于分布式的全样本均值估计）的计算复杂度高，时间成本昂贵。这个问题在没有强大的分布式计算系统作为支撑时尤其明显。然而，对于大多数实际数据分析而言，其对估计精度的实际需求是有限的，但时间支出预算是极其宝贵的。因此，在某种程度上牺牲统计效率以换取更少的时间成本对研究者而言可能更具吸引力。第二种情况是需要自动统计推断时。在这种情况下，如果使用基于分布式的全样本均值估计，则必须显示计算出估计量渐近分布的解析形式。此时最佳方案是存在一种全自动且一致的统计推断方法。而我们所提出的方法在这种情况下同样具有独特的优势。
2方法论
2.1模型和设定
首先给出符号说明和模型介绍。令为第个观测到的独立随机变量，其中且代表全样本量。令为全样本的指标集。令为的某个特定矩。简单起见，这里假设。本文的理论能够被容易地推广到多维矩以及估计等更一般的情形。令为目标参数，其中是一个已知的非线性函数，这里假设充分光滑。为估计，可采用一个简单的矩估计量 = ，其中。
方面起见，定义为全样本（whole sample，WS）估计量，因此代表该估计是基于全样本计算得出的。WS估计量的优势是它具有最优的统计效率。然而，当过大时，的计算较为困难，尤其是在研究者的计算资源非常有限的情况下。因此，亟需提出一种计算可行的方法。在这方面，大量卓越且可行的再抽样方法相继提出(Drineas et al., 2011; Ma et al., 2015, 2020; Mahoney, 2011; Wang, 2019; Wang et al., 2018; Yu et al., 2022)。
令为子样本量，通常远小于。令为抽取子样本的次数。记 为第个抽取到的子样本集，其中 （对任意 , ）通过简单有放回随机抽样，独立地从中抽取。换言之，给定的条件下，独立同分布，选中概率为 （对任意 ）。因此，基于的矩估计量能够被如下计算出：，其中。将所有的再抽样估计量取平均，我们就能得到一个更精确的估计量：。称 为再抽样单轮（subsample one-shot，SOS）估计。SOS估计量和分布式计算文献中提及的单轮（one-shot）估计 (Mcdonald et al., 2009; Zhang et al., 2013; Zinkevich et al., 2011)具有类似的思想。但二者最核心的区别是，分布式文献中的标准单轮估计涉及的子样本之间互不重叠。相反，在所提出的再抽样方法中，由于采用有放回抽样，不同子样本之间允许部分重叠。
2.2SOS估计的方差以及偏差分析
为阐述提出方法的动机，我们首先给出关于SOS估计量偏差以及方差的非正式分析。正式的理论结果将于2.4节中给出。具体而言，通过泰勒展式，可以对作如下近似展开
其中和分别是在处的一阶导和二阶导。随后可得如下式子由式(1) 可得，能够被的方差所近似。令，对符号进行重复使用，令同时代表全样本所包含的信息，即，所产生的-域。回顾一下，给定，对任意以及，间独立同分布。由此可得以及。假设的二阶矩有限，满足。我们有由式(2)，可发现能够被近似，其中。在的条件下，子样本估计量的方差可以进一步近似为，而这个量正是WS估计量的渐近方差。随后，我们研究的方差。对任意的估计量，定义的方差为 。则由式（1），我们有
由（3）可知的主项为，其中。不幸的是，这个量不会随着 增加而得到改善。注意到，因此，当与满足时，的偏差才会是一个可忽略的高阶小项。当使用分布式系统时 (Huang and Huo, 2019; Jordan et al., 2019)，假设较为合理。此时，是分布式设备的个数。因此，通常小于，这里是分配在每台设备上的样本量。然而，对于再抽样方法，假设可能不一定满足。此时，为抽取子样本的次数可能较大。相反，为计算方便，子样本量可能远比小。这使得式（3）产生的偏差可能是不可忽略的。为解决这个问题，考虑寻找一种更优的估计量使得它的偏差能被显著地减小。在这一点上，刀切法（jackknife）是一个非常著名的能够用于减小估计量偏差的方法 (Cameron and Trivedi, 2005; Efron and Stein, 1981)。然而，再抽样方法情形下的刀切法的表现和性质仍不清晰。这启发我们提出一个新颖的刀切法估计量。
2.3刀切法估计量
这一节的目标有两个。首先我们提出估计的刀切法纠偏再抽样（jackknife debiased subsample，JDS）估计量。其次，我们同时提出JDS估计量的刀切法标准误（jackknife standard error，JSE）估计量。
首先，给出能够修正SOS估计量偏差的JDS估计方法。为此，对第个子样本，定义刀切法估计量如下
通过与式（1）类似的分析方法，可得近似等于。随后有 ，以及 。这启发我们提出一个估计偏差的估计量，其形式如下：。随后，对第个子集，可以提出一个如下的修偏估计量：。进一步的，平均个不同子集下的修偏估计量，我们得到最终的JDS估计量。我们将在2.4小节中严格论证比的偏差更小。具体而言， 首先可得 以及 。随后，理论上能够证明，的渐近方差与WS估计量相同。因此，在假设足够大的前提下，基于一个很小的子样本量，就能达到卓越的统计效率。除了能够进行偏差修正，刀切法还能够用于得到标准误估计，也就是JDS估计量的标准差估计。我们接下来介绍其核心思想，首先注意到，由式（2），可得
因为，和是已知的，因此这里的关键因素就是去估计未知参数。此外，由刀切法的定义以及泰勒展示，可得对任意 和 有则立刻得知式（2.3）的结果和式（4）中的未知参数非常接近。注意到遍历不同的和，的样本均值是对的合理近似估计。这启发我们提出如下的JSE（jackknife standard error）估计量我们将在随后的内容中理论证明是的一致估计量。此外，我们还将验证。因此， 同时也是的一致估计量。2.4理论结果
在这一小节中，我们将研究并给出文中提到的三种估计量（SOS，JDS和JSE估计量）。为此，首先给出如下的标准技术性假设。
（C1）（亚高斯分布假设） 假设服从亚高斯（sub-Gaussian）分布，即，存在正的常数使得对任意有。
（C2）（光滑性假设） 定义为的第阶导数，假设对任意，是一个连续函数。
（C3）（再抽样条件） 随着，假设子样本量 。此外，假设 以及 。
条件（C1）是关于协变量的经典假设 (Jordan et al., 2019; Zhu et al., 2022)。。条件（C2） 要求函数需要足够充分地光滑，使得能够在附近进行高阶泰勒展开。值得一提的是，这里的八阶连续假设比传统的光滑性假设略强。这个假设是为得到渐近偏差更加精确的显示形式。如果仅保证估计量的渐近正态性，条件（C2）能够被进一步放宽到要求四阶连续可导 (Lehmann and Casella, 2006; Wu, 1986)。最后，条件 (C3) 陈述了， 和 之间的关系和限制。（C3）要求子样本量需要足够大，使得高阶项的分析较为容易。与此同时，要求以保证所有子样本的一致收敛性。而实际中，这个条件是非常容易满足的。
接下来考虑在有限矩不存在的情况下，如何分析各再抽样估计量的渐近行为。受Shao (2003)的渐近理论启发，这里我们采用一种泰勒展开方法来处理有限矩不存在的情况。具体而言，以为例，由泰勒展开，可得 ， 其中 ，，代表某些高阶小项。正如我们在第2.2小节初步讨论的结果，该式子表明的渐近行为能够被和完全决定。其中是无偏的，其主要决定方差，而的方差是可忽略的小项，其主要决定偏差。因此，能够根据来理解的渐近偏差，根据来得到的渐近方差。具体而言，有如下定理。
Theorem 1.假设条件（C1）--（C3）成立，则有如下式子：，其中， ， ， 且
根据定理1可得如下结论。首先，高阶项与和相比可以忽略不计。此外，的渐近偏差由决定，而其渐近方差由决定。其次，由式（5），可知的偏差受和的影响。和项分别代表由子样本误差和总体误差引起的渐近偏差。对于的方差主项还包括两个量：和项。其中第一项是由子样本误差引起的，而第二项是由总体误差引起的。回顾WS估计量，在前文中我们给出其渐近方差近似等于。因此，为了使SOS估计量达到与相同的渐近统计效率，必须要求。不幸的是，的再抽样误差的阶数是，与无关，因此不会随着时减少。因此，要使SOS估计量渐近偏差为，需要。否则，SOS估计量永远无法达到与WS估计量渐近相同的统计效率。类似于，也可以采用泰勒展开，得到，详细表达式在文章的附录中给出。定义，，以及。定理2给出了JDS估计量的性质。
Theorem 2.假设条件（C1）--（C3）成立，则有 。 其中 ， ， 以及
Theorem 3.定义 ， 进一步假设条件 （C1）--（C3） 成立。JSE估计量对于是比率一致的，即，其中“”表示“依概率收敛”。
最后，为了进行有效的渐近推断，需要研究JDS估计量和SOS估计量的渐近分布。因此，给出以下定理以建立和的渐近正态性。
Theorem 4.假设条件（C1）--（C3）成立，则JDS估计量具有渐近正态性，即， 其中“”表示“依分布收敛”。如果进一步加强条件，假设，则SOS估计量也具有渐近正态性，即。
3数值模拟及实际数据分析
3.1使用有放回抽样的原因
这一小节旨在设计一个基于GPU的算法，用于处理存放在硬盘上的数据，以实现前文介绍的再抽样方法。为此，首先需要明确硬盘上的抽样机制。具体而言，本小节将通过以下步骤详细说明不同抽样机制（即有放回和不放回简单随机抽样）之间的计算效率差异。
（1）首先，假设硬盘上有 个数据点（代表一个大型数据集），见图6左上角。数据集包含 列。第一列是样本 ID（），第二列是感兴趣的变量 。
（2）其次，为了进行随机抽样，可以在 和 之间随机生成一个整数，以确定应该抽取哪个数据行。不失一般性，假设抽中第 行，则将 读入内存。（需要注意的是，这个抽样过程是一个简化版本。实际上，研究者无法通过样本 ID 在硬盘上访问数据行。相反，我们根据样本点在硬盘上的物理地址来引用它，这不是一个非常直接的操作，需要较为复杂的技术）。然后，将指标集 从 更新为 。
（3）第三步将介绍如何进行无放回随机抽样。为此，从 到 随机且独立地生成另一个整数 。 可能是 中已经抽取的单元，这会导致重复抽样。为了避免重复抽样，需要将 与 中已经抽取的每个单元进行比较。若发现 已经存在于 中，则需重新生成 。否则将 更新为 ，并将 读入内存。
（4）假设需要生成 个大小为 的子样本。则索引集的大小约为 。为避免重复抽样，每个抽样单元都需要与 中的每个单元进行比较。这个过程中，平均意义下每个抽样单元的计算成本为。总计算成本为 。这显然是一个昂贵的成本。整个抽样流程图见图6。
（5）最后，如果进行有放回随机抽样，则可以\textbf{避免}：（a）不断更新并比较；或者（b）不断更新。这使得所提出的方法在计算上更加高效。
综上所述，与有放回抽样相比，大规模数据下的无放回抽样在实际操作上实现起来更加困难。因此，对于硬盘上的大规模数据集，更适合使用有放回的抽样方法。在这种情况下，整个过程不需要进行记录和比较操作。接下来，为了进一步说明这一点，我们开展了一项数值实验，比较了硬盘上有放回和无放回抽样的效果。为此，首先从标准二元正态分布中独立地生成，其中。目标参数是总体均值。为了估计，分别基于两种抽样策略计算样本均值。这里分别用和表示基于有放回和无放回抽样的估计量。重复进行次实验。对于第次重复，我们汇报两种抽样策略的平均均方误差（mean square error，MSE）和时间成本（time cost，TC）。表1给出所有模拟结果。
表 1: 不同抽样方式的随机模拟结果从表1中可得出以下结论。首先，两种抽样策略的均方误差（MSE）值是可比的，且随着或的增加而减小。然而，两种策略的时间成本（TC）相差很大。使用有放回抽样的方法比无放回抽样的方法快得多。当增加时，两种策略之间的差距显著增加。例如，考虑，的情形，使用有放回抽样方法完成整个过程只需要秒，而使用无放回抽样方法需要近秒。
3.2基于GPU的算法实现
接下来本小节提出一个基于GPU的快速计算算法。需要注意的是，前文所提出的方法在理论上和实际上都非常有用。从理论上而言，它保证了子样本的样本量较小时，所得到的的估计量的统计效率。从实际应用上讲，它是简单、自动且灵活的。然而，该方法的计算成本很高，因为它不仅需要进行次子抽样，还需要对每个子样本进行次自助法。因此，该方法的计算成本很高。所导致的后果是，它在中央处理器（central processing unit，CPU）上的实现可能是低效的，因为标准的CPU通常只有非常有限的核。为了解决这个问题，以MacBook Pro （13-inch, 2020）为例，它使用Intel Core i5处理器，只有四个核。相比之下，标准的GPU可以拥有成千上万的核。因此，GPU是适用于并行计算的强大工具。同时，我们提出的再抽样方法（尤其是刀切法部分）非常适合并行计算。这启发我们提出一个基于GPU的算法来实现再抽样方法。
一个标准的GPU系统应该具有两个独特的特征。为充分利用GPU的计算性能，需要将它们都纳入到考虑中。
GPU系统的第一个独特的特征是它存在两类通讯成本（communication cost）。见图1。第一类通讯成本是将硬盘（hard disk，HD）上的数据传输到CPU内存（CPU memory，CM）上的时间消耗。这个时间消耗是任意计算系统都存在的标准通讯成本。在我们所提出的算法中，这一类成本主要是由于再抽样而产生的。第二类通讯成本代表将CM上的数据传到GPU内存（GPU memory ，GM）上的时间成本。将数据从CM传输到GM的主要目的是为了准备数据，以便能够并行计算刀切法。因此，这部分通讯成本主要是由于刀切法造成的。需要注意的是，当前的GPU架构不允许GPU直接从HD中读取数据。因此，一个好的算法应同时尽量减少这两种类型的通讯成本，从而避免数据在HD、CM和GM之间进行多次通讯。
GPU系统的第二个独特特性是，它非常适合进行张量类型的并行计算。通过这种类型的计算，GPU系统的并行计算能力可得到充分利用。这表明，刀切法计算应被制定成一个张量类型的计算问题。具体而言，在这里，我们提出了一个三步算法来实现所提出的再抽样刀切法。该过程如图2所示。首先，从硬盘中获取第 个子样本 并将其放置在 CM 中。在这个小节中，我们对符号进行了重复使用：这里假设对于任何 ， 是一个 维向量。接下来，可将第 个子样本规范为一个 的矩阵格式，即 。然后，将 传递给 GM，并将 复制 次，以便可以构造一个三维张量 。接下来，定义一个函数来实现使用刀切法计算所需统计量。然后将这个函数映射到的不同通道上，其中每个表示的一个通道。通过这样的处理，GPU系统可以以并行方式执行刀切法。随后收集每个通道的计算结果，并将其整合为第个子样本的所需统计量和。这产生了最终的估计量。这就是整个GPU算法，具体细节见算法1。
3.3通讯成本和计算成本
为评估所提出方法的有限样本下表现，本节将给出丰富的数值模拟。首先，生成 的大规模数据作为全样本。具体而言，对于每个 ，独立地从二元正态分布中生成一个二维随机变量 ，其均值为，协方差为 ，其中 ，，。随后，定义目标参数为如下所示的相关系数
该参数是关于 的多种矩的复杂非线性变换函数。生成全样本后，它被放置在硬盘上作为单个文件，大约占 吉字节（GB）。可以看出，这个大小的数据几乎无法被读入 CM中。一旦数据被放置在硬盘上，它将在后续的数值模拟中固定不变。为了可靠地评估方法的表现，我们总共重复了 次模拟实验。所有计算都使用 TensorFlow 2.2.0， 在单个 GPU 设备上（NVIDIA Tesla P100）执行。
本小节关注时间成本的表现。具体而言，本节分析算法的通讯成本和计算成本。为此，进一步将通讯成本分为两部分。第一部分是从HD到CM的数据传输所需的时间成本。第二部分是从CM到GM所需的时间成本。接下来，将从变化到，将从变化到。随后，使用来表示第次模拟中得到的第个子样，使用来表示获得的时间成本。基于，可以得到矩阵。然后将从HD传输到CM上，其中相关的时间成本被记录为。计算和所需的计算成本为。因此，总时间成本由给出。它们的平均值由， ，和给出。随后，我们研究了它们与和的关系。
图3中给出了详细结果。从图3可以看出，所有类型的时间成本都随子样本数量 的增加而增加。其中，再抽样所需的通讯成本（即 ）远大于其他两类时间成本。相比之下，刀切法所需的通讯成本（即 ）最小。值得注意的是，仅使用CPU系统时，非常显著。然而，使用 GPU 系统时，相应的时间成本变得可以忽略不计。考虑 ，的情形，可得 ， 以及 。另外，图3中和图3右显示，对于固定的总子样本量 ， 和 随 的增加几乎不变。这个结果展示了 GPU 系统出色的并行能力。从中我们能够获得如下见解：只要计算机内存允许，就应该将子样本大小 设置得尽可能大，以获得更好的计算效率。
接下来，我们展示GPU系统的计算优势。为此，首先定义为第次模拟所需的总时间成本减去由再抽样产生的通讯成本（后者是任意计算系统都会产生的成本）。然后在仅使用CPU的系统上执行相同的算法（具体而言，采用TensorFlow 2.2.0的CPU-only指令）。这形成了仅使用CPU系统所需的总时间成本（减去由再抽样产生的通讯成本），将其记录为。随后计算第次模拟的时间之比，即。将平均比率定义为。我们在图4中报告不同组合的经过对数变换的值。从图4中可以看出，对数AR值始终小于0。这表明，基于GPU的系统所需的计算时间成本始终小于仅使用CPU系统的平均计算时间成本。实际上，所报告的对数AR值对于再抽样的抽样次数（即）似乎相当不敏感。此外，当固定时，随着子样本量的增加，对数AR值会降低。这是因为较大的需要更高的计算成本，由此能够更好地展示出GPU系统的并行计算能力。例如，考虑以及的情况，此时GPU系统的平均时间成本约为0.55秒，而CPU系统的时间成本约为4.83秒。相应的AR值为。这表明，GPU系统所需的计算成本仅约为CPU系统的1.1%。
3.4JSE估计量的模拟结果
本小节主要关注JSE估计量的有限样本表现。为此，我们采用和前一小节相同的模拟设定。值得一提的是，为节省实验时间，硬盘上的数据仅生成一次，一旦生成数据，则基于相同的全样本数据重复进行次实验。具体而言，对于第次实验，计算SOS估计量，JDS估计量和JSE估计量。将和（ ）的样本标准差分别定义为。因此，和分别衡量了给定全样本下，和的变异性。由于，和应该是和真实变异性的良好近似；见定理2和3。接下来，对于第次模拟，定义相对绝对误差为和。随后，我们在图5中给出它们的箱线图。
从图5的左图中可以看出，SOS 和 JDS 估计量的 值相似。它们随着 的增加而减小至 0。这表明在 固定的情况下，较大的 会得到更准确的 JSE 估计量。右图也观察到类似的趋势：我们发现，在 固定的情况下，较大的 会得到更准确的 JSE 估计。总而言之，图5中的两个箱线图表明，所提出的 JSE 估计器在 的情况下是一致的。
3.5JDS估计量的模拟结果
这一小节，我们评估了JDS估计量及其置信区间的有限样本性能。为了比较，我们也同时评估SOS估计量的性能。具体而言，按照与前一小节相同的模拟设置，我们基于硬盘上的固定全样本数据重复了次实验。对于第次实验，我们计算JDS估计量及其相应的JSE估计量。这总共产生个估计量。基于这些估计量，其平均偏差可以计算为，同时可以得到相应的标准误差（SE）。此外，对于每个估计量，我们构造的置信区间：，其中，为标准正态分布的下分位数。随后，定义经验覆盖率，其中是指示函数。SOS估计量的评估也类似。详细结果见表2。
从表2 中可得，SOS和JDS两个估计量在不同的 组合下，标准误（SE）的表现相似。然而，它们的偏差却相差较大。SOS 估计量 的偏差显著高于 。例如，考虑 和 的情况， 的偏差为 ，而 的偏差仅有 。可以看出，前者约为后者的四十倍。此外， 的偏差与其标准误相当。因此， 的置信区间表现较差，其 ECP 显著小于 。相比之下，JDS 估计量的置信区间表现较好，因为 相应的 ECP 值非常接近于 。
3.6实际数据分析
本小节中，我们研究一个真实数据集：美国航空公司数据集。该数据集可在美国统计协会（ASA）的官方网站上获取。航空公司数据集包含大约1.2亿条记录，占据硬盘大约12GB的空间。每个记录包含了美国国内从1987年10月到2008年4月所有商业航班的详细信息。数据集包含13个连续变量和16个分类变量。为方便起见，我们聚焦于13个连续变量。但部分连续变量的记录缺失率较高，因此仅考虑其中缺失率小于10％的变量，具体包括：实际飞行时间（）、计划飞行时间（）、航程（）、起飞延误（）和到达延误（）。如需更详细的变量信息，请参考ASA官方网站：http://stat-computing.org/dataexpo/2009。
接下来对原始数据进行预处理。具体而言，首先对每个变量进行符号对数变换： 。此变换仅用于展示的方便性，否则大量变量（例如 ）存在厚尾现象，因此其有限矩可能不存在。对于每个变换后的变量，我们研究以下参数：均值、标准差和峰度。它们的全样本估计量在表3中给出。这些全样本估计量被视为真实参数，由此我们可以执行与模拟实验中相同的实验。在此实际数据分析中，我们固定 ，并使用不同的 组合进行 次实验。详细结果见表4。
从表4中，可以得到如下结论。首先，注意到样本均值是均值的无偏估计。因此，SOS和JDS估计量都是无偏的。实际上，在这种情况下，这两个估计量完全相同。因此，它们在实验中表现出相近的模拟结果，对于所有五个变量，ECP值都非常接近95%。其次，对于其他两个参数（即标准差和峰度），样本估计量不再是无偏的。因此，SOS和JDS估计量不再相同。正如预期，两个估计量的标准误差（SE）相似。然而，它们的偏差结果显著不同。所有情况下，SOS估计量的偏差均显著大于JDS估计量的偏差。因此，SOS估计量的ECP值显著偏离其名义水平95%。相比之下，JDS估计量的ECP值仍非常接近95%。例如，考虑和时的峰度。SOS估计量的ECP值仅为35.6%。相比之下，JDS估计量的ECP值为95.7%。
图 3: 两类通讯成本的随机模拟结果。不同 (n, K) 组合的时间成本取对数后呈现在图中，其中 K =两类通讯成本的随机模拟结果]{两类通讯成本的随机模拟结果。不同 组合的时间成本取对数后呈现在图中，其中 ， 和 。左图汇报由于再抽样而导致的通讯成本 。中图汇报由于刀切法而导致的通讯成本 。计算成本 在右图中给出。图中展示的时间成本（以对数表示）是基于 次模拟的平均值。
图 4: GPU系统和CPU系统的计算效率比较]{ GPU系统和CPU系统的计算效率比较。其中，AR值以log尺度汇报，抽取的子样本数固定为，和，模拟重复次。
图 5: JDS（浅箱子）和SOS（深箱子）估计量RAE值的箱线图。左图对应的情况，右图对应的情况。模拟重复次。
表 2: SOS估计量和JDS估计量基于次模拟重复的随机模拟结果。
表 3: 使用符号对数变换后，基于整个航空数据集的5个连续变量的描述性统计。包括样本均值 （Mean）、样本标准差 （SD） 和样本峰度 （Kurt） 。
表 4: 航空数据集的模拟结果。模拟重复次。目标参数为经过符号对数变换的均值（Mean），标准差（SD），以及偏度（Kurt）。我们同时比较了SOS估计量和JDS估计量偏差，方差和经验覆盖率的表现，其中名义经验覆盖水平为95%。
参考文献
Cameron, A. C. and Trivedi, P. K. (2005), Microeconometrics: methods and applications, Cambridge university press.
Drineas, P., Mahoney, M. W., Muthukrishnan, S., and Sarlós, T. (2011), “Faster least squares approximation,” Numerische mathematik, 117, 219–249.
Efron, B. and Stein, C. (1981), “The jackknife estimate of variance,” The Annals of Statistics, 586–596.
Huang, C. and Huo, X. (2019), “A distributed one-step estimator,” Mathematical Programming, 174, 41–76.
Jordan, M. I., Lee, J. D., and Yang, Y. (2019), “Communication-efficient distributed statistical inference,” Journal of the American Statistical Association, 114, 668–681.
Lehmann, E. L. and Casella, G. (2006), Theory of point estimation, Springer Science & Business Media.
Ma, P., Mahoney, M. W., and Yu, B. (2015), “A statistical perspective on algorithmic leveraging,” The Journal of Machine Learning Research, 16, 861–911.
Ma, P., Zhang, X., Xing, X., Ma, J., and Mahoney, M. (2020), “Asymptotic analysis of sampling estimators for randomized numerical linear algebra algorithms,” in International Conference on Artificial Intelligence and Statistics, PMLR, pp. 1026–1035.
Mahoney, M. W. (2011), “Randomized algorithms for matrices and data,” Foundations and Trends® in Machine Learning, 3, 123–224.
Mcdonald, R., Mohri, M., Silberman, N., Walker, D., and Mann, G. S. (2009), “Efficient largescale distributed training of conditional maximum entropy models,” in Advances in neural information processing systems, pp. 1231–1239.
Shao, J. (2003), Mathematical statistics, Springer Science & Business Media.
Suresh, A. T., Felix, X. Y., Kumar, S., and McMahan, H. B. (2017), “Distributed mean estimation with limited communication,” in International Conference on Machine Learning, PMLR, pp. 3329–3337.
Wang, F., Zhu, Y., Huang, D., Qi, H., and Wang, H. (2021), “Distributed one-step upgraded estimation for non-uniformly and non-randomly distributed data,” Computational Statistics & Data Analysis, 162, 107265.
Wang, H. (2019), “More efficient estimation for logistic regression with optimal subsamples,” Journal of Machine Learning Research, 20, 1–59.
Wang, H., Zhu, R., and Ma, P. (2018), “Optimal subsampling for large sample logistic regression,” Journal of the American Statistical Association, 113, 829–844.
Wu, C.-F. J. (1986), “Jackknife, bootstrap and other resampling methods in regression analysis,” the Annals of Statistics, 14, 1261–1295.
Yu, J., Wang, H., Ai, M., and Zhang, H. (2022), “Optimal distributed subsampling for maximum quasi-likelihood estimators with massive data,” Journal of the American Statistical Association, 117, 265–276.
Zhang, Y., Duchi, J. C., and Wainwright, M. J. (2013), “Communication-efficient algorithms for statistical optimization,” The Journal of Machine Learning Research, 14, 3321–3363.
Zhu, X., Pan, R., Wu, S., and Wang, H. (2022), “Feature screening for massive data analysis by subsampling,” Journal of Business & Economic Statistics, 40, 1892–1903.
Zinkevich, M., Weimer, M., Smola, A. J., and Li, L. (2011), “Parallelized Stochastic Gradient Descent,” in Advances in Neural Information Processing Systems 23: Conference on Neural Information Processing Systems A Meeting Held December.