基于二元分割檢測(cè)分段線性趨勢(shì)中的多變點(diǎn)*

劉 偉, 胡 堯, 2**, 胡 倩

(1.貴州大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 貴陽(yáng) 550025; 2.貴州省公共大數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室, 貴陽(yáng) 550025)

0 引 言

變點(diǎn)檢測(cè)問(wèn)題一直是統(tǒng)計(jì)學(xué)中一個(gè)經(jīng)典的分支, 最初起源于Page在Biometrika上發(fā)表的一篇關(guān)于連續(xù)抽樣檢驗(yàn)的文章[1], 這是一篇關(guān)于質(zhì)量檢測(cè)問(wèn)題的理論研究, 從此開(kāi)啟了變點(diǎn)研究的篇章。在如今的大數(shù)據(jù)發(fā)展趨勢(shì)中, 掌握數(shù)據(jù)中的突變對(duì)于分析數(shù)據(jù), 挖掘其中的隱藏信息有著至關(guān)重要的作用, 所以越來(lái)越多的統(tǒng)計(jì)學(xué)者投入到這一研究領(lǐng)域[2-3]。近年來(lái), 變點(diǎn)檢測(cè)問(wèn)題已被廣泛應(yīng)用于各個(gè)領(lǐng)域, 在許多行業(yè)中, 都能看到變點(diǎn)的身影,例如自動(dòng)檢測(cè)云數(shù)據(jù)中的變點(diǎn)[4], 以保持應(yīng)用程序或網(wǎng)站的性能和可用性; 熱帶氣旋數(shù)據(jù)中的氣候變化檢測(cè)[5], 能及時(shí)預(yù)防重大自然災(zāi)害;根據(jù)光曲線數(shù)據(jù)的變化能檢測(cè)系外行星; DNA拷貝數(shù)的突變跟某些疾病的起因密切相關(guān)[6]; 對(duì)潛在協(xié)整股票價(jià)格的平穩(wěn)區(qū)間的估計(jì)可以降低損失的風(fēng)險(xiǎn)[7]等。

所謂的變點(diǎn), 就是在一個(gè)時(shí)間序列或過(guò)程中, 當(dāng)某個(gè)統(tǒng)計(jì)特性在某一時(shí)刻τ突然發(fā)生了變化, 就稱該時(shí)刻τ為變點(diǎn), 統(tǒng)計(jì)特性包括樣本的分布類型、分布參數(shù)、數(shù)字特征等, 變點(diǎn)檢測(cè)就是利用統(tǒng)計(jì)量或統(tǒng)計(jì)方法將該時(shí)刻τ估計(jì)出來(lái)。在數(shù)據(jù)被假定為分段恒定的情況下, 一類常用的方法是基于最小化成本函數(shù)的思想, 如Jackson等[8]提出的OP(Optimal Partitioning)算法, 該算法是在成本函數(shù)中引入懲罰項(xiàng), 將變點(diǎn)檢測(cè)轉(zhuǎn)化為成本函數(shù)懲罰最小化問(wèn)題, 但是在數(shù)據(jù)量比較大的情況下, 計(jì)算比較復(fù)雜；所以Killick等[9]提出基于不等式修剪的PELT(Pruned Exact Linear Time)算法, 它比OP更有效且計(jì)算簡(jiǎn)單；而Maidstone等[10]將PELT與pDPA(pruned Dynamic Programming Algorithm)相結(jié)合,提出一種更穩(wěn)健高效的FPOP(Functional Pruning Optimal Partitioning)算法等。而在數(shù)據(jù)具有線性趨勢(shì)變化的相關(guān)研究中, Bai和Perron[11]考慮通過(guò)最小二乘法估計(jì)具有多個(gè)結(jié)構(gòu)變化的線性模型, 并針對(duì)無(wú)變化的原假設(shè)提出Wald型檢驗(yàn); Kim等[12]和Tibshirani等[13]考慮了具有L1懲罰的“趨勢(shì)過(guò)濾”; Fearnhead和Maidstone等[14]通過(guò)動(dòng)態(tài)規(guī)劃算法用L0正則化來(lái)檢測(cè)斜率的變化; Spiriti等[15]研究了兩種優(yōu)化最小二乘和懲罰樣條中節(jié)點(diǎn)位置的算法; Anastasiou和Fryzlewicz提出了ID (Isolation-Detection)方法, 該方法不斷地搜索擴(kuò)展的數(shù)據(jù)段以檢測(cè)其中變化, 但正因如此, 會(huì)使得某些數(shù)據(jù)被多次重復(fù)計(jì)算, 而且每次擴(kuò)展的數(shù)據(jù)量只給出一個(gè)固定值λ=3, 并沒(méi)有說(shuō)明給出的原因, 并且在數(shù)據(jù)為長(zhǎng)時(shí)間的小跳躍情況下, 該方法比較乏力。

