毛文濤 蔣夢雪 李源 張仕光

在實際應用中,時間序列數(shù)據的變化往往受多種因素影響.傳統(tǒng)的時間序列預測算法主要針對影響事物變化的某一種因素進行單個時間變量的預測.而若對事物變化趨勢進行更加綜合、準確的評估,則需對多個相關因素進行同時預測,即多變量時間序列預測.相比于單變量時間序列預測,多變量時間序列預測可同時預測多個變量的走勢,同時可利用變量之間的相關信息提高動態(tài)預測的精度與穩(wěn)定性,已得到越來越多學者的關注.但是,傳統(tǒng)預測方法直接應用于多變量時間序列預測中,容易受變量之間冗余作用、誤差累積和缺乏關聯(lián)信息等特點的影響,無法取得令人滿意的預測效果.因此,選擇一個合適的建模方法對多變量進行更為準確地預測有著重要的理論價值和現(xiàn)實意義.

目前,多變量時間序列預測研究已取得一定的進展.傳統(tǒng)的預測方法將多變量分解為多個單變量,利用支持向量回歸機(Support vector machine,SVM)等方法[1]對每個變量單獨進行回歸建模.這種方法簡單明了,但是重復建模增加了計算量,同時未能有效利用變量之間的結構信息.一種典型的改進方法是利用時間序列本身的系統(tǒng)特性信息和數(shù)據特點進行預測.其中,張勇等[2]利用最大Lyapunov指數(shù)理論,選取多個鄰近參考點,實現(xiàn)了對多個混沌序列的同時預測;針對股指期貨價格預測問題,Sun等[3]利用模糊C均值對數(shù)據進行預處理,并結合粗糙集算法建立模糊邏輯關系組,實現(xiàn)對多個變量走勢的預測;韓敏等[4]利用主成分分析方法實現(xiàn)對輸入變量降維,再將動態(tài)儲備遲用作核函數(shù)充分映射多元混沌時間序列的動力學特性實現(xiàn)了避免過擬合,提高預測精度.另一種做法是圍繞著多變量輸出結構的算法改進.例如,Wang等[5]提出基于極限學習機的在線多變量時間序列預測方法,對多變量序列進行相空間重構后,建立極限學習機預測模型,實現(xiàn)在線預測;Han等[6]提出基于SCKF-γESN模型的在線多變量時間序列預測方法,利用平方根容積卡爾曼濾波算法更新γ回聲網絡參數(shù),以實現(xiàn)對多變量序列的在線預測,同時在濾波算法中加入異常值的檢測,使得預測結果更加穩(wěn)定;Chen等[7]利用K近鄰和互信息獲取多變量時序數(shù)據的重要性表示,并據此重構樣本,利用改進的加權LS-SVM進行預測.但是,這些方法大多通過多變量數(shù)據的簡單融合構建預測模型,缺乏在算法結構上針對多變量序列特點的有效改進.

由上述分析可知,對多變量時間序列預測來說,預測效果的好壞直接取決于數(shù)據中蘊含的有效信息量.而實際應用多為小樣本預測問題,若能有效利用變量之間的結構化信息,則可一定程度彌補樣本數(shù)不足帶來的信息缺失,有利于提高小樣本下預測模型的穩(wěn)定性和精度.因此,提升多變量時間序列預測效果的關鍵在于:1)如何有效挖掘序列之間的依賴關系等結構化信息?2)如何構建適用于多變量序列預測的預測模型?針對這兩點,本文提出一種基于異常序列剔除的多變量時間序列預測方法.該方法首先給出基于模糊熵的時間序列聚類算法,實現(xiàn)對相似序列的初步劃分,并引入主曲線,構建異常時間序列檢測算法,對其中的異常序列進行剔除,最后采用多維度支持向量回歸機(Multi-dimensional support vector regression,M-SVR)[8]對最終得到的相似序列進行多輸出時間序列預測.M-SVR是一種具有多輸出結構的支持向量機,不僅對小樣本有快速和準確的回歸預測效果,而且,利用超球損失函數(shù)度量多個輸出端的的風險損失,可有效利用輸出端之間的結構化信息,目前已在多步超前時間序列預測[9]等問題取得成功應用.但根據作者文獻調研,尚未發(fā)現(xiàn)M-SVR在多變量時間序列預測中的應用.鑒于此,本文采用M-SVR作為基礎建模算法,旨在利用M-SVR的結構化輸出特性,選擇具有相關性的序列同時進行預測,以達到更好的預測效果.此外,本文從理論上給出了異常序列剔除的信息損失上界和模型可靠度下界,從而證明所提算法的合理性.最后采用混沌時間序列數(shù)據與五個實際數(shù)據集數(shù)據進行仿真實驗,實驗結果驗證了所提算法的有效性.

1 相關理論

1.1 層次聚類方法

層次聚類方法是一種常用的聚類方法,分為凝聚層次聚類方法和分裂層次聚類算法.該類算法的核心思想是遞歸地采用自底向上或自頂向下策略對數(shù)據對象進行合并或分裂.分裂層次方法是把一個給定的數(shù)據對象族迭代地分裂,進而形成更小的數(shù)據對象.該類算法的優(yōu)點是能夠獲取到不同粒度的多層次聚類結構,但是算法的復雜度至少為O(n2).而凝聚層次方法從每一個數(shù)據對象開始,迭代地進行合并,形成更大的數(shù)據對象簇.目前使用的層次聚類方法多為凝聚層次聚類,中間的區(qū)別一般是定義的類間距離不同.此外,該類算法需要事先給定一個族合并或分裂的終止準則,并且某個族一旦完成合并或分裂步驟就無法撤消[10].層次聚類是一種嵌套聚類的方法,通用性強,適用于小數(shù)據集[11].

