韓冰凌，孫佳楠

（北京林業(yè)大學(xué) 理學(xué)院，北京 100083）

0 引言

變點(diǎn)檢測(cè)涉及的基礎(chǔ)學(xué)科涵蓋了數(shù)理統(tǒng)計(jì)、應(yīng)用數(shù)學(xué)、計(jì)算機(jī)科學(xué)等，并在金融學(xué)、經(jīng)濟(jì)學(xué)、氣象學(xué)、環(huán)境學(xué)等多個(gè)學(xué)科中廣泛應(yīng)用。例如，在金融學(xué)研究中，宿成建和陳潔[1]應(yīng)用變點(diǎn)模型研究了滬深股股市波動(dòng)性突變行為，并分析了1992—2002年上證和深證綜合指數(shù)的方差變點(diǎn)，對(duì)這些變點(diǎn)的經(jīng)濟(jì)意義進(jìn)行了解釋。在自然環(huán)境研究中，涂新軍和陳曉宏[2]基于變點(diǎn)原理，運(yùn)用似然比方法研究了存在變點(diǎn)的河川徑流量序列，并給出了一系列的結(jié)論。

注意到對(duì)泊松分布序列和指數(shù)分布序列的變點(diǎn)檢測(cè)研究不多見，但其具有較強(qiáng)的實(shí)際應(yīng)用價(jià)值。一些稀有事件如地震、煤礦災(zāi)難等的發(fā)生近似服從泊松分布，總結(jié)這些稀有事件的發(fā)生規(guī)律及發(fā)展過(guò)程中的規(guī)律突變，對(duì)于防范自然災(zāi)害等有重大意義，如對(duì)礦難發(fā)生次數(shù)的變點(diǎn)分析[3]、稀有事件變點(diǎn)問(wèn)題的分析[4]等。一方面，產(chǎn)品的壽命以及隨機(jī)服務(wù)系統(tǒng)的服務(wù)時(shí)長(zhǎng)等往往服從指數(shù)分布，檢測(cè)這些指數(shù)序列的變點(diǎn)，可以為提高生產(chǎn)質(zhì)量和改進(jìn)服務(wù)質(zhì)量提供科學(xué)依據(jù)，如黃志堅(jiān)和張志華[5]研究了可靠性數(shù)據(jù)在變點(diǎn)前后服從不同參數(shù)的指數(shù)分布產(chǎn)品的壽命分布，建立了產(chǎn)品故障分布的模型。

基于上述兩種分布的獨(dú)立序列數(shù)據(jù)，有必要通過(guò)模擬實(shí)驗(yàn)就不同變點(diǎn)檢測(cè)方法的檢測(cè)效果進(jìn)行比較研究，并給出能夠指導(dǎo)實(shí)際應(yīng)用的有價(jià)值的參考建議。R軟件中的Changepoint程序包[6]是近年開發(fā)的簡(jiǎn)單實(shí)用的變點(diǎn)檢測(cè)程序包，其中包含了經(jīng)典的僅一個(gè)變點(diǎn)（Atmost One Changepoint;AMOC）檢測(cè)法[6]和Binary Segmentation（BS）方法[7]，也包含了最近提出的Pruned Exact Linear Time（PELT）方法[8]。本文針對(duì)泊松和指數(shù)分布序列，使用該程序包下的上述三種變點(diǎn)檢測(cè)方法，分別對(duì)不同情境下的獨(dú)立泊松序列和指數(shù)序列進(jìn)行均值方差變點(diǎn)的檢測(cè)和比較。

1 三種均值方差變點(diǎn)檢測(cè)方法

1.1 變點(diǎn)問(wèn)題的提法

