基于二元分割的多變點(diǎn)估計

夏美美， 趙聯(lián)文， 范元靜， 楊 航

（西南交通大學(xué)數(shù)學(xué)學(xué)院，成都 611756）

變點(diǎn)問題一直是統(tǒng)計學(xué)中的焦點(diǎn)話題之一. 變點(diǎn)問題在1954年由Page提出，該問題的提出源于工業(yè)生產(chǎn)的質(zhì)量控制過程［1］，生產(chǎn)過程中產(chǎn)品質(zhì)量在相當(dāng)長的時間內(nèi)保持大致恒定的輸出，當(dāng)過程中某個階段出現(xiàn)故障時，產(chǎn)品質(zhì)量惡化，大部分產(chǎn)品質(zhì)量變得不可接受，此刻希望能夠停止生產(chǎn)或者發(fā)出預(yù)警，避免產(chǎn)出更多不合格的產(chǎn)品，該故障發(fā)生時刻稱為變點(diǎn). 陳希孺［2］給出了變點(diǎn)問題的一般提法，我們有一系列的觀察值（樣本），在多數(shù)情況下，這些觀察值按其出現(xiàn)的時間進(jìn)行先后排列，在某個未知的時刻，樣本的分布或其數(shù)字特征發(fā)生突然變化，這個時刻就是變點(diǎn). 由于變點(diǎn)發(fā)生時刻是未知的，因此變點(diǎn)問題研究的主要目的在于估計變點(diǎn)的個數(shù)和位置，同時圍繞變點(diǎn)估計量的相合性、漸進(jìn)分布以及變點(diǎn)的跳躍度、收斂速度等進(jìn)行理論分析.

隨著問題的深入研究，變點(diǎn)統(tǒng)計推斷已經(jīng)應(yīng)用到經(jīng)濟(jì)、金融、生物醫(yī)學(xué)、氣象學(xué)、地質(zhì)學(xué)等很多領(lǐng)域. 如金融領(lǐng)域中［3］，由于政治因素或潛在金融危機(jī)的影響，導(dǎo)致股市存在大的波動，通過對序列或過程中存在的變點(diǎn)做出相應(yīng)的估計，為金融危機(jī)等事件的監(jiān)測及預(yù)測提供有用的信息，從而幫助人們達(dá)到規(guī)避金融風(fēng)險的目的；在智能交通系統(tǒng)［4］中，通過對交通流數(shù)據(jù)的變點(diǎn)分析，實(shí)時監(jiān)控當(dāng)前和未來短時交通狀況，使車輛避開擁擠區(qū)域或使駕駛員降低行車速度，從而達(dá)到平穩(wěn)交通流、減輕區(qū)域交通壓力的目的；在地質(zhì)學(xué)中［5］，對地質(zhì)信號的變點(diǎn)分析，有助于檢測地下巖層中比較突出的斷層，從而能夠及時對地質(zhì)勘探以及地震預(yù)防等起到指導(dǎo)作用；在視頻監(jiān)控系統(tǒng)中［6］，將原始視頻分割成多個鏡頭，再對相似鏡頭進(jìn)行合并，將冗長的監(jiān)控視頻處理為簡短、準(zhǔn)確、包含視頻主要信息的視頻摘要，安保人員可以通過瀏覽視頻摘要快速做出防范，被廣泛應(yīng)用于電影編輯、道路交通、國家安防等各個方面. 因此，變點(diǎn)分析不僅有重要的理論意義，許多實(shí)際問題也可以通過轉(zhuǎn)換為變點(diǎn)分析問題進(jìn)行解決，這也是近年來變點(diǎn)問題受到統(tǒng)計學(xué)以及更多學(xué)科研究學(xué)者重視的原因.

對于單個變點(diǎn)問題的研究已經(jīng)比較成熟，但在互聯(lián)網(wǎng)時代下，數(shù)據(jù)維度、數(shù)量激增，單個變點(diǎn)問題的研究不能滿足實(shí)際應(yīng)用，因此多變點(diǎn)統(tǒng)計推斷問題亟須解決. 多變點(diǎn)估計問題不僅要估計變點(diǎn)的位置，也要確定變點(diǎn)的數(shù)目. 本文重點(diǎn)研究正態(tài)分布序列中均值多變點(diǎn)的估計問題，該問題可以簡述為：考慮模型Xi=μi+εi,i=1,2,…,T，其中μi是一維分段的恒定期望，εi～iidN(0,1)，其變點(diǎn)個數(shù)為k，變點(diǎn)位置為t1,t2,…,tk.

