采用密度比估計(jì)的多窗口變點(diǎn)檢測(cè)方法

張 曼，崔文泉

中國(guó)科學(xué)技術(shù)大學(xué) 管理學(xué)院 統(tǒng)計(jì)與金融系，合肥 230026

自新型冠狀病毒COVID-19爆發(fā)以來(lái)，全球疫情呈蔓延之勢(shì)，病毒具有較長(zhǎng)潛伏期和高度傳染性，在沒(méi)有政策干預(yù)的情況下傳播速度較快，會(huì)極大程度地干擾公共衛(wèi)生系統(tǒng)，給經(jīng)濟(jì)和社會(huì)帶來(lái)巨大的破壞。全球各國(guó)采取了多項(xiàng)緊急措施，如區(qū)域封鎖和大規(guī)模檢測(cè)等，一個(gè)自然的問(wèn)題是，這些干預(yù)措施是否以及在多大程度上有效地減緩疫情的傳播[1]，為此需要分析各個(gè)國(guó)家疫情的階段性特征。

COVID-19的傳播可能經(jīng)歷幾個(gè)階段，最初由于人群缺乏防范意識(shí)導(dǎo)致傳播速度較快，接著由于政府干預(yù)和公共衛(wèi)生措施傳播速度將逐漸放慢，這種階段性分析可以表述為疫情數(shù)據(jù)的變點(diǎn)檢測(cè)問(wèn)題。本文基于疫情背景，研究時(shí)間序列的變點(diǎn)檢測(cè)問(wèn)題。

時(shí)間序列是隨著時(shí)間變化對(duì)系統(tǒng)行為的描述，系統(tǒng)的行為會(huì)由于外部或內(nèi)部因素隨著時(shí)間而改變[2-3]。變點(diǎn)分析是統(tǒng)計(jì)學(xué)中的一個(gè)經(jīng)典分支，其基本定義是在一個(gè)序列或過(guò)程中，當(dāng)某個(gè)統(tǒng)計(jì)特性（分布類(lèi)型、分布參數(shù)）在某時(shí)間點(diǎn)受系統(tǒng)性因素而非偶然性因素影響發(fā)生變化，稱(chēng)該時(shí)間點(diǎn)為變點(diǎn)[4]。時(shí)間序列變點(diǎn)檢測(cè)問(wèn)題作為時(shí)間序列分析的一個(gè)重要研究方向，在經(jīng)濟(jì)、金融、醫(yī)學(xué)和氣象學(xué)等方面都有較廣的應(yīng)用[5]。

時(shí)間序列的變點(diǎn)檢測(cè)研究可追溯至Fox[6]在1972年提出的時(shí)間序列異常點(diǎn)檢測(cè)，此后時(shí)間序列變點(diǎn)檢測(cè)問(wèn)題逐漸受到更為廣泛的關(guān)注。當(dāng)前對(duì)于時(shí)間序列的變點(diǎn)檢測(cè)方法，可分為監(jiān)督式（Supervised）方法和無(wú)監(jiān)督（Unsupervised）方法[7]。

監(jiān)督式變點(diǎn)檢測(cè)方法是指收集足夠數(shù)量的不同狀態(tài)分布下的時(shí)間序列數(shù)據(jù)，結(jié)合這些樣本時(shí)間序列的類(lèi)別標(biāo)簽，對(duì)該序列集進(jìn)行監(jiān)督式訓(xùn)練，學(xué)習(xí)得到適合的模型，通過(guò)該模型對(duì)給定的時(shí)間序列進(jìn)行分類(lèi)識(shí)別，以確定其數(shù)據(jù)分布狀態(tài)[7]。監(jiān)督式方法包括基于支持向量機(jī)[8]、隨機(jī)森林[9-10]、神經(jīng)網(wǎng)絡(luò)[11]的檢測(cè)方法等。監(jiān)督學(xué)習(xí)需要事先大量地訓(xùn)練，在應(yīng)用中是一個(gè)挑戰(zhàn)[7]。

無(wú)監(jiān)督方法是指對(duì)無(wú)標(biāo)簽的時(shí)間序列直接進(jìn)行變點(diǎn)檢測(cè)，定量計(jì)算當(dāng)前時(shí)刻發(fā)生結(jié)構(gòu)變化的可能性，從而進(jìn)行變點(diǎn)分析[7]，與監(jiān)督變點(diǎn)檢測(cè)方法相比，無(wú)監(jiān)督的方法應(yīng)用更為廣泛。因此，本文對(duì)時(shí)間序列的無(wú)監(jiān)督變點(diǎn)檢測(cè)方法展開(kāi)研究。根據(jù)檢測(cè)延遲的不同，無(wú)監(jiān)督檢測(cè)方法可分為兩類(lèi)：實(shí)時(shí)檢測(cè)和回溯性檢測(cè)，其中回溯性檢測(cè)更穩(wěn)定，檢測(cè)效果一般會(huì)更準(zhǔn)確。

回溯性變點(diǎn)檢測(cè)方法需要待測(cè)時(shí)刻前后的數(shù)據(jù)，主要包括基于貝葉斯檢驗(yàn)[12]、似然比檢驗(yàn)[13]、累積和檢驗(yàn)[14]、密度估計(jì)的變點(diǎn)檢測(cè)等方法。

時(shí)間序列變點(diǎn)檢測(cè)問(wèn)題可描述如下[13]：

設(shè)y(t)∈?d是時(shí)間t時(shí)的d維時(shí)間序列，令：

