Tweedie復(fù)合泊松回歸模型的Bayes估計(jì)和影響分析

段星德,張 實(shí),羅露璐,張文專

(貴州財(cái)經(jīng)大學(xué)數(shù)統(tǒng)學(xué)院,貴州貴陽(yáng)550025)

§1 引 言

Tweedie復(fù)合泊松分布常用來(lái)分析由零和正的連續(xù)數(shù)據(jù)所構(gòu)成的半連續(xù)型數(shù)據(jù),這類數(shù)據(jù)廣泛存在于精算科學(xué)(比如:保險(xiǎn)公司對(duì)不同個(gè)體的賠償金額),環(huán)境科學(xué)(比如:一段時(shí)間某個(gè)地區(qū)的降雨量)等研究領(lǐng)域中.對(duì)半連續(xù)數(shù)據(jù)建模常見(jiàn)的有兩類方法,第一,兩部分模型,即混合Bernoulli分布和Gamma分布或者Lognormal分布,常被用來(lái)分別處理零數(shù)據(jù)和正的連續(xù)數(shù)據(jù)兩部分(Fernandes等[1],Piantadosi等[2]),然而這種分割處理方法常會(huì)破壞半連續(xù)數(shù)據(jù)的整體屬性(Hasan等[3]).第二,Hasan等[3]以及Hasan和Dunn[4-6]發(fā)展了Tweedie復(fù)合泊松模型來(lái)對(duì)降雨量和降雨發(fā)生率聯(lián)合建模.進(jìn)一步,針對(duì)Tweedie復(fù)合泊松模型的統(tǒng)計(jì)推斷問(wèn)題,Smyth和J?rgensen[7]研究了Tweedie復(fù)合泊松模型的均值和散度參數(shù)的聯(lián)合建模問(wèn)題以及約束極大似然估計(jì)方法;Dunn和Smyth[8-10]提出了不同的數(shù)值方法去逼近Tweedie復(fù)合泊松分布的密度函數(shù)并給出了R軟件包;Peters等[11]給出了Tweedie復(fù)合泊松模型的極大似然估計(jì)和基于Markov Chain Monte Carlo(簡(jiǎn)稱MCMC)模擬技術(shù)的Bayes估計(jì);Zhang[12]基于Monte Carlo EM和MCMC等算法研究了Tweedie復(fù)合泊松混合效應(yīng)模型的估計(jì)問(wèn)題;Qian等[13]基于分組彈性網(wǎng)估計(jì)技術(shù)研究Tweedie復(fù)合泊松模型的變量選擇問(wèn)題;Yang等[14]把梯度提升樹算法應(yīng)用到Tweedie復(fù)合泊松模型中并用來(lái)預(yù)測(cè)半連續(xù)保險(xiǎn)數(shù)據(jù);Bonat等[15]提出兩類擬似然和偽似然方法對(duì)Tweedie復(fù)合泊松模型進(jìn)行統(tǒng)計(jì)推斷.但據(jù)作者所知,基于Bayes數(shù)據(jù)刪除方法對(duì)Tweedie復(fù)合泊松模型的統(tǒng)計(jì)診斷分析還沒(méi)得到研究.

近幾年來(lái),基于Bayes數(shù)據(jù)刪除影響診斷以引起大量研究者的關(guān)注.特別地,Cho等[16]基于K-L差異統(tǒng)計(jì)量的Bayes數(shù)據(jù)刪除影響測(cè)度來(lái)研究生存模型的統(tǒng)計(jì)診斷問(wèn)題;Zhu等[17]系統(tǒng)研究了三種Bayes數(shù)據(jù)刪除影響測(cè)度的統(tǒng)計(jì)性質(zhì),并把這些統(tǒng)計(jì)診斷方法應(yīng)用于復(fù)雜統(tǒng)計(jì)模型中;Jackson等[18]提出兩類基于MCMC算法的方法去逼近Bayes數(shù)據(jù)刪除影響測(cè)度;Duan和Tang[19]對(duì)廣義部分線性混合效應(yīng)模型建立起一套基于Bayes數(shù)據(jù)刪除影響診斷方法來(lái)評(píng)價(jià)模型對(duì)于刪除一個(gè)數(shù)據(jù)點(diǎn)或數(shù)據(jù)組的敏感性.因此,本文將在Bayes估計(jì)的基礎(chǔ)上基于Bayes數(shù)據(jù)刪除影響方法來(lái)研究Tweedie復(fù)合泊松模型的統(tǒng)計(jì)診斷問(wèn)題.

