CVX软件包在统计实验教学中的应用

丁先文 陈雪平 陈建东 唐安民

摘 要：回归分析是高校统计学的专业必修课，关于模型的变量选择又是该门课程的重点内容。传统的变量选择方法具有很大的局限性。文章基于CVX凸优化包，给出了线性回归模型、分位数回归模型和复合分位数回归模型中变量选择的算法。通过模拟计算说明了该算法的可行性和有效性。

关键词：CVX； 变量选择； 教学

中图分类号：O212.2 文献标识码：A 文章编号：2095-7394（2017）02-0093-05

目前，许多开设统计学专业的高校都将模型的回归分析设为专业必修课，体现了该门课程在统计学中的重要地位。在该门课程的教学中，关于模型的变量选择问题是重点内容。现有的大部分教材都是介绍传统的变量选择方法，如向前法、向后法和逐步回归等。这些方法在回归分析中扮演着重要角色，然而，随着大数据时代的来临，在海量数据下，如何快速高效地进行变量选择面临着巨大挑战。笔者结合自身的教学实践，探索将目前流行的一些方法应用于具体的教学过程中。

近年来，关于模型的变量选择问题成为了统计学的热点研究课题。特别是随着大数据时代的来临，如何高效地处理和分析大数据对现有的统计方法提出了巨大的挑战。在一些实际问题中，虽然在一段时间内可以收集到海量数据，但并不是每一个变量都对兴趣变量都有显著影响，这就需要在建立模型時剔除一些与兴趣变量无关的变量，然后再进行统计分析，这正是统计学中的变量选择问题。采用传统的变量选择方法，需要分两步进行，首先要选择有显著影响的变量，其次再对模型进行统计分析。这类方法在大数据背景下很难实现，计算的效率也将受

到很大损失。Tibshirani[1]提出了一种压缩估计方法（LASSO），该方法的的一个显著优点就是可以将变量选择和参数估计同时进行，从而提高了计算效率。Fan and Li[2]针对惩罚函数提出了SCAD惩罚方法，并给出了估计量的Oracle性质。同时Fan and Li[2]指出，一个好的估计量应该具备Oracle性质，并说明了LASSO方法不具有Oracle性质。Zou[3]提出了自适应LASSO的变量选择方法，并证明了自适应LASSO方法具有Oracle性质。关于变量选择的详细介绍和研究进展，请参见王大荣和张忠占[4]。

在实施变量选择的过程中，由于目标函数或惩罚项的非光滑性，这给统计优化带来了极大的挑战。Fan and Li[2]提出了局部二次近似方法来优化目标函数，该方法依赖于初始值的选取且与阀值的选取较为敏感。Efron[5]针对线性回归模型提出了最小角回归算法，该方法的优点是收敛速度快且效果很好，该算法可以通过调用R中程序包来实现。但是该方法需要有一定的编程基础才能实现，这给教师的教学带来了一定的难度。目前，还没有一种通用的算法可以实现不同模型的变量选择问题，本文利用Matlab中的CVX软件包给出常见模型的变量选择的一般算法。

CVX（凸优化）是由Grant and Boyd[6] 基于Matlab软件编写的求解凸优化问题的软件包。该软件包采用的是一种规则化的编程语言来描述数学优化问题，与以往的优化软件包相比，它具有可读性和易用性等特点，教师在教学过程中可通过演示法让学生掌握该软件包的代码编写规则。牛佳[7]研究了基于CVX和非负矩阵分解的图像融合的问题；王芳， 陈勇， 叶志清等[8]研究了基于CVX工具箱的自适应波束形成实验。然而，基于CVX对模型的变量选择算法很少有学者研究。本文对线性回归模型、分位数回归模型和复合分位数回归模型给出基于CVX的变量选择算法。对其它的常见模型的变量选择可以作类似的推广应用。本文的方法可供统计学专业的教师在回归分析教学中借鉴使用。

1 线性回归模型的变量选择算法

考虑下面的线性回归模型

[Yi=XTiβ+εi，i=1，…，n，] （1）

