大样本Gamma回归的最优子抽样

Optimal Subsampling for Large-Sample Gamma Regression

摘要: 随着计算机行业的迅猛发展，人类社会逐渐迈入大数据时代。面对大规模右偏性和厚尾分布的数据，Gamma回归模型发挥着非常重要的作用。然而如何快速并准确估计出Gamma回归中感兴趣参数成为值得思考的热点问题。在本文中，我们提出两种两步算法分别有效地逼近Φ已知Gamma回归和Φ未知Gamma回归在全数据下的最大似然估计，从而解决了单参数与双参数大样本Gamma回归估计问题。首先在Φ已知情况下，可证明出在给定全数据下一般子抽样估计量渐近服从正态分布，推导出使估计量渐近均方误差最小的最优子抽样概率。为了进一步降低计算量，我们还提出了另一种最优子抽样概率。由于最优子抽样概率取决于未知参数，我们还提出了单参数两步算法。其次在Φ未知情况下，我们基于单参数两步算法提出了双参数两步算法。最后使用数值模拟表明两种算法的计算效率高，也证实了通过单参数两步算法得到的估计量与双参数两步算法得到的估计量差距不明显。

Abstract: With the rapid development of computer industry, human society is gradually moving into the era of big data. Gamma regression models play a very important role in the face of large-sample right-skewed and thick-tailed data. However, how to quickly and accurately estimate the parameters of interest in Gamma regression has become a hot issue to be considered. In this paper, we propose two two-step algorithms to efficiently approximate the maximum likelihood estimates of Φ known Gamma regression and Φ unknown Gamma regression under full data, respectively, thus solving the single-parameter and two-parameter large-sample Gamma regression estimation problems. Firstly, in the case where Φ is known, it can be shown that the general subsampling estimates asymptotically obey a normal distribution given the full data, and the optimal subsampling probability that minimizes the asymptotic mean square error of the estimates is derived. To further reduce the computational effort, we also propose an alternative optimal subsampling probability. Since the optimal subsampling probability depends on the unknown parameters, we also propose a single-parameter two-step algorithm. Secondly, in the case of Φ unknown, we propose a two-parameter two-step algorithm based on the one-parameter two-step algorithm. Finally, using numerical simulations, we show that the two algorithms are computationally efficient and also confirm that the difference between the estimates obtained by the one-parameter two-step algorithm and those obtained by the two-parameter two-step algorithm is not significant.

1. 引言

随着大数据兴起，科学家们面临的数据规模也逐渐庞大，前所未有的数据量为研究人员提供了前所未有的机遇和挑战。然而如何在有限的计算资源下处理规模巨大、结构复杂的数据越来越受到人们的关注。

针对线性回归模型，目前处理大规模数据比较流行的方法有分治方法(Lin等，2011) [1]、数据流在线更新方法(Schifano等，2016) [2] 以及子抽样方法(Politis等，1999) [3]。Ma和Sun (2015) [4] 在子抽样方法(Politis等) [3] 基础上提出了杠杆抽样方法，其方法成功的关键是利用不同权重构造非均匀抽样概率，使有影响的数据点以高概率抽样(Ma，Yu和Mahoney，2015) [5]。Ma等(2020) [6] 还研究了杠杆抽样估计的渐近正态性和渐近无偏性。还有很多将子抽样方法应用到其它问题上，Drineas等 [7] 针对L2回归模型提出了子抽样算法。Drineas等(2011) [8] 提出两种子抽样算法来解决线性回归中最小二乘逼近问题。Dhillon等(2013) [9] 针对线性回归中的普通最小二乘快速估计问题使用了子抽样算法。最优子抽样方法是子抽样方法的延伸之一，Zhu等(2015) [10] 针对线性回归模型提出了两种子抽样方法，即最优子抽样方法和预测长度的子抽样方法。Wang等 [11] 针对线性回归子抽样提出基于信息的最优子抽样方法(IBOSS)。

在广义线性模型下，Chen和Xie [12] 在广义线性模型(GLM)上利用分治方法处理海量数据。Wang等(2018) [13] 在logistic回归模型上提出最优子抽样方法来有效地逼近在全数据下的最大似然估计，该方法通过最小化估计量的渐近均方误差来确定最优子抽样概率。Yao等(2018) [14] 在Wang等 [13] 基础上研究了多分类logistics回归模型的最优子抽样算法。Wang和Ma (2021) [15] 研究了分位数回归的最优子抽样，推导出一般子抽样估计量的渐近性质以及两种最优子抽样概率。此外，李莉莉等(2022) [16] 针对岭回归模型提出最优子抽样算法，从而处理了变量间存在多重共线性问题，并证明了估计量的渐近性质。牛晓阳等(2022) [17] 对非参数局部多项式回归模型提出最优子抽样算法，提出OPT和PL两种抽样算法。

Gamma回归模型作为广义线性模型中的一个特例 [18]，它在很多方面有所应用，比如周围等(2021) [19] 将Gamma分布应用在海面目标检测中，并介绍了在Gamma分布杂波背景下的检测器；吴林林(2014) [20] 利用Gamma回归研究了雨滴谱对移动双偏振雷达的影响；王俊轶(2019) [21] 利用Gamma回归研究了图像质量的影响因素。针对Gamma回归估计问题，陈超(2015) [22] 总结了双参数Gamma分布参数估计方法，同时利用矩思想，证明了形状参数存在的相合性和渐近正态性。针对双参数Gamma回归多重共线性问题，Muhammad Qasim等(2018) [23] 使用了Liu估计器。根据模拟研究，表明Liu估计器总是优于最大似然估计。左卫兵和钱莉等 [24] 针对Gamma回归模型的多重共线性问题，提出了Gamma回归模型的几乎无偏岭估计方法，同时给出了该估计优于极大似然估计和一般岭估计的存在性定理。