二元分割方法[16](Binary Segmentation, BS)是多變點(diǎn)檢測(cè)的經(jīng)典方法之一, 與其他變點(diǎn)檢測(cè)方法相比,該方法檢測(cè)效果很好, 特別是對(duì)大量數(shù)據(jù), 長(zhǎng)期性數(shù)據(jù)的多變點(diǎn)檢測(cè), 很多單變點(diǎn)檢測(cè)方法都能跟二元分割相結(jié)合而轉(zhuǎn)化為多變點(diǎn)檢測(cè), 如Olshen A B等[17]的CBS(Circular Binary Segmentation)方法,Fryzlewicz的WBS(Wild Binary Segmentation)方法[18]和WBS2(Wild Binary Segmentation 2)方法等,但這兩種方法都是用于檢測(cè)均值變點(diǎn), 所以在數(shù)據(jù)存在異常值時(shí), 檢測(cè)結(jié)果會(huì)存在很大偏差。在二元分割方法中, 檢驗(yàn)統(tǒng)計(jì)量非常重要, 所以在本文中, 根據(jù)Baranowski等[19]提出的統(tǒng)計(jì)量作為的檢驗(yàn)統(tǒng)計(jì)量, 同時(shí)受到WBS理論對(duì)整個(gè)數(shù)據(jù)區(qū)域隨機(jī)“產(chǎn)生”區(qū)間以檢測(cè)變點(diǎn)的啟發(fā),也對(duì)整個(gè)數(shù)據(jù)序列隨機(jī)抽取檢測(cè)區(qū)間進(jìn)行變點(diǎn)檢測(cè)。

1 模型及檢測(cè)方法

1.1 模型介紹

在具體應(yīng)用過(guò)程中, 由于數(shù)據(jù)類型的多樣性, 不同的數(shù)據(jù), 其分布類型不能確定, 參數(shù)方法已經(jīng)無(wú)法滿足實(shí)際應(yīng)用的需求, 然而非參數(shù)方法對(duì)總體分布的假定要求低, 不會(huì)因?yàn)閷?duì)總體分布的假定不當(dāng)而導(dǎo)致重大問(wèn)題, 更能體現(xiàn)讓數(shù)據(jù)說(shuō)話的特點(diǎn), 具有很好的穩(wěn)健性, 所以基于非參數(shù)模型對(duì)變點(diǎn)進(jìn)行研究更具有通用性。

對(duì)于觀測(cè)的數(shù)據(jù)序列Y=(Y1,Y2,…,YT), 運(yùn)用如下經(jīng)典的單變量統(tǒng)計(jì)模型:

Yt=ft+σtεt,t=1,2,…,T

(1)

其中,Yt為單次觀測(cè)數(shù)據(jù),ft是確定的數(shù)據(jù)信號(hào),εt為獨(dú)立的隨機(jī)噪聲, 且εt～N(0,σt2), 在第二節(jié)數(shù)值研究設(shè)σt=1。假設(shè)時(shí)間序列數(shù)據(jù)Y有q個(gè)變點(diǎn), 則Y被分割為q+1個(gè)不同的區(qū)間段, 記變點(diǎn)的位置分別為0=τ0<τ1<,…,<τq<τq+1=T, 其中q的值和變點(diǎn)位置τ未知, 需要估計(jì), 第i個(gè)區(qū)間段的數(shù)據(jù)為Yτi-1+1:τi=(Yτi-1+1,Yτi-1+2,…,Yτi), 本文研究的是分段線性趨勢(shì)的數(shù)據(jù), 所以此處每個(gè)區(qū)間段的數(shù)據(jù)呈線性相關(guān), 研究的目的是估計(jì)出每個(gè)區(qū)間段的首尾端點(diǎn), 也就是變點(diǎn)。

在該模型中, 設(shè)σt=σ, 則對(duì)σ的估計(jì), 可以運(yùn)用中位數(shù)絕對(duì)偏差(Median Absolute Deviation, MAD)方法[20], 在εt為獨(dú)立同分布的高斯情況下, MAD定義為

其中,Φ-1(·)表示標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)函數(shù)。注意, MAD的估計(jì)值對(duì)ft中的任何變點(diǎn)都是穩(wěn)健的, 因?yàn)樗Y(jié)合了對(duì)差異數(shù)據(jù)的處理和對(duì)中位數(shù)的使用。

1.2 檢驗(yàn)統(tǒng)計(jì)量

檢驗(yàn)統(tǒng)計(jì)量作為變點(diǎn)識(shí)別的主要部分, 其檢測(cè)能力直接影響變點(diǎn)檢測(cè)的最終結(jié)果,選擇一個(gè)好的檢驗(yàn)統(tǒng)計(jì)量至關(guān)重要, 所以選擇Baranowski在2019年所提出的統(tǒng)計(jì)量[19]為本文的檢驗(yàn)統(tǒng)計(jì)量, 具體構(gòu)造如下。

對(duì)于區(qū)間(s,e], 在給定(Ys+1,Ys+2,…,Ye)的情況下, 設(shè)l(Ys+1,Ys+2,…,Ye;Θ)為Θ的似然,Θ為參數(shù)空間, 則廣義對(duì)數(shù)似然比定義如下:

而對(duì)比函數(shù)是由數(shù)據(jù)與對(duì)比向量的內(nèi)積所構(gòu)成, 定義對(duì)比函數(shù)為