§2 模型

本節(jié)首先介紹泊松-伽瑪復(fù)合分布,即Tweedie復(fù)合泊松分布,其次引入它們所對(duì)應(yīng)的廣義線性模型:Tweedie復(fù)合泊松回歸模型.

2.1 Tweedie復(fù)合泊松分布

如果一個(gè)隨機(jī)變量Y服從泊松-伽瑪復(fù)合分布,則它可以表示為如下的隨機(jī)表達(dá)式.

其中N為服從均值參數(shù)為λ的泊松分布的隨機(jī)變量;給定N的條件下,Xi(1≤i≤N)獨(dú)立同分布,服從均值為αγ,方差為αγ2的伽瑪分布;且N和Xi相互獨(dú)立.另外當(dāng)N=0時(shí),有Y=0.進(jìn)一步,P(Y=0)=P(N=0)=exp(?λ).

根據(jù)J?rgensen[20],Tweedie復(fù)合泊松分布是指數(shù)散度分布族的一類特殊分布,而指數(shù)散度分布族的概率密度函數(shù)可表示為

其中a(·)和κ(·)是已知函數(shù);θ為定義在R=(?∞,∞)上的參數(shù),φ>0為散度參數(shù).根據(jù)指數(shù)散度分布族的性質(zhì)有:

其中κ0(θ)和κ00(θ)是分別關(guān)于θ的一階和二階導(dǎo)數(shù).根據(jù)(2.1)和(2.2),可得兩組參數(shù)之間的關(guān)系:

顯然參數(shù)p的取值范圍為(1,2).因此把一個(gè)隨機(jī)變量Y服從泊松-伽瑪復(fù)合分布記為Y～Twp(μ,φ),其中E(Y)=μ,var(Y)=φμp,p∈(1,2)稱之為冪指標(biāo)參數(shù).

為了對(duì)模型進(jìn)行統(tǒng)計(jì)推斷,需要推導(dǎo)出泊松-伽瑪復(fù)合分布的概率密度函數(shù),即寫出(2.2)式中所對(duì)應(yīng)的具體表達(dá)式.盡管泊松-伽瑪復(fù)合分布的概率密度函數(shù)沒(méi)有顯示表達(dá),但涉及到的模型可以通過(guò)Y和N的聯(lián)合概率密度函數(shù)進(jìn)行統(tǒng)計(jì)推斷,其中N可視為潛變量.根據(jù)(2.1),Y和N的聯(lián)合概率密度函數(shù)可以表示為

2.2 Tweedie復(fù)合泊松回歸模型

設(shè)yi,i=1,2,···,m是響應(yīng)變量Yi的第i次觀測(cè)值,Y1,···,Ym相互獨(dú)立,則Tweedie復(fù)合泊松回歸模型表示為

其中xi=(xi1,xi2,···,xiq)T∈Rq是q維固定效應(yīng)協(xié)變量,β為對(duì)應(yīng)的q維待估的未知參數(shù).另外

§3 模型的Bayes分析

令

在統(tǒng)計(jì)推斷過(guò)程中,可把潛變量N視為缺失數(shù)據(jù),而Y仍然為觀測(cè)數(shù)據(jù),(Y,N)視為完全數(shù)據(jù).為了獲得模型參數(shù)的Bayes估計(jì),考慮參數(shù)θ的基于完全數(shù)據(jù)(Y,N)的后驗(yàn)分布.類似于Zhang[10],Y和N的聯(lián)合對(duì)數(shù)似然函數(shù)為

因此參數(shù)θ基于完全數(shù)據(jù)(Y,N)的后驗(yàn)分布為