在本文中，我们将最优子抽样方法应用到Gamma回归模型中。本文的剩余部分组织如下：在第二部分中，我们分别建立了单参数Gamma回归和双参数Gamma回归的最大似然估计。在第三部分中，我们首先考虑了在 $\varphi$ 已知情况下Gamma回归的一般子抽样算法。通过三条正则性假设，在给定全数据下，我们建立了一般子抽样算法估计量的相合性和渐近正态性。最后我们考虑了在 $\varphi$ 未知情况下Gamma回归的一般子抽样算法。在第四部分中，我们基于 $\varphi$ 已知情况下一般子抽样估计量的渐进性质，提出两种最优子抽样概率。由于最优子抽样概率取决于未知量，我们提出了单参数Gamma回归的两步算法来近似逼近全数据下的估计量。在 $\varphi$ 未知情况下，我们提出了双参数Gamma回归的两步算法。在第五部分中，我们进行了广泛的仿真研究，以验证我们的方法的有效性。在第六部分中，我们总结了本文。所有证明细节见附录。

2. Gamma回归

假设 ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^N$ 是一个独立且同分布的随机变量序列，当响应变量 $y_i$ 为正偏态且服从伽马分布时，Gamma回归模型非常有用。假设在给定协变量 $x_i\in R^d$ 条件下，响应变量 $y_i\in R$ 服从伽马分布，则下述模型

$\{\begin{array}{l}f\left(y_i;{\mu }_i,\varphi \right)=\exp \left(\frac{y_i/{\mu }_i-\ln \left(1/{\mu }_i\right)}{-\varphi }+\frac{1-\varphi }{\varphi }\ln y_i-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right),y_i>0.\\ g\left({\mu }_i\right)=x_i^{\text{T}}\beta .\end{array}$ (1)

称为Gamma回归模型。其中， ${\mu }_i$ 是响应变量 $y_i$ 的数学期望且有限，即 $E\left(y_i\right)={\mu }_i$ ；假设 $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_d\right)}^{\text{T}}$ 是d维感兴趣参数，属于 $R^d$ 的紧子集； $\Gamma \left(\alpha \right)={\displaystyle {\int }_0^{\infty }{x}^{\alpha -1}}{\text{e}}^{-x}\text{d}x$ 为Gamma函数； $g(\cdot )$ 是连接函数，它建立了样本总体均值与线性预测值之间的关系，在Gamma回归模型中我们一般采用对数连接函数，即 $g\left({\mu }_i\right)=\ln \left({\mu }_i\right)=x_i^{\text{T}}\beta$，从而 ${\mu }_i=\exp \left(x_i^{\text{T}}\beta \right)$ ； $\varphi$ 为离散参数。下面我们分别在 $\varphi$ 已知和 $\varphi$ 未知这两种情况下估计我们所感兴趣的参数 $\beta$。

2.1. $\varphi$ 已知的Gamma回归

我们先考虑 $\varphi$ 已知情况。由于 $\varphi$ 已知，根据(1)式我们很容易可以写出Gamma回归模型，即

$f\left(y_i|x_i;\beta \right)=\exp \left(\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{-\varphi }+\frac{1-\varphi }{\varphi }\ln y_i-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right),y_i>0.$

Gamma回归模型是广义线性模型的一种特殊情况。在Y关于 $\beta$ 的最大化对数似然函数的运算过程中，可以只考虑含参数 $\beta$ 的项，则Y的对数似然函数可以简洁的表达为

$l\left(\beta \right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{-\varphi }\right]}.$ (2)

估计Gamma回归中未知参数 $\beta$ 通常是最大化对数似然函数。然而， ${\hat{\beta }}_{MLE}$ 一般没有封闭解，通常采用迭代方法来进行数值求解，牛顿迭代法是一种常用的迭代方法。根据式(2)，我们可得到下面公式，

${\hat{\beta }}^{\left(t+1\right)}={\hat{\beta }}^{\left(t\right)}+{\left[{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{y_i{x}_i{x}_i^{\text{T}}}{\exp \left(x_i^{\text{T}}{\hat{\beta }}^{\left(t\right)}\right)}}\right]}^{-1}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{y_i{x}_i-\exp \left(x_i^{\text{T}}{\hat{\beta }}^{\left(t\right)}\right)x_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}^{\left(t\right)}\right)}}.$ (3)

2.2. $\varphi$ 未知的Gamma回归

现在我们考虑 $\varphi$ 未知情况。根据(1)式我们可以写出，在已知协变量 $x_i$ 情况下，响应变量 $y_i$ 的概率密度函数为

$f\left(y_i|x_i;\beta ,\varphi \right)=\exp \left(\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{-\varphi }+\frac{1-\varphi }{\varphi }\ln y_i-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right),y_i>0.$

由于 $\varphi$ 未知，所以回归中有两个未知参数 $\beta$ 与 $\varphi$ 待估计，记 $\theta ={\left(\beta ,\varphi \right)}^{\text{T}}$。估计 $\theta$ 通常是最大化对数似然函数，然而， ${\hat{\theta }}_{MLE}$ 一般没有封闭解，所以我们通常使用牛顿迭代算法。Y关于参数 $\theta$ 的对数似然函数可以表达为

$l\left(\theta \right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{-\varphi }+\frac{1-\varphi }{\varphi }\ln y_i-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right]}.$

根据牛顿迭代公式，我们可得到下面公式，

${\hat{\theta }}^{\left(t+1\right)}={\hat{\theta }}^{\left(t\right)}+{\left[I\left({\hat{\theta }}^{\left(t\right)}\right)\right]}^{-1}S\left({\hat{\theta }}^{\left(t\right)}\right).$ (4)

其中， $S\left(\theta \right)$ 、 $I\left(\theta \right)$ 分别是Gamma回归模型的得分函数和Fisher信息矩阵。

在(3)式中，我们发现每次迭代需要 $O\left(N{d}^2\right)$ 时间， $\zeta$ 为所需迭代次数，因此，整个优化过程需要 $O\left(\zeta N{d}^2\right)$ 时间。由于 $\theta$ 是 $d+1$ 维向量，所以(4)式每次迭代需要 $O\left(N{\left(d+1\right)}^2\right)$ 时间，整个优化过程需要 $O\left(\zeta N{\left(d+1\right)}^2\right)$ 时间。然而对于超大样本，(3)式单次运行的所需计算时间太长了，更别说计算(4)式了。所以我们使用子抽样算法来解决超大样本问题。