是長(zhǎng)度為k的時(shí)間t處的時(shí)間序列的子序列。變點(diǎn)檢測(cè)問(wèn)題即檢測(cè)兩個(gè)相鄰時(shí)間段（稱(chēng)為參考區(qū)間reference和測(cè)試區(qū)間test）之間是否存在變點(diǎn)。設(shè)tref和ttest(tref＜ttest)分別是參考區(qū)間和測(cè)試區(qū)間起始時(shí)間點(diǎn)，ptest(Y)和pref(Y)分別是參考區(qū)間和測(cè)試區(qū)間樣本的概率密度函數(shù)。假設(shè)參考區(qū)間和測(cè)試區(qū)間的樣本量分別為nref和ntest，則ttest=tref+nref，為簡(jiǎn)潔，記Yref(i)=Y(tref+i-1)，Ytest(i)=Y(ttest+i-1)，則第i個(gè)參考樣本記為Yref(i)，第i個(gè)測(cè)試樣本記為Ytest(i)，于是變點(diǎn)檢測(cè)問(wèn)題可表述為如下假設(shè)檢驗(yàn)問(wèn)題：

Brodsky[15]提出了基于密度估計(jì)的變點(diǎn)檢測(cè)方法，在高維的情形下因“維度災(zāi)難”導(dǎo)致檢測(cè)精度大大降低。

為了克服這一困難，直接估計(jì)概率密度比值的密度比估計(jì)方法被廣泛應(yīng)用起來(lái)，密度比框架在不進(jìn)行密度估計(jì)的情況下，直接估計(jì)密度比，可以高效、統(tǒng)一地解決很多機(jī)器學(xué)習(xí)問(wèn)題。因此研究如何更好地利用密度比來(lái)對(duì)時(shí)間序列進(jìn)行變點(diǎn)檢測(cè)具有重要研究意義和實(shí)際應(yīng)用價(jià)值。已有的密度比估計(jì)包含如下方法[16]：核均值匹配[17]、Logistic回歸方法[18]、Kullback-Leibler重要性估計(jì)程序[19（]KLIEP）、無(wú)約束最小二乘重要性估計(jì)[20]（unconstrained least squares importance fitting，uLSIF），相對(duì)無(wú)約束最小二乘重要性估計(jì)[21（]relative uLSIF，RuLSIF）。

Kawahara等[13]和Liu等[22]提出了基于密度比的時(shí)間序列變點(diǎn)檢測(cè)方法，主要是通過(guò)比較相鄰時(shí)間窗數(shù)據(jù)的分布差異來(lái)識(shí)別變點(diǎn)。本文的目標(biāo)是繼續(xù)推進(jìn)基于密度比的變點(diǎn)檢測(cè)這一研究路線，同時(shí)考慮結(jié)合COVID-19全球蔓延的應(yīng)用背景。然而該類(lèi)方法存在局限性，Liu等[22]提出，基于密度比的變點(diǎn)檢測(cè)方法的檢測(cè)效果受到超參數(shù)窗寬的影響。Kawahara等[13]指出，窗寬越大，密度比的估計(jì)就越準(zhǔn)確，但在變點(diǎn)個(gè)數(shù)較多的情形下，可能漏檢相鄰較近的一些變點(diǎn)。

鑒于本文的研究背景是基于疫情數(shù)據(jù)，政府干預(yù)、公共衛(wèi)生措施等諸多宏觀、微觀影響因素會(huì)導(dǎo)致疫情數(shù)據(jù)發(fā)生結(jié)構(gòu)性變化，由于影響因素的多元性和復(fù)雜性，致使疫情數(shù)據(jù)變點(diǎn)出現(xiàn)的時(shí)間尺度不一致。然而基于密度比的變點(diǎn)檢測(cè)方法[22]中不同時(shí)間窗檢測(cè)效果不盡相同，考慮對(duì)時(shí)間窗窗寬進(jìn)行超參數(shù)選優(yōu)，超參數(shù)選優(yōu)往往需要足夠數(shù)量的驗(yàn)證集更新超參數(shù)，當(dāng)COVID-19病例數(shù)據(jù)較少時(shí)，需要進(jìn)一步探究基于密度比的時(shí)間序列變點(diǎn)檢測(cè)方法中時(shí)間窗窗寬的選擇問(wèn)題。

Messer等[23]提出多重過(guò)濾算法（multiple filter algorithm，MFA），該方法按照一定的算法匯集不同時(shí)間窗檢測(cè)所得的備選變點(diǎn)集。

為了更好地對(duì)疫情數(shù)據(jù)等時(shí)間序列數(shù)據(jù)進(jìn)行變點(diǎn)分析，本文在基于密度比的變點(diǎn)檢測(cè)方法的基礎(chǔ)上，利用和發(fā)展動(dòng)態(tài)多重過(guò)濾MFA算法，提出一種新方法：基于密度比的時(shí)間序列多窗口變點(diǎn)檢測(cè)方法（mDRCPD），可同時(shí)應(yīng)用多個(gè)窗口進(jìn)行檢測(cè)。引入多窗口使得該方法更加適用于類(lèi)似于疫情數(shù)據(jù)等多變點(diǎn)時(shí)間序列數(shù)據(jù)，這是已有的時(shí)間序列變點(diǎn)檢測(cè)中所沒(méi)有出現(xiàn)過(guò)的方法。

接著采用絕對(duì)誤差、精確度、召回率和F1得分等指標(biāo)來(lái)驗(yàn)證mDRCPD方法的有效性，模擬結(jié)果表明該方法在各項(xiàng)性能指標(biāo)下表現(xiàn)優(yōu)于單窗口方法。

最后將mDRCPD變點(diǎn)檢測(cè)方法拓展到COVID-19傳播進(jìn)程的分析中，與Wang等[24]提出的疫情預(yù)測(cè)的生存卷積模型組合，通過(guò)對(duì)傳播率的分段建模來(lái)刻畫(huà)疫情的階段性特點(diǎn)，評(píng)估國(guó)家政策在控制疫情傳播上的效果?，F(xiàn)有文獻(xiàn)中關(guān)于疫情數(shù)據(jù)變點(diǎn)分析的研究數(shù)量相對(duì)不多，該分析能夠?yàn)檎疀Q策及疫情分析提供有效依據(jù)，在疫情大背景下具有重要的現(xiàn)實(shí)意義。

1 方法及算法

mDRCPD方法包括固定窗口下尋找變點(diǎn)集和變點(diǎn)集的匯集，下面將分別介紹這兩個(gè)步驟及其中涉及的方法原理。

