CN108171009A - 基于二维自回归模型参数估计的脑电信号间因果关系检测方法 - Google Patents

Abstract

本发明公开了一种基于二维自回归模型参数估计的脑电信号间因果关系检测方法，包括以下步骤：(1)自回归模型进行阶数估计；(2)对原始信号进行加窗处理；(3)结合OPS算法对自回归模型进行参数估计；(4)所得自回归模型应用于Wiener‑Granger因果方法检测脑电信号间的因果关系。该方法最大限度地减少自回归模型参数估计中的干扰项，并提高参数估计结果的精度。

Description

基于二维自回归模型参数估计的脑电信号间因果关系检测 方法

技术领域

本发明属于生物医学领域，具体涉及一种基于二维自回归模型参数估计的脑电信号间因果关系检测方法。

背景技术

由于大脑系统的复杂性，利用数学的方法对脑电信号进行精确建模比较困难，一般情况下脑电信号的动力学方程是未知的，为了对脑电信号进行分析，需要先对系统的模型进行辨识。在研究脑电信号间的线性关系时，将大脑系统看作一个线性系统，而满足一定条件的线性系统均可以用线性自回归模型来表示。自回归模型既包含了过去的输出信息，又包含了过去的输入信息，需要较少的参数集，利用自回归模型，可以有效地减少线性系统辨识的计算工作量。

对于一个平稳时间序列y，如果知道了y在某一时刻n之前的P个时刻的值，可以建立一个基于该时间序列过去P个值的自回归(Auto Regression，AR)模型来描述它，并且可对其将来的值进行预测：

其中,P为该AR模型的阶数；y_n和y_n-i为自回归模型中的“项”，y_n-i为信号y在时刻n延迟i个采样时间单位的采样值；a(i)为y_n-i项的系数；w_n为时间序列y真实值和由式(1)得到的预测值之间的误差项。

采集到的脑电信号都是离散的，线性自回归模型表示系统输入与输出之间的一个线性函数关系，该函数是一个线性的差分方程，可以用二维自回归模型来模拟两个脑电信号x和y的采样值在时域上的离散函数关系：

其中，x_n、y_n、x_n-i和y_n-i称为自回归模型中的“项”，x_n-i和y_n-i为信号x和y在时刻n延迟i个采样时间单位的采样值；q_ij是信号的阶数，即模型中每一项的最大延迟数；a_ij是各项系数，i,j∈{x,y}；w_n,x和w_n,y为时间序列x和y真实值与由二维自回归模型得到的预测值之间的误差项。

对于该模型可以采用最佳参数搜索算法(Optimal Parameter Search,OPS)对它进行阶数和系数的同时估计或采用赤池信息量准则(Akaike Information Criterion,AIC)方法进行阶数估计后再用最小二乘法进行系数估计。不过这两种方法都有缺陷。

1、OPS方法

以(2)式中的信号y为例：

其中，P和Q分别为x序列和y序列的阶数，a(i)和b(j)分别为y(n-i)和x(n-j)项的系数。对于该自回归模型，可以用OPS算法对系数a(i)，b(j)进行估计。OPS方法的主要内容分为以下几步：

第一步：对(3)式中所示模型，人工给定阶数的最大上限，例如：令P＝Q＝10；

第二步：利用筛选线性无关向量组的方法从模型内总共(P+Q)个项中逐个筛选出线性无关项。定义一个“投影距离”的概念，模型中的每一项，都可以计算出相应的投影距离。该投影距离表示了该项在模型中的重要程度，即投影距离的值越大，其所对应的项在模型中的重要度越高，也有利于减少预测误差。

第三步：根据投影距离值来对第二步结果中的项进行再次筛选，保留投影距离值较大的项而抛弃掉投影距离值较小的项。

第四步：对第三步结果保留下来的项利用最小二乘法进行系数估计。

OPS算法在执行时事先设定一个阶数的最大上限，然后在这个阶数范围内筛选出那些可以显著减少预测误差的项，再进行系数估计。而在实际应用中，为了尽可能完整地包含所有正确项，阶数上限的值往往需要设置的比较大，这样会导致其他无用的干扰项增多，筛选结果的精确度不高。另外，在数据有噪声的环境下，会产生无关项的投影距离值大于正确项的投影距离值的情况，会导致筛选过程选择了其他无关项或抛弃了正确项，从而造成结果不够精确。

