變系數(shù)模型的穩(wěn)健LS-SVR估計算法及數(shù)值分析

張輝國,張孟娟,胡錫健

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

1 引言

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

作為目前最流行且有效的方法之一,支持向量機(Support Vector Machines, SVM)將線性思想應用于非線性數(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)化為求解線性方程組的問題,降低了算法復雜度,求解速度較快,在各個領域中都得到一定應用[12-13]和進一步的研究發(fā)展[14-16]。但LSSVM存在兩個潛在的缺陷:一是解丟失了稀疏性[17-18];二是損失函數(shù)中采用誤差平方和度量損失,若訓練數(shù)據(jù)中存在異常值以及誤差不服從高斯分布時,LSSVM的穩(wěn)健性較差。為此,國內(nèi)外許多學者針對穩(wěn)健LSSVM進行了深入研究,增加LSSVM穩(wěn)健性的方法主要分為四個方面[19]:基于異常值剔除技術改進、基于加權函數(shù)改進、基于p-范數(shù)改進以及基于損失函數(shù)改進。

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

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

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

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

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

(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)構成的n×n核矩陣,⊙表示為點乘,即每個矩陣元素對應相乘。

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

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

(3)

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

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

(4)

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

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

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

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

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

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

(5)

(6)

=1.483MAD(i)

(7)

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

算法1:VC-WLS-SVR

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

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

算法2:VC-IRLS-SVR

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

Step5:設定k=k+1;

4 數(shù)值分析

在本節(jié)中,進行了一些數(shù)值實驗,并構造合適的評價指標以評估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ù)

β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% 。

當模型誤差取自case1中的正態(tài)分布N(0,0.52)時,可認為是數(shù)據(jù)中不存在異常值的情況。當模型誤差取自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重復產(chǎn)生N次誤差項,本文設置內(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

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

一次的計算可能會具有偶然性,因此本文設置外循環(huán)次數(shù)M=50,將x重復生成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之間的非線性關系,并使用廣義交叉驗證方法[20]確定正則化參數(shù)γ以及核參數(shù)σ。由于權重的范圍為0≤vi≤1,而本文提出的迭代重加權方法是權重之間相乘以達到對異常樣本多次降權的目的,當數(shù)據(jù)存在極端異常值時,VC-IRLS-SVR方法如果在進行迭代時以滿足停止條件停止迭代的話,可能會發(fā)生所有樣本被降權接近為0的情況,因此,除了設置停止條件外,還給定了最大迭代次數(shù),在數(shù)據(jù)正?；虍惓颖据^少時設定最大迭代次數(shù)k=15,在數(shù)據(jù)存在極端異常值時設定最大迭代次數(shù)k=3。本文數(shù)值實驗結果如表1-2及圖1-2所示。

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

表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

當誤差項來自柯西分布(即數(shù)據(jù)中存在極端異常值)時,VC-LS-SVR方法估計系數(shù)函數(shù)的Mmse(j)和SDmse(j)值非常大,是VC-WLS-SVR方法的幾千倍,這是由于最小二乘方法本身就不穩(wěn)健,回歸系數(shù)函數(shù)估計結果被極端異常值扭曲,即使可以調(diào)整學習參數(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)定的解。

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

當誤差項來自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也可以看出,當數(shù)據(jù)含有極端異常值時,VC-LS-SVR方法估計已失效,估計出的曲線受異常值影響劇烈波動,而所提出的兩種穩(wěn)健方法基本可以還原真實曲線趨勢。圖2展示了當數(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 結論

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