其中[Yi]与[Xi]分别表示响应变量及[p]维协变量，[β]是[p]维的回归系数，[εi]为独立同分布的随机误差项。假设模型（1）具有稀疏性，即参数[β]中的某些分量为0。参数[β]的最小二乘估计可以通过优化式（2）得到：

[β=argminβi=1n（Yi-XTiβ）2。] （2）

由（2）式得到的参数[β]的估计具有很多优良的性质[9]。然而，当模型中存在稀疏性时，由（2）式得到的参数估计结果往往不能将[β]中的不显著的分量估计为0，从而降低了估计的有效性。一个常用的办法是采用压缩估计法，即参数[β]的估计可通过优化式（3）得到：

[β=argminβi=1n（Yi-XTiβ）2+nj=1ppλ（βj），] （3）其中[nj=1ppλ（βj）]称为惩罚项，[pλ（.）]是惩罚函数，参数[λ]是调谐参数。通过选取不同的[λ]来调整惩罚程度的大小，从而达到压缩估计的目的。当[pλ（βj）=λβj]时，式（3）即为LASSO估计；当[pλ（βj）=λωjβj]时，式（3）即为自适应LASSO估计，特别地，若[ωj=1，j=1，…，p]，则自适应LASSO估计即为LASSO估计；当惩罚函数的导数满足

