基于篩選排序算法的多均值變點估計

李 揚， 吳密霞，， 胡 堯， 楊 超

(1. 北京工業(yè)大學(xué)理學(xué)部統(tǒng)計與數(shù)據(jù)科學(xué)系，北京 100124;2. 貴州大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，貴陽 550025; 3. 貴陽市第二中學(xué)，貴陽 550001)

0 引言

自上世紀70 年代以來，變點問題一直是統(tǒng)計中的一個熱門話題。它最早產(chǎn)生于工業(yè)質(zhì)量控制領(lǐng)域，目前在經(jīng)濟、金融、醫(yī)學(xué)、計算機等領(lǐng)域也有大量的應(yīng)用。關(guān)于單變點問題，已有一系列相當成熟的檢測方法和理論[1—7]。但這些方法較難推廣到多變點問題情形，因為多變點問題不但需要確定變點位置，更關(guān)鍵的是需要確定變點的個數(shù)。近年來，關(guān)于多變點的研究頗受到統(tǒng)計學(xué)者的廣泛關(guān)注。

關(guān)于均值多變點的研究，Yao[8]基于貝葉斯信息準則(Bayesian Information Criterion， BIC)提出了變點的數(shù)目和位置的估計方法，并證明了所得到的估計的相合性。Zhang 和Siegmund[9]基于帶有變化漂移的布朗運動模型提出修正的貝葉斯信息準則(Modified Bayesian Information Criterion， MBIC)，并給出了相應(yīng)的估計。與Yao[8]的方法相比，MBIC 通過修正懲罰項，其得到的變點個數(shù)和變點位置的估計精度更高。Harchaoui 和L′evy-Leduc[10]考慮了基于l1型懲罰準則的多變點診斷方法。另外，變點診斷的各方法的應(yīng)用離不開有效的算法。近年來，二元分割(Binary Segmentation， BS)算法[11]在多變點中廣泛應(yīng)用，其思想是遞歸地進行單變點檢測來確定所有的變點，該算法的一個缺點是檢測到的變點具有隨機性，在實際應(yīng)用中，停止判據(jù)不易計算。Fryzewicz[12]對BS 算法進行改進，提出了Wild 二元分割(Wild Binary Segmentation， WBS)算法，但WBS 算法對異常值較敏感，容易產(chǎn)生過估計問題。Niu 和Zhang[13]提出了一種快速檢測多變點的算法，即篩選排序(Screening and Ranking algorithm， SaRa)算法。該算法通過定義局部診斷函數(shù)降低運算復(fù)雜度，可適用于高維數(shù)據(jù)DNA 序列的變點檢測上。值得注意的是，SaRa 算法的閾值確定時，方差采用的是局部方差估計的最小值。因此，閾值易受局部數(shù)據(jù)的影響。結(jié)合非參數(shù)方法，本文提出了一個改進的SaRa 算法。

本文的結(jié)構(gòu)如下：在第1 節(jié)，我們主要介紹多均值模型和基于局部多項式擬合模型給出方差的全局估計，通過替換SaRa 算法中設(shè)定保守的閾值，給出了MBIC 準則下的多均值變點篩選的SaRa 改進算法；在第2 節(jié)，我們通過大量數(shù)據(jù)模擬所提出方法的有效性；在第3 節(jié)將該方法應(yīng)用于深圳市道路車流量數(shù)據(jù)的變點分析問題。

1 模型及算法

1.1 模型簡介

1.2 方差σ2 的全局估計

SaRa 算法是采用局部方差估計替換閾值表達式中的方差。下面我們基于非參的方法提供方差的全局估計。

將變點模型(1)當作是一個非參數(shù)模型

1.3 篩選排序算法

一般來講，變點周圍的數(shù)據(jù)應(yīng)該為該點是變點提供足夠的信息，而較遠的數(shù)據(jù)點并不能為該點是變點提供較多的信息?；诖?，文獻[13，15]定義了局部診斷函數(shù)來反映某位置是變點的概率，并通過計算所有位置的局部診斷函數(shù)值進行篩選排序，快速找到最有可能成為變點的候選點，具體步驟如下。

