基于蒙特卡洛方法的假設檢驗問題探討

魏艷華，王丙參，邢永忠

（天水師范學院a.數學與統計學院；b.電子信息與電氣工程學院，甘肅 天水 741001）

0 引言

在統計推斷中，由于目標統計量的真實分布往往很難獲得，這就對目標量的估計與檢驗造成了很大困難。假設檢驗的推斷基礎是樣本，是一種不完全歸納推斷，盡管在大多數情況下得到的結論是正確的，但也有可能犯錯：一類是棄真錯誤（第I型錯誤），一類是存偽錯誤（第II型錯誤）。通常，很難獲得統計量在已知條件（比如原假設成立或備擇假設成立條件）下的精確分布或極限分布，因此，也無法確定接受原假設的臨界值與犯錯的概率。這時，利用蒙特卡洛方法逼近統計量的分布不失為一種可行辦法，其主要思想是通過產生參考數據，構造新分布去逼近基于觀測數據的統計量分布[1-4]。本文將對此方法進行探討。

1 假設檢驗錯誤率的蒙特卡洛估計

設 (Ω，?X)為可測空間，樣本X取值于 (Ω，?X)，其分布族為{Pθ，θ∈Θ}，又因總體與樣本同分布，故總體分布族也為{Pθ，θ∈Θ}。一個統計假設是關于樣本X的函數，或者說，關于參數θ的一個命題：H0:θ∈Θ0，其中Θ0已知，是Θ的一個真子集。設非空集Θ1?Θ-Θ0，則命題H1:θ∈Θ1稱為H0的備擇假設或對立假設。注意，在一般情況下Θ1=Θ-Θ0，采用Θ1是Θ-Θ0的子集具有更大的靈活性。如果Θ1是Θ-Θ0的真子集，則H0成立未必H1不成立，這時，H1是本文最關心的非H0情況。基于樣本x而做出是否拒絕原假設H0，稱為對H0進行統計檢驗。這樣，一個統計檢驗就是把樣本空間Ω分解為兩個不交子集 ΩA與 ΩB，當x∈ΩA時接受H0，當x∈ΩB時拒絕H0，故ΩA與ΩB分別稱為檢驗的接受域與拒絕域。按照隨機化決策的思想，可引入更一般的檢驗方法：在樣本空間Ω上定義一個取值為[0，1]的?x可測的函數φ(x)，對于樣本x計算φ(x)，然后以概率φ(x)否定原假設H0，以1-φ(x)的概率接受原假設，這樣的φ(x)稱為檢驗函數，即檢驗模型。如果φ(x)∈(0，1)，檢驗稱為隨機化的。如果φ(x)取值只是0或1，則檢驗退化為一般的檢驗，即非隨機化檢驗，這時，檢驗的接受域為{x:φ(x)=0}，拒絕域為{x:φ(x)=1}。

對于正態總體，由于構造的檢驗統計量往往已知，所以可通過解析方法獲得兩類錯誤的表達式，但是，如果假設檢驗總體的正態性假設不滿足或者不能通過解析方法獲得兩類錯誤的表達式，此時，利用蒙特卡洛方法模擬兩類錯誤便是一個很好的選擇，這就是蒙特卡洛檢驗（MCT）[4]。MCT的關鍵在于利用蒙特卡洛方法逼近統計量的分布，其基本思想是利用產生的參考數據構造新的分布去逼近基于觀測數據的統計量分布，因此如何產生參考數據是至關重要的。當然，如果在原假設成立或備擇假設成立條件下總體的分布已知，則產生觀測數據是比較容易的，只需要生成已知總體的樣本觀測值即可。在參數的情況下，如果沒有討厭參數，MCT可能達到精確的顯著性水平，即使與一致最優檢驗相比，它的功效也是非常高的；在討厭參數存在情況下，MCT仍然可能達到精確的顯著性水平。在參數族情況下，不論是否有討厭參數，不論統計量的漸近分布是否樞軸，利用MCT逼近的誤差比統計量的漸近引起的誤差小，且可區分n-0.5的速度逼近備擇假設。這些結論進一步奠定了MCT法的理論依據。

利用蒙特卡洛方法模擬棄真錯誤α的步驟：

（1）確定原假設H0成立條件下總體的分布函數。

（2）抽取樣本容量為n的樣本的樣本觀測值，即生成n個總體分布隨機數。用臨界值方法進行假設檢驗，確定棄真錯誤是否發生，即觀測是否拒絕原假設H0，并記錄實驗結果：

（3）重復第（2）步N次。

同理，利用蒙特卡洛方法模擬存偽錯誤β的步驟：

（1）確定原假設H0不成立條件下總體的分布函數。

（2）抽取樣本容量為n的樣本的樣本觀測值，即生成n個總體分布隨機數。用顯著性水平α與相應臨界值進行假設檢驗，確定存偽錯誤是否發生，即觀測是否未拒絕原假設H0，并記錄實驗結果：

