基于差分進(jìn)化算法的正態(tài)分布均值變點(diǎn)檢測(cè)

朱偉業(yè)成，金良瓊，沈 婷

(貴州民族大學(xué) 數(shù)據(jù)科學(xué)與信息工程學(xué)院，貴陽(yáng) 550025)

變點(diǎn)檢測(cè)問題自20世紀(jì)50年代被提出以來(lái)，一直是統(tǒng)計(jì)領(lǐng)域當(dāng)中的一個(gè)熱點(diǎn)問題。變點(diǎn)指的是從某個(gè)時(shí)刻開始，樣本的分布或者數(shù)字特征發(fā)生了變化。變點(diǎn)檢測(cè)就是要檢測(cè)數(shù)據(jù)中是否存在變點(diǎn)及對(duì)變點(diǎn)的個(gè)數(shù)和位置進(jìn)行估計(jì)。隨著社會(huì)的發(fā)展，變點(diǎn)檢測(cè)理論已廣泛運(yùn)用在經(jīng)濟(jì)、醫(yī)學(xué)、氣象和圖像處理等領(lǐng)域。

在變點(diǎn)檢測(cè)問題中，均值變點(diǎn)是一類重要的變點(diǎn)類型，數(shù)據(jù)的均值發(fā)生變化往往引起人們的重視。針對(duì)正態(tài)分布序列的均值變點(diǎn)問題，中外許多學(xué)者都對(duì)此進(jìn)行過研究。Hawkins[1]應(yīng)用極大似然法研究正態(tài)序列均值變點(diǎn)問題。Sen等[2]利用貝葉斯方法給出了均值變點(diǎn)的精確和漸進(jìn)的分布函數(shù)。Kim等[3]采用一種退火隨機(jī)逼近蒙特卡洛法（ASAMC）對(duì)韓國(guó)傳染病數(shù)據(jù)進(jìn)行均值變點(diǎn)研究。趙江南等[4]基于ASAMC方法對(duì)正態(tài)分布序列均值變點(diǎn)問題進(jìn)行研究并將結(jié)果應(yīng)用于烏魯木齊市的溫度變點(diǎn)檢測(cè)問題。胡丹青等[5]在貝葉斯后驗(yàn)分布的基礎(chǔ)上應(yīng)用遺傳算法，研究了線性回歸模型中的多結(jié)構(gòu)變點(diǎn)問題。王曉原等[6]運(yùn)用加速遺傳算法結(jié)合均值變點(diǎn)模型對(duì)交通流變點(diǎn)問題進(jìn)行分析研究。郭衛(wèi)娟[7]采用二分法對(duì)方差已知且相等的正態(tài)分布均值多變點(diǎn)問題進(jìn)行研究，給出了變點(diǎn)位置的后驗(yàn)分布?？姲仄涞萚8]基于信息論準(zhǔn)則，研究了異方差情況下的多元正態(tài)分布均值變點(diǎn)問題。何朝兵等[9]采用馬爾科夫蒙特卡洛方法研究了左截?cái)嘤覄h失情況下的對(duì)數(shù)正態(tài)分布序列的變點(diǎn)檢測(cè)問題。楊豐凱等[10]采用一種非迭代抽樣方法（IBF）與迭代Gibbs抽樣算法做對(duì)比，研究了在無(wú)信息先驗(yàn)條件下的正態(tài)分布均值單變點(diǎn)問題。對(duì)于方差已知且相等的正態(tài)分布序列均值變點(diǎn)模型，貝葉斯法是一種常用的研究方法，但貝葉斯法存在的不足之處就是后驗(yàn)分布的計(jì)算量較大?；诖?，本文將差分進(jìn)化算法與貝葉斯法結(jié)合，對(duì)正態(tài)分布序列均值變點(diǎn)問題進(jìn)行研究，模擬結(jié)果表明該方法能快速有效地計(jì)算后驗(yàn)分布的最大值，并且檢測(cè)出均值變點(diǎn)的位置。

1 正態(tài)分布均值多變點(diǎn)模型

式中：μ1≠μ2，μ2≠μ3…μk≠μk+1，σ2為已知常數(shù)，則稱該模型為正態(tài)分布序列均值變點(diǎn)模型，k為變點(diǎn)個(gè)數(shù)，r1，r2，r3...rk為變點(diǎn)位置。當(dāng)k=1時(shí)，該模型為單變點(diǎn)模型，k>1時(shí)，為多變點(diǎn)模型。本文考慮的問題是在k已知的情況下，如何估計(jì)變點(diǎn)位置，即r1，r2，r3...rk的值。