2、AIC方法

给定一个零均值的m维自回归过程z_n＝[z_n(1),z_n(2),...,z_n(m)]^T

z_n＝Φ₁z_n-1+Φ₂z_n-2+…+Φ_pz_n-p+w_n (3)

其中Φ_i i＝1,2,...,p，是大小为m×m的系数矩阵，p是模型阶数，w_n是零均值独立同分布的向量过程，它的协方差矩阵为Σ。假定对于任意的变量i>0，w_n和z_n-i都是相互独立的。

假设观测值z₁,z₂,...,z_N都是由(4)式中的自回归过程所产生。为了求得阶数p的估计值q，设定一个估计阶数的最大上限K(K>1)，令q取1至K的所有整数值，对每一个q值均用以下最小二乘法对原模型进行估计，得到如下模型：

其中，是w_n的估计值，估计所得的系数矩阵可以通过求解以下方程组来得到：

ZA_j＝z_q+1,j；j＝1,2,...,m (5)

其中

且

Φ_i(j,:)表示Φ_i的第j行。利用最小二乘法对(6)式求解得到：

原模型的预测误差的协方差矩阵表示如下：

其中，

为了选择阶数q，AIC由下式给出：

令q取遍所有的1至K的整数值，使得AIC(q)值达到最小时的q值即为估计所得的阶数。得到估计所得的阶数q后，即可用最小二乘法对模型系数进行估计。

AIC方法是在认定(2)中不同信号具有相同阶数的基础上进行的，即q_xx＝q_yx＝q_xy＝q_yy。但实际情况中，不同信号的阶数往往不相同，AIC方法会使得估计模型比真实模型多出额外的项，导致接下来的系数估计运算量加大且精准度不高。

发明内容

发明目的：为了解决AIC方法不能对自回归模型中不同信号的阶数进行分别估计的困难，同时也为了解决OPS算法对自回归模型参数估计的不够稳定和不够准确的困难，本发明提供了一种基于滑动窗法的自回归模型参数估计方法来检测脑电信号间的因果关系，该方法最大限度地减少自回归模型参数估计中的干扰项，并提高参数估计结果的精度。

技术方案：本发明采用如下技术方案：

基于二维自回归模型参数估计的脑电信号间因果关系检测方法，所述二维自回归模型为：

其中x_n和y_n分别为两个脑电信号的采样值，q_ij是信号的阶数，即模型中每一项的最大延迟数；a_ij是各项系数，i,j∈{x,y}；w_n,x和w_n,y为时间序列x和y真实值与由二维自回归模型得到的预测值之间的误差项；

脑电信号间因果关系检测方法包括以下步骤：

(1)自回归模型进行阶数估计；

(2)对原始信号进行加窗处理；

(3)结合OPS算法对自回归模型进行参数估计；

(4)所得自回归模型应用于Wiener-Granger因果方法检测脑电信号间的因果关系。

所述步骤(1)具体为：

(1-1)定义广义赤池信息准则指数，其计算方式如下：

其中，是估计阶数的集合，N表示脑电信号长度，det(·)表示求方阵的行列式的值，表示预测误差的协方差矩阵；

(1-2)利用赤池信息准则方法估计得到原模型的阶数q_max，接下来的gAIC方法的模型阶数将在这个范围内选择，即i,j∈{x,y}；

(1-3)将第二个脑电信号y的阶数，即q_xy和q_yy，设定为q_max来计算第一个脑电信号x的阶数，即q_xx和q_yx；令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

(1-4)将第一个脑电信号x的阶数，即q_xx和q_yx，设定为最佳阶数和来计算第二个脑电信号y的阶数，即q_xy和q_yy；令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

所述步骤(2)具体为：

(2-1)设置长度为L的滑动窗和步长S，N为脑电信号总长度且S＜L＜N；

(2-2)利用所述滑动窗从信号x和y的头部开始，同时截取长度均为L且位置相同的信号，得到第一组信号；

(2-3)将滑动窗向后移动S个信号距离，再次截取得到一组信号；

(2-4)重复步骤(2-3)，直至滑动窗移动至信号尾部，此时共可得到组信号。

