基于WBS2的一種多維數據多變點檢驗方法及其應用

毛佳慧，施三支

（長春理工大學 理學院，長春 130022）

變點是指在該時間點上，樣本的分布或者數字特征在突然發生變化。研究變點問題可以判斷過程中某參數發生變化的時刻并有效控制該參數，也能夠分析系統的穩定性，從外部控制變量出發檢驗或預測形態發生的變化。變點檢驗就是估計變點的數量和位置，該研究現已廣泛用于工業質量控制[1]、醫學診斷[2]、交通流[3]和網絡安全[4]等許多領域。檢驗方法也從參數檢驗發展到非參數檢驗，檢驗對象也從一維擴展到多維甚至是高維數據。關于一維數據的多變點問題，許多學者給出了檢驗效果非常好的方法，如2012年Killick[5]考慮大型數據集的多變點問題，提出PELT法來找到變點可能的數量和位置的成本函數最小值，從而得到具有計算成本的變點的最佳數量和位置，其計算效率是O()n。Rigaill[6]于2015年提出的修剪動態規劃將最佳情況下的一些基于懲罰項的變點估計方法的計算耗時加速到線性時間，并提供快速實現。2014年Frick[7]針對指數族回歸中的變點提出了控制FWER的SMUCE法。該方法首先最小化α級的一個多尺度檢驗的接受域上的變點數量，從而估計未知步長函數，再通過構建漸進未知步長函數和變點的置信集，得到估計變點位置的指數界限。Li和 Munk[8]在 2016年提出了基于 FDR的相關方法。Fryzlewicz和 Subba Rao[9]于 2014 年以及 Cho 和 Fryzlewicz[10]于 2012 年利用二元分割對單變量時間序列數據進行分割，Cho和Fryzle?wicz[11]于2016年對多變量甚至高維時間序列數據進行分割。2007年Venkatraman和Olshen[12]提出的圓形二元分割，2014 年 Fryzlewicz[13]提 出 的WBS和 2018年 Baranowski等人[14]提出的最窄閾值法旨在提高二元分割的性能。2019年Fryzle?wicz[15]在克服了WBS算法缺點的基礎上提出了WBS2（Wild Binary Segmentation2），不同于 WBS中所有子分段都是預先分割好的，WBS2算法是數據自適應的，每一組子分段的位置都是由先前檢驗到的變點的位置決定的，并使用SDLL進行模型選擇。關于高維數據變點檢驗問題的可用工具較少，并且大多不具有普適性，2019年Fryzlewicz針對英國32個城區的房價指數，利用主成分分析將一個32維、長度為284的數據降維成兩組一維數據，再分別進行變點檢驗。

通過研究交通流多維數據的特點，給出了一種基于線性投影的多維數據變點檢驗方法，能夠同時檢驗多條線或者多個站點的突變點，并對其變點位置從發生時間進行估計。通過生成階梯狀模擬數據，利用不同的變點檢驗方法進行對比分析，發現WBS2.SDLL能得到較好的檢驗結果，由于檢測的交通數據具有高度相關性，先采用主成分分析和線性投影將多維數據降為一維數據后，再利用WBS2.SDLL進行變點檢驗，檢驗結果同每個站點單獨進行多變點檢驗相比，具有快速準確的效果。再高效地找出交通高峰時間段，并推斷其發生的原因。

1 數據降維方法

首先利用主成分分析將分量相關的原隨機向量通過正交變換轉化成分量不相關的新隨機變量。具體步驟如下：

假設進行主成分分析的指標變量有m個：x1,x2,…,xm，共有n個評價對象，第i個評價對象的第j個指標的取值為aij。將各指標aij轉換成標準化 指 標

其中i=1,2,…,n;j=1,2,…,m；且

為標準化指標變量。再計算相關系數矩陣R=(rij)m×m：

式中，rii=1；rij=rji；rij是第i個指標與第j個指標的相關系數。再計算相關系數矩陣R的特征值λ1≥λ2≥…≥λm≥0，及對應的特征向量u1,u2,…,um，其中uj=(u1j,u2j,…,umj)T，由特征向量組成m個新的指標變量：

式中，y1是第1主成分；y2是第2主成分，…，ym是第m主成分。再選擇p(p≤m)個主成分，計算綜合評價值。首先計算特征值λj(j=1,2,…,m)的信息貢獻率和累積貢獻率。主成分yj的信息貢獻率為：

