張輝國,張孟娟,胡錫健

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆 烏魯木齊 830046)

1 引言

變系數(shù)模型作為模擬變量間回歸關(guān)系動態(tài)變化的重要工具,是經(jīng)典線性回歸模型的有效擴展,在社會科學(xué)和自然科學(xué)領(lǐng)域受到了廣泛的關(guān)注[1,2]。該模型中的回歸系數(shù)允許隨其它協(xié)變量的值光滑變化,而不是設(shè)為固定常數(shù),有效解決了非參數(shù)回歸研究中的“維數(shù)災(zāi)難”問題,并且繼承了經(jīng)典線性回歸模型的簡單性和易解釋性。其中系數(shù)函數(shù)的估計是變系數(shù)模型的關(guān)鍵問題,人們一直致力于為其發(fā)展有效的估計方法[3-5]。

作為目前最流行且有效的方法之一,支持向量機(Support Vector Machines, SVM)將線性思想應(yīng)用于非線性數(shù)據(jù),主要解決分類和非線性函數(shù)估計問題,由Vapnik等人[6]于1995年提出后被許多人進一步研究[7-10]。Suykens等人[11]在1999年提出了最小二乘支持向量機(Least Squares Support Vector Machine, LSSVM),該方法是一種基于SVM的改進算法,在繼承SVM優(yōu)點的同時,采用等式約束替代SVM中的不等式約束,并使用誤差的2-范數(shù)替代SVM中的ε-不敏感損失函數(shù),從而將求解SVM的凸二次規(guī)劃問題轉(zhuǎn)化為求解線性方程組的問題,降低了算法復(fù)雜度,求解速度較快,在各個領(lǐng)域中都得到一定應(yīng)用[12-13]和進一步的研究發(fā)展[14-16]。但LSSVM存在兩個潛在的缺陷:一是解丟失了稀疏性[17-18];二是損失函數(shù)中采用誤差平方和度量損失,若訓(xùn)練數(shù)據(jù)中存在異常值以及誤差不服從高斯分布時,LSSVM的穩(wěn)健性較差。為此,國內(nèi)外許多學(xué)者針對穩(wěn)健LSSVM進行了深入研究,增加LSSVM穩(wěn)健性的方法主要分為四個方面[19]:基于異常值剔除技術(shù)改進、基于加權(quán)函數(shù)改進、基于p-范數(shù)改進以及基于損失函數(shù)改進。

用于回歸的LSSVM問題一般稱為最小二乘支持向量回歸機(LS-SVR),基于LS-SVR對于非線性數(shù)據(jù)的適應(yīng)性以及計算的簡單性,Shim和Hwang[20]在2015年提出了利用LS-SVR技術(shù)擬合變系數(shù)模型的方法VC-LS-SVR,給出了廣義交叉驗證法來選擇算法中的超參數(shù),并提供了構(gòu)造估計系數(shù)函數(shù)置信區(qū)間的方法。該方法被證明估計系數(shù)函數(shù)的表現(xiàn)優(yōu)于常用的局部多項式擬合方法,是一種簡單且高效的新方法。由于該方法把LS-SVR直接應(yīng)用于變系數(shù)模型,所以也具有其面對異常值不穩(wěn)健的缺點,當(dāng)存在異常值時,可能會導(dǎo)致系數(shù)函數(shù)的估計失效,因此,本文將在原始的VC-LS-SVR方法框架下,基于加權(quán)函數(shù)提出兩種變系數(shù)模型的穩(wěn)健最小二乘支持向量回歸估計方法,預(yù)期在數(shù)據(jù)包含異常值時,兩種穩(wěn)健估計方法在估計系數(shù)函數(shù)方面比VC-LS-SVR方法有更好的表現(xiàn)。最后還做了數(shù)值實驗,以評估所提出方法在恢復(fù)真實回歸系數(shù)曲線方面的穩(wěn)健性,并對三種方法估計系數(shù)函數(shù)的性能進行了全面的比較。

2 VC-LS-SVR方法的介紹

(1)

其中xij是xi的第j個分量,βj(·),j=1,2,…,p是需要估計的未知系數(shù)函數(shù),Var(yi)=σ2(ti,xi)>0,εi是均值為0、方差為1的獨立同分布的隨機變量。

其中σ>0和d是需預(yù)先指定的核參數(shù)。

則式(1)中回歸函數(shù)f(·)可重寫為

基于LS-SVR的思想,VC-LS-SVR的優(yōu)化問題可定義為