1.1 基于密度比的時(shí)間序列變點(diǎn)檢測(cè)方法

設(shè)y(t)∈?d是時(shí)間t時(shí)的d維時(shí)間序列樣本，令

是長(zhǎng)度為k的時(shí)間t處的時(shí)間序列的子序列。本文以Y(t)為樣本點(diǎn)，為了能夠利用密度比建模方法處理具有相依性結(jié)構(gòu)的時(shí)間序列數(shù)據(jù)，類(lèi)似于文獻(xiàn)[13，22]中的處理方式，需要進(jìn)行“工作獨(dú)立假定”，視窗寬為n的樣本是獨(dú)立同分布的[25]，對(duì)原始時(shí)間序列的相依性是由作為樣本點(diǎn)的子序列Y(t)進(jìn)行刻畫(huà)的。不妨設(shè)新樣本容量為T(mén)。

稱(chēng)Ytref和Yttest分別為t時(shí)刻窗口為nref和ntest的相鄰時(shí)間窗，在工作獨(dú)立假定下，Ytref和Yttest分別為樣本容量為nref和ntest的簡(jiǎn)單隨機(jī)樣本。

基于密度比的時(shí)間序列變點(diǎn)檢測(cè)方法步驟如下：

（1）考慮兩個(gè)相鄰時(shí)間窗Ytref和Yttest，然后計(jì)算Ytref和Yttest間的某種差異性度量作為點(diǎn)t的檢測(cè)得分[26]，記為score(t)，檢測(cè)得分越高，差異性越大，點(diǎn)t越可能是變點(diǎn)。

（2）計(jì)算每個(gè)時(shí)間t的檢測(cè)分?jǐn)?shù)：

檢測(cè)得分隨t變化的曲線的超過(guò)預(yù)先設(shè)定的某個(gè)閾值的峰值點(diǎn)[22]處即為估計(jì)的變點(diǎn)。

這里給出t點(diǎn)處閾值的選取方法：仿照Sugiyama等[27]提出的最小二乘兩樣本置換檢驗(yàn)，針對(duì)樣本集Ytref和Yttest，對(duì)Ytref∪Yttest進(jìn)行隨機(jī)置換，置換后的前nref個(gè)樣本記為，剩余ntest個(gè)樣本稱(chēng)為，再對(duì)和計(jì)算檢測(cè)分?jǐn)?shù)，如此重復(fù)M0次，得到置換后樣本的M0組檢測(cè)分?jǐn)?shù)的上α1分位數(shù)即為t處閾值。

由于nref和ntest表示時(shí)間窗的長(zhǎng)度，故稱(chēng)之為窗寬，不妨稱(chēng)nref為左窗寬，ntest為右窗寬，由于這里的左右窗寬是固定不變的值，為與本文提出的方法區(qū)分，稱(chēng)這里的方法為基于密度比的單窗口變點(diǎn)檢測(cè)方法（single window density-ratio change point detection，sDRCPD）。

方法使用對(duì)稱(chēng)散度來(lái)度量分布的差異性：

其中，Ptref和Pttest分別表示Ytref和Yttest所對(duì)應(yīng)總體的分布，PEα(P||P′)表示α相對(duì)Pearson（PE）散度[21]：

其中，p(Y)和p′(Y)分別是P和P′的概率密度函數(shù)，0≤α＜1，p′α(Y)=αp(Y)+(1-α)p′(Y)稱(chēng)為α混合密度，rα(Y)=p(Y)/(αp(Y)+(1-α)p′(Y))稱(chēng)為α相對(duì)密度比。

本文根據(jù)基于RuLSIF的密度比估計(jì)方法對(duì)PE散度進(jìn)行估計(jì)[21]，Sugiyama等[16]給出以下核模型對(duì)密度比p(Y)/p′α(Y)進(jìn)行建模：

其中，θ:=(θ1,θ2,…,θn)T是從樣本中學(xué)習(xí)的參數(shù)，K(Y,Y′)是基函數(shù)，本文使用高斯核：

其中，σ(＞0)表示核寬，在平方損失J(Y)最優(yōu)的準(zhǔn)則下對(duì)密度比進(jìn)行估計(jì)：

第一項(xiàng)與θ無(wú)關(guān)，故只需考慮最后三項(xiàng)。利用式（6），并對(duì)目標(biāo)函數(shù)（8）進(jìn)行經(jīng)驗(yàn)化處理，給出了RuLSIF優(yōu)化問(wèn)題：

其中，λθTθ/2為懲罰項(xiàng)，λ(≥0)表示正則化參數(shù)，是n×n矩陣，第(?,?′)元素為：

式（9）的解為：

其中，In表示n維單位矩陣，密度比估計(jì)為：

由于缺乏非負(fù)限制，估計(jì)的密度比可能是負(fù)的，為了解決這個(gè)問(wèn)題，負(fù)的輸出可以這樣處理：

利用α相對(duì)密度比的估計(jì)量，α相對(duì)PE散度的可估計(jì)為[21]：

1.2 mDRCPD方法

在基于密度比的單窗口變點(diǎn)檢測(cè)方法中，不同的窗寬檢測(cè)效果不盡相同，大窗寬時(shí)密度比估計(jì)更準(zhǔn)確，小窗寬則更適用于檢測(cè)變化較快的點(diǎn)[13]，但同一段時(shí)間序列數(shù)據(jù)包含的變點(diǎn)時(shí)間尺度多樣，僅僅用單窗口滑動(dòng)來(lái)檢測(cè)變點(diǎn)有一定的局限性。為了進(jìn)一步驗(yàn)證引入多窗口方法的必要性，本節(jié)實(shí)驗(yàn)分析了改進(jìn)后的多窗口方法和未改進(jìn)的單窗口方法在檢測(cè)不同時(shí)間尺度變點(diǎn)上的區(qū)別。

