高维协方差矩阵估计方法的比较
李小雪;明瑞星
【摘 要】The differences between the thresholding estimation and the shrinking estimation are reported by a series of simulations,and the proper estimation is proposed within these two estimations in practice. The simulations show that if the population covariance matrix is a sparse matrix,the thresholding estimation is better than that of the shrinking estimation,and vice versa.%通过模拟比较门限估计方法和收缩估计方法之间的差异，得出2种方法在实际应用中的使用范围。由模拟结果可知，若有确切的证据表明总体协方差矩阵是稀疏矩阵，则采用门限估计方法，否则，采用稳健的收缩估计方法比较恰当。
【期刊名称】《江西师范大学学报（自然科学版）》
【年(卷),期】2015(000)006
正则化一个5 5随机矩阵
【总页数】6页(P599-604)
【关键词】高维协方差矩阵;稀疏矩阵;非稀疏矩阵;门限估计;收缩估计
【作 者】李小雪;明瑞星
【作者单位】浙江工商大学统计与数学学院，浙江杭州 310018;浙江工商大学统计与数学学院，浙江杭州 310018
【正文语种】中 文
【中图分类】F224
0 引言
设Xk=(Xk1，Xk2，…，Xkp)T，k=1，2，…，n 为来自p 元正态总体Np(μ，∑)的独立样本.以 =
表示样本协方差矩阵.众所周知，当p 很小时，S 是∑的一个“好”的估计［1］.当p 较大时，S 不是∑的一个“好”的估计［2-3］.特别地，当p 大于n 时(此时的协方差矩阵估计称为高维协方差矩阵估计)，S 是一个病态矩阵，它显然不是∑的一个好的估计.要想得到∑的好的估计，需要对S 进行正则化处理，如thresholding 方法［4-8］，tapering 方法［9-15］和banding 方
法［16-20］.事实上，上述提及的各种方法都需要对协方差矩阵的结构做出相应的假设，在因子模型的协方差矩阵估计中，文献［21］假定误差协方差矩阵是一对角矩阵，文献［22］假定误差协方差矩阵是一稀疏矩阵，但是在实际应用中，如何验证误差协方差矩阵是对角矩阵还是稀疏矩阵是一件非常困难的事情.其次，对总体协方差矩阵的这些假设一旦与事实相违背，则得到的最终估计不是一个好的估计，其均方误差(MSE)比用S 去估计∑还要大.自然地，当对协方差矩阵的认识不够充分时，会选择一种稳健的估计.收缩估计正是一种稳健的估计，它可以不依赖表示样本均值，于对协方差矩阵的任何假设.收缩估计最初由C.Stein 提出［22］，在高维情形下，O.Ledoit 等［23］建议用
作为协方差矩阵∑的估计，其中I 是单位矩阵.如果对协方差矩阵有着先验的认识(这种先验认识用矩阵T 表示)，J.Schafer 等［3］建议用T 取代(1)式中的I，即
本质上，对协方差矩阵的各种假设是基于对其的一些先验认识，从这种意义上讲，收缩估计中T 的选取也是对协方差矩阵做出某种假设的表现.相对于对协方差矩阵做出的假设而言，收缩估计中T 的选取具有很大的灵活性，具体参见文献［3］.
本文对门限估计和收缩估计进行模拟研究，为应用统计工作者如何使用门限估计和收缩估计
提出一些建议.
1 门限估计
P.J.Bickel 等［4］提出用硬门限估计方法来估计总体协方差矩阵∑=(σij)p×p，其基本原理是:设定某个门槛值θ，当 时，σij 的估计值为0;否则，σij 的估计值为sij.A.J.Rothman 等［5］在文献［4］的基础上把门限估计中的硬门限扩展到了一般门限函数，即∀θ ≥0，定义函数sθ:R →R，∀z ∈R，sθ 满足(i))对于θ，sθ(z)=0;(iii).文献［5］提出的广义门限估计包含许多门限估计，如硬门限估计、软门限估计、SCAD 门限估计和adaptive LASSO 门限估计等.这几种门限估计的原理相似，下面以软门限估计为例来简单阐述其估计原理.
门限估计方法中，最为关键的一个假定是总体协方差矩阵是稀疏矩阵，即总体协方差矩阵的许多非对角元素为0 或者接近0.在软门限估计中，门限函数sθ(z)=sign(z)，其中(·)+表示对于实数a，(a)+=max(a，0)，则有
通常把θ 称为门限参数，文献［5］采用文献［4］中提出的交叉验证法来确定θ，具体步骤:把n 个样本随机分成2 部分，第1 部分的样本个数为n1 =nn/logn 和第2 部分的样本个数为n2
=n-n1，重复这种分割N 次，令S1，v 和S2，v 分别表示当第v 次分割时由第1 部分的样本和第2 部分的样本得到的样本协方差矩阵，其中v=1，2，…，N，从而有
则总体协方差矩阵∑的软门限估计为
2 收缩估计
O.Ledoit 等［23］和J.Schafer 等［3］提出用(2)式作为总体协方差矩阵的收缩估计矩阵.在(2)式中最为关键的就是λ 的确定和矩阵T 的选取.下面先阐述如何确定λ.考虑 的均方误差(MSE)，其均方误差可以分解为
如果 是总体协方差矩阵∑的1 个无偏估计，即，则减小均方误差就等价于减小.注意到MSE是关于λ 的一个函数，不妨记为R(λ)，此外记T=(tij)p×p，这样，
将上式展开，
对R(λ)求导，得
因为样本协方差阵S 是总体协方差矩阵的1 个无偏估计，即Bias(sij)=0，所以有
由于λ 的取值范围为［0，1］，所以其极值点或者端点就是取最值的点，所以
收缩估计的另外一个关键问题是如何选取目标矩阵T.文献［3］给出了T 的6 种情形.文献［24］也给出了T 的9 种选法.其实，目标矩阵如何选取没有统一的定论，需要根据实际情况选取.一般而言，目标矩阵应该是1 个低维的矩阵，其包含较少的参数.本文仅选取3 种常见的目标矩阵进行研究，其它情形可类似地进行研究.
(i)T 为单位矩阵.记Ta=I，总体协方差矩阵的收缩估计为 =λTa+(1-λ)S.
(ii)T 为对角矩阵.记
总体协方差矩阵的收缩估计为
(iii)T 为非对角矩阵.记
其中表示样本相关系数矩阵的平均值，总体协方差矩阵的收缩估计为
3 模拟
前面分别给出了门限估计和收缩估计的理论基础，本部分通过模拟来具体的比较这2 种估计方法的优劣.具体的模拟过程如下:先从均值为零向量、协方差矩阵为∑0 的多元正态分布中产生50 个样本(即n=50)，再根据前述的理论确定门限方法的估计矩阵和收缩方法的估计矩阵和，将上述过程重复30 次.用矩阵的Frobenius范数度量估计矩阵和总体协方差矩阵∑0 之间的差异.设定∑0 具有下面3 种结构:∑0 为稀疏矩阵、非稀疏矩阵或者随机矩阵(此时∑0 的元素是随机产生的，无法确定其稀疏性).
3.1 ∑0 为稀疏矩阵
类似文献［6］中模型1 产生如下的∑0:
其中h1，…，hp 均是取自均匀分布U(1，5)的随机数.模拟结果见表1.
表1 当∑0 为稀疏矩阵时的估计效果p ‖S-∑0‖F ‖∑∧t-∑0‖F ‖∑∧sa-∑0‖F ‖∑∧sb-∑0‖F ‖∑∧sc-∑0‖F 2 0.749 706 0.746 830 1.016 662 0.768 169 1.623 049 5 2.465 439 2.049 045 2.708 103 1.733 708 2.992 748 10 4.633 239 3.009 409 4.288 424 2.407 053 3.878 976 20 9.344 986 5.033 152 7.100 850 3.524 470 4.634 618 50 25.076 830 12.471 950 16.570 070 9.95
6 072 10.471 230 100 49.126 630 22.899 820 30.979 550 24.073 720 24.239 240 200 98.360 930 43.904 870 60.218 750 53.079 700 53.164 350 300 152.928 000 67.918 060 93.124 660 85.590 600 85.640 410 400 203.239 700 90.482 560 122.596 800 115.098 800 115.135 800 500 247.420 700 103.907 100 147.899 600 140.525 300 140.551 500
若总体协方差矩阵为稀疏矩阵，则由表1 可知，当维数特别低时(相对于样本量而言，如n=2)，样本协方差矩阵是一个好的估计.当维数比较低时(p ＜n)，门限估计方法比收缩估计方法的表现要好，而收缩估计方法比样本协方差矩阵的表现要好.当维数p 大于样本量n 时，门限估计方法也比其他方法表现的要好.另外，在3 种收缩估计中，T 为单位矩阵时表现的比其他2 种收缩估计要差，而T 为对角矩阵和非对角矩阵的表现很接近，几乎相同.图1 更加直观地体现了这些结论.
图1 当∑0 为稀疏矩阵时各种估计方法的表现
由图1 可以看出，由于T 为对角矩阵和非对角矩阵的表现很接近，所以在图1 中2 条线合成1 条线，而门限估计比其他方法的表现都要好.所以，当总体协方差矩阵是稀疏矩阵时，选用门限估计方法最为合适.
3.2 ∑0 为非稀疏矩阵
首先叙述产生非稀疏矩阵的过程.从正态分布N(0，5)中产生∑0 的对角元素，从均匀分布U(2，5)中产生∑0 的非对角元素，得出∑0 后进行谱分解，其中λi 为∑0 的特征值，μi为λi 对应的特征向量.如果λi 为负值或0，则用1 来代替，即，此时，，上面产生∑0 用到的正态分布和均匀分布的参数以及λi 的替代值都可以进行修改，但是要注意一点，产生的∑0 的非对角元素的值不能太小也不能为0，这是为了保证∑0 为一非稀疏矩阵.模拟结果见表2.
若总体协方差矩阵为非稀疏矩阵，则由表2 可知，当维数特别低时(相对于样本量而言，如n=2)，样本协方差矩阵比其他估计方法表现的要好，是一个好的估计.当维数比较低时(p ＜n)，门限估计方法比样本协方差矩阵表现的差，但是差别不是太大.当维数p 大于样本量n 时，门限估计方法比样本协方差矩阵的表现差很多，且样本协方差矩阵比3 种收缩估计方法的表现差很多.这表明在总体协方差矩阵是非稀疏矩阵时，门限估计方法不再适用，而收缩估计方法比样本协方差矩阵的估计效果要好的多.具体如图2 所示.
表2 当∑0 为非稀疏矩阵时的估计效果p ‖S-∑0‖F ‖∑∧t-∑0‖F ‖∑∧sa-∑0‖F ‖∑∧sb-∑0‖F ‖∑∧sc-∑0‖F 2 0.837 686 0.907 051 0.991 598 1.142 594 1.180 436 5 3.339 797 3.664 142
3.384 116 3.500 350 3.695 759 10 8.137 748 8.281 267 7.964 004 8.013 193 9.001 249 20 13.608 750 15.352 020 13.751 460 13.796 020 13.540 260 50 42.649 590 44.765 930 42.112 170 42.145 620 40.930 300 100 85.958 390 91.417 720 84.919 140 84.946 990 76.713 390 200 180.306 900 193.087 200 176.470 900 176.464 500 158.032 900 300 334.022 400 339.773 100 324.450 500 324.439 600 281.302 100 400 432.133 200 468.268 000 425.445 800 425.440 900 332.726 900 500 588.745 800 604.403 000 569.064 500 569.001 300 449.255 400

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。