监控profile数据的AMEWMA控制图

AMEWMA Control Chart for Monitoring Profile Data

摘要: 在线性轮廓数据的监控中，多变量指数加权移动平均(MEWMA)控制图因其对中、小漂移具有良好检测性能而被广泛研究。该控制图的监控性能依赖光滑参数的选取。然而，实际应用中过程漂移的大小通常是未知的，固定光滑参数往往难以同时兼顾不同幅度漂移的检测性能。因此，提出了一种自适应MEWMA控制图(AMEWMA^*)，基于控制图统计量构建漂移估计量${\widehat{\delta }}^{\ast }{}_t$ ，将光滑参数定义为权重函数$g({\widehat{\delta }}^{\ast }{}_t)$ ，从而实现光滑参数的自适应调整，引入RMI作为性能评价准则。通过统计模拟，对比了AMEWMA^*与AMEWMA、MEWMA及T²控制图的平均运行长度和RMI值，统计结果显示，AMEWMA^*的RMI值低于对比方法，说明其在小、中、大不同漂移场景下具有较快的响应速度。

Abstract: In the monitoring of linear profile data, the Multivariate Exponentially Weighted Moving Average (MEWMA) control chart has been extensively studied due to its strong detection capability for small to moderate process shifts. The monitoring performance of this control chart depends on the selection of the smoothing parameter. However, in practical applications, the magnitude of process shifts is typically unknown, and a fixed smoothing parameter often fails to balance detection performance across shifts of varying magnitudes. To address this limitation, an adaptive MEWMA control chart (AMEWMA^*) is proposed. By constructing a shift estimator ${\widehat{\delta }}^{\ast }{}_t$ based on the control chart statistic, the smoothing parameter is defined as a weight function $g\left({\widehat{\delta }}^{\ast }{}_t\right)$ , enabling adaptive adjustment of the smoothing parameter. Additionally, the Relative Monitoring Index (RMI) is introduced as a performance evaluation criterion. Through statistical simulations, AMEWMA is compared with AMEWMA, MEWMA, and T² control charts in terms of average run length and RMI values. The results demonstrate that AMEWMA^* achieves lower RMI values than the benchmark methods, indicating its enhanced responsiveness across small, moderate, and large shift scenarios.

1. 引言

随着现代制造业对产品质量与过程稳定性的要求不断提高，统计过程控制(SPC)技术的应用领域不断拓展，已从传统的单变量与多变量质量特征监控，延伸至“轮廓数据”(profiles)的过程监控。在大多数SPC应用中，通常假设过程或产品的质量，能够通过单个或多个质量特征的单变量、多变量分布进行充分描述。然而，在实际生产中，某些质量特征更适合由响应变量与一个或多个解释变量之间的函数关系来表示，这种函数关系被称为“轮廓”。轮廓数据的监控也成为近年来SPC领域的重要研究方向。

简单线性轮廓是常见的轮廓类型之一，因其模型结构相对简单被广泛应用。在对简单线性轮廓数据的监控方面，Kang和Albin [1] (2000)设计了两种控制图用于监控轮廓参数：一种是用多元$T^2$ 控制图对截距和斜率进行监控；另一种是利用EWMA控制图和R控制图监控误差的均值和方差。结果表明，EWMA/R控制图在监控小漂移方面表现出色。在此基础上，Kim等[2] (2003)提出对自变量$x$ 进行中心化，使截距和斜率的估计值相互独立。并提出使用三个独立的$EWMA$ 控制图来监控简单线性轮廓的截距、斜率和标准差，该控制图被称为EWMA-3图。该方法提高了控制图对不同参数变化的识别能力，但由于需分别绘制和更新多个控制图，其实现过程较为繁琐，运算效率较低。为解决上述问题，Zou等[3] (2007)提出了MEWMA控制图。通过构建一个联合监控截距、斜率和标准差的多元控制图，实现了参数估计与控制过程的统一简化。该方法不仅在中小幅度漂移的检测上优于EWMA-3控制图，而且提升了整体监控效率。