實(shí)驗(yàn)仿照Takeuchi and Yamanish（i2006）[28]中一維自回歸模型下生成1 000個(gè)樣本（即t=1,2,…,1 000）：

其中，εt是均值為μ和標(biāo)準(zhǔn)差σ=1.5的高斯噪聲，初始值y(1)=y(2)=0，通過(guò)設(shè)置μ插入變點(diǎn)。變點(diǎn)分別為C1=400，C2=420，C3=440，C4=460，C5=480，C6=500，C7=600，C8=700，C9=800，C10=900，其中Ci表示第i個(gè)變點(diǎn)，人為設(shè)置前5個(gè)變點(diǎn)變化頻率快，變化程度大，后5個(gè)變點(diǎn)變化頻率慢，變化程度小。

實(shí)驗(yàn)結(jié)果如表1，表1顯示不同窗寬下基于密度比的變點(diǎn)檢測(cè)方法分別檢測(cè)到了不同類(lèi)型變點(diǎn)，小窗寬檢測(cè)方法對(duì)前5個(gè)變點(diǎn)較為敏感，大窗寬則對(duì)后5個(gè)變點(diǎn)檢測(cè)效果好。實(shí)驗(yàn)結(jié)果說(shuō)明了“小窗寬不易漏檢變化較快的點(diǎn)，大窗寬則對(duì)點(diǎn)識(shí)別更準(zhǔn)確”，主要原因在于大窗寬時(shí)間窗包含的數(shù)據(jù)信息多，密度比估計(jì)更為準(zhǔn)確，故變點(diǎn)識(shí)別更為準(zhǔn)確，但在變點(diǎn)距離較近時(shí)，容易受到同一個(gè)時(shí)間窗內(nèi)存在不止一個(gè)變點(diǎn)的影響從而降低檢測(cè)效果。

表1 不同時(shí)間尺度變點(diǎn)檢測(cè)結(jié)果Table 1 Detection results of change points on different time scales

在本次實(shí)驗(yàn)中，各窗寬檢測(cè)效果差異較大，因此有必要引入多窗口方法，從而改善單窗口方法受窗寬影響，檢測(cè)效果不足的問(wèn)題。為了進(jìn)一步提高檢測(cè)效果，本文引入Messer等提出的多重過(guò)濾算法MFA[23]，即同時(shí)應(yīng)用多個(gè)窗口進(jìn)行檢測(cè)，并將結(jié)果匯集。將引入MFA算法后的基于密度比的時(shí)間序列變點(diǎn)檢測(cè)方法稱(chēng)為多窗口方法，即mDRCPD，mDRCPD算法步驟如下：

（1）首先在每個(gè)固定的窗寬nref和ntest下尋找各自變點(diǎn)集，從觀測(cè)時(shí)間點(diǎn)t∈{nref,nref+1,…,T-ntest}開(kāi)始，得到時(shí)間點(diǎn)t的相對(duì)PE散度，即檢測(cè)分?jǐn)?shù)。

（2）若檢測(cè)分?jǐn)?shù)的最大值（設(shè)為M）大于閾值，則認(rèn)為最大值所在的t是估計(jì)的變點(diǎn)。

（3）通過(guò)去除上一個(gè)估計(jì)變點(diǎn)的(t-nref,t+ntest)鄰域來(lái)更新候選點(diǎn)并繼續(xù)尋找下一個(gè)最大值點(diǎn)。

（4）這種算法重復(fù)直到M低于閾值停止。

（6）對(duì)于任意點(diǎn)v∈P2，僅當(dāng)該點(diǎn)不屬于集合任何現(xiàn)有變點(diǎn)的左鄰域，右鄰域時(shí)，點(diǎn)v才被添加到集合P中。

（7）該過(guò)程一直向前推進(jìn)，直到左右窗寬為和停止。

算法1和算法2具體描述mDRCPD方法的代碼實(shí)現(xiàn)過(guò)程：

算法1基于密度比的單窗口變點(diǎn)檢測(cè)（sDRCPD）

2 模擬分析

為了驗(yàn)證提出的mDRCPD方法相較于單窗口變點(diǎn)檢測(cè)方法的優(yōu)越性，本文模擬分析了兩種方法在四類(lèi)多變點(diǎn)時(shí)間序列數(shù)據(jù)集上的表現(xiàn)。對(duì)于所有實(shí)驗(yàn)，RuLSIF方法中的α=0.1，k=10[22]，閾值選取中M0=100，α1=0.01，窗寬取nref=ntest=n，niref=ntiest=ni，i=1,2,…,|H|，σ和λ的參數(shù)調(diào)優(yōu)采用5折交叉驗(yàn)證，兩組候選參數(shù)如下：

（1）σ=0.6dmed,0.8dmed,dmed,1.2dmed,1.4dmed

（2）λ=10-3,10-2,10-1,100,101

其中，dmed表示樣本之間的距離中位數(shù)。

2.1 模擬數(shù)據(jù)集

模擬使用的時(shí)間序列數(shù)據(jù)集具體介紹如下，其中包含人工插入的變點(diǎn)：

數(shù)據(jù)集1（均值變化）：采用Takeuchi and Yamanishi（2006）[28]中一維自回歸模型下生成1 000個(gè)樣本（即t=1,2,…,1 000）：

其中，εt是均值為μ和標(biāo)準(zhǔn)差σ=1.5的高斯噪聲，初始值y(1)=y(2)=0，通過(guò)設(shè)置μ在每100個(gè)時(shí)間點(diǎn)插入一個(gè)變點(diǎn)：

其中，M是使得100(M-1)+1≤t≤100M的自然數(shù)。

數(shù)據(jù)集2（方差變化）使用與數(shù)據(jù)集1相同的自回歸模型生成1 000個(gè)樣本，但通過(guò)設(shè)置σ在每100個(gè)時(shí)間點(diǎn)中插入一個(gè)變點(diǎn)：

