二項(xiàng)分布下基于鞍點(diǎn)逼近的總體成數(shù)置信區(qū)間的構(gòu)造

夏麗麗，田茂再，朱 鈺

(1.北京工業(yè)大學(xué) 應(yīng)用數(shù)理學(xué)院，北京 100872；2.中國(guó)人民大學(xué) a.應(yīng)用統(tǒng)計(jì)科學(xué)研究中心；b.統(tǒng)計(jì)學(xué)院，北京 100872；3.西安財(cái)經(jīng)大學(xué) 統(tǒng)計(jì)學(xué)院，陜西 西安 710100)

一、引 言

當(dāng)今社會(huì)發(fā)展迅速，人們的生活也越來(lái)越富裕，衛(wèi)生保健意識(shí)也在不斷增強(qiáng)，因此就醫(yī)的人數(shù)就大大增加，這也就增加了醫(yī)療事故的發(fā)生率。雖然政府已經(jīng)出臺(tái)了很多相關(guān)法律政策來(lái)制約醫(yī)療事故的發(fā)生，但醫(yī)療事故的發(fā)生率呈現(xiàn)上升的趨勢(shì)，特別是一些醫(yī)療機(jī)構(gòu)，出于商業(yè)目的盲目追求效益而造成的醫(yī)療糾紛必須要重視。通過(guò)對(duì)醫(yī)療機(jī)構(gòu)醫(yī)療事故的發(fā)生率進(jìn)行研究，可以讓更多的人認(rèn)識(shí)到問(wèn)題的嚴(yán)重性，加強(qiáng)監(jiān)管力度，從而找到醫(yī)療事故發(fā)生率比較低的醫(yī)療機(jī)構(gòu)進(jìn)行就醫(yī)，希望一些醫(yī)院能夠認(rèn)識(shí)到問(wèn)題的嚴(yán)重性，從源頭上進(jìn)行治理，因此對(duì)醫(yī)療事故發(fā)生率的區(qū)間估計(jì)尤為重要。通過(guò)分析發(fā)現(xiàn)，醫(yī)療事故的發(fā)生率和統(tǒng)計(jì)學(xué)中的成數(shù)非常相似，這就讓我們有了一個(gè)新的角度對(duì)醫(yī)療事故的發(fā)生率進(jìn)行估計(jì)?？墒菍?duì)醫(yī)療事故發(fā)生率進(jìn)行點(diǎn)估計(jì)不夠精確和可靠，經(jīng)常用區(qū)間估計(jì)來(lái)彌補(bǔ)點(diǎn)估計(jì)的不足，這樣一來(lái)，問(wèn)題就變成了從研究總體成數(shù)的置信區(qū)間入手研究醫(yī)療事故發(fā)生率的置信區(qū)間。

抽樣調(diào)查研究的基本問(wèn)題之一就是對(duì)總體成數(shù)進(jìn)行區(qū)間估計(jì)，因此選擇較好的區(qū)間估計(jì)方法顯得尤為重要。曲衛(wèi)彬在1993年闡述了總體成數(shù)估計(jì)中兩個(gè)不容忽視的問(wèn)題，既要明確總體成數(shù)的含義，又要弄明白成數(shù)估計(jì)中應(yīng)用正態(tài)近似的條件，這對(duì)我們接下來(lái)的研究也很有啟發(fā)，要首先清楚成數(shù)的具體含義和條件才能更好的研究[1]。2012年黃良文等在《統(tǒng)計(jì)學(xué)》一書中介紹了正態(tài)近似法(Z估計(jì)法)對(duì)總體成數(shù)的區(qū)間估計(jì)方法，但此方法的使用必須滿足大樣本條件，而在小樣本的情況此方法便無(wú)能為力了，這就促使我們尋找適應(yīng)性更好的區(qū)間估計(jì)方法[2]。朱惠健通過(guò)抽樣分布理論，用T估計(jì)法對(duì)總體成數(shù)的區(qū)間估計(jì)方法(Z估計(jì)法)進(jìn)行修正[3]。我們發(fā)現(xiàn)這些方法都是基于大樣本理論基礎(chǔ)上的正態(tài)近似方法，針對(duì)區(qū)間估計(jì)，我們可以從構(gòu)造新的統(tǒng)計(jì)量入手，但是此方法的難度是相當(dāng)大的，于是想從小樣本理論出發(fā)，找出更理想的方法對(duì)總體成數(shù)進(jìn)行區(qū)間估計(jì)。

