姜云盧，袁晶晶

(1．中山大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院，廣東廣州510275;2．暨南大學(xué)經(jīng)濟(jì)學(xué)院，廣東廣州510632)

對于線性回歸模型，最小二乘方法是估計回歸系數(shù)的一種最常見的方法，描述了自變量對于因變量的均值影響．如果線性回歸模型中的隨機(jī)誤差項服從均值為零且同方差的正態(tài)分布，則回歸系數(shù)的最小二乘估計為極大似然估計．然而，在實際應(yīng)用中，上述的正態(tài)假設(shè)經(jīng)常不能被滿足，例如數(shù)據(jù)經(jīng)常出現(xiàn)尖峰或重尾的分布，或存在顯著的異方差等情形．而且，最小二乘估計對離群點特別敏感，一個“壞”的數(shù)據(jù)點就可能破壞整個回歸線的估計．

1978 年，KOENKER 和 BASSETT［1］首次提出了分位數(shù)回歸(Quantile Regression)的思想．假設(shè)隨機(jī)變量Y具有分布函數(shù)為F(y)=P(Y≤y)，則Y的τ(0≤τ≤1)分位數(shù)定義為:Q(τ)=inf{y:F(y)≥τ};特別地，Q(0．5)就是中位數(shù)．指出了隨機(jī)變量Y的τ分位數(shù)可定義為

其中 ρτ(r)=τr-rΙ(r＜0)稱為檢驗函數(shù)(check function)．進(jìn)一步，他們將檢驗函數(shù)推廣到分位數(shù)回歸的情形．分位數(shù)回歸的假設(shè)與均值回歸不同，對誤差項方差齊性與正態(tài)性沒有要求，同時善于捕捉分布的尾部特征．注意到把檢驗函數(shù)替換成ρτ(r)=r2即得到傳統(tǒng)的均值回歸估計，對于異常值其二次增長明顯比一次增長受影響更大，所以在穩(wěn)健性方面分位數(shù)回歸優(yōu)于均值回歸．另外，分位數(shù)回歸能給出響應(yīng)變量Y更為全面的描述，τ取不同的值可以反映解釋變量X對不同水平的Y的影響．還討論了分位數(shù)回歸估計的一致性與漸進(jìn)正態(tài)性．

有關(guān)分位數(shù)回歸估計的計算，1993年，KOENKER和OREY［2］把線性規(guī)劃中常用的單純形法擴(kuò)展到分位數(shù)回歸中，該算法所得到的估計雖然具有很好的穩(wěn)定性，但是，在處理大型數(shù)據(jù)時效率低下．1996年，KOENKER和PARK［3］用線性規(guī)劃的內(nèi)點算法得到了非線性分位數(shù)回歸的估計．雖然在大型數(shù)據(jù)下比單純形算法效率要高，但是，當(dāng)自變量比較多時效率仍然很低．2000年，HUNTER 和 LANGE［4］首次提出MM(Majorization-Minimization)算法，它概念簡單，數(shù)值穩(wěn)定，且易于執(zhí)行，并將其用于非線性分位數(shù)回歸中．他們的數(shù)值結(jié)果表明，在許多非線性分位數(shù)回歸問題中，MM-算法優(yōu)于內(nèi)點方法．

分位數(shù)回歸的理論研究和實際應(yīng)用發(fā)展迅速，其中理論已經(jīng)拓展到二元響應(yīng)模型、非參數(shù)模型、時間序列、極值理論、面板數(shù)據(jù)、久期模型、非線性模型和多元分析等諸多方面［5］．在其他學(xué)科領(lǐng)域包括環(huán)境科學(xué)［6］、生態(tài)學(xué)［7］、生存分析［5］、醫(yī)學(xué)［8］、經(jīng)濟(jì)學(xué)［9］、金融學(xué)［10］等已得到廣泛應(yīng)用．

本文主要考慮非參數(shù)分位數(shù)回歸模型Y=mτ(X)+ ετ，其中 ετ是隨機(jī)誤差，其分布滿足 τ分位數(shù)為0．此模型關(guān)注于在給定X=x條件下，Y的條件 τ分位數(shù)函數(shù) mτ(x)=argmin E{ρτ(Y-mτ(X))|X=x}(0≤τ≤1)．在非參數(shù)分位數(shù)回歸模型下，常用的回歸估計方法包括局部常數(shù)核估計和局部線性核估計．YU和JONES［11］比較了這2種估計的內(nèi)部點漸進(jìn)均方誤(MSE)和邊界點漸進(jìn)均方誤(MSE)，指出在內(nèi)點情形下，兩者的實際差距并不顯著．然而，對于非參數(shù)分位數(shù)回歸，檢驗函數(shù)并不光滑，使得非參數(shù)分位數(shù)回歸的計算成為一個重要的挑戰(zhàn)．使用局部線性的方法，YU 和 JONES［12］利用 MM-算法對分位數(shù)回歸逐點進(jìn)行估計．但是，這種方法有2個缺點:(1)計算量繁重;(2)即使分位數(shù)回歸函數(shù)是光滑的，基于局部線性方法的估計仍可能不光滑甚至不連續(xù)．