數(shù)據(jù)集3（協(xié)方差變化）從原點(diǎn)為中心的正態(tài)分布中生成1 000個(gè)二維樣本，通過(guò)設(shè)置時(shí)間t處的協(xié)方差矩陣Σ在每100個(gè)時(shí)間點(diǎn)處插入一個(gè)變點(diǎn)：

數(shù)據(jù)集4（頻率變化）通過(guò)以下模型生成1 000個(gè)樣本：

其中，εt是以原點(diǎn)為中心的高斯噪聲，標(biāo)準(zhǔn)差為0.8，通過(guò)改變時(shí)間t的頻率ω在每100個(gè)時(shí)間點(diǎn)插入一個(gè)變點(diǎn)：

上述數(shù)據(jù)集中后一個(gè)變點(diǎn)均比前一個(gè)變點(diǎn)更重要。

模擬實(shí)驗(yàn)中統(tǒng)計(jì)了100次蒙特卡洛模擬中的平均指標(biāo)數(shù)據(jù)，采用的評(píng)價(jià)指標(biāo)為：絕對(duì)誤差、精確度、召回率、TP（True Positive）和F1得分。其中絕對(duì)誤差、精確度、召回率、TP和F1得分都取為100次模擬的均值，2.2節(jié)將介紹所使用的評(píng)價(jià)指標(biāo)。

2.2 評(píng)價(jià)指標(biāo)

真實(shí)的變點(diǎn)集用T*={t1*,t2*,…,tK**}表示，估計(jì)的變點(diǎn)集用表示[30]。絕對(duì)誤差記為ΔAE，表示預(yù)測(cè)的變點(diǎn)數(shù)和真實(shí)的變點(diǎn)數(shù)|T*|(=K*)之間的差異：

精確度表示估計(jì)的變點(diǎn)中是真實(shí)變點(diǎn)的比例。召回率表示真正的變點(diǎn)中被估計(jì)出的比例。類(lèi)似Messer等[23]的做法，對(duì)于單窗口檢測(cè)方法來(lái)說(shuō)，一個(gè)估計(jì)的變點(diǎn)若它的n2鄰域包含了真實(shí)的變點(diǎn)，就認(rèn)為它是正確檢測(cè)的；而對(duì)于多窗口檢測(cè)方法，一個(gè)估計(jì)的變點(diǎn)若它的ni2鄰域包含了真實(shí)的變點(diǎn)，方認(rèn)為它是正確檢測(cè)的，這里的ni指的是該點(diǎn)被添加到最終變點(diǎn)集前所屬的窗寬ni。TP（true positive）表示被正確檢測(cè)的變點(diǎn)數(shù)，即：

精確度（Precision）和召回率（Recall）定義為：

F1得分是精度和召回率的調(diào)和平均值，記為ΔF1：

F1得分最優(yōu)為1，最差為0。

2.3 模擬結(jié)果

表2～5分別展示了四種數(shù)據(jù)集在單窗口檢測(cè)方法和本文提出的多窗口檢測(cè)方法下的評(píng)價(jià)結(jié)果，其中n1=100，n2=90，n3=80，n4=70，n5=60，n6=50。

在表2～表5中，絕對(duì)誤差指標(biāo)越小表現(xiàn)越優(yōu)，而TP、Presicion、Recall和F1得分四列指標(biāo)均是越大表現(xiàn)越優(yōu)。

表2 數(shù)據(jù)集1評(píng)價(jià)指標(biāo)對(duì)比結(jié)果Table 2 Comparison of evaluation index on dataset 1

表5 數(shù)據(jù)集4評(píng)價(jià)指標(biāo)對(duì)比結(jié)果Table 5 Comparison of evaluation index on dataset 4

表3 數(shù)據(jù)集2評(píng)價(jià)指標(biāo)對(duì)比結(jié)果Table 3 Comparison of evaluation index on dataset 2

表4 數(shù)據(jù)集3評(píng)價(jià)指標(biāo)對(duì)比結(jié)果Table 4 Comparison of evaluation index on dataset 3

針對(duì)數(shù)據(jù)集1，mDRCPD方法相對(duì)于表現(xiàn)較好的單窗口方法（加粗顯示，下同），絕對(duì)誤差降低了1個(gè)單位，TP提高了1.4個(gè)單位，精確度提高了3.76%，召回率提高了29.23%，F(xiàn)1得分提高了19.91%。

針對(duì)數(shù)據(jù)集2，mDRCPD方法相對(duì)于表現(xiàn)較好的單窗口方法，絕對(duì)誤差降低了1.46個(gè)單位，TP提高了1.5個(gè)單位，召回率提高了17.08%，F(xiàn)1得分提高了7.08%。

針對(duì)數(shù)據(jù)集3，mDRCPD方法相對(duì)于表現(xiàn)較好的單窗口方法，絕對(duì)誤差降低了4.12個(gè)單位，TP提高了2.94個(gè)單位，召回率提高了66.03%，F(xiàn)1得分提高了29.68%。

針對(duì)數(shù)據(jù)集4，mDRCPD方法相對(duì)于表現(xiàn)較好的單窗口方法，絕對(duì)誤差降低了0.5個(gè)單位，TP提高了0.3個(gè)單位，召回率提高了4.76%，F(xiàn)1得分提高了2.19%。模擬結(jié)果表明了提出的mDRCPD方法在檢測(cè)多變點(diǎn)時(shí)間序列時(shí)的優(yōu)越性。

3 實(shí)證分析

3.1 基于嬰兒心率數(shù)據(jù)集

為進(jìn)一步論證mDRCPD多窗口變點(diǎn)檢測(cè)方法在真實(shí)數(shù)據(jù)集上仍然取得較好效果，利用R包wavethresh中提供的醫(yī)學(xué)時(shí)間序列數(shù)據(jù)集BabyECG和BabySS進(jìn)行實(shí)證分析。BabyECG是某嬰兒心率的記錄（以每分鐘心跳次數(shù)為單位），每16秒采樣一次，記錄時(shí)間為21：17：59到06：27：18，包含2 048個(gè)觀測(cè)值。BabySS是嬰兒睡眠狀態(tài)的記錄，由經(jīng)過(guò)培訓(xùn)的專(zhuān)家監(jiān)測(cè)EEG（大腦）和EOG（眼動(dòng)）確定，等級(jí)為1到4。睡眠狀態(tài)代碼是1=安靜睡眠，2=安靜和活躍睡眠之間，3=活躍睡眠，4=清醒。

