銀炳皓，施三支，齊德全

（長(zhǎng)春理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，長(zhǎng)春 130022）

變點(diǎn)估計(jì)和檢驗(yàn)問題是統(tǒng)計(jì)學(xué)中的一個(gè)經(jīng)典問題，用于識(shí)別時(shí)間序列中的突變，發(fā)生這種突然變化的時(shí)間點(diǎn)稱為“變點(diǎn)”。通過估計(jì)變點(diǎn)的位置，可以將原始數(shù)據(jù)集劃分為較短的段，然后使用平穩(wěn)時(shí)間序列的方法進(jìn)行分析，解決各個(gè)領(lǐng)域相關(guān)的問題。變點(diǎn)問題在金融行業(yè)[1]、互聯(lián)網(wǎng)安全[2]、生物醫(yī)學(xué)[3]等領(lǐng)域都有著廣泛的應(yīng)用，例如，變點(diǎn)估計(jì)可用于發(fā)出有關(guān)異常金融事件的警報(bào)，檢測(cè)網(wǎng)絡(luò)上的分布式拒絕服務(wù)攻擊，精確定位人腦中腦波活動(dòng)的開始等。隨著時(shí)代的發(fā)展和科技的進(jìn)步，現(xiàn)實(shí)生活中的數(shù)據(jù)多以高維的形式出現(xiàn)，高維數(shù)據(jù)流導(dǎo)致了數(shù)據(jù)量的激增，這就要求變點(diǎn)分析方法能有所進(jìn)步，并且在運(yùn)算時(shí)間和方法準(zhǔn)確性上達(dá)到一個(gè)良好的平衡效果。

經(jīng)典的變點(diǎn)檢驗(yàn)關(guān)注的是單變量時(shí)間序列，例如2012年Killick[4]提出的PELT方法是通過找到變點(diǎn)可能的數(shù)量和位置的成本函數(shù)最小值，來推斷變點(diǎn)的最佳數(shù)量和位置。2014年Frick和Munk等人[5]提出的SMUCE算法針對(duì)指數(shù)族回歸中的變點(diǎn)問題，引入了一種新的估計(jì)量；2014年Fryzlewicz[6]提出的 WBS 算法和 2018 年 Baranowski等人[7]提出的最窄超閾值法（NOT）提高了二元分割算法的性能。然而，經(jīng)典的單變量變點(diǎn)檢驗(yàn)方法通常不直接適用于現(xiàn)代應(yīng)用中經(jīng)常遇到的高維數(shù)據(jù)集。因此，學(xué)者們提出了幾種方法來估計(jì)和檢驗(yàn)高維數(shù)據(jù)中的變點(diǎn)，并且為了解決維數(shù)災(zāi)難的問題，現(xiàn)有的幾種高維變點(diǎn)檢驗(yàn)方法假設(shè)了變化信號(hào)具有某種形式的稀疏性，稀疏性假設(shè)大大降低了原始高維問題的復(fù)雜性。例如，在 2015年 Jirak[8]使用 CUSUM 統(tǒng)計(jì)量對(duì)不同子部分進(jìn)行L∞聚合，對(duì)于稀疏變點(diǎn)的效果很好；2015 年 Cho 和 Fryzlowicz[9]的稀疏二元分割算法，也考慮了稀疏性；2018年Wang和Samworth[10]基于投影的方法中假設(shè)了一個(gè)變點(diǎn)前后的均值差異只在小部分維數(shù)坐標(biāo)中是非零的；2019 年 Enikeeva和 Harchaoui[11]提出了基于具有稀疏性的掃描統(tǒng)計(jì)量的方法。關(guān)于高維數(shù)據(jù)的變點(diǎn)檢驗(yàn)問題的工具仍在不斷創(chuàng)新和發(fā)展，但大多不具有普適性，且運(yùn)算時(shí)間較長(zhǎng)。