本文利用MM-算法對局部常數(shù)分位數(shù)回歸進(jìn)行估計，所得到的目標(biāo)函數(shù)具有MM-算法的下降性質(zhì)，保證了所得估計的數(shù)值穩(wěn)定性．同時，所提出的算法不僅成功地減輕了計算負(fù)擔(dān)，而且，在適當(dāng)?shù)臈l件下，得到的估計保持了回歸函數(shù)連續(xù)光滑的性質(zhì)．

1 局部常數(shù)核加權(quán)分位數(shù)回歸

1．1 局部常數(shù)估計

假設(shè)一組隨機(jī)樣本{(Xi，Yi)(i=1，2，…，n)}來自如下的非參數(shù)分位數(shù)回歸模型

Y=mτ(X)+ ετ(0≤τ≤1)，

其中ετ是隨機(jī)誤差項，滿足第 τ分位數(shù)為0，協(xié)變量X和ετ是獨立的．在分位數(shù)回歸中，我們的任務(wù)是在給定X=x的條件下，估計響應(yīng)變量Y的τ條件分位數(shù)函數(shù) mτ(x)．

在Y的τ分位數(shù)處取得最小值．類似地，在非參數(shù)分位數(shù)回歸中，通過最小化目標(biāo)函數(shù)

來估計mτ(x)．為了使估計具有光滑性，我們采用格點x附近的局部常數(shù)估計來近似條件分位數(shù)函數(shù)，即mτ(X)≈β0，然后最小化核加權(quán)的目標(biāo)函數(shù)

此處Kh(Xi-x)=h-1K((Xi-x)/h)是一個核函數(shù)．

對一個給定的格點x，式(1)是非光滑的，其最小優(yōu)化可通過一個MM迭代算法實現(xiàn)，即構(gòu)造凸函數(shù)包，將不光滑函數(shù)的優(yōu)化問題轉(zhuǎn)化為光滑凸函數(shù)的優(yōu)化．現(xiàn)給定β(l)，并定義ri(β)=Yi-β，這些函數(shù)有以下幾個性質(zhì)．

性質(zhì)1 局部加權(quán)的二次函數(shù)

在β(l)處是l(β)=Σni=1ρτ(Yi- β)Kh(Xi-x)的一個優(yōu)超函數(shù)．即在給定 β(l)下，對所有的 β，有 Qτ，并且 Q

證明 令

因為

MM-算法的第l+1次迭代值β(l+1)是二次函數(shù)的最小值點，通過微分，得到

性質(zhì)2 由式(2)定義的序列{β(l)}滿足l(β(l+1))≤l(β(l))，等式成立當(dāng)且僅當(dāng)

證明 事實上，由性質(zhì)1可知若 Qτ()，即此時二次函數(shù)))在β(l+1)和β(l)處均取得最小值，也就是說 β(l+1)=β(l)，于是

反過來，由式(3)，結(jié)論顯然成立．

性質(zhì)3 假設(shè)MM-算法在有限次迭代后收斂，并且在每次迭代中，對所有i{1，2，．．．，n}都有ri(β(l))≠0，則得到函數(shù)mτ(x)的估計連續(xù)且光滑．

由于ri(β(l))≠0，Kh(Xi-x)=h-1K((Xi-x)/h)是x的連續(xù)且光滑函數(shù)，所以對任一次迭代β(l+1)(x)也是x的連續(xù)且光滑函數(shù)．

1．2 具有擾動的局部常數(shù)估計

在性質(zhì)3中，需要在每次迭代中都有ri(β(l))≠0的條件，才能保證回歸估計連續(xù)且光滑．然而，這個條件在實際數(shù)據(jù)處理時很難保證．

受 HUNTER 和 LANGE［4］的啟發(fā)，在優(yōu)超函數(shù))的二次項系數(shù)的分母中增加一個很小的擾動項ε．于是，定義改進(jìn)的分位數(shù)目標(biāo)函數(shù)為

