基于高斯随机扰动的渐近协方差估计
张疏影，中国人民大学统计学院2023级专业硕士研究生。
今天和大家分享Jing Zhou, Wei Lan, and Hansheng Wang (2022). Asymptotic Covariance Estimation by Gaussian Random Perturbation. 该文章于2022年发表在Computational Statistics & Data Analysis上。
背景
-估计量是一大类估计量，通过最小化（或最大化）一个适当定义的损失函数来获得(Lehmann and Casella, 1983)。这个大类中有许多不同的估计量。例如，传统线性回归模型的损失函数可以定义为残差平方和(Casella et al., 2017)，由此得到的普通最小二乘估计量就是一个 -估计量。除此之外，广义最小二乘估计量也是一个 -估计量。事实上，-估计量也被称为最大似然估计量的泛化形式，所有的最大似然估计量都是 -估计量。此外，通过添加一些适当选择的惩罚函数，各种收缩估计量（例如 LASSO 和 SCAD）也是 -估计量。最后，如果参数通过一组矩条件来确定，广泛使用的广义矩估计法可以被转化为一个最小化问题。因此，它也是一个 -估计量(Wooldridge, 2001)。
在大多数情况下，经过精心定义的-估计在适当的条件下是一致的且渐近正态的。-估计的渐近协方差矩阵通常采用三明治形式(Shao, 2003)。在一些情况下（例如，最大似然估计），三明治形式可以进一步简化为非三明治形式。对于渐近有效的统计推断（例如假设检验和置信区间），需要估计三明治型渐近协方差矩阵。解决这个问题的一种方法是获得其解析公式，然后用适当的估计替换未知参数。如果渐近协方差矩阵的公式具有简单形式，则此方法是有效的。例如，在大多数情况下，各种最小二乘估计和最大似然估计都有简单的解析解。
然而，渐近协方差矩阵涉及损失函数的一阶和二阶导数的计算，如果损失函数很复杂，那么这将变得很繁琐。一个典型的例子是带有缺失数据的回归模型。在这种情况下，完整的似然函数涉及一些未知的干扰参数，特别是这些干扰参数通常与缺失机制有关并需要积分消除(Shao and Wang, 2002; Wang and Dai, 2008; Lin et al., 2019; Zhou et al., 2020)。因此，由此得到的渐近协方差矩阵涉及一个具有复杂积分的目标函数的一阶和二阶导数，因此很难计算和估计(Chen et al., 2015; Zhao and Shao, 2015)。在不知道其解析公式的情况下估计渐近协方差矩阵是一个很重要的问题。为了缓解这个问题，重采样类型的方法如自助法和交叉验证法已被提出并广泛使用(Efron and Stein, 1981; Efron and Gong, 1983; Efron and Tibshirani, 1986; Efron, 1994; Jiao and Han, 2020)。它们可以在不知道解析公式的情况下一致地估计渐近协方差矩阵，避免了计算某些复杂积分函数的导数。然而，这种重采样方法也受到计算复杂度的影响。对于传统的数据分析，当样本量不是很大且数据维度相对较低时，这可能不是一个问题。然而，却可能成为大规模数据集的一个严重负担。在这种情况下，计算-估计本身已经有很大计算量，更进一步的计算实际上是不可行的。
为了解决上述问题，我们提出了一种新方法，称为高斯随机扰动。其关键思想可以概括如下。首先，对于给定的损失函数及其-估计，我们在-估计周围生成一个小的随机扰动。随机扰动是从具有微小变异性的多元正态分布中生成的。因此，随机扰动的-估计在局部仍非常接近原始的-估计。其次，我们重新评估这些局部随机扰动的-估计上的损失函数。通过使用泰勒级数展开的严格数学推导，我们发现在-估计处评估的损失函数的一阶和二阶导数可以进一步近似为两个组成部分，分别是在-估计处评估的损失函数和在随机扰动的-估计处重新评估的损失函数。这表明，三明治型渐近协方差矩阵中的元素可以通过使用损失函数进行近似估计，而不是计算导数。这是一种新的渐近协方差矩阵估计量，我们随后严格证明了得到的协方差估计量是统计一致的。
值得注意的是，上述提出的协方差估计量具有两个重要特点。首先，它不涉及在-估计处计算损失函数的一阶和二阶导数。因此，可以在不知道-估计量的解析公式的情况下自动一致地估计渐近协方差矩阵。其次，所提出的协方差估计量可以用向量形式表示。向量可以自然地分解为不同的元素，这些元素可以分别由不同的计算机同时处理，这使得向量形式更适合于并行计算。通过使用并行策略，我们可以将大规模计算问题分解成许多小问题，然后以并行方式解决(Battey et al., 2018; Jordan et al., 2018; Fan et al., 2019; Li et al., 2020)。这使得该方法在进行大规模数据分析时特别有吸引力，尤其是在计算复杂度和隐私保护非常重要的情况下。
研究方法
为了说明所提出的高斯随机扰动方法的实用性，我们考虑两种类型的-估计量的渐近协方差。第一种是传统的-估计量，在第2.1小节中介绍。第二种是带有未知干扰参数的-估计量，我们将在第2.2小节中讨论。
M-估计量的渐近协方差
设 是从第 个（）受试者收集的独立同分布观测数据。其中 表示 维实数向量空间。我们的主要关注点是基于观测数据 （）对未知参数 进行推断，其中 。这里， 是参数空间，是 的开集。相应的损失函数定义为 。然后，提出了一个 估计量 ，其中 。因此，我们有，其中 表示 相对于 的一阶导数。设 为真实值， 和 分别为 的一阶和二阶导数。在适当的正则性条件下，应该有
并且 。这里，两个未知矩阵分别由 和 给出。需要注意的是，当 是在正确指定的似然函数下获得的最大似然估计量时，应该有 ，其中 是某个正定矩阵。因此， 的渐近协方差变为 。可以看出，这里的关键问题是估计 和 。
接下来我们考虑如何在不知道它们的解析公式的情况下估计和。具体来说，对于给定的损失函数在相应的-估计处进行评估时，我们在周围生成一个小的随机扰动。随机扰动是从一个均值为0，协方差为的多元正态分布中生成的，其中，是预先指定的重复次数。这里代表的单位矩阵。方差是一个精心选择的小正数。应该尽可能选择小的值，以使随机扰动的损失函数保持接近原始损失函数。
首先考虑估计。的一个自然估计量是 。通过泰勒级数展开，可以验证
该式近似成立是因为选择的足够小。在(2.1)两边同时乘以，我们得到。利用的事实，我们将替换为其期望形式，并期望可以通过小偏差地估计。这激励我们构建一个的初始估计量我们有
不过，上述初始估计量是有偏的。它的渐近偏差由以下命题1给出。
命题 1.定义，则有。
基于上述命题，我们可以得到
其中包含一个不可忽略的偏差项。因此，需要进行偏差校正。从(2.2)可以知道，，因此可以用 来近似 ，将这个表达式代入(2.2)，得到以下偏差校正估计量：
偏差校正的估计量 几乎是无偏的。但其方差仍然很大。为了减少方差，我们对 进行不同重复次数的平均，最终得到 的估计量为
这便是 的一致估计量。其严格的理论证明，请参见本小节末提出的定理1。
接下来我们考虑如何估计 。其自然估计量为 。类似于 ，我们使用泰勒级数展开，得到
在两边都乘上 ，我们得到-+。注意到 ，由下面的条件(C2)，有。我们期望能够通过来估计。这启发我们构造一个初步的 估计量， 。 然后，通过泰勒级数展开，我们有
容易验证 ，因为 的定义使得 。因此，。与 类似， 是有偏的，其渐近偏差由命题2给出。
命题 2.定义，那么我们有。
基于命题2的结果，我们有
类似地，为了矫正偏差项，我们定义偏差矫正估计量为通过对所有的重复实验取平均，我们得到的最终估计量为
在适当的正则条件下，我们可以理论上严格地证明 是 的一致无偏估计量。与 相似，计算 不需要知道其解析公式，也不需要计算导数。在提供 和 的理论结果之前，我们首先考虑以下正则条件。
(C1) 假设 至少有 阶连续导数，其中 。此外， 和 在 在 的小邻域内是一致有界的。进一步假设 ，且 是 中的开集。
(C2) 假设 。此外，对于任意的 ， 和 都是正定矩阵，其特征值有界。
条件 (C1) 和 (C2) 是标准的，类似的条件在文献中被广泛使用（例如，Shao (2003)）。基于这两个条件， 和 的理论结果在定理1中给出。
定理 1. 在条件(C1)-(C2)下， 当 时，我们有 和 ，其中 。
上述定理意味着，通过将 设置为当 时足够小， 和 都是一致的。我们的模拟结果表明， 的效果令人满意。
需要注意的是， 和 的公式可以用向量形式表示。例如， 可以重写为 ，其中 ，，。这种向量形式非常适合并行计算（例如，Vegh (2018); Masliah (2019)）。这是因为向量可以自然地分解成不同的元素，这些元素可以由不同的计算机分别同时处理，这使得向量形式方便进行并行计算。
带有干扰参数的M-估计量的渐近协方差
我们接下来考虑具有干扰参数的 估计量。相应的损失函数定义为 ，其中 （）是感兴趣的目标参数，（）是干扰参数，并且具有一致的 估计量。一个典型的例子是具有异方差误差方差的线性回归模型(Greene, 1997; Wooldridge, 2015)。在这种情况下，与方差相关的参数是干扰参数。另一个例子是具有缺失数据的回归模型，在该模型中假设倾向函数具有参数形式；例如，Huang et al. (2005) 和 Ibrahim and Molenberghs (2009)。在这种情况下，与倾向函数相关的参数是干扰参数。基于初步估计量 ，提出了一个两步 估计量 ，其中 。设 和 是真实参数。然后，在下面的条件 (C3) 和 (C4) 下，我们有
因此，可以得到 ，其中
类似于第2.1节中的讨论， 和 可以通过以下方法估计：
其中 ，
以及
其中。
这里的 是在第2.1节中定义的随机扰动。在提供 和 的理论结果之前，我们首先考虑以下正则条件。
(C3) 假设 对于某个 具有至少 阶连续导数。此外，对于 在 的一个小邻域内， 和 均有一致的有界性，同时对于 在 的一个小邻域内也是如此。进一步地，假设 ，且 是 中的开集。
(C4) 假设 且 。此外，对于任意 ， 和 都是正定矩阵，其特征值有界。
定理2表明，当 足够小，且 时， 和 都是 和 的一致估计量。
定理 2. 在条件 (C3)-(C4) 下，当 且 时，我们有 和 。
实际数据分析
为了说明我们的方法的实际用途，我们给出一个实际数据的例子，用于分析客户流失。该研究旨在了解哪些因素影响移动通信行业中客户的流失。数据来自中国的一家移动通信公司，包含 名客户的通话记录信息。对于每个客户 ，我们定义一个二元响应变量，其中 表示该客户停止使用服务，否则 。对应的流失率为 1.17 。文献表明，除传统因素（例如续约时间、开销等）之外，社交因素对客户流失具有重要影响(Nitzan and Libai, 2011)。为了研究可能影响客户流失的因素，我们考虑了五个协变量，分别称为tenure, expense, degree, tightness和 entropy. 具体而言，tenure被定义为客户使用服务的时间长度。Expense被定义为客户在一段时间内使用移动电话的平均成本。
为了解释接下来的三个变量，我们需要定义一个辅助变量，称为邻接矩阵。具体而言，我们假设 个客户的网络结构由邻接矩阵 来捕捉，其中 表示节点 与节点 （）之间存在通话（或被通话），否则 。我们定义 以保证完整性，且对于所有 ，有。然后，degree被定义为 ，表示一个焦点客户在自己的网络中涉及的联系人数量。Tightness被定义为 ，其中 是客户 与其连接成员之间的总通话时间。Entropy被定义为 ，其中 ， 是 和 之间的总通话时间。通常，较大的熵表示平均通话时间更分散。
表 1: 实际数据案例的估计结果
为了研究所提出的因素对客户流失率的影响，我们进行了标准逻辑回归分析。所有变量都已经被标准化，使得均值为0，方差为1。表1报告了通过最大似然估计得出的系数，标准误差（即，），以及使用传统方法和我们提出的方法估计的相应-值（即，）。从表中我们可以看出，所提出的随机扰动方法估计的标准误与传统方法得到的非常相似。这表明我们的方法在实际数据中具有鲁棒性。此外，我们发现所有提出的因素在0.1 的水平上都是显著的，这意味着它们都对解释客户流失行为有显著的影响。
总结
我们在本篇文章中提出了一种高斯随机扰动方法，用于估计一般-估计量的渐近协方差矩阵。其关键思想是在局部-估计量周围生成一个小的随机扰动。通过重新计算在随机扰动的-估计量处的损失函数，我们通过泰勒级数展开获得了损失函数的一阶和二阶导数的估计量。这便得到了一种新的渐近协方差矩阵估计量。然后，我们严格证明了在适当的正则条件下，所得到的协方差估计量是统计一致的。该方法不需要计算损失函数的导数，并且适合并行计算。通过模拟和实际数据分析，我们进一步展示了该方法的实际用处。
为了推广所提出的高斯随机扰动方法的实用性，我们在这里提供了两个可能的未来研究方向。首先，将所提出的方法推广到适应具有不连续一阶导数的损失函数，如分位数回归，具有重要意义。其次，我们的方法可以扩展到非参数回归模型。我们相信这些努力可以显著提高高斯随机扰动概念的价值。
参考文献
Battey, H. , Fan, J. , Liu, H. , Lu, J. , and Zhu, Z. . (2018). “ Distributed testing and estimation under sparse high dimensional models,” The Annals of Statistics. 46, 1352-1382.
Casella, G., Fienberg, S., Olkin, I., (2015), “An Introduction to Statistical Learning with Applications in R,” Seventh Edition, Springer New York Heidelberg Dordrecht London.
Chen, K., Guo, S., Lin, Y., and Ying, Z. (2010). “Least absolute relative error estimation,” Journal of the American Statistical Association, 105(491), 1104- 1112.
Chen, X., Wan, A. and Zhou, Y. (2015). “Efficient quantile regression analysis with missing observations,” Journal of the American Statistical Association. 110, 723–741.
Efron, B. & Stein, C. (1981). “The jackknife estimate of variance,” The Annals of Statistics. 9, 586-596.
Efron, B. & Tibshirani, R. (1986). “Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy,” Statistical Science. 1, 54-75.
Efron, B. & Gong, G. (1983). “A leisurely look at the bootstrap, the jackknife, and cross-validation,” The American Statistician 37, 36-48.
Efron, B. (1994). “Missing data, imputation, and the bootstrap,” Journal of the American Statistical Association. 89, 463-475. 11 Fan, J., Wang, D., Wang, K. & Zhu, Z. (2019). “Distributed estimation of principal eigenspaces,” The Annals of Statistics. 47, 3009-3031.
Greene, W. H., (1997), “Econometric Analysis,” Third Edition., Prentice-Hall , Inc.
Huang, R., Liang, Y., and Carriere, KC. (2005), “The role of proxy information in missing data analysis,” Statistical Methods in Medical Research. 14(5), 457.
Ibrahim, J., and Molenberghs, G. (2009), “Rejoinder on: missing data methods in longitudinal studies: a review,” TEST: An Official Journal of the Spanish Society of Statistics and Operations Research. 18.
Jiao, J., & Han, Y. (2020). “Bias correction with jackknife, bootstrap, and Taylor series,” IEEE Transactions on Information Theory 66, 4392-4418.
Jordan, M. I., Lee, J. D. & Yang, Y. (2018). “Communication-efficient distributed statistical inference,” Journal of the American Statistical Association. 114, 1-14.
Lehmann, E. L, & Casella, G. (1983). Theory of point estimation, Wiley. Lin, H., Liu, W. & Lan, W. (2019). “Regression analysis with individual-specific patterns of missing covariates,” Journal of Business Economic and Statistics.
In Press. Li, X., Li, R., Xia, Z. & Xu, C. (2020). “Distributed feature screening via component wise debiasing,” Journal of Machine Learning Research. 21, 1-32.
Maslian, I., Abdelfattah, A., Haidar, A., Tomov, S., Baboulin, M., & Falcou, J., et al. (2019). “ Algorithms and optimization techniques for high-performance matrix-matrix multiplications of very small matrices,” Parallel Computing. 81, 1-21.
Nitzan, I. & Libai, B. (2011). “Social effects on customer retention,” Journal of Marketing. 75, 24–38.
Shao, J. (2003), Mathematical Statistics, Second Edition, Springer-Verlag New York.
Shao, J. & Wang, H. (2002), “Sample correlation coefficients based on survey data under regression imputation,” Journal of the American Statistical Association.
97, 544-552. Wang, Q. & Dai, P. (2008), “Semiparametric model-based inference in the presence of missing responses,”Biometrika. 95, 721–734.
Vegh, J. (2018). “Introducing the explicitly many-processor approach,” Parallel Computing 75, 28-40.
Wooldridge, J. M. (2001). “Applications of generalized method of moments estimation,” Journal of Economic Perspectives. 15, 87-100.
Wooldridge, J. M. (2015), “Introductory Econometrics A Modern Approach,” Sixth Edition, Cengage Learing.
Zhao, J. & Shao, J. (2015), “Semiparametric pseudo-likelihoods in generalized linear models with nonignorable missing data,” Journal of the American Statistical Association. 110, 1577–1590.
Zhou, J., Liu, J., Wang, F. & Wang, H. (2020). “Autoregressive model with spatial dependence and missing data,” Journal of Business Economic and Statistics. forthcoming.