CVX軟件包在統(tǒng)計實驗教學中的應用

丁先文 陳雪平 陳建東 唐安民

摘 要：回歸分析是高校統(tǒng)計學的專業(yè)必修課，關于模型的變量選擇又是該門課程的重點內(nèi)容。傳統(tǒng)的變量選擇方法具有很大的局限性。文章基于CVX凸優(yōu)化包，給出了線性回歸模型、分位數(shù)回歸模型和復合分位數(shù)回歸模型中變量選擇的算法。通過模擬計算說明了該算法的可行性和有效性。

關鍵詞：CVX； 變量選擇； 教學

中圖分類號：O212.2 文獻標識碼：A 文章編號：2095-7394（2017）02-0093-05

目前，許多開設統(tǒng)計學專業(yè)的高校都將模型的回歸分析設為專業(yè)必修課，體現(xiàn)了該門課程在統(tǒng)計學中的重要地位。在該門課程的教學中，關于模型的變量選擇問題是重點內(nèi)容?，F(xiàn)有的大部分教材都是介紹傳統(tǒng)的變量選擇方法，如向前法、向后法和逐步回歸等。這些方法在回歸分析中扮演著重要角色，然而，隨著大數(shù)據(jù)時代的來臨，在海量數(shù)據(jù)下，如何快速高效地進行變量選擇面臨著巨大挑戰(zhàn)。筆者結合自身的教學實踐，探索將目前流行的一些方法應用于具體的教學過程中。

近年來，關于模型的變量選擇問題成為了統(tǒng)計學的熱點研究課題。特別是隨著大數(shù)據(jù)時代的來臨，如何高效地處理和分析大數(shù)據(jù)對現(xiàn)有的統(tǒng)計方法提出了巨大的挑戰(zhàn)。在一些實際問題中，雖然在一段時間內(nèi)可以收集到海量數(shù)據(jù)，但并不是每一個變量都對興趣變量都有顯著影響，這就需要在建立模型時剔除一些與興趣變量無關的變量，然后再進行統(tǒng)計分析，這正是統(tǒng)計學中的變量選擇問題。采用傳統(tǒng)的變量選擇方法，需要分兩步進行，首先要選擇有顯著影響的變量，其次再對模型進行統(tǒng)計分析。這類方法在大數(shù)據(jù)背景下很難實現(xiàn)，計算的效率也將受

到很大損失。Tibshirani[1]提出了一種壓縮估計方法（LASSO），該方法的的一個顯著優(yōu)點就是可以將變量選擇和參數(shù)估計同時進行，從而提高了計算效率。Fan and Li[2]針對懲罰函數(shù)提出了SCAD懲罰方法，并給出了估計量的Oracle性質(zhì)。同時Fan and Li[2]指出，一個好的估計量應該具備Oracle性質(zhì)，并說明了LASSO方法不具有Oracle性質(zhì)。Zou[3]提出了自適應LASSO的變量選擇方法，并證明了自適應LASSO方法具有Oracle性質(zhì)。關于變量選擇的詳細介紹和研究進展，請參見王大榮和張忠占[4]。

在實施變量選擇的過程中，由于目標函數(shù)或懲罰項的非光滑性，這給統(tǒng)計優(yōu)化帶來了極大的挑戰(zhàn)。Fan and Li[2]提出了局部二次近似方法來優(yōu)化目標函數(shù)，該方法依賴于初始值的選取且與閥值的選取較為敏感。Efron[5]針對線性回歸模型提出了最小角回歸算法，該方法的優(yōu)點是收斂速度快且效果很好，該算法可以通過調(diào)用R中程序包來實現(xiàn)。但是該方法需要有一定的編程基礎才能實現(xiàn)，這給教師的教學帶來了一定的難度。目前，還沒有一種通用的算法可以實現(xiàn)不同模型的變量選擇問題，本文利用Matlab中的CVX軟件包給出常見模型的變量選擇的一般算法。