變點(diǎn)一般是指觀察序列中統(tǒng)計(jì)性質(zhì)發(fā)生變化的點(diǎn)的位置，統(tǒng)計(jì)學(xué)變點(diǎn)檢測(cè)問(wèn)題是對(duì)該位置的估計(jì)問(wèn)題。設(shè)一個(gè)按時(shí)間順序排列的觀察值序列記為 y1:n=(y1，...，yn)，若存在一個(gè)時(shí)間點(diǎn)τ∈{1，...，n-1}，使得這個(gè)時(shí)間點(diǎn)之前的序列{y1，...，yτ}和這個(gè)時(shí)間點(diǎn)之后的序列{yτ+1，...，yn}具有某方面不同的統(tǒng)計(jì)性質(zhì)，那么該時(shí)間點(diǎn)τ稱為一個(gè)變點(diǎn)。當(dāng)這兩個(gè)子序列的均值參數(shù)變化，τ稱為均值變點(diǎn)；當(dāng)這兩個(gè)子序列的均值和方差參數(shù)都變化，τ稱為均值方差變點(diǎn)。如果序列只存在一個(gè)變點(diǎn)，稱為單變點(diǎn)；如果變點(diǎn)數(shù)量為m，即存在不只一個(gè)變點(diǎn)，稱 τ1:m=(τ1，...，τm)為多變點(diǎn)[9]。

1.2 單變點(diǎn)問(wèn)題

單變點(diǎn)檢測(cè)可以轉(zhuǎn)化為假設(shè)檢驗(yàn)問(wèn)題，原假設(shè)是觀察值序列無(wú)變點(diǎn)，備擇假設(shè)為存在一個(gè)變點(diǎn)。該檢驗(yàn)問(wèn)題可通過(guò)似然比檢驗(yàn)實(shí)現(xiàn)，具體參見正態(tài)分布下均值單變點(diǎn)的檢測(cè)研究[10]和正態(tài)方差單變點(diǎn)的檢測(cè)研究[11]。

1.3 多變點(diǎn)問(wèn)題

對(duì)于多變點(diǎn)問(wèn)題，常將變點(diǎn)檢測(cè)問(wèn)題轉(zhuǎn)換為目標(biāo)函數(shù)的優(yōu)化問(wèn)題其中，C為損失函數(shù)，可以采用負(fù)對(duì)數(shù)似然函數(shù)。βf(m)為懲罰函數(shù)，可以采用 AIC[12]、BIC[13]的懲罰形式。

具體地，使用BS方法[7]優(yōu)化上述目標(biāo)函數(shù)的思想：第一步，在觀察值序列中只檢測(cè)一個(gè)變點(diǎn)的位置，如果序列中存在一個(gè)τ滿足 C(y1:τ)+C(y(τ+1):n)+β＜C(y1:n)，則認(rèn)為發(fā)現(xiàn)了一個(gè)變點(diǎn)；第二步，針對(duì)yτ分得的兩個(gè)子序列，分別進(jìn)行單變點(diǎn)檢測(cè)……直到每個(gè)子序列中不再檢測(cè)出變點(diǎn)。若第一步找不到單變點(diǎn)，則認(rèn)為此序列沒(méi)有變點(diǎn)。BS方法是將單變點(diǎn)檢測(cè)的思路應(yīng)用于多變點(diǎn)檢測(cè)問(wèn)題，方法中常取 f(m)=m。BS方法具有運(yùn)算效率高的優(yōu)點(diǎn)，但不能保證檢測(cè)出的變點(diǎn)是目標(biāo)函數(shù)優(yōu)化的全局最優(yōu)解。

若使用PELT方法[8]優(yōu)化上述目標(biāo)函數(shù)，則需以O(shè)ptimal Partitioning(OP)算法為基礎(chǔ)。OP算法的思想是采用遞歸的方式優(yōu)化目標(biāo)函數(shù)。記F(s)=min{F(t)+C(y(t+1):n) +β} ，其中，F(xiàn)(t)表示數(shù)據(jù) y1:t中函數(shù)最小值。OP算法沒(méi)有BS方法的求解效率高，于是PELT方法在OP算法的基礎(chǔ)上增加了一個(gè)剪枝[8]過(guò)程，通過(guò)剪枝操作來(lái)提高運(yùn)算效率，剪枝的本質(zhì)是去掉每次迭代過(guò)程中不能起到減小F(t)作用的τ。

2 模擬研究

針對(duì)獨(dú)立泊松分布和指數(shù)分布序列中的變點(diǎn)檢測(cè)問(wèn)題，分別應(yīng)用AMOC、PELT、BS方法進(jìn)行模擬實(shí)驗(yàn)并比較其效果，從而給出觀察值序列服從兩種不同分布下的方法選擇建議。

2.1 研究設(shè)計(jì)