本文通過研究具有稀疏變點(diǎn)的高維數(shù)據(jù)的特點(diǎn)，提出了一種基于稀疏最優(yōu)投影的多變點(diǎn)檢驗(yàn)方法，SOP-SeedBS方法能夠檢驗(yàn)出數(shù)據(jù)的維度p在小于或大于等于數(shù)據(jù)量n的情況下變點(diǎn)的位置及個(gè)數(shù)，并研究了變點(diǎn)間隔大小對(duì)該檢驗(yàn)方法穩(wěn)定性的影響。通過仿真模擬，在多種情況下與Wang和Samworth的Inspect方法（2018）進(jìn)行了比較，結(jié)果表明，SOP-SeedBS方法具有較高的準(zhǔn)確度，并在計(jì)算時(shí)間方面具有明顯的優(yōu)勢(shì)。

1 SOP-SeedBS高維稀疏多變點(diǎn)檢驗(yàn)

1.1 稀疏變點(diǎn)模型

設(shè)X1,…,Xn是n個(gè)相互獨(dú)立的p維隨機(jī)向量，Xt～Np(μt,σ2Ip)，1 ≤t≤n，σ2是方差，記

其中，p可能與序列的長(zhǎng)度n相等，甚至大于序列的長(zhǎng)度n。

假設(shè)存在υ個(gè)變點(diǎn)，變點(diǎn)檢驗(yàn)的問題能被表示為經(jīng)典的假設(shè)檢驗(yàn)問題

假設(shè)備擇假設(shè)成立，數(shù)據(jù)可劃分成υ+1段的分段常數(shù)均值結(jié)構(gòu)，換句話說，如果有υ個(gè)變點(diǎn) 1≤z1＜z2＜…＜zυ，即?0 ≤i≤υ。記z0=0,zυ+1=n， 記θ(i)=μ(i)-μ(i-1),i=1,…,υ。

本文假設(shè)均值的變化是稀疏的，即一個(gè)變點(diǎn)前后的均值差異只在小部分維度（p）中是非零的。記‖θ(i)‖0表示向量θ(i)中非零元素的個(gè)數(shù)，則‖θ(i)‖0≤k，k∈{1,…,p}且k?p。本文只考慮參數(shù)μ(i)的檢驗(yàn)問題，其他參數(shù)皆視為討厭參數(shù)。假設(shè)變點(diǎn)位置z1，z2，…，zυ滿足：

其中，τ是變點(diǎn)的最小間隔，0＜τ＜1。均值變化的大小滿足：

其中，?是均值變化大小的閾值。

1.2 稀疏最優(yōu)投影（SOP）

Wang等人[10]（2018）通過對(duì)高維數(shù)據(jù)進(jìn)行高維CUSUM變換進(jìn)而推導(dǎo)出數(shù)據(jù)矩陣的最優(yōu)投影方向，并解決了計(jì)算第一個(gè)左奇異向量引出的凸優(yōu)化問題。因此對(duì)于高維稀疏變點(diǎn)數(shù)據(jù)，本文根據(jù)類似的思想，并基于奇異值分解和凸松弛優(yōu)化方法將高維數(shù)據(jù)投影成一維數(shù)據(jù)，由于本文研究的高維數(shù)據(jù)中均值變化是稀疏的，從而稱該方法為稀疏最優(yōu)投影（Sparse Optimal Projection）。

對(duì)于任意的a∈Rp，滿足 ‖a‖2=1，其中

選擇a=θ/‖θ‖2來最大限度地?cái)U(kuò)大兩個(gè)數(shù)據(jù)片段之間的均值差異的大小。然而θ在實(shí)踐中通常是未知的，本文的想法是尋找一個(gè)最優(yōu)的投影方向v=θ/‖θ‖2，文獻(xiàn)[10]啟發(fā)了本文對(duì)數(shù)據(jù)矩陣X進(jìn)行分解，再利用奇異值分解得到一個(gè)合適的投影方向后，將數(shù)據(jù)投影到一維空間。

