基于LS-SVM的時空變系數回歸模型算法及仿真

蘇磊·乃比，張輝國，胡錫健

(新疆大學數學與系統科學學院，新疆烏魯木齊市 830046)

1 引言

SVM首先由Vapnik等[7]提出，其主要特點為可以通過引入核函數估計較為復雜的非線性變化趨勢，且比起Bagging，嶺回歸及回歸樹，SVM具有稀疏化運算，避免維數災難的優點。Suykens等[8]提出的最小二乘支持向量機(least squares support vector machine，LS-SVM)則是將SVM的二次規劃問題轉化成線性方程組，大幅度提高了求解的可行性。Quan等[9]建立了一類針對非線性時間序列的加權LS-SVM，并通過仿真模擬證實了其精度高于普通的LS-SVM。Shim和Hwang[10]提出了變系數模型的LS-SVM回歸估計方法，并通過三組仿真模擬及一項實證分析證實其處理帶有不同變化趨勢的回歸系數及誤差項具有精度高，偏差小的優點。隨后Shim和Hwang[11]基于GTWAR模型提出了一類基于SVM核函數的時空加權自回歸模型，并通過深圳房價數據將提出的模型與普通最小二乘回歸，空間自回歸，GWR，GTWAR等模型比較分析，驗證了其具有更低的誤差，更高的F統計量以更高的決策系數R2。

然而絕大部分空間或時空模型的研究很少涉及將傳統回歸框架與新興的機器學習框架結合并通過仿真驗證其還原真實系數曲面的能力。本文通過將LS-SVM框架應用于時空變系數回歸模型，在保持傳統模型的可解釋性特征的前提下提高模型精度及靈活性。

2 GTWR模型的介紹

(1)

其中，βk(ui，vi，ti)，k=1，…，p為 (ui，vi，ti)點處的系數函數，εi(i=1，…，n)為均值為0，方差為σ2的獨立同分布隨機變量。模型允許通過規定xi1=(1，1，…，1)T引入時空變系數截距項β1(ui，vi，ti)。系數函數βj(ui，vi，ti)的估計值可由下式給出：

其中X=(x1，…，xp)，Y=(y1，…，yn)T，時空權重矩陣W(ui，vi，ti)=diag(wi1，…，win)，其對角元素wij(j=1，…，n)為樣本點之間的權重函數，可以通過引入以下時空距離確定該權重函數：

(2)

相應的

(3)

GTWR模型中的超參數參數一般可通過交叉驗證(CV)或赤池信息準則(AIC)等方法選取最有超參數。

3 LSSVM-STVC回歸模型

3.1 LSSVM-STVC回歸模型及其系數估計

基于文獻[10]中的思想，本文假設模型(1)中的時空變系數函數βk(ui，vi，ti)，k=1，…，p與時空坐標(ui，vi，ti)有如下非線性關系

(4)

其中，Φ(ui，vi，ti)為非線性特征映射，其作用是將輸入空間映射到更高維(甚至無窮維)特征空間；wj為與特征映射Φ(ui，vi，ti)相應維度的權重向量。于是對于給定樣本點(ui，vi，ti，xi)，回歸函數可以表示為

根據Mercer[12]，特征映射的內積可以由輸入空間的一類核函數表示：

K((ui，vi，ti)，(uj，vj，tj))=Φ(ui，vi，ti)TΦ(uj，vj，tj)

(5)

最常用的核函數之一稱之為高斯核函數，其形式如下

Kij=K((ui，vi，ti)，(uj，vj，tj))

基于經典LS-SVM模型，模型(1)對應的優化問題可定義為

i=1，…，n

其中，γ為正則化參數，ei,i=1，…，n為均值為0，方差Var(e)<∞的獨立同分布隨機誤差。通過引入一組拉格朗日乘子αi，i=1，…，n可得以下拉格朗日函數

+ei-yi)，i=1，…，n

(6)

于是，原優化問題的KKT條件可以表示為

k=1，…，p，

