基于篩選排序算法的多均值變點(diǎn)估計

李 揚(yáng)， 吳密霞，， 胡 堯， 楊 超

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

0 引言

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

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

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

1 模型及算法

1.1 模型簡介

1.2 方差σ2 的全局估計

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

將變點(diǎn)模型(1)當(dāng)作是一個非參數(shù)模型

1.3 篩選排序算法

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

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

2 數(shù)值模擬

其中n= 500，誤差εi～N(0，σ2)，噪聲參數(shù)σ= 0.1，存在6 個均值變點(diǎn)，μ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 的間距較小，使得變點(diǎn)137、298、307 較難檢測到。

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

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

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

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

圖3 交叉驗(yàn)證CV(h)圖

圖4 局部多項(xiàng)式擬合圖

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

圖5 不同趨勢下的局部多項(xiàng)式的標(biāo)準(zhǔn)差估計

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

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

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

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

3 實(shí)例分析

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

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

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

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

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

4 結(jié)論