近年來(lái)國(guó)內(nèi)外已經(jīng)開(kāi)始研究小樣本理論，此理論方法的應(yīng)用可以大大提高效率，減少誤差。自從Daniels提出鞍點(diǎn)逼近以來(lái)，小樣本情況得到了進(jìn)一步發(fā)展，此方法在樣本量很小的情況，鞍點(diǎn)逼近的效果明顯要高于正態(tài)逼近[4]。1994年，Barndorff-Nielsen和Cox推導(dǎo)出極大似然估計(jì)密度函數(shù)的鞍點(diǎn)逼近式，推動(dòng)了鞍點(diǎn)逼近方法在統(tǒng)計(jì)學(xué)中的廣泛應(yīng)用[5]。Casella使用傅立葉反演公式和埃奇沃斯展開(kāi)，推導(dǎo)出鞍點(diǎn)近似于單個(gè)隨機(jī)變量的密度，并且還展示了用于近似指數(shù)族中最大似然估計(jì)密度的技術(shù)[6]。Butler系統(tǒng)地介紹了鞍點(diǎn)逼近的方法和應(yīng)用，并給出了高階漸進(jìn)推理證明[7]。

羅玉波等利用鞍點(diǎn)逼近的方法來(lái)近似廣義卡方型混合隨機(jī)變量的密度函數(shù)和分布函數(shù)，并通過(guò)蒙特卡洛模擬將得到的結(jié)果與正態(tài)近似法得出的結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)鞍點(diǎn)逼近方法整體表現(xiàn)效果出色，尤其是密度函數(shù)尾部區(qū)域逼近效果更為突出[8]。王玉琢將鞍點(diǎn)逼近法應(yīng)用到Pareto分布和非中心卡方分布中，給出Pareto分布極大似然估計(jì)的密度函數(shù)，以及非中心卡方分布的密度函數(shù)和分布函數(shù)的鞍點(diǎn)逼近式[9]。錢政超等人是在二項(xiàng)采樣下利用鞍點(diǎn)逼近方法構(gòu)造相對(duì)風(fēng)險(xiǎn)的置信區(qū)間，其研究方法是用兩組相互獨(dú)立的隨機(jī)變量且服從二項(xiàng)分布的情況下，用其分布參數(shù)比構(gòu)成相對(duì)風(fēng)險(xiǎn)(即兩種條件下一個(gè)事件分別發(fā)生的概率之比)，然后用鞍點(diǎn)逼近對(duì)相對(duì)風(fēng)險(xiǎn)的置信區(qū)間進(jìn)行構(gòu)造，對(duì)比并證明了小樣本情況下此方法的優(yōu)越性[10]。孟令賓等人主要是做了二項(xiàng)抽樣下優(yōu)勢(shì)比(兩種條件下各自的優(yōu)勢(shì)之比)的研究，并將其應(yīng)用到流行病學(xué)的相關(guān)領(lǐng)域中，證明鞍點(diǎn)逼近法(特別是在小樣本情況下)對(duì)優(yōu)勢(shì)比的區(qū)間估計(jì)較其他方法精確[11]。

