馮巖, 宋珊, 徐常青

(蘇州科技大學(xué)數(shù)學(xué)科學(xué)學(xué)院,江蘇 蘇州 215009)

1.引言

Rayleigh分布是一種連續(xù)的概率分布,它最早于1880年在聲學(xué)和光學(xué)領(lǐng)域由Lord Rayleigh提出并以他的名字命名.由于Rayleigh分布有良好的單調(diào)失效率性, 因此被廣泛地應(yīng)用到通訊工程、物理學(xué)、醫(yī)學(xué)、可再生能源、機(jī)械設(shè)計(jì)和航空航天等領(lǐng)域中.在通信理論中, 對(duì)到達(dá)接收器的密集散射信號(hào)的多條路徑進(jìn)行建模.在物理學(xué)中模擬風(fēng)速, 波高、聲輻射和光輻射.在醫(yī)學(xué)成像科學(xué)中, 要對(duì)磁共振成像中的噪聲方差建模.

Rayleigh分布在可靠性壽命試驗(yàn)方面也有許多實(shí)際應(yīng)用, 并且得出設(shè)備的生命周期取決于其年齡的結(jié)論.Rayleigh分布通常用于不同領(lǐng)域的物理建模過程, 如海浪的高度, 聲音和光線輻射, 無線電信號(hào)和風(fēng)能, 超聲圖像建模[1]等.它也被用于模型壽命只有數(shù)小時(shí)的管, 電阻, 網(wǎng)絡(luò), 晶體, 旋鈕, 變壓器, 飛機(jī)雷達(dá)組中的繼電器和電容器等.Rayleigh分布還用于研究風(fēng)力渦輪機(jī)站點(diǎn)一年以上的風(fēng)速和日平均風(fēng)速.在所有上述應(yīng)用中, 通常認(rèn)為特定設(shè)備的的壽命依賴于其年齡.由于基于壽命分布的廣義模型在建模過程中更具有直觀的適用性和吸引力, Surles和Padgett[2]提出了雙參數(shù)Burr Type X分布, 也稱為廣義Rayleigh(GR)分布.Surles等[3]和Al-khedhari等[4]討論了使用不同技術(shù)對(duì)該分布的參數(shù)估計(jì).

本文引入了一種新的二參數(shù)廣義Rayleigh (BVGR)分布, 其邊緣分布是GR分布.BVGR分布的邊緣分布的風(fēng)險(xiǎn)率函數(shù)可以是遞增或遞減、常數(shù)或浴盆形, 這一性質(zhì)豐富了BVGR分布的應(yīng)用.另外本文推導(dǎo)出了BVGR分布的矩生成函數(shù)和k階矩的表達(dá)式, 通過對(duì)樣本點(diǎn)進(jìn)行歸一化, 得到了它的近似累積分布函數(shù)(CDF）、概率密度函數(shù)(PDF)、生存函數(shù)(SF)、風(fēng)險(xiǎn)率函數(shù)(HRF)等.該模型未知參數(shù)的最大似然估計(jì)(MLEs)并不總是存在, 因此本文在對(duì)樣本點(diǎn)進(jìn)行歸一化的基礎(chǔ)上用EM算法非常有效地計(jì)算了參數(shù)的MLEs, 而且極大地簡(jiǎn)化了估計(jì)參數(shù)的表示形式.另外, 利用文[5]中Louis的思想, 根據(jù)EM算法獲得了觀測(cè)Fisher信息矩陣, 最后根據(jù)參數(shù)的MLEs漸進(jìn)服從正態(tài)分布, 構(gòu)造了未知參數(shù)的漸近置信區(qū)間.

2.二參數(shù)廣義Rayleigh分布

設(shè)Z1,Z2是相互獨(dú)立的隨機(jī)變量,它們都服從正態(tài)分布N(0,σ2),則隨機(jī)變量的CDF和PDF分別為

和

其中z ≥0, 稱Z服從參數(shù)為σ(σ >0)的Rayleigh分布.

參照Surles和Padgett[2]對(duì)GR分布的定義, 對(duì)Rayleigh分布進(jìn)行改進(jìn), 令=α, 并引入?yún)?shù)β, 則當(dāng)x>0時(shí), GR分布的CDF和PDF分別定義為

和

其中α>1,β >0.當(dāng)β=1,它與Rayleigh分布一致.為簡(jiǎn)便起見,把分布函數(shù)為式(2.3)的GR分布用GR(α,β)表示.

現(xiàn)在假設(shè)

是相互獨(dú)立的.定義X1= min(U0,U1),X2= min(U0,U2), 那么二元向量X= (X1,,X2)T的聯(lián)合分布稱之為二參數(shù)廣義Rayleigh 分布, 用BVGR(α,β0,β1,β2)表示, 其參數(shù)為α,β0,β1,β2.

如果X ～BVGR(α,β0,β1,β2), 對(duì)于x0=max{x1,x2}, 則它們的CDF用下式表示[6]:

即

X的聯(lián)合PDF可以寫成

其中

由于矩生成函數(shù)和矩函數(shù)在顏色傳輸、故障診斷、大型反射面天線仿真、信號(hào)和噪聲傳播、編程與算法設(shè)計(jì)和分析等方面發(fā)揮著重要的作用, 接下來我們將引入X的矩生成函數(shù)以及k階矩.令t=(t1,t2)T, 那么二元隨機(jī)向量X的矩生成函數(shù)為MX(t)=E[etTX], 即

其中

特別地, 當(dāng)α=2時(shí), 利用二項(xiàng)式級(jí)數(shù)展開上式可以寫為

參考文[7]中k階矩的概率, 可以推導(dǎo)出X的k階矩, 即

同樣為計(jì)算方便，我們令α=2, 則利用二項(xiàng)式級(jí)數(shù)展開上式可以寫為

利用Taylor展開, (2.5)可近似改為:

這里, 我們對(duì)樣本點(diǎn)采用球面歸一化, 即令+= 1.球面歸一化使(2.11)產(chǎn)生的誤差達(dá)到最小.注意到其他歸一化方法, 如1-范數(shù)歸一化:‖x‖1=|x1|+|x2|= 1, 或最大范數(shù)歸一化,即‖x‖∞≡max{|x1|,|x2|}=1, 對(duì)應(yīng)的相對(duì)誤差均大于球面歸一化情形.

接下來的討論中我們都假設(shè)對(duì)樣本進(jìn)行歸一化處理, 可以很容易地驗(yàn)證Xi,i=1,2的邊際分布遵循GR(α,β0+βi),即Xi的邊際CDF近似為

因此,Xi的邊際生存函數(shù)(SF)近似為

利用(2.12),Xi的邊際PDF近似為

而Xi的邊際失效函數(shù)(HRF)近似為

Raqab和Kundu[8]觀察到: 1)當(dāng)β0+βi ≤1/2時(shí), GR(α,β0+βi)的PDF呈下降趨勢(shì), HRF呈浴盆狀; 2)當(dāng)β0+βi >1/2時(shí), PDF呈右偏單峰趨勢(shì), HRF 呈上升趨勢(shì).與Marshal-Olkin(MO)分布和雙變量廣義指數(shù)(BVGE)分布相比, BVGR分布的邊緣分布的風(fēng)險(xiǎn)率函數(shù)可以是遞增或遞減、常數(shù)或浴盆形, 這一性質(zhì)豐富了BVGR分布的應(yīng)用.

若X ～BVGR(α,β0,β1,β2), 則min{X1,X2}的聯(lián)合CDF近似是

該結(jié)果表明, min{X1,X2}的SF近似為

根據(jù)式(2.5)、(2.12)以及X的聯(lián)合SF與聯(lián)合CDF的關(guān)系, 有

我們可以推出X的聯(lián)合SF近似為

利用聯(lián)合PDF和X的聯(lián)合SF之間的關(guān)系, 可以下式得出X的聯(lián)合HRF近似

3.BVGR分布的EM算法

在這一節(jié)中, 我們用EM算法來計(jì)算二參數(shù)廣義Rayleigh分布的未知參數(shù)的MLEs.假設(shè){(x11,x21),...,(x1n,x2n)}是來自BVGR(α,β0,β1,β2)的隨機(jī)二元樣本, 參照文[9], 引入下面的符號(hào)：

這里‖Ij‖對(duì)于j=0,1,2表示集合Ij中元素的個(gè)數(shù).當(dāng)i ∈I0時(shí), 其對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)近似為

當(dāng)i ∈I1時(shí), 其對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)近似為

當(dāng)i ∈I2時(shí), 其對(duì)應(yīng)的對(duì)數(shù)似然函數(shù)近似為

文[10]指出, 當(dāng)有一個(gè)ni= 0時(shí), MLEs不存在.如果n0>0,n1>0,n2>0, 則MLEs存在, 分別為:β0,β1,β2, 但沒有顯式表達(dá)式.在本文中我們使用EM算法來計(jì)算BVGR模型中未知參數(shù)的MLEs.因此, 估計(jì)的BVGR可以視為一個(gè)缺失數(shù)據(jù)的問題.假設(shè)對(duì)于二元隨機(jī)向量(X1,X2), 存在一個(gè)相關(guān)的隨機(jī)向量(δ1,δ2), 其中(δ1,δ2)定義為:

和

這里的Ui與第2節(jié)開始時(shí)定義的相同.可以看出, 在這種情況下, 即使我們知道(X1,X2), 但是對(duì)應(yīng)的(δ1,δ2)可能并不總是已知的.例如X1=X2, 則δ1=δ2= 0已知；但是, 如果X12,則(δ1,δ2)未知.如果(x1,x2)∈I1, 則(δ1,δ2)的可能值為(1,0)或(1,2), 同理, 如果(x1,x2)∈I2,則(δ1,δ2)的可能值為(0,2)或(1,2), 概率不為零.

