單指標分位回歸模型估計的MM算法*

郭媛媛，楊雪梅，孫志華,3?

(1 中國科學院大學數(shù)學科學學院， 北京 100049； 2 華北電力大學數(shù)理學院， 北京 102206；3 中國科學院大數(shù)據(jù)挖掘與知識管理重點實驗室， 北京 100190)(2019年5月14日收稿； 2019年10月9日收修改稿)

分位回歸模型具有穩(wěn)健的特點，并且能夠?qū)憫兞康姆植甲龀鼍毜拿枋觯?因此獲得很多學者的關(guān)注。 目前， 分位回歸方法已經(jīng)成為分析數(shù)據(jù)的一個非常重要的工具，廣泛地應用在金融、醫(yī)療、生物研究等領(lǐng)域。 分位回歸模型的研究和應用可參考文獻[1-4]及相關(guān)文獻。 參數(shù)分位回歸模型是在實際中應用廣泛的一類回歸模型[5]。但很多時候無法對參數(shù)分位回歸模型進行正確的設定?；谡`定的參數(shù)分位回歸的統(tǒng)計推斷結(jié)果經(jīng)常是不可信的， 文獻[6-7] 闡述了參數(shù)分位回歸需要進行檢驗的必要性。非參數(shù)分位回歸不存在誤定的風險，但當樣本量比較小且協(xié)變量比較多時，非參數(shù)分位回歸方法可能會受到維數(shù)禍根的問題的困擾。

對分位回歸，構(gòu)建目標函數(shù)時用到的損失函數(shù)ρτ(r)=τr-I{r<0}具有不光滑的特點, 從而使得求解目標函數(shù)的最小值比較困難，且可能出現(xiàn)多個最小值點的情況，參見文獻[5,8-11]。一種解決上面問題的方法是將分位回歸模型的求解問題轉(zhuǎn)化為線性規(guī)劃問題，再利用單純形法或內(nèi)點法進行計算。不管是單純形法，還是內(nèi)點法，運算效率都不能令人滿意。2000 年，Hunte和Lange[12]提出一種新的用于求解分位回歸問題的算法，即MM算法。MM算法概念簡單，易于執(zhí)行，且數(shù)值穩(wěn)定，比內(nèi)點法擁有更強的數(shù)值計算能力。文獻[11]對4種求解分位回歸問題的算法，即內(nèi)點法、MM算法、坐標下降法和ADMM算法進行比較研究，驗證了MM 算法具有數(shù)值穩(wěn)定和計算效率高的特點。

單指標分位回歸模型具有降維的效果，同時保持了非參數(shù)分位回歸的穩(wěn)健性，其估計問題的研究吸引了很多研究者的興趣。文獻[13-14] 提出基于兩步迭代的估計方法，文獻[15-16]進一步提出不需要迭代的估計方法，文獻[17-18]又提出基于貝葉斯方法的估計方法，文獻[19-21]探討單指標分位回歸模型的變量選擇以及加權(quán)復合單指標分位回歸模型的估計。然而這些文獻所提的估計方法都基于內(nèi)點法來實現(xiàn)，內(nèi)點法在計算分位回歸模型時，計算效率低、耗時久，尤其在樣本量較大的情況下，這種缺點更為明顯。MM 算法在求解分位回歸模型的估計時比較高效和便捷，這在文獻[22-24]中均有體現(xiàn)，但是沒有文獻研究單指標分位回歸模型的MM算法，故本文研究單指標分位回歸模型估計的MM算法。

我們借鑒文獻[12] 的方法，對單指標分位回歸模型的每一步迭代程序中目標函數(shù)構(gòu)建其替代函數(shù)，從而將復雜的優(yōu)化問題簡單化。然后，基于優(yōu)化函數(shù)再進行求解計算得到估計值。我們構(gòu)建的優(yōu)化函數(shù)是光滑的，并能夠保證每次迭代目標函數(shù)是下降的。數(shù)值模擬和實例分析結(jié)果表明基于MM 算法的估計程序具有較好的穩(wěn)定性，能夠得到比較準確的估計結(jié)果，并且相較于傳統(tǒng)的內(nèi)點算法具有更強的數(shù)值計算能力，用時更短。

1 單指標分位回歸模型的估計介紹

對于給定的分位數(shù)τ∈(0,1)，在給定x的條件下，響應變量y的τ分位數(shù)θτ(x)與協(xié)變量x之間的關(guān)系如下：

θτ(x)=g(xTγ),

其中x∈d是d維協(xié)變量，g(·)表示未知的一元聯(lián)系函數(shù)。另外γ=(γ1,…,γd)T為未知的單指標向量，滿足‖γ‖=1且γ1>0,‖·‖表示Euclidean范數(shù)。這個約束條件是為了模型的可識別性[25]，已廣泛應用在有關(guān)單指標模型的文獻中。

本文采用局部線性方法對γ和g(·)進行估計，詳細內(nèi)容可參考文獻[13]，具體算法如下：

(1)

(2)

step 4重復step 1和step 2直至收斂。

(3)

2 估計程序的MM算法

2.1 MM算法介紹

下面介紹MM算法的基本思想。假設需要最小化的目標函數(shù)為L(θ):p→，θk為第k步的迭代值。MM算法每次迭代分兩步來進行。首先，構(gòu)造目標函數(shù)的優(yōu)化函數(shù)Q(θ|θk):p×p→ 滿足

Q(θk|θk)=L(θk),

Q(θ|θk)≥L(θ)?θ.

(4)

然后, 對優(yōu)化函數(shù)Q(θ|θk)進行最小化, 得到下一步的迭代值θk+1，則有

Q(θk+1|θk)≤Q(θk|θk).

(5)

綜合式(4)和式(5)，可知L(θk+1)≤L(θk)。這種下降趨勢保證了MM算法具有顯著的數(shù)值穩(wěn)定性。

觀察目標函數(shù)式(1)～式(3)，可以發(fā)現(xiàn)它們均為非光滑函數(shù)， 因而不易得到最優(yōu)解，故借用文獻[12] 提出的MM 算法的思想來處理這個問題。其主要的處理方式如下：首先給ρτ(r)加一擾動ε，得到其近似函數(shù)

2.2 非參數(shù)部分的估計

首先，將式(1)中的目標函數(shù)修正為

(6)

(7)

由此，可以將MM算法總結(jié)為如下步驟：

1) 選擇迭代初始值θ0和一個較小的正常數(shù)ε，置k=0；

3) 令k=k+1，判斷是否滿足收斂準則，若滿足收斂準則，即： 當

時，可終止迭代，其中δ是預先取定的足夠小的數(shù)。否則返回2)繼續(xù)迭代，直到滿足收斂準則。

2.3 單指標向量的估計

首先定義式(2)中目標函數(shù)的近似函數(shù)：

在γk處的優(yōu)化函數(shù)可以構(gòu)建為

(8)

最后，可將基于MM算法的單指標模型的估計總結(jié)為如下步驟：

1)參考第1節(jié)step 1所提供的方法，得到γ的初始估計；

4)重復2)、3)步驟，直至收斂。

2.4 聯(lián)系函數(shù)g(·)的估計

對于聯(lián)系函數(shù)的估計，式(6)中目標函數(shù)式的近似函數(shù)可定義為

3 數(shù)值模擬

3.1 模擬1

借鑒文獻[13]模擬1的模型設置，考慮模型

在不同的分位數(shù)下，計算上述評價指標值，所得結(jié)果列于表1和表2，可以看出，用MM算法計算單指標分位回歸模型，無論是單指標向量的估計還是聯(lián)系函數(shù)的估計，都有比較好的結(jié)果，且在較小和較大的分位數(shù)下依然有良好的表現(xiàn)，這表明本文所提出的計算方法是有效的。將該估計結(jié)果與文獻[13] 中關(guān)于該模型的模擬結(jié)果進行對比，可以看出，用MM 算法計算單指標分位回歸模型，與內(nèi)點法相比，估計的精度是相似的，估計系數(shù)的偏差大小都在10-3～10-2，估計系數(shù)的樣本標準差數(shù)量級均為10-2，但在計算效率上，本文所提出的方法大大優(yōu)于內(nèi)點法，這將在模擬3 中展示和說明。

表1 模擬1中不同的τ下，估計的Bias、SE、MSE、CI、CPTable 1 The Bias、SE、MSE、CI、CP of under different choices of τ in simulation 1

表2 模擬1中不同τ的選擇下，的ASE、AAE的Mean、SETable 2 Outcomes of for the models under different choices of τ in simulation 1

3.2 模擬2

借鑒文獻[13]模擬2的模型設置，考慮模型

驅(qū)動層介于系統(tǒng)層和硬件設備之間，提供應用程序訪問硬件設備資源的接口，同時也為存儲軟件提供了基礎(chǔ)環(huán)境和接口。linux系統(tǒng)將設備分為3類：字符設備、塊設備、網(wǎng)絡設備。考慮到安全存儲的硬件加解密設備與塊設備、網(wǎng)絡設備特性的差異，而與面向流的字符設備類似，驅(qū)動層在采用字符設備的基礎(chǔ)上提供相應的庫文件供應用程序?qū)崿F(xiàn)加解密功能[15]。

3.3 與內(nèi)點法比較

由前兩個模擬可以發(fā)現(xiàn)，MM算法在單指標分位回歸估計的計算問題中表現(xiàn)良好，接下來比較MM算法與內(nèi)點法的計算效率，考慮如下3個模型：

表3 模擬2中不同的τ下，估計的Bias、SE、MSE、CI、CP Table 3 The Bias、SE、MSE、CI、CP of under different choices of τ in simulation 2

表4 模擬2中不同τ的選擇下，的ASE、AAE的Mean、SE Table 4 Outcomes of for the models under different choices of τ in simulation 2

表5 n=100， 模型1、2、3的估計結(jié)果比較 Table 5 Estimation comparison among models 1， 2， and 3 with n=100

可以發(fā)現(xiàn)MM算法所用的時間遠遠少于內(nèi)點法，且隨著樣本量的增大，這種計算效率上的優(yōu)勢更加明顯。這是由于用內(nèi)點法解決分位回歸問題，是將目標函數(shù)及約束條件轉(zhuǎn)化為線性規(guī)劃問題，再用內(nèi)點法來求解該問題，但轉(zhuǎn)化之后的線性規(guī)劃問題，協(xié)變量維數(shù)與樣本量的大小有著正相關(guān)的關(guān)系，這種方法極大地增加了算法的計算量與所用時間。本文第3.2節(jié) 中的問題轉(zhuǎn)化為線性規(guī)劃后，協(xié)變量的維數(shù)為2n+p，第3.3 節(jié) 中的問題轉(zhuǎn)化為線性規(guī)劃后，維數(shù)為2n2+p，具體轉(zhuǎn)化方法及維數(shù)的增加量可見文獻[11]。而MM 算法只需對p維矩陣做運算，故兩種方法的計算效率隨樣本量的增加會產(chǎn)生越來越大的差距。

表6 n=200，模型1、2、3的估計結(jié)果比較Table 6 Estimation comparison among models 1， 2， and 3 with n=200

4 實例分析

θτ(MEDV|RM,TAX,PTRATIO,LSTAT)=

g(γ1RM+γ2log(TAX)+γ3PTRATIO)+

γ4log(LSTAT).

用本文提出的方法對該問題進行估計，計算在不同分位數(shù)下系數(shù)的估計值，并采用bootstrap方法估計標準差，方法如下，具體細節(jié)可參考文獻[13,29]。

重復模擬100次計算標準差，所得結(jié)果列于表7。從表7可以發(fā)現(xiàn)，RM的系數(shù)在不同的分位數(shù)下皆為正，這表明每棟房屋的房間數(shù)量越多，房價就越高且收入越多的家庭更加在意每棟房屋的房間數(shù)量；log(TAX)的系數(shù)為負且隨分位數(shù)逐漸變大，這表明不動產(chǎn)的稅率越高，房價越低且收入較低的家庭更加在意不動產(chǎn)稅率的大??；PTRATIO 的系數(shù)為負且隨分位數(shù)變化較小，這表明學生與教師的比例越大，即教師資源越匱乏，房價越低且低收入家庭與高收入家庭對教育的重視程度是同樣大的；log(LSTAT) 系數(shù)為負且隨分位數(shù)逐漸變小，這表明一個地區(qū)低收入人群所占的百分比越高，房價越低且收入較高的家庭更加在意一個地區(qū)的低收入人群比例。

表7 波士頓數(shù)據(jù)集在單指標分位回歸模型下的系數(shù)估計及標準差估計Table 7 Coefficient estimation and standard deviationestimation of Boston data set under the single-indexquantile regression mode

圖1 聯(lián)系函數(shù)g(u)及在不同的τ下,的估計Fig.1 The link function g(u) and the estimation of under different choices of τ

5 結(jié)論

本文研究單指標分位回歸模型估計方法的MM算法。相比于內(nèi)點法，MM算法極大地縮短了計算時間，提高了運算效率。此外，本文給出單指標分位回歸模型在MM 算法下的參數(shù)估計公式，在每次迭代過程中，將協(xié)變量與響應變量的觀測值直接代入公式，即可得到參數(shù)的估計值，避免了每次迭代都要優(yōu)化目標函數(shù)的麻煩。