沙雪云， 周菊玲， 董翠玲

(新疆師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 新疆 烏魯木齊 830017)

0 引言

變點(diǎn)問題是近幾年統(tǒng)計(jì)學(xué)的熱點(diǎn)研究問題，變點(diǎn)模型則是研究變點(diǎn)問題的一種非常重要的統(tǒng)計(jì)模型，其應(yīng)用廣泛，研究方法多樣. 常用的研究方法有貝葉斯方法、極大似然法和最小二乘法等，其中貝葉斯方法在解決變點(diǎn)問題上綜合了先驗(yàn)信息以及樣本信息，使得判斷更為準(zhǔn)確. MCMC算法是一種高效的貝葉斯方法，將Gibbs抽樣與M-H抽樣相結(jié)合的算法，根據(jù)參數(shù)的滿條件分布形式來選取Gibbs抽樣和M-H抽樣，進(jìn)而得到參數(shù)的Gibbs樣本，最終把Gibbs樣本的均值作為各參數(shù)的貝葉斯估計(jì). Lomax分布是一種非常重要的壽命分布，具有單調(diào)的失效率，所以在壽命試驗(yàn)的數(shù)據(jù)處理中起著至關(guān)重要的作用. 姚惠等人分別研究了基于沒有任何數(shù)據(jù)缺失的情況下，Lomax分布在熵?fù)p失函數(shù)、Linex損失函數(shù)等不同的損失函數(shù)下參數(shù)的bayes估計(jì)[1-4]. 龍兵等人針對在數(shù)據(jù)不完整的情況下對Lomax分布進(jìn)行各參數(shù)估計(jì)[5-8]. 岑泰林討論了在完全數(shù)據(jù)和隨機(jī)刪失數(shù)據(jù)不同情況下Lomax分布的參數(shù)估計(jì)問題[9]. 但是，目前對Lomax分布單變點(diǎn)問題的研究較少. 本文將研究在尺度參數(shù)已知的情況下，對Lomax分布的形狀參數(shù)及變點(diǎn)位置進(jìn)行貝葉斯估計(jì)，并運(yùn)用MCMC算法進(jìn)行隨機(jī)模擬，結(jié)果表明，各參數(shù)的估計(jì)值與真實(shí)值的MC誤差較小，說明各參數(shù)估計(jì)值在較高水平上是行之有效的.

1 Lomax分布及變點(diǎn)模型

設(shè)隨機(jī)變量X服從參數(shù)為λ,θ的Lomax分布，則分布函數(shù)和密度函數(shù)分別為

其中尺度參數(shù)λ>0，形狀參數(shù)θ>0.

假設(shè)樣本yi(i=1,2,…,n)值服從Lomax分布，若在序列{y1,y2,…,yn}存在變點(diǎn)，則在該序列中存在一個(gè)時(shí)間點(diǎn)，在該時(shí)間點(diǎn)的前后樣本值yi所服從的Lomax分布的尺度參數(shù)θ也將會發(fā)生變化，即在這一時(shí)間點(diǎn)之前的序列服從參數(shù)θ1為的Lomax分布，在這一時(shí)間點(diǎn)之后的序列服從參數(shù)為θ2的Lomax分布，稱該時(shí)間點(diǎn)為序列中的一個(gè)變點(diǎn)，當(dāng)在序列中存在兩個(gè)及以上的變點(diǎn)時(shí)，則稱此模型為多變點(diǎn)模型. 含有k個(gè)未知變點(diǎn)的Lomax分布的模型為

(1)

其中m1,m2,…,mk是所要求得未知變點(diǎn)參數(shù)，θ1,θ2,…,θk分別為不含變點(diǎn)區(qū)間的尺度參數(shù). 針對含有多個(gè)未知變點(diǎn)的情況，可以采用二分分段法和貝葉斯因子進(jìn)行逐段尋找，其中用貝葉斯因子來判斷序列是否存在單個(gè)變點(diǎn)可以借助于Kass和Raftery在1995年提出的貝葉斯因子模型選擇的一般準(zhǔn)則[10]. 下面給出用二分分段法及貝葉斯因子解決多變點(diǎn)問題的具體思路：

在處理多變點(diǎn)問題時(shí)，首先需要確定變點(diǎn)是否存在，若存在變點(diǎn)，則要在給定變點(diǎn)個(gè)數(shù)的情況下求出變點(diǎn)的所在位置，借助二分分段法及貝葉斯因子先來確定變點(diǎn)的個(gè)數(shù)，然后判斷在一個(gè)序列集中是否存在一個(gè)變點(diǎn)及這個(gè)變點(diǎn)所在的位置.

第1步：在整個(gè)序列集S中，利用貝葉斯因子來檢測是否存在單個(gè)變點(diǎn). 若貝葉斯因子檢測得不存在變點(diǎn)，就表示整個(gè)序列集不存在變點(diǎn)，則終止程序；否則就能夠檢測出這個(gè)序列里的第一個(gè)變點(diǎn)；

第2步：基于檢測出的第一個(gè)變點(diǎn)，將整個(gè)序列集S分為兩個(gè)子序列，對每個(gè)子序列，用第1步分別尋找出一個(gè)變點(diǎn)，一直持續(xù)這個(gè)過程直到在每個(gè)子序列中找不到變點(diǎn)為止.