假設(shè)均值變化只出現(xiàn)在數(shù)據(jù)矩陣X的k個(gè)列中，且在沒有發(fā)生均值變化的p-k個(gè)列中的元素均值為μt=0，t∈{1,…,n}，k∈{1,…,p} 。對(duì)數(shù)據(jù)矩陣X進(jìn)行分解，得到X=μ+W，其中μ=(μ1,…,μn)T∈Rn×p是均值矩陣，W=(W1,…,Wn)T∈Rn×p，W1,…,Wn是獨(dú)立的Np(0,σ2Ip)的隨機(jī)噪聲。當(dāng)σ已知，零假設(shè)成立時(shí)，μ為恒定的常數(shù)向量，備擇假設(shè)下只存在一個(gè)變點(diǎn)z時(shí)，矩陣X的元素均值可以計(jì)算出：

其中，Xr,j為矩陣X中第r行，第j列的元素；t∈{1,…,n} ；j∈ {1,…,p}。記前z行向量元素相同，后n-z行向量元素相同。

這里St,j為第j列中變點(diǎn)z前后分段數(shù)據(jù)之和，t∈{1,…,n}，j∈{1,…,p} ，即：

則有：

文獻(xiàn)[10]中，計(jì)算矩陣XT的第一個(gè)左奇異向量是一個(gè)組合優(yōu)化問題。公式（10）是一個(gè)計(jì)算復(fù)雜度很高的凸優(yōu)化問題，解決這個(gè)問題的一種方法是使用優(yōu)化問題中的凸松弛得到M?來代替X，具體過程如下：

懲罰函數(shù)可以將有約束的目標(biāo)函數(shù)轉(zhuǎn)化為無約束的目標(biāo)函數(shù)，即當(dāng)?shù)玫降腗不滿足約束條件時(shí)，估計(jì)的?可能會(huì)遠(yuǎn)離最優(yōu)解。通過添加一個(gè)懲罰項(xiàng) -λ‖M‖1，優(yōu)化問題變?yōu)橐韵滦问剑?/p>

其中，λ為軟閾值。則有：

記是的第一個(gè)左奇異向量，即：

將高維數(shù)據(jù)X沿著最優(yōu)方向投影，得到一維數(shù)據(jù)，記為Y。

如果數(shù)據(jù)存在變點(diǎn)，投影后的一維數(shù)據(jù)呈現(xiàn)出了分段常數(shù)的特征，適用于此種數(shù)據(jù)的幾種經(jīng)典的單變量多變點(diǎn)檢驗(yàn)方法，包括2012年Killick[4]給 出 的 PELT 方 法 、2014 年 Fryzlewicz[6]提出的 WBS方法、2018年 Baranowski等人[7]給出的最窄超閾值方法（NOT）以及2020年Kovács等人[12]提出的SeedBS方法。PELT方法是通過找到變點(diǎn)可能的數(shù)量和位置的成本函數(shù)最小值來推斷變點(diǎn)的最佳數(shù)量和位置，計(jì)算復(fù)雜度較高。WBS方法和NOT方法都是基于二元分割思想的改進(jìn)算法，NOT方法中隨機(jī)抽取間隔的想法與WBS方法類似，但NOT方法更關(guān)注數(shù)據(jù)中可能存在變點(diǎn)的最小間隔。WBS和NOT方法中用于搜索變點(diǎn)的間隔都是通過隨機(jī)抽樣構(gòu)造的，在數(shù)據(jù)量較大的情況下，需要構(gòu)造的隨機(jī)間隔數(shù)量會(huì)很多，且會(huì)出現(xiàn)許多長(zhǎng)度與原始序列相同或者重復(fù)出現(xiàn)的間隔，這增加了算法復(fù)雜度和運(yùn)算時(shí)間，特別是高維稀疏變點(diǎn)數(shù)據(jù)中變點(diǎn)較少的情況下。2020 年 Kovács等人[12]提出的種子二元分割（SeedBS）是一種解決大規(guī)模數(shù)據(jù)變點(diǎn)檢驗(yàn)問題的方法，該方法能夠有效地降低計(jì)算復(fù)雜度，提高計(jì)算效率。因此，本文選擇SeedBS算法來對(duì)投影后的一維數(shù)據(jù)進(jìn)行多變點(diǎn)檢驗(yàn)。