1.2 時間序列的相似性度量

時間序列的相似性度量是衡量兩個時間序列相似的標準,其有效性直接關系到后續(xù)工作的性能.目前對時間序列相似性的度量主要有兩種方法.

1)歐氏距離通過序列中對應值計算得到對應的相似性,其計算公式為

其中,i為序列中的第i個值[12].時間序列的歐氏距離需要兩個序列等長,且要求序列的值一一對應.

2)動態(tài)時間彎曲(Dynamic time warping distance,DTW)距離[13]是一種通過彎曲時間軸來更好地對時間序列形態(tài)進行匹配映射的相似性度量方法.DTW 距離在兩個時間序列之間尋找最小的映射路徑,允許一個序列中的點對應于另一個序列中的多個相鄰的點.DTW 不需要利用領域知識,只是假設相近的序列存在低消耗的平移匹配,并在時間序列分類上取得了廣泛應用[14].DTW 距離是一種基于距離測量指標的算法,其基本假設前提為數(shù)據符合正態(tài)分布.設序列s的長度是M,序列c的長度是L,兩個時間序列的DTW距離計算過程如下:首先計算矩陣中的(i,j)兩點之間的歐氏距離:d(si,cj)=‖si?cj‖2,得到一個M×L的距離矩陣,然后將兩個序列的規(guī)整路徑定義為W=(w1,···,wk),其中wk=(i,j)k=d(si,cj).DTW 的目的是找到一個從(1,1)到(M,L)使取到最小值的單調增長路徑[13].最佳路徑是γi,j=d(si,cj)+min{γi?1,j?1+γi?1,j+γi,j?1}.

1.3 主曲線

主曲線是通過數(shù)據集“中間”的光滑無參數(shù)曲線,基于一定概率分布下曲線的“自相合”特性,能將數(shù)據信息保持好,有效勾勒出原始信息的輪廓[15],目的是找到一條通過數(shù)據分布“中央”并能夠真實反映數(shù)據的形態(tài)的曲線,這意味著此曲線是數(shù)據集合的“骨架”,數(shù)據集合是這個曲線的“云”,所以主曲線對數(shù)據的信息保持性好.具體步驟如下[16]:

步驟1.令初始曲線f(j)(λ)為X的第一主成分,設j=0.

步驟2.對所有的x∈Rd,求投影指標

步驟3.定義在x上f的投影點為

步驟4.如果1?Δ(f(j+1))/Δ(f(j))小于某個閾值,則停止(其中Δ(f(j))表示點x到曲線f的歐氏平方距離),否則令j=j+1,轉步驟2.

2 基于異常時間序列檢測的多變量時間序列預測

為提高多變量時間序列預測的準確度和穩(wěn)定性,本文提出一種基于異常時間序列檢測的多變量時間序列預測方法,步驟如下:

步驟1.計算各條時間序列的模糊熵,通過比較模糊熵值判斷得到序列之間的相似度,選擇層次聚類算法進行聚類以得到初步的相似序列;

步驟2.對各類別中的序列計算主曲線,根據各個序列到主曲線的距離求得異常因子,并據此進一步剔除聚類結果中相似性相對較弱的序列;

步驟3.將得到的相似序列作為輸入,利用MSVR進行多輸出建模.

整體算法流程如圖1所示.

圖1 算法流程圖Fig.1 The flowchart of the algorithm

2.1 定義

對已知序列定義如下:{Xi},i=1,···,M,其中,Xi為{xj},j=1,···,N.M表示序列個數(shù),N代表樣本數(shù).

定義1.剔除序列過程中序列權重:構建已知序列集D的主曲線,則D中每條序列xi對應的序列權重定義為

其中,fi為序列xi到主曲線的DTW 距離.由式(1)可知,序列xi到主曲線的距離越小,其對應的序列權重越大,反之fi越大,wi越小.

定義2.序列的CP-距離:通過距離參數(shù)衡量各序列在整個序列挑選過程中的相似度.

其中,fi為序列xi到主曲線的DTW距離,oi為聚類中心,DFuzz(xi,oi)為序列到聚類中心的距離.由式(2)可知,當序列到聚類中心和主曲線的距離越小,其對應的CP-距離越小.

定義3.整個序列挑選過程中序列權重:針對聚類和主曲線的內容可知,序列權重跟序列到主曲線的DTW 距離成反比,與序列到聚類中心的距離成反比,定義如下:

其中,N和分別為序列數(shù)和序列到主曲線距離之和.

2.2 基于模糊熵的時間序列聚類算法

為有效識別序列之間的相似程度,本文首先對時間序列進行初步的相似性劃分.基本思想是利用時間序列的模糊熵的大小作為衡量序列之間相似性的標準,對時間序列進行凝聚層次聚類,建立初始模型.時間序列凝聚層次聚類將初始序列的每個對象看作一類,根據相似性度量準則計算類間距離進行合并,直至滿足終止條件.該方法不需首先確定聚類中心,可以充分利用類間距離將相對相似的序列聚成一類.模糊熵采用模糊隸屬度函數(shù)計算向量之間相似度,從而使得熵值連續(xù)平穩(wěn),對于不同的參數(shù)結果穩(wěn)定,抗噪性強[17],因此適合用做層次聚類的類間距.該算法步驟如下:

步驟1.將已知各序列視為一類,計算每兩個序列之間模糊熵的差值作為類間距離.

步驟2.根據步驟1中得到的類間距離將距離最近的兩類合并成一類.

步驟3.用average-linkage衡量兩類A與B之間的距離:

