0–1 膨脹幾何分布的客觀貝葉斯分析

肖 翔，古 晞

（1.上海工程技術(shù)大學(xué) 數(shù)理與統(tǒng)計(jì)學(xué)院，上海 201620；2.同濟(jì)大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，上海 200092）

在醫(yī)療衛(wèi)生、保險(xiǎn)金融及可靠性等許多實(shí)際應(yīng)用領(lǐng)域，樣本數(shù)據(jù)不僅會(huì)出現(xiàn)零過(guò)多的情況，也會(huì)出現(xiàn)一過(guò)多的情況.例如，在新型冠狀病毒（COVID–19）大流行時(shí)，個(gè)體感染COVID–19 后，自身就會(huì)產(chǎn)生抗體，使其感染次數(shù)最多可能一次.又如，在商場(chǎng)買衣服時(shí)，出于貨比三家的心理，很多顧客沒(méi)有購(gòu)買衣服或者只購(gòu)買一件衣服.

近年來(lái)，國(guó)內(nèi)外很多文獻(xiàn)對(duì)0–1 膨脹泊松分布模型進(jìn)行了深入研究，取得豐富的研究成果.田震[1]基于數(shù)據(jù)刪失和加權(quán)擾動(dòng)模型對(duì)0–1 膨脹泊松分布模型進(jìn)行統(tǒng)計(jì)推斷.Tang 等[2]構(gòu)造0–1 膨脹泊松分布模型的等價(jià)表達(dá)式，采用極大似然估計(jì)與貝葉斯方法對(duì)新加坡軍團(tuán)菌感染病例數(shù)據(jù)進(jìn)行研究.Liu 等[3]利用廣義最大期望（EM）算法對(duì)0–1 膨脹泊松分布回歸模型中的參數(shù)進(jìn)行估計(jì)，對(duì)美國(guó)底特律城市交通事故死亡數(shù)據(jù)進(jìn)行擬合.夏麗麗等[4]使用局部多項(xiàng)式核回歸法對(duì)0–1 膨脹泊松分布模型進(jìn)行參數(shù)估計(jì)，通過(guò)對(duì)北京市糖尿病患者數(shù)據(jù)的分析，驗(yàn)證了局部多項(xiàng)式核回歸方法的有效性.

對(duì)于0–1 膨脹泊松分布模型，當(dāng)數(shù)據(jù)存在較大變異時(shí)，即樣本均值與樣本方差不相等時(shí)，如果仍然用模型進(jìn)行擬合，效果往往不好.而0–1 膨脹幾何分布模型，不僅可以用于處理樣本數(shù)據(jù)的變異，也適應(yīng)于樣本尾部數(shù)據(jù)退化較慢的情形.對(duì)于0–1 膨脹幾何分布及其回歸模型，目前研究文獻(xiàn)較少，肖翔[5]利用貝葉斯方法對(duì)0–1 膨脹幾何分布回歸模型進(jìn)行參數(shù)估計(jì)，Xiao 等[6]基于Polya-Gamma 潛變量設(shè)計(jì)0–1膨脹幾何分布回歸模型中后驗(yàn)樣本的抽樣機(jī)制.本研究對(duì)0–1 膨脹幾何分布模型進(jìn)行參數(shù)變換，計(jì)算出客觀貝葉斯先驗(yàn)，以期得到更好的擬合效果.

1 0–1 膨脹幾何分布模型

本研究提出0–1 膨脹幾何分布(簡(jiǎn)稱為ZOIGE）模型，即一個(gè)非負(fù)的0-1 膨脹幾何分布的隨機(jī)變量Y，可以表示為Y=V(1?B1)+B1(1?B2).其 中，B1、B2、V相互獨(dú)立，B1為一個(gè)試驗(yàn)成功概率為p1的伯努利隨機(jī)變量；B2為一個(gè)試驗(yàn)成功概率為p2的伯努利隨機(jī)變量；V為一個(gè)服從于試驗(yàn)成功概率為 θ的幾何分布隨機(jī)變量，即P(V=k)=θk(1?θ)，k=0,1,···.隨機(jī)變量Y的分布律為

式中：0 ≤p1≤1，0 ≤p2≤1，0 ≤θ ≤1.可以看出，0–1 膨脹幾何分布是由伯努利分布與幾何分布按照比例p1和1?p1組成的混合分布.當(dāng)p2=1時(shí)，ZOIGE變成零膨脹幾何分布(ZIGE)[7?8]，當(dāng)p1=0時(shí)，ZOIGE 退化成幾何分布.

進(jìn)行參數(shù)變換，令

可得

通過(guò)上述重參數(shù)化，式(1)變?yōu)?/p>

式中：q1≥0，q2≥0，q1+q2≤1，0 ≤θ ≤1.

設(shè)Y=(Y1,Y2,···,Yn)為取自0–1 膨脹幾何分布的觀測(cè)值，由式(4)得出似然函數(shù)公式為

式中：S0=#{i:Yi=0}為集合{i:Yi=0}中包含元素的個(gè)數(shù)；為集合{i:Yi=1}中包含元素的個(gè)數(shù)；

