基于Bootstrap思想的非參數(shù)檢驗p值的估計方法

呂書龍, 劉文麗

(福州大學數(shù)學與計算機科學學院,福建 福州 350116)

0 前言

參數(shù)型假設(shè)檢驗特別是基于正態(tài)總體的假設(shè)檢驗,無論是從理論上還是實踐上都已經(jīng)成為經(jīng)典. 參數(shù)型假設(shè)檢驗的前提要求總體服從的分布形式或分布族是已知的,這對構(gòu)造檢驗統(tǒng)計量并確定其分布和精確計算檢驗p值都提供了最有利的條件. 但在實際問題中,總體的分布通常是未知的,于是,人們就希望在總體分布不假定的前提下,基于數(shù)據(jù)本身來構(gòu)造反映總體特征及問題的統(tǒng)計量或信息,從而達到統(tǒng)計推斷的目的,這也是非參數(shù)統(tǒng)計的宗旨. 眾所周知,采用檢驗p值完成假設(shè)檢驗的推斷已成為一種最基本的模式. 作為非參數(shù)統(tǒng)計的一個主要分支,非參數(shù)檢驗在醫(yī)學、 生物學、 信息學、 金融、 管理、 教育等眾多領(lǐng)域應(yīng)用廣泛.

如何構(gòu)造檢驗統(tǒng)計量并計算檢驗p值,是非參數(shù)檢驗的一個關(guān)鍵問題,特別是在樣本量較小的情況下,這個問題尤為突出. 但本研究注意到由于非參數(shù)方法不假定總體分布,對原始數(shù)據(jù)的使用也不像參數(shù)型方法,這可能會導(dǎo)致非參數(shù)方法的有效性低于參數(shù)方法. 而通過重抽樣Bootstrap方法[1-2]可以很好地給出總體特征的估計及精度表示[3]、 區(qū)間估計以及假設(shè)檢驗等,甚至能發(fā)現(xiàn)常規(guī)檢驗方法中可能存在的不穩(wěn)定問題. 檢驗p值的估計是非參數(shù)假設(shè)檢驗的基本問題,Nilsson[4]提出檢驗p值的數(shù)據(jù)量縮減型Monte-Carlo非參數(shù)估計方法,并給出其在生物信息學上的應(yīng)用. Silverman[5]探討了基于Bootstrap的漸近樞軸統(tǒng)計量的經(jīng)驗分布函數(shù)的核方法及檢驗p值的估計. 本研究嘗試探討寬泛的檢驗p值估計方法及其適用范圍.

1 問題描述

設(shè)總體X～F(x),F(x)未知,X1,X2, …,Xn為其獨立同分布樣本,x1,x2, …,xn為樣本觀測值. 為推斷與總體相關(guān)的某些問題,通常要構(gòu)造樣本函數(shù)T=T(X1,X2, …,Xn),即統(tǒng)計量. 但由于F(x)未知,此時統(tǒng)計量T的研究基本上可歸入非參數(shù)領(lǐng)域.

假設(shè)檢驗應(yīng)用廣泛,很多的實際問題都可以轉(zhuǎn)化成假設(shè)檢驗而得到較好解決. 在作假設(shè)檢驗時,是否拒絕原假設(shè)H0,基本上都是通過先假定H0成立,再計算其檢驗p值并根據(jù)檢驗p值的大小來進行推斷. 各種統(tǒng)計軟件僅根據(jù)統(tǒng)計量T的分布給出檢驗p值,但是否拒絕原假設(shè)全部交給使用者自行裁定. 如果統(tǒng)計量T的分布難以確定,就會導(dǎo)致檢驗p值的計算變得無章可循,實際問題也就難以得到有效的解決.