所述步骤(3)具体为：

(3-1)取出步骤(2-1)到(2-4)中所得到的一组信号，假设为和利用序列构建个列向量，利用序列构建个列向量，方法如下：

将上述列向量顺序排列，组成维数为的矩阵

(3-2)从矩阵H的列向量中顺次筛选出它的一个极大线性无关组，组成维数为的W矩阵，记作W＝[w₁w₂…w_R]，其中w_k(1≤k≤R)表示W矩阵的列向量，R是线性无关向量的个数，即W矩阵的列数；

(3-3)利用最小二乘法对W矩阵所对应的模型中的项的系数进行估计，得到系数向量Θ；

(3-4)利用系数向量Θ对W矩阵中每一个列向量w_k都计算出它的投影距离，对于那些存在于H矩阵中但不存在于W矩阵中的列向量，将他们的投影距离设置为0；所以H矩阵中的每一个列向量均有一个属于自己的投影距离值，将它们按H矩阵中列向量的顺序也写成向量形式，即投影距离向量C＝(C₁,C₂,...C_k)且

(3-5)对步骤(2-4)中所得的的组信号，分别利用步骤(3-1)至(3-4)的过程计算每组信号的投影距离向量C；

(3-6)对每一次求得的投影距离向量C相加并求均值，得到向量将中的元素从大到小依次相加，当累积和大于阈值t时即停止，0<t<1；

(3-7)对于步骤(3-6)中被相加的C_i，在W矩阵中保留其对应的列向量，更新W矩阵；

(3-8)对W矩阵中列向量所对应的项，利用最小二乘法对模型系数进行估计。

所述步骤(4)具体为：

(4-1)Wiener-Granger因果影响索引指数的定义为：

其中分别为一维自回归模型中对脑电信号x和y的预测误差的方差；分别为二维自回归模型中对脑电信号x和y的预测误差的方差；

(4-2)利用步骤(1)、(2)和(3)估计所得的系数来计算模型的预测误差，最后由预测误差的方差得到上述WGCI的值。

有益效果：本发明对自回归模型利用gAIC方法进行阶数估计，并对原始信号增加一个滑动窗，多次计算自回归模型中每一项的投影距离并求平均值，然后再对模型中的项数进行约减。与原有的最佳参数搜索方法相比，本发明的有益效果是：有效减少了自回归模型参数估计过程中的干扰项，并且使筛选结果更加准确和稳定，提高了脑电信号间因果关系检测的精度。

附图说明

图1为本发明公开的检测方法流程图；

图2为实施例中的结果对比图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

如图1所示，为本发明公开的基于二维自回归模型参数估计的脑电信号间因果关系检测方法的流程图，所述二维自回归模型为：

其中x_n和y_n分别为两个脑电信号的采样值，q_ij是信号的阶数，即模型中每一项的最大延迟数；a_ij是各项系数，i,j∈{x,y}；w_n,x和w_n,y为时间序列x和y真实值与由二维自回归模型得到的预测值之间的误差项；

对不同信号的阶数分别进行估计，并且对原始信号增加一个滑动窗后进行分组，结合OPS算法，计算每组中每一项的投影距离值并求平均值，然后再对原自回归模型中的项数进行约减和系数估计，具体包括以下步骤：

(1)自回归模型进行阶数估计；

(1-1)定义广义赤池信息准则指数，其计算方式如下：

其中，是估计阶数的集合，N表示脑电信号长度，det(·)表示求方阵的行列式的值，表示预测误差的协方差矩阵；

(1-2)利用赤池信息准则方法估计得到原模型的阶数q_max，接下来的gAIC方法的模型阶数将在这个范围内选择，即i,j∈{x,y}；

(1-3)将第二个脑电信号y的阶数，即q_xy和q_yy，设定为q_max来计算第一个脑电信号x的阶数，即q_xx和q_yx；利用(6)到(8)式中的过程和(11)式的定义，令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

(1-4)将第一个脑电信号x的阶数，即q_xx和q_yx，设定为最佳阶数和来计算第二个脑电信号y的阶数，即q_xy和q_yy；利用(6)到(8)式中的过程和(11)式的定义，令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

(2)对原始信号进行加窗处理；具体步骤为：

(2-1)设置长度为L的滑动窗和步长S，N为脑电信号总长度且S＜L＜N；