對(duì)總體成數(shù)的區(qū)間估計(jì)通常是用正態(tài)近似法，但此方法在小樣本情況下估計(jì)的效果差強(qiáng)人意。所以本文引入樞軸量法和鞍點(diǎn)逼近方法對(duì)總體成數(shù)進(jìn)行區(qū)間估計(jì)，基于鞍點(diǎn)逼近法的小樣本特點(diǎn)，本文將其應(yīng)用在總體成數(shù)置信區(qū)間的估計(jì)，并將其與正態(tài)近似法、樞軸量法進(jìn)行對(duì)比，然后進(jìn)行模擬和實(shí)例分析得出：鞍點(diǎn)逼近法較其他方法對(duì)總體成數(shù)的區(qū)間估計(jì)效果較優(yōu)，特別是小樣本條件下，估計(jì)的區(qū)間長(zhǎng)度和覆蓋率整體表現(xiàn)最佳。相比錢政超和孟令賓等人的研究，本文的研究對(duì)象發(fā)生了變化，并將該方法應(yīng)用到對(duì)醫(yī)療事故發(fā)生率的研究中，本研究具有理論和現(xiàn)實(shí)意義。

二、總體成數(shù)置信區(qū)間的構(gòu)造方法

(1)

(2)

已知總體成數(shù)服從二項(xiàng)分布，沿著這條主線再結(jié)合二項(xiàng)分布的性質(zhì)和特點(diǎn)可以對(duì)總體成數(shù)的置信區(qū)間進(jìn)行構(gòu)造。本文分別從樞軸量累計(jì)分布法、正態(tài)近似法(Z估計(jì)法和T估計(jì)法)和鞍點(diǎn)逼近法對(duì)成數(shù)的置信區(qū)間進(jìn)行構(gòu)造，在進(jìn)行構(gòu)造前先介紹本文中用到的兩個(gè)重要引理。

引理1[12]如果G(t，θ)是θ的連續(xù)、嚴(yán)格減函數(shù)，那么

(1)關(guān)于θ的方程G(T-0，θ)=1-α的解是θ的置信水平為1-α的置信下限；

(2)關(guān)于θ的方程G(T，θ)=α的解是θ的置信水平為1-α的置信上限；

引理2[13]設(shè)隨機(jī)變量X～b(n，p)，Y～b(m，p)，且X與Y相互獨(dú)立，證明Z=X+Y～b(m+n，p)。

(一)樞軸量累積分布法

由式(1)可得每個(gè)樣本點(diǎn)服從兩點(diǎn)分布，設(shè)x1，x2，…，xn是來(lái)自參數(shù)為p且服從伯努利分布的獨(dú)立隨機(jī)樣本。結(jié)合引理1可求總體成數(shù)p的置信區(qū)間，首先定義p的充分統(tǒng)計(jì)量為：

因此，不難驗(yàn)證下列等式成立。

(3)

=Be(p|k，n-k+1)

(4)

同理可得置信上限為：

(二)正態(tài)近似法

假設(shè)xi相互獨(dú)立，在大樣本下(np>5，n(1-p)>5)，利用中心極限定理，則可把二項(xiàng)分布問(wèn)題轉(zhuǎn)化成正態(tài)分布問(wèn)題近似地去求解，但是目前這種正態(tài)近似法并不能保證對(duì)置信區(qū)間的可信度和準(zhǔn)確性，在實(shí)際計(jì)算中很可能造成難以接受的誤差，而為了減少誤差我們目前唯一的辦法只能使樣本量盡可能的大。

其服從標(biāo)準(zhǔn)正態(tài)分布，由于p未知，可用樣本方差代替總體方差，則總體成數(shù)p的置信水平為1-α的置信區(qū)間，為：

2.T估計(jì)法[3]。在大樣本下，由正態(tài)分布的可加性和中心極限定理可得：

(5)

(6)

將式(6)帶入式(5)可得總體成數(shù)p的置信水平1-α的置信區(qū)間為：

(三)鞍點(diǎn)逼近下總體成數(shù)置信區(qū)間的構(gòu)造