1.3 種子二元分割（SeedBS）

種子二元分割方法是通過構(gòu)造一個(gè)具有確定性，且可以預(yù)先計(jì)算的搜索區(qū)間集來搜索變點(diǎn)。區(qū)間集的構(gòu)造獨(dú)立于潛在的變點(diǎn)數(shù)量，并能達(dá)到搜索間隔的總長(zhǎng)度盡可能小的效果。

其中，μt是t時(shí)刻的信號(hào)；σ是恒定標(biāo)準(zhǔn)差；εt是正態(tài)噪聲，記有υ個(gè)變點(diǎn)1≤z1＜z2＜ ...＜zυ，則μt可以被分為υ+1段，其中第i段中。

記a∈(1/2,1]為一個(gè)給定的衰減參數(shù)，對(duì)于任意的q：1≤q≤log1/a(n)，將第q層定義為由nq個(gè)區(qū)間組成，初始長(zhǎng)度為lq且進(jìn)行確定性位移sq的集合Lq：

例如當(dāng)a=1/2時(shí)，得到一個(gè)重疊的二元分隔區(qū)間。來自前一層的每個(gè)間隔都被分成左、右和重疊的中間間隔，每個(gè)間隔的大小都是前一層的一半。當(dāng)間隔的長(zhǎng)度在層上迅速衰減時(shí)，隨著間隔長(zhǎng)度變小，間隔的數(shù)量迅速增加，還保證了在同一層內(nèi)的連續(xù)間隔之間有足夠大的重疊。生成的間隔如圖1所示[12]。

圖1 a=1/2時(shí)生成搜索間隔的可視化圖

用于種子二元分割的檢驗(yàn)統(tǒng)計(jì)量是定義一段內(nèi)分割點(diǎn)在s的CUSUM統(tǒng)計(jì)量：

其中，0≤l＜r≤n和h=r-l為整數(shù)，達(dá)到CUSUM統(tǒng)計(jì)量絕對(duì)值最大值的位置s記為最好的分割點(diǎn)。因此，計(jì)算m個(gè)種子間隔的檢驗(yàn)統(tǒng)計(jì)量，得到一個(gè)候選變點(diǎn)集，記為，如果候選變點(diǎn)集合中對(duì)應(yīng)間隔的檢驗(yàn)統(tǒng)計(jì)量超過了閾值ζn，就認(rèn)為得到一個(gè)最終的變點(diǎn)估計(jì)集合O(0,n]，并將聲明為一個(gè)變點(diǎn)：

其中：

閾值ζn可以通過強(qiáng)化施瓦茨信息準(zhǔn)則計(jì)算：

其中，α≥ 1；是模型（2）分段參數(shù)μ(i)的最大似然估計(jì)；nυ表示需要估計(jì)參數(shù)的總個(gè)數(shù)，包括變點(diǎn)位置。另外，當(dāng)α=1時(shí)相當(dāng)于施瓦茨信息準(zhǔn)則。

在具有固定的方差σ2和分段常數(shù)均值的單變量高斯信號(hào)模型中，也可以通過貝葉斯信息準(zhǔn)則（BIC）來計(jì)算閾值：

其中，σ2是方差；n為樣本量；T={z1,…,zυ} ? {1,…,n} ；|T|為集合T的基數(shù)。通過BIC準(zhǔn)則計(jì)算出的閾值篩選出的變點(diǎn)個(gè)數(shù)要多于強(qiáng)化施瓦茨信息準(zhǔn)則篩選出的變點(diǎn)個(gè)數(shù)。

1.4 SOP-SeedBS方法

