一类改进广义岭估计的实例研究

Improved Generalized Ridge Regression and Its Application

摘要: 本文对广义岭回归方法进行了改进，并在真实数据中进行了实例验证。改进广义岭估计主要对广义岭估计在线性回归模型中当存在若干较大异常值影响模型精度的情况进行了修正，加入适当的修正参数，使模型达到对数据更精确的拟合以及预测作用，并对中俄贸易数据进行了实例验证。

Abstract: In this paper, the generalized ridge regression method is improved and verified by real data. The improved generalized ridge regression mainly corrects the generalized ridge estimation in the linear regression model when there are some large outliers that affect the accuracy of the model. Appropriate correction parameters are added to make the model achieve more accurate data fitting and prediction effects. Finally, we verified it in trade data.

1. 引言

本文主要提出改进广义岭估计。广义岭估计是Horel和Kennard提出的对于一般岭估计岭参数设置改进的一种方法；Hawkins和Yin提出了针对实际情况中样本量n小于未知参数的p形成降秩矩阵的岭回归的快速迭代算法；由于实际研究中还没有一种能够确定岭参数的完美的方法，或多或少会造成估计误差不稳定，游、王和刘(2002) [1] 提出了广义岭估计的直接解法，跳过岭参数的估计，直接求得具有最小均方误差的解；凌和叶对平衡损失下的双h岭估计的优良性进行了验证，给出了风险函数的表达式。改进岭估计在基于前人对该领域研究的基础上，适当改进，对于实际数据若干异常值较大影响回归结果的问题提出新的方法予以解决，将大大提高模型准确度。

2. 模型建立及数据预处理

2.1. 模型建立

在将中蒙俄2008~2018的贸易数据整理，所有变量重新命名。定义 $Y_1$ 和 $Y_2$ 为因变量， $X_1$ 、 $X_2$ 、 $X_3$ 、 $X_4$ 、 $X_5$ 、 $X_7$ 、 $X_8$ 、 $X_{10}$ 、 $X_{12}$ 和 $X_{14}$ 是对应 $Y_1$ 的自变量； $X_1$ 、 $X_2$ 、 $X_3$ 、 $X_4$ 、 $X_6$ 、 $X_7$ 、 $X_9$ 、 $X_{11}$ 、 $X_{13}$ 和 $X_{15}$ 是对应 $Y_2$ 的自变量。所有数据均能在统计局官网获得，这里不予介绍。建立如下线性模型：

$Y_1=X^{\prime }{\beta }_1+\epsilon ,\epsilon ~N\left(0,{\sigma }^{2}I\right),$

$Y_2=X^{″}{\beta }_2+\epsilon ,\epsilon ~N\left(0,{\sigma }^{2}I\right),$

其中： $X^{\prime }=\left[1X_1~X_{14}\right]$， $X^{″}=\left[1X_1~X_{15}\right]$， ${\beta }_1,{\beta }_2\in R^{11}$。

2.2. 多重共线性诊断

如果数据设计阵存在多重共线性问题，就不能简单的用最小二乘估计对其进行参数估计，虽然看似是对数据进行了最小均方误差的拟合，但是有时可能导致拟合的结果完全和事实背离，造成数据分析无效等后果。为了使数据得到更好的处理，我们将引用岭估计、主成分回归和逐步回归等方法才能对数据进行更加合理的处理，所以进行多重共线性诊断是十分必要的，这关系到模型的合理性以及对数据的解释程度。

对Y₁、Y₂设计阵进行诊断

对矩阵 $Z_1$ 和 $Z_2$ 进行中心化和标准化：

$Z_1=\left[{\left(\frac{X_1{}_{ij}-{\mu }_j}{{\sigma }_j}\right)}_{ij}\right],Z_2=\left[{\left(\frac{X_2{}_{ij}-{\mu }_j}{{\sigma }_j}\right)}_{ij}\right],$

其中 ${\mu }_j$ 和 ${\sigma }_j$ 分别表示矩阵Z的第j列的均值和方差。计算 ${Z^{\prime }}_1{Z}_1$ 和 ${Z^{\prime }}_2{Z}_2$ 的特征值。