3. 一般子抽样

对于海量数据，为了缩短其计算时间，我们利用抽样概率 ${\pi }_i\left(i=1,2,\cdots ,N\right)$，从全部数据中有放回随机抽取 $n\left(\ll N\right)$ 个子样本，其中 ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{\pi }_i}=1$， ${\pi }_i$ 既取决于完整数据 $D_N={\left\{\left(x_i,y_i\right)\right\}}_{i=1}^N$，也取决于响应变量。我们需要使用 ${\pi }_i$ 的倒数作为子样本对数似然函数的权重，否则，估计结果往往是有偏差的。我们首先考虑在 $\varphi$ 已知情况下的Gamma回归模型的子抽样。

3.1. $\varphi$ 已知的Gamma回归的子抽样

子样本表示为 $\left(x_i^{*},y_i^{*}\right)$，子抽样概率是 ${\pi }_i^{*}\left(i=1,2,\cdots ,n\right)$， $\tilde{\beta }$ 为子样本估计器，通过最大化(5)式的加权对数似然函数，可得到基于子样本的估计值。

$l^{*}\left(\beta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\pi }_i^{*}}\left[\frac{y_i^{*}/\exp \left({\left(x_i^{*}\right)}^{\text{T}}\beta \right)+{\left(x_i^{*}\right)}^{\text{T}}\beta }{\varphi }\right]}.$ (5)

最大化(5)式以获得基于子样本的估计值 $\tilde{\beta }$。由于 $l^{*}\left(\beta \right)$ 的凸性，最大化可以用牛顿迭代法法来实现，反复应用以下公式，直到 ${\tilde{\beta }}^{\left(t+1\right)}$ 和 ${\tilde{\beta }}^{\left(t\right)}$ 足够接近，

${\tilde{\beta }}^{\left(t+1\right)}={\tilde{\beta }}^{\left(t\right)}+{\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{y_i^{\text{*}}{x}_i^{\text{*}}{\left(x_i^{\text{*}}\right)}^{\text{T}}}{{\pi }_i^{\text{*}}\exp \left({\left(x_i^{\text{*}}\right)}^{\text{T}}{\tilde{\beta }}^{\left(t\right)}\right)}}\right]}^{-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{y_i^{\text{*}}{x}_i^{\text{*}}-\exp \left({\left(x_i^{\text{*}}\right)}^{\text{T}}{\tilde{\beta }}^{\left(t\right)}\right)x_i^{\text{*}}}{{\pi }_i^{\text{*}}\exp \left({\left(x_i^{\text{*}}\right)}^{\text{T}}{\tilde{\beta }}^{\left(t\right)}\right)}}.$ (6)

在 $\varphi$ 已知情况下，我们研究子样本估计量 $\tilde{\beta }$ 的渐近性质。在整篇论文中， $O_p\left(1\right)$ 指依概率有界。 $\Vert v\Vert$ 表示向量v的欧几里德范数，即 $\Vert v\Vert ={\langle v^{\text{T}}v\rangle }^{1/2}$。下面我们列出一些正则性假设。

假设1. 当 $N\to \infty$ 时，有 $G_N=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{y_i{x}_i{x}_i^{\text{T}}}{\varphi \exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}}$ 依概率收敛于正定矩阵，并且 $\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left[\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}\right]}^4}=O_P\left(1\right)$。

假设2. $\frac{1}{N^3}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\Vert x_i\Vert }^8}{{\pi }_i^2}}=O_P\left(1\right)$。如果 $\rho \ge 0$，有 $\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left[\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right]}^{4+2\rho }}=O_P\left(1\right)$， $\frac{1}{N^{3+2\rho }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\pi }_i^{-2-2\rho }{\Vert x_i\Vert }^{4+2\rho }}=O_p\left(1\right)$。

假设3. $\frac{1}{N^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left[y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta \right]}^2}{{\pi }_i}}=O_P\left(1\right)$， $\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{y_i}{\exp \left(x_i^{\text{T}}\beta \right)}{\Vert x_i\Vert }^3}=O_p\left(1\right)$。

假设1中 $l\left(\beta \right)$ 关于 $\beta$ 的二阶导函数是渐近非奇异的，而且 $y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)$ 的四阶矩是依概率有界的，这个假设是很易实现的。在假设2中，当 $\rho \ge 0$ 时， $y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1$ 的四阶矩是依概率有界的，同时，我们还对子抽样概率引入了一些条件，如果子抽样概率选为均匀子抽样概率 ${\pi }_i=N^{-1}$，我们不难发现，协变量 $x_i$ 的 $4+2\rho$ 阶矩依概率有界。在假设1和假设2成立的条件下，假设3对指数分布族均成立，这个假设是合理的。

定理1. 如果假设1、2和3均成立，当 $N\to \infty$ 和 $n\to \infty$ 时，在给定全数据 $D_N$ 的条件下，近似误差 $\tilde{\beta }-{\hat{\beta }}_{MLE}$ 依概率收敛于0，同时收敛速度是 $n^{-1/2}$，也就是说，对任意 $\epsilon >0$，存在一个有限的 ${\Delta }_{\epsilon }$ 和 $n_{\epsilon }$，使得

$P\left(\Vert \tilde{\beta }-{\hat{\beta }}_{MLE}\Vert \ge n^{-1/2}{\Delta }_{\epsilon }|D_N\right)<\epsilon .$

对任意 $n>n_{\epsilon }$ 均成立。

定理2. 如果假设1、2和3均成立，当 $N\to \infty$ 和 $n\to \infty$ 时，在给定全数据 $D_N$ 的条件下，近似误差 $\tilde{\beta }-{\hat{\beta }}_{MLE}$ 的条件分布渐近服从正态分布，即

$V_N^{-\frac{1}{2}}\left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)\to N\left(0,I\right).$

其中

$V_N=G_N^{-1}{V}_n{G}_N^{-1}$

$V_n=\frac{1}{n{N}^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left\{\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right\}}^2{x}_i{x}_i^{\text{T}}}{{\varphi }^2{\pi }_i}}$