其中,1≤i,j≤M,ij,|A|和|B|分別代表類A及類B的大小.

步驟4.重復步驟2和步驟3,直到滿足事先給定的聚類類別數(shù),即聚類完成.

由上述過程可知,該時間序列聚類算法采用了模糊熵作為類間距離的衡量方式,可在整體上有效度量時間序列的相似度,從而實現(xiàn)對相似序列的初步劃分.

2.3 基于主曲線的異常序列的檢測

從實驗結果可以看出,在第2.2節(jié)得到的相似序列中,往往會因為人為設置類別初始數(shù)目等原因,將一些相似性較弱的序列強制聚到某一類中.本文將這類序列定義為聚類結果中的異常序列.異常序列的存在會減弱序列之間的相似度,從而降低M-SVR的結構化預測效果.因此,有必要對初步的聚類結果進行進一步的篩選,剔除相似性相對較弱的異常序列.

根據第1.3節(jié)分析,主曲線具有可有效勾勒出原始信息的輪廓和對數(shù)據信息保持性好的優(yōu)點[15].因此,本文在初始階段聚類的結果上,通過計算類中序列數(shù)據的均值,找到距離均值距離最近的兩條序列,求得這兩條序列的主曲線作為本類的主曲線.

序列集合中,x序列到主曲線S的可達距離定義為

其中,x(k)表示樣本中距離序列x的k近鄰樣本,‖x?S‖指從x到S的DTW 距離.基于式(4),x的局部可達密度為

x局部異常因子可定義為x(i)的局部可達密度的平均值與x局部可達密度的比,即

LOFk(x)的值越大,x的異常度就越大.當x(i)的周圍密度較高,而x的周圍密度較低時,局部異常因子變大,x被視為異常序列.相反,當x(i)的周圍密度較低,而x的周圍密度較高時,局部異常因子減小,x被視為正常序列.因此,可以根據LOFk(x)的值,從聚類結果中將異常序列進行剔除.檢測異常值在統(tǒng)計學上已有多種方法,本文借鑒置信區(qū)間的概念,采用如下方法:若LOFk(x)＞mean(X)+n×std(X),則視為異常序列,其中mean(X)為序列集合X的局部異常因子的均值,std(X)為其標準差.由于此處要剔除的通常為局部異常因子明顯較大的序列,因此在本文實驗中,n直接設為0.5即可有效識別各數(shù)據集中的異常序列.

2.4 M-SVR建模

通過以上步驟得到具有依賴關系的相似序列,根據這些序列構建輸入樣本集.

其中,n=1,2,···.相應序列的下一時刻值作為輸出Yi=(x(n+mτ)),其中m為嵌入維,τ為時延,利用M-SVR實現(xiàn)對相似序列的多輸出回歸建模.其中M-SVR優(yōu)化目標為[18]

利用一階泰勒展開式可得到其二次近似[18]:

其中k為迭代次數(shù),CT為獨立的高階項.

通過利用加權迭代最小二乘法(Iteratively reweighted least squares,IRWLS)連續(xù)迭代過程中線性搜索變量的下降方向,并收斂到最優(yōu)值.由于以上求解可轉換為對參數(shù)β和b的求解,過程如下[18]:

步驟1.設置初值k=0,βs=0,bk=0,計算和ai.

步驟2.計算βs和bs.

步驟3.采用線性搜索,計算βk+1和bk+1,并得到和ai,判斷是否收斂至最優(yōu)值.若不收斂,返回步驟2.

利用上述過程得到輸出端的參數(shù),可構建得到相應的回歸模型.

3 理論分析

本文所提方法旨在利用M-SVR的結構特性,選擇具有相似性的多變量時間序列同時進行預測,關鍵在于根據序列的異常因子大小進行異常序列的剔除,得到最具相關性的序列.為說明該方法的合理性,本文從信息熵的角度證明在異常序列檢測剔除的過程中的信息損失存在上界[19].同時,為證明本算法的有效性,本文從模型可靠度[20]入手證明在整個序列挑選的過程中模型的可靠度存在下界.

3.1 異常序列剔除過程信息損失

設剔除序列集合為Ψ={(xj,tj),j=1,2,···,M?m},其中xj對應的序列權重為wj,則剔除序列集Ψ的總體序列權重之和為

易知,表示序列集Φd中所有序列到主曲線的DTW 距離之和,對已知序列來說,為定值,故令則Ψ的總體序列損失權重之和為

定理1.令H(Ψ)表示異常序列剔除過程中的整體信息損失,那么有

其中,整體信息損失大小H(Ψ)的上界僅與剔除序列集Ψ中所有序列到主曲線的DTW距離之和有關.

證明.根據熵的定義,有

根據最大熵原理,當每一個wi都取相同的值時,H(Ψ)達到最大值.則有

由式(12)可以看出,H(Ψ)上界僅和剔除的序列到主曲線的DTW距離有關,dk越大,該界越小.

定理1從熵的角度給出異常序列剔除過程中的信息損失存在上界,在理論上證明了根據序列到主曲線的DTW 距離剔除異常序列樣本的可行性與有效性.考慮極端情況,若剔除的序列到主曲線的DTW 距離趨近于無窮大(即極端不相似),則對應的信息損失上界趨近于無窮小,這意味著該序列對整體信息幾近可忽略不計,進一步證明了本文所提異常序列檢測方法的合理性.

3.2 模型可靠度

根據上述對算法的描述可知,假設已知序列{Xi},i=1,···,M,其中,Xi為{xj},j=1,···,N.M表示序列個數(shù),N代表樣本數(shù).挑選之后序列從原來的{Xi},i=1,···,M到選取得到的{Xj},i=1,···,n.整體摒棄的序列集為{Xt},t=1,···,M?n.已知序列xi的權重為wi,損失序列集合的總權重之和為