本文應(yīng)用變點(diǎn)檢測(cè)模型來(lái)檢測(cè)該嬰兒睡眠狀態(tài)的變化時(shí)刻，而B(niǎo)abySS數(shù)據(jù)中睡眠狀態(tài)切換的時(shí)刻記為真實(shí)變點(diǎn)。表6展示了嬰兒心率記錄數(shù)據(jù)集在單窗口檢測(cè)方法和提出的多窗口檢測(cè)方法下的評(píng)價(jià)結(jié)果，其中n1=100，n2=90，n3=80，n4=70，n5=60，n6=50。

表6 嬰兒心率記錄數(shù)據(jù)集評(píng)價(jià)結(jié)果Table 6 Comparison of evaluation index in infant’s heart rate record data

從表6中可以看出，mDRCPD方法相對(duì)于表現(xiàn)較好的單窗口方法（加粗顯示），絕對(duì)誤差降低了3個(gè)單位，TP提高了2個(gè)單位，召回率提高了13.34%，F(xiàn)1得分提高了6.53%，該實(shí)證分析結(jié)果也表明mDRCPD方法相較于單窗口方法在實(shí)際時(shí)間序列數(shù)據(jù)變點(diǎn)檢測(cè)中的優(yōu)越性。

3.2 基于COVID-19數(shù)據(jù)集

3.2.1 COVID-19數(shù)據(jù)變點(diǎn)檢測(cè)及政策效果評(píng)估

全球各國(guó)已經(jīng)采取了前所未有的措施來(lái)緩解COVID-19的蔓延，一個(gè)自然的問(wèn)題是各國(guó)的干預(yù)措施是否有效遏制了疫情傳播的速度。

Wang等[24]提出了生存-卷積模型來(lái)評(píng)估國(guó)家政策對(duì)降低疫情傳播速度的效果和預(yù)測(cè)COVID-19病例數(shù)據(jù)。該方法對(duì)傳播率函數(shù)分段線性建模，當(dāng)大規(guī)模公共衛(wèi)生干預(yù)在特定日期實(shí)施時(shí)，引入額外的線性函數(shù)和新的斜率參數(shù)，則干預(yù)前后斜率的變化反映了政策對(duì)降低疾病傳播速度的影響。接著假設(shè)未來(lái)一段時(shí)期傳播率斜率保持不變，將歷史數(shù)據(jù)代入生存卷積模型中，預(yù)測(cè)未來(lái)數(shù)據(jù)。本節(jié)將提出的mDPCPD變點(diǎn)檢測(cè)方法與生存卷積模型組合，應(yīng)用到COVID-19傳播進(jìn)程的分析中。

本小節(jié)實(shí)證分析的數(shù)據(jù)為COVID-19傳播過(guò)程中每日新增病例（含境外輸入）時(shí)間序列數(shù)據(jù)，來(lái)源是https：//ourworldindata.org/coronavirus-source-data，分析的國(guó)家包括：美國(guó)、巴西、俄羅斯、英國(guó)、德國(guó)、墨西哥、印度、韓國(guó)，數(shù)據(jù)對(duì)應(yīng)時(shí)間從每日新增確診病例數(shù)超過(guò)3的日期開(kāi)始，截止日期為2021年4月6日，各國(guó)實(shí)證分析中使用的時(shí)間序列數(shù)據(jù)樣本量在表7中給出。在分析過(guò)程中，選取時(shí)間序列中2021年3月18日至2021年4月6日的20個(gè)樣本作為測(cè)試集，其余數(shù)據(jù)作為訓(xùn)練集。

利用mDRCPD變點(diǎn)檢測(cè)模型估計(jì)變點(diǎn)，以變點(diǎn)為節(jié)點(diǎn)，對(duì)傳播率函數(shù)分段線性建模，變點(diǎn)前后傳播率函數(shù)斜率的變化反映了變點(diǎn)處政策對(duì)降低疫情傳播速度的影響。

表7總結(jié)了8個(gè)代表性國(guó)家的基于mDRCPD方法的檢測(cè)結(jié)果，其中包括檢測(cè)時(shí)間序列數(shù)據(jù)的開(kāi)始日期，時(shí)間序列數(shù)據(jù)長(zhǎng)度(m)，估計(jì)的變點(diǎn)個(gè)數(shù)（NO.CP），前3個(gè)變點(diǎn)對(duì)應(yīng)的日期（1stCP，2ndCP，3rdCP），ai表示第i個(gè)變點(diǎn)到第i+1個(gè)變點(diǎn)之間傳播率函數(shù)a(t)的斜率。

表7 基于mDRCPD的疫情數(shù)據(jù)變點(diǎn)細(xì)節(jié)展示Table 7 Detail of change point detection on epidemic data based on mDRCPD

將部分變點(diǎn)與國(guó)家政策和估計(jì)得到的傳播率聯(lián)系，對(duì)政策效果進(jìn)行評(píng)估。對(duì)于美國(guó)，變點(diǎn)出現(xiàn)在2020年3月24日，距離2020年3月13日“國(guó)家宣布全國(guó)進(jìn)入緊急狀態(tài)”有11天，差不多是一個(gè)正常的病毒潛伏期，且從3月20日起，“美國(guó)開(kāi)始禁止在過(guò)去14天內(nèi)曾去往28個(gè)歐洲國(guó)家旅行的外國(guó)公民入境”，變點(diǎn)的發(fā)生與此也對(duì)應(yīng)。

對(duì)于巴西，變點(diǎn)是3月27日與3月30日，“感染者逼近一萬(wàn)例時(shí)，巴西政府陸續(xù)關(guān)閉邊境，并呼吁民眾戴口罩出行”事件對(duì)應(yīng)。