本文結(jié)合了類似 Wang等人[10]（2018）的投影思想和 Kovács等人[12]（2020）提出的種子二元分割（SeedBS）方法，進(jìn)行高維數(shù)據(jù)的稀疏變點(diǎn)檢驗(yàn)，先將高維稀疏變點(diǎn)數(shù)據(jù)沿著最優(yōu)的投影方向投影到一維，再利用SeedBS方法對(duì)一維數(shù)據(jù)進(jìn)行多變點(diǎn)檢驗(yàn)，SOP-SeedBS方法具體步驟如下：

輸入：數(shù)據(jù)矩陣X=(X1,...,Xn)T∈Rn×p，軟閾值λ，衰減參數(shù)足夠大的變點(diǎn)上界K，max閾值ζn。

（1）對(duì)數(shù)據(jù)矩陣X進(jìn)行分解并推導(dǎo)出最優(yōu)投影方向。

（3）對(duì)進(jìn)行SVD分解，得到第一個(gè)左奇異向量，計(jì)算投影數(shù)據(jù)

（4）根據(jù)衰減參數(shù)a構(gòu)造包含m個(gè)種子區(qū)間的區(qū)間集。

（5）計(jì)算m個(gè)種子間隔的檢驗(yàn)統(tǒng)計(jì)量，并根據(jù)閾值ζn確定變點(diǎn)位置：

2 數(shù)值模擬

構(gòu)造由模型1生成的具有均值變點(diǎn)的高維數(shù)據(jù)時(shí)，設(shè)置參數(shù)σ2=1，數(shù)據(jù)量n分別取200，500，1 000，2 000，維數(shù)p分別為200，500，1 200，將高維數(shù)據(jù)噪聲表示為ε，分別取0.1，0.5，0.7，1，變點(diǎn)位置滿足z/n={0 .25,0.5,0.75}，定義?(i)為第i個(gè)變點(diǎn)的均值跳躍大小，=(?(1),?(2),?(3))=(0.4,0.8,1.2)，?(0)=0，，其中k控制著發(fā)生變化維度數(shù)量的稀疏程度，取值分別為10，40。根據(jù)文獻(xiàn)[10]的建議，。根據(jù)文獻(xiàn)[12]的建議，設(shè)置衰減參數(shù)a= 2。另外，下面三種情形分別考慮了均值結(jié)構(gòu)中發(fā)生變化的維數(shù)坐標(biāo)(p)上的重疊問題，高維數(shù)據(jù)中均值結(jié)構(gòu)發(fā)生變化的三種不同重疊情形說明如下：

情形1（S1）：在沒有重疊的情況下，對(duì)于相應(yīng)的變點(diǎn)，發(fā)生變化的k個(gè)維數(shù)坐標(biāo)是來自完全不相交的坐標(biāo)集；

情形3（S3）：在完全重疊的情況下，則對(duì)于相應(yīng)的變點(diǎn)，變化的維數(shù)發(fā)生在相同的k個(gè)坐標(biāo)上。

本文先給出了SOP-SeedBS方法在情形3下，參數(shù)設(shè)置為σ2=1，n=2 000，p=200，k=40，??=(?(1),?(2),?(3))=(0.4,0.8,1.2)，變點(diǎn)位置z/n={0.25,0.5,0.75}時(shí)，不同噪聲下的變點(diǎn)估計(jì)結(jié)果可視化，用ε表示高維數(shù)據(jù)的噪聲，不同ε取值對(duì)于投影后一維數(shù)據(jù)的變點(diǎn)檢驗(yàn)結(jié)果如圖2所示，圖2的結(jié)果顯示了當(dāng)ε的取值越小時(shí)，噪聲越小，分段特征越明顯。另外，當(dāng)均值跳躍越大，分段特征也越明顯。后續(xù)的仿真實(shí)驗(yàn)中，噪聲ε皆設(shè)置為1。

圖2 不同 ε下的投影后一維數(shù)據(jù)的變點(diǎn)檢驗(yàn)結(jié)果