基于上述，替換wk,k=1，…，p及ei,i=1，…，n，可得以下線性方程組

(7)

其中，X=(x1，…，xp)，Y=(y1，…，yn)T，α=(α1，…，αn)T，b=(b1，…，bj)T，核函數矩陣K由Kij=K((ui，vi，ti)，(uj，vj，tj))組成，Op×p表示p維全0方陣，0p×1為p維全0向量，⊙表示矩陣對應元素相乘。

給定樣本點(u0，v0，t0，x0)，LSSVM-STVC回歸模型的系數估計為

對應的回歸函數的估計值為

3.2 LSSVM-STVC回歸模型的參數選擇

基于文獻[10]中模型選擇的思想，本文給出一類適用于LSSVM-STVC回歸模型的一類廣義交叉驗證(GCV)方法。該方法相比較于傳統的交叉驗證(CV)可大幅度提高運算效率。

假設λ是LSSVM-STVC回歸模型的一組超參數(包括核參數κ，τ以及正則化參數γ)，且

=((u1，v1，t1，x1)|λ)，…，(un，vn，tn，xn)|λ))T=HY

其中，H=(XXT⊙K)H0稱之為帽子矩陣，且有

及Z=XXT⊙K+(1γ)In，于是可得相應的GCV函數

最優超參數組λ亦選取恰當使得上式取最小值即可。

4 仿真與結果分析

基于本文模型，本節分別設計樣本量為n=245及n=2205情形下的仿真，并通過構造合適的統計量對擬合效果進行檢驗。

給定時空變系數模型(1)，令p=2，xi1≡1 及xi2=xi，i=1，…，n，于是仿真模擬數據可由下式生成

yi=β1(ui，vi，ti)+β2(ui，vi，ti)xi+εi，

i=1，…，n

(8)

其中，解釋變量xi，i=1，…，n通過均勻分布U(0，1)隨機生成，隨機誤差項εi，i=1，…，n由正態分布N(0，0.52)隨機生成，回歸系數βk(u，v，t)，k=1，2分別取定為以下時空非線性及時空線性趨勢：

β1(u，v，t)=2sin(πu)exp(-(4-t)/40)

β2(u，v，t)=(u+v)exp(-(4-t)/40)，

其中，對應于每一個時間點t，均有1×1矩形區域內的m×m個空間坐標u和v。根據上式，回歸系數βk(u，v，t)，k=1，2值均限定在0到2之間。由上述步驟生成的兩個回歸系數的初始曲面如圖1所示。

為了評價本文模型的擬合優度，以下引入幾類統計量。假設響應變量的估計值為

其中H為上一節定義的帽子矩陣。于是殘差的樣本方差可定義為：

(9)

(10)

圖1 初始回歸系數曲面

對于給定時空坐標點(ui，vi，ti)，第k個回歸系數βk(ui，vi，ti)的估計值為：

(11)

(12)

雖然實際問題中真實的βk(ui，vi，ti)往往不是已知的，進而導致MMSEβk無法計算，但是本文設計的仿真模擬實驗可以有效地展示在一組同時包含時空線性以及非線性趨勢的回歸系數下的擬合效果。

實驗分別選取樣本量為n=245及n=2205的情形來展示本文模型在樣本量較小及較大情形下的擬合效果。對于每一個樣本量，解釋變量x重復隨機生成50次，且對于每一次生成的x，誤差項隨機生成N=100次并計算(8)-(11)式中各統計量。最終結果如表1及圖2、3所示。

圖2 重復試驗50次MMSEβ1(左)與MMSEβ2(右)散點圖

圖3 重復試驗50次散點圖

表1 50次重復試驗MMSEβk(k=1，2)以及的均值及標準差

圖4 n=245下的m(1 (ui，vi，ti))(左)及m(2 (ui，vi，ti))(左)曲面

圖5 n=2205下的m(1 (ui，vi，ti))(左)及m(2 (ui，vi，ti))(左)曲面

5 結論