針對(duì)變點(diǎn)問題的研究方法主要有（加權(quán)）最小二乘法、極大似然法、非參數(shù)方法和貝葉斯法等。貝葉斯法利用貝葉斯理論，對(duì)包括變點(diǎn)位置在內(nèi)的未知參數(shù)選取合適的先驗(yàn)分布，利用似然函數(shù)推導(dǎo)變點(diǎn)位置的后驗(yàn)分布，結(jié)合樣本數(shù)據(jù)計(jì)算分布的最大值，對(duì)變點(diǎn)的位置進(jìn)行估計(jì)。運(yùn)用貝葉斯法，如何選取參數(shù)的先驗(yàn)分布是關(guān)鍵，本文選擇的是無(wú)信息先驗(yàn)分布。

2 貝葉斯估計(jì)

2.1 變點(diǎn)位置的無(wú)信息先驗(yàn)分布

正態(tài)分布序列均值變點(diǎn)模型中的未知參數(shù)為變點(diǎn)位置r1，r2，r3…rk和分布的均值μ1，μ2…μk+1。無(wú)信息先驗(yàn)分布可以理解為對(duì)參數(shù)的任何可能值既無(wú)偏愛，又同等無(wú)知[11]。在1，2…n-1上每個(gè)變點(diǎn)序列出現(xiàn)的概率都是相等的，因此本文選取r1，r2，r3…rk的聯(lián)合先驗(yàn)分布為

2.2 正態(tài)均值的無(wú)信息先驗(yàn)分布[11]

正態(tài)分布中的均值μ為位置參數(shù)，總體X的密度函數(shù)具有形式p（x-μ）。為求均值μ的無(wú)信息先驗(yàn)分布，對(duì)X做平移變換，得到Y(jié)＝X＋c，對(duì)參數(shù)μ也做同樣的變換，得到η＝μ+c。設(shè)分布Y有密度函數(shù)p（y-η），此時(shí)，μ和η應(yīng)有相同的先驗(yàn)分布。即

π（·）和π*（·）分別為μ和η的無(wú)信息先驗(yàn)分布。此外，由變換η＝μ＋c可得到η的無(wú)信息先驗(yàn)分布為

比較式（3）和式（4）可得π（η）=π（η－c），由于η與c的任意性，可得均值μ的無(wú)信息先驗(yàn)分布為

這是正態(tài)均值μ的一個(gè)廣義先驗(yàn)分布。

2.3 變點(diǎn)位置的后驗(yàn)分布

設(shè)參數(shù)θ=（r1，r2，…rk，μ1，μ2…μk+1），假設(shè)變點(diǎn)位置與均值參數(shù)互相獨(dú)立，即

利用貝葉斯公式

可得θ的后驗(yàn)分布為

對(duì)式（8）中的μ1，μ2…μk+1積分，可得r1，r2，…rk的后 驗(yàn)分布為

為計(jì)算方便取對(duì)數(shù)，得

在多變點(diǎn)問題中，變點(diǎn)個(gè)數(shù)越多，變點(diǎn)位置的參數(shù)空間越復(fù)雜，計(jì)算后驗(yàn)分布最大值的計(jì)算量也變得越大。為快速找到該后驗(yàn)分布的最大值，本文引用差分進(jìn)化算法進(jìn)行分析。接下來(lái)介紹差分進(jìn)化算法的原理和具體步驟。

3 差分進(jìn)化算法介紹

差分進(jìn)化算法（DE）是一種高效的全局優(yōu)化算法。該算法模擬生物進(jìn)化的過程，首先隨機(jī)生成初代種群，然后通過變異、交叉操作生成子代種群，根據(jù)“優(yōu)勝劣汰”原則，淘汰適應(yīng)度低的種群，保留適應(yīng)度高的種群。經(jīng)過不斷迭代進(jìn)化，最終收斂到最優(yōu)結(jié)果。差分進(jìn)化算法具備結(jié)構(gòu)簡(jiǎn)單、收斂快和自適應(yīng)能力強(qiáng)等優(yōu)點(diǎn)，被廣泛運(yùn)用在各個(gè)領(lǐng)域。具體算法流程如下[12-13]。

3.1 初始化種群

首先確定各個(gè)參數(shù)，包括迭代次數(shù)G、種群大小NP、種群維數(shù)D、變異算子F、交叉算子CR，以及搜索空間的上下限xmax和xmin。然后隨機(jī)生成NP個(gè)長(zhǎng)度為D的解向量xi，j，i=1，2…NP，j=1，2…D，xmin≤xi，j≤xmax。種群大小NP一般選取[50，200]。

3.2 變異操作

式中：i，k1，k2，k3為兩兩互不相等的隨機(jī)整數(shù)表示第k個(gè)種群的第g代。在變異操作中，變異算子F取值范圍為[0，2]，F(xiàn)過小可能陷入局部最優(yōu)，F(xiàn)過大則不容易收斂，一般取[0.4，1]居多。另外，如果變異以后的值νi，j超出了邊界，可以隨機(jī)在范圍內(nèi)再選擇1個(gè)數(shù)，或者直接取邊界值。

