基于貝塔分布的最優(yōu)置信區(qū)間研究

溫利民,張良超,章 溢,劉 蔚

(1.江西師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,江西 南昌 330022;2.江西師范大學(xué)財政金融學(xué)院,江西 南昌 330022)

0 引言

在概率論與數(shù)理統(tǒng)計中，區(qū)間估計是參數(shù)估計的重要內(nèi)容.設(shè)θ是總體的一個未知參數(shù)，在總體為連續(xù)型分布的情況下，通過樞軸量法可得到參數(shù)的置信區(qū)間.然而當(dāng)樞軸量的分布為單峰非對稱時，利用傳統(tǒng)方法構(gòu)造的區(qū)間是等尾置信區(qū)間，而不是最優(yōu)置信區(qū)間.

關(guān)于最優(yōu)置信區(qū)間的定義，常見的有2種:一種是在給定置信水平的區(qū)間估計下要求平均區(qū)間長度最短，另一種是在給定平均區(qū)間長度下要求置信度盡可能大或精確度盡可能高.本文主要考慮第1種定義，即在給定置信度水平下求解平均值區(qū)間長度最短的區(qū)間估計.在數(shù)理統(tǒng)計中，關(guān)于置信區(qū)間的最優(yōu)性的研究較多.夏樂天等[1]討論了指數(shù)分布參數(shù)的最短區(qū)間估計;袁長迎等[2]在伽瑪分布形狀參數(shù)已知時研究了尺度參數(shù)的最短區(qū)間估計;徐美萍等[3]研究了在威布爾分布中尺度參數(shù)的最短區(qū)間估計;王秀麗[4]研究了均勻分布參數(shù)的最短置信區(qū)間;薛峰等[5]利用粒子群優(yōu)化算法研究了貝塔分布參數(shù)的最短置信區(qū)間.在區(qū)間估計問題中，樞軸量G的分布通常是單峰分布，如正態(tài)分布、t分布、卡方分布、F分布等.由于正態(tài)分布和t分布是單峰對稱分布，且未知參數(shù)θ在樞軸量G的分子上，所以用傳統(tǒng)方法構(gòu)造的區(qū)間就是最短置信區(qū)間.然而由于卡方分布和F分布為單峰非對稱分布，所以傳統(tǒng)方法構(gòu)造的區(qū)間不是最短置信區(qū)間.孫鵬哲等[6]研究了卡方分布的最短置信區(qū)間，得出在各種自由度和常用的置信水平下最優(yōu)左側(cè)尾概率分配統(tǒng)計表.李廣正[7]運用拉格朗日乘數(shù)法和利用Mathematica軟件給出了F分布的最短區(qū)間估計用表.上述文獻(xiàn)對于F分布和卡方分布的最短置信區(qū)間討論較多，而對于貝塔分布的研究較少.

在已有的研究基礎(chǔ)上，本文給出在一類形狀參數(shù)下基于貝塔分布的最短置信區(qū)間，這可以適用于當(dāng)樞軸量服從貝塔分布時求解最短置信區(qū)間.利用具有單峰的貝塔分布的密度函數(shù)，可得到待估參數(shù)的最短置信區(qū)間.此外，在貝葉斯框架下，當(dāng)參數(shù)的后驗分布服從貝塔分布時，可得到參數(shù)的最短后驗區(qū)間估計.通過模擬分析，驗證貝塔分布的最短置信區(qū)間的優(yōu)越性.

1 最短置信區(qū)間的概念界定

1.1 區(qū)間估計的定義及說明

定義1設(shè)總體X具有概率函數(shù)f(x;θ),θ為未知參數(shù).X1,X2,…,Xn是取自總體X的一個樣本，若對于事先給定的α(0<α<1),存在2個統(tǒng)計量

T1=T1(X1,X2,…,Xn)與T2=T2(X1,X2,…,Xn)

