聶浩巍，李志強(qiáng)

（北京化工大學(xué) 數(shù)理學(xué)院，北京 100029）

0 引言

分位數(shù)回歸（Quantile Regession，QR）由Koenker 和Basset（t1978）[1]提出，相比于傳統(tǒng)的均值回歸，分位數(shù)回歸可以研究不同分位數(shù)下協(xié)變量對結(jié)果的影響，而不需要對誤差作出任何假設(shè)，因此更加靈活和穩(wěn)健。從分位數(shù)回歸被正式提出至今，學(xué)者們不斷地研究其各種參數(shù)估計，并成功將其應(yīng)用于計量經(jīng)濟(jì)學(xué)、醫(yī)學(xué)等不同的領(lǐng)域中。

如今，海量數(shù)據(jù)集常見于各大研究領(lǐng)域，有時數(shù)據(jù)集甚至以流的形式出現(xiàn)。然而，傳統(tǒng)分位數(shù)回歸需要同時處理整個數(shù)據(jù)集，而海量數(shù)據(jù)集由于內(nèi)存限制很難由單獨的一臺計算機(jī)進(jìn)行處理。為了解決這個問題，已經(jīng)開發(fā)出了許多基于分治（Divide-and-Conquer，DC）的估計算法。他們大致分為兩大類：一類是基于多輪通信的迭代算法，它通過子機(jī)器與主機(jī)器間多輪信息傳輸進(jìn)行迭代以達(dá)到處理全數(shù)據(jù)集的目的[1]，從而快速得到有效的估計量[2—4]。該方法的估計效率較高，但除Chen等（2019）[2]基于核光滑估計方程的迭代算法外，他們均無法處理流數(shù)據(jù)。另一類是只需要一輪通信（One-shot）的分治算法，首先從各個子數(shù)據(jù)集中得到局部估計量，然后通過簡單平均或加權(quán)平均進(jìn)行聚合，最終得到聚合估計量[5,6]。其中，Lin 和X（i2011）[5]通過展開估計方程得到了一種行之有效的聚合估計算法（Aggregated Estimating Equation Estimation,AEEE），但AEEE要求估計方程可微。然而眾所周知，分位數(shù)回歸的估計方程是不可微的，因此AEEE不能直接用于分位數(shù)回歸。Chen和Zhou（2020）[6]改進(jìn)AEEE并成功將其應(yīng)用于分位數(shù)回歸中。然而，他們的方法需要通過使用重采樣方法獲得權(quán)重矩陣，這無疑降低了計算速度。因此，有必要開發(fā)一種計算效率高且適用于流數(shù)據(jù)的算法，用于海量數(shù)據(jù)分位數(shù)回歸的參數(shù)估計。

對此，本文建議使用Fernandes等（2021）[7]提出的光滑方法，將分位數(shù)回歸的求解問題光滑化，從而滿足AEEE中的可微條件，由此提出一種計算高效的海量數(shù)據(jù)下光滑分位數(shù)回歸聚合估計（Divide-and-Conquer Smoothing Quantile Regession，DCSQR）算法。具體而言，本文先計算每個數(shù)據(jù)塊的核卷積光滑估計方程估計量和對應(yīng)的Hessian矩陣，并只需要保留每個數(shù)據(jù)塊的這兩個統(tǒng)計信息。若數(shù)據(jù)是以流的形式接收，則可以不斷計算和保存相應(yīng)統(tǒng)計信息并丟棄原數(shù)據(jù)集。最后，通過AEEE得到原數(shù)據(jù)集的有效估計量。本文將通過詳細(xì)的理論證明給出該估計量的漸近正態(tài)性，并通過模擬研究和實證分析證實該方法的有效性。

1.1 光滑分位數(shù)回歸模型