定理1和定理2表明，使用足够大的子样本，可以使近似误差 $\tilde{\beta }-{\hat{\beta }}_{MLE}$ 尽可能小，并且其渐近服从正态分布。在这里当n和N都趋于无穷，对它们的相对阶数并没有限制。即使n大于N，这个定理仍然有效，但是过抽样在计算上没有好处，所以在实践中 $n\ll N$ 是很典型的。下面我们考虑在 $\varphi$ 未知情况下的Gamma回归模型的子抽样。

3.2. $\varphi$ 未知的Gamma回归的子抽样

在 $\varphi$ 未知下，假设子样本为 $\left(x_i^{*},y_i^{*}\right)$，子抽样概率是 ${\pi }_i^{*}\left(i=1,2,\cdots ,n\right)$， $\tilde{\theta }$ 为子样本估计器，通过最大化(7)式的加权对数似然函数，可得到基于子样本的估计值。

$l^{*}\left(\theta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{{\pi }_i^{*}}\left[\frac{y_i^{*}/\exp \left({\left(x_i^{*}\right)}^{\text{T}}\beta \right)+{\left(x_i^{*}\right)}^{\text{T}}\beta }{\varphi }+\frac{1-\varphi }{\varphi }\ln y_i^{*}-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right]}.$ (7)

我们假设 ${\delta }_i$ 表示是否选择了第i个数据点的指标变量，即当 $\left(x_i,y_i\right)$ 在子样本时， ${\delta }_i=1$，否则 ${\delta }_i=0$。我们假设 ${\delta }_i$ 的分布是带参数的伯努利分布，即

$P\left({\delta }_i=1|x_i,y_i\right)={\pi }_i,i=1,2,\cdots ,N.$

(7)式可以等价地写成

$l^{*}\left(\theta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\delta }_i}{{\pi }_i^{*}}\left[\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{\varphi }+\frac{1-\varphi }{\varphi }\ln y_i-\frac{\ln \varphi }{\varphi }-\ln \Gamma \left(\frac{1}{\varphi }\right)\right]}.$

最大化(7)式以获得基于子样本的估计值 $\tilde{\theta }$。由于 $l^{\text{*}}\left(\theta \right)$ 的凸性，最大化可以用牛顿迭代法法来实现，反复应用以下公式，直到 ${\tilde{\theta }}^{\left(t+1\right)}$ 和 ${\tilde{\theta }}^{\left(t\right)}$ 足够接近，

${\tilde{\theta }}^{\left(t+1\right)}={\tilde{\theta }}^{\left(t\right)}+{\left[I^{\text{*}}\left({\tilde{\theta }}^{\left(t\right)}\right)\right]}^{-1}{S}^{\text{*}}\left({\tilde{\theta }}^{\left(t\right)}\right).$ (8)

其中 $S^{\text{*}}\left(\theta \right)I^{\text{*}}\left(\theta \right)$ 分别是(7)式的得分函数和Fisher信息矩阵。

4. 最优子抽样

我们发现无论在(6)式，还是在(8)式，最关键的是选择一个最佳的子抽样概率 ${\pi }_i$。使用均匀子抽样概率 $\pi ={\left\{{\pi }_i=N^{-1}\right\}}_{i=1}^N$ 的算法往往不是“最优”的。从定理2可以看出，渐近方差矩阵 $V_N$ 取决于子抽样概率。因此，我们可以导出使渐近方差矩阵最小的最优子抽样概率。由于 $V_N$ 是矩阵，“最小化”含义需要定义，这里我们采用A-最优性和L-最优性准则(Kiefer, 1959 [25] )。在 $\varphi$ 已知情况下，通过采用两种不同的最优性准则，下面的定理分别给出了两种最优子抽样概率，同时这两种最优子抽样概率分别最小化了 $\tilde{\beta }$ 和 $G_N\tilde{\beta }$ 的渐近均方误差。

定理3. 使迹 $tr\left(V_N\right)$ 最小的最优子抽样概率为

${\pi }_i^{A\text{opt}}=\frac{\left|\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right|\cdot \Vert G_N^{-1}{x}_i\Vert }{{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right|\cdot \Vert G_N^{-1}{x}_j\Vert }}\text{,}i=1,2,\cdots ,N.$ (9)

使迹 $tr\left(V_n\right)$ 最小的最优子抽样概率为

${\pi }_i^{L\text{opt}}=\frac{\left|\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right|\cdot \Vert x_i\Vert }{{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right|\cdot \Vert x_j\Vert }}\text{,}i=1,2,\cdots ,N.$ (10)

最优子抽样概率 ${\pi }_i^{A\text{opt}}$ 需要计算 $\Vert G_N^{-1}{x}_i\Vert$，消耗 $O\left(N{d}^2\right)$ 时间。而对于 ${\pi }_i^{L\text{opt}}$ 是需要计算 $\Vert x_i\Vert$，其在计算上只需要消耗 $O\left(Nd\right)$ 时间，所以 ${\pi }_i^{L\text{opt}}$ 在计算上优于 ${\pi }_i^{A\text{opt}}$，同时 ${\pi }_i^{L\text{opt}}$ 也避免对矩阵求逆和减少出现奇异值的现象。要想使用 ${\pi }_i^{L\text{opt}}$ 和 ${\pi }_i^{A\text{opt}}$，我们发现(9)和(10)式的抽样概率均取决于 ${\hat{\beta }}_{MLE}$，这正是需要完整数据来估计参数所得到的，因此最优子抽样概率不能直接使用。为了更加清晰的表达，下面我们描述了在 $\varphi$ 已知情况下Gamma回归估计的两步算法过程。

5. 模拟研究

在本节中，我们将使用第4节中提出的算法1以及算法2在R上进行数值模拟实验，分别评估其性能。假设 $\beta$ 是7维向量，将参数 $\beta$ 的真实值设为 ${\beta }_0={\left(-0.5,-0.5,-0.5,-0.5,0.5,0.5,0.5\right)}^{\text{T}}$。我们考虑样本大小为 $N=10000$ 的完整数据集，其中X分布情况如下，