CVX（凸優(yōu)化）是由Grant and Boyd[6] 基于Matlab軟件編寫的求解凸優(yōu)化問題的軟件包。該軟件包采用的是一種規(guī)則化的編程語言來描述數(shù)學優(yōu)化問題，與以往的優(yōu)化軟件包相比，它具有可讀性和易用性等特點，教師在教學過程中可通過演示法讓學生掌握該軟件包的代碼編寫規(guī)則。牛佳[7]研究了基于CVX和非負矩陣分解的圖像融合的問題；王芳， 陳勇， 葉志清等[8]研究了基于CVX工具箱的自適應波束形成實驗。然而，基于CVX對模型的變量選擇算法很少有學者研究。本文對線性回歸模型、分位數(shù)回歸模型和復合分位數(shù)回歸模型給出基于CVX的變量選擇算法。對其它的常見模型的變量選擇可以作類似的推廣應用。本文的方法可供統(tǒng)計學專業(yè)的教師在回歸分析教學中借鑒使用。

1 線性回歸模型的變量選擇算法

考慮下面的線性回歸模型

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

其中[Yi]與[Xi]分別表示響應變量及[p]維協(xié)變量，[β]是[p]維的回歸系數(shù)，[εi]為獨立同分布的隨機誤差項。假設模型（1）具有稀疏性，即參數(shù)[β]中的某些分量為0。參數(shù)[β]的最小二乘估計可以通過優(yōu)化式（2）得到：

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

由（2）式得到的參數(shù)[β]的估計具有很多優(yōu)良的性質(zhì)[9]。然而，當模型中存在稀疏性時，由（2）式得到的參數(shù)估計結果往往不能將[β]中的不顯著的分量估計為0，從而降低了估計的有效性。一個常用的辦法是采用壓縮估計法，即參數(shù)[β]的估計可通過優(yōu)化式（3）得到：

[β=argminβi=1n（Yi-XTiβ）2+nj=1ppλ（βj），] （3）其中[nj=1ppλ（βj）]稱為懲罰項，[pλ（.）]是懲罰函數(shù)，參數(shù)[λ]是調(diào)諧參數(shù)。通過選取不同的[λ]來調(diào)整懲罰程度的大小，從而達到壓縮估計的目的。當[pλ（βj）=λβj]時，式（3）即為LASSO估計；當[pλ（βj）=λωjβj]時，式（3）即為自適應LASSO估計，特別地，若[ωj=1，j=1，…，p]，則自適應LASSO估計即為LASSO估計；當懲罰函數(shù)的導數(shù)滿足

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

時，其中[α>0，θ>0]，式（3）即為SCAD估計。

注意到，式（3）的第二項在原點不可導，普通的通過梯度法尋求（3）式的最優(yōu)值不可行。然而利用關系式[βj=β+j+β-j]，[βj=β+j-β-j]，其中[β+j=βI（β>0）]和[β-j=βI（β<0）]，可以將式（3）轉化為凸線性規(guī)劃問題來解決。以下以懲罰項為自適應LASSO為例，給出基于CVX的優(yōu)化式（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表示調(diào)諧參數(shù)[λ]。對于懲罰函數(shù)為SCAD情形，也可類似運用以上代碼進行變量選擇，這時需要對SCAD懲罰函數(shù)采用一步近似方法。

2 分位數(shù)回歸模型的變量選擇算法

作為對普通最小二乘方法的一種替代方法，Koenker and Bassett （1978） 提出了分位數(shù)回歸模型。通過估計不同的條件分位數(shù)函數(shù)，分位數(shù)回歸可以系統(tǒng)地刻畫協(xié)變量對響應分布的影響。此外，分位數(shù)回歸模型對誤差分布不作任何假設，這使得分位數(shù)回歸模型得到了許多研究者的深入研究并在各領域得到了廣泛應用。關于分位數(shù)回歸模型的研究進展和詳細介紹，請參見 Koenker[10]。

考慮下面的線性回歸模型

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

其中[Yi]與[Xi]分別表示響應變量及[p]維協(xié)變量，[β]是[p]維的回歸系數(shù)，[εi]為具有未知分布函數(shù)的隨機誤差項。在給定[Xi]的條件下，令[Yi]的[τ]條件分位數(shù)為[Qτ（Yi][Xi）=XTiβτ]且滿足[P（YiXTiβτXi）=τ，]其中[0<τ<1]。當模型（4）中存在稀疏性時，可通過優(yōu)化（5）式得到參數(shù)的估計

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