其中p(θ)為參數(shù)θ的先驗(yàn)概率密度函數(shù).進(jìn)一步,根據(jù)Zhang[10],θ的先驗(yàn)分布可假定為

其中β0,Σ0,φ0為事先給定的超參數(shù),U(·)表示均勻分布.

3.1 條件分布和Gibbs抽樣

為了進(jìn)行Bayes推斷,用Gibbs抽樣技術(shù)從聯(lián)合后驗(yàn)分布p(θ,N|Y,X)中抽取所需的隨機(jī)樣本,并基于此隨機(jī)樣本對(duì)參數(shù)θ進(jìn)行Bayes推斷.

由(3.1)-(3.3)式可知,給定θ,Y,X下N的對(duì)數(shù)條件分布可以表示為

于是有

其中I{·}表示示性函數(shù).由于{yi=0}與{ni=0}是兩個(gè)等價(jià)事件組,即當(dāng)yi=0時(shí)隱含著ni=0,因此基于條件分布(3.4)式進(jìn)行抽樣時(shí)只需對(duì)ni>0的情形進(jìn)行抽樣.

由(3.1)-(3.3)式可知,給定N,Y,X,β,φ下β的對(duì)數(shù)條件分布可以表示為

利用Gibbs抽樣技術(shù)對(duì)上述條件分布進(jìn)行抽樣時(shí)涉及到的條件分布都是一些非標(biāo)準(zhǔn)分布且非常復(fù)雜,因此將選擇一些合適的,有針對(duì)性的MH算法來(lái)進(jìn)行隨機(jī)抽樣.

3.2 MH算法

為了實(shí)施MH算法,假設(shè)潛變量ni和參數(shù)θ的當(dāng)前迭代值為和θ(t),新的隨機(jī)樣本和θ?服從如下建議分布:

3.3 Bayes估計(jì)

基于上述Gibbs抽樣和MH算法的混合算法獲得隨機(jī)樣本,就可以得到θ=(β,φ,p)的聯(lián)合Bayes估計(jì)及其對(duì)應(yīng)的標(biāo)準(zhǔn)差.

設(shè)

是來(lái)自于聯(lián)合后驗(yàn)分布p(θ,N|Y,X)的隨機(jī)樣本,則θ=(β,φ,p)的Bayes相合估計(jì)可以表示為

類似地,參數(shù)θ的后驗(yàn)協(xié)方差矩陣var(θ|Y,X)的相合估計(jì)可以通過(guò)它們的隨機(jī)樣本的樣本協(xié)方差矩陣得到,即

另外可以通過(guò)隨機(jī)樣本序列的樣本協(xié)方差矩陣的對(duì)角線元素來(lái)估計(jì)對(duì)應(yīng)的標(biāo)準(zhǔn)差.

§4 Bayes數(shù)據(jù)刪除影響測(cè)度

在貝葉斯框架下,將基于Tweedie復(fù)合泊松回歸模型引入兩類Bayes數(shù)據(jù)刪除影響測(cè)度.首先引入記號(hào),用{yi,xi}表示第i個(gè)數(shù)據(jù)點(diǎn),D={Y,X}為完全數(shù)據(jù)集,表示刪除第i個(gè)數(shù)據(jù)點(diǎn)后剩余的數(shù)據(jù)集;令參數(shù)向量θ=(β,φ,p).

為了度量數(shù)據(jù)點(diǎn){yi,xi}對(duì)參數(shù)θ的聯(lián)合后驗(yàn)分布的影響程度,根據(jù)Cho等[16]和Zhu等[17]的工作,首先引入如下的第一類型Bayes數(shù)據(jù)刪除影響測(cè)度:

其中p(θ|D)和p(θ|D[i])是θ關(guān)于數(shù)據(jù)集D和D[i]的后驗(yàn)分布.該影響測(cè)度Dφ(i)用來(lái)度量數(shù)據(jù)點(diǎn){yi,xi}刪除前后的影響程度,若它的估計(jì)值比較大,則對(duì)應(yīng)的數(shù)據(jù)點(diǎn){yi,xi}可能是強(qiáng)影響點(diǎn)或異常點(diǎn).

