基于Bootstrap思想的非参数检验p值的估计方法

吕书龙, 刘文丽

(福州大学数学与计算机科学学院,福建 福州 350116)

0 前言

参数型假设检验特别是基于正态总体的假设检验,无论是从理论上还是实践上都已经成为经典. 参数型假设检验的前提要求总体服从的分布形式或分布族是已知的,这对构造检验统计量并确定其分布和精确计算检验p值都提供了最有利的条件. 但在实际问题中,总体的分布通常是未知的,于是,人们就希望在总体分布不假定的前提下,基于数据本身来构造反映总体特征及问题的统计量或信息,从而达到统计推断的目的,这也是非参数统计的宗旨. 众所周知,采用检验p值完成假设检验的推断已成为一种最基本的模式. 作为非参数统计的一个主要分支,非参数检验在医学、 生物学、 信息学、 金融、 管理、 教育等众多领域应用广泛.

如何构造检验统计量并计算检验p值,是非参数检验的一个关键问题,特别是在样本量较小的情况下,这个问题尤为突出. 但本研究注意到由于非参数方法不假定总体分布,对原始数据的使用也不像参数型方法,这可能会导致非参数方法的有效性低于参数方法. 而通过重抽样Bootstrap方法[1-2]可以很好地给出总体特征的估计及精度表示[3]、 区间估计以及假设检验等,甚至能发现常规检验方法中可能存在的不稳定问题. 检验p值的估计是非参数假设检验的基本问题,Nilsson[4]提出检验p值的数据量缩减型Monte-Carlo非参数估计方法,并给出其在生物信息学上的应用. Silverman[5]探讨了基于Bootstrap的渐近枢轴统计量的经验分布函数的核方法及检验p值的估计. 本研究尝试探讨宽泛的检验p值估计方法及其适用范围.

1 问题描述

设总体X～F(x),F(x)未知,X1,X2, …,Xn为其独立同分布样本,x1,x2, …,xn为样本观测值. 为推断与总体相关的某些问题,通常要构造样本函数T=T(X1,X2, …,Xn),即统计量. 但由于F(x)未知,此时统计量T的研究基本上可归入非参数领域.

假设检验应用广泛,很多的实际问题都可以转化成假设检验而得到较好解决. 在作假设检验时,是否拒绝原假设H0,基本上都是通过先假定H0成立,再计算其检验p值并根据检验p值的大小来进行推断. 各种统计软件仅根据统计量T的分布给出检验p值,但是否拒绝原假设全部交给使用者自行裁定. 如果统计量T的分布难以确定,就会导致检验p值的计算变得无章可循,实际问题也就难以得到有效的解决.

本研究在总体分布未知且小样本条件下,避开统计量的理论分布探讨,单纯从模拟仿真的角度,通过Bootstrap法产生的自助样本,构造“检验统计量”并结合随机模拟、 经验分布和核密度估计[6]等方法实现原假设下的检验p值的估计,从而在一定程度上解决了上述假定条件下假设检验问题的推断.

2 主要理论和工具

1) 格列汶科(Glivenko)定理给出用样本的经验分布来估计总体的分布函数的理论,即若设总体X的分布函数为F(x),经验分布函数为Fn(x),则有

(1)

3) 设K(x)为定义在(-∞, +∞)上的一个Borel可测函数,hn>0为常数,则称

(2)

为总体密度函数f(x)的一个核估计. 其中, K(x)称为核函数； hn称为窗宽. 利用核密度函数估计法可对T(X1, X2, …, Xn; Fn)的密度进行估计.

3 检验p值估计的基本思路和算法

先确定假设检验问题和类型,构造检验统计量； 然后利用Bootstrap重抽样得到检验统计量的模拟抽样序列； 在抽样序列的基础上,利用Bootstrap法、 经验分布法和核密度估计法计算不同检验问题和类型对应的检验p值.

3.1 三类检验问题

假定检验统计量为T(X1,X2, …,Xn;Fn),其一次试验的值为Tx(x1,x2, …,xn),检验p值的直观意义表示为发生比Tx(x1,x2, …,xn)更极端事件的概率. 对于三类假设检验问题,下面给出检验p值的形式表示.

1) 双侧检验,p1=P(T≥Tx(x1,x2, …,xn)),如果p1<=0.5,则检验p值=2p1,否则检验p值=2(1-p1)；

2) 右侧检验,检验p值=P(T≥Tx(x1,x2, …,xn))；

3) 左侧检验,检验p值=P(T≤Tx(x1,x2, …,xn)).

3.2 检验p值的模拟估计

3.2.1Bootstrap法估计检验p值

设T∈(Ti,Ti+1),则

(3)

3) 基于Bootstrap的百分位数法,给出三类检验问题的检验p值公式：

(4)

3.2.2 经验分布法估计检验p值

(5)

通过经验分布函数法得到的检验p值本质上和Bootstrap法得到的检验p值是一致的,只是表达方式不同而已.

3.2.3 核密度法估计检验p值

p双侧=2p1(p1≤0.5)或2p2(p2≤0.5)；p右侧=p2；p左侧=p1

(6)

4 模拟分析与R实现

为了更好地说明本研究提出的算法,下面给出正态型参数检验和非参数检验的两个例子,并与常规的假设检验进行比较.

例1 设总体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, 试给出关于E(X)=100的三类检验问题的检验p值.

表1 常规算法与本研究算法的检验p值Tab.1 p value of conventional algorithm and our algorithm