首先，定義局部診斷函數(shù)為yi的加權(quán)和函數(shù)

2 數(shù)值模擬

其中n= 500，誤差εi～N(0，σ2)，噪聲參數(shù)σ= 0.1，存在6 個均值變點，μ0=－0.18，均值向量μ= (0.18，1.07，－0.53，0.16，－0.69，－0.16)T，位置向量τ= (137，224，241，298，307，331)T，則跳躍度向量δ= (0.36，0.89，－1.6，0.69，－0.85，0.53)T，正弦趨勢參數(shù)a ∈{0，0.01，0.025}分別代表無趨勢、短趨勢和長趨勢。由于位置137 的跳躍度較小，而298 和307 的間距較小，使得變點137、298、307 較難檢測到。

圖1 是根據(jù)(15)式隨機產(chǎn)生a= 0 的模擬數(shù)據(jù)，圖2 是對應(yīng)模擬數(shù)據(jù)取d= 7 時的局部診斷函數(shù)D(x)，可以看出，在真實變點處的|D(x)|較大。這表明若在某點處局部函數(shù)的絕對值越大，此點越有可能成為變點。因此，需要通過設(shè)定一個初始閾值對樣本點進行篩選。

圖1 模擬數(shù)據(jù)及均值函數(shù)圖

圖2 局部診斷函數(shù)及變點分布

通過(13)式得到閾值λ，需先估計出?σ。圖3 展示了交叉驗證CV(h)圖，其中圖中“o”處表示最優(yōu)窗寬h0，即最小的CV值所對應(yīng)的窗寬，進而將最優(yōu)窗寬代入(4)式得到擬合值。圖4 展示了局部多項式?yi(i= 1，2，···，n)擬合圖，可以看出擬合效果較好，將所得的擬合值代入(6)式開方后得到?σ。

圖3 交叉驗證CV(h)圖

圖4 局部多項式擬合圖

對σ= 0.1 的不同趨勢的數(shù)據(jù)各模擬1 000 次，并用局部多項式結(jié)合交叉驗證得到σ的估計值?σ，將1 000 次標準差的估計值畫為箱線圖。圖5 分別展示無趨勢a=0、短趨勢a= 0.01 和長趨勢a= 0.025 下的標準差估計值。結(jié)果顯示：無論在什么趨勢下，局部多項式估計出的?σ都在0.1 附近波動，模擬效果較好。因此，可將估計的?σ代入(13)式得出初始閾值λ對數(shù)據(jù)進行篩選。

圖5 不同趨勢下的局部多項式的標準差估計

將得到的?σ代入(13)式得出初始閾值λ(令C= 2，d= 7)，若局部最大值D(x，d)＞λ，則x被選入候選池。對候選池中的點運用MBIC 中的(14)式進行最佳子集選擇，運用BS 法、WBS 方法、加入初始閾值的SaRa 方法(其中的?σ是由變點分段得出的局部常數(shù)估計得到)、本文提出的改進的SaRa 方法分別對(15)式產(chǎn)生數(shù)據(jù)進行變點檢測，并模擬1 000 次，模擬結(jié)果見表1 和表2。

表1 四種方法檢測出的變點數(shù)目(1 000 次模擬)