模擬數(shù)據(jù)分別來(lái)自獨(dú)立泊松分布和獨(dú)立指數(shù)分布。每種分布下分別設(shè)計(jì)觀察值序列的樣本量為1000、1500；當(dāng)變點(diǎn)個(gè)數(shù)設(shè)計(jì)一個(gè)變點(diǎn)時(shí)，分布參數(shù)的變化范圍為由3變?yōu)?，或由2變?yōu)?.5；當(dāng)變點(diǎn)個(gè)數(shù)為兩個(gè)時(shí)，分布參數(shù)的變化范圍為由3變?yōu)?再變?yōu)?，或由2變?yōu)?.5再變?yōu)?；故共8種情境。每種模擬情境生成5組數(shù)據(jù)來(lái)進(jìn)行重復(fù)實(shí)驗(yàn)。目標(biāo)函數(shù)中分別采用AIC、BIC兩種信息準(zhǔn)則作懲罰項(xiàng)。研究中使用R Changepoint程序包的不同變點(diǎn)檢測(cè)方法來(lái)檢測(cè)泊松和指數(shù)分布中的變點(diǎn)。泊松分布與指數(shù)分布有一個(gè)共同的特點(diǎn)：均值參數(shù)和方差參數(shù)同時(shí)變化。因此使用均值方差變點(diǎn)命令cpt.meanvar進(jìn)行檢測(cè)。該程序包可以選擇檢測(cè)變點(diǎn)的懲罰項(xiàng)類型如AIC、BIC。

2.2 研究結(jié)果

模擬結(jié)果從以下角度分析：檢測(cè)的變點(diǎn)數(shù)、變點(diǎn)位置、輸出的負(fù)對(duì)數(shù)似然值的情況。影響結(jié)果的變量為觀察值序列的樣本量、變點(diǎn)個(gè)數(shù)、分布的參數(shù)、變點(diǎn)檢測(cè)方法的選取、懲罰函數(shù)類型。用N表示樣本個(gè)數(shù)，n表示變點(diǎn)個(gè)數(shù)，λ表示泊松或指數(shù)分布的參數(shù)。

2.2.1 泊松分布序列的變點(diǎn)檢測(cè)結(jié)果

在表1中，數(shù)字代表正確識(shí)別的變點(diǎn)數(shù)，“-”代表變點(diǎn)個(gè)數(shù)為2時(shí)不再使用AMOC方法。作為判定變點(diǎn)檢測(cè)效果的標(biāo)準(zhǔn)，此處著重考察每種方法得到的變點(diǎn)中，是否包含變點(diǎn)真值，即變點(diǎn)的準(zhǔn)確位置。具體地，從檢測(cè)到的變點(diǎn)中，首先選出距離真實(shí)變點(diǎn)最近的位置，再判定其是否距離真實(shí)變點(diǎn)在三個(gè)時(shí)間點(diǎn)以內(nèi)；若是，則視為檢測(cè)正確。由表1看出：當(dāng)變點(diǎn)數(shù)為一個(gè)時(shí)，使用三種方法正確檢測(cè)的變點(diǎn)數(shù)的均值相同，此時(shí)這三種方法沒(méi)有太大差異。當(dāng)變點(diǎn)數(shù)為兩個(gè)時(shí)，PELT方法正確檢測(cè)的變點(diǎn)數(shù)平均而言多于BS方法。比較兩種懲罰類型，BIC懲罰下正確檢測(cè)的變點(diǎn)數(shù)平均而言多于AIC懲罰。

表1 不同方法正確檢測(cè)泊松分布序列變點(diǎn)的個(gè)數(shù)

由表2看出：從檢測(cè)到的變點(diǎn)個(gè)數(shù)看，不管選取的樣本量及參數(shù)如何變化，當(dāng)真實(shí)情況存在一個(gè)變點(diǎn)時(shí)，顯然AMOC檢測(cè)的變點(diǎn)總數(shù)一定準(zhǔn)確，而其他兩種方法的變點(diǎn)數(shù)在使用AIC類型的懲罰項(xiàng)時(shí)均大于一個(gè)，使用BIC類型的懲罰項(xiàng)時(shí)表現(xiàn)較好。若樣本量不同，其他條件相同，使用AIC懲罰項(xiàng)的PELT方法找出的變點(diǎn)數(shù)會(huì)隨著樣本量的增大而增大，而其他情況檢測(cè)出的變點(diǎn)數(shù)量與樣本量變化無(wú)關(guān)。當(dāng)真實(shí)情況存在兩個(gè)變點(diǎn)時(shí)，顯然AMOC不再適用，使用AIC懲罰項(xiàng)的PELT方法找出的變點(diǎn)數(shù)會(huì)隨著樣本量的增大而增大?？傮w而言，針對(duì)懲罰類型選取的不同，BIC懲罰明顯優(yōu)于AIC懲罰下的變點(diǎn)識(shí)別效果。