使得P(T1≤θ≤T2)=1-α成立，則稱區(qū)間[T1,T2]為參數(shù)θ的置信水平為1-α的置信區(qū)間,其中T1和T2分別被稱為置信水平1-α的置信下限和置信上限.

由定義可以看出,T1和T2都是不依賴于未知參數(shù)的隨機(jī)變量，因此置信區(qū)間[T1,T2]是隨機(jī)區(qū)間.Pθ(T1≤θ≤T2)=1-α表示:對樣本X1,X2,…,Xn觀測多次，得到許多不同的區(qū)間[T1,T2],在這些確定的區(qū)間中，大約有(1-α)×100%的比例包含了未知參數(shù)θ的真值，而約有α×100%的比例不包含其真值.需特別注意，對于一次抽樣所得到的一個區(qū)間，決不能理解為“不等式T1≤θ≤T2成立的概率為1-α”.因為在給定樣本下的T1和T2是2個確定的數(shù)，從而只有2種可能:要么這個區(qū)間包含θ；要么這個區(qū)間不包含θ.因此，定義說明區(qū)間[T1(X1,X2,…,Xn),T2(X1,X2,…,Xn)]屬于包含未知參數(shù)θ的區(qū)間類的置信水平是1-α，這也說明置信水平與概率是有所區(qū)別的，不可混淆.

1.2 區(qū)間估計的可靠度和精確度及其關(guān)系

當(dāng)參數(shù)真值為θ時，自然希望隨機(jī)區(qū)間[T1,T2]包含θ的概率Pθ(T1≤θ≤T2)要大.因此,一個好的區(qū)間估計應(yīng)該對所有屬于參數(shù)空間Θ的θ，概率Pθ(T1≤θ≤T2)都相當(dāng)大.

若一個區(qū)間估計的置信系數(shù)越大，則該區(qū)間估計的可靠度越高.但是，構(gòu)造一個置信系數(shù)很大的區(qū)間估計并不是一件難事.如將明天中午12點的氣溫估計在-10～50 ℃之間，這個估計的可靠度很高，但由于它的范圍太大，很不精確，所以一個好的區(qū)間估計還有一個精確度的要求.

1.3 最短置信區(qū)間的確定

區(qū)間估計的精確度的標(biāo)準(zhǔn)不止一個，常用的標(biāo)準(zhǔn)有2個:

1)區(qū)間[T1(X1,X2,…,Xn),T2(X1,X2,…,Xn)]的平均長度Eθ[T2-T1]要短，即區(qū)間的范圍不能太大，這是符合實際的;

2)設(shè)參數(shù)真值為θ，在θ*≠θ時，自然希望區(qū)間[T1,T2]包含θ*的概率要小，即區(qū)間[T1,T2]包含非真值的情況出現(xiàn)越少越好.

在給定樣本容量n后，可靠度與精確度是相互制約著的.為了提高可靠度，可以通過增大區(qū)間范圍來實現(xiàn)，但是會降低精確度.反過來，為了提高精確度，可通過減小區(qū)間范圍來實現(xiàn)，但是會降低可靠度.為此本文采用J. Neyman建議的某種折中方案:在使得置信系數(shù)達(dá)到一定要求的前提下，尋找精確度盡可能高的區(qū)間估計，也就是要求區(qū)間平均長度盡可能短，或者區(qū)間包含非真值的概率盡可能小，這2個要求可能同時達(dá)到，也可能不同時達(dá)到.

下面介紹構(gòu)造置信區(qū)間的常用方法，即樞軸量法.可按下列3個步驟構(gòu)造η=g(θ)的置信區(qū)間.

1)構(gòu)造樣本(X1,X2,…,Xn)和未知參數(shù)η的一個函數(shù)G=G(X1,X2,…,Xn;η),要求G的分布與未知參數(shù)η無關(guān)，稱具有這種性質(zhì)的函數(shù)為樞軸量.

2)對給定的α(0<α<1),選取2個常數(shù)c和d(c

Pθ(c≤G(X1,X2,…,Xn;η)≤d)=1-α,?θ∈Θ.

3)若不等式c≤G(X1,X2,…,Xn;η)≤d可等價地變換為