另外為了度量數(shù)據(jù)點(diǎn){yi,xi}對(duì)θ的后驗(yàn)均值的影響程度,根據(jù)Zhu等[17]的工作,定義如下的第二類型Bayes數(shù)據(jù)刪除影響測(cè)度.

在實(shí)際計(jì)算中,通常使用參數(shù)θ的樣本后驗(yàn)協(xié)方差矩陣的逆矩陣來(lái)估計(jì)正定陣Wθ.同理若CM(i)的估計(jì)值比較大,那么所對(duì)應(yīng)的數(shù)據(jù)點(diǎn){yi,xi}可能對(duì)參數(shù)θ的后驗(yàn)均值的估計(jì)產(chǎn)生較大的影響;因此數(shù)據(jù)點(diǎn){yi,xi}可能被診斷為強(qiáng)影響點(diǎn)或異常點(diǎn).

在計(jì)算這兩類影響測(cè)度Dφ(i)和CM(i)的過(guò)程中,發(fā)現(xiàn)需要計(jì)算所涉及到的后驗(yàn)分布p(θ|D)和p(θ|D[i]),但所考慮的模型涉及到潛變量N的處理,這給計(jì)算帶來(lái)一定的困難.為了有效的處理潛變量N,發(fā)展了以下一些簡(jiǎn)單有效的計(jì)算公式.

首先根據(jù)Zhu等[17]的計(jì)算方法,令

pi(θ)=p(Y|X,θ)/p(Y[i]|X[i],θ)=p(Yi|Xi,θ).

通過(guò)一些計(jì)算可得到如下的表達(dá)式:

因此第一類型Bayes數(shù)據(jù)刪除影響測(cè)度可以簡(jiǎn)化為

Dφ(i)=logEθ|D{pi(θ)}?1+Eθ|D{logpi(θ)},

這種影響測(cè)度也稱之為K-L差異測(cè)度.進(jìn)一步,可基于上述的MCMC算法從從后驗(yàn)分布p(θ|D)中進(jìn)行抽取樣本來(lái)計(jì)算該影響測(cè)度.

同樣利用上述的MCMC算法從后驗(yàn)分布p(θ|D)中抽取樣本來(lái)計(jì)算后驗(yàn)均值以及θ的后驗(yàn)協(xié)方差矩陣Wθ,從而可計(jì)算出影響測(cè)度CM(i)的值.總之,計(jì)算上述兩類Bayes數(shù)據(jù)刪除影響測(cè)度的關(guān)鍵所在是計(jì)算pi(θ)的值.根據(jù)pi(θ)的定義,可得

這就意味著對(duì)存在潛變量的模型中計(jì)算pi(θ)是繁瑣的.盡管這給計(jì)算影響測(cè)度Dφ(i)和CM(i)帶來(lái)困難,但可利用下面的數(shù)值方法來(lái)計(jì)算pi(θ).具體地,可以利用Monte Carlo方法來(lái)逼近pi(θ),即

其中若yi=0時(shí),所對(duì)應(yīng)的潛變量ni的取值為0;若yi>0時(shí),所對(duì)應(yīng)的潛變量ni的隨機(jī)樣本

從零截?cái)嗖此煞植糵(ni;λ|ni>0)中產(chǎn)生,這里λ是所對(duì)應(yīng)泊松分布的均值參數(shù),它可通過(guò)潛變量N的樣本均值來(lái)進(jìn)行估計(jì).在下一節(jié)的模擬研究和實(shí)例分析中,取J=200或更大的數(shù)值.

§5 模擬研究

在第一個(gè)模擬研究中,假設(shè)響應(yīng)變量{yi,i=1,2,···,m}服從Tweedie復(fù)合泊松分布,即Yi～Twp(μi,φ);模型的系統(tǒng)部分表示為:

log(μi)=β1x1+β2x2+β3x3+ ···+β10x10,

其中樣本量m=150,協(xié)變量xi～N(0,1),i=1,2和xi～N(0,0.25),i=3,4,···,10,且這10個(gè)協(xié)變量是相互獨(dú)立的,回歸參數(shù)向量β的真值為