表2 不同方法檢測(cè)的泊松分布序列的變點(diǎn)總數(shù)

再?gòu)呢?fù)對(duì)數(shù)似然值的大小來(lái)看（由于篇幅所限，不展示負(fù)對(duì)數(shù)似然值的表格）：若檢測(cè)變點(diǎn)的方法選取不同（不再考慮AMOC方法），使用AIC懲罰下的PELT方法有時(shí)會(huì)出現(xiàn)NAN的情況。針對(duì)不同的懲罰類型進(jìn)行比較，BIC懲罰下負(fù)對(duì)數(shù)似然值小于AIC懲罰下的負(fù)對(duì)數(shù)似然值。

2.2.2 指數(shù)分布序列的變點(diǎn)檢測(cè)結(jié)果

由表3，當(dāng)變點(diǎn)數(shù)為一個(gè)時(shí)，使用三種檢測(cè)方法正確檢測(cè)變點(diǎn)的平均數(shù)相差不大，此時(shí)這三種方法沒(méi)有太大差異。PELT方法相比于另兩種方法正確檢測(cè)的變點(diǎn)數(shù)略多。當(dāng)變點(diǎn)數(shù)為兩個(gè)時(shí)，比較BS方法和PELT方法，PELT方法正確檢測(cè)的變點(diǎn)數(shù)的均值與BS方法無(wú)明顯差異。BIC懲罰下正確檢測(cè)的變點(diǎn)數(shù)平均而言與AIC懲罰相似。

表3 不同方法正確檢測(cè)指數(shù)分布序列變點(diǎn)的個(gè)數(shù)

由表4看出，類似于泊松序列的研究結(jié)果，當(dāng)真實(shí)情況存在一個(gè)變點(diǎn)時(shí)，顯然AMOC檢測(cè)的變點(diǎn)總數(shù)一定準(zhǔn)確，而其他兩種方法的變點(diǎn)數(shù)在使用BIC懲罰項(xiàng)比AIC好。若樣本量不同，其他條件相同，使用AIC懲罰項(xiàng)的PELT方法找出的變點(diǎn)數(shù)會(huì)隨著樣本量的增大而增大，而其他情況檢測(cè)出的變點(diǎn)數(shù)量與樣本量變化無(wú)關(guān)。當(dāng)真實(shí)情況存在兩個(gè)變點(diǎn)時(shí)，使用AIC懲罰項(xiàng)的PELT方法明顯比BS方法差，但適用BIC懲罰時(shí)二者表現(xiàn)相似。

從負(fù)對(duì)數(shù)似然值的大小來(lái)看（由于篇幅所限，不展示負(fù)對(duì)數(shù)似然值的表格），N=1500時(shí)的負(fù)對(duì)數(shù)似然值要大于N=1000時(shí)的負(fù)對(duì)數(shù)似然值；若懲罰類型選取不同，其他變量均相同，使用PELT方法時(shí)，BIC懲罰下負(fù)對(duì)數(shù)似然值小于AIC懲罰下的負(fù)對(duì)數(shù)似然值；若使用BS方法，兩種懲罰方式下負(fù)對(duì)數(shù)似然值相同。

表4 不同方法檢測(cè)的指數(shù)分布序列的變點(diǎn)總數(shù)

3 實(shí)證

Carlin等（1992）[3]針對(duì)1851—1962年這 112年間英國(guó)每年發(fā)生煤礦災(zāi)難次數(shù)的數(shù)據(jù)，使用貝葉斯方法進(jìn)行變點(diǎn)檢測(cè)并找到一個(gè)變點(diǎn)k=41，其對(duì)應(yīng)年份為1891年；每年發(fā)生礦難的平均數(shù)由1891年之前的3.10下降到1891年之后的0.90。圖1為1852—1962年英國(guó)煤礦災(zāi)難每年的發(fā)生次數(shù)時(shí)序圖。