模擬中真實變點為6 個，由表1 可看出，對異常值敏感的BS、WBS 方法往往造成變點數(shù)目過高的估計，SaRa 算法容易造成變點數(shù)目過低的估計，而我們提出的方法無論正弦趨勢參數(shù)a設(shè)為0、0.01 還是0.025，都比BS、WBS 以及SaRa 這三種方法更精準，變點數(shù)目的估計更理想。在變點位置估計方面，分別統(tǒng)計出三種方法在1 000 次模擬中變點檢測數(shù)目為6 的各變點位置，與真實位置τ= (137，224，241，298，307，331)T作對比。這里規(guī)定若估計的位置與真實位置的距離小于2，則認為檢測正確，否則認為檢測錯誤。表2 列出三種方法檢測數(shù)目為6 的數(shù)據(jù)中的各變點檢測率和平均錯誤發(fā)現(xiàn)數(shù)(Average Falsely Discovered， AFD)。由表2 可看出，在變點檢測率與AFD 方面，雖然SP 方法與BS、WBS 這兩種方法相比效果稍差，與SaRa 方法相差不多。但SP 方法在每個變點的檢測率都能達到95%以上，且AFD 都小于1。因此，認為SP 方法能夠較精準的估計變點位置。綜上，改進的SaRa 方法在檢測變點數(shù)目和位置優(yōu)于現(xiàn)有的方法，是一種較為理想的均值變點檢測方法。

表2 四種方法的各變點檢測率及平均錯誤發(fā)現(xiàn)數(shù)(AFD)

3 實例分析

實例數(shù)據(jù)來源于深圳市道路車流量數(shù)據(jù)(http://m2ct.org/)，以北環(huán)大道新洲立交東往西方向的卡口為例，分別選取2018 年3 月14 日(周三，工作日時間)和2018 年3 月18 日(周日，非工作日時間)這兩天的車流量數(shù)據(jù)。數(shù)據(jù)結(jié)構(gòu)為每天00:00～22:00 的每兩分鐘(共660 個)過車數(shù)。

采用基于MBIC 的篩選排序算法進行車流量的均值變點估計。分別運用SP 方法和WBS 方法進行變點檢測，變點位置估計結(jié)果如圖6 和圖7 所示。由于本文的SP 方法需要估計方差得出?σ，這里運用局部多項式對原始數(shù)據(jù)進行擬合。由圖6 和圖7 可看出，局部多項式擬合效果較好，為確定初始閾值奠定良好的基礎(chǔ)，便于篩選候選點。然后，經(jīng)過MBIC 的最優(yōu)子集選擇得出最終變點。

圖6 北環(huán)大道新洲立交東往西方向2018 年3 月14 日車流量變點位置估計

圖7 北環(huán)大道新洲立交東往西方向2018 年3 月18 日車流量變點位置估計

在2018 年3 月14 日(周三，工作日時間)的變點估計中，SP 方法和WBS 方法檢測結(jié)果相近，分別在早高峰(SP 方法6:56，WBS 方法6:50)和晚高峰(WBS 方法18:40，WBS方法18:56)檢測出均值變化，與實際相吻合。在2018 年3 月18 日(周日，非工作日時間)的變點位置估計中，早高峰(兩種方法檢測都為7:38)比3 月14 日(工作日)稍晚，這一點與人們想在休息日多休息一下的實際情況相符，二者晚高峰檢測效果也接近(SP 方法18:42，WBS 方法18:36)，但二者在早晚高峰期之間各檢測出一個變點(SP 方法13:52，WBS 方法9:32)，由于研究卡口附近有公園、商業(yè)街等人流量較大的區(qū)域，在周末的13:52 左右人們可能去購物或游玩，因而出現(xiàn)新的“小高峰”。事實上，圖6 的13:52 也能看出存在均值跳躍，而在9:32 并未有均值的跳躍。因此，認為9:32 檢測效果不理想。可能是WBS 方法對異常值敏感造成的過估計，與表2 的模擬結(jié)果相符。

綜上，研究區(qū)域在工作日和非工作日的交通流狀態(tài)不盡相同，人們可以由此合理規(guī)劃出行時間，盡量避開高峰時間。檢測結(jié)果也進一步驗證深圳市對外地車限制行駛政策，即工作日早晚高峰(早7:00～9:00，晚17:30～19:30)時間禁止外地車通行。此外，本文的SP 方法檢測均值變點可應(yīng)用于城市中不同路段，相關(guān)部門可根據(jù)不同方向和路段發(fā)生變點的時刻進行及時調(diào)控，實時為出行人群提供合理的選擇和建議，在一定程度上利于交通管理，避免造成較大的擁堵。

4 結(jié)論