本研究在總體分布未知且小樣本條件下,避開統(tǒng)計量的理論分布探討,單純從模擬仿真的角度,通過Bootstrap法產(chǎn)生的自助樣本,構(gòu)造“檢驗統(tǒng)計量”并結(jié)合隨機模擬、 經(jīng)驗分布和核密度估計[6]等方法實現(xiàn)原假設(shè)下的檢驗p值的估計,從而在一定程度上解決了上述假定條件下假設(shè)檢驗問題的推斷.

2 主要理論和工具

1) 格列汶科(Glivenko)定理給出用樣本的經(jīng)驗分布來估計總體的分布函數(shù)的理論,即若設(shè)總體X的分布函數(shù)為F(x),經(jīng)驗分布函數(shù)為Fn(x),則有

(1)

3) 設(shè)K(x)為定義在(-∞, +∞)上的一個Borel可測函數(shù),hn>0為常數(shù),則稱

(2)

為總體密度函數(shù)f(x)的一個核估計. 其中, K(x)稱為核函數(shù)； hn稱為窗寬. 利用核密度函數(shù)估計法可對T(X1, X2, …, Xn; Fn)的密度進行估計.

3 檢驗p值估計的基本思路和算法

先確定假設(shè)檢驗問題和類型,構(gòu)造檢驗統(tǒng)計量； 然后利用Bootstrap重抽樣得到檢驗統(tǒng)計量的模擬抽樣序列； 在抽樣序列的基礎(chǔ)上,利用Bootstrap法、 經(jīng)驗分布法和核密度估計法計算不同檢驗問題和類型對應(yīng)的檢驗p值.

3.1 三類檢驗問題

假定檢驗統(tǒng)計量為T(X1,X2, …,Xn;Fn),其一次試驗的值為Tx(x1,x2, …,xn),檢驗p值的直觀意義表示為發(fā)生比Tx(x1,x2, …,xn)更極端事件的概率. 對于三類假設(shè)檢驗問題,下面給出檢驗p值的形式表示.

1) 雙側(cè)檢驗,p1=P(T≥Tx(x1,x2, …,xn)),如果p1<=0.5,則檢驗p值=2p1,否則檢驗p值=2(1-p1)；

2) 右側(cè)檢驗,檢驗p值=P(T≥Tx(x1,x2, …,xn))；

3) 左側(cè)檢驗,檢驗p值=P(T≤Tx(x1,x2, …,xn)).

3.2 檢驗p值的模擬估計

3.2.1Bootstrap法估計檢驗p值

設(shè)T∈(Ti,Ti+1),則

(3)

3) 基于Bootstrap的百分位數(shù)法,給出三類檢驗問題的檢驗p值公式：

(4)

3.2.2 經(jīng)驗分布法估計檢驗p值

(5)

通過經(jīng)驗分布函數(shù)法得到的檢驗p值本質(zhì)上和Bootstrap法得到的檢驗p值是一致的,只是表達方式不同而已.

3.2.3 核密度法估計檢驗p值

p雙側(cè)=2p1(p1≤0.5)或2p2(p2≤0.5)；p右側(cè)=p2；p左側(cè)=p1

(6)

4 模擬分析與R實現(xiàn)

為了更好地說明本研究提出的算法,下面給出正態(tài)型參數(shù)檢驗和非參數(shù)檢驗的兩個例子,并與常規(guī)的假設(shè)檢驗進行比較.

例1 設(shè)總體X～N(100,σ2),其中σ2未知(模擬取σ2=16),其樣本為X1,X2,…,X20, 觀測值為100.09, 100.53, 101.84, 94.25, 101.35, 101.88, 100.96, 96.45, 102.87, 98.82, 99.88, 96.97, 105.99, 103.23, 100.60, 95.67, 101.89, 100.65, 97.48, 99.13, 試給出關(guān)于E(X)=100的三類檢驗問題的檢驗p值.

表1 常規(guī)算法與本研究算法的檢驗p值Tab.1 p value of conventional algorithm and our algorithm