則預測值相較于真實值的偏離率

預測值與真實值的差值

已知閾值θ的衡量預測效果,則預測值與真實值的差值E超出θ時為預測效果不理想的一類L,綜上可知,L與p呈正相關關系.

定理2.由上述描述可知,此時模型的結果服從二項分布,在給定置信度α的條件下,模型可靠度的求解公式為

則模型的可靠度RL存在下限,且僅與序列CP-距離有關.

證明.由于已知序列個數(shù)固定,聚類之后聚類中心固定,序列到對應主曲線的投影距離和為定值,已知

根據可靠度的定義式(14)可知,當α確定時,L與RL呈負相關,即p與RL呈負相關,p存在上限的同時可靠度RL存在下限,且結合式(15)可知,可靠度RL僅與序列的CP-距離有關,該值越小,L越小,可靠度RL越高,即若序列與聚類中心、主曲線的距離越近,則對應的模型可靠度越高,進一步表明本文所提算法的有效性.□

4 仿真與實驗分析

為驗證本文所提方法的有效性,分別引入混沌時間序列數(shù)據與五個實際數(shù)據集數(shù)據進行對比實驗.為方便起見,本文所提方法簡稱為OE-MSVR(Outlier eliminating multi-dimensional SVR).

為進行定量比較,本文引入均方根誤差(Root mean square error,RMSE)和平均絕對值誤差(Mean absolute error,MAE)衡量預測結果.對應的表達式為

4.1 混沌時間序列預測

一般認為,混沌時間序列廣泛存在于氣象、天文和水文等領域,具有微觀無序、宏觀有序的特點.混沌預測方法就是在相空間中找到一個非線性模型去逼近系統(tǒng)動態(tài)特性.現(xiàn)有的預測方法多為只針對單獨序列進行預測.若能利用混沌序列之間的結構特性實現(xiàn)對多條序列的同時預測,則可以提高預測效率,增強預測效果,更好地把握數(shù)據的下一步走勢.

對比算法采用超限學習機(Extreme learning machine,ELM)和LibSVM.ELM是一種單層前饋式神經網絡算法,具有良好的多輸出回歸能力[21];LibSVM是目前廣泛使用的SVM工具包,需對每個輸出單獨建模[22].實驗中,OE-MSVR和LibSVM均使用RBF核函數(shù),核參數(shù)設為5,正則化參數(shù)C均為100,松弛變量ε為0.01;ELM隱神經元個數(shù)為400,激活函數(shù)為hardlim.實驗開始前,所有的樣本均歸一化至[?1,1].

4.1.1 混沌時間序列數(shù)據生成

本文采用三種典型的混沌時間序列Lorenz[23]、Mackey-Glass[24]和Henon[25]系統(tǒng)進行實驗.

利用上述三種混沌序列構造8條時間序列.其中Lorenz序列2條,序號為1,2,分別采用參數(shù)σ=10,r=28,b=8/3,初始值為y=[5,5,15],如圖2所示.Mackey-Glass序列4條,序號為3,4,5,6,分別采用參數(shù)a=0.2,b=0.1,TAU=25,如圖3所示.Henon序列2條,序號為7,8,分別采用參數(shù)a=1.4,x0=0.03,y0=0.02,如圖4所示.

圖2 Lorenz序列數(shù)據Fig.2 The Lorenz sequences

圖3 Mackey-Glass序列數(shù)據Fig.3 The Mackey-Glass sequences

圖4 Henon序列數(shù)據Fig.4 The Henon sequences

從幾何形態(tài)看,Henon序列與其他兩種序列差別較為明顯.每一條時間序列分別產生1500個樣本,其中前1000個樣本訓練,后500個樣本用作測試.為使得實驗場景更接近真實數(shù)據,這里在數(shù)據中加入均值為0、方差為1的白噪聲.

4.1.2 實驗結果與分析

首先,通過聚類方法得到初步的相似序列集.考慮到已知序列的數(shù)量較少,為模擬實際應用場景,將聚類的初始類別數(shù)設置成兩類.利用提出的基于模糊熵的層次聚類算法進行聚類,結果如圖5所示.可以看出,序列7和序列8被聚為一類,稱為A類,其他6條序列為一類,稱為B類,與幾何圖形的相似程度一致,由此可知本文所提時間序列聚類算法結果與預期一致,可較好區(qū)分不同類型的混沌時間序列.

其次,選擇B類序列1～6進行預測.由上述描述可知,序列本身添加了白噪聲,因此可選擇主曲線描述數(shù)據的分布特性.此處采用基于主曲線的異常序列檢測算法進一步剔除序列中的異常序列(序號1,2).具體過程:計算B類序列的均值,并找到與其DTW距離最近的一條序列,定義該序列的主曲線作為B類序列的主曲線,如圖6所示.由圖6可知,主曲線能很好地反映數(shù)據分布的特性.

圖5 初始序列聚類結果Fig.5 The result of clustering on original sequences

圖6 B類序列的主曲線Fig.6 The principal curve of B class

再次,利用式(5)計算B類各序列基于主曲線的異常因子.由圖7可以看出,序列1和序列2的異常因子明顯高于其他四條序列,序列3～6的異常因子基本趨于一致.因此可將序列1和序列2從該類剔除,與生成數(shù)據時的設置完全一致.

圖7 B類各序列的異常因子Fig.7 Abnormal factor of every sequence in B class