由于EWMA和MEWMA控制图中的光滑参数是固定的，导致其在监控特定幅度漂移时表现良好。然而现实应用中发生的漂移一般都是未知的，固定参数的控制图不能兼顾小漂移与大漂移的检测性能。因此研究者开始引入自适应机制，通过动态调整光滑参数提升控制图的监控性能。Capizzi等[4] (2010)提出了AEWMA-3控制图，用三个自适应EWMA控制图分别监控线性轮廓的截距、斜率与标准差，相比传统的EWMA-3图，AEWMA-3图能更有效的检测不同幅度的过程漂移，在多种漂移情形下均显示出更优的性能表现。Haq和Khoo [5] (2019)提出了一种自适应MEWMA (AMEWMA)控制图，用于监测多元正态过程中均值向量的未知漂移。引入自适应的权重函数$g({\widehat{\delta }}^{\ast }{}_t)$ ，使光滑参数根据监控过程中漂移的大小动态调整，显著提升了对未知均值漂移的检测能力。在此基础上，Haq [6] (2022)将相同的权重函数应用于轮廓监控领域，构建了自适应MEWMA (AMEWMA)控制图监控简单线性轮廓数据。数值模拟结果表明，AMEWMA控制图的性能明显优于现有控制图。

虽然Haq [6] (2022)提出的AMEWMA控制图，在多种漂移情形下的监控效果优于固定光滑参数的控制图。但是，AMEWMA图的权重函数对过程变化的反应不够连续，影响控制图的稳定性和灵敏度。为提升监控方法对过程变化的适应能力，本文定义了新的权重函数。通过引入两个调节参数𝜆和k，实现光滑参数值的动态调整，从而提升了对多种类型过程漂移的响应速度。为优化控制图性能，引入相对效率值的平均值RMI，通过实验确定使控制图监控性能最优的参数值k。模拟结果表明，本文提出的自适应MEWMA控制图RMI值更小，在监控效果上优于Haq [6] (2022)提出的AMEWMA控制图。

2. 变点模型

假设在时间$t$ 抽取一个大小为$n$ 的简单随机样本$\left(X_i,Y_{it}\right)$ ，$i=1,2,\dots ,n$ ，$t\ge 1$ ，其中$Y_{it}$ 和$X_i$ 分别是响应变量$Y$ 和解释变量$X$ 在$t$ 个样本中的第$i$ 个值。当过程处于受控状态时，假设响应变量和解释变量之间的关系可以用简单线性回归模型来表示：

$Y_{it}=A_0+A_1{X}_i+{\epsilon }_{it},i=1,2,\dots ,n$

其中$A_0$ 和$A_1$ 分别为回归模型的截距和斜率，假设误差项${\epsilon }_{it}$ ，服从均值为零方差为${\sigma }^2$ 的正态分布，即${\epsilon }_{it}\sim N\left(0,{\sigma }^2\right)$ 。

假设在$\tau$ 时刻过程发生了未知漂移，即监控如下变点模型：

