戴明鋒，金勇進，孫婕

（中國人民大學a．統(tǒng)計學院；b．應用統(tǒng)計科學研究中心，北京 100872）

一、引 言

在抽樣調查中，目標參數(shù)的估計通常有兩種途徑，一種是傳統(tǒng)的基于頻率學派的隨機化抽樣估計。Neyman提出了基于設計概率的抽樣估計理論［1］7。這種理論認為觀測樣本的取值是固定的，沒有誤差的，唯一的隨機性來源于抽取的樣本。設計概率在頻率學派的理論中起著非常重要的作用，很多估計量、比如總體均值、方差等的估計都是基于設計概率計算的，這種理論對目標量的統(tǒng)計分布不做專門的假設，因此在大樣本的情況下，具有很好的性質。比如估計量的無偏性、穩(wěn)健性、最小方差性等。π估計就是基于設計概率得出的計算總體估計量的一般公式。另一種目標量的估計基于模型的抽樣估計，這種方法將總體取值視為隨機過程，認為有限總體是從超總體模型中抽取的一個隨機樣本，是一個或多個隨機過程的一次實現(xiàn)［2］。polya后驗方法是基于模型估計的一種特殊形式，屬于無信息貝葉斯方法。這種方法首先要對總體的分布進行一個假定，認為在樣本抽取出來以后，后驗分布即是在已知觀測樣本的條件下，沒有觀測樣本的條件分布，抽樣估計不應該依賴于設計概率。當總體只有少量或者幾乎沒有先驗信息的情況下，polya后驗方法能夠得到很好地應用［3］。Rubin認為，polya后驗方法類似于Bootstrap方法［4］。在國外，polya后驗方法在分位數(shù)估計、小域估計中都得到了很好地應用［5－6］；而在國內，目前從事這方面的研究還比較少。本文對polya后驗估計方法進行模擬研究，探討了polya后驗估計在有限總體抽樣中的應用，并和運用比較成熟的Bootstrap方法進行比較，研究運用polya后驗估計方法的一些估計效果。

二、polya后驗抽樣

（一）polya后驗抽樣過程

假設有兩個盒子，第一個盒子中裝有n個相同的球，每個球的觀測值已知，記為y1，y2，…，yn，設這n個球為觀測的一個樣本。第二個盒子裝有N－n個球（n＜＜N），觀測值未知，現(xiàn)在分別從兩個盒子中隨機地各取一個球，假設從第一個盒子中取出的球的觀測值為yi，令從第二個盒子中取出的球的值和從第一個盒子中取出的球的值相等，都為yi，然后將這兩個球都放回第一個盒子中。再分別從兩個盒子中隨機地各取一個球，令從第二個盒子中取出的球的值和從第一個盒子中取出的球的值相同。將這個過程重復進行，直到這N－n個球都賦給了觀測值。因此在給定觀測樣本的情況下，運用這種polya抽樣就產生了未觀測值的一個分布，這個分布稱為polya后驗分布。當這個過程完成后，就得到了總體的一個完全觀測值。抽樣過程如圖1所示。

圖1 polya后驗抽樣過程示意圖

設m（y1，y2，…，yn）為邊際密度函數(shù)，則給定一個樣本統(tǒng)計量Z＝（s，y（s））和樣本單元值y（s）＝（y1，y2，…yn）后，參數(shù)θ的后驗密度為：

（二）polya后驗估計過程

由于polya后驗抽樣估計屬于無信息bayes估計，在進行點估計時，通常有以下步驟：1．確定參數(shù)的無信息先驗分布；2．根據(jù)觀測的樣本信息，結合參數(shù)的無信息先驗分布計算出未觀測的樣本的條件分布；3．定義一個損失函數(shù)；4．在損失函數(shù)下根據(jù)條件分布計算出相應的估計量。

設U＝｛1，2，…N｝表示容量為N 的有限總體，s＝｛1，2，…n｝為n的一個樣本，樣本單元數(shù)計為n（s），y（s）＝（y1，y2，…yn）為樣本單元值，向量Y＝（y1，y2，…，yN）T是參數(shù)空間 Θ中的一個樣本點，樣本關于參數(shù)的所有信息包含在統(tǒng)計量Z＝（s，y（s））中。設B為參數(shù)空間Θ上的一個Borelσ代數(shù)，p為可測空間（Θ，B）上的一個測度。Bayes學派認為參數(shù)空間Θ上相應的先驗分布屬于先驗分布族Π（θ）。