為進一步分析本文所提方法的性能,本文對初始序列、聚類階段及異常序列剔除三個階段分別進行預測.限于篇幅,以序列3為例,預測效果如圖8所示.嵌入維m為4,時延τ為1.其中,圖8(a)是針對初始8條序列同時回歸建模的預測效果圖,圖8(b)是對已知序列進行聚類初步劃分得到的1～6序列同時回歸建模的預測效果圖,圖8(c)是在聚類初步劃分后進行異常序列剔除挑選得到的序列3～6同時回歸建模的預測效果圖.從圖8可以看出,隨著相似序列的一步步篩選,目標序列的預測效果明顯變好,同時驗證了M-SVR建模效果取決于輸出端的相似性,越是相似的序列,利用M-SVR進行預測效果越好.

表1給出了3,4,5,6四條序列在三個階段的預測效果.由表1可知,隨著各個階段序列相似程度的提高,預測精度明顯提高,表明本文所提算法可以有效提高多變量時間序列的預測精度和數(shù)值穩(wěn)定性,進一步驗證了M-SVR可有效利用多個輸出端之間的結構化特性,從而增加模型的信息含量,提高建模的精度及穩(wěn)定性.

以序列3為例,圖9～11分別給出本文所提方法,SVR和ELM的預測效果及對應的局部放大圖.

表1 各個階段預測結果性能指標對比Table 1 Prediction performance parameters of capillary of three stages

圖8 OE-MSVR在三個階段的預測效果Fig.8 The prediction of three stages with OE-MSVR algorithms

由圖9～11可知,OE-MSVR在整體上預測值與真實值最接近,除了少數(shù)的波峰和波谷外,大部分時間的預測曲線均貼近真實曲線;SVR算法的預測值與真實值的走勢基本相同,但大部分時間點的預測值均與真實值存在偏差;ELM算法的預測結果不太穩(wěn)定,相比其他兩種算法與真實值的偏差相對較大.

從表1可以看出,隨著剔除相似度較低的序列,OE-MSVR預測效果在不斷提高,表明異常序列的存在影響著數(shù)據集整體內在的關聯(lián)性,當篩選得到具有依賴關系較強相似性序列時,結構化數(shù)據中蘊含的有價值的信息得到了有效利用.結合圖9～11可知,在篩選得到的序列上進行操作時,OE-MSVR預測效果明顯好于其他兩種算法,進一步表明了OE-MSVR可以有效挖掘多變量時間序列的內在結構信息,使得預測精度更高,數(shù)值穩(wěn)定性更好.

使用OE-MSVR,SVR,ELM 對序列3～6等四條序列的預測效果對比如圖12所示.由圖12可知,本文所提算法在兩種誤差指標上明顯比另外兩種算法小,預測效果明顯較好,再次表明本文方法對多變量時間序列的預測效果好,驗證了此方法的有效性與穩(wěn)定性.對于其他類別序列,本文方法具有類似的對比效果.

4.2 實際數(shù)據集時間序列預測

為驗證算法的性能,選擇不同規(guī)模的數(shù)據集對算法進行測試,并將本文所提算法與SVM[22],ELM[21],Fol-BP[26]及AR[27]算法進行對比.其中,Fol-BP是侯公羽等于2014年提出的多變量混沌時間序列預測算法,通過RMSE和MAE進行評估.實驗設置與第4.1節(jié)相同.為消除ELM和Fol-BP算法的隨機性,其結果為重復10次的平均值.與第4.1節(jié)不同,對于實際采集到的數(shù)據,事先并沒有其內在結構的先驗知識.

實驗中,OE-MSVR均使用RBF核函數(shù),在澳門氣象數(shù)據數(shù)據集上的核參數(shù)設為2,正則化參數(shù)C為2?1,松弛變量ε為0.01;在A monitor system數(shù)據集上的核參數(shù)設為23,正則化參數(shù)C為25,松弛變量ε為0.01;在Italian air quality數(shù)據集上的核參數(shù)設為2,正則化參數(shù)C為25,松弛變量ε為0.01;在Istanbul stock exchange數(shù)據集上的核參數(shù)設為1,正則化參數(shù)C為25,松弛變量ε為0.01;在Gas sensor array drift數(shù)據集上的核參數(shù)設為22,正則化參數(shù)C為22,松弛變量ε為0.01;Fol-BP,LibSVM和ELM參數(shù)均為網絡搜索得到最優(yōu)值.實驗開始前,所有的樣本均歸一化至[?1,1].

4.2.1 數(shù)據集

選定澳門氣象數(shù)據,A monitor system,Italian air quality,Istanbul stock exchange和Gas sensor array drift五種真實的時間序列數(shù)據集對OEMSVR算法的性能進行驗證.使用的數(shù)據集可在澳門氣象局官網、UCI公共數(shù)據集下載得到.五個真實數(shù)據集信息見表2.

圖9 序列3的OE-MSVR預測效果圖Fig.9 The prediction of the third sequence with OE-MSVR algorithm

圖10 序列3的SVR預測效果圖Fig.10 The prediction of the third sequence with SVR algorithm

圖11 序列3的ELM預測效果圖Fig.11 The prediction of the third sequence with ELM algorithm

圖12 三種方法下預測誤差對比圖Fig.12 The errors of three algorithms

表2 實際數(shù)據集信息Table 2 Real datasets

本文首先選擇澳門氣象局官網提供的空氣質量數(shù)據集做進一步對比實驗.該數(shù)據包括影響空氣質量的多種因素:PM10,SO2,NO2,O3,伽馬射線(γ射線),氣壓(P),氣溫(T),濕度(H),風速(W),雨量(R),日照量(S),因此是一種典型的多變量時間序列數(shù)據.對這些指標的同時預測,具有明確的工程需求.具體為2002年～2006年1823組大譚山的氣象數(shù)據.