[p'λ（θ）=λ（I（θλ）+（αλ-θ）+（α-1）λI（θ>λ））]

时，其中[α>0，θ>0]，式（3）即为SCAD估计。

注意到，式（3）的第二项在原点不可导，普通的通过梯度法寻求（3）式的最优值不可行。然而利用关系式[βj=β+j+β-j]，[βj=β+j-β-j]，其中[β+j=βI（β>0）]和[β-j=βI（β<0）]，可以将式（3）转化为凸线性规划问题来解决。以下以惩罚项为自适应LASSO为例，给出基于CVX的优化式（3）的代码。

cvx_begin quiet

variable s（p）

variable t（p）

minimize（（y-x?（s-t））?（y-x?（s-t））+ n?lambda?weight?（s+t））

subject to

s>=0；

t>=0；

cvx_end

在以上代码中，y为n维的响应变量，[X]为[n×p]的设计矩阵，weight表示自适应权重[ω=（ω1，…，ωp）T]，在计算时可令[ωj=（β0j-2]，s表示[β+j]，t表示[β-j]，lambda表示调谐参数[λ]。对于惩罚函数为SCAD情形，也可类似运用以上代码进行变量选择，这时需要对SCAD惩罚函数采用一步近似方法。

2 分位数回归模型的变量选择算法

作为对普通最小二乘方法的一种替代方法，Koenker and Bassett （1978） 提出了分位数回归模型。通过估计不同的条件分位数函数，分位数回归可以系统地刻画协变量对响应分布的影响。此外，分位数回归模型对误差分布不作任何假设，这使得分位数回归模型得到了许多研究者的深入研究并在各领域得到了广泛应用。关于分位数回归模型的研究进展和详细介绍，请参见 Koenker[10]。

考虑下面的线性回归模型

[Yi=XTiβ+εi，i=1，…，n，] （4）

其中[Yi]與[Xi]分别表示响应变量及[p]维协变量，[β]是[p]维的回归系数，[εi]为具有未知分布函数的随机误差项。在给定[Xi]的条件下，令[Yi]的[τ]条件分位数为[Qτ（Yi][Xi）=XTiβτ]且满足[P（YiXTiβτXi）=τ，]其中[0<τ<1]。当模型（4）中存在稀疏性时，可通过优化（5）式得到参数的估计

[βτ=argminβ{i=1nρτ（Yi-XTiβ）+nj=1ppλ（βj）}，] （5）

其中[ρτ（t）=（τ-I（t0））]为检查函数，[I（.）]为示性函数。由于式（5）中的两项在原点均不可导，因此无法通过普通的梯度方法来优化。注意到检查函数[pτ（t）]满足[pτ（t）=τt++（1-τ）t-]，其中[t+=tI（t>0）]，[t-=tI（t<0）]，t=[t++t-]。可以将式（5）转化为凸线性规划问题来解决。具体地，以惩罚项为自适应LASSO为例，优化式（5）等价于

[mint+i，t+i，η+i，η+i{i=1nτt+i+（1-τ）t-i+nλj=1pωj（η+i+η-i）}，]

满足的约束条件为：

[t+i-t-i=Yi-XTi（η+-η-）；t+i0；t-i0；η+j0；η-j0；i=1，…，n；j=1，…，p，]

其中[η+=（η+1，…，η+p）T，η-=（η-1，…，η-p）T，]。由此可以得到参数[β]的估计[βr=η+-η-]。下面给出基于CVX的优化式（5）的执行代码。

cvx_begin quiet

variable t1（n）

variable t2（n）

variable eta1（p）

variable eta2（p）

minimize（sum（tau?s+（1-tau）?t）+n?lamb da?weight?（eta1+eta2））

subject to

t1-t2==y-x?（eta1-eta2）；

t1>=0；t2>=0；eta1>=0；eta2>=0；

cvx_end

3 复合分位数回归模型的变量选择算法

分位数估计只考虑了在某个给定的分位点上的估计，这可能对许多可能感兴趣的分布无效。Zou and Yuan[11]提出了复合分位数回归模型，其思想是通过极小化来自不同分位数回归模型中的目标函数的一个混合结构，是一种稳健的统计方法。基于复合分位数回归模型进行变量选择会产生稳健的结果。

考虑下面的线性回归模型

[Yi=XTiβ+εi，i=1，…，n，] （6）

其中[Yi]与[Xi]分别表示响应变量及[p]维协变量，[β]是[p]维的回归系数，[εi]为具有未知分布函数的随机误差项。假设有K个分位点[τk，k=1，…，K]，则模型（6）中的参数估计可以通过优化下面的复合分位数目标函数得到

[βargminβ{k=1Ki=1npτk（Yi-XTiβ-bτk）}，]

其中[0<τk<1]是给定的K个分位点。若模型（6）中存在稀疏性，可通过优化（7）式得到参数[β]的估计 [ β=argminβ{k=1Ki=1npτk（Yi-XTiβ-bτk）+nj=1ppλ（βj）}，]（7）

其中[pr（t）=t（τ-I（t0））]为检查函数，[I（.）]为示性函数。利用类似于式（5）的方法，可以将（7）式转化为线性规划问题

[mint+ik，t+ik，η+i，η+i{k=1Ki=1nτkt+ik+（1-τk）t-ik+nλj=1pωj（η+i+η-i）}，]

满足的约束条件为：

[t+ik-t-ik=Yi-XTi（η+-η-）-bτk；t+ik0；t-ik0；η+j0；η-j0；i=1，…，n；j=1，…，p；k=1，…，K，]

其中[η+=（η+1，…，η+p）T，η-=（η-1，…，η-p）T，]。由此可以得到参数[β]的估计[βr=η+-η-]。下面给出基于CVX的优化式（7）的执行代码。

cvx_begin quiet

variable t1（n，K）

variable t2（n，K）

variable eta1（p）

variable eta2（p）

variable btau（K）

minimize（sum（sum（（repmat（tauseq，n，1））. ?t1+（repmat（1-tauseq，n，1））.?t2））+n?lamb da?weight'?（eta1+eta2））

subject to

t1-t2==repmat（y-x?（eta1-eta2），1，K）-rep mat（btau，n，1）；

t1>=0；t2>=0；eta1>=0；eta2>=0；

cvx_end

在上述代码中，tauseq表示事先给定的分位数序列，其他符号的含义可参见优化式（3）的代码。

4 模拟计算

为实施模拟，本文从以下模型中产生数据

[Yi=XTiβ+εi，i=1，…，100，]

其中[β=（1，2，3，0，0，0，0，0）T]为待估参数向量，对应的[Xi]的每一个分量都独立同分布于标准正态分布[N（0，1）]，[Yi]根据模型产生，模型误差服从以下分布：M1：标准正态分布[N（0，1）]；M2：自由度为3的t分布[t（3）]；M3：混合正态分布[0.1N（0，1）+0.9N（0，10）]；M4：混合拉普拉斯分布[0.1Lap（0，1）+0.9Lap（0，10）]。为了便于比较，分位数回归模型中取分位点为[τ=0.5]。复合分位数回归中从区间[0.1，0.9]上均匀选取9点分位点。

在模拟计算中，调谐参数根据BIC准则选取。将模拟实验重复进行1 000次，结果如表1所示。表1中LSE表示基于最小二乘方法得到的结果，QR表示基于分位數回归得到的结果，CQR表示基于复合分位数得到的结果。“C”表示在1 000次模拟试验中，回归系数中5个为0的系数估计为0的平均个数，“I”表示在1 000次模拟试验中，回归系数中三个非零系数估计为0的平均个数。GMSE（广义均方误差）根据以下公式计算

[ GMSE（β）=（β-β）TE（XXT）（β-β）]。

通过比较GMSE的大小可以判断参数估计的好坏。

从表1可以看出：三种方法的计算结果都较好，能够很好地对模型进行变量选择，这说明文中给出的基于CVX的变量选择算法是有效的。

5 结语

本文基于CVX软件包对线性回归模型、分位数回归模型和复合分位数回归模型的变量选择算法进行了探讨，给出了Matlab代码，解决了一类回归模型中的变量选择算法问题。此方法可以推广到更多的统计模型，这需要在以后的教学中继续完善和推广，也可为回归分析的教学提供参考。

参考文献：

[1] TIBSHIRANI R. Regression Shrinkage and Selection via the Lasso：a retrospective[J]. Journal of the Royal Statistical Society， 1994， 58（1）：267-288.

[2] FAN J， LI R. Variable selection via nonconvave penalized likelihood and its oracle properties[J].Journal of the American Statistical Association， 2001， 96（456）：1 348-1 360.

[3] ZOU H. The Adaptive Lasso and Its Oracle Properties[J]. Journal of the American Statistical Association， 2006， 101（476）：1 418-1 429.

[4] 王大荣， 张忠占. 线性回归模型中变量选择方法综述[J]. 数理统计与管理， 2010， 29（4）：615-627.

[5] EFRON B，HASTIE T. Least angle regression[J]. Mathematics， 2004， 32（2）：407-451.

[6] GRANT M， BOYD S P. CVX： MATLAB software for disciplined convex programming[J]. Global Optimization， 2014：155-210.

[7] 牛佳. 基于CVX和非负矩阵分解的图像融合研究[J]. 计算机工程与设计， 2008， 29（20）：5 311-5 313.

[8] 王芳， 陈勇， 叶志清，等. 基于CVX工具箱的自适应波束形成实验[J]. 电气电子教学学报， 2016， 38（2）：136-139.

[9] 唐年胜， 李会琼. 应用回归分析[M]. 北京：科学出版社， 2014.

[10] KOENKER R. Quantile regression[M]. Cambridge Massachusetts：Cambridge university press， 2005.

[11] ZOU H， YUAN M. Composite quantile regression and the Oracle model selection theory [J]. The Annals of Statistics， 2008，36（3）：1 108-1 126.

Application of CVX Software Package in Statistical Experiment Teaching

DING Xian-wen1，CHEN Xue-ping1 ， CHEN Jian-dong1， TANG An-min2

（1.School of Mathematics and Physics， Jiangsu University of Technology， Changzhou 213001， China；

2.Department of Statistics， Yunnan University， Kunming 65000， China）

Abstract： Regression analysis is a compulsory subject of statistics in college and the variable selection of model is the key content of this course. The traditional variable selection method has a lot of limitations. Based on the software package of CVX in Matlab， we propose an optimization algorithm of variable selection in linear regression model， quantile regression model and composite quantile regression model. The simulation study presents the feasibility and validity of the proposed algorithm.

Key words： CVX； variable selection； teaching

责任编辑 祁秀春