αp為主成分y1,y2,…,yp的累積貢獻率：

當αp接近于1(αp=0.85,0.90,0.95 )時，則選擇前p個指標變量y1,y2,…,yp作為p個主成分，代替原來m個指標變量。文章提取出累積貢獻率達到95%的前三個主成分，再把每個主成分的載荷占三個主成分載荷的和的比例值作為這三個主成分的投影系數，從而能夠通過結合主成分分析和線性投影的方法將多維數據降成一維數據。

2 WBS2.SDLL

考慮模型：

其中，ft是一個確定的一維分段持續信號，它的變點數量N和位置η1,…,ηN未知。隨機序列是相互獨立且服從均值為0且方差為σ2的正態分布，序列是有界的，即對于t=1,…,T，有|ft|0時，δT=δT，間隔長度，滿足，故對于一個足夠大的常數C，δT和B由要求。Fryzlewicz給出了用于識別候選變點位置的主要定位統計量是累積和，即對于數據(Xs,…,Xe) ：

其中，s≤b≤e，且n=e-s+1。

Fryzlewicz[15]在 2019 年提出了一種遞歸算法：WBS2。它在第一次遍歷數據時，抽取少量M個子區間樣本，起點和終點隨機選擇，并均勻且獨立地從集合{1,…,T} 中替換，WBS2將M個CUSUM的絕對值的argmax記為第一個候選變點，然后利用此候選變點將集合{1,…,T} 劃分成兩個，并再次遞歸地在候選變點的左側和右側抽取M個子區間樣本，依此類推，在短子域上所有可能子區間樣本將小于M，在這種情況下，WBS2將繪制每個這樣的子域的所有子區間樣本。

對于檢驗的置信水平為90%和95%的情況，將兩個算法分別稱為WBS2.90和WBS2.95。對降維后的一維數據，利用找到(fs,…,fe)中最有可能的一個變點位置。WBS2僅在子域長度為1時停止。完成該步驟后，WBS2按照CUSUM絕對值大小對候選變點進行降序排列。WBS2方案在計算時是數據自適應的，每個下一批M子區間的子域都是由先前檢驗到的候選變點的位置確定，由于該過程在{1,…,T}不斷變短的子域上操作，使得WBS2的計算速度加快，且其只需要設置M這一個參數值。

接下來將對模擬數據進行對比實驗，再利用WBS2對降維后的真實數據進行變點檢驗，相關算法流程圖如圖1所示。

圖1 算法流程圖

3 模擬實驗與實證分析

先利用多變點檢驗法PELT、SMUCE和WBS作為對比對象對階梯狀模擬數據進行檢驗，并從中選出表現最好的方法，再對交通行人流數據進行變點檢驗。

3.1 模擬計算

使用兩組階梯狀模擬數據，其中一組是1 000個帶有10個不同均值相同方差的服從正態分布的模擬數據，其均值分別是 1，3，5，7.5，8，6，5，5.5，6，4.5，方差取 1，其數據量均為 100，模擬次數分別為100，200，500和1 000次。第二組模擬數據與第一組的唯一區別是，其數據量是從50，140，120，75，40，160，125，80，60，150中隨機無放回地抽取。模擬次數分別是500和1 000。

令N是數據集的真實變點數量，?是使用算法估計變點數量，則表示變點估計數量的絕對平均值，它與估計的變點位置的許多精度測量值強烈相關。令?為一個分段常數函數，其每對連續估計變點之間的值是這對變點界定的區間內的數據的平均值。令為模擬中估計ft的均方誤差，作為變點位置估計誤差的度量。選用作為變點估計的準確率的評價標準，T表示模擬次數，模擬結果如表1、表2和表3所示。

表1 多維數據變點檢驗模擬結果1

表2 多維數據變點檢驗模擬結果2

表3 各組數據量不等時數據變點檢驗模擬結果

表1和表2在第一組模擬數據的基礎上得到上述結果，而表3使用的是第二組模擬數據，僅考慮模擬次數分別為500和1 000的情況。由表中結果可知，利用WBS2.95做變點檢驗得到的變點數量和變點位置的偏差都是最小的，而隨著模擬次數的增加誤差逐漸越來越小。第二組數據的檢驗誤差相較于第一組的要大，說明隨機抽取間隔數使得其中某幾組的數據量偏小，會直接影響檢驗的準確率。