表1給出了σ2=1時(shí)，SOP-SeedBS方法和Inspect方法的比較結(jié)果，兩種方法在四組不同的參數(shù)設(shè)置M1、M2、M3、M4下，且每組參數(shù)設(shè)置都具有三種不同的情形，通過100次運(yùn)行算法（兩種方法通過相同參數(shù)設(shè)置下生成的高維數(shù)據(jù)進(jìn)行比較）檢測(cè)到的變點(diǎn)數(shù)量的頻數(shù)和平均Hausdorff距離dH分別來度量估計(jì)變點(diǎn)的數(shù)量和估計(jì)變點(diǎn)位置的準(zhǔn)確度。M1的參數(shù)設(shè)置為n=2 000，p=200，k=10，=(0.8,1.6,2.4)；M2的參數(shù)設(shè)置為n=2 000，p=200，k=40，=(0.4,0.8,1.2)；M3的參數(shù)設(shè)置為n=500，p=500，k=40，=(0.4,0.8,1.2)；M4的參數(shù)設(shè)置為n=200，p=500，k=40，=(0.4,0.8,1.2)。表中Z1，Z2，Z3，Z4，Z5表示估計(jì)到的變點(diǎn)個(gè)數(shù)，取值分別為1，2，3，4，5。

表1 兩種方法在不同參數(shù)設(shè)置下進(jìn)行100次模擬的估計(jì)變點(diǎn)頻數(shù)和平均Hausdorff距離的比較

從表1可以看出，在σ2=1，n=2 000，p=200，k=10的參數(shù)設(shè)置下，SOP-SeedBS方法在三種不同均值變化重疊情況下的結(jié)果都優(yōu)于Inspect方法，進(jìn)行100次模擬的平均Hausdorff距離分別為 6.87，2.32，0.47，而 Inspect方法分別為12.72，26.01，13.91；第二組實(shí)驗(yàn)將第一組實(shí)驗(yàn)中的發(fā)生變化的維數(shù)k=10替換為k=40，第一個(gè)變點(diǎn)的均值跳躍大小由0.8變?yōu)?.4，SOPSeedBS方法在三種情形下進(jìn)行100次模擬得到的平均Hausdorff距離分別為4.48，2.42，0.34。雖然兩組實(shí)驗(yàn)在變點(diǎn)頻數(shù)上同樣具有良好的檢驗(yàn)結(jié)果，但第二組結(jié)果與第一組相比，在均值跳躍變小的情況下，第二組結(jié)果中情形1和情形3的平均Hausdorff距離要比第一組更小，說明k的大小也會(huì)對(duì)實(shí)驗(yàn)結(jié)果產(chǎn)生影響。

從SOP-SeedBS方法的第三組和第四組實(shí)驗(yàn)結(jié)果可以看出，在其他參數(shù)固定的情況下，n越大，對(duì)比實(shí)驗(yàn)結(jié)果顯示n較大時(shí)檢驗(yàn)結(jié)果較好，說明由于樣本量的增加，適當(dāng)增加種子區(qū)間的個(gè)數(shù)有利于進(jìn)行變點(diǎn)檢驗(yàn)。本文還通過模擬實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)其他參數(shù)固定時(shí)，維度p變大，對(duì)比實(shí)驗(yàn)的結(jié)果表現(xiàn)出檢驗(yàn)效果相近，因此維度的增加對(duì)于變點(diǎn)檢驗(yàn)效果影響不大。

圖3給出了SOP-SeedBS方法與Inspect方法在其他參數(shù)不變的情況下，不同的n和p設(shè)置下算法平均運(yùn)算時(shí)間的比較，可以看出SOP-SeedBS方法在算法運(yùn)算速度上具有較大的優(yōu)勢(shì)。圖3（a）給出了當(dāng)p=200，n分別為 500、1 000、2 000時(shí)，兩種方法100次模擬的平均運(yùn)算時(shí)間；圖3（b）給出了當(dāng)p=500，n分別為 200、500、1 000時(shí)，兩種方法進(jìn)行100次模擬的平均運(yùn)算時(shí)間。