其中γ為正則化參數(shù),ei是獨立同分布的隨機變量(均值為0,方差Var(e)<∞)。

通過構(gòu)造拉格朗日函數(shù)和KKT最優(yōu)條件求解該優(yōu)化問題,結(jié)果由下列線性方程組給出:

(2)

其中X=(x1,x2,…,xn)t,Y=(y1,y2,…,yn)t,α=(α1,α2,…,αn)t,b=(b1,b2,…,bj)t,In表示一個維度為n的單位矩陣,Op×p表示一個p維的零矩陣,0p×1表示一個p維的零向量,K是由元素Kij=K(ti,tj)構(gòu)成的n×n核矩陣,⊙表示為點乘,即每個矩陣元素對應(yīng)相乘。

給定一個數(shù)據(jù)點(t0,x0),系數(shù)函數(shù)的VC-LS-SVR估計形式為

j(t0)=∑ni=1xijK(t0,ti)i+j,j=1,2,…,

(3)

對應(yīng)的回歸函數(shù)估計值為

(t0,x0)=∑pj=1x0j(∑ni=1xijK(t0,ti)i+j)

(4)

3 變系數(shù)模型的穩(wěn)健LS-SVR估計

上節(jié)中提到的VC-LS-SVR方法將LS-SVR直接應(yīng)用于變系數(shù)模型,因此也具有LS-SVR面對異常值不穩(wěn)健的缺點,為了在之前的VC-LS-SVR解的基礎(chǔ)上獲得穩(wěn)健估計,在后續(xù)的討論中,將利用WLS-SVR[22]的思想,通過加權(quán)因子vi對每個樣本點的誤差變量ei=αi/γ進行加權(quán),則有如下優(yōu)化問題

其中γ為正則化參數(shù),且ei為均值為0,同方差Var(e)<∞的獨立同分布的隨機變量。

構(gòu)造拉格朗日函數(shù)

根據(jù)KKT(Karush-Kuhn-Tucker)最優(yōu)條件可得

將w、ei的顯式解代入第4個偏導(dǎo)式,寫成矩陣形式可得

(5)

(6)

=1.483MAD(i)

(7)

其中MAD代表絕對中位差,常數(shù)c1、c2通常被選擇為c1=2.5和c2=3。由此可給出以下算法:

算法1:VC-WLS-SVR

Step3:根據(jù)式(6)確定權(quán)重vi;

在文獻[23]中,Brabanter等人證實了在帶有重尾的非高斯噪聲分布下,加權(quán)方法會失效。因此為了抵御極度異常情況的存在,本文又提出了一種迭代重加權(quán)的方案,這種方法簡稱為VC-IRLS-SVR,其基本思想是對那些極度異常的樣本點多次降權(quán),使其賦權(quán)結(jié)果趨于穩(wěn)定。由此引出如下算法:

算法2:VC-IRLS-SVR

設(shè)定迭代次數(shù)k=0;

Step5:設(shè)定k=k+1;

4 數(shù)值分析

在本節(jié)中,進行了一些數(shù)值實驗,并構(gòu)造合適的評價指標以評估VC-LS-SVR、VC-WLS-SVR以及VC-IRLS-SVR方法在數(shù)據(jù)有無異常值的情況下得到的系數(shù)估計值的準確性。

實驗數(shù)據(jù)由以下變系數(shù)模型生成:

yi=β0(ti)+β1(ti)xi+εi

(8)

其中ti～U(0,1),xi～U(0,2),系數(shù)函數(shù)設(shè)定為復(fù)雜的非線性周期函數(shù)

β0(t)=cos(7πt);β1(t)=sin(7πt)

式(8)中模型的誤差項εi(i=1,2,…,n)獨立的從下列分布中抽取:

case1:正態(tài)分布N(0,0.52);

case2:位置參數(shù)為0,尺度參數(shù)為1的柯西分布C(0,1);

case3:污染正態(tài)分布(1-δ)N(0,0.52)+δN(0,82),其中δ為污染率,在模擬中依次取值為10%、20%、30%、40%、50%;

case4:另一種污染正態(tài)分布(1-δ)N(0,0.52)+δLaplace(0,2,2),也稱Huber分布,同樣地δ為污染率,依次取值為10%、20%、30%、40%、50% 。