βT=(β1,β2,β3,β4,β5,β6,β7,β8,β9,β10)=(0.85,0.85,1,1,1,1,1,1,1,1),

散度參數(shù)φ的真值為φ=2,冪指標(biāo)參數(shù)p的真值為p=1.6.

為了研究參數(shù)分布的不同先驗(yàn)信息對(duì)Bayes估計(jì)的影響,將考慮如下三種不同的超參數(shù)設(shè)置.

類型I:參數(shù)β的先驗(yàn)分布的超參數(shù)β0的取值為β的真值,即

β0T=βT=(0.85,0.85,1,···,1),

協(xié)方差陣Σ0=0.25I10,I10表示10階單位陣.這種情形可視為參數(shù)的先驗(yàn)分布具有良好的先驗(yàn)信息.

類型Ⅱ:設(shè)β0T=1.5×βT=1.5×(0.85,0.85,1,···,1),協(xié)方差陣Σ0=0.75I10.這種情形可視為具有不準(zhǔn)確的先驗(yàn)信息.

類型Ⅲ:設(shè)β0T=(0,0,0,···,0),協(xié)方差陣Σ0=100I10.這種情形可視為具有無(wú)先驗(yàn)信息.

基于上述設(shè)計(jì)產(chǎn)生的100個(gè)模擬數(shù)據(jù),用前面提出的Bayes混合算法進(jìn)行100次重復(fù)試驗(yàn)得到參數(shù)的聯(lián)合Bayes估計(jì).這里用計(jì)算出的所有參數(shù)的EPSR(estimated potential scale reduction)值來(lái)判斷混合算法中所產(chǎn)生的Markov鏈的收斂性.在每一次模擬測(cè)試中,利用每一個(gè)參數(shù)的3組不同初值分別產(chǎn)生3條平行的Markov鏈,并基于每一次迭代值計(jì)算出所有參數(shù)的EPSR值.在所有模擬測(cè)試中,發(fā)現(xiàn)所有參數(shù)的EPSR值在迭代次數(shù)為1000次時(shí)均小于1.2,說(shuō)明模擬研究中的混合算法是收斂的.在本次模擬研究中,利用隨機(jī)抽取一次模擬測(cè)試中12個(gè)參數(shù)的迭代值進(jìn)行計(jì)算相應(yīng)的EPSR值,其計(jì)算結(jié)果見(jiàn)圖(1左圖).因此在每次重復(fù)實(shí)驗(yàn)中迭代10000次,為了保證算法的收斂性,丟掉了所有參數(shù)的前5000次迭代值,保留后5000次迭代值用來(lái)計(jì)算參數(shù)的Bayes估計(jì).另外在實(shí)施MH算法和隨機(jī)游走M(jìn)H算法中,選擇調(diào)節(jié)參數(shù)

讓它們所對(duì)應(yīng)的參數(shù)抽樣的平均接受率分別為25.7%,29.1%,24.7%,27.7%.表1給出了所有參數(shù)θ的相對(duì)偏差值(記為RB)、標(biāo)準(zhǔn)差和RMS,其中

在第二個(gè)模擬研究中,為了檢驗(yàn)協(xié)變量出現(xiàn)多重共線性情形時(shí)對(duì)模型的參數(shù)估計(jì)的影響情況,在第一個(gè)模擬研究的基礎(chǔ)上,對(duì)協(xié)變量進(jìn)行如下兩種假定:

圖1 模擬研究中所有參數(shù)的EPSR值(左圖);實(shí)例分析中所有參數(shù)的EPSR值(右圖)

表1 隨機(jī)模擬研究中未知參數(shù)的Bayes估計(jì)

(1)協(xié)變量

(x1,x2)～N2(μ0,Σ0);

(其中μ0=(0,0),協(xié)方差陣Σ0的元素為0.5|i?j|,i,j=1,2)和xi～N(0,0.25),i=3,4,···,10;

(2)協(xié)變量

(x1,x2,x3)～N3(μ0,Σ0).

(其中μ0=(0,0,0),協(xié)方差陣Σ0的元素為0.5|i?j|,i,j=1,2,3)和xi～N(0,0.25),i=4,5,···,10.其它假定均和第一個(gè)模擬研究中的一致.