利用二分分段法，只需要檢測并估計(jì)出單個(gè)變點(diǎn)的位置，接下來，利用貝葉斯方法來討論Lomax分布的單變點(diǎn)問題.

假設(shè)變點(diǎn)個(gè)數(shù)有且只有一個(gè)，Lomax分布形狀參數(shù)單變點(diǎn)模型為

(2)

其中x>0,θ1>0,θ2>0,λ=c(c為常數(shù))，且m,θ1,θ2均為未知. 當(dāng)θ1≠θ2時(shí)，m=2,3,…,n-1就是所要研究的變點(diǎn)，同時(shí)還需要對變點(diǎn)m以及參數(shù)θ1,θ2進(jìn)行貝葉斯估計(jì).

2 單變點(diǎn)模型參數(shù)的貝葉斯估計(jì)

當(dāng)θ1≠θ2時(shí)，設(shè)m為變點(diǎn)，記α=(m,θ1,θ2)，則此變點(diǎn)問題的似然函數(shù)為

1)θ1,θ2通過用Fisher信息陣確定無信息先驗(yàn)分布.

首先, 寫出樣本似然函數(shù)為

通過樣本似然函數(shù)可得到與其相關(guān)的樣本的信息矩陣,

進(jìn)而得到θ1,θ2的無信息先驗(yàn)矩陣為,

最后θ1,θ2的無信息先驗(yàn)分布為

由貝葉斯公式求得聯(lián)合后驗(yàn)分布為,

π(α|x)∝L(α|x)π(α)=

各參數(shù)滿條件分布為,

2) 取θ1,θ2的共軛先驗(yàn)分布為Ga(bi,ci)，即

且m,θ1,θ2相互獨(dú)立，由貝葉斯公式得聯(lián)合后驗(yàn)分布為,

π(α|x)∝L((α|x)π(m)π(θ1)π(θ2))∝

可得各參數(shù)滿條件分布為,

3 隨機(jī)模擬

下面進(jìn)行隨機(jī)模擬估計(jì)，通過上述分別得到了在無信息先驗(yàn)分布和共軛先驗(yàn)分布下的Lomax分布變點(diǎn)位置及各參數(shù)的滿條件分布，可以利用MCMC算法對變點(diǎn)位置及各參數(shù)進(jìn)行隨機(jī)模擬估計(jì)，由于參數(shù)θ1,θ2,m的滿條件分布都較為復(fù)雜，用Gibbs方法對樣本進(jìn)行抽樣較為困難，所以利用M-H方法對各參數(shù)θ1,θ2,m滿條件分布進(jìn)行抽樣.

取隨機(jī)產(chǎn)生n=183個(gè)數(shù)據(jù)作為樣本，參數(shù)(m,θ1,θ2,λ)的真實(shí)值取(100,3,18,5)，即，

(3)

利用得到的各參數(shù)滿條件分布，運(yùn)用MATLAB軟件來進(jìn)行MCMC模擬，在模擬過程中檢驗(yàn)?zāi)Ｐ蛥?shù)的收斂性，通過輸入多組初始值，形成多層Markov鏈，再用MATLAB軟件對參數(shù)進(jìn)行多層Markov鏈迭代分析，當(dāng)參數(shù)m,θ1,θ2收斂時(shí)，Markov鏈迭代圖將趨于重合. 為確保參數(shù)的收斂性，先進(jìn)行10 000次的預(yù)迭代的基礎(chǔ)上在進(jìn)行20 000次迭代，即從第10 001次開始至30 000次開始迭代，可得MATLAB的運(yùn)行結(jié)果，如表1所示.

當(dāng)θ1,θ2的先驗(yàn)分布為無信息先驗(yàn)分布時(shí):

表1 無信息先驗(yàn)分布下，參數(shù)m,θ1,θ2的貝葉斯估計(jì)

當(dāng)θ1,θ2取共軛先驗(yàn)分布時(shí)，θ1～Ga(15,3.2)，θ2～Ga(27,1.2)，

表2 共軛先驗(yàn)分布下，參數(shù)m,θ1,θ2的貝葉斯估計(jì)

對模擬結(jié)果進(jìn)行分析，首先，通過表1、表2可以觀察到對參數(shù)選取不同的先驗(yàn)分布后，分別對它們進(jìn)行隨機(jī)模擬，得到的估計(jì)值與真實(shí)值的MC誤差均不超過2%，說明各參數(shù)的估計(jì)值在較高水平上是有效的；其次，可取[2.5%分位數(shù)，97.5%分位數(shù)]作為參數(shù)置信水平為0.95的置信區(qū)間，從表1、表2可以看出置信區(qū)間較窄，區(qū)間估計(jì)效果良好；最后，通過圖2、圖4可以看出變點(diǎn)m的兩條Markov鏈圖像趨于重合，收斂性較好. 綜上分析，利用MCMC算法得到的參數(shù)估計(jì)及變點(diǎn)位置估計(jì)的效果都較為理想，因此使用該方法處理Lomax分布的變點(diǎn)問題是可行且有效的.