T1(X1,X2,…,Xn)≤η≤T2(X1,X2,…,Xn),

則

Pθ(T1(X1,X2,…,Xn)≤η≤T2(X1,X2,…,Xn))=1-α,?θ∈Θ,

從而[T1(X1,X2,…,Xn),T2(X1,X2,…,Xn)]是η的一個置信水平為1-α的置信區(qū)間.當(dāng)G(X1,X2,…,Xn;η)是η的連續(xù)嚴(yán)格單調(diào)函數(shù)時，這2個不等式的等價關(guān)系總是可以做到的.當(dāng)g(θ)=θ時，[T1(X1,X2,…,Xn),T2(X1,X2,…,Xn)]是θ的一個置信水平為1-α的置信區(qū)間.

一般來講，滿足要求的c和d是不唯一的，若有可能，應(yīng)選在平均區(qū)間長度Eθ[T2-T1]達(dá)到最短時的c與d,則此時所求得的置信區(qū)間被稱為置信水平1-α的最短置信區(qū)間.由于區(qū)間平均長度與所構(gòu)造的樞軸量密切相關(guān)，所以接下來考慮在2類樞軸量形式下的最短置信區(qū)間估計.以下討論都是在樞軸量服從單峰分布情況下進(jìn)行的.

1)樞軸量G具有如下形式:

G=T(X1,X2,…,Xn)(θ+U(X1,X2,…,Xn)),

其中T(X1,X2,…,Xn)>0.

由Pθ(c≤G≤d)=1-α(即Pθ(c≤T(θ+U)≤d)=1-α)得到參數(shù)θ的置信區(qū)間為[cT-1-U,dT-1-U],平均區(qū)間長度為(d-c)Eθ(T-1).考慮在平均區(qū)間長度最短下的區(qū)間估計，即為求如下條件極值問題:

(1)

運用拉格朗日乘數(shù)法，令L=d-c+λ(F(d)-F(c)-1+α),對L關(guān)于c、d分別求偏導(dǎo)并令其為0得

(2)

其中F(·)、f(·)分別表示G的分布函數(shù)與密度函數(shù).由式(2)可知,f(c)=f(d).所以條件極值問題(1)可轉(zhuǎn)化為如下所示的2元方程組求解問題:

(3)

當(dāng)f(x)為單峰對稱密度函數(shù)(比如正態(tài)分布、t分布的密度函數(shù))時，由式(3)容易看出，此時等尾置信區(qū)間即為最短置信區(qū)間.當(dāng)f(x)為單峰非對稱密度函數(shù)(如卡方分布和F分布的密度函數(shù))時，只要求解式(3)就可得最短置信區(qū)間，文獻(xiàn)[8]證明了式(3)有唯一解，具體求解可利用求根法或黃金分割法[9]得到.

2)樞軸量G具有如下形式:

G=T(X1,X2,…,Xn)(θ+U(X1,X2,…,Xn))-1,

其中T(X1,X2,…,Xn)≥0.

由Pθ(c≤G≤d)=1-α(即Pθ(c≤T(θ+U)-1≤d)=1-α,其中c、d同號)得到參數(shù)θ的置信區(qū)間為[d-1T-U,c-1T-U],平均區(qū)間長度為(c-1-d-1)Eθ(T).考慮在平均區(qū)間長度最短下的區(qū)間估計，即為求如下條件極值問題:

(4)

同理運用拉格朗日乘數(shù)法，把條件極值問題(4)可轉(zhuǎn)化為如下所示的2元方程組求解問題:

(5)

當(dāng)f(x)為單峰對稱密度函數(shù)(如正態(tài)分布、t分布的密度函數(shù))時，由式(5)容易看出，此時等尾置信區(qū)間對應(yīng)的c、d滿足f(c)=f(d),然而c2f(c)≠d2f(d),從而等尾置信區(qū)間不是最短置信區(qū)間，且文獻(xiàn)[8]證明了式(5)有唯一解.