在具體計算條件分布時，首先需要確定先驗分布和似然函數(shù)。無信息先驗分布的分布函數(shù)可以選擇π（θ）＝1或者選擇Jeffreys先驗分布，也可以選擇其他的形式［7］。定義似然函數(shù)為：

記μ（y）為總體均值，在平方損失的情況下，總體均值的平均值為：

總體均值的方差估計為：

其中vs為樣本方差，計算公式為：

三、polya后驗方法特點的模擬研究

在復雜的統(tǒng)計問題中，有時候進行推導運算是很困難的，這時候可以借助統(tǒng)計模擬很好地解決問題。目前，統(tǒng)計模擬方法廣泛用于很多統(tǒng)計問題的研究中［8］。本文利用R軟件進行統(tǒng)計模擬，主要計算均值估計量、均值的置信區(qū)間、真值（總體均值）覆蓋率。

在polya后驗估計中，給定一個樣本，產生一個完整的模擬總體，重復這個過程多次，進行點估計時，取每次估計值的平均值作為最終點估計值。在構造置信區(qū)間時，這里選擇根據(jù)分位數(shù)構造置信區(qū)間。比如要構建均值的95%的置信區(qū)間，根據(jù)點估計得到一系列的均值，按照從小到大的順序排序。取其下側0．025分位數(shù)作為下側置信區(qū)間，取其上側0．975分位數(shù)作為上側置信區(qū)間。

設有兩個總體，記為pop1和pop2，各包含10 000個單元。pop1服從均值為100方差為20的正態(tài)分布；pop2服從位置參數(shù)為20尺度參數(shù)為5的伽馬分布。在每個總體中，分別抽取樣本量為1 000、2 000、3 000、4 000的四個樣本，每種情況模擬500次。計算均值的估計值、置信區(qū)間以及真值（總體均值）覆蓋率。模擬結果如表1。

表1 polya后驗估計的模擬結果

在設定好隨機種子后，pop1中，本文中利用總體10 000個樣本單元直接計算出的總體均值是99．984；pop2中，利用總體10 000個樣本單元直接計算出的總體均值是3．984。對比表1中顯示的估計結果，可以看出采用polya后驗估計可以精確地估計出總體的均值。當樣本量增加時，均值的平均置信區(qū)間會變短，但是真值覆蓋率會增加。pop1中，可以看到當樣本容量由3 000增加到4 000后，平均區(qū)間長度和真值覆蓋率并沒有發(fā)生變化，說明當樣本量增加到一定程度后，采用這種方法估計的精確度就會穩(wěn)定在一定水平上。pop2中，當樣本容量由3 000增加到4 000后，平均區(qū)間長度和真值覆蓋率雖然發(fā)生了一定的變化，但是變化幅度已經(jīng)非常小了。

四、polya后驗估計與Bootstrap方法估計的模擬比較

（一）Bootstrap方法的基本原理

Bootstrap方法是一類非參數(shù)蒙特卡洛方法，它通過再抽樣對總體分布進行估計，再抽樣方法將觀測到的樣本視為一個有限總體，從中進行隨機再抽樣來估計總體的特征以及對總體參數(shù)作出統(tǒng)計推斷［9］。

假設y＝（y1，y2…，yn）為從一個總體分布F（y）中觀測的樣本，Y＊為從y中隨機選擇的一個樣本，則P（Y＊＝y(tǒng)i）＝，i＝1，2，…，n，從y中有放回的再抽樣得到隨機樣本Y1＊，Y2＊，…，Yn＊，則隨機變量Y1＊，Y2＊，…，Yn＊為獨立同分布的隨機變量，服從｛y1，y2…，yn｝上的均勻分布。

由于經(jīng)驗分布函數(shù)Fn（y）是F（y）的估計，并且其本身是｛y1，y2…，yn｝上的均勻分布隨機變量Y＊的分布函數(shù)。因此Bootstrap重復下的經(jīng)驗分布函數(shù)Fn＊是Fn的逼近。從y中再抽樣，等價于從Fn中產生隨機樣本。

采用Bootstrap方法構造置信區(qū)間可以采用標準Bootstrap、百分位數(shù)Bootstrap、t百分位數(shù)Bootstrap和修正偏差后的百分位Bootstrap四種方法來估計。本文采用百分位數(shù)Bootstrap方法。