(2-2)利用所述滑动窗从信号x和y的头部开始，同时截取长度均为L且位置相同的信号，得到第一组信号；

(2-3)将滑动窗向后移动S个信号距离，再次截取得到一组信号；

(2-4)重复步骤(2-3)，直至滑动窗移动至信号尾部，此时共可得到组信号。

(3)结合OPS算法对自回归模型进行参数估计；具体步骤为：

(3-1)取出步骤(2-1)到(2-4)中所得到的一组信号，假设为和利用序列构建个列向量，利用序列构建个列向量，方法如下：

将上述列向量顺序排列，组成维数为的矩阵H_Y和维数为的矩阵H_X，具体如下：

可以看出，根据(10)式中的二维自回归模型，利用矩阵H_Y的第i行和矩阵H_X的第i行可以对进行估计，可以看作是矩阵H_Y和H_X中的每一个列向量都对应于模型中的某一项。即，H_Y的第i列对应(10)式模型中的y(n-i)项；H_X的第i列对应(14)式模型中的x(n-i)项。

(3-2)从矩阵H的列向量中顺次筛选出它的一个极大线性无关组，组成维数为的W矩阵，记作W＝[w₁w₂…w_R]，其中w_k(1≤k≤R)表示W矩阵的列向量，R是线性无关向量的个数，即W矩阵的列数；具体包括如下步骤：

(3-2-1)考察向量H₁和H₂之间是否线性无关。即利用最小二乘法，用H₁来拟合H₂，假设它们线性相关，则应该存在一个实数θ，使得H₁θ＝H₂。

(3-2-2)根据最小二乘法，计算

(3-2-3)在数据有噪声的情况下，需要给定一个拟合误差的阈值η(一般此阈值较小，例如η＝0.001)，如果有(“|| ||₂”表示求二范数)，则判定H₁与H₂线性相关。然后继续用同样的方法考察H₁与H₃的线性无关性，直到找到一个H_i与H₁是线性无关的。

(3-2-4)将H₁与H_i合并成矩阵W＝[H₁ H_i]，并利用上述方法继续考察W与H_i+1的线性无关性。每当找到与W线性无关的向量时，便将该向量与W矩阵合并，以用与判断与后面向量的线性无关性。直到所有向量均被考察完毕。此时将W记作W＝[w₁ w₂ … w_R]。其中w_k(1≤k≤R)表示W矩阵的列向量，R是线性无关向量的个数，即W矩阵的列数。

(3-3)利用最小二乘法对W矩阵所对应的模型中的项的系数进行估计，得到系数向量Θ；具体方法如下：

(3-3-1)用W矩阵和序列可以构建等式：

其中，向量表示预测误差。Θ＝(Θ(1),Θ(2),...,Θ(R))T表示原双变量AR模型中和W矩阵列向量所对应的项的系数所构成的向量。

(3-3-2)对上式进行最小二乘法的求解：

其中，表示对Θ用最小二乘估计所得结果，上标“T”表示矩阵的转置，上标“-1”表示求逆矩阵。

(3-4)利用系数向量Θ对W矩阵中每一个列向量w_k都计算出它的投影距离，对于那些存在于H矩阵中但不存在于W矩阵中的列向量，将他们的投影距离设置为0；所以H矩阵中的每一个列向量均有一个属于自己的投影距离值，将它们按H矩阵中列向量的顺序也写成向量形式，即投影距离向量C＝(C₁,C₂,...C_k)且投影距离的计算公式为：

(3-5)对步骤(2-4)中所得的的组信号，分别利用步骤(3-1)至(3-4)的过程计算每组信号的投影距离向量C；

(3-6)对每一次求得的投影距离向量C相加并求均值，得到向量将中的元素从大到小依次相加，当累积和大于阈值t时即停止，0<t<1；

(3-7)对于步骤(3-6)中被相加的C_i，在W矩阵中保留其对应的列向量，更新W矩阵；

(3-8)对W矩阵中列向量所对应的项，利用最小二乘法对模型系数进行估计。

(4)所得自回归模型应用于Wiener-Granger因果方法检测脑电信号间的因果关系，包括如下步骤：