現(xiàn)在，構(gòu)造Lε(β)在 β(l)處的優(yōu)超函數(shù)

性質(zhì)4 β(0)是R R上任意一點，定義緊集

證明 由ri(β)=Yi-β的連續(xù)性和Ω集合的緊致性可知，其上界d是有限的．如果ε+|r|≤1，則有．如果 ε+＞1，由于 ε＜d-1，那么由d的定義可得 ln(ε+＜ln d＜-lnε．因此，

2 窗寬的選取和數(shù)值模擬

2．1 窗寬選取

對于非參數(shù)核光滑方法來說，窗寬h的選取至關(guān)重要．交叉驗證(Cross-Validation)是選擇光滑參數(shù)最流行的一種方法，它將整個數(shù)據(jù)集做多次劃分，每次劃分將數(shù)據(jù)分為2個部分:訓(xùn)練集與檢驗集，訓(xùn)練集數(shù)據(jù)用來做估計，檢驗集數(shù)據(jù)用來做驗證．本節(jié)的數(shù)據(jù)模擬采用缺一交叉驗證(Leave-One-Out Cross-Validation)的方法來選擇窗寬，也就是每次驗證時檢驗數(shù)據(jù)集中只有一個樣本點，并且把所有樣本點都充當(dāng)一次檢驗數(shù)據(jù)．缺一交叉驗證標(biāo)準(zhǔn)為

2．2 數(shù)值模擬

模擬的數(shù)據(jù)來自模型:Y=sin(2πX)+ε，其中X～U［0，1］，誤差項 ε 與 X 獨立．誤差項 ε 的分布取標(biāo)準(zhǔn)正態(tài) N(0，1)和 t1分布．對 τ=0．25、τ=0．5、τ=0．75這3種不同的分位數(shù)水平，基于我們提出的MM迭代算法，計算了局部常數(shù)估計、最小二乘局部常數(shù)估計，并且基于估計積分平方誤差(ISE)的中位數(shù)進(jìn)行了比較．在樣本數(shù)n=100、n=200、n=500的情形下分別模擬100次，每次采用交叉驗證法選出最優(yōu)窗寬．所得到的結(jié)果見表1，其中LS、QR0．5、QR0．25、QR0．75分別表示最小二乘核估計、MM-算法下的中位數(shù)、1/4、3/4分位數(shù)局部常數(shù)估計．誤差項ε在標(biāo)準(zhǔn)正態(tài)N(0，1)和t1分布下的均值、中位數(shù)均為0，這2種模型下LS和QR0．5均為真實模型的估計，而 QR0．25和 QR0．75是真實模型的一個平移估計，預(yù)計 QR0．25和 QR0．75的 ISE 要高于 QR0．5．

表1 各估計積分平方誤差的中位數(shù)Table 1 Median ISE of these estimators

結(jié)果顯示，在正態(tài)誤差情形下，QR0．5的ISE略高于LS的ISE，兩者差別不大．但是，在重尾的t1誤差分布下，前者的ISE遠(yuǎn)低于后者，這說明我們提出的方法得到的估計在重尾分布的誤差下是穩(wěn)健的．另外，在t1誤差分布下，n=500時，我們畫出了中位數(shù)局部常數(shù)分位數(shù)回歸(MLCQR)的估計擬合曲線、最小二乘局部常數(shù)(LSLC)的估計曲線以及真實模型(RM)．從圖1看到在重尾誤差的情況下，提出的方法能刻畫真實模型．

圖1 n=500時中位數(shù)回歸估計曲線及LS估計曲線Figure 1 Median regression estimated curve and LS estimated curve when n=500

2．3 實證研究

應(yīng)用所提的方法去分析一份青少年骨密度的數(shù)據(jù)．該數(shù)據(jù)最早出現(xiàn)于文獻(xiàn)［13］．隨后，在文獻(xiàn)［14］中也被作為一個范例進(jìn)行分析．?dāng)?shù)據(jù)下載網(wǎng)址為:http://www-stat．stanford．edu/～ tibs/ElemStatLearn/，包含了4個變量:idnum、age(測量時孩子的平均年齡)、gender(男或女)和spnbmd(相對的脊椎骨密度測量值)，共有458個觀測值．本文所關(guān)注的是相對的脊椎骨密度測量值隨年齡的變化情況．

由于事先并不知道相對的脊椎骨密度測量值隨年齡的變化情況，所以采用非參數(shù)回歸模型去擬合其變化情況．運用本文所提出的基于MM算法的局部常數(shù)核加權(quán)分位數(shù)回歸去估計非參數(shù)回歸函數(shù)．這里年齡為自變量X，相對的脊椎骨密度測量值為因變量Y．