2 貝塔分布的最短置信區(qū)間

現(xiàn)有文獻(xiàn)對于F分布和卡方分布的最短置信區(qū)間討論較多，而對于貝塔分布的研究較少.事實上，貝塔分布在一類形狀參數(shù)下也是單峰非對稱分布.接下來研究在樞軸量服從貝塔分布時的最短置信區(qū)間，首先給出一個引理.

引理1[10]設(shè)總體X的密度函數(shù)為p(x),分布函數(shù)為F(x),X1,X2,…,Xn為樣本,則樣本極差Z=X(n)-X(1)的分布函數(shù)為

2.1 單峰貝塔分布的最短置信區(qū)間

定理1設(shè)X1,X2,…,Xn是來自均勻分布U(θ1,θ2)的一個樣本，則θ2-θ1的一個無偏估計為

證令Yi=(Xi-θ1)/(θ2-θ1),i=1,2,…,n,則Yi獨立同分布于U(0,1),由引理1可知

Y(n)-Y(1)～Beta(n-1,2).

由貝塔分布的數(shù)學(xué)期望公式有

E(Y(n)-Y(1))=E((X(n)-X(1))/(θ2-θ1))=(n-1)/(n+1),

由定理1的證明過程知，事實上可構(gòu)造樞軸量:

G=(X(n)-X(1))/(θ2-θ1)～Beta(n-1,2).

對給定的置信水平1-α,若c、d滿足

P(c≤G≤d)=1-α,

可得θ2-θ1的置信區(qū)間為

[(X(n)-X(1))d-1,(X(n)-X(1))c-1].

(6)

注意到，要使得式(6)的平均區(qū)間長度最短等價于求解式(5)，其中F(·)為Beta(n-1,2)的分布函數(shù),f(·)為Beta(n-1,2)的密度函數(shù).

2.2 與等尾置信區(qū)間的比較

接下來討論給定置信水平0.90和0.95以及樣本容量n,比較θ2-θ1的最短置信區(qū)間與等尾置信區(qū)間.首先由式(6)知，平均區(qū)間長度正比于c-1-d-1,因此記等尾置信區(qū)間長度為

L1=(Betaα/2(n-1,2))-1-(Beta1-α/2(n-1,2))-1.

其中c、d是在式(6)中當(dāng)平均區(qū)間長度達(dá)到最短時的取值，具體可通過Matlab軟件中的解方程組命令fsolve求解，因此記最短置信區(qū)間長度為

L2=c-1-d-1.

2者的相對差異記為e(n)=(L1-L2)/L2.對不同的樣本容量，得到如表1和表2所示的結(jié)果.

在表1、表2中的c、d是在式(6)中當(dāng)平均區(qū)間長度達(dá)到最短時的取值.從表1、表2可以看出:在給定置信水平情況下，2種方法求得的置信區(qū)間長度都隨著樣本容量n的增加而變短,而且c、d也隨著樣本容量增加而變大.這是由于貝塔分布的形狀參數(shù)n-1增加，其密度函數(shù)越呈現(xiàn)“尖峰左偏”形狀，樣本的集中趨勢越來越明顯.

表1 在置信水平0.90下最短置信區(qū)間與等尾置信區(qū)間的比較

表2 在置信水平0.95下最短置信區(qū)間與等尾置信區(qū)間的比較

通過對比可以看出:在每一個置信水平和樣本容量的組合下，本文計算的最短置信區(qū)間要優(yōu)于等尾置信區(qū)間，隨著樣本容量的增大，2種置信區(qū)間的相對差異逐漸減小，但即使樣本容量取到26，在置信水平0.95的情況下相對差異仍有15.16%.在樣本容量n≤10時，2者的相對差異更大，大多數(shù)達(dá)到20%以上.由于給定樣本容量，在置信水平0.95情況下的2種置信區(qū)間的相對差異比在置信水平0.90的情況下的更小,所以,在小樣本情況下最短置信區(qū)間優(yōu)勢明顯.