（3）重復第（2）步N次。

注意，由于原假設H0不成立具有多種可能，因此在分析存偽錯誤時應指明選取的備擇假設。

設φ(x)是檢驗H0的檢驗函數，則函數：

gφ(θ)=Pθ(否定H0)

稱為檢驗方法φ(x)的功效函數。通俗地說，樣本觀測值落入拒絕域內的概率稱為該檢驗的功效函數，即：

gφ(θ)=Pθ(否定H0)=Pθ(X∈W)，θ∈Θ

其中X為樣本觀測值，W為拒絕域。顯然，當θ∈Θ0時，gφ(θ)=α=α(θ)，當θ∈Θ1時，gφ(θ)=1-β=1-β(θ)。可見，求功效函數就是求檢驗犯兩類錯誤的概率。如果利用蒙特卡洛方法估計出兩類錯誤的概率，也就估計出了功效函數。

例 1：假定X～P(λ) ，利用樣本觀測值檢驗假設：H0:λ=5VSH1:λ＞5 。令，機械地利用中心極限定理會得出：當Z≥1.645時拒絕原假設H0。顯然，對于確定的n（尤其當n較小時），統計量的精確分布是未知的，但利用MCT方法進行統計推斷。下面利用MCT方法模擬棄真錯誤的概率α。每次模擬α時生成N=1000個參考值z，總共模擬1000個α值，直方圖見圖1。

圖1 利用MCT方法模擬棄真錯誤的概率（左圖n=10，右圖n=100）

當n=10，棄真錯誤α的均值為0.0551，標準差為0.0071，利用分位數方法估計的95%的置信區間為（0.0420,0.0690）；當n=100，棄真錯誤α的均值為0.0523，標準差為0.0070，利用分位數方法[5]估計的95%的置信區間為（0.0390,0.0660）。如果統計量服從標準正態分布，棄真錯誤的概率應為0.05。隨著n的增大，Z越來越近似于標準正態分布，棄真錯誤的概率也越來越可能接近0.05，方差也會相應的減少，模擬結果也驗證了此結論。注意，這些結論是在概率意義下成立的，在實際模擬中也會出現例外，尤其當n相差不大時。

圖2 λ∈[5.1,6.5]時的功效函數

對于λ∈[5.1,6.5]，利用蒙特卡洛方法估計存偽錯誤的概率，進而得出相應區間的功效函數曲線，如圖2所示，其中n=100，每次模擬存偽概率β時生成N=1000個參考值z。λ從5.1～6.2（間隔為0.1）的功效函數取值分別為：

（0.1070 0.2350 0.3820 0.5580 0.7210 0.8300 0.9260 0.9660 0.9870 0.9990 0.9990 1.0000）

可見，隨著λ的增大，功效函數遞增，從λ=6.2開始，功效函數取值幾乎為1，即當λ較大時利用中心極限定理進行近似檢驗的效率非常高，但是，λ較小時，效率則是非常低的，比如λ=5.1時，功效函數才為0.1070。

例2：（非平穩序列虛假回歸的隨機模擬）為正確理解虛假回歸，本文考慮一元動態線性回歸模型：yt=β0+β1xt+εt。為了檢驗模型的顯著性，需要檢驗：H0:β1=0 VSH1:β1≠0 。如果序列相互獨立，理論上應該接受原假設H0。如果檢驗結果恰好與理論相反，就會犯拒真錯誤。由于樣本的隨機性，拒真錯誤會一直存在。通常采用統計量進行顯著性檢驗，如果都平穩，則該統計量服從t(n)，其中n為樣本容量。如果不平穩，則檢驗統計量不再服從t(n)分布，且樣本方差遠遠大于t(n)分布的方差。如果仍用t(n)分布臨界值進行檢驗，拒真概率就會大大增加。這就會導致無法控制拒真概率，容易接受回歸模型顯著的錯誤結論，即產生虛假回歸現象。假如擬合兩個隨機游走模型序列：

圖3 檢驗統計量的頻數直方圖（左）與頻率直方圖（右）

圖4 檢驗統計量的頻數直方圖（左）與頻率直方圖（右）

拒真錯誤概率才0.0650，這說明在此非平穩序列場合，虛假回歸現象幾乎不存在，原因主要在隨機游動系數的選擇上。如果保持其他參數不變，令a=0.5，b=0.5，拒真錯誤概率為0.11700；令a=0.1，b=0.1，拒真錯誤概率為0.0500。可見，隨機游動系數越小，拒真錯誤概率越低，甚至會接近正常值，不出現虛假回歸現象。

2 檢驗p值的蒙特卡洛估計

由于假設檢驗的結論很簡單，在給定顯著性水平α下，要么接受原假設，要么拒絕原假設。然而，會出現這樣情況：在一個較大的顯著性水平得到拒絕原假設的結論，但是一個較小的顯著性水平下得到拒絕原假設的結論，因為顯著性水平變小，會使得檢驗的拒絕域變小。為此，本文引入檢驗p值，即利用樣本所能夠拒絕原假設的最小顯著性水平[6，7]。

