1. 引言
函数型数据最早是1982年Ramsay [1],1991年Ramsay和Dalzell [2]提出的,其为定义在某个集合T上的函数,表示为
。例如随时间变化的股票交易价格,随位置变化的温度等。假设检验作为两大统计推断问题之一,历来受到研究学者的极大关注。近年来,函数型数据中的假设检验问题也受到了关注。例如1998年Fan [3]研究了函数型数据的显著水平;针对函数型数据,2012年Ferraty [4]对其条件分布,Bugni [5]对缺失数据进行了检验;针对协方差算子,2013年Fremdt [6]对其是否相等以及Jaruvskova [7]对其是否存在变点进行了检验;针对两样本均值函数相等性的检验问题,2009年Horváth [8]等利用FPCA给出了协方差算子特征根和特征函数的估计,然后基于此估计构造了渐近分布和渐近混合分布两个检验统计量。2010年Zhang [9]等研究了FDA中两样本Behrens-Fisher检验问题,对两个高斯过程均值函数相等性进行了检验,提出了混合卡方渐近统计量。2014年Portugue [10]提出的针对响应变量是实值的函数线性模型的拟合优度检验;2020年Hu [11]等对主成分结构不明显的函数数据,构造了一个4阶的U统计量,并证明了该检验统计量是渐近正态的。
长期以来,样本量与统计功效总是相伴出现的。统计功效由显著性水平、效应量和样本量三者共同决定。在显著性水平和效应量不变的情况下,样本量越大,功效就越高。所以,统计功效实际上取决于样本量的大小,而提高统计功效的主要方法是增加样本量。然而如果样本量太小,研究结果的可重复性及代表性就欠佳可能得出错误的结论。但是如果样本量过大,试验过程中对研究对象的潜在伤害越大,比如一些动物实验还可能存在一定的伦理问题,而且研究所需的经费和资源就越多,项目执行的难度也越大。样本量在实验结果、动物伦理和经济效应方面具有的突出、核心和重要作用,这要求我们在假设检验中必须提前给出合适、准确的样本量。
本文余下内容结构如下:第2节介绍预备知识,提出检验统计量,证明其理论性质。第3节对检验统计量在不同功效下做Monte Carlo模拟计算最低样本量。第4节验证该样本量确定方法在气象数据的应用。第5节总结全文。
2. 相依函数型数据均值检验
2.1. 可逼近的函数时间序列
引理1 我们假设
是;L2-m-可逼近的函数时间序列,因此也是平稳的,并且
满足在
中
,以及
有界即
,那么。
其中
是一个高斯过程,
,
(1)
核函数
是
上的平方可积函数。
不失为一般性,我们假设
是非零的均值函数,
(2)
其中随机误差函数
满足引理1的假设。
核函数
满足以下条件:
,
是连续有界函数,
且如果
,那么对于某些
(3)
接下来我们定义经验协方差函数,
其中
。
的估计量为
(4)
其中
是满足以下条件的平滑窗宽:
当
时,
(5)
引理2 假设函数型时间序列满足(1)中的模型,在条件(3) (5)下
其中
是(1)式中的定义,
是(4)式中定义。
引理1和引理2的证明Horváth,Kokoszka,Reeder [12]已给出。这里我们直接引用结论。
2.2. 均值函数的检验
我们知道两个样本假设检验过程中,样本量的大小一致时检验效果会更好,同时样本量数量相同也符合很多情况下我们设计实验的习惯。所以在该问题下我们第一个假设就是训练样本和测试样本的样本量数量相同。
我们考虑两个样本数量相同的训练样本
和测试样本
。我们假设这两个样本分别来自以下模型:
其中
。
我们假设随机误差函数
和
都满足所提出的条件,特别注意的是他们的长期协方差核定义如下:
.
我们假设随机误差函数
和
是相互独立的,对于这个假设检验问题我们得:
我们知道在
空间里,原假设
意味着
,备择假设
意味着
。因为统计推断是关于样本观测值的平均函数,并且我们知道
和
分别是
和
的无偏估计量,因此如果
很大,我们就有足够的我们有充分的理由拒绝原假设。
定理1:如果
,以及引理1和引理2均成立,并且当
时,我们有:
其中,
是一个均值为零的高斯过程并且具有协方差
证明:为了证明以下定理我们先说明一个事实:当
我们有:
(6)
其中
,这样的事实我们在经验协方差算子和总体协方差算子的特征值和特征函数的估计里已经广泛应用。
引理1和引理2表明,当
时我们有
,
,
其中
和
是
上相互独立的均值为零的高斯过程并且它们的协方差分别为
,
,
并且当
成立时我们知道还有
(7)
其中
。定理1得证。
定理2:如果
和定理1成立,,特别地我们还知道
。
证明:通过希尔伯特空间中的Ergodic Theorem我们知道,
以及
。定理2得证。
定理1中的核函数
所定义出来的协方差算子
,记
中的特征值满足
,通过Karhunen-Loève展开式(K-L展开),我们可以得到
其中
是相互独立的标准正态随机变量。
由于特征值
未知,所以无法直接用上述公式来模拟
的分布,下面我们来解释如何估计
。
假设
是
在
中的一致估计量,当
,
(8)
我们设想满足上式关系式的估计量
的构造。首先对于估计量
不管
和
是否成立上式都成立,其次估计量
也不依赖于
和
。并且我们还知道(6)式在
和
下都成立,因此我们可以在备择假设
估计
的分布。
假设是
的特征值,我们有
其中,
是相应的特征函数且。然后选择合适的
使得累计方差比例
足够大,我们可以用
的分布来近似
的分布。
我们的统计推断是基于
的差值,当
,
因此,
是
的差值的渐近协方差核。我们使用
在
上的投影,定义投影为
从而得到向量
。
我们还能从定理3的证明过程知道:
(9)
其中
表示均值为0的
元正态随机变量,并且协方差矩阵
由于协方差算子
未知,我们无法直接计算出
。但是我们知道满足协方差算子
的估计量也可以用来计算
的估计量。我们假设
是上式中定义的特征函数且
根据上式我们提出下面这个检验统计量
定理3 假设
,定理1和(8)式成立且
,那么
。其中
是一个自由度为
的卡方随机变量。
证明:在定理3的假设下,我们有
其中
和
是相互独立的
元正态随机变量。
因为
其中
。
我们知道矩阵
满足
其中
是一个克罗内克函数,
是正交的特征函数,
是相应的特征值。
由(6)式我们有
再由(7)式我们有
再由(9)式以及矩阵Q的对角性。定理3得证。
参考Jacob Cohen [13]效应量的计算方式,
从而就能计算特定显著性水平
下要检测到特定效应量Effect size和检验功效power的最低样本量。同样地,我们规定:0.2或更小的w是较小的效应量;0.5附近的w是中等效应大小;0.8或更大的w是较大的效应量。
3. 数值模拟
为了计算最低样本量,本节将使用蒙特卡洛模拟的方法进行数值模拟。由定理3可知检验统计量在原假设下渐近服从
,其中d是所选取的主成分的个数。主成分个数d选取是函数型数据主成分分析经常遇到的问题。常见的方法有交叉验证法、碎石图方法、pseudo-AIC、以及累计方差占比(CPV)方法,本文采用的是累计方差占比(CPV)方法,即
选取合适的
使得
不小于95%。
Horváth,Kokoszka,Reeder [12]通过他们提出的方法进行数值模拟时只研究了检验水平(size)和检验功效(power),而没有给出高功效下的最低样本量表。参考Horváth et al.的数据生成方式,不失为一般性,设定均值函数
时,原假设成立,备择假设下,设定
以及
。随机误差函数
和
都是独立同分布的布朗桥。
核函数选择的是Politis和Romano [14]的flap top核,即
其中窗宽为
。
在
区间上等间距取301个观测点,c从0.5变化到1,在显著性水平
为0.01、0.05、0.10三种情形下做1000次重复模拟结果如表1所示。本文所有程序都是用R完成的。
Table 1. Sample size and empirical power of mean function equality test for dependency functional data
表1. 相依性均值函数相等检验的样本量和经验功效
|
Power = 0.8 |
|
Power = 0.9 |
|
α = 0.01 |
α = 0.05 |
α = 0.10 |
|
α = 0.01 |
α = 0.05 |
α = 0.10 |
c |
n |
EP |
n |
EP |
n |
EP |
n |
EP |
n |
EP |
n |
EP |
0.5 |
320 |
0.801 |
270 |
0.801 |
228 |
0.802 |
413 |
0.906 |
323 |
0.901 |
302 |
0.900 |
0.6 |
227 |
0.800 |
182 |
0.801 |
160 |
0.804 |
301 |
0.914 |
232 |
0.903 |
215 |
0.918 |
0.7 |
165 |
0.808 |
118 |
0.803 |
106 |
0.809 |
213 |
0.915 |
170 |
0.903 |
159 |
0.900 |
0.8 |
131 |
0.803 |
102 |
0.806 |
93 |
0.803 |
173 |
0.903 |
135 |
0.901 |
113 |
0.902 |
0.9 |
106 |
0.803 |
85 |
0.801 |
75 |
0.806 |
132 |
0.901 |
109 |
0.901 |
95 |
0.907 |
1.0 |
85 |
0.809 |
66 |
0.803 |
58 |
0.800 |
102 |
0.900 |
87 |
0.912 |
75 |
0.903 |
4. 气象数据实例分析
改革开放以来,随着我国工业化和城市化进程的不断加速,大气污染问题日益凸显。2013年我国遭遇了有观测以来最严重的雾霾天气。空气质量问题是关乎国家生态安全和人民福祉的重大课题。因此持续关注和研究空气质量数据有重要意义。空气质量指数(AQI)和六大空气主要污染物:细颗粒物(PM2.5)、可吸入颗粒物(PM10)、二氧化硫(SO2)、二氧化氮(NO2)、臭氧(O3)、一氧化碳(CO)不仅是天然的函数型数据,另一个重要特征是这些数据还是相依型函数型数据。以PM2.5浓度为例,可以看出本文选取的数据具有明显的函数曲线特征,因此可以用函数型数据中的平滑方法对数据进行处理。同时,这些曲线并不是具有相同分布的独立曲线,它们不是简单的随机曲线样本。一方面,某一天的PM2.5浓度曲线会受到前几天PM2.5浓度曲线的影响,曲线之间具有记忆性;另一方面,某一时刻的PM2.5浓度也总受到前面时刻PM2.5浓度的影响。因此,本文选取的样本符合假设条件,可以进行实例分析。
为了验证我们提出的均值函数假设检验问题的样本量确定方法的有效性,本章将选取北京市东城区和西城区空气质量历史数据进行研究。此时我们关心的是自2020年9月中国明确提出实现“双碳”目标以来,北京市东城区和西城区的空气质量指数(AQI)和六大空气主要污染物PM2.5、PM10、SO2、NO2、O3、CO的浓度是否相同。数据均来自北京市环境保护检测中心。
因为检验的思路基本一致,所以我们以空气质量指数(AQI)为例。首先,对采集的每天24个空气质量指数(AQI)数据补齐缺失值,利用缺失时刻前后三小时的平均值代替缺失值。在拟合之间我们首先压缩一下数据范围,我们将空气质量指数(AQI)数据对数化处理,然后通过数据变换,将每天时间
变换到
上面。然后对于这些离散的观测值我们再使用50个傅里叶基光滑拟合。对2024年3月1日至3月30日北京市东城区和西城区的空气质量指数AQI、六大污染物PM2.5、PM10、SO2、NO2、O3、CO的浓度数据进行光滑处理如图1所示。
Figure 1. AQI、PM2.5、PM10、SO2、NO2、O3、CO concentration curves smoothed in Dongcheng District and Xicheng District
图1. 东城区和西城区平滑处理后的AQI、PM2.5、PM10、SO2、NO2、O3、CO浓度曲线
使用本文提出的方法估计出特征值和特征函数再按照CPV ≥ 95%,均选取5个主成分,计算统计量
,利用公式计算效应量大小
设置显著性水平
,功效power = 0.9,利用R语言的pwr包我们计算得到样本量大小为35。选择北京市东城区和西城区2024年3月1日至4月4日这35天的空气质量指数AQI计算本文提出检验的p值,p值等于2.60981e−07。因此,显然有充分的理由拒绝原假设。同样的方法运用于东城区和西城区六大空气主要污染物PM2.5、PM10、SO2、NO2、O3、CO的浓度是否相同检验问题时,计算得到所需的样本量分别为:53、28、29、29、25、56。选择相应天数进行检验,计算得到p值均接近零,如表2所示。这意味着该样本量的确定方法是有效的。
Table 2. Sample size and p-value for mean test of air quality data
表2. 空气质量数据均值检验的样本量及p-value
空气质量数据 |
Effect Size |
样本量大小 |
p-value |
AQI |
0.788989 |
35 |
2.60981e−07 |
PM2.5 |
0.64600 |
53 |
0.00066 |
PM10 |
0.89394 |
28 |
0 |
SO2 |
0.87843 |
29 |
0 |
NO2 |
0.86958 |
29 |
0 |
O3 |
0.93815 |
25 |
0 |
CO |
0.62870 |
56 |
0 |
5. 结论
本文针对相依性函数型数据,假设训练样本和测试样本在样本数量相同条件下提出检验统计量,并对该检验统计量的理论性质进行研究,得到了原假设下的渐近分布,并提出了效应量的度量方式。在显著性水平
、检验功效power、效应量Effect size确定的情况下计算最低样本量。最后将该方法应用到空气质量指数(AQI)和六大空气主要污染物PM2.5、PM10、SO2、NO2、O3、CO的浓度数据证明了该方法的有效性。