如果研究的問題是： 設(shè)Y=X2,求關(guān)于E(Y)=10 016的三種假設(shè)檢驗問題的檢驗p值. 顯然使用常規(guī)的方法難以計算檢驗p值,而本研究算法卻可以很輕松地估計出檢驗p值,只需將統(tǒng)計量函數(shù)定義為

(注:EY=[E(X)]2+σ2=μ2+σ2)

例2 對于中風病人與健康成人血液中尿酸濃度數(shù)據(jù),見表2[7],問兩類人血液中的尿酸濃度的變異是否存在顯著差異?

表2 中風病人與健康成人血液中尿酸濃度數(shù)據(jù)

Tab.2 Data of uric acid concentration in the blood of stroke patients and healthy adults (g·L-1)

本例的假設(shè)檢驗問題可表示為

(7)

文獻[7]采用Moses檢驗法,并給出拒絕H0的推斷. Moses檢驗的作法是將兩樣本分組,保證每組長度都一樣,如各分成m1,m2組； 然后計算兩樣本各小組的離差平方和SSAi, SSBj,i=1, 2, …,m1；j=1, 2, …,m2； 再將SSAi, SSBj合并后計算樣本1的m1個組的秩和S,從而構(gòu)造Moses統(tǒng)計量TM

(8)

王星[7]認為如果兩組數(shù)據(jù)的方差存在很大差異,從平均看,其中一組的離差平方和比另一組的離差平方和要小. 文獻[7]以每組3個樣本點依據(jù)前后順序?qū)颖?分成4組,樣本2分成3組且棄用最后1個樣本點,并給出Mann-Whitney的Wα表法進行查表檢驗. 文獻[7]介紹了這種方法,并給出具體的計算過程. 本研究對其中的三個細節(jié)問題作深入探討,一是如何分組,二是不正好分組時棄用哪些樣本點,三是檢驗的結(jié)論是否會隨分組不同而不同. 本研究認為引入隨機分組和隨機棄用,結(jié)合大量模擬來進行檢驗更符合實際情況. 為了說明上述問題,設(shè)計R程序(見附錄2)并整理運行結(jié)果,如表3所示.

表3 文獻[7]方法的某1 000次隨機模擬的統(tǒng)計結(jié)果Tab.3 Statistical results by 1 000 stochastic simulations based on the method in literature [7]

若按照文獻[7]提供的方法,得到顯著性水平為0.05時的臨界值0和12,上述程序模擬1 000次,合計出現(xiàn)了383次接受原假設(shè)的情況(如表3所示),且經(jīng)過多次模擬表明,出現(xiàn)上述檢驗錯誤的頻率大致都在38%左右. 以上說明不同的隨機分組和取點會造成檢驗結(jié)論的不同,即上述三個細節(jié)問題必須慎重考慮. 依據(jù)文獻[7]構(gòu)造的隨機模擬表明該方法可能呈現(xiàn)出相當大的不穩(wěn)定性.

為此有必要充分考慮上述三個細節(jié)問題,并探討更加穩(wěn)定的算法來解決上述方差相等與否的檢驗問題. 對數(shù)據(jù)先作一下基本分析,計算其樣本均值、 樣本方差和樣本極差,直觀判斷兩組數(shù)據(jù)的差異性.

表4為基本統(tǒng)計量,從表4看出這兩組數(shù)據(jù)的各項統(tǒng)計量指標都存在較大的差異,檢驗的結(jié)論應(yīng)該拒絕原假設(shè)才合理. 為了避免出現(xiàn)上述隨機分組、 取點的干擾,引入刀切法思想,即構(gòu)造原樣本的舍一樣本方差：

； i=1, 2, …, n

(9)

(10)

若式(7)中原假設(shè)成立,則上述兩個TM值既不太大也不太小,否則可拒絕原假設(shè),程序見附錄3.

表5為檢驗p值的比較,從表5看出修正后的Moses檢驗的檢驗p值遠小于文獻[7]方法得到的檢驗p值,相對于文獻[7]而言有更強的理由拒絕原假設(shè).