(4-1)自回归模型的每个信号不仅依赖于自己过去的信息，还依赖于另一个信号过去的信息。在一维模型中，只包含信号Y，对于信号Y的预测效果可以用预测误差的方差来表示，即其中，Y^-表示Y的过去。在二维模型中，包含信号X和Y，对于信号Y预测误差的方差为如果X因果Y，则应当有小于Wiener-Granger因果影响索引指数的定义为：

其中分别为一维自回归模型中对脑电信号x和y的预测误差的方差；分别为二维自回归模型中对脑电信号x和y的预测误差的方差；

(4-2)利用步骤(1)、(2)和(3)估计所得的系数来计算模型的预测误差，最后由预测误差的方差得到上述WGCI的值。

本实施例以自回归模型生成的模拟信号为例说明本发明的步骤。

本实验共进行1000次，每次实验的对象为长度为1000的模拟信号，由以下模型生成：

其中，序列e_x(n)和e_y(n)为高斯白噪声。本实验中主要对第二个信号y(n)进行参数估计。

实验步骤：

采用具体实施方式中的步骤(1)，利用gAIC方法对实验模型的阶数进行估计。在本实验中，估计所得的信号y中y的阶数为6，x的阶数为3，即估计得原模型中共有9项。

采用具体实施方式中的步骤(2)，对上述模型所生成的信号进行加窗处理。在本试验中，滑动窗长度为800，步长为40。所以原始信号可以被分成6组，分别为原始信号的1至800点、41至840点、81至880点、121至920点、161至960点、201至1000点。

采用具体实施方式中的步骤(3)，其中(3-3-3)中的拟合误差阈值η＝0.001，(3-6)中的阈值t＝0.92。对每组信号分别计算，得到模型中每一项所对应的投影值向量C＝(C₁,C₂,...C₉)，其中C_i,i＝1,2,...,9表示第i项的投影值。并对6组投影值向量求均值向量根据来对模型中的项进行筛选，并对保留下来的项用最小二乘法进行系数估计。

采用具体实施方式中的步骤(4)，重新生成信号并重复进行实验1000次，保留每次实验对模型中项的筛选结果和最终的系数估计结果。

图2中列出了1000次实验中，模型项的筛选结果。在这里，正确项表示的是第二个信号y中的y(n-5)、y(n-6)、x(n-1)和x(n-3)项；其余的项均为错误的项，即y(n-4)、y(n-3)、y(n-2)、y(n-1)和x(n-2)项。某次实验结果正确指的是刚好可以筛选出了正确的四项；结果多余指的是除了选出这四项之外还包含了其他多余的干扰项；结果丢失指的是丢失了正确项中的一项或者多项，即结果中没有完整包含所有正确项。

从图2的结果中能看出，一方面，本发明所提出的方法在阈值大于0.8之后的正确率都要大于原始方法；另一方面，本发明所提出的方法的最佳正确率(在阈值为0.92时所取得)要高于原始方法的正确率(在阈值为0.75时所取得)，证明了本发明方法的有效性。

在本实验中，Wiener-Granger因果影响索引指数的定义为：

在原始的WGCI方法中，首先利用AIC准则来确定模型阶数，然后利用最小二乘方法对模型系数进行估计并计算预测误差，最后由预测误差的方差来得到WGCI；在新的方法中，利用估计所得的系数来计算模型的预测误差及其方差，进而算得WGCI。在本次实验中，进行两种方法的对比，并将结果的均值展示在表1和表2中。

表1两种方法的WGCI_X→Y对比

表2两种方法的WGCI_Y→X对比

从表1中可以看出，两种方法计算出的WGCI_X→Y分别为0.7940与0.8072，WGCI_X→Y分别为0.1970与0.4826，这里WGCI大于0表明了在对X(或Y)的预测过程中，通过引入Y(或X)可以增加预测的准确度，即不论从X到Y还是Y到X均是有Granger因果的，这个结果与实验模型也是相吻合的。而在对WGCI_X→Y的计算过程中，新的方法计算出的WGCI要明显大于原始WGCI方法，这主要是由于新的方法在引入Y的前后能检测出更大的差别(即明显大于)，证明了在对WGCI_X→Y的识别中，新的方法是更加鲁棒有效的。

Claims (5)

1.基于二维自回归模型参数估计的脑电信号间因果关系检测方法，所述二维自回归模型为：