3.3 交叉操作

式中：rand（0，1）是1個(gè)0到1的隨機(jī)數(shù)，變異算子CR用來(lái)控制子代向量值是變異向量值還是原向量值，jrand是1個(gè)0到D的隨機(jī)整數(shù)，其保證了子代向量至少有1個(gè)元素與原向量不同。

3.4 選擇操作

選擇操作的原理是計(jì)算子代向量和原向量各自的適應(yīng)度，如果子代向量的適應(yīng)度更高，則將子代向量保留，替代原向量進(jìn)入下一代循環(huán)，否則保留原向量。具體表達(dá)式如下

式中：f（·）為向量對(duì)應(yīng)的適應(yīng)度函數(shù)。

3.5 退出循環(huán)

當(dāng)最后的解向量滿足所需要的精度，或循環(huán)次數(shù)達(dá)到預(yù)設(shè)值G后，退出循環(huán)，否則重復(fù)步驟3.2到步驟3.4。

4 數(shù)值模擬

在變點(diǎn)個(gè)數(shù)已知的情況下，對(duì)正態(tài)分布序列均值變點(diǎn)的位置進(jìn)行檢測(cè)?？紤]變點(diǎn)個(gè)數(shù)為2的情況，模型如下

根據(jù)上述模型，分別生成樣本量為300、400、500、600的數(shù)據(jù)作為檢測(cè)的樣本。對(duì)應(yīng)的變點(diǎn)位置分別為100和200、150和300、200和350、200和400。采用2種差分進(jìn)化算法進(jìn)行估計(jì)：差分進(jìn)化算法和自適應(yīng)差分進(jìn)化算法。兩者的不同之處在于自適應(yīng)差分進(jìn)化算法將變異算子F設(shè)為1個(gè)隨迭代次數(shù)變化的值，在迭代初期，F(xiàn)較大，能保證種群多樣性；迭代后期，F(xiàn)較小，能提高搜索效率。算法的參數(shù)選擇：種群數(shù)NP＝50，最大迭代次數(shù)Gm＝20，變異算子F0＝0.5，交叉算子CR=0.4，自適應(yīng)變異算子由式（15）產(chǎn)生

式中：G表示當(dāng)前迭代的次數(shù)。使用Matlab軟件進(jìn)行數(shù)值模擬1 000次，結(jié)果見表1、表2。

表1 不同樣本量下差分進(jìn)化算法的r1、r2識(shí)別結(jié)果

表2 不同樣本量下自適應(yīng)差分進(jìn)化算法的r1、r2識(shí)別結(jié)果

從表1和表2可以看出，在樣本量分別為300、400、500、600的情況下，迭代20次，得到的變點(diǎn)位置r1、r2的估計(jì)值與真實(shí)值的相對(duì)誤差均小于0.05%，標(biāo)準(zhǔn)差均小于1.3。說明2種方法均能夠較好地識(shí)別出變點(diǎn)位置r1、r2。對(duì)比2種方法估計(jì)結(jié)果的相對(duì)誤差與標(biāo)準(zhǔn)差，除了在樣本量為400的情況下，差分進(jìn)化算法對(duì)于r1的估計(jì)相對(duì)誤差為0.007%，大于自適應(yīng)差分進(jìn)化算法對(duì)應(yīng)的相對(duì)誤差0.005%以外，其余情況差分進(jìn)化算法的估計(jì)結(jié)果的相對(duì)誤差與標(biāo)準(zhǔn)差都小于自適應(yīng)差分進(jìn)化算法，整體估計(jì)效果更好。

在樣本量為300、400、500、600的情況下，統(tǒng)計(jì)2種方法正確識(shí)別變點(diǎn)位置r1、r2及正負(fù)1個(gè)或2個(gè)單位偏差的次數(shù)，對(duì)比評(píng)估2種方法的識(shí)別率，結(jié)果見表3、表4。

表3 不同樣本量下基本差分進(jìn)化算法的r1、r2識(shí)別率

表4 不同樣本量下自適應(yīng)差分進(jìn)化算法的r1、r2識(shí)別率

由表3和表4可以看出，除了在樣本量n=300的情況下，差分進(jìn)化算法對(duì)第二個(gè)變點(diǎn)正負(fù)1個(gè)單位偏差的識(shí)別率為96.7%，低于自適應(yīng)差分進(jìn)化算法的97%。其余情況下，差分進(jìn)化算法的識(shí)別率均大于等于自適應(yīng)差分進(jìn)化算法，進(jìn)一步說明前者的識(shí)別效果更好。從表中還看出r2的識(shí)別率大于r1的識(shí)別率，這是因?yàn)閞2前后的均值跳躍度要大于r2前后的均值跳躍度，2組數(shù)據(jù)的差異更明顯，因此更容易識(shí)別出來(lái)，符合常理。