(2)

(3)

其中,

γ(s,e]=(γ(s,e](1),γ(s,e](2),…,γ(s,e](T))T,

以及常數(shù)向量

1(s,e]=1(s,e](1),1(s,e](2),…,1(s,e](T))T,

在向量

1.3 變點(diǎn)檢測(cè)方法

本文的主要思想是將數(shù)據(jù)Y=(Y1,Y2,…,YT)的局部與全局處理相結(jié)合, 將統(tǒng)計(jì)量與二元分割相結(jié)合用以檢測(cè)數(shù)據(jù)中的多變點(diǎn), 首先在全局階段,隨機(jī)繪制一些子區(qū)間(s,e], 此區(qū)間的數(shù)據(jù)為子樣本(Ye+1,Ye+2…,Ys), 1≤s

最后，運(yùn)用強(qiáng)化型施瓦茨信息準(zhǔn)則(Strengthened Schwarz Information Criterion, sSIC)對(duì)集合S集進(jìn)行最優(yōu)篩選, 得出最終變點(diǎn)。

1.4 強(qiáng)化型施瓦茨信息準(zhǔn)則

1.5 計(jì)算復(fù)雜度分析

2 模擬研究

在本節(jié)中,根據(jù)模型式(1), 產(chǎn)生兩組模擬數(shù)據(jù)以檢測(cè)本文方法,如圖1的兩幅圖所示,其中波動(dòng)較大的灰色實(shí)線(有噪音)表示所模擬的原始數(shù)據(jù)Yt,中間的白色實(shí)線表示無(wú)噪音的分段線性數(shù)據(jù)ft,與白色實(shí)線相近的黑色虛線表示的是根據(jù)本文的方法,擬合出來(lái)的分段線性數(shù)據(jù)。從圖1中可以明顯看出,本文的方法對(duì)數(shù)據(jù)轉(zhuǎn)折點(diǎn)(變點(diǎn))處的識(shí)別很精準(zhǔn),使得運(yùn)用線性回歸來(lái)擬合兩變點(diǎn)之間數(shù)據(jù)有著不錯(cuò)的效果,中間的白色實(shí)線與黑色虛線幾乎重疊。

3 實(shí)例分析

城市道路的交通狀態(tài)是交通運(yùn)行系統(tǒng)中通行能力的體現(xiàn), 交通流變點(diǎn)就是由于某種變化而導(dǎo)致的, 這種變化有可能是交通事故、自然災(zāi)害、交通管控等, 有效及時(shí)地分析出交通流的突變情況對(duì)提升道路交通的通行能力有很大幫助。

選取深圳市北環(huán)大道新洲立交的交通流卡口數(shù)據(jù)作為研究對(duì)象, 以2018-03-17(周六)和2018-03-27日(周二)00∶00—22∶00的數(shù)據(jù)為例, 每日共660個(gè)數(shù)據(jù), 對(duì)道路卡口每2 min的車流量進(jìn)行變點(diǎn)檢測(cè)(數(shù)據(jù)來(lái)源于2018年深圳杯競(jìng)賽D題)。

從圖2可知, 工作日和周六的數(shù)據(jù)分布以及變化在某些時(shí)段不太一致, 在工作日(周二), 檢測(cè)出的變點(diǎn)分別為06∶00,07∶08,13∶38,14∶24,16∶44, 對(duì)照?qǐng)D2(a), 可以得出在時(shí)間段06∶00—07∶08期間, 車流量的增量非常明顯, 在13∶38—14∶24時(shí)間段, 有短暫的增加, 在16∶44以后, 道路車流量也增加地比較明顯, 可以得出, 在工作日, 該道路的早中晚高峰期比較明顯, 這與實(shí)際情況完全吻合。

在休息日(周六), 檢測(cè)出的變點(diǎn)分別為06∶00,08∶44,12∶58,14∶42,相比工作日, 變點(diǎn)減少了一個(gè), 而且明顯可以看出, 在06∶00—08∶44和12∶58—14∶42這時(shí)間段的車流量增長(zhǎng)比工作日要緩慢許多,這是由于周末大家都沒(méi)有上班, 車流量的增長(zhǎng)速度相對(duì)工作日就比較緩慢, 而在下午卻沒(méi)有晚高峰,由于下午沒(méi)有下班的高峰期, 所以下午就沒(méi)有變點(diǎn)的出現(xiàn)。

綜上可得, 在工作日和非工作日的交通流量, 在增長(zhǎng)速度方面的差異比較大, 工作日的變化情況比較明顯, 該方法能很好地檢測(cè)出交通流中的變化情況, 如果出現(xiàn)了交通事故等問(wèn)題, 交通流量情況一定會(huì)在數(shù)據(jù)中體現(xiàn)出來(lái), 只需運(yùn)用此方法就可知道交通中出現(xiàn)的事故等,并且出行人可以合理調(diào)整自己的出行時(shí)間, 避開(kāi)出行的高峰期, 交管部門也可根據(jù)此實(shí)際情況對(duì)交通進(jìn)行有效調(diào)控。

4 結(jié)束語(yǔ)