其中x_n和y_n分别为两个脑电信号的采样值，q_ij是信号的阶数，即模型中每一项的最大延迟数；a_ij是各项系数，i,j∈{x,y}；w_n,x和w_n,y为时间序列x和y真实值与由二维自回归模型得到的预测值之间的误差项；

其特征在于，包括以下步骤：

(1)自回归模型进行阶数估计；

(2)对原始信号进行加窗处理；

(3)结合OPS算法对自回归模型进行参数估计；

(4)所得自回归模型应用于Wiener-Granger因果方法检测脑电信号间的因果关系。

2.根据权利要求1所述的脑电信号间因果关系检测方法，其特征在于，所述步骤(1)具体为：

(1-1)定义广义赤池信息准则指数，其计算方式如下：

其中，是估计阶数的集合，N表示脑电信号长度，det(·)表示求方阵的行列式的值，表示预测误差的协方差矩阵；

(1-2)利用赤池信息准则方法估计得到原模型的阶数q_max，接下来的gAIC方法的模型阶数将在这个范围内选择，即

(1-3)将第二个脑电信号y的阶数，即q_xy和q_yy，设定为q_max来计算第一个脑电信号x的阶数，即q_xx和q_yx；令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

(1-4)将第一个脑电信号x的阶数，即q_xx和q_yx，设定为最佳阶数和来计算第二个脑电信号y的阶数，即q_xy和q_yy；令取遍所有可能值，当达到最小值时即可得到最佳阶数估计和

3.根据权利要求1所述的脑电信号间因果关系检测方法，其特征在于，所述步骤(2)具体为：

(2-1)设置长度为L的滑动窗和步长S，N为脑电信号总长度且S＜L＜N；

(2-2)利用所述滑动窗从信号x和y的头部开始，同时截取长度均为L且位置相同的信号，得到第一组信号；

(2-3)将滑动窗向后移动S个信号距离，再次截取得到一组信号；

(2-4)重复步骤(2-3)，直至滑动窗移动至信号尾部，此时共可得到组信号。

4.根据权利要求3所述的脑电信号间因果关系检测方法，其特征在于，所述步骤(3)具体为：

(3-1)取出步骤(2-1)到(2-4)中所得到的一组信号，假设为和利用序列构建个列向量，利用序列构建个列向量，方法如下：

将上述列向量顺序排列，组成维数为的矩阵

(3-2)从矩阵H的列向量中顺次筛选出它的一个极大线性无关组，组成维数为的W矩阵，记作W＝[w₁w₂…w_R]，其中w_k(1≤k≤R)表示W矩阵的列向量，R是线性无关向量的个数，即W矩阵的列数；

(3-3)利用最小二乘法对W矩阵所对应的模型中的项的系数进行估计，得到系数向量Θ；

(3-4)利用系数向量Θ对W矩阵中每一个列向量w_k都计算出它的投影距离，对于那些存在于H矩阵中但不存在于W矩阵中的列向量，将他们的投影距离设置为0；所以H矩阵中的每一个列向量均有一个属于自己的投影距离值，将它们按H矩阵中列向量的顺序也写成向量形式，即投影距离向量C＝(C₁,C₂,...C_k)且

(3-5)对步骤(2-4)中所得的的组信号，分别利用步骤(3-1)至(3-4)的过程计算每组信号的投影距离向量C；

(3-6)对每一次求得的投影距离向量C相加并求均值，得到向量将中的元素从大到小依次相加，当累积和大于阈值t时即停止，0<t<1；

(3-7)对于步骤(3-6)中被相加的C_i，在W矩阵中保留其对应的列向量，更新W矩阵；

(3-8)对W矩阵中列向量所对应的项，利用最小二乘法对模型系数进行估计。

5.根据权利要求4所述的脑电信号间因果关系检测方法，其特征在于，所述步骤(4)具体为：

(4-1)Wiener-Granger因果影响索引指数的定义为：

其中分别为一维自回归模型中对脑电信号x和y的预测误差的方差；分别为二维自回归模型中对脑电信号x和y的预测误差的方差；

(4-2)利用步骤(1)、(2)和(3)估计所得的系数来计算模型的预测误差，最后由预测误差的方差得到上述WGCI的值。