a) 协变量X服从均值为0，方差为 $\Sigma$ 的多元正态分布，即 $N\left(0,\Sigma \right)$，其中 $\Sigma$ 是一个矩阵，所有对角元素等于1，非对角元素等于0.5；

b) 协变量X服从均值非零，方差为 $\Sigma$ 的多元正态分布，即 $N\left(1,\Sigma \right)$，其中 $\Sigma$ 与a)相同；

c) $X~1/2N\left(1,\Sigma \right)+1/2N\left(-1,\Sigma \right)$，其中 $\Sigma$ 与情况a)相同；

d) 协变量X的分量是独立的，且每个分量都服从指数分布，速率参数为2。

5.1. $\varphi$ 已知的模拟研究

在 $\varphi$ 已知情况下，我们假设 $\varphi =2$，在这一节中我们研究了 $\stackrel{⌣}{\beta }$ 的均方误差和CPU运行时间等，来评估算法1的性能。

5.1.1. 数据集分布、子样本量与 $\stackrel{⌣}{\beta }$ 的MSE

为了评估在算法1的性能，我们对上述四种情况分别进行探讨，研究数据集分布、子样本量与MSE的关系。我们固定第一步子样本量 $n_0=200$，第二步子样本大小 $n=500,1000,1500,2000,2500,3000$，重复执行 $S=1000$ 次，每次循环要计算估计值 $\stackrel{⌣}{\beta }$ 的均方误差。为了进行比较，我们提供了子样本大小为 $n_0+n$ 的均匀子抽样的均方误差结果，并使用1000个bootstrap样本计算了全数据的估计值。图1结果表明，对于四种不同的分布，算法1比均匀子抽样方法更有优势，并且基于抽样概率 ${\pi }_i^{A\text{opt}}$ 的算法1优于基于 ${\pi }_i^{L\text{opt}}$ 的算法1。

Figure 1. n₀ = 200, MSE of different data sets with n not fixed and different sampling methods

图1. n₀ = 200，不同数据集在n不固定和不同抽样方法下的MSE

5.1.2. 子样本量与CPU运行时间

在第一步样本量 $n_0=200$ 和第二步子样本量为 $n=100,200,300,500,700,1000$ 下，我们分别记录了CPU运行时间。经模拟研究，我们在全数据下使用了1000次bootstrap抽样，其CPU运行时间为15.761秒。表1给出了基于抽样概率为 ${\pi }_i^{A\text{opt}}$ 和 ${\pi }_i^{L\text{opt}}$ 的算法1以及均匀子抽样在情况(a)上的结果。每种算法我们重复循环1000次，并记录所使用的CPU时间。由于其它三种数据集的结果相似，因此在这省略它们以节省空间。我们不难发现，均匀子抽样算法计算时间最短。其次，基于 ${\pi }_i^{A\text{opt}}$ 的算法1比基于 ${\pi }_i^{L\text{opt}}$ 的算法需要的计算时间更长，同时这三种算法均比不上在全数据下的bootstrap抽样花费的时间长，这是符合常理的。

Table 1. n0 = 200, CPU times for case (a) with n not fixed and different sampling methods

表1. n₀ = 200，情况(a)在n不固定和不同抽样方法下的CPU时间

5.1.3. 样本量与CPU时间

为了进一步研究在大样本高维情况下算法1的性能，我们将维数扩大至9维，将参数 $\beta$ 的真实值设为 ${\beta }_0={\left(-0.5,-0.5,-0.5,-0.5,0.5,0.5,0.5,0.5,0.5\right)}^{\text{T}}$，并考虑样本量分别为 $N={10}^4,{10}^5,{10}^6$。表2我们探讨了在不同抽样方法下的样本量与CPU时间的关系。在数据情况(a)下，我们每次抽样方法重复循环200次，同时固定第一步和第二步样本大小分别为 $n_0=200$ 、 $n=1000$。

结果表明，随着样本量N的增加，每种抽样方法所使用的CPU时间都有所增加；在相同样本量下，最优子抽样方法相对于全数据方法大大减少了计算时间，计算效率明显提高，同时，基于 ${\pi }_i^{L\text{opt}}$ 的算法1在计算上占有很大优势，这个实验结论和表1的是一样的。

Table 2. n0 = 200, n = 1000 and d = 10, case (a) CPU times under different N and different sampling methods

表2. n₀ = 200，n = 1000和d = 10，情况(a)在N不同和不同抽样方法下的CPU时间

5.2. $\varphi$ 未知的模拟研究

在 $\varphi$ 未知情况下，我们假设初始值 ${\varphi }_0=2$。我们利用算法2分别研究 $\stackrel{⌢}{\theta }$ 和 $\stackrel{⌢}{\beta }$ 的均方误差，进而评估算法2的性能，最后还与算法1进行比较。

5.2.1. 数据集分布、子样本量与 $\stackrel{⌢}{\beta }$ 的MSE

在本节中，我们研究数据集分布、子样本量与 $\stackrel{⌢}{\beta }$ 均方误差的关系，进而评估算法2的性能。我们选取第一步子样本量 $n_0=200$，第二步子样本大小为 $n=500,1000,1500,2000,2500$。为了节省计算时间，本节重复执行 $S=500$ 次，利用算法2计算估计值并分别求出 $\stackrel{⌢}{\beta }$ 的均方误差 ${\text{MSE}}_{\stackrel{⌢}{\beta }}$。为了进行比较，我们提供了子样本大小为 $n_0+n$ 的均匀子抽样的均方误差结果，并使用500次bootstrap抽样计算了在全数据下估计值的均方误差。图2结果表明，对于四种不同的分布，使用最优子抽样算法比均匀子抽样方法更有优势，并且算法2在基于抽样概率 ${\pi }_i^{A\text{opt}}$ 的MSE小于基于 ${\pi }_i^{L\text{opt}}$ 的MSE。值得注意的是，对于EXP后尾数据集，算法2优于均匀抽样，同时基于抽样概率 ${\pi }_i^{A\text{opt}}$ 的算法2与基于 ${\pi }_i^{L\text{opt}}$ 的性能差距不明显。

Figure 2. n₀ = 200, MSE of $\stackrel{⌢}{\beta }$ for different data sets with n not fixed and different sampling methods