當(dāng)模型誤差取自case1中的正態(tài)分布N(0,0.52)時,可認為是數(shù)據(jù)中不存在異常值的情況。當(dāng)模型誤差取自case2中的柯西分布C(0,1)時,由于柯西分布是厚尾分布,其均值和方差都不存在,因此取自這個分布的誤差一般都會包括一些極端異常的值。而case3中的污染正態(tài)分布,其主體是正態(tài)分布,尾部是均值相同、方差更大的正態(tài)分布,由于尾部方差與主體方差差別很大,可能會產(chǎn)生一些過大或過小的數(shù)值。同樣地,case4中的污染正態(tài)分布其主體是正態(tài)分布,尾部是比正態(tài)分布厚尾的拉普拉斯分布,更易于產(chǎn)生較大或較小的數(shù)值,由于想在真實曲線上下產(chǎn)生數(shù)目不對稱的異常點,來評估穩(wěn)健方法的有效性,因此這里采用非對稱拉普拉斯分布。

為了評價所提出方法估計系數(shù)函數(shù)的準確性,在下文中定義了一些評價指標來衡量每種方法在估計系數(shù)時的表現(xiàn)。

首先在數(shù)值實驗過程中,為減少抽樣誤差,對生成的x重復(fù)產(chǎn)生N次誤差項,本文設(shè)置內(nèi)循環(huán)次數(shù)N=100,即可得到給定點ti處第j個回歸系數(shù)的平均估計值

m(j(ti))=1N∑Nk=1(k)j(ti)

mse(j)=1n∑ni=1(βj(ti)-m(j(ti)))2

雖然在實際應(yīng)用中βj(ti)的真實值是未知的,因此mse(j)是無法計算出來的,但通過數(shù)值實驗,可以知道真實的βj(ti),從而證明所提出方法的有效性。

一次的計算可能會具有偶然性,因此本文設(shè)置外循環(huán)次數(shù)M=50,將x重復(fù)生成M次,則可得到M個mse(j),最終可計算出mse(j)的均值及標準差

Mmse(j)=1M∑Ml=1msel(j)

SDmse(j)=1M∑Ml=1(msel(j)-Mmse(j))2

在本次數(shù)值實驗中,使用高斯核描述系數(shù)函數(shù)βj(t)與光滑變量t之間的非線性關(guān)系,并使用廣義交叉驗證方法[20]確定正則化參數(shù)γ以及核參數(shù)σ。由于權(quán)重的范圍為0≤vi≤1,而本文提出的迭代重加權(quán)方法是權(quán)重之間相乘以達到對異常樣本多次降權(quán)的目的,當(dāng)數(shù)據(jù)存在極端異常值時,VC-IRLS-SVR方法如果在進行迭代時以滿足停止條件停止迭代的話,可能會發(fā)生所有樣本被降權(quán)接近為0的情況,因此,除了設(shè)置停止條件外,還給定了最大迭代次數(shù),在數(shù)據(jù)正常或異常樣本較少時設(shè)定最大迭代次數(shù)k=15,在數(shù)據(jù)存在極端異常值時設(shè)定最大迭代次數(shù)k=3。本文數(shù)值實驗結(jié)果如表1-2及圖1-2所示。

從表1和表2中可以看出,當(dāng)誤差項來自正態(tài)分布時,三種方法估計β0(t)與β1(t)的準確性與穩(wěn)定性相差不大,考慮到此時數(shù)據(jù)無異常,VC-WLS-SVR方法會產(chǎn)生幾乎全等于1的權(quán)重矩陣,實際上也就是VC-LS-SVR方法,且VC-IRLS-SVR方法由于權(quán)重矩陣快速穩(wěn)定此時也不會進行多次迭代,因此三種方法的性能大致相同,這一結(jié)果符合預(yù)期。

表1 0 (t)的Mmse(j)和SDmse(j)

表1 0 (t)的Mmse(j)和SDmse(j)

誤差分布評價指標VC-LS-SVRVC-WLS-SVRVC-IRLS-SVRN(0,0.52)Mmse(^βj)0.0170 0.0182 0.0182 SDmse(^βj)0.0038 0.0039 0.0039 C(0,1)Mmse(^βj)2218.1360 0.4047 0.0482 SDmse(^βj)14838.3500 2.0304 0.0183 (1-δ)N(0,0.52 )+δN(0,82 )δ=10%Mmse(^βj)0.0283 0.0205 0.0232 SDmse(^βj)0.0107 0.0046 0.0051 δ=20%Mmse(^βj)0.0411 0.0235 0.0306 SDmse(^βj)0.0172 0.0065 0.0069 δ=30%Mmse(^βj)0.0457 0.0264 0.0353 SDmse(^βj)0.0179 0.0067 0.0055 δ=40%Mmse(^βj)0.0614 0.0276 0.0106 SDmse(^βj)0.0309 0.0206 0.0122 δ=50%Mmse(^βj)0.0864 0.0637 0.0335 SDmse(^βj)0.0445 0.0444 0.0303 (1-δ)N(0,0.52 )+δLaplace(0,2,2)δ=10%Mmse(^βj)0.1421 0.0254 0.0226 SDmse(^βj)0.0684 0.0082 0.0053 δ=20%Mmse(^βj)0.3183 0.0506 0.0280 SDmse(^βj)0.1054 0.0227 0.0079 δ=30%Mmse(^βj)0.5620 0.1002 0.0398 SDmse(^βj)0.1981 0.0454 0.0141 δ=40%Mmse(^βj)0.9959 0.2667 0.0819 SDmse(^βj)0.2673 0.1221 0.0673 δ=50%Mmse(^βj)1.4259 0.5044 0.2172 SDmse(^βj)0.3656 0.1889 0.1078