鞍點(diǎn)逼近是做漸進(jìn)分析的一個(gè)經(jīng)常用到的工具，它具有公式簡(jiǎn)潔、計(jì)算效率高和逼近效果更加準(zhǔn)確的優(yōu)點(diǎn)。1999年Goutis系統(tǒng)地介紹了鞍點(diǎn)逼近的理論推導(dǎo)，這里就簡(jiǎn)單介紹本文用到的鞍點(diǎn)逼近的相關(guān)理論基礎(chǔ)。假設(shè)存在一個(gè)非負(fù)函數(shù)f(x)，我們想要在某點(diǎn)x0處逼近它的值。這里需要用到泰勒展開(kāi)式思想，但是我們不是針對(duì)f(x)進(jìn)行展開(kāi)，而是在f(x)為正的情況下，使其滿足h(x)=lnf(x)，即f(x)=eh(x)，并選擇x0作為展開(kāi)點(diǎn)，我們得到：

f(x)≈

(7)

(8)

f(x)

(9)

借用公式(9)，根據(jù)Fourier逆轉(zhuǎn)公式，通過(guò)適當(dāng)?shù)姆e分變換，可以得到密度函數(shù)f(x)的Deniels的鞍點(diǎn)逼近公式為：

采用了Lugannani等人介紹的鞍點(diǎn)逼近公式[14]，可以得到累積分布函數(shù)的鞍點(diǎn)逼近公式為：

(10)

其中，Z=n1-np。根據(jù)累積生成函數(shù)的逼近方法，首先求出Z的矩生成函數(shù)：

φz(θ)=E(eθz)=E(eθ(n1-np))

=E(eθn1)E(e-θnp)

=(peθ+(1-p))ne-θnp

(11)

由式(11)可以求出關(guān)于隨機(jī)變量Z的累積生成函數(shù)為：

Kz(θ)=lnφz(θ)=nln(peθ+(1-p))-θnp

因此，可以求出隨機(jī)變量Z的一階、二階導(dǎo)數(shù)為：

peθ(n-n1)-n1(1-p)=0

三、蒙特卡洛模擬

接下來(lái)，我們將通過(guò)蒙特卡洛模擬比較樞軸量法、正態(tài)近似法(Z估計(jì)法和T估計(jì)法)和鞍點(diǎn)逼近法對(duì)總體成數(shù)置信區(qū)間的估計(jì)效果。關(guān)于對(duì)置信區(qū)間的估計(jì)好壞，我們通常用區(qū)間長(zhǎng)度和覆蓋率進(jìn)行評(píng)判，希望得到的覆蓋率越高越好，區(qū)間長(zhǎng)度越短越好。本文主要考慮如下情形：當(dāng)樣本量n=30，50，100，300的情況下，分別對(duì)p=0.05，0.25，0.5，0.75的情況進(jìn)行比較分析，具體模擬結(jié)果見(jiàn)表1和表2。

表1 不同的估計(jì)方法對(duì)總體成數(shù)P<0.5時(shí)區(qū)間估計(jì)的模擬結(jié)果表

表2 不同的估計(jì)方法對(duì)總體成數(shù)P≥0.5時(shí)區(qū)間估計(jì)的模擬結(jié)果表

對(duì)蒙特卡洛模擬的結(jié)果進(jìn)行分析可以得到：小樣本情況下，即n=30，50的情形，鞍點(diǎn)逼近法較Z估計(jì)法和T估計(jì)法的區(qū)間長(zhǎng)度短，并且覆蓋率近次于樞軸量法。Z估計(jì)法雖然較T估計(jì)法區(qū)間長(zhǎng)度短，但有時(shí)候覆蓋率表現(xiàn)并不是太理想。較大樣本情況下，即n=100，n=300的情形，隨著p的增大，鞍點(diǎn)逼近的區(qū)間長(zhǎng)度略有增加，但相比其他三種方法仍然是最短置信區(qū)間長(zhǎng)度，整體表現(xiàn)最優(yōu)。此外，這幾種方法普遍更適用于小概率事件。

四、實(shí)例分析