同理用上述設(shè)置產(chǎn)生的兩組模擬數(shù)據(jù)集在模擬一中考慮的三種先驗(yàn)情形下分別進(jìn)行模型參數(shù)估計(jì),為了節(jié)約空間,只列出無(wú)先驗(yàn)信息情形下參數(shù)的RB計(jì)算結(jié)果.具體為:

假定(1)中回歸參數(shù)向量β的RB計(jì)算值分別為:

0.0008,?0.0087,0.0182,?0.0267,?0.0355,0.0186,?0.0216,?0.0078,?0.0418,?0.0434,

冪指標(biāo)參數(shù)p的RB計(jì)算值為:0.0080,散度參數(shù)φ的RB計(jì)算值為?0.1445;

假定(2)中β的RB計(jì)算值為:

?0.0404,0.0117,?0.0231,?0.0221,?0.0167,?0.0016,0.0374,?0.0137,?0.0177,?0.0186,

冪指標(biāo)參數(shù)p的RB計(jì)算值為:?0.0157,散度參數(shù)φ的RB計(jì)算值為?0.4314.

從上述計(jì)算結(jié)果發(fā)現(xiàn),在模擬設(shè)置中,隨著協(xié)變量多重共線性的嚴(yán)重程度增加時(shí),模型的散度參數(shù)φ的估計(jì)偏差增大,而對(duì)參數(shù)向量β和冪指標(biāo)參數(shù)p的估計(jì)影響較小.

在第三個(gè)模擬研究中,假定模型的系統(tǒng)部分表為:

log(μi)=β0+β1xi,

其中樣本量m=150,協(xié)變量xi來(lái)自于均勻分布U(?1,1),回歸參數(shù)向量β、散度參數(shù)φ、指標(biāo)參數(shù)p的真值分別為:βT=(β0,β1)=(0.5,0.5)、φ=2、p=1.6.基于上述模型在模擬數(shù)據(jù)集中通過(guò)如下兩種情形分別產(chǎn)生異常點(diǎn).

情形一,把協(xié)變量重置為:

x1=x1+8,x75=x75+8,x150=x150+7;

情形二,把響應(yīng)變量重置為

yi=yi+10,i=1,75,150.

并基于上述的混合算法產(chǎn)生的參數(shù)迭代值計(jì)算兩類Bayes數(shù)據(jù)刪除影響測(cè)度Dφ(i)和CM(i),只列出情形二下這兩類影響測(cè)度的散點(diǎn)圖(見(jiàn)圖2).從圖2可看出,正如所預(yù)期的一樣,數(shù)據(jù)點(diǎn)1,75,150被識(shí)別為異常點(diǎn).

圖2 模擬研究中兩類Bayes數(shù)據(jù)刪除影響測(cè)度Dφ(i)和CM(i)的散點(diǎn)圖

§6 實(shí)例分析

在這個(gè)部分里,分析的一組包含240份汽車保單的索賠數(shù)據(jù)來(lái)自Ismail[21],該數(shù)據(jù)集包括車年數(shù),索賠次數(shù),索賠強(qiáng)度,其中索賠強(qiáng)度等于平均每次索賠的賠款金額.模型的響應(yīng)變量是累計(jì)損失,表示在保險(xiǎn)期間累計(jì)的索賠金額,其數(shù)值等于索賠強(qiáng)度與索賠次數(shù)的乘積.另外,數(shù)據(jù)集包括5個(gè)因子解釋變量及其對(duì)應(yīng)的水平:保單類型(綜合險(xiǎn),非綜合險(xiǎn)),汽車產(chǎn)地(國(guó)產(chǎn),進(jìn)口),汽車用途和駕駛?cè)诵詣e(商業(yè),私人-女性,私人-男性),車齡和行駛區(qū)域(0-1年,2-3年,4-5年,6年以上),行駛區(qū)域(中部,東部,東馬來(lái)西亞,北部,南部).在數(shù)據(jù)預(yù)處理過(guò)程中,發(fā)現(xiàn)有7份保單中車年數(shù)均為0,因此只分析其它233份保單的數(shù)據(jù).對(duì)響應(yīng)變量做變換,令y=y/100000,變換后的數(shù)據(jù)樣本均值為1.7029,樣本方差為19.154,零數(shù)據(jù)的比例為56.7%.這里,用上述發(fā)展的Tweedie回歸模型對(duì)該數(shù)據(jù)集建立模型,其模型系統(tǒng)部分表示為:

log(μi)=log(di)+β0+β1xi1+β2xi2+ ···+β11xi,11,

其中車年數(shù)的對(duì)數(shù)log(di)是抵消項(xiàng),虛擬變量x1,x2,x3,x4,x5,x6,x7,x8,x9,x10和x11分別表示非綜合險(xiǎn)水平,進(jìn)口水平,私人-女性水平,私人-男性水平,2-3年水平,4-5年水平,6年以上水平,東部水平,東馬來(lái)西亞水平,北部水平,南部水平;而把綜合險(xiǎn)水平,國(guó)產(chǎn)水平,商業(yè)水平,0-1年水平,中部水平分別作為各個(gè)分類變量的基準(zhǔn)水平.在利用前面提出的Bayes混合算法對(duì)參數(shù)進(jìn)行估計(jì)時(shí),由于不知道未知參數(shù)的真實(shí)值,因此采用回歸參數(shù)的無(wú)先驗(yàn)信息,即令

β0T=(0,···,0)12×1,

協(xié)方差陣

Σ0=100I12.

為了檢驗(yàn)算法的收斂性,用模擬研究中所引進(jìn)的方法來(lái)計(jì)算14個(gè)參數(shù)的EPSR值,其計(jì)算結(jié)果見(jiàn)圖1(右圖).從該圖可以看出,當(dāng)?shù)螖?shù)約為500次時(shí),所對(duì)應(yīng)的14個(gè)參數(shù)的EPSR值均小于1.2,說(shuō)明所使用的混合算法在500次左右時(shí)就已收斂.因此,在進(jìn)行Bayes估計(jì)時(shí),丟掉前面5000次迭代值而保留后面5000次迭代值用來(lái)進(jìn)行計(jì)算.通過(guò)計(jì)算得到散度參數(shù)φ的估計(jì)值(標(biāo)準(zhǔn)差)為0.6381(0.0431),冪指標(biāo)參數(shù)的估計(jì)值(標(biāo)準(zhǔn)差)為1.4193(0.0217);另外模型回歸參數(shù)的估計(jì)值(標(biāo)準(zhǔn)差)分別為-12.3324(1.1789),-0.5805(0.1947),0.0092(0.1294),6.6860(1.1684),7.3796(1.1651),-0.5943(0.1639),-0.9383(0.1685),-0.9665(0.1590),-0.6654(0.2168),-0.5250(0.2115),-0.2655(0.1490),-0.0662(0.1650).

類似于模擬研究,用上述的混合算法產(chǎn)生一系列隨機(jī)樣本,丟掉所有參數(shù)的前5000次迭代值并利用后5000次迭代值來(lái)計(jì)算Bayes數(shù)據(jù)刪除影響測(cè)度Dφ(i)和CM(i)的值,然后作出它們的散點(diǎn)圖,結(jié)果見(jiàn)圖3.從圖3可看出:數(shù)據(jù)點(diǎn)1和6被Dφ(i)和CM(i)測(cè)度同時(shí)診斷為異常點(diǎn).

§7 結(jié) 論

本文給出了基于Gibbs抽樣和MH算法的混合算法的Tweedie復(fù)合泊松回歸模型的Bayes聯(lián)合估計(jì),并在估計(jì)的基礎(chǔ)上發(fā)展了兩類Bayes數(shù)據(jù)刪除影響測(cè)度;最后通過(guò)隨機(jī)模擬研究和一個(gè)實(shí)例分析說(shuō)明了提出的方法是可行和有效的.

圖3 實(shí)例分析中兩類Bayes數(shù)據(jù)刪除影響測(cè)度Dφ(i)和CM(i)的散點(diǎn)圖