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

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

(1. 北京工業大學理學部統計與數據科學系，北京 100124;2. 貴州大學數學與統計學院，貴陽 550025; 3. 貴陽市第二中學，貴陽 550001)

0 引言

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

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

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

1 模型及算法

1.1 模型簡介

1.2 方差σ2 的全局估計

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

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

1.3 篩選排序算法

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

首先，定義局部診斷函數為yi的加權和函數

2 數值模擬

其中n= 500，誤差εi～N(0，σ2)，噪聲參數σ= 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，正弦趨勢參數a ∈{0，0.01，0.025}分別代表無趨勢、短趨勢和長趨勢。由于位置137 的跳躍度較小，而298 和307 的間距較小，使得變點137、298、307 較難檢測到。

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

圖1 模擬數據及均值函數圖

圖2 局部診斷函數及變點分布

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

圖3 交叉驗證CV(h)圖

圖4 局部多項式擬合圖

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

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

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

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

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

表2 四種方法的各變點檢測率及平均錯誤發現數(AFD)

3 實例分析

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

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

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

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

在2018 年3 月14 日(周三，工作日時間)的變點估計中，SP 方法和WBS 方法檢測結果相近，分別在早高峰(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)，由于研究卡口附近有公園、商業街等人流量較大的區域，在周末的13:52 左右人們可能去購物或游玩，因而出現新的“小高峰”。事實上，圖6 的13:52 也能看出存在均值跳躍，而在9:32 并未有均值的跳躍。因此，認為9:32 檢測效果不理想。可能是WBS 方法對異常值敏感造成的過估計，與表2 的模擬結果相符。

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

4 結論