給定Y∈? 為單變量響應(yīng)變量，X=(x1，…，xp)T∈?p為p維協(xié)變量向量,其中x1≡1。假設(shè)數(shù)據(jù)集D={Yi，中含有來自(Y，X)的N個i.i.d.的樣本，在給定分位數(shù)水平τ∈(0，1)下，本文考慮線性分位數(shù)回歸模型為：

其中，β0(τ)為關(guān)于τ的p維回歸參數(shù)真值向量，εi滿足P[εi≤0|Xi]=τ。為簡單起見，下文將省略τ。

分位數(shù)回歸估計[1]可通過求解如下最小化問題得到：

其中，ρτ(u)=u(τ-I(u＜0)) 是檢查損失函數(shù)（check loss function），而I(·)是示性函數(shù)。根據(jù)Buchinsky（1998）[8]的研究，可通過求解以下估計方程來獲得β0的經(jīng)典估計方程估計量

其中，ψτ(u)=τ-I(u＜0)為檢查函數(shù)。

然而,由于估計方程（3）不可微，因此Lin 和X（i2011）[5]的方法不能推廣到分位數(shù)回歸中。為了避免估計方程的不可微性，本文使用Fernandes等（2021）[7]所提出的核卷積光滑（Kernel Convolution Smoothing）方法，最小化以下光滑分位數(shù)回歸（Smoothing Quantile Regession，SQR）的目標(biāo)函數(shù)來求解模型（1）中β0的分位數(shù)回歸估計量

估計問題式（3）轉(zhuǎn)化為求解以下光滑估計方程：

1.2 海量數(shù)據(jù)下光滑分位數(shù)回歸聚合估計

當(dāng)樣本量N過大時，由于單臺計算機(jī)內(nèi)存有限，直接解決式（6）中的估計方程并不可行，因此本文考慮使用AEEE方法解決上述問題。將數(shù)據(jù)集D隨機(jī)分為K塊，每塊含有n個數(shù)據(jù)，各塊數(shù)據(jù)集分別記為…，K，其中N=nK，以保證每個塊都可以存儲在計算機(jī)的內(nèi)存中。對于每塊數(shù)據(jù)集Dk，其對應(yīng)的光滑估計方程為：

當(dāng)Rk足夠小時，通過簡單的推導(dǎo)就可以得到一個式（6）的閉式近似解

1.3 理論性質(zhì)

A1：參數(shù)空間?是Rp的緊子集，參數(shù)向量β0是?的內(nèi)點。

A2：Xki有有界支撐，且Σ0=E[Xi XiT]非奇異。

A3：對于所有0的鄰域內(nèi)的u和幾乎所有的x，f(u|x)存在并遠(yuǎn)離0和∞，且r階對u連續(xù)可微。

A5：窗寬? 滿足當(dāng)n→∞時：（a）N?2r→0 ；（b）N?/lgN→∞。

A6：D0=E[Xi XiTf(0|Xi)]正定且有界。

下面的定理給出了估計量的漸近性質(zhì)。為了證明定理1，先給出引理1。

利用分部積分公式可得：

對F(-?v|Xi)在0處進(jìn)行r階泰勒展開，可得：

（b）本文僅證明第一個式子,第二個式子同理。利用分部積分公式，可得：

（c）對F(-?v|Xi)在0處進(jìn)行r階泰勒展開，可得：

引理1證畢。

定理1：假設(shè)條件A1至A4和A5（a）成立，則有：

根據(jù)Lindeberg中心極限定理，?ε＞0，都有：

定理1給出了核卷積光滑估計量的漸近性質(zhì)。對于各塊數(shù)據(jù)下的局部核卷積光滑估計量，該定理也同樣適用。

證明：由條件A2 和A6 可以得到Ak是正定的。根據(jù)式（10）可得：

定理2表明，當(dāng)K以慢于子數(shù)據(jù)集大小n的速度趨于無窮大時，是β0的相合估計量。

為了證明定理3，先證明引理2。

在β0的η鄰域內(nèi)，對GN(β)使用中值定理：

引理2證畢。

1.4 估計算法

由于SQR具有優(yōu)良的性質(zhì)，因此本文可以使用高效的Newton-Raphson 迭代算法來估計，并避免了對討厭參數(shù)的額外估計，從而降低了計算成本。為了進(jìn)一步降低計算成本，本文選擇使用第一塊子數(shù)據(jù)集的標(biāo)準(zhǔn)QR 估計量作為每一塊數(shù)據(jù)的迭代初值。具體算法如下：

步驟1：參數(shù)設(shè)置：給定窗寬h與核函數(shù)K(·)。

步驟2：將數(shù)據(jù)集D分割成K塊,并將各小塊數(shù)據(jù)集Dk分別發(fā)送給各子節(jié)點。

2 模擬研究

本文使用蒙特卡羅模擬來檢驗所提出算法在線性模型下的有限樣本性能。所有程序都是用Python編寫的，并在搭建好的Spark集群上運行，集群包含3 臺內(nèi)存為8G 的計算機(jī)，其中一臺為主節(jié)點，另外兩臺為子節(jié)點。在模擬實證中，統(tǒng)一設(shè)置核函數(shù)為標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)最優(yōu)窗寬的選擇可參考Fernandes 等（2021）[7]或He 等（2023）[9]研究中的最小化漸近均方誤差（Asymptotic Mean Square Error，AMSE）。由于在模擬研究中結(jié)果對窗寬不敏感，因此為簡單起見，將窗寬固定為?=1.5N-1/3。

Case 1：同方差正態(tài)分布,?i～N(0，1)。

Case 2：異方差正態(tài)分布,?i～N(0，(1+0.1Xi1)2)。

Case 3：指數(shù)分布，?i～Exp(1)。

因此,對于任何給定的分位數(shù)水平τ,給定X的Y的τ條件分位數(shù)分別為：