近年來(lái)，醫(yī)療事故頻繁發(fā)生，其發(fā)生率是一個(gè)重要的指標(biāo)，它反映了醫(yī)療機(jī)構(gòu)的工作人員對(duì)工作認(rèn)真負(fù)責(zé)的態(tài)度。這里列舉兩個(gè)例子，對(duì)醫(yī)療機(jī)構(gòu)醫(yī)療事故的發(fā)生率進(jìn)行實(shí)例研究，讓更多的人認(rèn)識(shí)到問(wèn)題的嚴(yán)重性，政府部門應(yīng)該加強(qiáng)監(jiān)管力度，患者應(yīng)到醫(yī)療事故發(fā)生率比較低的醫(yī)療機(jī)構(gòu)就醫(yī)。為了更加精確地算出相關(guān)醫(yī)療事故的發(fā)生率，這里比較了上述四種置信區(qū)間的構(gòu)造方法，具體應(yīng)用如下。

通過(guò)上面的三種方法分別構(gòu)造p的95%的置信區(qū)間見(jiàn)表3。表3表明：小樣本下，樞軸量累積分布法求出的置信區(qū)間比較可靠，但區(qū)間長(zhǎng)度太大;Z估計(jì)法和T估計(jì)法的區(qū)間長(zhǎng)度很接近，但相比鞍點(diǎn)逼近方法，鞍點(diǎn)逼近法求出的p的覆蓋率較高，區(qū)間長(zhǎng)度最短。

表3 不同估計(jì)方法對(duì)婦產(chǎn)科醫(yī)療事故發(fā)生率的區(qū)間估計(jì)結(jié)果表

通過(guò)上面的三種方法分別構(gòu)造p的95%的置信區(qū)間，見(jiàn)表4。由表4可知，大樣本情況下，鞍點(diǎn)逼近在區(qū)間長(zhǎng)度上有明顯的優(yōu)勢(shì)，覆蓋率僅次于樞軸量法。

表4 不同估計(jì)方法對(duì)醫(yī)療機(jī)構(gòu)責(zé)任比率的區(qū)間估計(jì)結(jié)果表

五、結(jié) 論

本文重點(diǎn)介紹了用樞軸量法和鞍點(diǎn)逼近法構(gòu)造的總體成數(shù)的置信區(qū)間，并通過(guò)模擬得出如下結(jié)論：樞軸量法的覆蓋率較高，但區(qū)間長(zhǎng)度較長(zhǎng)，在實(shí)際應(yīng)用中并不是較好的選擇；相比正態(tài)近似法對(duì)總體成數(shù)的區(qū)間估計(jì)方法，鞍點(diǎn)逼近法無(wú)論是在平均區(qū)間長(zhǎng)度，還是在覆蓋率上都表現(xiàn)出明顯的優(yōu)勢(shì)，尤其是在小樣本情況下，鞍點(diǎn)逼近法對(duì)總體成數(shù)置信區(qū)間的估計(jì)精度明顯提高，是一種較好的區(qū)間構(gòu)造方法。另外，研究結(jié)果表明，這幾種方法普遍適用于小概率事件，特別是正態(tài)近似法和鞍點(diǎn)逼近法。通過(guò)對(duì)兩個(gè)醫(yī)療事故發(fā)生率的實(shí)例分析，我們認(rèn)為，根據(jù)抽樣分布理論對(duì)醫(yī)療事故發(fā)生率的研究，可以轉(zhuǎn)化到對(duì)總體成數(shù)的估計(jì)問(wèn)題，研究結(jié)果表明，鞍點(diǎn)逼近法構(gòu)造醫(yī)療事故的發(fā)生率有較好區(qū)間估計(jì)效果，特別是對(duì)一些小概率、小樣本事件的研究，此方法具有普遍適用性。

本文的創(chuàng)新之處主要是運(yùn)用鞍點(diǎn)逼近法對(duì)總體成數(shù)的置信區(qū)間進(jìn)行估計(jì)，并且通過(guò)蒙特卡洛模擬和實(shí)例分析，將其與正態(tài)近似法、樞軸量法進(jìn)行比較。本文的研究不足之處是在模擬過(guò)程中，當(dāng)抽取的樣本較少并且成數(shù)的值也較小時(shí)，容易出現(xiàn)抽不到所需樣本的情況。