图2. n₀ = 200，不同数据集在n不固定和不同抽样方法下的 $\stackrel{⌢}{\beta }$ MSE

5.2.2. 数据集分布、子样本量与 $\stackrel{⌢}{\theta }$ 的MSE

针对情况(c)和情况(d)数据集，我们还考虑了在不同抽样方法和不同第二步子样本量下对算法2的影响。为了进行比较，我们使用 $\stackrel{⌢}{\theta }$ 的均方误差 ${\text{MSE}}_{\stackrel{⌢}{\theta }}$ 作为指标。第一步子样本量、第二步子样本量、迭代次数以及抽样方法与(5.2.1)节一样。图3结果表明，针对情况(c)和情况(d)的数据集，算法2均优于均匀子抽样方法。这个实验结论与前面一致。

5.2.3. $\stackrel{⌢}{\beta }$ 与 $\stackrel{⌣}{\beta }$ 的MSE比较

我们通过比较 $\stackrel{⌢}{\beta }$ 与 $\stackrel{⌣}{\beta }$ 的均方误差，进而比较算法1与算法2的性能。在情况(a)和情况(c)下分别研究子样本量与MSE的关系。固定第一步子样本量 $n_0=200$，第二步子样本量n = 300, 500, 800, 1500, 2000, 2500, 3000, 3200, 3300, 3500为了缩短计算时间，我们重复执行 $S=500$ 次。在图4中，我们发现算法1与算法2的性能差距不大，随着子样本量的逐渐增加，几乎没有差距。

Figure 3. n₀ = 200, MSE of $\stackrel{⌢}{\theta }$ for different data sets with n not fixed and different sampling methods

图3. n₀ = 200，不同数据集在n不固定和不同抽样方法下的 $\stackrel{⌢}{\theta }$ MSE

Figure 4. n₀ = 200, MSE for different data sets with n not fixed and different sampling probabilities

图4. n₀ = 200，不同数据集在n不固定和不同抽样概率下的MSE

6. 结论

在本文中，我们提出了两种算法分别有效地逼近 $\varphi$ 已知和未知的Gamma回归在全数据下的最大似然估计。在 $\varphi$ 已知下，我们首先建立了一般子抽样估计量的渐近分布，利用其渐近性质导出使渐近均方误差最小的最优子抽样概率。根据L-最优性准则，我们还提出了另一种最优子抽样概率，从而缩短了计算时间。由于最优子抽样概率取决于全数据估计量，我们提出单参数两步算法有效地解决了此问题。在 $\varphi$ 未知下，我们使用在 $\varphi$ 已知情况下所提出的两种子抽样概率，由于子抽样概率取决于参数 $\varphi$ 以及全数据估计量，因此我们提出了双参数两步算法解决此问题。最后使用数值模拟分别评估这两种两步算法性能，数值结果均表明了两种算法在大样本高维下提取有价值信息具有巨大潜力。

附录

引理1. 如果假设1和2成立，则

$G_N^{*}-G_N=O_{P|D_N}\left(n^{-\frac{1}{2}}\right),$ (11)

$\frac{1}{N}\frac{\partial l^{*}\left(\beta \right)}{\partial \beta }=O_{P|D_N}\left(n^{-\frac{1}{2}}\right),$ (12)

其中 $G_N^{*}=\frac{1}{Nn}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{y_i^{*}{x}_i^{*}{\left(x_i^{*}\right)}^{\text{T}}}{\varphi {\pi }_i^{*}\exp \left({\left(x_i^{*}\right)}^{\text{T}}{\hat{\beta }}_{MLE}\right)}}.$。

证明(引理1)：通过计算，可直接得出

$E\left(G_N^{*}|D_N\right)=G_N$ (13)

根据赫尔德不等式、假设2以及假设1，我们有

$\begin{array}{c}\mathrm{var}\left(G_N^{*}|D_N\right)=\frac{1}{n{\varphi }^2{N}^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right]}^2{\Vert x_i\Vert }^4}{{\pi }_i}}-\frac{1}{n}{\left(G_N\right)}^2\\ \le \frac{1}{n{\varphi }^2}{\left[\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left[\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}\right]}^4}\right]}^{\frac{1}{2}}{\left[\frac{1}{N^3}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\Vert x_i\Vert }^8}{{\pi }_i^2}}\right]}^{\frac{1}{2}}-\frac{1}{n}{\left(G_N\right)}^2\\ =O_P\left(\frac{1}{n}\right).\end{array}$ (14)

根据(13)式和(14)式结果，并利用切比雪夫不等式，可证明出(11)式。下面证明(12)式，我们通过计算很容易得到

$E\left[\frac{1}{N}\frac{\partial l^{*}\left({\hat{\beta }}_{MLE}\right)}{\partial \beta }|D_N\right]=\frac{1}{N}\frac{\partial l\left({\hat{\beta }}_{MLE}\right)}{\partial \beta }=0.$ (15)

利用假设2和赫尔德不等式，我们有

$\begin{array}{c}\mathrm{var}\left\{\frac{1}{N}\frac{\partial l^{*}\left({\hat{\beta }}_{MLE}\right)}{\partial \beta }|D_N\right\}=\frac{1}{{\varphi }^2{N}^{2}n}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1\right]}^2{\Vert x_i\Vert }^2}{{\pi }_i}}\\ \le \frac{1}{{\varphi }^{2}n}{\left\{\frac{1}{N^3}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\Vert x_i\Vert }^4}{{\pi }_i^2}}\right\}}^{\frac{1}{2}}\cdot {\left\{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1\right]}^4}\right\}}^{\frac{1}{2}}\\ =O_p\left(n^{-1}\right).\end{array}$ (16)

通过(15) (16)式和切比雪夫不等式，可证出(12)式。