表4 基本統(tǒng)計量Tab.4 Basic statistics

表5 檢驗p值的比較Tab.5 Comparison of p values

5 實例分析

風險價值方法 (value at risk, VaR)作為金融風險度量的一種重要方法,在金融風險管理中應(yīng)用廣泛. 金融資產(chǎn)價格變化的分布經(jīng)常出現(xiàn)尖峰厚尾的現(xiàn)象,若使用模型設(shè)定法(如價格變化服從t分布或正態(tài)分布等),估計往往不精確甚至不符合實際. 本研究以2016年上證指數(shù)前一日收盤價與后一日開盤價(共計244個數(shù)據(jù))的日對數(shù)收益率rxt=log(開盤價t)-log(收盤價t-1),t=2, 3, …, 244(合計243個數(shù)據(jù))為例,構(gòu)造這樣一個風險檢驗問題： 在置信度為95%時,H0:VaR=-3.40%,即以95%的置信度認為日對數(shù)收益率不低于-3.40%. 圖1(a)給出日對數(shù)收益率的核密度估計圖,該圖顯示收益率的分布呈現(xiàn)尖峰厚尾的現(xiàn)象. 這個檢驗問題(左側(cè)檢驗)用常規(guī)的方法難以處理,因為不知道VaR的分布,而用本研究方法卻能方便地給出結(jié)果,過程如下：

1) 利用核估計對243個對數(shù)收益率求得置信度為95%的VaR值(該值等于-3.40%)；

2) 利用Bootstrap方法,求得置信度為95%的VaR的n個(不妨取n=10 000)模擬值,并再次利用核估計給出VaR的估計分布； 然后基于該分布計算VaR的檢驗p值,給出結(jié)論. 圖1(b)給出VaR的模擬分布密度. 計算程序見附錄4,主要結(jié)果見表6,由于檢驗p值都較大,不妨接受VaR=-3.40%的假設(shè).

圖1 日對數(shù)收益率和VaR的核密度估計Fig.1 Kernel density estimations of daily log returns and VaR表6 VaR的檢驗p值Tab.6 p value of VaR

Bootstrap法經(jīng)驗分布函數(shù)法核密度法0.48700.48690.4990

6 結(jié)語

例1說明了使用Bootstrap進行檢驗p值估計的可行性； 例2說明Bootstrap方法可用來解釋一些不當?shù)淖鞣? 實例分析表明,基于Bootstrap法的檢驗p值可方便用來解決實際問題. 總的來說,使用Bootstrap方法估計檢驗p值,關(guān)鍵在于構(gòu)造合適的統(tǒng)計量,確定統(tǒng)計量的類型(連續(xù)或離散),最后通過自助樣本的經(jīng)驗分布序列得到統(tǒng)計量的近似分布,進而給出p值估計.

[1] EFRON B, TIBSHIRANI R J. An introduction to bootstrap[M]. New York: Chapman and Hall, 1993.

[2] EFRON B. Bootstrap methods: another look at the Jackknife[J]. Ann Statist, 1979, 7(1): 1-26.

[3] KULESA A, KRZYWINSKI M, BLAINEY P,etal. Points of siginificance: sampling distribution and the bootstrap[J]. Nature methods, 2015, 12(6): 477-478.

[4] NILSSON B. A compression algorithm for pre-simulated Monte Carlop-value functions: application to the ontological analysis of microarray studies[J]. 2008, 29(6): 768-772.

[5] SILVERMAN B W. Density estimation for statistics and data analysis[M]. London: Chapman and Hall, 1986.

[6] RACINEA J S, MACKINNON J G. Inference via kernel smoothing of bootstrappvalues[J]. 2006, 51(12): 5949-5957.

[7] 王星. 非參數(shù)統(tǒng)計[M]. 北京： 清華大學出版社, 2009.