傳統(tǒng)的多變點(diǎn)估計問題首先對變點(diǎn)個數(shù)進(jìn)行估計，得到變點(diǎn)個數(shù)估計后，對所有可能的變點(diǎn)位置組合進(jìn)行搜索，選取擬合優(yōu)度最好的組合作為變點(diǎn)的位置估計，通常利用誤差平方和來衡量擬合優(yōu)度［7］. 在不考慮計算復(fù)雜度的情況下，傳統(tǒng)的多變點(diǎn)估計方法有較高的準(zhǔn)確性. 但當(dāng)樣本量較大時，算法有可能難以實(shí)現(xiàn). 對多變點(diǎn)估計常見的有兩種解決辦法：其一，利用變量選擇，例如文獻(xiàn)［8］將變點(diǎn)問題轉(zhuǎn)化為變量選擇問題；其二，利用二元分割思想［9］，首先估計序列中的第一個變點(diǎn)，該變點(diǎn)將序列分為變點(diǎn)前后兩段，再在前后兩段中分別重復(fù)該過程，直至序列不可繼續(xù)分段.

二元分割思想在檢測多變點(diǎn)問題中應(yīng)用廣泛，在每個階段可以看作是對單個變點(diǎn)的估計，執(zhí)行較為簡單，能同時估計出變點(diǎn)個數(shù)及其位置. 另一方面，BS 按照順序每個階段只涉及搜索單個變點(diǎn)，不會返回執(zhí)行，這意味著BS 可能不適用于兩個相鄰變點(diǎn)之間的間距較小的情況. Fryzlewicz 指出，當(dāng)任意兩個相鄰變點(diǎn)之間的最小間距大于n3/4時（n 是樣本量），BS 才是有效的. 2004 年Olshen 等改進(jìn)二元分割過程，提出了circular binary segmentation（CBS），稱為環(huán)形二元分割［10］，首先將序列的首尾連接為一個環(huán)，再對環(huán)形序列進(jìn)行二元分割，并將其應(yīng)用到染色體數(shù)據(jù)的異?；蚪M區(qū)域識別中；2014年Fryzlewicz 擴(kuò)展了二元分割方法，提出了wild binary segmentation（WBS），稱為野生二元分割，依據(jù)一定的隨機(jī)抽取機(jī)制抽取子區(qū)間，能夠改善相鄰變點(diǎn)間距較小時無法檢測的情況，使得多變點(diǎn)檢測更具有自適應(yīng)性［11］. 更多變點(diǎn)問題的統(tǒng)計研究可見參考文獻(xiàn)［12-13］.

對于正態(tài)分布的均值多變點(diǎn)模型，運(yùn)用二元分割的思想，不斷最優(yōu)化每個分段的擬合優(yōu)度，用最小二乘估計衡量擬合誤差，將經(jīng)典的最小二乘估計與二元分割思想結(jié)合，對模型中的多個變點(diǎn)進(jìn)行估計，稱其為最小二乘二元分割算法. 該改進(jìn)算法不僅保持了最小二乘估計的準(zhǔn)確性，而且能夠有效估計間隔較小的兩個相鄰變點(diǎn).

1 介紹

1.1 傳統(tǒng)變點(diǎn)估計的最小二乘法

考慮離散情況的均值變點(diǎn)模型，Xi=μi+εi,i=1,2,…,T，假定隨機(jī)誤差εi獨(dú)立同分布，且E(εi)=0，Var(εi)=σ2＜+∞，并且有μ1=…=μk1-1=b1，μk1=…=μk2-1=b2，…，μkq=…=μn=bq+1，其中1 ＜k1＜k2＜…＜kq＜T. 如果bj+1≠bj，則ki就是一個變點(diǎn)，q 為變點(diǎn)個數(shù). 具體過程如下：

第一，考慮目標(biāo)函數(shù)：

此處約定k0=1，kq+1=n+1，其中令

即觀察值的算數(shù)平均值，用來代替（1）式中的bj.

第二，得到如下只依賴于k1,…,kq的量

任意取一組初始值k1,…,kq，滿足1 ＜k1＜k2＜…＜kq＜T .

第五，這樣繼續(xù)下去得到一組新值k′1,…,k′q，將其作為初始值，回到第二步，繼續(xù)下去得到一組新值k″1,…,k″q，再回到第二步，這個過程繼續(xù)下去直到無可調(diào)整時為止. 最后所得值k?1,…,k?q，可作為k1,…,kq的估計值.

全局遍歷最小二乘中考慮變點(diǎn)個數(shù)q 的確定，有如下結(jié)論：式（3）中取q=1,2,…，則有T1≥T2≥…. 若只有q 個變點(diǎn)，則Tq與Tq+1不應(yīng)相差太多. 因此，考慮兩者之比Tq/Tq+1不小于1，可以設(shè)定一個略大于1的數(shù)，例如1.1，取最大的q，使Tq/Tq+1≥1.1，這個q 就是變點(diǎn)個數(shù)的估計.