事實上，當(dāng)參數(shù)m>1、n>1時，貝塔分布Beta(m,n)都是單峰分布.如m>n>1是左偏單峰分布，n>m>1是右偏單峰分布，n=m>1是單峰對稱分布.當(dāng)樞軸量G具有如下形式:

G=T(X1,X2,…,Xn)(θ+U(X1,X2,…,Xn))-1

且服從貝塔分布時,本文構(gòu)造了在置信水平0.90和0.95下基于Beta(m,n)的最短區(qū)間估計用表，整理在表3～表6中.

表3 在置信水平0.90下最短置信區(qū)間的左側(cè)端點值

表4 在置信水平0.90下最短置信區(qū)間的右側(cè)端點值

表6 在置信水平0.95下最短置信區(qū)間的右側(cè)端點值

從表3～表6可以看出:在給定置信水平情況下，n固定，m越大，最短置信區(qū)間的左側(cè)端點值越大,最短置信區(qū)間的右側(cè)端點值也越大，但增大幅度比左側(cè)端點值更小.在給定置信水平情況下，m固定，n越大，最短置信區(qū)間的左側(cè)端點值越小,最短置信區(qū)間的右側(cè)端點值也越小，但減小幅度比左側(cè)端點值更大.m和n都固定，在置信水平更高情況下的左側(cè)端點值更小以及右側(cè)端點值更大，這導(dǎo)致置信區(qū)間長度增大，這體現(xiàn)了可靠度與精確度其實是相互制約的.

2.3 Beta(a,1)型的最短置信區(qū)間

前面的分析討論都是在基于樞軸量的分布為單峰分布的假設(shè)下得到的.但當(dāng)a>1時，Beta(a,1)分布的密度函數(shù)是單調(diào)遞增的.接下來討論當(dāng)樞軸量服從Beta(a,1)分布時在2類樞軸量形式下的最短置信區(qū)間估計.

1)樞軸量G具有如下形式:

G=T(X1,X2,…,Xn)(θ+U(X1,X2,…,Xn)),

其中T(X1,X2,…,Xn)>0.

由Pθ(c≤G≤d)=1-α(即Pθ(c≤T(θ+U)≤d)=1-α)得到參數(shù)θ的置信區(qū)間為[cT-1-U,dT-1-U],平均區(qū)間長度為(d-c)Eθ(T-1).此外，G的密度函數(shù)為

則尋求c、d使得

(7)

根據(jù)密度函數(shù)單調(diào)遞增的特點，對于給定的置信水平，在平均區(qū)間長度最短下，應(yīng)當(dāng)選取g(x)取值較大的部分，即應(yīng)選取c0,使得

(8)

其中c0=α1/a,得到參數(shù)θ的置信區(qū)間為[c0T-1-U,T-1-U],平均區(qū)間長度為(1-c0)Eθ(T-1).由于密度函數(shù)嚴(yán)格單增，所以觀察式(7)和式(8),顯然有d-c≥1-c0.綜上所述，此時的最短置信區(qū)間為[c0T-1-U,T-1-U].

2)樞軸量G具有如下形式:

G=T(X1,X2,…,Xn)(θ+U(X1,X2,…,Xn))-1,

其中T(X1,X2,…,Xn)≥0.

若Pθ(c≤G≤d)=1-α(即Pθ(c≤T(θ+U)-1≤d)=1-α)得到參數(shù)θ的置信區(qū)間為[d-1T-U,c-1T-U],平均區(qū)間長度為(c-1-d-1)Eθ(T).令

所以(d-c)/(cd)≥(1-c0)/c0.

2.4 最短后驗置信區(qū)間