A monitor system,Italian air quality,Istanbul stock exchange和Gas sensor array drift數(shù)據集均來源于UCI公共數(shù)據集.其中A monitor system數(shù)據是安裝在多人室中的監(jiān)視器系統(tǒng)收集的數(shù)據,包括19個連續(xù)屬性;Italian air quality數(shù)據是在意大利部署的氣體多傳感器收集到的空氣質量數(shù)據,包括13個空氣質量屬性;Istanbul stock exchange數(shù)據集包括伊斯坦布爾證券交易所與七個其他國際指數(shù)的回報;Gas sensor array drift數(shù)據是描述在不同濃度情況下氣體傳感器的漂移數(shù)據,包括30個連續(xù)屬性.它們均屬于典型的多變量時間序列數(shù)據.

4.2.2 實驗結果

以澳門氣象數(shù)據為例,采用基于模糊熵的層次聚類對初始序列進行初步劃分,其結果如圖13所示.

從圖13可以看出,聚類將序列1,3,4,7,8,9聚為一類,稱作A類;序列2,5,6,10聚為一類,稱為B類;序列11單獨為一類,稱作C類.

圖13 初始序列聚類結果Fig.13 The result of clustering on original sequences

以A類為例,為得到相似度更高的序列,計算各序列的異常因子.構建A類序列的主曲線作為各序列異常因子檢測的基準,如圖14(a)所示;并利用式(5)計算A類各個序列基于主曲線的異常因子,如圖14(b)所示.由圖14可知,序列7和序列8的異常因子明顯較大,故而將序列7和序列8從該類剔除.

利用篩選得到的序列構建模型,用前1276組數(shù)據訓練,后547組數(shù)據進行測試驗證.圖15給出了五種算法的預測效果對比.由圖15可知,雖然所提算法在較少序列上存在欠缺,但在整體上OEMSVR所提算法預測效果最好,表明在處理實際多變量預測的問題時,OE-MSVR有效提高了多條序列的預測精度,進一步驗證了所提方法對多變量時間序列預測的有效性與穩(wěn)定性,對實際多變量預測更具有實際的應用價值.對于B類序列也有類似效果.

圖14 異常序列檢測過程圖Fig.14 The detection of abnormal sequences

圖15 五種方法下預測誤差對比圖Fig.15 The errors of five algorithms

針對A monitor system,Italian air quality,Istanbul stock exchange和Gas sensor array drift四個UCI實際數(shù)據集,其初始序列集在聚類階段的結果如圖16所示.由圖16可知,A monitor system數(shù)據集中當聚類數(shù)目設置為4類時,其中序列1,2,4,5,6,7,10,11,17,18,19為一類,稱作A類;序列3,8,9,14,15,16為一類,稱作B類;序列12和序列13分別為兩類,成為C類和D類.Italian air quality數(shù)據集中當聚類數(shù)目設置為2類時,其中序列1,3,4,11,12,13為一類,稱作A類;序列2,5,6,7,8,9,10為一類,稱作B類.Istanbul stock exchange數(shù)據集中當聚類數(shù)目設置為2類時,其中序列1,2,5,6,7為一類,稱作A類;序列3,4,8,9為一類,稱作B類;Gas sensor array drift數(shù)據集中當聚類數(shù)目設置為4類時,其中序列1,18,19,20,21,26,28,29為一類,稱作A類;序列2,3,4,5,6,10,11,12,13,14,22,27,30為一類,稱作B類;序列7,8,9,15,16,17,23,24為一類,稱作C類,25為一類,稱作D類.

針對各數(shù)據集的初始聚類結果,本文以A monitor system數(shù)據集的A類,Italian air quality數(shù)據集的B類,Istanbul stock exchange數(shù)據集的A類、Gas sensor array drift數(shù)據集的B類為例進行操作.通過構建各類序列的主曲線以衡量序列集內部各序列的異常因子,從而選取相似性高的序列集.圖17為四個數(shù)據集相應類的主曲線.

利用式(5)得到各序列基于主曲線的異常因子.從圖18可以看出,A monitor system,Italian air quality和Gas sensor array drift數(shù)據集的各序列異常因子相對明顯,Istanbul stock exchange數(shù)據的異常因子相對不明顯.

完成對各個數(shù)據集中相似序列的篩選后,利用選擇得到的序列集進行回歸建模.選擇各數(shù)據集的前70%組數(shù)據進行訓練,后30%組數(shù)據測試驗證,具體樣本數(shù)見表2.其中,OE-MSVR,Fol-BPELM和SVR四種方法的嵌入維m為2,時延τ為1,AR算法的階數(shù)為1.表3～6分別給出了A monitor system,Italian air quality,Istanbul stock exchange和Gas sensor array drift四個實際數(shù)據集在5種算法下的預測效果對比.

圖16 四個數(shù)據集的聚類結果圖Fig.16 The results of clustering on four datasets

圖17 四個數(shù)據集對應類的主曲線結果圖Fig.17 The principal curve of classes on four datasets