對(duì)于俄羅斯，變點(diǎn)出現(xiàn)在2020年4月2日與“3月18日至5月1日期間，限制外國(guó)人和無(wú)國(guó)籍人員進(jìn)入俄羅斯；3月27日起，完全停止與所有國(guó)家的包機(jī)和常規(guī)航班”的政策對(duì)應(yīng)。

對(duì)于德國(guó)，變點(diǎn)是2020年3月16日，當(dāng)時(shí)“3月16日德國(guó)封鎖了邊境，17日，德國(guó)措施開(kāi)始加碼，聯(lián)邦政府提議關(guān)閉大量商店”，變點(diǎn)日期與之呼應(yīng)。

對(duì)于墨西哥，變點(diǎn)日期為2020年4月14日與“4月4日，墨西哥政府宣布即日起至4月30日全國(guó)進(jìn)入公共衛(wèi)生緊急狀態(tài)，在此期間公共和私營(yíng)部門(mén)的所有非必要活動(dòng)全部暫?！笔录鄬?duì)應(yīng)。

對(duì)于印度，變點(diǎn)發(fā)生時(shí)間為2020年4月3日與“3月25日，印度開(kāi)始進(jìn)行為期21天的全國(guó)封鎖；印度鐵路停駛?cè)珖?guó)客運(yùn)列車(chē)至4月14日”事件對(duì)應(yīng)。

對(duì)于韓國(guó)，變點(diǎn)出現(xiàn)在2020年2月29日與“3月3日，韓國(guó)總統(tǒng)文在寅在國(guó)務(wù)會(huì)議發(fā)言中強(qiáng)調(diào)，所有政府機(jī)構(gòu)要啟動(dòng)24小時(shí)應(yīng)急體制”事件對(duì)應(yīng)。

上面所分析的變點(diǎn)對(duì)應(yīng)的內(nèi)在表現(xiàn)均為變點(diǎn)后一個(gè)階段傳播率的較快下降，說(shuō)明了政策的效果和公共衛(wèi)生干預(yù)對(duì)降低疫情傳播速度的重要性。

3.2.2 基于mDRCPD的COVID-19兩階段預(yù)測(cè)

對(duì)COVID-19的每日新增確診病例數(shù)的預(yù)測(cè)對(duì)于公共衛(wèi)生決策至關(guān)重要，可以幫助政府以數(shù)據(jù)作為決策依據(jù)。提出一種兩階段方法來(lái)進(jìn)行每日新增病例數(shù)短期預(yù)測(cè)，第一階段利用mDRCPD變點(diǎn)檢測(cè)模型估計(jì)變點(diǎn)，第二階段將變點(diǎn)作為傳播率的分段節(jié)點(diǎn)，假設(shè)傳播率短期的線性趨勢(shì)保持不變（即斜率保持不變）進(jìn)行預(yù)測(cè)，預(yù)測(cè)結(jié)果部分展示如圖1，其中虛線表示訓(xùn)練數(shù)據(jù)的估計(jì)變點(diǎn)。

對(duì)于預(yù)測(cè)結(jié)果的評(píng)價(jià)，仿照Wang等[24]的方式構(gòu)造預(yù)測(cè)結(jié)果的置信區(qū)間。設(shè)θ表示模型中的參數(shù)，將每日新增病例數(shù)劃分為訓(xùn)練集和測(cè)試集，訓(xùn)練集表示為：

其中，t1和t2分別表示訓(xùn)練集的開(kāi)始時(shí)間和結(jié)束時(shí)間。模型為，其中ε(t)表示殘差項(xiàng)，Y(t;θ)表示生存卷積模型預(yù)測(cè)的每日新增病例數(shù)。估計(jì)置信區(qū)間時(shí)，假設(shè)殘差值是可交換的，因此可以使用置換方法。首先擬合模型，得到殘差，再置換殘差，并利用置換的殘差重構(gòu)觀察到的每日新增病例數(shù)，重新擬合模型得到新的預(yù)測(cè)結(jié)果，如此重復(fù)這個(gè)過(guò)程100次，以產(chǎn)生一組新的預(yù)測(cè)值Y(t;θ)，最后估計(jì)100組預(yù)測(cè)值的經(jīng)驗(yàn)分位數(shù)獲得了95%的置信區(qū)間。

圖1顯示，真實(shí)值大部分都包含在預(yù)測(cè)值的95%置信區(qū)間內(nèi)，說(shuō)明預(yù)測(cè)是有效的。

圖1 新增確診病例觀測(cè)值和預(yù)測(cè)值Fig.1 Observed and predicted values of newly confirmed cases

4 算法分析

針對(duì)本文提出的mDRCPD檢測(cè)方法，本文給出了算法的性能分析和與其他變點(diǎn)檢測(cè)方法的比較分析，據(jù)此，來(lái)驗(yàn)證本文算法的有效性和可行性。

4.1 性能分析

本節(jié)對(duì)改進(jìn)后方法的性能分析分為兩個(gè)階段，一是算法完成前的理論分析，即復(fù)雜度分析，可以從時(shí)間和空間兩個(gè)維度去衡量，二是算法完成后的實(shí)驗(yàn)分析，即靈敏度分析。

4.1.1 算法復(fù)雜度分析

任何模型都不是兩全其美的，在基于密度比的變點(diǎn)檢測(cè)中引入多窗口算法MFA的同時(shí)，可能帶來(lái)了算法復(fù)雜度的增加，這里進(jìn)行復(fù)雜度分析。假設(shè)所有其他因素，如處理器的速度等是恒定的，對(duì)算法的實(shí)現(xiàn)沒(méi)有影響。本文的算法復(fù)雜度從以下兩個(gè)方面來(lái)說(shuō)明。