$\begin{array}{l}{H}_0:A_{0t}=A_0,A_{1t}=A_1,{\epsilon }_{it}~N\left(0,{\sigma }^2\right),t=1,2,\cdots \\ H_1:\{\begin{array}{l}{A}_{0t}=A_0,A_{1t}=A_1,{\epsilon }_{it}~N\left(0,{\sigma }^2\right),t=1,2,\cdots ,\tau -1,\\ A_{0t}={A^{\prime }}_0,A_{1t}={A^{\prime }}_1,{\epsilon }_{it}~N\left(\mu ,{\sigma }_1^2\right),t=\tau ,\tau +1,\cdots .\end{array}\end{array}$

3. 自适应MEWMA控制图

首先，分别对截距、斜率和标准差进行最小二乘估计，得到t时刻截距、斜率和误差方差的最小二乘估计量分别为：

${\widehat{A}}_{0t}^{}={\overline{y}}_t-{\widehat{A}}_{1t}^{}\overline{x}$ ，

${\widehat{A}}_{1t}^{}=\frac{{\displaystyle {\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}},$

${\stackrel{\wedge }{MSE}}_t=\frac{{\displaystyle {\sum }_{i=1}^n{\left(y_{it}-{\widehat{y}}_{it}\right)}^2}}{n-2}.$

将上述估计量${\widehat{A}}_{0t}$ 、${\widehat{A}}_{0t}^{}$ 和${\stackrel{\wedge }{MSE}}_t$ 进行标准正态化，使其服从标准正态分布，得到的标准化结果如下：

$Z_{1t}=\frac{\sqrt{n}\left({\widehat{A}}_{0t}^{}-A_0\right)}{\sigma }\sim N\left(0,1\right)$

$Z_{2t}=\frac{\sqrt{S_{XX,t}}\left({\widehat{A}}_{1t}^{}-A_1\right)}{\sigma }\sim N\left(0,1\right)$

$Z_{3t}={\text{Φ}}^{-1}\left(F_{{\chi }_{n-2}^2}\left(\frac{\left(n-2\right){\stackrel{\wedge }{MSE}}_t}{{\sigma }^2}\right)\right)\sim N\left(0,1\right)$

其中，$S_{XX}={\displaystyle {\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}$ ，${\text{Φ}}^{-1}(\cdot )$ 是标准正态分布的分布函数的逆函数，$F\left(\cdot ;v\right)$ 是自由度为$v$ 的卡方分布的分布函数。令$Z_t={\left[Z_{1t},Z_{2t},Z_{3t}\right]}^{\prime }$ ，则$Z_t$ 服从多元标准正态分布，即$Z_t\sim N_3\left(0,I\right)$ ，其中，0为$3\times 1$ 的零向量，$I$ 为$3\times 3$ 的单位矩阵。

基于$Z_t$ 建立AMEWMA控制图对截距、斜率和标准差同时进行监控。计算AMEWMA统计量：

${\text{J}}_t=g\left({\widehat{\delta }}_t^{\ast }\right){\text{Z}}_t+\left(1-g\left({\widehat{\delta }}_t^{\ast }\right)\right){\text{J}}_{t-1}$

其中$J_0=0$ 是3维零向量，${\widehat{\delta }}_t^{*}$ 是漂移的估计量，$g\left({\widehat{\delta }}^{\ast }{}_t\right)$ 是权重函数。当${J^{\prime }}_t{J}_t$ 大于控制线$h$ 时，控制图发出过程失控的警报。

权重函数$g\left({\widehat{\delta }}^{\ast }{}_t\right)$ 定义如下：

$g\left({\widehat{\delta }}^{\ast }{}_t\right)=\{\begin{array}{l}\lambda /4,{\widehat{\delta }}^{\ast }{}_t\le \text{2}k\\ \lambda /2,{\widehat{\delta }}^{\ast }{}_t\le \text{3}k\\ \text{2}\lambda ,{\widehat{\delta }}^{\ast }{}_t\le \text{4}k\\ \text{4}\lambda ,{\widehat{\delta }}^{\ast }{}_t>\text{4}k\end{array}$

其中$\lambda$ 类似于EWMA控制图中的光滑参数。这里当$\lambda$ 取值0.2时，$\lambda /4=0.05$ 适合监控小漂移；$\lambda /2=0.1$ 适合监控中等漂移，$4\lambda =0.8$ 适合监控大漂移。参数k可以理解为漂移的临界值，选取使控制图的RMI (Han和Tsung [7] (2006))达到最小的临界值，从而保障该控制图针对多种的漂移量综合表现最好。

漂移${\delta }_t$ 的估计量${\widehat{\delta }}_t^{*}$ 采用Dai等[8] (2011)和Haq [6] (2022)的方法。首先，求MEWMA统计量

$A_t=\theta {\text{Z}}_t+\left(1-\theta \right)A_{t-1}$

其中$A_0=0$ ，$\theta$ 是光滑参数($\theta \in (0,1]$ )。通过统计量$A_t$ 计算${{\delta }^{\prime }}_t{\delta }_t^{}$ 的无偏估计量(无偏性的证明见附录)

${\widehat{\delta }}_t^2=\frac{1}{{\left(1-{\left(1-\theta \right)}^t\right)}^2}\left[A_t^{\prime }{A}_t-\frac{3\theta \left(1-{\left(1-\theta \right)}^{2t}\right)}{\left(2-\theta \right)}\right]$

基于${\widehat{\delta }}_t^2$ 构造漂移${\delta }_t$ 的估计量

${\widehat{\delta }}^{\ast }{}_t=+\sqrt{\left|{\widehat{\delta }}_t^2\right|}.$

4. 数值模拟

在接下来的数值模拟中，对比了本文提出的自适应MEWMA控制图(简记为AMEWMA^*)与Haq [6] (2022)提出的AMEWMA控制图、传统MEWMA控制图以及$T^2$ 控制图在不同漂移下的监控效果。为了讨论漂移估计方法对控制图的影响，本文参考廉惠然等[9] (2023)通过马氏距离估计均值漂移量的方法，得到漂移估计量${\widehat{\delta }}^{\ast }{}_t=\sqrt{Z_t\text{'}{Z}_t}$ ，然后利用上一节的权重函数建立自适应控图(简记为AME)，将其作为另一个对比方案。

比较控制图的监控性能时，平均运行长度(ARL)和相对平均指数(RMI)是两个重要的指标。其中，受控状态下的ARL (记为$AR{L}_0$ )越大，说明控制图的虚警率越低；失控状态下的ARL (记为$AR{L}_1$ )越小，说明控制图对漂移的敏感度和检测变化的能力越强。RMI值是所有相对效率值的平均值，RMI越小，表明控制图在监测漂移方面的整体性能越好(Zhou等[10] (2012))。

参考Kang和Albin [1] (2000)、Zou，Tsung和Wang [3] (2007)所采用的线性profile模型$y_{it}={\text{3+2x}}_i+{\epsilon }_{it}$ ，其中${\epsilon }_{it}~N\left(0,1\right)$ ，解释变量$x_i\left(i=1,2,3,4\right)$ 的固定取值为2、4、6和8。对$x_i$ 进行中心化处理$x_i^{\ast }=\left(x_i-\overline{x}\right)$ ，使得截距和斜率的最小二乘估计量是相互独立的，简化分析。得到新的profile模型为$y_{it}=1{\text{3+2x}}_i^{\ast }+{\epsilon }_{it}$ ，其中${\epsilon }_{it}~N\left(0,1\right)$ 。在模拟过程中，设定控制图的$AR{L}_0$ 值约为370，光滑参数λ值取为0.2，控制图的运行结果均是基于100,000次重复实验获得的。

为了评估不同控制图的监控效果，分别在截距、斜率和标准差上引入不同大小的漂移来模拟线性轮廓参数的变化。定义截距、斜率和误差标准差的漂移形式分别为：$A_0+{\eta }_1$ 、$A_1+{\eta }_2$ 和${\eta }_{\text{3}}\sigma$ 。其中，${\eta }_{\text{1}}、{\eta }_{\text{2}}和{\eta }_{\text{3}}$ 表示漂移的大小，$\sigma$ 为标准差。当${\eta }_{\text{1}}\text{=}{\eta }_{\text{2}}\text{=0}且{\eta }_{\text{3}}\text{=1}$ 时，过程处于受控状态。表1、表2和表3分别给出了回归模型的$A_0$ 、$A_1$ 和$\sigma$ 发生不同漂移时，各控制图的$AR{L}_1$ 和RMI结果。

Table 1. Comparison of ARL and RMI values for different control charts when $A_0$ drifts

表1. $A_0$ 发生漂移时不同控制图的ARL和RMI值对比

${\eta }_{\text{1}}$	$AMEWM{A}^{\ast }$	$AME$	$AMEWMA$	$MEWMA$	$T^2$
	h = 0.3915	h = 9.5	h = 1.4734	h = 1.4796	h = 14.172
0	370.96	370.1284	370.92	370.0698	370.8666
0.0425	300.0374	358.7396	319.4	332.9862	365.4987
0.079	200.0066	335.7322	234.13	264.4465	351.7919
0.1	152.9324	315.4782	187.37	222.2761	342.9611
0.133	100.8644	284.2703	129.78	164.8279	319.8381
0.177	60.9539	237.7636	79.34	109.7806	287.0489
0.25	30.9163	170.8379	38.71	57.674	230.1347
0.43	10.3236	70.4531	11.2	18.4764	114.8224
0.61	5.0689	30.0412	5.28	9.418	54.3605
1	1.9325	7.1659	2.03	4.4728	12.4
2	1.0207	1.4053	1.03	2.1167	1.4476
RMI	0	2.5078	0.1487	0.7222	4.0225

Table 2. Comparison of ARL and RMI values for different control charts when $A_1$ drifts

表2. $A_1$ 发生漂移时不同控制图的ARL和RMI值对比

${\eta }_{\text{2}}$	$AMEWM{A}^{\ast }$	$AME$	$AMEWMA$	$MEWMA$	$T^2$
	h = 0.3915	h = 9.5	h = 1.4734	h = 1.4796	h = 14.172
0	370.96	370.1284	370.92	370.0698	370.8666
0.0077	300.389	358.0455	320.15	332.7309	365.5008
0.0144	200.6735	335.265	234.21	264.7534	351.6405
0.0185	150.3728	315.2835	185.36	220.3473	339.9255
0.0243	100.8644	284.3571	129.63	166.3119	319.5782
0.0325	60.1469	238.6802	78.8	108.9222	286.0575
0.046	30.4346	169.5377	37.96	56.967	227.8045
0.08	9.8976	67.2651	10.74	17.6298	111.7506
0.11	5.2113	31.0347	5.42	9.6105	56.231
0.175	2.0853	8.12	2.18	4.733	14.2292
0.3	1.0952	2.0257	1.12	2.5247	2.3483
RMI	0	2.5812	0.1501	0.7424	4.1688

Table 3. Comparison of ARL and RMI values for different control charts when $\sigma$ drifts

表3. $\sigma$ 发生漂移时不同控制图的ARL和RMI值对比

${\eta }_{\text{2}}$	$AMEWM{A}^{\ast }$	$AME$	$AMEWMA$	$MEWMA$	$T^2$
	h = 0.3915	h = 9.5	h = 1.4734	h = 1.4796	h = 14.172
1	370.96	370.1284	370.92	370.0698	370.8666
1.016	299.192	305.8158	300.54	313.4657	309.8672
1.043	201.9391	224.0548	202.53	230.5034	229.3764
1.062	150.2859	181.8772	152.73	184.7385	186.7205
1.089	99.8699	137.1698	102.9	135.2346	141.6334
1.126	59.9954	94.6169	61.41	90.1909	99.3484
1.185	30.3204	55.9436	30.64	51.4607	59.4493
1.32	10.2058	21.3717	9.94	20.8522	23.0222
1.45	5.1003	11.0836	5.02	12.0757	11.7928
1.8	1.9375	3.9093	2	5.4436	3.9655
3	1.0807	1.3624	1.09	2.027	1.3369
RMI	0.0043	0.5774	0.0129	0.7162	0.6398

根据表1和表2的实验结果可以看出，当截距和斜率发生小、中和大漂移时，AMEWMA^*控制图的$AR{L}_1$ 值和RMI值均为最小。例如当${\eta }_{\text{1}}\text{=0}\text{.1}$ 时，AMEWMA^*与AME、AMEWMA、MEWMA和$T^2$ 控制图的$AR{L}_1$ 分别为152.9324、315.4782、187.37、222.2716和342.9611；当${\eta }_{\text{2}}\text{=0}\text{.0144}$ 时，AMEWMA^*与AME、AMEWMA、MEWMA和$T^2$ 控制图的$AR{L}_1$ 分别为200.6735、335.265、234.21、264.7534和351.6405。截距发生漂移时，各控制图的RMI值分别为0、2.5078、0.1487、0.7222和4.0225；斜率发生漂移时，各控制图的RMI值分别为0、2.5812、0.1501、0.7424和4.1688。因此AMEWMA^*控制图监控截距和斜率漂移时的监控效果更好。当误差标准差$\sigma$ 发生漂移时，AMEWMA^*控制图与AMEWMA控制图的$AR{L}_1$ 值很接近，当${\eta }_3\text{=1}\text{.089}$ 时AMEWMA^*与AME、AMEWMA、MEWMA和$T^2$ 控制图的$AR{L}_1$ 分别为99.8699、137.1698、102.9、135.2346和141.6334；当${\eta }_3\text{=1}\text{.32}$ 时AMEWMA^*与AME、AMEWMA、MEWMA和$T^2$ 控制图的$AR{L}_1$ 分别为10.2058、21.3717、9.94、20.8522和23.0222。为了进一步评估这些控制图的整体性能，对RMI值进行对比，根据表3的结果可以看出，AMEWMA^*控制图的RMI值最小，整体监控效果最好。综上，与AME控制图相比，AMEWMA^*控制图的漂移估计方法更有效(实证分析时，不再考虑AME控制图)，本文提出的AMEWMA^*控制图在监控过程漂移时表现更好。

5. 实证分析

为了说明所提出方法的有效性，选取Amiri等[11] (2011)提供的皮革行业鞋革染色过程实际数据进行分析。该数据集反映了染色过程中一个关键质量特性——颜色流出量与温度之间的关系。作者把样本分成五份，在25、32、39、46、53℃五个不同温度下对颜色流出量进行检测，并连续记录11天的数据。在此期间，染色过程处于受控状态。其中，温度(℃)为解释变量，颜色流出量为响应变量。颜色流出量与温度之间的关系可以用简单的线性回归模型拟合，定义为$y_j=-0.0509+0.0034x_i$ 。首先，计算出每个控制图的控制线h，$AR{L}_0$ 的取值约为370。模拟生成30个受控的轮廓数据，在第20个轮廓之后，在参数$A_0$ 和$A_1$ 上引入漂移。通过监测控制图统计量是否超出控制线来判断过程是否失控。当统计量超过控制线时，控制图发出失控信号。各控制图的监控效果如图1所示。

Figure 1. Comparison of results

图1. 结果对比

由图1可知，AMEWMA^*控制图在第1个失控样本出现时发出了失控信号；AMEWMA控制图在第3个失控样本出现时发出了失控信号；MEWMA控制图在第5个失控样本出现时发出了失控信号；$T^2$ 控制图在第9个失控样本出现时发出了失控信号。过程失控时，报警越早的控制图表现越好。综上所述，AMEWMA^*控制图具有较好的表现。

6. 结论

本文提出了一种改进的自适应MEWMA控制图(AMEWMA^*)以提高监控简单线性profile数据过程漂移的综合性能。根据漂移大小，动态调整光滑参数的值。发生大漂移时，选择较大的光滑参数；发生中等漂移时，选择中等的光滑参数；发生较小漂移时，选择较小的光滑参数。本方法较容易推广至监控一般线性profile数据等情形。数值模拟结果表明，当截距和斜率发生漂移的时候，AMEWMA^*控制图的性能一致好于AMEWMA、MEWMA和$T^2$ 控制图，尤其在监控过程中的小到中等漂移时表现出更好的敏感度；误差标准差发生漂移时，AMEWMA^*与AMEWMA控制图相比，优势不够明显，性能相近，但是AMEWMA^*控制图具有更小的RMI值，综合表现最好。

基金项目

国家自然科学基金面上项目(12271271)。

附 录：${\widehat{\delta }}_t^2$ 统计特性的理论证明

性质：${\widehat{\delta }}_t^2$ 是${{\delta }^{\prime }}_t{\delta }_t^{}$ 的无偏估计量，即$E\left({\widehat{\delta }}_t^2\right)={\delta }_t^{\text{'}}{\delta }_t$ 。

证明：当过程失控时，$Z_t\sim N_p\left({\delta }_t,I_p\right)$ ，其中${\delta }_t$ 是漂移量，维数$p=3$ ，$I_p$ 为$p\times p$ 的单位矩阵。

因为$A_t=\theta {\text{Z}}_t+\left(1-\theta \right)A_{t-1}=\theta Z_t+\theta \left(1-\theta \right)Z_{t-1}+\theta {\left(1-\theta \right)}^2{Z}_{t-2}+\cdots +\theta {\left(1-\theta \right)}^{t-1}{Z}_1$ ，所以

$E\left(A_t\right)=\theta {\delta }_t{\displaystyle \sum _{i=0}^{t-1}{\left(1-\theta \right)}^i}={\delta }_t\cdot \left(1-{\left(1-\theta \right)}^t\right)$ ,

$Cov\left(A_t\right)={\theta }^2{\displaystyle \sum _{i=0}^{t-1}{\left(1-\theta \right)}^{2i}\cdot I_p}={\theta }^2\cdot \frac{1-{\left(1-\theta \right)}^{2t}}{1-{\left(1-\theta \right)}^2}\cdot I_p.$

化简得

$Cov\left(A_t\right)=\frac{\theta }{2-\theta }\left[1-{\left(1-\theta \right)}^{2t}\right]\cdot I_p$ .

于是有

$tr\left(Cov\left(A_t\right)\right)=\frac{\theta }{2-\theta }\left[1-{\left(1-\theta \right)}^{2t}\right]\cdot P$ .

从而有

$E\left(A_t\text{'}{A}_t\right)=\left[E\left(A_t\right)\right]\text{'}E\left(A_t\right)+tr\left(Cov\left(A_t\right)\right)={\delta }_t^{\text{'}}{\delta }_t{\left(1-{\left(1-\theta \right)}^t\right)}^2+\frac{\theta }{2-\theta }\left(1-{\left(1-\theta \right)}^{2t}\right)\cdot p$ .

进一步得

$\begin{array}{c}E\left({\widehat{\delta }}_t^2\right)=\frac{1}{{\left(1-{\left(1-\theta \right)}^t\right)}^2}\left[E\left(A_t\text{'}{A}_t\right)-\frac{p\theta \left(1-{\left(1-\theta \right)}^{2t}\right)}{2-\theta }\right]\\ =\frac{1}{{\left(1-{\left(1-\theta \right)}^t\right)}^2}\left[{\delta }_t^{\text{'}}{\delta }_t{\left(1-{\left(1-\theta \right)}^t\right)}^2+\frac{\theta }{2-\theta }\left(1-{\left(1-\theta \right)}^{2t}\right)\cdot p-\frac{p\theta \left(1-{\left(1-\theta \right)}^{2t}\right)}{2-\theta }\right]\\ =\frac{1}{{\left(1-{\left(1-\theta \right)}^t\right)}^2}\left({\delta }_t^{\text{'}}{\delta }_t{\left(1-{\left(1-\theta \right)}^t\right)}^2\right)\\ ={\delta }_t^{\text{'}}{\delta }_t.\end{array}$

NOTES

^*通讯作者。