引入条件数公式 $k={\lambda }_1/{\lambda }_{p-1}$，k的值越大，说明存在的共线性越严重。

$k_1=\frac{{\lambda }_1}{{\lambda }_{p-1}}=\frac{1313.9}{2.9}\approx 453.81>1000,$

$k_2=\frac{{\lambda }_1}{{\lambda }_{p-1}}=\frac{1301.2}{0.9}\approx 1408.03>1000,$

从 $k_1$ 和 $k_2$ 的取值来看，不管是 $Y_1$ 还是 $Y_2$ 的设计阵，都存在略严重的多重共线性。此时如果要对 $Y_1$ 和 $Y_2$ 分别进行线性拟合，就不能再使用最小二乘估计了，处理此类问题常见的方法有常见的一般岭估计、逐步回归分析、主成分回归。

3. 理论推导

3.1. 改进广义岭估计理论推导

线性模型拟合中，对于参数的估计最常用的的是最小二乘估计，但是在遇到设计阵存在多重共线性的时候，往往最小二乘估计不是最好的，所以A.E.Hoerl提出了岭估计，引入岭参数 $k{I}_{p\times p}$ 来解决多重共线性问题，后来广义岭估计的提出又定义了岭参数 $k_1~k_p$ 可以不相同来减小拟合结果的均方误差。叶、朱(1998) [2] 提出了奇异值分解的解法。本节在前人建立的广义岭估计的基础上又提出对于若干异常值较大情况的一种广义岭估计的修正方法，改进广义岭估计。先对广义岭估计进行回顾，

岭估计参数 $\hat{\beta }$ 的估计方法：

$\hat{\beta }\left(k\right)={\left(X^{\prime }X+k{I}_{p\times p}\right)}^{-1}{X}^{\prime }Y={\left(X^{\prime }X+k{I}_{p\times p}\right)}^{-1}\left(X^{\prime }X\right)\hat{\beta },$

建立典则方程：

$Y={\alpha }_{0}1+Z\alpha +\epsilon ,$

典则回归系数的岭估计：

$\hat{\alpha }\left(k\right)={\left(\Lambda +k{I}_{p\times p}\right)}^{-1}{Z}^{\prime }Y={\left(\Lambda +kI\right)}^{-1}\Lambda \hat{\alpha },$

其中 $\hat{\alpha }={\left(\text{Λ}\right)}^{-1}{Z}^{\prime }Y$， $\Lambda ={\varphi }^{\prime }\left(X^{\prime }X\right)\varphi =\text{diag}\left({\lambda }_1,\cdots ,{\lambda }_p\right)$， $\varphi =\left({\phi }_1,\cdots ,{\phi }_p\right)$ 是 ${\lambda }_i$ 对应的特征根， $Z=X\varphi$， $\hat{\alpha }={\varphi }^{\prime }\hat{\beta }$。定义 $K=\text{diag}\left(k_1,\cdots ,k_p\right)$， $k_i$ 是不必全相等的岭参数，于是定义该模型广义岭估计的典则参数和平差参数估计如下：

$\begin{array}{l}\hat{\alpha }\left(K\right)={\left(Z^{\prime }Z+K\right)}^{-1}{Z}^{\prime }Y={\left(\Lambda +K\right)}^{-1}{Z}^{\prime }Y={\left(\Lambda +K\right)}^{-1}\Lambda \hat{\alpha }={\left(\Lambda +K\right)}^{-1}C,\\ \hat{\beta }\left(K\right)=\varphi \hat{\alpha }\left(K\right)=\varphi {\left(\Lambda +K\right)}^{-1}\Lambda {{\varphi }^{\prime }}^{\hat{\beta }}=\varphi {\left(\Lambda +K\right)}^{-1}{\varphi }^{\prime }{X}^{\prime }Y=\left(X^{\prime }X+\varphi K{\varphi }^{\prime }\right)X^{\prime }Y,\end{array}$ (A)

当 $K=k{I}_{p\times p}$ 时，广义岭估计就简化成为了一般的岭估计。

由于实际数据或许存在若干较大异常值情况，故提出修正参数对广义岭估计进行改进：