表2 1 (t)的Mmse(j)和SDmse(j)

誤差分布評價指標VC-LS-SVRVC-WLS-SVRVC-IRLS-SVRN(0,0.52)Mmse(^βj)0.0117 0.0125 0.0125 SDmse(^βj)0.0031 0.0031 0.0032 C(0,1)Mmse(^βj)836.5237 0.2760 0.0423 SDmse(^βj)5143.3330 1.1785 0.0421 (1-δ)N(0,0.52 )+δN(0,82 )δ=10%Mmse(^βj)0.0214 0.0144 0.0162 SDmse(^βj)0.0069 0.0034 0.0037 δ=20%Mmse(^βj)0.0314 0.0166 0.0219 SDmse(^βj)0.0142 0.0051 0.0051 δ=30%Mmse(^βj)0.0346 0.0189 0.0252 SDmse(^βj)0.0158 0.0051 0.0079 δ=40%Mmse(^βj)0.0439 0.0196 0.0068 SDmse(^βj)0.0202 0.0092 0.0045 δ=50%Mmse(^βj)0.0619 0.0463 0.0238 SDmse(^βj)0.0293 0.0252 0.0153 (1-δ)N(0,0.52 )+δLaplace(0,2,2)δ=10%Mmse(^βj)0.0720 0.0171 0.0155 SDmse(^βj)0.0365 0.0063 0.0046 δ=20%Mmse(^βj)0.0955 0.0285 0.0197 SDmse(^βj)0.0432 0.0116 0.0053 δ=30%Mmse(^βj)0.1347 0.0468 0.0251 SDmse(^βj)0.0672 0.0224 0.0094 δ=40%Mmse(^βj)0.1969 0.0921 0.0396 SDmse(^βj)0.0951 0.0534 0.0335 δ=50%Mmse(^βj)0.2452 0.1504 0.0927 SDmse(^βj)0.1252 0.0781 0.0515

當(dāng)誤差項來自柯西分布(即數(shù)據(jù)中存在極端異常值)時,VC-LS-SVR方法估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值非常大,是VC-WLS-SVR方法的幾千倍,這是由于最小二乘方法本身就不穩(wěn)健,回歸系數(shù)函數(shù)估計結(jié)果被極端異常值扭曲,即使可以調(diào)整學(xué)習(xí)參數(shù)使VC-LS-SVR方法獲得穩(wěn)健解,但在這種含有極端異常值的情況下,它已經(jīng)完全失效。相比之下,VC-WLS-SVR方法表現(xiàn)得就較為良好,VC-IRLS-SVR方法更為占優(yōu)勢,其估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值比前兩種方法要小的多,能產(chǎn)生更為準確和穩(wěn)定的解。

當(dāng)誤差項來自case3中的污染正態(tài)分布時,三種方法在估計系數(shù)的準確性與穩(wěn)定性方面雖然有差異,但都表現(xiàn)出了一定的穩(wěn)健性。首先隨著污染率的增加,VC-LS-SVR方法可通過設(shè)置更高的正則化參數(shù)γ去懲罰誤差平方和損失函數(shù)以增強穩(wěn)健性,因此其估計系數(shù)函數(shù)的性能即使在污染率極高的情況下,也表現(xiàn)的較為良好。但在相同污染率下VC-WLS-SVR方法與VC-IRLS-SVR方法更為穩(wěn)健,在污染率不高(10%～30%)時,兩者估計的性能不相上下,當(dāng)污染率高達40%、50%時,VC-IRLS-SVR方法估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值就比前兩個方法小得多,估計系數(shù)函數(shù)更為準確和穩(wěn)定。