证明(定理1)：我们记 $t_i\left(\beta \right)=\frac{y_i/\exp \left(x_i^{\text{T}}\beta \right)+x_i^{\text{T}}\beta }{\varphi }$， $t_i^{*}\left(\beta \right)=\frac{y_i^{*}/\exp \left({\left(x_i^{*}\right)}^{\text{T}}\beta \right)+{\left(x_i^{*}\right)}^{\text{T}}\beta }{\varphi }$。则 $l^{*}\left(\beta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{t_i^{*}\left(\beta \right)}{{\pi }_i^{*}}}$， $l\left(\beta \right)={\displaystyle \underset{i=1}{\overset{N}{\sum }}{t}_i\left(\beta \right)}$。通过计算，我们可以直接得出

$E\left[\frac{1}{N}{l}^{*}\left(\beta \right)|D_N\right]=\frac{1}{N}l\left(\beta \right).$ (17)

$\begin{array}{c}\mathrm{var}\left[\frac{1}{N}{l}^{*}\left(\beta \right)|D_N\right]=\frac{1}{n}\left[\frac{1}{N^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{t_i^2\left(\beta \right)}{{\pi }_i}}-{\left(\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{t}_i\left(\beta \right)}\right)}^2\right]\\ \le \frac{1}{n{N}^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{t_i^2\left(\beta \right)}{{\pi }_i}}\\ =O_P\left(\frac{1}{n}\right).\end{array}$ (18)

最后一个等式是根据假设3得出的。根据(17)式和(18)式以及切比雪夫不等式，在给定全数据 $D_N$ 的条件下，我们有 $1/N{l}^{*}\left(\beta \right)-1/Nl\left(\beta \right)\to 0$。根据定理5.9 (Van der Vaart, 1998 [26] )，我们得到

$\Vert \tilde{\beta }-{\hat{\beta }}_{MLE}\Vert =o_{P|D_N}\left(1\right).$ (19)

使用泰勒定理(Ferguson, 1996 [27] )，我们有

$0=\frac{{\stackrel{\dot{}}{l}}_j^{*}\left(\tilde{\beta }\right)}{N}=\frac{{\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{N}+\frac{1}{N}\frac{\partial {\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{\partial {\beta }^T}\left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)+\frac{1}{N}{R}_j.$ (20)

其中 ${\stackrel{\dot{}}{l}}_j^{*}\left(\tilde{\beta }\right)$ 是 $l^{*}\left(\beta \right)$ 关于 ${\beta }_j$ 的偏导，并且

$R_j={\left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)}^{\text{T}}{\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\frac{{\partial }^2{\stackrel{\dot{}}{l}}_j^{*}\left\{{\hat{\beta }}_{MLE}+\mu \nu \left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)\right\}}{\partial \beta \partial {\beta }^{\text{T}}}\nu \text{d}\mu \text{d}\nu }}\left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right).$ (21)

利用假设3和马尔可夫不等式，当 $\tau \to \infty$ 时，

$P\left(\frac{1}{\varphi nN}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{y_i^{*}/\exp \left({\left(x_i^{*}\right)}^{\text{T}}\beta \right)}{{\pi }_i^{*}}{\Vert x_i^{*}\Vert }^3}\ge \tau |D_N\right)\le \frac{1}{\varphi N\tau }{\displaystyle \underset{i=1}{\overset{N}{\sum }}{y}_i/\exp \left(x_i^{\text{T}}\beta \right){\Vert x_i\Vert }^3}\to 0.$

因此，

$\begin{array}{l}{\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\frac{{\partial }^2{\stackrel{\dot{}}{l}}_j^{*}\left\{{\hat{\beta }}_{MLE}+\mu \nu \left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)\right\}}{\partial \beta \partial {\beta }^{\text{T}}}\nu \text{d}\mu \text{d}\nu }}\\ \le {\displaystyle {\int }_0^1{\displaystyle {\int }_0^1\nu \Vert \frac{{\partial }^2{\stackrel{\dot{}}{l}}_j^{*}\left\{{\hat{\beta }}_{MLE}+\mu \nu \left(\tilde{\beta }-{\hat{\beta }}_{MLE}\right)\right\}}{\partial \beta \partial {\beta }^{\text{T}}}\Vert \text{d}\mu \text{d}\nu }}\\ =O_{P|D_N}\left(N\right).\end{array}$ (22)

将(19)式和(22)式带入(21)式中，于是

$1/N{R}_j=o_{P|D_N}\left(1\right).$ (23)

通过(20)式、(23)式以及引理1，得出

$\begin{array}{c}\tilde{\beta }-{\hat{\beta }}_{MLE}=-{\left[\frac{1}{N}\frac{\partial {\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{\partial {\beta }^{\text{T}}}\right]}^{-1}\left\{\frac{{\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{N}+o_{P|D_N}\left(1\right)\right\}\\ ={\left(G_N^{*}\right)}^{-1}\left\{\frac{{\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{N}+o_{P|D_N}\left(1\right)\right\}\\ =O_{P|D_N}\left(n^{-\frac{1}{2}}\right).\end{array}$ (24)

证明(定理2)：

$\frac{{\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)}{N}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{\left[y_i^{*}/\exp \left({\left(x_i^{*}\right)}^{\text{T}}\beta \right)-1\right]x_i^{*}}{N\varphi {\pi }_i^{*}}}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\eta }_i}.$ (25)

给定 $D_N$， ${\eta }_1,{\eta }_2,\cdots ,{\eta }_n$ 是独立同分布的。其均值

$E\left({\eta }_i|D_N\right)=0.$ (26)

通过假设2和赫尔德不等式，可计算出其方差为

$\mathrm{var}\left({\eta }_i|D_N\right)=n{V}_n=\frac{1}{{\varphi }^2{N}^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1\right]}^2{\Vert x_i\Vert }^2}{{\pi }_i}}=O_P\left(1\right).$ (27)