為更加直觀地觀察估計(jì)效果，以樣本量為600的實(shí)驗(yàn)結(jié)果為例，繪制r1、r2的散點(diǎn)圖如圖1—圖4所示。

圖1 差分進(jìn)化中r1的散點(diǎn)圖

圖2 自適應(yīng)差分進(jìn)化中r1的散點(diǎn)圖

圖3 差分進(jìn)化中r2的散點(diǎn)圖

圖4 自適應(yīng)差分進(jìn)化中r2的散點(diǎn)圖

由圖1—圖4可以明顯看出，2種方法對(duì)r1、r2的識(shí)別結(jié)果集中在各自對(duì)應(yīng)的真值附近。其中，對(duì)r2的識(shí)別結(jié)果比r1更為集中，效果更好。

為了進(jìn)一步對(duì)比2種方法的優(yōu)劣，繪制迭代圖如圖5、圖6所示。

圖5 差分進(jìn)化算法迭代圖

圖6 自適應(yīng)差分進(jìn)化算法迭代圖

從圖5和圖6可以看出，差分進(jìn)化算法在迭代到第16代的時(shí)候就已經(jīng)搜索到了最優(yōu)解，但自適應(yīng)差分進(jìn)化算法到第20代的時(shí)候才搜索到最優(yōu)解，這說明差分進(jìn)化算法能夠更快地收斂到函數(shù)的最大值。為探究迭代次數(shù)對(duì)2種方法估計(jì)效果的影響。下面將最大迭代次數(shù)Gm增加到100次，以樣本量為600的變點(diǎn)模型為例，進(jìn)行數(shù)值模擬，結(jié)果見表5、圖7、圖8。

表5 最大迭代次數(shù)Gm=100下的2種算法的識(shí)別結(jié)果

圖7 Gm=100時(shí)差分進(jìn)化算法的收斂圖

圖8 Gm=100時(shí)自適應(yīng)差分進(jìn)化算法的收斂圖

由表5、圖7、圖8可知，在迭代次數(shù)足夠大的情況下，2種算法對(duì)變點(diǎn)位置的檢測(cè)效果相同。

由此可知：在迭代次數(shù)較小時(shí)，差分進(jìn)化算法的檢測(cè)效果要優(yōu)于自適應(yīng)差分進(jìn)化算法。當(dāng)?shù)螖?shù)達(dá)到一定的程度時(shí)，2個(gè)方法的檢測(cè)效果一致。

為了觀察變點(diǎn)檢測(cè)前后序列的均值特征，繪制變點(diǎn)檢測(cè)前后的均值特征圖如圖9、圖10所示。

圖9 變點(diǎn)檢測(cè)前的均值特征

圖10 變點(diǎn)檢測(cè)后的均值特征

圖9的虛線表示無(wú)變點(diǎn)模型的數(shù)據(jù)均值特征，圖10的虛線表示在變點(diǎn)模型下的數(shù)據(jù)均值特征?？梢钥闯鲎凕c(diǎn)模型將數(shù)據(jù)分成3段，每段分別刻畫數(shù)據(jù)的均值。更符合數(shù)據(jù)的真實(shí)情況。這說明了本文引用的差分進(jìn)化算法能夠檢測(cè)出變點(diǎn)的位置從而更好地解釋正態(tài)分布均值的特征。

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

本文在貝葉斯理論的基礎(chǔ)上，對(duì)已知變點(diǎn)個(gè)數(shù)條件下的正態(tài)分布序列均值變點(diǎn)模型進(jìn)行研究。首先選擇無(wú)信息先驗(yàn)作為變點(diǎn)位置和正態(tài)均值的先驗(yàn)分布，利用貝葉斯公式得到變點(diǎn)位置的后驗(yàn)分布。為快速計(jì)算后驗(yàn)分布的最大值，引入2種差分進(jìn)化算法結(jié)合后驗(yàn)分布進(jìn)行優(yōu)化并進(jìn)行數(shù)值模擬。模擬結(jié)果表明：2種算法均能夠有效識(shí)別出變點(diǎn)的位置。其中，差分進(jìn)化算法的效果要優(yōu)于自適應(yīng)差分進(jìn)化算法，因此在選擇方法時(shí)，應(yīng)選擇差分進(jìn)化算法估計(jì)變點(diǎn)。綜上，使用差分進(jìn)化算法結(jié)合貝葉斯法能夠更有效地檢測(cè)變點(diǎn)個(gè)數(shù)已知下的正態(tài)分布序列均值變點(diǎn)模型中的變點(diǎn)位置。