圖1 1851—1962年英國(guó)每年發(fā)生煤礦災(zāi)難次數(shù)的時(shí)序圖

根據(jù)該數(shù)據(jù)的產(chǎn)生背景，不妨假設(shè)序列中各隨機(jī)變量相互獨(dú)立并服從泊松分布。這里分別使用AMOC、BS和PELT方法分析該數(shù)據(jù)，觀察不同檢測(cè)變點(diǎn)方法及懲罰方式檢測(cè)變點(diǎn)的效果，并與Carlin等（1992）的研究結(jié)果進(jìn)行比較。由表5看出，AMOC方法可以較準(zhǔn)確地檢測(cè)到變點(diǎn)。PELT方法采用AIC懲罰時(shí)過(guò)于敏感，檢測(cè)出的變點(diǎn)數(shù)較多；采用BIC懲罰時(shí)可以減輕這種情況。BS方法表現(xiàn)較好，適用BIC懲罰時(shí)比AIC懲罰表現(xiàn)更好。

表5 三種方法對(duì)英國(guó)礦難數(shù)據(jù)的變點(diǎn)檢測(cè)結(jié)果

進(jìn)一步，如果對(duì)BS方法約束檢測(cè)到的變點(diǎn)個(gè)數(shù)為一個(gè)，BS方法檢測(cè)到的變點(diǎn)也是準(zhǔn)確的，結(jié)果見下頁(yè)表6。

表6 約束BS方法僅檢測(cè)一個(gè)變點(diǎn)的結(jié)果

4 結(jié)論

本文的研究得到以下結(jié)論：首先，對(duì)模擬研究，當(dāng)觀測(cè)值服從泊松分布且只有一個(gè)變點(diǎn)時(shí)，AMOC方法一定可以檢測(cè)到準(zhǔn)確的變點(diǎn)數(shù)量和位置，并且此時(shí)具有較小的負(fù)對(duì)數(shù)似然值，因此相對(duì)其他兩種方法較優(yōu)。若使用另外兩種方法，最好選擇BIC懲罰類型，不論從檢測(cè)出的變點(diǎn)個(gè)數(shù)準(zhǔn)確度或負(fù)對(duì)數(shù)似然值來(lái)看，使用BIC懲罰要優(yōu)于使用AIC懲罰。相對(duì)而言，BS方法要優(yōu)于PELT方法。當(dāng)泊松序列中存在兩個(gè)變點(diǎn)時(shí)，比較BS方法與PELT方法，看出兩種方法均在BIC懲罰下可以得到較準(zhǔn)確的變點(diǎn)數(shù)量和較小的負(fù)對(duì)數(shù)似然值。其次，指數(shù)分布的結(jié)果與泊松分布類似，當(dāng)變點(diǎn)個(gè)數(shù)為一個(gè)時(shí)，使用AMOC方法可以檢測(cè)到準(zhǔn)確的變點(diǎn)位置，并且此時(shí)具有較小的負(fù)對(duì)數(shù)似然值，相對(duì)其他兩種方法較優(yōu)。當(dāng)變點(diǎn)數(shù)為兩個(gè)時(shí)，使用BS方法與BIC懲罰結(jié)合使用、PELT方法與BIC懲罰結(jié)合使用得到的結(jié)果是類似的。再有，通過(guò)對(duì)實(shí)證研究中的變點(diǎn)檢測(cè)并與前人研究結(jié)果對(duì)比，發(fā)現(xiàn)AMOC、BS、PELT三種檢測(cè)方法的檢測(cè)效果優(yōu)劣與模擬結(jié)果類似。總之，泊松序列或指數(shù)序列存在一個(gè)變點(diǎn)時(shí)，使用均值方差同時(shí)變化的AMOC方法相比另外兩種方法更優(yōu)；對(duì)存在兩個(gè)變點(diǎn)的情況，BS或PELT結(jié)合BIC懲罰均較好，前者略優(yōu)于后者。本文的結(jié)果對(duì)于泊松和指數(shù)分布序列如何選擇三種方法來(lái)檢測(cè)變點(diǎn)具有較好的指導(dǎo)意義，未來(lái)研究還可探索對(duì)隨機(jī)變量序列服從其他分布類型時(shí)上述三種方法的變點(diǎn)檢測(cè)效果的比較。