從表3、表4和表6可以看出,在大多數(shù)序列上,OE-MSVR均明顯優(yōu)于其他對比算法,尤其在表3的序列1,2;表4的序列2,9,10和表6的序列4,5,12,13,14上,本文所提方法均取得了較低的RMSE預測誤差.不難發(fā)現(xiàn),Fol-BP算法在數(shù)據集Italian air quality上取得了較好的整體預測效果,在序列2,7,9,10上的RMSE預測誤差均小于ELM,SVR和AR,表明較好利用多變量序列之間的結構化信息有助于降低整體的預測誤差,但由于該算法僅采用一階局域法對多變量序列進行建模,并沒有深入考慮其間的結構化信息,因此預測誤差低于OE-MSVR.我們留意到AR在表3的序列5,19等序列取得了最低的RMSE誤差,根據這些序列的圖形走勢可以看到,該序列較為平緩,波動不劇烈,AR易于取得較好的預測效果,而在序列波動相對較劇烈的表4中,AR的預測誤差則要明顯高于本文所提算法和Fol-BP,也略高于ELM,與文獻[27]的觀測結果一致.同時觀察到,ELM在少數(shù)序列上(例如表3中序列6,19和表4中序列5)的RMSE預測誤差低于OE-MSVR,但這種結果來自于大量的網格搜索后取的最優(yōu)值,同時ELM 本身也帶有多輸出的網絡結構,相比較而言,本文所提算法參數(shù)為直接指定,未做模型選擇.SVR的預測效果在所有算法中相對較差,盡管其采用了網格搜索,但由于SVR本身只能做單輸出的預測,因此在大多數(shù)序列上預測誤差較高,從另一個方面驗證了多變量時間序列預測的必要性.

圖18 四個數(shù)據集的異常因子結果圖Fig.18 The abnormal factors of four datasets

表3 A monitor system數(shù)據集預測結果性能指標對比Table 3 Prediction performance parameters of capillary of A monitor system dataset

表4 Italian air quality數(shù)據集預測結果性能指標對比Table 4 Prediction performance parameters of capillary of Italian air quality dataset

表5 Istanbul stock exchange數(shù)據集預測結果性能指標對比Table 5 Prediction performance parameters of capillary of Istanbul stock exchange dataset

表6 Gas sensor array drift數(shù)據集預測結果性能指標對比Table 6 Prediction performance parameters of capillary of Gas sensor array drift dataset

表5的結果驗證了異常序列剔除的作用.由圖18(c)可知,與澳門氣象數(shù)據(圖14(b))和Italian air quality數(shù)據(圖18(b))的異常因子分析對比,Istanbul stock exchange數(shù)據集各序列異常因子區(qū)別并不顯著.當剔除具有相對較高異常因子的序列6和序列7后,OE-MSVR雖然取得了最低的預測誤差,但是這種提高并不顯著.經過大量的網格搜索,SVR和ELM 也取得了較好的預測效果.綜合表3～5的結果,我們發(fā)現(xiàn),當異常因子存在較大差異時(如圖14(b)和圖18(b)),剔除異常序列后預測誤差有顯著下降(如圖15和表4所示);而當異常序列并不顯著時(如圖18(c)),所提算法的整體預測效果與經過模型選擇后的SVR和ELM等方法相仿,提高幅度并不明顯,而這恰恰表明了異常序列在多變量時間序列預測中的負面影響.

5 結論

基于結構化輸出的多變量時間序列預測可通過挖掘變量間蘊含的領域信息同時提高多個變量序列的預測效果.其中的關鍵問題在于如何提取變量間的依賴關系.本文提出了一種基于異常序列剔除的多變量時間序列預測方法.該方法利用基于模糊熵的層次聚類對時間序列進行初步劃分,提出了基于主曲線的異常序列檢測算法,進一步檢測并剔除異常序列,最終引入多輸出SVR進行建模和預測,同時在理論上證明了該算法的可行性與合理性,最終利用混沌時間序列數(shù)據與實際數(shù)據集數(shù)據驗證了算法的有效性.下一步的工作將集中在算法的泛化性理論分析和不同類型的變量間結構特性的建模.

References

1 Sch?lkopf B B,Smola A J.Learning with Kernels.Cambridge,Britain:MIT Press,2002,3:2165?2176

2 Zhang Yong,Guan Wei.Predication of multivariable chaotic time series based on maximal Lyapunov exponent.Acta Physica Sinica,2009,58(2):756?763(張勇,關偉.基于最大Lyapunov指數(shù)的多變量混沌時間序列預測.物理學報,2009,58(2):756?763)

3 Sun B Q,Guo H F,Karimi H R,Ge Y J,Xiong S.Prediction of stock index futures prices based on fuzzy sets and multivariate fuzzy time series.Neurocomputing,2015,151:1528?1536

4 Han Min,Xu Mei-Ling,Ren Wei-Jie.Research on multivariate chaotic time series prediction using mRSM model.Acta Automatica Sinica,2014,40(5):822?829(韓敏,許美玲,任偉杰.多元混沌時間序列的相關狀態(tài)機預測模型研究.自動化學報,2014,40(5):822?829)

5 Wang X Y,Han M.Improved extreme learning machine for multivariate time series online sequential prediction.Engineering Applications of Arti ficial Intelligence,2015,40:28?36

6 Han M,Xu M L,Liu X X,Wang X Y.Online multivariate time series prediction using SCKF-γESN model.Neurocomputing,2015,147:315?323

7 Chen T T,Lee S J.A weighted LS-SVM based learning system for time series forecasting.Information Sciences,2015,299:99?116

8 Sanchez-FernandezM,de-Prado-CumplidoM,Arenas-Garcia J,Perez-Cruz F.SVM multiregression for nonlinear channel estimation in multiple-input multiple-output systems.IEEE Transactions on Signal Processing,2005,52(8):2298?2307

9 Bao Y K,Xiong T,Hu Z Y.Multi-step-ahead time series prediction using multiple-output support vector regression.Neurocomputing,2014,129:482?493