全局遍歷搜索最小二乘法對均值變點(diǎn)的估計較為準(zhǔn)確，但需要預(yù)先假定變點(diǎn)的個數(shù)，變點(diǎn)個數(shù)的確定較為復(fù)雜.

1.2 二元分割方法

二元分割方法的步驟如下.

步驟1：檢驗無變點(diǎn)與單個變點(diǎn)，如果不存在變點(diǎn)，則停止；否則，估計第一個變點(diǎn)位置，它將整個序列分為兩部分.

步驟2：分別檢驗這兩個分段，以進(jìn)行進(jìn)一步細(xì)分.

步驟3：對每個分段重復(fù)該過程，直到無法進(jìn)一步細(xì)分為止.

2 最小二乘二元分割

2.1 模型概述

二元分割的思想即在整段時間序列中首先找到一個變點(diǎn)，變點(diǎn)將時間序列一分為二，再在兩段子序列中重復(fù)搜索過程直至結(jié)束，相當(dāng)于在每段中分別找單變點(diǎn). 考慮到最小二乘估計方法的準(zhǔn)確性，故提出最小二乘二元分割方法. 考慮離散情況的均值變點(diǎn)模型，Xi=μi+εi,i=1,2,…,T ，假定隨機(jī)誤差εi獨(dú)立同分布，且E(εi)=0，Var(εi)=σ2＜+∞，并且有μ1=…=μk1-1=b1，μk1=…=μk2-1=b2，…，μkq=…=μn=bq+1，其中1 ＜k1＜k2＜…＜kq＜T . 如果bj+1≠bj，則ki就是一個變點(diǎn)，q 為變點(diǎn)個數(shù)，這里我們假定樣本序列中至少有一個變點(diǎn)存在，即有q ≥1，考慮如下統(tǒng)計量

接下來討論有無變點(diǎn)存在的假設(shè)檢驗問題，針對上述模型給出假設(shè)檢驗如下：

第四，將得到的變點(diǎn)估計b1,b2,…,bq按照其大小順序排列后即為樣本X1,…,XT的變點(diǎn)估計值.

2.2 閾值C 的選取

假定X1,…,XT獨(dú)立同分布，有非0有限方差σ2，且有高于二階的矩存在，對于給定的S 和S*=min{S2,…,Sn} ，當(dāng)樣本大小T →∞時，有以下極限定理［14］：

2.3 統(tǒng)計量的漸進(jìn)分布

對于均值變點(diǎn)模型，若假定方差為1，則在不存在變點(diǎn)的情況下［15］，有

其中：an=(2 ln lnn)12，dn=2 ln lnn+(ln ln ln n)/2.

3 隨機(jī)模擬及實(shí)證分析

3.1 隨機(jī)模擬

為驗證最小二乘二元分割方法對正態(tài)分布多均值變點(diǎn)估計的有效性，利用統(tǒng)計軟件R 進(jìn)行模擬仿真，并將估計結(jié)果與二元分割方法估計進(jìn)行對比. 首先生成帶有四個變點(diǎn)的時間序列，變點(diǎn)位置分別為11、23、38、61，隨機(jī)數(shù)據(jù)生成如下：x=c（rnorm（10，2，1），rnorm（12，3，1），rnorm（15，5，1），rnorm（23，7，1），rnorm（20，4，1）），其中得到x 的數(shù)據(jù)圖如圖1. 對上述隨機(jī)數(shù)據(jù)重復(fù)多次實(shí)驗，針對相同的數(shù)據(jù)，分別使用二元分割方法以及最小二乘二元分割方法進(jìn)行估計，對比結(jié)果如表1所示.

圖1 模擬仿真的正態(tài)分布均值變點(diǎn)數(shù)據(jù)圖Fig.1 Data diagram of the normal distribution mean change point

表1 傳統(tǒng)二元分割估計值以及最小二乘二元分割估計值與真實(shí)變點(diǎn)的對比表Tab.1 Comparison of traditional binary segmentation estimates and least squares binary segmentation estimates with true change points

表1中b 是已知的真實(shí)變點(diǎn)，b?代表兩種方法的估計變點(diǎn)，進(jìn)而 ||b-b? 代表真實(shí)變點(diǎn)與估計變點(diǎn)的距離，距離越小則代表估計值越準(zhǔn)確. 由表1可以看出，最小二乘二元分割方法比傳統(tǒng)二元分割的估計值更準(zhǔn)確，同時在一定程度上能夠減少均方誤差.