考慮總體為F(.)的簡單隨機樣本假定原假設為H0:F(.)=G(.，θ)，其中θ是未知參數，G(.)為已知函數。對于這個假設檢驗的任何統計量蒙特卡洛檢驗方法就是從中獨立產生參考數據并計算作為統計量的參考值，其中?是θ的估計值。不妨令=T0，T1，…，TN表示由蒙特卡洛方法得到的N個參考值，如果T值較大，拒絕原假設，則檢驗p值的估計為，其中k為T0，T1，…，TN大于等于T0的個數。如果?小于等于給定的顯著性水平α，則拒絕原假設；反之則否。同理，對于其他檢驗情況（比如雙邊檢驗）不難進行相應調整。值得注意的是，參數蒙特卡羅檢驗的具體步驟非常類似20世紀80年代發展的參數自助近似[8]。

對于非參數或半參數的情形，在原假設成立情況下很難模擬參考數據，進而也無法計算統計量對應的MCT的條件統計量。這是因為即使在原假設成立下，模型不能用含有幾個未知數的具體模型刻畫。自助法是解決此問題方法之一，它是完全的非參數統計方法，對模型結構和數據結構的限制很少。但是，如果模型不是非參數的，而是半參數，自助法就效率不高，這時可改進蒙特卡洛逼近，充分利用數據提供的信息，這也是目前研究的熱點之一。如果隨機向量X獨立可分解，即X=Y·Z依分布獨立，Y，Z獨立且已知，Y·Z表示Y，Z點乘，這時可采用非參數蒙特卡洛檢驗（NMCT）。記表示總體X的簡單隨機樣本，如果xi在原假設下可獨立分解為xi=yi·zi，則檢驗統計量NMCT的具體做法為：給定z1，…，zn，從分布Y中獨立生成一組參考數據于是相應統計量的值為假設T值較大，拒絕原假設；雙邊檢驗不難做相應調整。如果原始數據對應的T值為T0，通過蒙特卡洛方法生成N組參考數據，則檢驗p值的估計為其中k為T0，T1，…，TN大于等于T0的個數。給定顯著性水平α，只要就拒絕原假設。因為同分布于，而且在給定下，它們有相同的條件分布，因此檢驗可能精確有效。原因如下：

如果Ti之間有結，且k≤[α(N+1)]，那么T0至少比Ti中第個大的元素大，因此也有

對zi求積分，可得這說明在變量可獨立分解時，NMCT是可以精確有效的，而自助法就沒有這個優點。

例3：（分布族的卡方擬合優度檢驗）考慮檢驗H0:F(.)=G(.，θ)，其中θ是r維未知參數，在H0下X的所有可能取值Ω可分為k個互不相交的子集其中k＞r+1，記fi為n個樣本觀測值落入集合Ai的次數，則樣本觀測值落入集合Ai頻率為另一方面，利用假設分布函數G(.，θ)可計算：

某地區在夏季一個月中由100個氣象站報告的暴雨次數如表1所示，要求檢驗假設H0:F(x)=P(θ)。

表1 某地區100次觀測的暴雨次數

下面利用MCT方法，每次生成100個參考數據，共模擬10000次，結果如圖5所示，p值為0.8336。如果每次生成10000個參考數據，共模擬1000次，結果如圖6所示。注意，此時分組應變為Ai=P(X=i-1)，i=1，…，8，并把 ≥8作為一類，記為A9。由于最初樣本值只有100個，所以此時不能得到原檢驗的p值，只能模擬檢驗統計量的分布。

圖5 檢驗統計量的頻率分布圖（n=100）與χ2(3)密度曲線

顯然，當n=100時，檢驗統計量與漸近分布χ2(3)相比，右偏，厚尾，這會導致低估p值，故利用MCT可提高檢驗效率。由于p值0.8336顯著大于0.05，故可認為接受原假設，即某地區暴雨次數分布服從泊松分布。

圖6 檢驗統計量的頻率分布圖（n=10000）與χ2(7)密度曲線

顯然，當n=10000時，檢驗統計量與漸近分布χ2(7)吻合的非常好，這也驗證了檢驗統計量的漸近分布為χ2(k-r-1)。故在進行分布族的卡方擬合優度檢驗時，樣本觀測值要盡量多（建議大于2000），否則，建議使用MCT方法。

3 結論

本文利用蒙特卡洛方法產生參考數據去逼近統計量的分布，模擬假設檢驗的兩類錯誤發生概率與p值。在非平穩序列場合構造一元動態回歸模型，探討虛假回歸現象出現的原因與影響因素。對于分布族的卡方擬合優度檢驗，得出：當樣本量較少時，使用漸進分布誤差較大，建議使用蒙特卡洛檢驗方法。