圖3 SOP-SeedBS方法與Inspect方法計(jì)算時(shí)間比較

圖3中的結(jié)果顯示，當(dāng)p=200，n=2 000時(shí)，Inspect方法進(jìn)行100次模擬的平均運(yùn)算時(shí)間達(dá)到了13.7 s，而SOP-SeedBS方法的平均運(yùn)算時(shí)間僅為0.14 s，約是Inspect方法運(yùn)算時(shí)間的1%。當(dāng)p=500，n=500，Inspect方法進(jìn)行 100次模擬的平均運(yùn)算時(shí)間是13.6 s，而SOP-SeedBS方法的平均運(yùn)算時(shí)間僅為0.15 s，約是Inspect方法運(yùn)算時(shí)間的1%，由此能體現(xiàn)SOP-SeedBS方法在運(yùn)算時(shí)間上優(yōu)于Inspect方法。

為了研究變點(diǎn)間隔對(duì)實(shí)驗(yàn)結(jié)果準(zhǔn)確性的影響，從前面的對(duì)比實(shí)驗(yàn)表1中發(fā)現(xiàn)，當(dāng)均值變化為情形1時(shí)，SOP-SeedBS方法的效果不如情形2和情形3的情況，這也符合不同變點(diǎn)在不同稀疏維數(shù)坐標(biāo)中發(fā)生變化時(shí)，信號(hào)強(qiáng)度會(huì)減弱的預(yù)期效果。因此，表2顯示的是在情形1下的模擬結(jié)果，即發(fā)生變化的k個(gè)維數(shù)坐標(biāo)是來自完全不相交的坐標(biāo)集。表中Z1、Z2、Z3、Z4表示估計(jì)到的變點(diǎn)個(gè)數(shù)，取值分別為1，2，3，4。

表2 不同變點(diǎn)間隔在p=200，σ2=1時(shí)四組不同參數(shù)設(shè)置下100次模擬的結(jié)果

表2中有四組不同的參數(shù)設(shè)置B1,B2,B3,B4，其中B1的參數(shù)設(shè)置為zi+1-zi=20，表示相鄰變點(diǎn)之間間隔20個(gè)數(shù)據(jù)，信號(hào)的均值跳躍大小為=(0.8,1.6,2.4)；B2的參數(shù)設(shè)置為zi+1-zi=30,=(0.8,1.6,2.4)；B3參數(shù)設(shè)置為zi+1-zi=40,=(0.8,1.6,2.4)；B4的參數(shù)設(shè)置為zi+1-zi=50,=(0.8,1.6,2.4)。從表2的結(jié)果可以看出，通過模型1生成的四組不同參數(shù)設(shè)置下的變點(diǎn)數(shù)據(jù)，在σ2=1，p=200，n=2 000，k=40，變點(diǎn)間隔為20時(shí)，真實(shí)變點(diǎn)個(gè)數(shù)為3，三個(gè)變點(diǎn)的均值跳躍大小分別為 0.6、1.2、1.8，SOP-SeedBS方法在100次模擬中估計(jì)正確變點(diǎn)個(gè)數(shù)次數(shù)僅為 31 次，而在σ2=1、p=200、n=2 000、k=60時(shí)，100次模擬中估計(jì)正確變點(diǎn)個(gè)數(shù)次數(shù)僅有61次。由最后一組實(shí)驗(yàn)可以看出當(dāng)變點(diǎn)之間的間隔達(dá)到50個(gè)數(shù)據(jù)點(diǎn)及以上時(shí)，SOP-SeedBS方法不僅能穩(wěn)定地檢測(cè)到正確的變點(diǎn)個(gè)數(shù)，在準(zhǔn)確性上也能達(dá)到很好的效果。

3 實(shí)證分析

