毛佳慧,施三支
(長春理工大學(xué) 理學(xué)院,長春 130022)
變點(diǎn)是指在該時(shí)間點(diǎn)上,樣本的分布或者數(shù)字特征在突然發(fā)生變化。研究變點(diǎn)問題可以判斷過程中某參數(shù)發(fā)生變化的時(shí)刻并有效控制該參數(shù),也能夠分析系統(tǒng)的穩(wěn)定性,從外部控制變量出發(fā)檢驗(yàn)或預(yù)測形態(tài)發(fā)生的變化。變點(diǎn)檢驗(yàn)就是估計(jì)變點(diǎn)的數(shù)量和位置,該研究現(xiàn)已廣泛用于工業(yè)質(zhì)量控制[1]、醫(yī)學(xué)診斷[2]、交通流[3]和網(wǎng)絡(luò)安全[4]等許多領(lǐng)域。檢驗(yàn)方法也從參數(shù)檢驗(yàn)發(fā)展到非參數(shù)檢驗(yàn),檢驗(yàn)對象也從一維擴(kuò)展到多維甚至是高維數(shù)據(jù)。關(guān)于一維數(shù)據(jù)的多變點(diǎn)問題,許多學(xué)者給出了檢驗(yàn)效果非常好的方法,如2012年Killick[5]考慮大型數(shù)據(jù)集的多變點(diǎn)問題,提出PELT法來找到變點(diǎn)可能的數(shù)量和位置的成本函數(shù)最小值,從而得到具有計(jì)算成本的變點(diǎn)的最佳數(shù)量和位置,其計(jì)算效率是O()n。Rigaill[6]于2015年提出的修剪動(dòng)態(tài)規(guī)劃將最佳情況下的一些基于懲罰項(xiàng)的變點(diǎn)估計(jì)方法的計(jì)算耗時(shí)加速到線性時(shí)間,并提供快速實(shí)現(xiàn)。2014年Frick[7]針對指數(shù)族回歸中的變點(diǎn)提出了控制FWER的SMUCE法。該方法首先最小化α級的一個(gè)多尺度檢驗(yàn)的接受域上的變點(diǎn)數(shù)量,從而估計(jì)未知步長函數(shù),再通過構(gòu)建漸進(jìn)未知步長函數(shù)和變點(diǎn)的置信集,得到估計(jì)變點(diǎn)位置的指數(shù)界限。Li和 Munk[8]在 2016年提出了基于 FDR的相關(guān)方法。Fryzlewicz和 Subba Rao[9]于 2014 年以及 Cho 和 Fryzlewicz[10]于 2012 年利用二元分割對單變量時(shí)間序列數(shù)據(jù)進(jìn)行分割,Cho和Fryzle?wicz[11]于2016年對多變量甚至高維時(shí)間序列數(shù)據(jù)進(jìn)行分割。2007年Venkatraman和Olshen[12]提出的圓形二元分割,2014 年 Fryzlewicz[13]提 出 的WBS和 2018年 Baranowski等人[14]提出的最窄閾值法旨在提高二元分割的性能。2019年Fryzle?wicz[15]在克服了WBS算法缺點(diǎn)的基礎(chǔ)上提出了WBS2(Wild Binary Segmentation2),不同于 WBS中所有子分段都是預(yù)先分割好的,WBS2算法是數(shù)據(jù)自適應(yīng)的,每一組子分段的位置都是由先前檢驗(yàn)到的變點(diǎn)的位置決定的,并使用SDLL進(jìn)行模型選擇。關(guān)于高維數(shù)據(jù)變點(diǎn)檢驗(yàn)問題的可用工具較少,并且大多不具有普適性,2019年Fryzlewicz針對英國32個(gè)城區(qū)的房價(jià)指數(shù),利用主成分分析將一個(gè)32維、長度為284的數(shù)據(jù)降維成兩組一維數(shù)據(jù),再分別進(jìn)行變點(diǎn)檢驗(yàn)。
通過研究交通流多維數(shù)據(jù)的特點(diǎn),給出了一種基于線性投影的多維數(shù)據(jù)變點(diǎn)檢驗(yàn)方法,能夠同時(shí)檢驗(yàn)多條線或者多個(gè)站點(diǎn)的突變點(diǎn),并對其變點(diǎn)位置從發(fā)生時(shí)間進(jìn)行估計(jì)。通過生成階梯狀模擬數(shù)據(jù),利用不同的變點(diǎn)檢驗(yàn)方法進(jìn)行對比分析,發(fā)現(xiàn)WBS2.SDLL能得到較好的檢驗(yàn)結(jié)果,由于檢測的交通數(shù)據(jù)具有高度相關(guān)性,先采用主成分分析和線性投影將多維數(shù)據(jù)降為一維數(shù)據(jù)后,再利用WBS2.SDLL進(jìn)行變點(diǎn)檢驗(yàn),檢驗(yàn)結(jié)果同每個(gè)站點(diǎn)單獨(dú)進(jìn)行多變點(diǎn)檢驗(yàn)相比,具有快速準(zhǔn)確的效果。再高效地找出交通高峰時(shí)間段,并推斷其發(fā)生的原因。
首先利用主成分分析將分量相關(guān)的原隨機(jī)向量通過正交變換轉(zhuǎn)化成分量不相關(guān)的新隨機(jī)變量。具體步驟如下:
假設(shè)進(jìn)行主成分分析的指標(biāo)變量有m個(gè):x1,x2,…,xm,共有n個(gè)評價(jià)對象,第i個(gè)評價(jià)對象的第j個(gè)指標(biāo)的取值為aij。將各指標(biāo)aij轉(zhuǎn)換成標(biāo)準(zhǔn)化 指 標(biāo)
其中i=1,2,…,n;j=1,2,…,m;且
為標(biāo)準(zhǔn)化指標(biāo)變量。再計(jì)算相關(guān)系數(shù)矩陣R=(rij)m×m:
式中,rii=1;rij=rji;rij是第i個(gè)指標(biāo)與第j個(gè)指標(biāo)的相關(guān)系數(shù)。再計(jì)算相關(guān)系數(shù)矩陣R的特征值λ1≥λ2≥…≥λm≥0,及對應(yīng)的特征向量u1,u2,…,um,其中uj=(u1j,u2j,…,umj)T,由特征向量組成m個(gè)新的指標(biāo)變量:
式中,y1是第1主成分;y2是第2主成分,…,ym是第m主成分。再選擇p(p≤m)個(gè)主成分,計(jì)算綜合評價(jià)值。首先計(jì)算特征值λj(j=1,2,…,m)的信息貢獻(xiàn)率和累積貢獻(xiàn)率。主成分yj的信息貢獻(xiàn)率為:
αp為主成分y1,y2,…,yp的累積貢獻(xiàn)率:
當(dāng)αp接近于1(αp=0.85,0.90,0.95 )時(shí),則選擇前p個(gè)指標(biāo)變量y1,y2,…,yp作為p個(gè)主成分,代替原來m個(gè)指標(biāo)變量。文章提取出累積貢獻(xiàn)率達(dá)到95%的前三個(gè)主成分,再把每個(gè)主成分的載荷占三個(gè)主成分載荷的和的比例值作為這三個(gè)主成分的投影系數(shù),從而能夠通過結(jié)合主成分分析和線性投影的方法將多維數(shù)據(jù)降成一維數(shù)據(jù)。
考慮模型:
其中,ft是一個(gè)確定的一維分段持續(xù)信號,它的變點(diǎn)數(shù)量N和位置η1,…,ηN未知。隨機(jī)序列是相互獨(dú)立且服從均值為0且方差為σ2的正態(tài)分布,序列是有界的,即對于t=1,…,T,有|ft|0時(shí),δT=δT,間隔長度,滿足,故對于一個(gè)足夠大的常數(shù)C,δT和B由要求。Fryzlewicz給出了用于識別候選變點(diǎn)位置的主要定位統(tǒng)計(jì)量是累積和,即對于數(shù)據(jù)(Xs,…,Xe) :
其中,s≤b≤e,且n=e-s+1。
Fryzlewicz[15]在 2019 年提出了一種遞歸算法:WBS2。它在第一次遍歷數(shù)據(jù)時(shí),抽取少量M個(gè)子區(qū)間樣本,起點(diǎn)和終點(diǎn)隨機(jī)選擇,并均勻且獨(dú)立地從集合{1,…,T} 中替換,WBS2將M個(gè)CUSUM的絕對值的argmax記為第一個(gè)候選變點(diǎn),然后利用此候選變點(diǎn)將集合{1,…,T} 劃分成兩個(gè),并再次遞歸地在候選變點(diǎn)的左側(cè)和右側(cè)抽取M個(gè)子區(qū)間樣本,依此類推,在短子域上所有可能子區(qū)間樣本將小于M,在這種情況下,WBS2將繪制每個(gè)這樣的子域的所有子區(qū)間樣本。
對于檢驗(yàn)的置信水平為90%和95%的情況,將兩個(gè)算法分別稱為WBS2.90和WBS2.95。對降維后的一維數(shù)據(jù),利用找到(fs,…,fe)中最有可能的一個(gè)變點(diǎn)位置。WBS2僅在子域長度為1時(shí)停止。完成該步驟后,WBS2按照CUSUM絕對值大小對候選變點(diǎn)進(jìn)行降序排列。WBS2方案在計(jì)算時(shí)是數(shù)據(jù)自適應(yīng)的,每個(gè)下一批M子區(qū)間的子域都是由先前檢驗(yàn)到的候選變點(diǎn)的位置確定,由于該過程在{1,…,T}不斷變短的子域上操作,使得WBS2的計(jì)算速度加快,且其只需要設(shè)置M這一個(gè)參數(shù)值。
接下來將對模擬數(shù)據(jù)進(jìn)行對比實(shí)驗(yàn),再利用WBS2對降維后的真實(shí)數(shù)據(jù)進(jìn)行變點(diǎn)檢驗(yàn),相關(guān)算法流程圖如圖1所示。
圖1 算法流程圖
先利用多變點(diǎn)檢驗(yàn)法PELT、SMUCE和WBS作為對比對象對階梯狀模擬數(shù)據(jù)進(jìn)行檢驗(yàn),并從中選出表現(xiàn)最好的方法,再對交通行人流數(shù)據(jù)進(jìn)行變點(diǎn)檢驗(yàn)。
使用兩組階梯狀模擬數(shù)據(jù),其中一組是1 000個(gè)帶有10個(gè)不同均值相同方差的服從正態(tài)分布的模擬數(shù)據(jù),其均值分別是 1,3,5,7.5,8,6,5,5.5,6,4.5,方差取 1,其數(shù)據(jù)量均為 100,模擬次數(shù)分別為100,200,500和1 000次。第二組模擬數(shù)據(jù)與第一組的唯一區(qū)別是,其數(shù)據(jù)量是從50,140,120,75,40,160,125,80,60,150中隨機(jī)無放回地抽取。模擬次數(shù)分別是500和1 000。
令N是數(shù)據(jù)集的真實(shí)變點(diǎn)數(shù)量,?是使用算法估計(jì)變點(diǎn)數(shù)量,則表示變點(diǎn)估計(jì)數(shù)量的絕對平均值,它與估計(jì)的變點(diǎn)位置的許多精度測量值強(qiáng)烈相關(guān)。令?為一個(gè)分段常數(shù)函數(shù),其每對連續(xù)估計(jì)變點(diǎn)之間的值是這對變點(diǎn)界定的區(qū)間內(nèi)的數(shù)據(jù)的平均值。令為模擬中估計(jì)ft的均方誤差,作為變點(diǎn)位置估計(jì)誤差的度量。選用作為變點(diǎn)估計(jì)的準(zhǔn)確率的評價(jià)標(biāo)準(zhǔn),T表示模擬次數(shù),模擬結(jié)果如表1、表2和表3所示。
表1 多維數(shù)據(jù)變點(diǎn)檢驗(yàn)?zāi)M結(jié)果1
表2 多維數(shù)據(jù)變點(diǎn)檢驗(yàn)?zāi)M結(jié)果2
表3 各組數(shù)據(jù)量不等時(shí)數(shù)據(jù)變點(diǎn)檢驗(yàn)?zāi)M結(jié)果
表1和表2在第一組模擬數(shù)據(jù)的基礎(chǔ)上得到上述結(jié)果,而表3使用的是第二組模擬數(shù)據(jù),僅考慮模擬次數(shù)分別為500和1 000的情況。由表中結(jié)果可知,利用WBS2.95做變點(diǎn)檢驗(yàn)得到的變點(diǎn)數(shù)量和變點(diǎn)位置的偏差都是最小的,而隨著模擬次數(shù)的增加誤差逐漸越來越小。第二組數(shù)據(jù)的檢驗(yàn)誤差相較于第一組的要大,說明隨機(jī)抽取間隔數(shù)使得其中某幾組的數(shù)據(jù)量偏小,會(huì)直接影響檢驗(yàn)的準(zhǔn)確率。
為了更直觀地展示檢驗(yàn)結(jié)果和所用算法的優(yōu)劣性,以隨機(jī)產(chǎn)生的其中一組模擬數(shù)據(jù)為例,將數(shù)據(jù)和利用不同算法檢驗(yàn)得到的變點(diǎn)個(gè)數(shù)和位置繪制在圖2中。
圖2 不同算法得到的變點(diǎn)個(gè)數(shù)和位置
圖2從上至下使用的變點(diǎn)檢驗(yàn)方法分別是WBS2.95,WBS2.90,WBS,PELT和 SMUCE 法。實(shí)線表示其中一個(gè)含有10組服從不同均值相同方差的正態(tài)分布,且數(shù)據(jù)量均為100的模擬數(shù)據(jù)集,豎直的虛線表示算法檢驗(yàn)出的變點(diǎn)位置。
圖2中使用的模擬數(shù)據(jù)集的真實(shí)變點(diǎn)個(gè)數(shù)為 9,變點(diǎn)位置分別是 100,200,300,400,500,600,700,800和900。由圖2中的豎直虛線的個(gè)數(shù)和位置可知使用的5個(gè)方法中只有WBS2.95法檢驗(yàn)得到的變點(diǎn)個(gè)數(shù)為9,其他方法得到的變點(diǎn)個(gè)數(shù)均為8,盡管利用五個(gè)方法檢驗(yàn)得到的變點(diǎn)位置與實(shí)際變點(diǎn)位置均有一定的偏差,但是WBS2.95法的偏差相對較小。通過圖2與表1和表2的對比結(jié)果,采用WBS2.95法對降維后的地鐵進(jìn)出閘機(jī)口行人流數(shù)據(jù)進(jìn)行變點(diǎn)檢驗(yàn)。
選取2015年4月1日上海市地鐵一號線上25個(gè)站點(diǎn)(富錦路,友誼西路,寶安公路,共富新村,呼蘭路,通河新村,共康路,彭浦新村,汶水路,上海馬戲城,延長路中山北路,上?;疖囌?,漢中路,新闡路,人民廣場,黃陂南路,陜西南路,常熟路,衡山路徐家匯,上海體育館,漕寶路,上海南站,錦江樂園)的閘機(jī)進(jìn)出口的行人流數(shù)據(jù),每五分鐘作為一個(gè)計(jì)數(shù)節(jié)點(diǎn)。為保證各維度數(shù)據(jù)量一致,且通過檢驗(yàn)發(fā)現(xiàn)每天6:30之前和22:00之后地鐵閘機(jī)進(jìn)出口的人數(shù)較少,不影響之后的檢驗(yàn)結(jié)果,故將行人流數(shù)據(jù)的時(shí)間段截取為6:20-22:00,從而得到一個(gè)25維長度為188的數(shù)據(jù)集。通過運(yùn)用主成分分析和線性投影將多維的數(shù)據(jù)降成1維數(shù)據(jù)后,再利用WBS2.SDLL對地鐵閘機(jī)口進(jìn)站和出站的人數(shù)進(jìn)行變點(diǎn)檢驗(yàn),得到行人進(jìn)出的高峰期,通過算法檢驗(yàn)得到高峰時(shí)間段的行人流均值如表4所示。
表4 多個(gè)閘機(jī)進(jìn)出口高峰時(shí)間段
由表4可看出在時(shí)間段7:00-9:00和17:30-18:00內(nèi),進(jìn)站閘機(jī)口人流量達(dá)到最高峰,在時(shí)間段 8:05-9:30和 17:40-18:30內(nèi),出站閘機(jī)口人流量達(dá)到最高峰。結(jié)合實(shí)際與人們通勤時(shí)間相符合。同時(shí),為驗(yàn)證上述結(jié)果的可信度,接下來將對單個(gè)站點(diǎn)閘機(jī)進(jìn)出口行人流進(jìn)行變點(diǎn)檢驗(yàn),選取一號線中的一個(gè)換乘站:漢中路站(可換乘12和13線)作為表4中變點(diǎn)檢驗(yàn)結(jié)果的對比驗(yàn)證對象,結(jié)果如表5所示。
表5 漢中路站閘機(jī)進(jìn)出口高峰時(shí)間段
由表5可看出在時(shí)間段8:00-8:30和17:20-18:00內(nèi),漢口站進(jìn)站閘機(jī)口人流量達(dá)到最高峰,在時(shí)間段8:30-9:00內(nèi),該出站閘機(jī)口人流量達(dá)到最高峰。上述時(shí)間段均被涵蓋于多個(gè)站點(diǎn)檢驗(yàn)得出的高峰時(shí)間段內(nèi),說明通過結(jié)合主成分分析和線性投影來將數(shù)據(jù)進(jìn)行降維,再對降維后的數(shù)據(jù)進(jìn)行變點(diǎn)檢驗(yàn)的方法對地鐵行人流數(shù)據(jù)是有效的,檢驗(yàn)結(jié)果也是可信的。
通過對上海市一號線25個(gè)地鐵站的閘機(jī)進(jìn)出口行人數(shù)量的變點(diǎn)檢驗(yàn)分析研究,根據(jù)算法可以計(jì)算出行人流發(fā)生突變的次數(shù)和突變發(fā)生的具體時(shí)間,從而給交通管理部門對特定時(shí)間段的進(jìn)出閘機(jī)口數(shù)量的合理開放提供相關(guān)數(shù)據(jù)支持。
首先對大量階梯狀的模擬數(shù)據(jù)利用不同的的多變點(diǎn)檢驗(yàn)算法進(jìn)行變點(diǎn)估計(jì),得到表現(xiàn)較好的WBS2.95。再利用主成分分析和線性投影將一個(gè)25維的上海市一號線地鐵站閘機(jī)進(jìn)口的行人流數(shù)據(jù)有效地降成一維數(shù)據(jù),再利用WBS2.95進(jìn)行變點(diǎn)檢驗(yàn),結(jié)果表明7:00-8:30和17:20-19:00為進(jìn)站高峰期,其中 7:46-8:20人流量最大,8:05-9:00和 17:45-20:45為出站高峰期,其中18:15-18:55人流量最大。
將WBS2.SDLL方法運(yùn)用到行人交通樞紐中,同時(shí)考慮多個(gè)站點(diǎn)行人流的變點(diǎn)問題,為深入研究多個(gè)站點(diǎn)或多條線上的行人流突變問題提供一種新方法,但是WBS2.SDLL只適用于一維數(shù)據(jù),之后還將嘗試研究一種不需要事先進(jìn)行數(shù)據(jù)降維就能夠直接用于檢驗(yàn)多維數(shù)據(jù)的方法,為多維甚至高維數(shù)據(jù)的變點(diǎn)檢驗(yàn)問題提供參考。