10 Han J W,Kamber M,Pei J[Author],Fan Ming,Meng Xiao-Feng[Translator].Data Mining Concepts and Techniques(3rd edition)(Computer Science Series).Beijing:China Machine Press,2012.297?301(Han J W,Kamber M,Pei J[著],范明,孟小峰[譯].數(shù)據挖掘概念與技術(第3版)(計算機科學叢書).北京:機械工業(yè)出版社,2012.297?301)

11 Han Zhong-Ming,Chen Ni,Le Jia-Jin,Duan Da-Gao,Sun Jian-Zhi.An efficient and effective clustering algorithm for time series of hot topics.Chinese Journal of Computers,2012,35(11):2337?2347(韓忠明,陳妮,樂嘉錦,段大高,孫踐知.面向熱點話題時間序列的有效聚類算法研究.計算機學報,2012,35(11):2337?2347)

12 Lee H.The Euclidean distance degree of Fermat hypersurfaces.Journal of Symbolic Computation,2017,80:502?510

13 Hautamaki V,Nykanen P,Franti P.Time-series clustering by approximate prototypes.In:Proceedings of the 19th International Conference on Pattern Recognition.Tampa,FL,USA:IEEE,2008.1?4

14 Yang Yi-Ming,Pan Rong,Pan Jia-Lin,Yang Qiang,Li Lei.A comparative study on time series classi fication.Chinese Journal of Computers,2007,30(8):1259?1266(楊一鳴,潘嶸,潘嘉林,楊強,李磊.時間序列分類問題的算法比較.計算機學報,2007,30(8):1259?1266

15 Zhang Jun-Ping,Wang Yu.An overview of principal curves.Chinese Journal of Computers,2003,26(2):129?146(張軍平,王玨.主曲線研究綜述.計算機學報,2003,26(2):129?146)

16 Mao Wen-Tao,Wang Jin-Wan,He Ling,Yuan Pei-Yan.Hybrid sampling extreme learning machine for sequential imbalanced data.Journal of Computer Application,2015,35(8):2221?2226(毛文濤,王金婉,何玲,袁培燕.面向貫序不均衡數(shù)據的混合采樣極限學習機.計算機應用,2015,35(8):2221?2226)

17 Sun Ke-Hui,He Shao-Bo,Yin Lin-Zi,A Di-Li·Duo Li-Kun.Application of fuzzyen algorithm to the analysis of complexity of chaotic sequence.Acta Physica Sinica,2012,61(13):130507(孫克輝,賀少波,尹林子,阿地力·多力坤.模糊熵算法在混沌序列復雜度分析中的應用.物理學報,2012,61(13):130507)

18 Mao Wen-Tao,Zhao Sheng-Jie,Zhang Jun-Na.Multi-inputmulti-output support vector machine based on principal curve.Journal of Computer Application,2013,33(5):1281?1284,1293(毛文濤,趙勝杰,張俊娜.基于主曲線的多輸入多輸出支持向量機算法.計算機應用,2013,33(5):1281?1284,1293)

19 Yuan P Y,Ma H D,Fu H Y.Hotspot-entropy based data forwarding in opportunistic social networks.Pervasive and Mobile Computing,2015,16:136?154

20 Tang Li-Dong,Song Bao-Wei,Li Zheng,Zheng Ke.A fuzzy reliability evaluation method for sub-sample products based on information entropy theory.Journal of Projectiles,Rockets,Missiles and Guidance,2005,25(S1):214?216(湯禮東,宋保維,李正,鄭珂.基于信息熵理論的小子樣模糊可靠性評定方法.彈箭與制導學報,2005,25(S1):214?216)

21 Mao W T,Zhao S J,Mu X X,Wang H C.Multi-dimensional extreme learning machine.Neurocomputing,2015,149:160?170

22 Chang C C,Lin C J.LIBSVM:a library for support vector machines[Online],available:http://www.csie.ntu.edu.tw/～cjlin/libsvm/index.html.February 1,2016

23 Xu Mei-Ling,Han Min.Factor echo state network for multivariate chaotic time series prediction.Acta Automatica Sinica,2015,41(5):1042?1046(許美玲,韓敏.多元混沌時間序列的因子回聲狀態(tài)網絡預測模型.自動化學報,2015,41(5):1042?1046)

24 Li Jun,Li Da-Chao.Wind power time series prediction using optimized kernel extreme learning machine method.Acta Physica Sinica,2016,65(13):33?42(李軍,李大超.基于優(yōu)化核極限學習機的風電功率時間序列預測.物理學報,2016,65(13):33?42)

25 Ma Qian-Li,Zheng Qi-Lun,Peng Hong,Qin Jiang-Wei.Chaotic time series prediction based on fuzzy boundary modular neural networks.Acta Physica Sinica,2009,58(3):1410?1419(馬千里,鄭啟倫,彭宏,覃姜維.基于模糊邊界模塊化神經網絡的混沌時間序列預測.物理學報,2009,58(3):1410?1419)

26 Hou Gong-Yu,Liang Rong,Sun Lei,Liu Lin,Gong Yan-Fen.Risk analysis on long inclined-shaft construction in coalmine by TBM techniques based on multiple variables chaotic time series.Acta Physica Sinica,2014,63(9):90505(侯公羽,梁榮,孫磊,劉琳,龔硯芬.基于多變量混沌時間序列的煤礦斜井TBM施工動態(tài)風險預測.物理學報,2014,63(9):90505)

27 Liu C H,Shang Y L,Duan L,Chen S P,Liu C C,Chen J.Optimizing workload category for adaptive workload prediction in service clouds.Service-Oriented Computing.Lecture Notes in Computer Science.Berlin,Heidelberg,Germany:Springer,2015.87?104