本節(jié)研究了R中ecp包的膀胱腫瘤微陣列數(shù)據(jù)集（Bladder Tumor Micro-Array Data），目的是檢驗(yàn)樣本目標(biāo)區(qū)域染色體拷貝數(shù)的變化，常常用于醫(yī)學(xué)中檢驗(yàn)癌癥的研究，染色體拷貝數(shù)變化反映了DNA的擴(kuò)增和驟減。從醫(yī)學(xué)的角度來說，擴(kuò)增的片段可能存在致癌基因，驟減的片段可能存在抑癌基因。本文對(duì)數(shù)據(jù)集中43個(gè)個(gè)體（p）的基因組上的2 215個(gè)不同基因位點(diǎn)（n）的對(duì)數(shù)強(qiáng)度比測(cè)量值進(jìn)行檢驗(yàn)，圖4繪制了前10個(gè)個(gè)體在前200個(gè)基因位點(diǎn)上的對(duì)數(shù)強(qiáng)度比散點(diǎn)圖，并采用SOP-SeedBS方法估計(jì)變點(diǎn)位置（其中虛線表示檢測(cè)到的變點(diǎn)位置），同時(shí)給出了文獻(xiàn)[10]和文獻(xiàn)[14]中Inspect方法和 Adapt-WBS方法的檢驗(yàn)結(jié)果，來驗(yàn)證SOP-SeedBS方法的有效性。

圖4 前10個(gè)個(gè)體在前200個(gè)位點(diǎn)上的ACGH數(shù)據(jù)散點(diǎn)圖

從圖4中可以看出，SOP-SeedBS方法檢測(cè)到的變點(diǎn)結(jié)果為：15，26，28，52（56），73，91，102，134，147（146），177（174），183（180），185（186），188（191），與文獻(xiàn)[10]和文獻(xiàn)[14]中檢驗(yàn)到的變點(diǎn)相近，其他無法檢測(cè)到的變點(diǎn)，例如位點(diǎn)30到40之間的數(shù)據(jù)出現(xiàn)的分段特征不明顯，許多個(gè)體在位點(diǎn)110～130之間測(cè)量值分布較為分散，沒有呈現(xiàn)出明顯的分段特征，位點(diǎn)155附近的測(cè)量值分布非常分散。SOP-SeedBS算法的結(jié)果與Inspect的結(jié)果的相似度約為80%。

4 結(jié)論與展望

本文提出SOP-SeedBS方法結(jié)合了奇異值分解和投影的思想，通過分析高維稀疏變點(diǎn)數(shù)據(jù)的特點(diǎn)，對(duì)數(shù)據(jù)矩陣進(jìn)行奇異值分解找到一個(gè)最優(yōu)的投影方向，將存在稀疏變點(diǎn)的高維數(shù)據(jù)沿著最優(yōu)的投影方向投影成一維數(shù)據(jù)，再利用SeedBS方法對(duì)一維數(shù)據(jù)進(jìn)行多變點(diǎn)檢驗(yàn)。通過仿真實(shí)驗(yàn)，分別用變點(diǎn)檢驗(yàn)頻數(shù)和平均Hausdorff距離評(píng)估了變點(diǎn)檢驗(yàn)的準(zhǔn)確度。本文提出的SOP-SeedBS方法在運(yùn)算時(shí)間、算法準(zhǔn)確度上都優(yōu)于Inspect方法，特別是在時(shí)間上具有較大的優(yōu)勢(shì)。最后將SOP-SeedBS方法應(yīng)用于ecp包中的膀胱腫瘤微陣列數(shù)據(jù)集中，檢驗(yàn)結(jié)果與文獻(xiàn)[14]中給出的Inspect方法的檢驗(yàn)結(jié)果有約80%的相似度。

本文是先將高維數(shù)據(jù)投影到一維，再進(jìn)行變點(diǎn)檢驗(yàn)。進(jìn)一步還可以研究能否在保持?jǐn)?shù)據(jù)維度不變的情況下，通過某種方式，例如引入一種針對(duì)高維數(shù)據(jù)的檢驗(yàn)統(tǒng)計(jì)量，將SeedBS方法直接用于高維數(shù)據(jù)的變點(diǎn)檢驗(yàn)問題中。