Case 1：同方差正態(tài)分布，θ(τ)=θ0+Φ-1(τ)(1，0，0，0，0)T。

Case 2：異方差正態(tài)分布，θ(τ)=θ0+Φ-1(τ)(1，0，0，0，0)T+0.1Φ-1(τ)(0，1，0，0，0)T。

Case 3：指數(shù)分布，θ(τ)=θ0+Fexp-1(τ)。

其中，Φ和Fexp分別為服從標(biāo)準(zhǔn)正態(tài)分布和均值為1的指數(shù)分布的向量。

為了證實本文方法的有效性，將總樣本量固定為N=1000000，令K在{10，50，100，200，500，1000}內(nèi)取值，并分別在以上3種不同隨機(jī)誤差下重復(fù)模擬實驗100次。

本文給出了在分位數(shù)水平τ=0.25，0.5，0.75 下，估計量的平均均方誤差（Mean Squared Error，MSE）MSE=和計算時間，并將結(jié)果與Chen 和Zhou（2020）[6]的算法（Divide-and-Conquer Quantile Regession，DCQR）進(jìn)行對比，用以證明本文算法的性能。模擬結(jié)果基于100次模擬重復(fù)實現(xiàn)。

從表1 中可以看到，DCSQR 比DCQR 花費的時間更少，這是因為DCSQR不需要額外估計權(quán)重矩陣Ak。這證明了DCSQR在計算速度上的優(yōu)越性。

表1 不同環(huán)境下DCSQR與DCQR計算時間對比（單位：秒）

而從圖1 中可以看出，當(dāng)K≤200 時，DCSQR 的MSE曲線與DCQR相近且變化幅度更??；當(dāng)K＞200 時，在多數(shù)情況下DCQR 的MSE 小于DCSQR。這說明當(dāng)K≤200 時DCSQR 的MSE 與DCQR 的穩(wěn)健性相當(dāng)。注意到，即使本文模擬研究中使用的分塊K的數(shù)量超過了定理4 中的理論限制(N=1000000 時分塊理論上限K≈31)，DCSQR 在K=200 時也仍然表現(xiàn)良好，當(dāng)K＞200 以后MSE才快速增大，這意味著關(guān)于K的理論條件可以進(jìn)一步放寬。

3 實證分析

本文將所提出的算法應(yīng)用于UCI 機(jī)器學(xué)習(xí)存儲庫報告的溫室氣體（GHG）觀測網(wǎng)絡(luò)數(shù)據(jù)集。該數(shù)據(jù)集由955167 個觀測值組成。響應(yīng)變量是合成觀測的GHG 濃度，共有15 個預(yù)測因子。這些預(yù)測因子是加利福尼亞州14個不同空間區(qū)域和加利福尼亞州以外一個區(qū)域排放的示蹤物的GHG濃度（記為Reg1-Reg15）。

本文先評估了DCSQR 方法的預(yù)測精度，并將其與DCQR 方法進(jìn)行比較。數(shù)據(jù)集被分為訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集，其中訓(xùn)練數(shù)據(jù)集含有900000個數(shù)據(jù)，而測試數(shù)據(jù)集含有55167個數(shù)據(jù)。然后，本文將訓(xùn)練數(shù)據(jù)隨機(jī)分成K塊（K∈{5，10，20，50，100，200}），并分別通過DCSQR 和DCQR兩種算法來估計回歸系數(shù)。

表2 給出了兩種算法在預(yù)測精度和計算成本（總秒數(shù)）兩個方面的比較結(jié)果。從表2中可以看到，DCQR方法較為穩(wěn)定，當(dāng)K≤50 時，DCSQR 幾乎與DCQR 相當(dāng)；當(dāng)K≥100 時DCQR 比DCSQR 稍好，這證明了當(dāng)K≤50 時DCSQR 的預(yù)測穩(wěn)健性。此外，該表還比較了這兩種方法的時間成本。顯然，從表2 中可以看到，不論分多少塊，DCSQR的用時都要遠(yuǎn)遠(yuǎn)小于DCQR的用時。

表2 GHG數(shù)據(jù)集下兩種算法的MAPE和計算總秒數(shù)對比

4 結(jié)論

本文提出了一種基于光滑估計方程的聚合估計算法DCSQR，用于解決海量數(shù)據(jù)下分位數(shù)回歸的參數(shù)估計問題。理論研究證明，當(dāng)K以慢于n的一定速度趨于無窮大時，聚合估計量具有和經(jīng)典分位數(shù)回歸相同的漸近正態(tài)性。模擬實證表明，由于避免了對討厭參數(shù)的估計，DCSQR 算法在保持原有估計精度的基礎(chǔ)上，和Chen 和Zhou（2020）[6]基于經(jīng)典分位數(shù)回歸的DCQR 算法相比，計算效率顯著提高，這證明了DCSQR的有效性。