另外，仿真模擬數(shù)據(jù)中共有400個變點(diǎn)，而傳統(tǒng)二元分割只能檢測出298個變點(diǎn)，低估了變點(diǎn)的數(shù)量，針對這一現(xiàn)象進(jìn)行詳細(xì)分析發(fā)現(xiàn)，當(dāng)在均值變化幅度較小以及變化持續(xù)較短的情況下，傳統(tǒng)二元分割方法無法準(zhǔn)確地識別變點(diǎn)，如圖2所示. 而改進(jìn)的最小二乘二元分割能夠較好地估計出此變點(diǎn)，如圖3所示. 由此可見，改進(jìn)的最小二乘二元分割能夠在均值變化幅度較小以及變化持續(xù)較短的情況下有更穩(wěn)定的估計性質(zhì).

圖2 傳統(tǒng)二元分割算法得到的變點(diǎn)分段常數(shù)擬合圖Fig.2 Variable point piecewise constant fitting graph obtained by traditional binary segmentation algorithm

圖3 最小二乘二元分割算法得到的變點(diǎn)分段常數(shù)擬合圖Fig.3 Fitting graph of the change point piecewise constant obtained by the least squares binary segmentation algorithm

3.2 實(shí)證分析

近年來，隨著社會經(jīng)濟(jì)的發(fā)展，人們的生活越來越富裕，出行的頻率也越來越高，可供出行的公共交通也越來越多. 以火車為例，為了更好地滿足旅客的出行需求，有必要對鐵路客運(yùn)量進(jìn)行分析，研究其潛在的規(guī)律，從而能夠更及時地規(guī)劃部署工作. 收集到某站點(diǎn)一年內(nèi)每天發(fā)送的旅客量數(shù)據(jù)，從歷史數(shù)據(jù)中觀察得出，每一年的數(shù)據(jù)有大致相同的變化規(guī)律. 旅客數(shù)量一般能夠維持在平穩(wěn)狀態(tài)內(nèi)，遇到節(jié)假日或者地域性活動時會打破該平穩(wěn)狀態(tài)，在該段時間內(nèi)旅客數(shù)量維持另一個平穩(wěn)狀態(tài)，該平穩(wěn)狀態(tài)的持續(xù)時間與打破平穩(wěn)狀態(tài)的原因有關(guān). 為了刻畫此規(guī)律，選取其中一年的數(shù)據(jù)進(jìn)行研究. 在整個一年數(shù)據(jù)中，工作日及周末休息日等對出行有一定的影響，因此，為了消除其對整個數(shù)據(jù)規(guī)律的影響，將全年數(shù)據(jù)分成周一到周日七個子序列，針對其中一個子序列研究其變化規(guī)律，這里選取的是周四的旅客量子序列. 將數(shù)據(jù)用散點(diǎn)圖描述，如圖4所示.

運(yùn)用最小二乘二元分割方法估計變點(diǎn)位置，得到的變點(diǎn)位置估計分別為9，24，25，27，36，37，43，44，45，51. 計算得出的各分段的均值分別為：414 647，32 471，36 235，49 911，36 343，32 007，29 623，28 011，26 612，25 277，32 263. 得到的分段常數(shù)擬合圖如圖5所示，其中細(xì)線為旅客數(shù)量變化的折線圖，粗線為使用最小二乘法計算得到的變點(diǎn)以及均值的分段常數(shù)擬合.

圖4 某站點(diǎn)周四的客運(yùn)量數(shù)據(jù)Fig.4 Passenger traffic data for a station on Thursday

圖5 最小二乘二元分割算法得到客流量分段常數(shù)擬合圖Fig.5 Fitting chart of passenger flow segmentation constant obtained by the least squares binary segmentation algorithm

4 結(jié)語

大多數(shù)檢測多變點(diǎn)的方法都是在已知變點(diǎn)個數(shù)的假設(shè)下估計變點(diǎn)位置，但是確定變點(diǎn)個數(shù)是比較困難的. 本文結(jié)合二元分割的思想，在不假定變點(diǎn)個數(shù)的情況下，實(shí)現(xiàn)了對變點(diǎn)個數(shù)以及位置的估計.

對均值多變點(diǎn)檢測，提出了改進(jìn)的二元分割最小二乘檢測方法，通過模擬以及實(shí)證，驗證了模型的合理性以及有效性.

基于二元分割最小二乘方法有如下結(jié)論：

1）得到了對于均值多變點(diǎn)模型的統(tǒng)計量；

2）給出了統(tǒng)計量的漸進(jìn)分布；

3）結(jié)合實(shí)例，進(jìn)一步驗證該方法能更容易并且快速地同時檢測出變點(diǎn)的個數(shù)及位置.

本文只考慮了將最小二乘與二元分割結(jié)合，在后續(xù)的研究中可以考慮將野二元分割與最小二乘相結(jié)合.