同时，对任何 $\epsilon >0$，存在 $\rho >0$，

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{n}{\sum }}E\left\{{\Vert n^{-\frac{1}{2}}{\eta }_i\Vert }^{2}I\left(\Vert {\eta }_i\Vert >n^{\frac{1}{2}}\epsilon \right)|D_N\right\}}\\ \le \frac{1}{n^{1+\frac{\rho }{2}}{\epsilon }^{\rho }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}E\left\{{\Vert {\eta }_i\Vert }^{2+\rho }I\left(\Vert {\eta }_i\Vert >n^{\frac{1}{2}}\epsilon \right)|D_N\right\}}\\ \le \frac{1}{n^{1+\frac{\rho }{2}}{\epsilon }^{\rho }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}E\left\{{\Vert {\eta }_i\Vert }^{2+\rho }|D_N\right\}}\\ =\frac{1}{n^{\frac{\rho }{2}}}\frac{1}{{\epsilon }^{\rho }}\frac{1}{{\varphi }^{2+\rho }}\frac{1}{N^{2+\rho }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1\right]}^{2+\rho }{\Vert x_i\Vert }^{2+\rho }}{{\pi }_i^{1+\rho }}}\\ \le \frac{1}{n^{\frac{\rho }{2}}}\frac{1}{{\epsilon }^{\rho }}\frac{1}{{\varphi }^{2+\rho }}{\left\{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left[y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)-1\right]}^{4+2\rho }}\right\}}^{\frac{1}{2}}{\left\{\frac{1}{N^{3+2\rho }}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{{\Vert x_i\Vert }^{4+2\rho }}{{\pi }_i^{2+2\rho }}}\right\}}^{\frac{1}{2}}\\ \to 0.\end{array}$

最后一个等式是利用假设2以及赫尔德不等式得出。根据(25) (26) (27)式以及Lindeber-Feller中心极限定理，在 $D_N$ 条件下，

$\frac{1}{N}{V}_n^{-\frac{1}{2}}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)=n^{-\frac{1}{2}}{\left\{\mathrm{var}\left({\eta }_i|D_N\right)\right\}}^{-\frac{1}{2}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\eta }_i}\to N\left(0,I\right).$

从(24)式，得出

$\tilde{\beta }-{\hat{\beta }}_{MLE}=\frac{1}{N}{\left(G_N^{*}\right)}^{-1}{\stackrel{\dot{}}{l}}_j^{*}\left({\hat{\beta }}_{MLE}\right)+o_{P|D_N}\left(1\right).$ (28)

利用引理1结果，我们可得到

$\begin{array}{c}{\left(G_N^{*}\right)}^{-1}-{\left(G_N\right)}^{-1}=-{\left(G_N\right)}^{-1}\left(G_N^{*}-G_N\right){\left(G_N^{*}\right)}^{-1}\\ =O_{P|D_N}\left(n^{-\frac{1}{2}}\right).\end{array}$ (29)

根据假设1和(27)式，经验证

$V_N=G_N^{-1}{V}_n{G}_N^{-1}=\frac{1}{n}{G}_N^{-1}\left(n{V}_n\right)G_N^{-1}=O_P\left(n^{-1}\right).$ (30)

因此，通过(28) (29) (30)式，我们可得到

$\begin{array}{c}{V}_N^{-\frac{1}{2}}\left(\tilde{\beta }-{\hat{\beta }}_{MIE}\right)=V_N^{-\frac{1}{2}}{N}^{-1}{\left(G_N^{*}\right)}^{-1}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)+o_{P|D_N}\left(1\right)\\ =V_N^{-\frac{1}{2}}{G}_N^{-1}{N}^{-1}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)+V_N^{-\frac{1}{2}}\left[{\left(G_N^{*}\right)}^{-1}-G_N^{-1}\right]N^{-1}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)+o_{P|D_N}\left(1\right)\\ =V_N^{-\frac{1}{2}}{G}_N^{-1}{V}_n^{\frac{1}{2}}{V}_n^{-\frac{1}{2}}{N}^{-1}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)+O_{P|D_N}\left(n^{-\frac{1}{2}}\right).\end{array}$

由于在 $D_N$ 条件下， $N^{-1}{V}_n^{-\frac{1}{2}}{\stackrel{\dot{}}{l}}^{*}\left({\hat{\beta }}_{MIE}\right)\to N\left(0,I\right)$。因此

$V_N^{-\frac{1}{2}}\left(\tilde{\beta }-{\hat{\beta }}_{MIE}\right)\to N\left(0,I\right).$

即证明出定理2。下面利用定理2的结果来证明定理3。

证明(定理3)：利用迹的性质以及柯西–施瓦茨不等式可知，

$\begin{array}{c}tr\left(V_N\right)=tr\left(G_N^{-1}{V}_n{G}_N^{-1}\right)\\ =\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}tr\left\{\frac{{\left\{\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right\}}^2}{{\pi }_i}{G}_N^{-1}{x}_i{x}_i^{\text{T}}{G}_N^{-1}\right\}}\\ =\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left\{\frac{1}{{\pi }_i}{\left[\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right]}^2{\Vert G_N^{-1}{x}_i\Vert }^2\right\}}\\ =\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{1}{{\pi }_i}}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left\{\frac{1}{{\pi }_i}{\left[\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right]}^2{\Vert G_N^{-1}{x}_i\Vert }^2\right\}}\\ \ge \frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left\{\left|\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right|\cdot \Vert G_N^{-1}{x}_i\Vert \right\}}^2}.\end{array}$

当且仅当 ${\pi }_i\propto \left|\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right|\cdot \Vert G_N^{-1}{x}_i\Vert$ 等号成立，即证明出(9)式。

$\begin{array}{c}tr\left(V_n\right)=\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}tr\left\{\frac{{\left\{\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right\}}^2}{{\pi }_i}{x}_i{x}_i^{\text{T}}\right\}}\\ =\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left\{\frac{1}{{\pi }_i}{\left[\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right]}^2{\Vert x_i\Vert }^2\right\}}\\ =\frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{1}{{\pi }_i}}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left\{\frac{1}{{\pi }_i}{\left[\left(y_i/\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)\right)-1\right]}^2{\Vert x_i\Vert }^2\right\}}\\ \ge \frac{1}{n\varphi }{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left\{\left|\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right|\cdot \Vert x_i\Vert \right\}}^2}.\end{array}$

当且仅当 ${\pi }_i\propto \left|\frac{y_i}{\exp \left(x_i^{\text{T}}{\hat{\beta }}_{MLE}\right)}-1\right|\cdot \Vert x_i\Vert$ 等号成立，即证明出(10)式。