${\gamma }_{ij}=1-I\left\{\left|x_{ij}-{\bar{x}}_j\right|>Q_{ja}\right\}\left(x_{ij}-\overline{x_j}\right),$

其中 $I\left\{x_{ij}>Q_{ja}\right\}$ 为示性函数， $Q_{ja}$ 为对应设计阵X的第j列进行中心化并求绝对值后的a%分位数， $\overline{x_j}$ 为对应设计阵X的第j列的均值。提出修正参数 ${\gamma }_{ij}$，将消除由于自变量中存在若干较大异常值对模型拟合的影响。

相应的，加入修正参数后的典则方程和平差方程如下：

$\hat{\alpha }\left(K\right)={\left[{\left(Z{\gamma }_{ij}\right)}^{\prime }\left(Z{\gamma }_{ij}\right)+K\right]}^{-1}{\left(Z{\gamma }_{ij}\right)}^{\prime }Y,$

$\hat{\beta }\left(K\right)={\left[{\left(X{\gamma }_{ij}\right)}^{\prime }\left(X{\gamma }_{ij}\right)+\varphi K{\varphi }^{\prime }\right]}^{-1}{\left(X{\gamma }_{ij}\right)}^{\prime }Y,$

令 $ZR=Z{\gamma }_{ij},XR=X{\gamma }_{ij}$ 得：

$\hat{\alpha }\left(K\right)={\left[Z{R}^{\prime }ZR+K\right]}^{-1}Z{R}^{\prime }Y,$

$\hat{\beta }\left(K\right)={\left[X{R}^{\prime }XR+\varphi K{\varphi }^{\prime }\right]}^{-1}X{R}^{\prime }Y,$

易证 $\text{MSE}\left[\hat{\alpha }\left(K\right)\right]=\text{MSE}\left[\hat{\beta }\left(K\right)\right]$，这里略去证明过程，于是可以通过 $\hat{\alpha }\left(K\right)$ 来作为评价 $\hat{\beta }\left(K\right)$ 估计值好坏的指标。

$\text{MSE}\left[\hat{\alpha }\left(K\right)\right]=\text{MSE}\left[\hat{\beta }\left(K\right)\right]={\sigma }^2{\displaystyle \underset{i=1}{\overset{p-1}{\sum }}\frac{{\lambda }_i}{{\left({\lambda }_i+k_i\right)}^2}}+{\displaystyle \underset{i=1}{\overset{p-1}{\sum }}\frac{k_i^2{\alpha }_i^2}{{\left({\lambda }_i+k_i\right)}^2}}.$

3.2. 改进广义岭估计解法介绍

对于岭参数的估计，将是建立合理的线性模型的至关重要的一点，虽然前人已经提出了各种确定岭参数的方法和准则，但是到目前为止，并没有一种公认的非常稳健的方法，所以本文将采用游、王和刘(2002) [1] 提出的广义岭估计的直接解法DSGRE (derct solution to generalized ridge estimate)，其内容如下：

由于 ${\lambda }_i$ 和 ${\alpha }_i$ 与 $k_i$ 没有关系且 ${\alpha }_i$ 是非随机的，为了求得最小的 $\text{MSE}\left[\hat{\alpha }\left(K\right)\right]$，故对 $\text{MSE}\left[\hat{\alpha }\left(K\right)\right]$ 关于 $k_i$ 进行一阶求导：

$f^{\prime }\left(k_i\right)=\frac{\text{dMSE}\left[\hat{\alpha }\left(K\right)\right]}{\text{d}{k}_i}=-2{\sigma }^2{\displaystyle \underset{i=1}{\overset{p-1}{\sum }}\frac{{\lambda }_i}{{\left({\lambda }_i+k_i\right)}^3}}+2{\displaystyle \underset{i=1}{\overset{p-1}{\sum }}\frac{k_i{\lambda }_i{\alpha }_i^2}{{\left({\lambda }_i+k_i\right)}^3}},$

令 $f^{\prime }\left(k_i\right)=0$,

解得：

$k_i=\frac{{\sigma }^2}{{\alpha }_i^2},i=1,2,\cdots ,p-1,$