其中[ρτ（t）=（τ-I（t0））]為檢查函數(shù)，[I（.）]為示性函數(shù)。由于式（5）中的兩項在原點均不可導，因此無法通過普通的梯度方法來優(yōu)化。注意到檢查函數(shù)[pτ（t）]滿足[pτ（t）=τt++（1-τ）t-]，其中[t+=tI（t>0）]，[t-=tI（t<0）]，t=[t++t-]。可以將式（5）轉化為凸線性規(guī)劃問題來解決。具體地，以懲罰項為自適應LASSO為例，優(yōu)化式（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，]。由此可以得到參數(shù)[β]的估計[βr=η+-η-]。下面給出基于CVX的優(yōu)化式（5）的執(zhí)行代碼。

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 復合分位數(shù)回歸模型的變量選擇算法

分位數(shù)估計只考慮了在某個給定的分位點上的估計，這可能對許多可能感興趣的分布無效。Zou and Yuan[11]提出了復合分位數(shù)回歸模型，其思想是通過極小化來自不同分位數(shù)回歸模型中的目標函數(shù)的一個混合結構，是一種穩(wěn)健的統(tǒng)計方法?；趶秃戏治粩?shù)回歸模型進行變量選擇會產(chǎn)生穩(wěn)健的結果。

考慮下面的線性回歸模型

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

其中[Yi]與[Xi]分別表示響應變量及[p]維協(xié)變量，[β]是[p]維的回歸系數(shù)，[εi]為具有未知分布函數(shù)的隨機誤差項。假設有K個分位點[τk，k=1，…，K]，則模型（6）中的參數(shù)估計可以通過優(yōu)化下面的復合分位數(shù)目標函數(shù)得到

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

其中[0<τk<1]是給定的K個分位點。若模型（6）中存在稀疏性，可通過優(yōu)化（7）式得到參數(shù)[β]的估計 [ β=argminβ{k=1Ki=1npτk（Yi-XTiβ-bτk）+nj=1ppλ（βj）}，]（7）

其中[pr（t）=t（τ-I（t0））]為檢查函數(shù)，[I（.）]為示性函數(shù)。利用類似于式（5）的方法，可以將（7）式轉化為線性規(guī)劃問題

[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，]。由此可以得到參數(shù)[β]的估計[βr=η+-η-]。下面給出基于CVX的優(yōu)化式（7）的執(zhí)行代碼。

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表示事先給定的分位數(shù)序列，其他符號的含義可參見優(yōu)化式（3）的代碼。

4 模擬計算

為實施模擬，本文從以下模型中產(chǎn)生數(shù)據(jù)

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

其中[β=（1，2，3，0，0，0，0，0）T]為待估參數(shù)向量，對應的[Xi]的每一個分量都獨立同分布于標準正態(tài)分布[N（0，1）]，[Yi]根據(jù)模型產(chǎn)生，模型誤差服從以下分布：M1：標準正態(tài)分布[N（0，1）]；M2：自由度為3的t分布[t（3）]；M3：混合正態(tài)分布[0.1N（0，1）+0.9N（0，10）]；M4：混合拉普拉斯分布[0.1Lap（0，1）+0.9Lap（0，10）]。為了便于比較，分位數(shù)回歸模型中取分位點為[τ=0.5]。復合分位數(shù)回歸中從區(qū)間[0.1，0.9]上均勻選取9點分位點。

在模擬計算中，調(diào)諧參數(shù)根據(jù)BIC準則選取。將模擬實驗重復進行1 000次，結果如表1所示。表1中LSE表示基于最小二乘方法得到的結果，QR表示基于分位數(shù)回歸得到的結果，CQR表示基于復合分位數(shù)得到的結果?！癈”表示在1 000次模擬試驗中，回歸系數(shù)中5個為0的系數(shù)估計為0的平均個數(shù)，“I”表示在1 000次模擬試驗中，回歸系數(shù)中三個非零系數(shù)估計為0的平均個數(shù)。GMSE（廣義均方誤差）根據(jù)以下公式計算

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

通過比較GMSE的大小可以判斷參數(shù)估計的好壞。

從表1可以看出：三種方法的計算結果都較好，能夠很好地對模型進行變量選擇，這說明文中給出的基于CVX的變量選擇算法是有效的。

5 結語

本文基于CVX軟件包對線性回歸模型、分位數(shù)回歸模型和復合分位數(shù)回歸模型的變量選擇算法進行了探討，給出了Matlab代碼，解決了一類回歸模型中的變量選擇算法問題。此方法可以推廣到更多的統(tǒng)計模型，這需要在以后的教學中繼續(xù)完善和推廣，也可為回歸分析的教學提供參考。

參考文獻：

[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]. 數(shù)理統(tǒng)計與管理， 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

責任編輯 祁秀春