式(5)兩邊取對(duì)數(shù)，得到對(duì)數(shù)似然函數(shù)為

2 客觀先驗(yàn)

2.1 Fisher 信息矩陣

計(jì)算隨機(jī)變量Y,S0,S1,S的期望為

計(jì)算對(duì)數(shù)似然函數(shù)式(6)的一階偏導(dǎo)數(shù)為

計(jì)算對(duì)數(shù)似然函數(shù)式(6)的二階偏導(dǎo)數(shù)為

進(jìn)一步計(jì)算二階偏導(dǎo)數(shù)期望的相反數(shù)，它們是Fisher 信息陣的組成元素.表達(dá)式為

因此，(q1,q2,θ)的Fisher 信息陣為

2.2 Jeffreys 先驗(yàn)

與Laplace 先驗(yàn)比較，Jeffreys 先驗(yàn)?zāi)軌蛟趨?shù)變換下保持不變性，比Laplace 先驗(yàn)具有更廣泛的應(yīng)用場(chǎng)合[9].參數(shù)(q1,q2,θ)的Jeffreys 先驗(yàn)與Fisher信息矩陣行列式的平方根成正比，通過(guò)式(7)可以計(jì)算(q1,q2,θ)的Jeffreys 先驗(yàn)，公式為

2.3 Reference 先驗(yàn)

對(duì)于參數(shù)組合{(∑q1,q2),θ}，(q1,q2)為感興趣的參數(shù)，F(xiàn)isher 信息矩陣(q1,q2,θ)可寫成

其中

根據(jù)文獻(xiàn)[10]，reference 先驗(yàn)求解過(guò)程中，先求出h1和h2，公式為

再完成以下4 個(gè)步驟.

步驟1選取參數(shù)空間的一組緊子集為?i=?12×?3i={(q1,q2)|0

步驟2當(dāng)(q1,q2)給定時(shí)，θ的條件先驗(yàn)為

步驟3結(jié)合式(8)，(q1,q2)的邊緣先驗(yàn)為

步驟4Φ的reference 先驗(yàn)為

對(duì)于參數(shù)組∑合{θ,(q1,q2)}，θ為感興趣的參數(shù)，F(xiàn)isher 信息矩陣(q1,q2,θ)可寫成

步驟1選取與{(q1,q2),θ}參數(shù)空間中相同的一組緊子集.

步驟2當(dāng) θ給定時(shí)，結(jié)合式(8)，(q1,q2)的條件先驗(yàn)為

步驟3結(jié)合式(8)，θ的邊緣先驗(yàn)為

步驟4Φ的reference 先驗(yàn)為

3 后驗(yàn)分析

基于先驗(yàn)分布 πJ，πR1和πR2，分別得到它們的后驗(yàn)分布，通過(guò)R 軟件進(jìn)行抽樣，獲取后驗(yàn)樣本.以πR1為例，(q1,q2,θ)的后驗(yàn)分布為

式中：Y=(Y1,Y2,···,Yn)為觀測(cè)數(shù)據(jù).式(9)的具體形式為

從式(10)可以看出，(q1,q2)的后驗(yàn)邊緣分布為

4 數(shù)值模擬

本節(jié)基于Jeffreys 先驗(yàn)和reference 先驗(yàn)，通過(guò)數(shù)值模擬對(duì)ZOIGE 模型的參數(shù)進(jìn)行估計(jì).樣本容量分別設(shè)為n=20和n=50，θ值設(shè)為 0.8，q1的值分別設(shè)為 0.3和0.7，q2值分別設(shè)為 0.4和 0.6，所有模擬重復(fù)2 000 次，計(jì)算出參數(shù)估計(jì)量的均值和均方誤差見(jiàn)表1 和表2.從表中可以看出，隨著樣本容量的增加，3 種客觀貝葉斯先驗(yàn)下的估計(jì)值越來(lái)越接近真值，均方誤差也越來(lái)越小.對(duì)于q1和q2的估計(jì)，πR1、πR2比 πJ表現(xiàn)更好，這是因?yàn)樵讦蠷1、πR2中包含q1和q2的信息更加豐富.對(duì)于 θ的估計(jì)，πR2比πR1、πJ表現(xiàn)更好，這是因?yàn)樵讦蠷2中 θ為感興趣參數(shù)，集中了更多的樣本信息.

表1 θ=0.8下參數(shù)估計(jì)量的均值Table 1 Mean of parameter estimators whenθ=0.8

表2 θ=0.8下參數(shù)估計(jì)量的均方誤差Table 2 Mean squared error of parameter estimators whenθ=0.8

5 結(jié)語(yǔ)

本研究對(duì)0–1 膨脹幾何分布模型進(jìn)行客觀貝葉斯分析，巧妙地進(jìn)行重參數(shù)化，寫出具有分塊對(duì)角形式的Fisher 信息矩陣.因而，較容易推導(dǎo)出參數(shù)的Jeffreys 先驗(yàn)和reference 先驗(yàn)，這種方法和技巧可以推廣到其他形式的0–1 膨脹分布模型中去.