由于真实的 ${\sigma }^2$ 和 ${\alpha }_i^2$ 是未知的，故用 ${\sigma }^2$ 的估计值 $\widehat{{\sigma }^2}$ 以及 ${\alpha }_i$ 的迭代值 ${\widehat{{\alpha }_i}}^{\left(j\right)}$ 来对 $k_i$ 进行迭代求解：

${\widehat{k_i}}^{\left(j\right)}=\frac{\widehat{{\sigma }^2}}{{\widehat{{\alpha }_i}}^{\left(j\right)}{\widehat{{\alpha }_i}}^{\left(j\right)}},$ (B)

那么就可以用第j次 ${\widehat{{\alpha }_i}}^{\left(j\right)}$ 的迭代值来估计与之相同次数的岭参数 $k_i$，由式(A)可知，第 $j+1$ 次估计为：

${\widehat{{\alpha }_i}}^{\left(j+1\right)}=\frac{c_i}{{\lambda }_i+{\widehat{k_i}}^{\left(j\right)}},$

带入式(B)可知：

${\widehat{{\alpha }_i}}^{\left(j+1\right)}=\frac{c_i}{{\lambda }_i+\frac{\widehat{{\sigma }^2}}{{\widehat{{\alpha }_i}}^{\left(j\right)}{\widehat{{\alpha }_i}}^{\left(j\right)}}}=\frac{c_i{\widehat{{\alpha }_i}}^{\left(j\right)}{\widehat{{\alpha }_i}}^{\left(j\right)}}{{\lambda }_i{\widehat{{\alpha }_i}}^{\left(j\right)}{\widehat{{\alpha }_i}}^{\left(j\right)}+{\sigma }^2}.$ (C)

3.3. 解的存在性证明

薛、梁(2002) [3] 给出了参数的迭代算法。结合已有的证明过程，对改进广义岭估计的解的存在性予以证明。如果 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 是一致收敛的，那么经过 $j+1$ 次迭代过程，必然能够求解出 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 的数值解，故下证 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 的一致收敛性：

证明：若 $c_i>0$，根据 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 的形式知： ${\widehat{{\alpha }_i}}^{\left(0\right)}=\frac{c_i}{{\lambda }_i}>0$，由于 ${\sigma }^2>0\Rightarrow {\widehat{{\alpha }_i}}^{\left(0\right)}>{\widehat{{\alpha }_i}}^{\left(1\right)}>\cdots >{\widehat{{\alpha }_i}}^{\left(\infty \right)}>0$，故