百分位數(shù)Bootstrap方法利用Bootstrap經(jīng)驗分布的和1－分位點分別是在1－α置信水平下統(tǒng)計量T的置信區(qū)間的上、下限。具體如下：假設取B＝1 000時，可以得到1 000個Bootstrap隨機樣本，將每個樣本的均值按由小到大的順序依次排列，可以得到順序統(tǒng)計量（i），（i＝1，2，…B），則和 1－分位點分別是在1－α置信水平下統(tǒng)計量T的置信區(qū)間的上下限。

根據(jù)伯努利大數(shù)定理，，當抽樣次數(shù)充分大時，這些區(qū)間中包含θ的真值的頻率接近于置信度（即概率）1－α，即在這些區(qū)間中包含θ的真值的區(qū)間大約有100（1－α）%個，不包含θ的真值的區(qū)間大約有100α%個。

（二）polya后驗估計與自Bootstrap方法估計的模擬比較

設有四個總體，記為pop1、pop2、pop3和pop4，各包含1 000個單元。pop1服從區(qū)間（－1，1）上的均勻分布，pop2服從均值為10方差為1的的正態(tài)分布，pop3服從位置參數(shù)為5尺度參數(shù)為1的伽馬分布，pop4服從參數(shù)為4的指數(shù)分布。在每個總體中，分別抽取樣本量為100、200的兩個樣本，每種情況模擬500次。計算在Bootstrap方法和polya后驗方法下均值的估計值、置信區(qū)間以及真值（總體均值）覆蓋率。模擬結果如表2。

表2顯示：總體來看，在pop1、pop2、pop3、pop4四個分布不同的總體中，采用Bootstrap方法和polya后驗估計估計時，在估計總體均值時，二者沒有什么大的區(qū)別，樣本量發(fā)生變化后，采用兩種方法估計的總體均值差異也差不多。在構造的平均置信區(qū)間方面，采用polya估計方法構造的置信區(qū)間區(qū)間長度相對較短，而采用Bootstrap方法構造的置信區(qū)間區(qū)間長度相對較長一些。在真值覆蓋率方面，采用Bootstra方法估計時，真值覆蓋率要高一些，但是采用polya后驗方法也能以較高的概率覆蓋真值。

表2 polya后驗估計與Bootstrap估計的結果比較表

五、小 結

本文研究了polya后驗方法在有限總體抽樣中的應用，通過模擬研究，發(fā)現(xiàn)polya后驗估計可以精確地估計出總體的均值，據(jù)此構造的置信區(qū)間相比較根據(jù)Bootstrap方法構造的置信區(qū)間長度較短一些，置信區(qū)間能夠以較高的概率覆蓋真值。對于polya后驗估計在估計方差時效果如何，polya后驗估計的使用范圍，以及polya后驗方法在實際問題中的應用等一些問題，還需要進一步的研究。

［1］ Cochran W G．抽樣技術［M］．張堯廷，吳輝，譯．北京：中國統(tǒng)計出版社，1985．

［2］ 金勇進，賀本嵐．復雜抽樣推斷方法體系的比較研究［J］．統(tǒng)計與信息論壇，2011（10）．

［3］ Glen Meeden．A Noninformative Bsyesian Approach to Domain Estimation［J］．Journal of Statistical Planning and Inference，2005（12）．

［4］ Rubin D．The Bayesian Bootstrap［J］．Annals of statistics，1981（9）．

［5］ David Nelson，Meeden Glen．Noninformative Nonparametric Quantile Estimation for Simple Random Samples［J］．Journal of Journal of Statistical Planning and Inference，2006（9）．

［6］ Xiaoyin Wang，Chong Z He，Dongchu Sun．Bayesian Population Estimation for Small Sample Capture-Recapture Data Using Noninformative Priors［J］．Journal of Statistical Planning and Inference，2007（4）．

［7］ 茆詩松，王靜龍，濮曉龍．高等數(shù)理統(tǒng)計 ［M］．北京：高等教育出版社，1998．

［8］ 楊貴軍，劉艷林，王清．捕獲在捕獲抽樣估計量的模擬研究 ［J］．統(tǒng)計與信息論壇，2011（3）．

［9］ 謝益輝，朱鈺．Bootstrap方法的歷史發(fā)展和前沿研究［J］．統(tǒng)計與信息論壇，2008（2）．