在MM-算法中，首先采用核回歸獲得均值函數(shù)的估計，即

在常數(shù)方差的假設(shè)下，進(jìn)一步計算方差的估計:

由正態(tài)近似，取MM算法中的迭代初始值為

其中Φ-1(τ)是標(biāo)準(zhǔn)正態(tài)經(jīng)驗分布函數(shù)的反函數(shù)．MM算法中的迭代終止條件是前后2次迭代估計之差小于10-6．在 τ=0．1、τ=0．25、τ=0．5、τ=0．75 和 τ=0．9的分位數(shù)水平下，通過缺一交叉驗證法則選出的最優(yōu)窗寬分別是 h0．1=0．59、h0．25=0．59、h0．5=0．60、h0．75=0．52 和 h0．9=0．55．其估計的結(jié)果見圖 2．

圖2 局部常數(shù)核加權(quán)分位數(shù)回歸的估計曲線Figure 2 Local constant kernel-weighted quantile vegression estimated curves

圖2表明，所提出的方法得到了光滑的回歸曲線．在12～13歲之間相對的脊椎骨密度達(dá)到最大．在此之前，相對的脊椎骨密度隨著年齡的增大而增加;在此之后，隨著年齡的增加反而降低．同時，隨著分位點τ的增加，相對應(yīng)的回歸曲線起伏越大．這說明，對不同水平的相對的脊椎骨密度，所反映的趨勢各不相同，也就是說，不同的相對的脊椎骨密度隨年齡的變化不一樣．特別地，相對的脊椎骨密度的分位數(shù)為10%左右的人群，其相對的脊椎骨密度隨年齡的變化不大;而分位數(shù)為90%左右的人群，其相對的脊椎骨密度隨年齡的變化幅度比較大．

3 總結(jié)

本文針對非參數(shù)分位數(shù)回歸模型，提出了一種新的基于MM-算法的局部常數(shù)核加權(quán)計算法．在某些條件下，由此計算出的非參數(shù)估計是連續(xù)光滑的．同時，利用缺一交叉驗證法則來選擇協(xié)調(diào)參數(shù)．雖然本文只考慮了非參數(shù)回歸模型，但是所提出的方法同樣適用于半?yún)?shù)回歸模型．

［1］KOENKER R，BASSETT G W．Regression quantiles［J］．Econometrica，1987，46(1):33-50．

［2］KOENKER R，OREY D．Algorithm AS 229:Computing regression quantiles［J］．Appl Stat，1987，36(3):383-393．

［3］KOENKER R，PARK B J．An interior point algorithm for nonlinear quantiles regression［J］． J Econometrics，1996，71(1):265-283．

［4］HUNTERD R，LANGEK．Quantile regression via an MM algorithm［J］．JComput Graph Stat，2000，9(1):60-77．

［5］KOENKER R，HALLOCK K．Quantile regression:An introduction［J］．JEcon Perspect，2001，15:143-156．

［6］CHOCK D P，WINKLER S L，CHEN C．A study of the association between daily mortality and ambient air pollutant concentrations in Pittsburgh，Pennsylvania［J］．Journal of the Air＆ Waste Management Association，2000，50(8):1481-1500．

［7］KOENKER R，SCHORFHEIDE F．Quantile splinemodels for global temperature change［J］．Climate Change，1994，28(4):395-404．

［8］COLE T J，GREEN P J．Smoothing reference centile curves:The LMSmethod and penalized likelihood［J］．Stat Med，1992，11(10):1305-1319．

［9］GOSLING A，MACHIN S，MEGHIR C．What has happened to men's wages since the mid-1960s［J］Fiscal Studies，1994，15(4):63-87．

［10］BASSETT G，CHEN H．Quantile style:Return-based attribution using regression quantiles［J］．Empir Econ，2001，26(1):7-40．

［11］YU K，JONESM C．A comparison of local constant and local linear regression quantile estimators［J］．Comput Stat Data An，1997，25(2):159-166．

［12］YU K，JONESM C．Local linear quantile regression［J］．JAmer Stat Assoc，1998，93(441):228-237．

［13］BACHRACH L K，HASTIE T，WANFM C．Bonemineral acquisition in healthy Asian，Hispanic，Black，and Caucasian youth:A longitudinal study［J］．JClin Endocr Metab，84(12):4702-4712．

［14］HASTIE T，TIBSHIRANL R，F(xiàn)RIEDMAN J．The elements of statistical leavning［M］．New York:Springer，2001．