當(dāng)誤差項來自case4中的污染正態(tài)分布時,由于其主體是正態(tài)分布,尾部是比正態(tài)分布厚尾且為非對稱的拉普拉斯分布,因此會在真實曲線上下產(chǎn)生數(shù)目不對稱的異常點。在這種誤差分布下,同一污染率中VC-IRLS-SVR方法估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值總是最小的,VC-WLS-SVR方法估計系數(shù)函數(shù)的準確性和穩(wěn)定性略差一些,而VC-LS-SVR方法與前兩者相比,估計性能相比差異較大。盡管隨著污染率的增加三種方法估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值都在增加,但VC-IRLS-SVR方法在高污染率下的表現(xiàn)也依舊保持了準確與穩(wěn)定。

從圖1也可以看出,當(dāng)數(shù)據(jù)含有極端異常值時,VC-LS-SVR方法估計已失效,估計出的曲線受異常值影響劇烈波動,而所提出的兩種穩(wěn)健方法基本可以還原真實曲線趨勢。圖2展示了當(dāng)數(shù)據(jù)含有非對稱異常值,污染率δ高達50%時,三種方法對截距項β0(t)的估計更易受到非對稱異常值的影響,估計曲線會被整體拉向異常值數(shù)目多的方向,其中 VC-LS-SVR方法估計出來的曲線被大幅度拉開遠離真實系數(shù)曲線,相比之下,其它兩種方法估計出來的曲線遠離真實系數(shù)曲線的幅度就較為小,尤其是VC-IRLS-SVR方法表現(xiàn)得最好。有趣的是,即使非對稱異常值較多的時候,三種方法基本上也都還原了系數(shù)函數(shù)的非線性趨勢,并且對β1(t)的估計都較為準確,VC-IRLS-SVR方法得到的β1(t)曲線在邊界處估計得比其它兩種方法更為準確。

圖1 誤差服從case2柯西分布時β0 (t)和β1 (t)通過三種方法估計的系數(shù)曲線

圖2 誤差服從case4非對稱污染正態(tài)分布時β0 (t)和β1 (t)通過三種方法估計的系數(shù)曲線(其中污染率δ=50%)

5 結(jié)論

本文研究了基于VC-LS-SVR應(yīng)用加權(quán)方法及迭代重加權(quán)方法獲得穩(wěn)健估計,提出了VC-WLS-SVR方法及VC-IRLS-SVR方法,既保持VC-LS-SVR模型對非線性數(shù)據(jù)的適應(yīng)性以及計算的高效性,又具有穩(wěn)健性。數(shù)值實驗結(jié)果表明,當(dāng)數(shù)據(jù)中不存在異常值時,三種估計方法的表現(xiàn)幾乎一樣好,它們都能給出相對準確和穩(wěn)定的系數(shù)估計值。當(dāng)數(shù)據(jù)中含有異常值時,VC-LS-SVR方法的Mmse(j)和SDmse(j)較大,可見回歸估計準確性和穩(wěn)定性下降,說明VC-LS-SVR方法對異常值非常敏感,即使該方法通過調(diào)節(jié)學(xué)習(xí)參數(shù)能夠增強一定的穩(wěn)健性,可以抵御少量一般異常點的影響,但當(dāng)數(shù)據(jù)包含異常值的比例增加時,其估計結(jié)果的扭曲程度迅速增大。相比較,本文提出的VC-WLS-SVR方法對異常樣本降權(quán)獲得了較高的穩(wěn)健性,即使在較多異常數(shù)據(jù)的情況下,系數(shù)估計值的均方誤差的均值和標準差都較小,而VC-IRLS-SVR方法通過使用迭代方法對極度異常的樣本點多次降權(quán),數(shù)值實驗顯示該方法更穩(wěn)健,無論數(shù)據(jù)包含一般異常值、高比例污染、存在極端異常值,還是非對稱的情況,其估計精度高,并且穩(wěn)定?？紤]到VC-IRLS-SVR方法進行迭代,需要循環(huán)重復(fù)學(xué)習(xí)樣本權(quán)重,因此比VC-WLS-SVR方法更為耗時。根據(jù)本文的數(shù)值實驗,當(dāng)異常樣本值少時,選擇VC-WLS-SVR方法就可以獲得準確且穩(wěn)健的系數(shù)函數(shù)估計,當(dāng)異常樣本數(shù)比例較高或存在極端異常值時,可選擇VC-IRLS-SVR方法擬合變系數(shù)模型,用較多計算時間獲得更為精確和穩(wěn)定的估計結(jié)果。