為了更直觀地展示檢驗結果和所用算法的優劣性，以隨機產生的其中一組模擬數據為例，將數據和利用不同算法檢驗得到的變點個數和位置繪制在圖2中。

圖2 不同算法得到的變點個數和位置

圖2從上至下使用的變點檢驗方法分別是WBS2.95，WBS2.90，WBS，PELT和 SMUCE 法。實線表示其中一個含有10組服從不同均值相同方差的正態分布，且數據量均為100的模擬數據集，豎直的虛線表示算法檢驗出的變點位置。

圖2中使用的模擬數據集的真實變點個數為 9，變點位置分別是 100，200，300，400，500，600，700，800和900。由圖2中的豎直虛線的個數和位置可知使用的5個方法中只有WBS2.95法檢驗得到的變點個數為9，其他方法得到的變點個數均為8，盡管利用五個方法檢驗得到的變點位置與實際變點位置均有一定的偏差，但是WBS2.95法的偏差相對較小。通過圖2與表1和表2的對比結果，采用WBS2.95法對降維后的地鐵進出閘機口行人流數據進行變點檢驗。

3.2 實證分析

選取2015年4月1日上海市地鐵一號線上25個站點（富錦路，友誼西路，寶安公路，共富新村，呼蘭路，通河新村，共康路，彭浦新村，汶水路，上海馬戲城，延長路中山北路，上?；疖囌?，漢中路，新闡路，人民廣場，黃陂南路，陜西南路，常熟路，衡山路徐家匯，上海體育館，漕寶路，上海南站，錦江樂園）的閘機進出口的行人流數據，每五分鐘作為一個計數節點。為保證各維度數據量一致，且通過檢驗發現每天6：30之前和22：00之后地鐵閘機進出口的人數較少，不影響之后的檢驗結果，故將行人流數據的時間段截取為6：20-22：00，從而得到一個25維長度為188的數據集。通過運用主成分分析和線性投影將多維的數據降成1維數據后，再利用WBS2.SDLL對地鐵閘機口進站和出站的人數進行變點檢驗，得到行人進出的高峰期，通過算法檢驗得到高峰時間段的行人流均值如表4所示。

表4 多個閘機進出口高峰時間段

由表4可看出在時間段7：00-9：00和17：30-18：00內，進站閘機口人流量達到最高峰，在時間段 8：05-9：30和 17：40-18：30內，出站閘機口人流量達到最高峰。結合實際與人們通勤時間相符合。同時，為驗證上述結果的可信度，接下來將對單個站點閘機進出口行人流進行變點檢驗，選取一號線中的一個換乘站：漢中路站（可換乘12和13線）作為表4中變點檢驗結果的對比驗證對象，結果如表5所示。

表5 漢中路站閘機進出口高峰時間段

由表5可看出在時間段8：00-8：30和17：20-18：00內，漢口站進站閘機口人流量達到最高峰，在時間段8：30-9：00內，該出站閘機口人流量達到最高峰。上述時間段均被涵蓋于多個站點檢驗得出的高峰時間段內，說明通過結合主成分分析和線性投影來將數據進行降維，再對降維后的數據進行變點檢驗的方法對地鐵行人流數據是有效的，檢驗結果也是可信的。

通過對上海市一號線25個地鐵站的閘機進出口行人數量的變點檢驗分析研究，根據算法可以計算出行人流發生突變的次數和突變發生的具體時間，從而給交通管理部門對特定時間段的進出閘機口數量的合理開放提供相關數據支持。

4 結論

首先對大量階梯狀的模擬數據利用不同的的多變點檢驗算法進行變點估計，得到表現較好的WBS2.95。再利用主成分分析和線性投影將一個25維的上海市一號線地鐵站閘機進口的行人流數據有效地降成一維數據，再利用WBS2.95進行變點檢驗，結果表明7：00-8：30和17：20-19：00為進站高峰期，其中 7：46-8：20人流量最大，8：05-9：00和 17：45-20：45為出站高峰期，其中18：15-18：55人流量最大。

將WBS2.SDLL方法運用到行人交通樞紐中，同時考慮多個站點行人流的變點問題，為深入研究多個站點或多條線上的行人流突變問題提供一種新方法，但是WBS2.SDLL只適用于一維數據，之后還將嘗試研究一種不需要事先進行數據降維就能夠直接用于檢驗多維數據的方法，為多維甚至高維數據的變點檢驗問題提供參考。