如果研究的问题是： 设Y=X2,求关于E(Y)=10 016的三种假设检验问题的检验p值. 显然使用常规的方法难以计算检验p值,而本研究算法却可以很轻松地估计出检验p值,只需将统计量函数定义为

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

例2 对于中风病人与健康成人血液中尿酸浓度数据,见表2[7],问两类人血液中的尿酸浓度的变异是否存在显著差异?

表2 中风病人与健康成人血液中尿酸浓度数据

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

本例的假设检验问题可表示为

(7)

文献[7]采用Moses检验法,并给出拒绝H0的推断. Moses检验的作法是将两样本分组,保证每组长度都一样,如各分成m1,m2组； 然后计算两样本各小组的离差平方和SSAi, SSBj,i=1, 2, …,m1；j=1, 2, …,m2； 再将SSAi, SSBj合并后计算样本1的m1个组的秩和S,从而构造Moses统计量TM

(8)

王星[7]认为如果两组数据的方差存在很大差异,从平均看,其中一组的离差平方和比另一组的离差平方和要小. 文献[7]以每组3个样本点依据前后顺序将样本1分成4组,样本2分成3组且弃用最后1个样本点,并给出Mann-Whitney的Wα表法进行查表检验. 文献[7]介绍了这种方法,并给出具体的计算过程. 本研究对其中的三个细节问题作深入探讨,一是如何分组,二是不正好分组时弃用哪些样本点,三是检验的结论是否会随分组不同而不同. 本研究认为引入随机分组和随机弃用,结合大量模拟来进行检验更符合实际情况. 为了说明上述问题,设计R程序(见附录2)并整理运行结果,如表3所示.

表3 文献[7]方法的某1 000次随机模拟的统计结果Tab.3 Statistical results by 1 000 stochastic simulations based on the method in literature [7]

若按照文献[7]提供的方法,得到显著性水平为0.05时的临界值0和12,上述程序模拟1 000次,合计出现了383次接受原假设的情况(如表3所示),且经过多次模拟表明,出现上述检验错误的频率大致都在38%左右. 以上说明不同的随机分组和取点会造成检验结论的不同,即上述三个细节问题必须慎重考虑. 依据文献[7]构造的随机模拟表明该方法可能呈现出相当大的不稳定性.

为此有必要充分考虑上述三个细节问题,并探讨更加稳定的算法来解决上述方差相等与否的检验问题. 对数据先作一下基本分析,计算其样本均值、 样本方差和样本极差,直观判断两组数据的差异性.

表4为基本统计量,从表4看出这两组数据的各项统计量指标都存在较大的差异,检验的结论应该拒绝原假设才合理. 为了避免出现上述随机分组、 取点的干扰,引入刀切法思想,即构造原样本的舍一样本方差：

； i=1, 2, …, n

(9)

(10)

若式(7)中原假设成立,则上述两个TM值既不太大也不太小,否则可拒绝原假设,程序见附录3.

表5为检验p值的比较,从表5看出修正后的Moses检验的检验p值远小于文献[7]方法得到的检验p值,相对于文献[7]而言有更强的理由拒绝原假设.

表4 基本统计量Tab.4 Basic statistics

表5 检验p值的比较Tab.5 Comparison of p values

5 实例分析

风险价值方法 (value at risk, VaR)作为金融风险度量的一种重要方法,在金融风险管理中应用广泛. 金融资产价格变化的分布经常出现尖峰厚尾的现象,若使用模型设定法(如价格变化服从t分布或正态分布等),估计往往不精确甚至不符合实际. 本研究以2016年上证指数前一日收盘价与后一日开盘价(共计244个数据)的日对数收益率rxt=log(开盘价t)-log(收盘价t-1),t=2, 3, …, 244(合计243个数据)为例,构造这样一个风险检验问题： 在置信度为95%时,H0:VaR=-3.40%,即以95%的置信度认为日对数收益率不低于-3.40%. 图1(a)给出日对数收益率的核密度估计图,该图显示收益率的分布呈现尖峰厚尾的现象. 这个检验问题(左侧检验)用常规的方法难以处理,因为不知道VaR的分布,而用本研究方法却能方便地给出结果,过程如下：

1) 利用核估计对243个对数收益率求得置信度为95%的VaR值(该值等于-3.40%)；

2) 利用Bootstrap方法,求得置信度为95%的VaR的n个(不妨取n=10 000)模拟值,并再次利用核估计给出VaR的估计分布； 然后基于该分布计算VaR的检验p值,给出结论. 图1(b)给出VaR的模拟分布密度. 计算程序见附录4,主要结果见表6,由于检验p值都较大,不妨接受VaR=-3.40%的假设.

图1 日对数收益率和VaR的核密度估计Fig.1 Kernel density estimations of daily log returns and VaR表6 VaR的检验p值Tab.6 p value of VaR

Bootstrap法经验分布函数法核密度法0.48700.48690.4990

6 结语

例1说明了使用Bootstrap进行检验p值估计的可行性； 例2说明Bootstrap方法可用来解释一些不当的作法. 实例分析表明,基于Bootstrap法的检验p值可方便用来解决实际问题. 总的来说,使用Bootstrap方法估计检验p值,关键在于构造合适的统计量,确定统计量的类型(连续或离散),最后通过自助样本的经验分布序列得到统计量的近似分布,进而给出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] 王星. 非参数统计[M]. 北京： 清华大学出版社, 2009.