${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 在 $j\in \left(0,\infty \right)$ 上市有界并且单调的，有Dieichlet判别法可知， ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 在 $j\in \left(0,\infty \right)$ 上是一致收敛的。

若 $c_i<0$，根据 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 的形式知： ${\widehat{{\alpha }_i}}^{\left(0\right)}=\frac{c_i}{{\lambda }_i}<0$，由于 ${\sigma }^2>0\Rightarrow {\widehat{{\alpha }_i}}^{\left(0\right)}<{\widehat{{\alpha }_i}}^{\left(1\right)}<\cdots <{\widehat{{\alpha }_i}}^{\left(\infty \right)}<0$，故 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 在

$j\in \left(0,\infty \right)$ 上是有界并且单调的，有Dieichlet判别法可知， ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 在 $j\in \left(0,\infty \right)$ 上是一致收敛的。

综上所知： ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 在 $j\in \left(0,\infty \right)$ 上是一致收敛的，迭代式 ${\widehat{{\alpha }_i}}^{\left(j+1\right)}$ 必然有解。

于是，令 ${\widehat{{\alpha }_i}}^{\left(\infty \right)}=\widehat{{\alpha }_i}$，则式(C)可以写成如下形式：

$\widehat{{\alpha }_i}=\frac{c_i{\widehat{{\alpha }_i}}^2}{{\lambda }_i{\widehat{{\alpha }_i}}^2+{\sigma }^2},$ (D)

解式(D)得：

$\widehat{{\alpha }_i}=\{\begin{array}{l}0,c_i=0或c_i^2-4{\lambda }_i{\sigma }^2\ge 0,\\ \frac{c_i+\sqrt{c_i^2-4{\lambda }_i{\sigma }^2}}{2{\lambda }_i},c_i>0,c_i^2-4{\lambda }_i{\sigma }^2\ge 0,\\ \frac{c_i+\sqrt{c_i^2-4{\lambda }_i{\sigma }^2}}{2{\lambda }_i},c_i<0,c_i^2-4{\lambda }_i{\sigma }^2\ge 0.\end{array}$

又因为 $\hat{\beta }=\varphi \hat{\alpha }$，故该方法直接给出了直接解改进广义岭估计参数的方式。本节将游、王和刘(2002) [1] 提出的DSGRE在改进广义岭估计中应用，推导过程验证，也是适用的。

4. 实例探究

基于第3节中对于改进广义岭估计的研究，本节将采用该方法对于中蒙俄贸易数据进行回归分析，在广义岭估计的基础上引入修正参数形成改进广义岭估计，训练相应的模型对总体数据进行更加精确地拟合。

4.1. 算法

① 对X的每一列进行中心化后再求绝对值得到 $A={\left(\left|x_{ij}-\overline{x_j}\right|\right)}_{ij}$ ；

② 对 $\left|x_{ij}-\overline{x_j}\right|$ 按列进行排序，得到新的排序矩阵B；

③ 找到排序矩阵每一列的a%分位数 $Q_{ja}$，这里不妨取a = 95；

④ 通过判别函数将每一列小于对应这一列a%分位数 $Q_{ja}$ 的元素重置为0，大于 $Q_{ja}$ 的元素则保留得到矩阵E；

⑤ A中每个元素取绝对值得到矩阵C，用A点除C得到矩阵D；

⑥ 所以 $ZR=Z-E.*D$， $XR=X-E.*D$，其中的.*为MATLAB中的点乘语句；

⑦ 最后将ZR和XR进行标准中心化，任用ZR和XR表示。

4.2. 改进广义岭估计参数 $\hat{\beta }$

利用上述算法对相关参数直接求解：

$c{1}_i={\left(51.0118.2147.7103.8-113.4126.745.815.122.014.6\right)}^{\text{T}},$

$c{2}_i={\left(59.2136.1142.3125.9-48.8125.937.413.2-28.535.9\right)}^{\text{T}},$

$\lambda 1_i={\left(5.044.780.6101.2196.0270.6309.3330.7623.81194.1\right)}^{\text{T}},$

$\lambda 2_i={\left(5.26.345.584.5263.6274.4320.5344.8621.71193.2\right)}^{\text{T}}.$.

最后被估的回归系数为：

$\widehat{\beta 1}=\varphi \hat{\alpha }={\left(0-0.17-0.05-0.080.021-0.250.220.05000\right)}^{\text{T}},$

$\widehat{\beta 2}=\varphi \hat{\alpha }={\left(0-0.040.17-0.080.0210.240.12-0.160.10.04000\right)}^{\text{T}}.$

从回归系数来看，对于响应变量 $Y_1$ 来说， $X_1$ 、 $X_{10}$ 、 $X_{12}$ 、 $X_{14}$ 的系数为0，与 $Y_1$ 不存在线性关系，其余变量在线性回归中对 $Y_1$ 起着决定性作用。对于响应变量 $Y_2$ 来说 $X_1$ 、 $X_{11}$ 、 $X_{13}$ 、 $X_{15}$ 与 $Y_2$ 几乎没有线性关系，其余变量决定 $Y_2$ 的线性回归。

5. 结论

本文我们对广义岭估计进行了适当的改进，加入了修正参数，这使得岭估计在当数据存在部分异常点时，具有一定的稳健性。同时我们将中蒙俄的贸易数据引入到本文中，对改进广义岭估计进行了实例探究，从结果发现，改进广义岭估计能够估出线性回归系数。我们还将其结果与现有的岭估计、主成分回归、最小二乘法进行了比较，改进广义岭估计估计的结果具有最小的均方误差。岭估计、主成分回归、逐步回归是较为成熟的方法，其估计结果我们没有在文中给出。