時(shí)間復(fù)雜度是指執(zhí)行這個(gè)算法所需要的計(jì)算工作量，其復(fù)雜度反映了程序執(zhí)行時(shí)間隨輸入規(guī)模增長(zhǎng)而增長(zhǎng)的量級(jí)，在很大程度上能很好地反映出算法的優(yōu)劣與否。在已有的單窗口方法中，見(jiàn)1.2節(jié)中算法1，只需要在固定窗口下遍歷集合C={nref,nref+1,…,T-ntest}中所有的點(diǎn)，計(jì)算每個(gè)時(shí)間點(diǎn)的基于密度比的檢測(cè)得分，時(shí)間復(fù)雜度記為O(T)，其中T為樣本容量。

在提出的多窗口mDRCPD方法中，見(jiàn)1.2節(jié)算法2，需要在窗寬序列集n={n1,n2,…n|H|}中每個(gè)窗寬下遍歷集合C，時(shí)間復(fù)雜度變成O(|H|×T)，然而在本文的模擬和實(shí)際應(yīng)用中，|H|往往取為3～6的較小常數(shù)。故相較于單窗口方法，多窗口方法并未給時(shí)間復(fù)雜度帶來(lái)較大增長(zhǎng)，仍可記為O(T)。

空間復(fù)雜度是指一個(gè)算法在運(yùn)行過(guò)程中臨時(shí)占用存儲(chǔ)空間大小的量度，已有的單窗口方法和改進(jìn)后的多窗口方法空間復(fù)雜度均為O(d×k×T)，由于d和k為常數(shù)，故記為O(T)。

在模擬分析和實(shí)際數(shù)據(jù)分析中已論證多窗口方法在檢測(cè)效果上的顯著提升，故權(quán)衡的來(lái)說(shuō)，在實(shí)際應(yīng)用中，往往愿意犧牲計(jì)算時(shí)間上的稍微代價(jià)，換取更好的檢測(cè)效果。

4.1.2 靈敏度分析

針對(duì)本文提出的mDRCPD檢測(cè)方法，這里給出了基于窗寬序列集n={n1,n2,…,nl}和調(diào)整參數(shù)k的靈敏度分析。采取不同的窗寬序列集n={n1,n2,…,nl}和調(diào)整參數(shù)k進(jìn)行模擬實(shí)驗(yàn)，表8給出了mDRCPD方法和單窗口方法在不同窗寬序列集和不同k下的F1得分，其中黑體加粗部分表示基于mDRCPD方法的F1得分，未加粗部分表示各單窗口下的最優(yōu)F1得分。

表8 靈敏度分析結(jié)果Table 8 Sensitivity analysis results

表8顯示，隨著窗寬序列集和調(diào)整參數(shù)k的變化，本文提出的多窗口方法性能始終優(yōu)于單窗口方法，該結(jié)果并沒(méi)有因?yàn)楦鲄?shù)的選取不同而產(chǎn)生較大差異，可見(jiàn)提出的模型是穩(wěn)定有效的。

4.2 與現(xiàn)有方法的對(duì)比分析

本文提出的方法是引入多重過(guò)濾算法MFA，同時(shí)對(duì)時(shí)間序列多窗口分割，在文獻(xiàn)[30]中總結(jié)了變點(diǎn)檢測(cè)的主要分割算法，典型的有二進(jìn)制分割（binary segmentatio，Binseg）和自下而上分割（bottom-up segmentation，BottomUp）。Binseg是一種順序方法，首先在完整的輸入序列中檢測(cè)到一個(gè)變點(diǎn)，然后據(jù)此變點(diǎn)將序列拆分，并對(duì)拆分后的子序列重復(fù)上述操作。BottomUp方法是先將序列沿著規(guī)則網(wǎng)格分割成許多子序列，后根據(jù)相似程度不斷合并連續(xù)的段，即先從許多變點(diǎn)開(kāi)始，然后不斷刪除不太重要的變點(diǎn)。

為了進(jìn)一步驗(yàn)證引入多重過(guò)濾算法對(duì)時(shí)間序列同時(shí)多窗口分割的必要性，本節(jié)進(jìn)行了與其他兩種主要分割方法的對(duì)比分析實(shí)驗(yàn)，實(shí)驗(yàn)中mDRCPD方法各參數(shù)取值同第2章模擬分析，實(shí)驗(yàn)結(jié)果如表9。

表9 與現(xiàn)有方法比較結(jié)果Table 9 Comparison of results with existing methods

從表9可知，在數(shù)據(jù)集1～4中，本文提出的多窗口方法的表現(xiàn)均優(yōu)于另外兩種分割過(guò)濾算法。

5 結(jié)束語(yǔ)

本文在基于密度比的變點(diǎn)檢測(cè)方法基礎(chǔ)上進(jìn)行了改進(jìn)，引入了MFA算法，提出了多窗口mDRCPD方法，并進(jìn)行了模擬實(shí)驗(yàn)和實(shí)證分析，結(jié)果均表明該方法在各項(xiàng)性能指標(biāo)下表現(xiàn)優(yōu)于單窗口方法，因此，本文的研究具有應(yīng)用和實(shí)用價(jià)值。

其次將變點(diǎn)檢測(cè)mDRCPD算法和疫情預(yù)測(cè)的生存卷積模型組合，對(duì)疫情數(shù)據(jù)進(jìn)行統(tǒng)計(jì)學(xué)變點(diǎn)檢測(cè)并與國(guó)家政策聯(lián)系，從而進(jìn)行政策效果的評(píng)估以及短期預(yù)測(cè)，最后利用置換法構(gòu)造預(yù)測(cè)值的95%置信區(qū)間，說(shuō)明了預(yù)測(cè)結(jié)果的有效性。

本文提出的多窗口mDRCPD方法適用于生物醫(yī)學(xué)、經(jīng)濟(jì)學(xué)、金融學(xué)等領(lǐng)域時(shí)間序列數(shù)據(jù)的變點(diǎn)分析問(wèn)題，有效改善了單窗口方法受時(shí)間窗窗寬影響，檢測(cè)變點(diǎn)時(shí)間尺度單一，檢測(cè)效果不足的問(wèn)題。