對于區(qū)間估計問題，在上述討論中把參數(shù)看成一個常數(shù)，在求置信區(qū)間時要構(gòu)造一個樞軸量，這一點技巧性較強(qiáng)，有時是比較困難的.并且在理解置信水平和置信區(qū)間時也會產(chǎn)生困難，而貝葉斯方法具有處理方便和含義清晰的優(yōu)點.貝葉斯統(tǒng)計方法是英國統(tǒng)計學(xué)家托馬斯·貝葉斯(Thomas Bayes)提出的一種方法,其主要的核心思想是將未知參數(shù)看成隨機(jī)變量.這使得統(tǒng)計學(xué)的區(qū)間估計得到了更好的解釋.貝葉斯統(tǒng)計方法已成為現(xiàn)代統(tǒng)計學(xué)不可或缺的重要內(nèi)容,在數(shù)理統(tǒng)計、生物統(tǒng)計、醫(yī)學(xué)統(tǒng)計、環(huán)境統(tǒng)計、金融統(tǒng)計與精算等領(lǐng)域[11-15]中都有廣泛的應(yīng)用.

若參數(shù)θ的先驗分布為π(θ),樣本分布函數(shù)為F(x;θ),由貝葉斯定理可得參數(shù)θ的后驗分布π*(θ|x).若給定概率1-α,找到一個區(qū)間[c,d],使得P(c≤θ≤d|x)=1-α成立，這樣求得的區(qū)間就是參數(shù)θ的貝葉斯置信區(qū)間，稱1-α為置信水平.注意到，在貝葉斯統(tǒng)計中，把參數(shù)θ看成是隨機(jī)變量，直接從后驗分布中推導(dǎo)得出置信區(qū)間，并且把置信水平1-α很自然地解釋為參數(shù)落入這一區(qū)間的概率.

置信水平和平均區(qū)間長度是評價貝葉斯區(qū)間估計的2個標(biāo)準(zhǔn)，在置信水平給定的情況下，希望平均區(qū)間長度越短越好.接下來，考慮在參數(shù)θ的后驗分布為貝塔分布時的最短后驗區(qū)間估計.

設(shè)X1,X2,…,Xn是來自負(fù)二項分布NB(m,θ)的樣本，其分布函數(shù)為

給定置信水平1-α，若P(c≤θ≤d|x)=F(d)-F(c)=1-α,則得到θ的貝葉斯置信區(qū)間為[c,d],區(qū)間長度為d-c.欲求最短置信區(qū)間，即為求如下條件極值問題:

(9)

運用拉格朗日乘數(shù)法，把條件極值問題(9)可轉(zhuǎn)化為如下所示的2元方程組求解問題:

其中F(·)、f(·)分別表示Beta(a*,b*)的分布函數(shù)與密度函數(shù).

3 數(shù)值例子

假設(shè)某企業(yè)生產(chǎn)的圓盤直徑服從均勻分布U(c,d),實際統(tǒng)計14個該種圓盤直徑，數(shù)據(jù)如下:8.022 2,7.965 0,8.016 0,8.001 9,8.047 3,8.014 9,8.030 0,7.995 4,7.993 2,8.032 5,7.958 3,7.963 3,7.967 3,7.989 1.

由定理1得到d-c的估計值為0.102 7.當(dāng)置信水平為0.90時，經(jīng)過簡單計算，等尾置信區(qū)間為[0.091 4,0.126 5];查表1可知，最短置信區(qū)間為[0.089 4,0.119 0].同理計算,當(dāng)置信水平為0.95時，等尾置信區(qū)間為[0.090 6,0.134 6],最短置信區(qū)間為[0.089 2,0.126 7].以上結(jié)果表明：當(dāng)置信水平為0.90和0.95時，本文方法得到的置信區(qū)間長度比等尾置信區(qū)間長度更短，且它們的長度之比分別為0.843 3和0.852 3，由此可見這種誤差是不可忽略的.

4 總結(jié)

通過上述分析可以看出:研究在樞軸量服從貝塔分布時的最短置信區(qū)間是十分有必要且有意義的，尤其是當(dāng)樣本容量較小時.本文介紹的方法在理論上可以計算在更廣泛的參數(shù)組合下的貝塔分布的最短置信區(qū)間.