如果(X1,X2)和相關(guān)的(δ1,δ2)對(duì)于所有的觀測(cè)都是已知的, 那么通過求解一維優(yōu)化問題,可以很容易地得到.但是對(duì)于所有的觀測(cè)結(jié)果(δ1,δ2)并不都是已知的.因此我們可以用EM算法計(jì)算未知參數(shù)的MLEs, 首先是“E”步, 如果對(duì)應(yīng)的(δ1,δ2)缺失, 則可以對(duì)對(duì)數(shù)似然函數(shù)(3.1) 求期望值, 形成偽對(duì)數(shù)似然函數(shù).“M”步驟是通過最大化偽對(duì)數(shù)似然函數(shù)得到未知參數(shù), 步驟如下.

在“E”步驟中, 當(dāng)觀測(cè)值屬于I0時(shí), 對(duì)應(yīng)的(δ1,δ2)是完全已知的, 因此可以保持所有觀測(cè)的對(duì)數(shù)似然函數(shù)不變.如果觀測(cè)值屬于I1或I2, 我們將其視為缺失觀測(cè)值.如果(x1,x2)∈I1, 我們形成偽觀測(cè),通過分式(x1,x2)分為兩個(gè)部分完整的偽觀測(cè),形式為(x1,x2,u1(θ))和(x1,x2,u2(θ)).這里, 給定X1< X2,θ= (α,β0,β1,β2),u1(θ)和u2(θ)分別為(δ1,δ2)取值(1,0)和(1,2)的條件概率.同理,若(x1,x2)∈I2,則偽觀測(cè)的形式為(x1,x2,v1(θ))和(x1,x2,v2(θ)),其中v1(θ)和v2(θ)分別為(δ1,δ2)取值(0,2)或(1,2)的條件概率, 給定X1>X2.因?yàn)?/p>

且

因此根據(jù)條件概率公式有

同理

現(xiàn)在令μ1=μ1(θ),μ2=μ2(θ),ν1=ν1(θ),ν2=ν2(θ), 則“偽數(shù)據(jù)”的對(duì)數(shù)似然函數(shù)可以寫成

通過簡(jiǎn)化可以得到

“M”步驟是關(guān)于α,β0,β1,β2最大化(3.4).首先固定α, 關(guān)于β0,β1,β2最大化(3.4)可以得到

可以得到

4.觀測(cè)Fisher信息矩陣及其置信區(qū)間

假設(shè)文[11]中的規(guī)則性條件成立, 這保證了Fisher信息矩陣的存在性.記M為完全數(shù)據(jù),X為觀測(cè)數(shù)據(jù),IM(θ)為完全信息矩陣,IX(θ)為觀測(cè)信息矩陣,IM|X(θ)為遺失信息矩陣, 根據(jù)文[2]提出的遺失信息原則, 得到觀測(cè)信息矩陣為

要計(jì)算EM中的Fisher信息矩陣, 可令S = (Si)為偽對(duì)數(shù)似然函數(shù)關(guān)于θ的梯度向量, B =(Bij)是其關(guān)于θ的負(fù)的二階導(dǎo)數(shù)矩陣.則

因此得到觀測(cè)Fisher信息矩陣為

其中矩陣B的元素為

向量S的元素為

注意, 所有這些條件期望值都可以在EM算法中僅使用S和B來計(jì)算, 僅需要在EM過程進(jìn)行最后一次迭代時(shí)對(duì)參數(shù)進(jìn)行估計(jì).由于極大似然估計(jì)是漸近服從正態(tài)分布, 期望為(α,β0,β1,β2), 協(xié)方差矩陣為, 這里是觀測(cè)信息矩陣的逆矩陣, 于是有

因此, 參數(shù)α,β0,β1,β2的置信水平為100(1?γ)%的漸近置信區(qū)間分別為

其中U11,U22,U33,U44分別是協(xié)方差矩陣的主對(duì)角線的元素:zγ/2是標(biāo)準(zhǔn)正態(tài)分布的γ/2上分位數(shù).

5.結(jié)論

本文在廣義Rayleigh分布的基礎(chǔ)上引入了BVGR分布, 并且推導(dǎo)出了其矩生成函數(shù)和k階矩的表達(dá)式, 在對(duì)樣本點(diǎn)進(jìn)行歸一化的前提下, 得到了它的近似CDF、PDF、SF、HRF等.在此基礎(chǔ)上, 用EM算法來計(jì)算未知參數(shù)的MLEs, 極大地簡(jiǎn)化了估計(jì)參數(shù)的表示形式.同時(shí), 還得到了觀測(cè)Fisher信息矩陣, 最后根據(jù)參數(shù)MLEs漸進(jìn)服從正態(tài)分布, 構(gòu)造了置信水平為100(1?γ)%的漸近置信區(qū)間.