函数型数据异常值检验研究范文

《统计与信息论坛杂志》2014年第六期

一、函数型数据图形描述方法综述

数据平滑方法是函数型数据图形化描述的基本工具。对于一组面板数据，当使用数据平滑方法在同一个坐标系内绘出同一个指标不同时间点的拟合曲线，而且这些曲线依照时间顺序以不同的颜色排列，就形成了类似雨后彩虹一样的图形，称这个图形为彩虹图。Hyndman使用惩罚样条插值法，对1899年到2005年间法国特定年龄段男性的死亡率曲线形成的彩虹图数据进行了平滑处理，得到一个形似彩虹的图形［7］。本文采用1991—2013年的中国上证指数统计数据，使用惩罚样条插值法对每个交易日的收盘价和总的交易金额作彩虹图，从1991年开始到2013年9月25日的数据依照不同颜色顺序排列形成了一组彩虹图。由于2007年和2010年中国股市的两次暴涨暴跌，使图形看起来没有形成一个完整的彩虹形状，但色谱同样逐次排列的光滑曲线在时间不断重复的条件下，仍形成了函数型数据的彩虹图（见图1）。同时，还可基于数据深度、数据密度或其他特性顺序作彩虹图。由于曲线在很多取值上重合，很难识别平均曲线的位置或者大多数曲线的下降位置，当异常值被曲线其他特性混淆时（例如某一段上的曲线形状与其他部分很不相同）则很难识别。对于单变量数据，通常用箱线图来解决上述问题，而本文的目的是以箱线图的形式定义函数型数据的变动，这种箱线图可以给出偏离曲线、一条中心线和一个包含曲线中间50％部分的区域。图1是以上证指数收盘价和成交金额等自然指标进行的排序，这样的排序在函数型数据的分析上缺乏新意，很难从中发现统计规律。在函数型数据研究的文献中，很多研究者采用不同的思路进行数据排序的尝试，如可利用稳健主成分方法计算出每个时间点的主成分得分，并以此排序形成彩虹图。几乎所有针对函数型数据的作图方法都涉及到对函数型数据的排序。图1中的数据是基于时间排序，然而对于很多数据集以数据本身潜在的价值来排序，似乎是更好的选择。关于函数型数据的很多排序方法都使用了数据深度或者数据密度的概念，这些方法揭示了给定的函数型数据的观测值及其潜在分布的“深度”或“密度”的情形［8－9］。一般情况下，一个深度函数或密度函数的轮廓图可以用来展现多变量数据的可以看得见的形状和结构特征。

（一）函数型数据深度的测度方法Febrero等人提出了一种基于函数型深度概念的异常值检测方法，其基本定义为：其中对于给定的x值，D（yi（x））是对其深度的度量函数。在这个定义下，通过一个不断增加的o｛i｝序列来定义曲线的顺序，因此接近x轴的第1条曲线的函数型数据深度最小，而最后一条的最大。

（二）二元主成分得分深度的测度对多变量函数型数据进行主成分分析，设φk（x｛｝）为主成分，zi｛，k｝是对函数型数据进行主成分分解后的主成分得分。原始数据yi（x｛｝）中的大量信息可以由少数的前几个主成分及其得分反映出来。大多数情形下，对于一些经济或自然观察数据而言，少数几个主成分得分往往可以积累超过80％的方差贡献率，这也正是主成分分析的优点之一。因此，将考虑前两个得分向量（z1，1，z2，1，…，zn，1）和（z1，2，z2，1…，zn，2），并考虑将这两个向量应用到深度函数的方法中，还可将二维平面上的点（zi，1，zi，2）看作zi。Tukey还提出了二元得分可以利用半空间位置深度排序的方法，即用d（θ，z）表示，θ∈R2对应于二元数据区域z＝｛zi；i＝1，2，…，n｝［8］。Tukey深度函数被定义为：当θ存在于封闭半平面边界上时，该半平面内全部数据点的最小值可以按照距离d（zi，Z）以升序排列，这种顺序下的第一条曲线可认为是平均线，而最后一条曲线是在样本曲线中离中心最远的曲线。

（三）数据密度方法Scott提出的数据密度方法是通过每个观测值上的二元核密度估计值排序［10］209－210。设oi＝＾f（zi），＾f（zi）是由所有二元主成分得分计算得到的核密度估计值，这样函数型数据就可以按照o｛i｝的值以升序排序。因此，有最高密度的曲线是第一个观测值，而最后一条曲线是最低密度值；第一条曲线被认为是模板曲线（这里可以理解为基本的参照曲线），而最后一条曲线被认为是最不同寻常的曲线，实际上也最可能是异常值。应注意到，这种排序下的最后一条曲线取值与其他曲线差异可能不大，其二元得分也可能并不在散点图（zi，1，zi，2）的边缘，可能的情形是有一点在散点图内，但该点附近再无其他点，这样该曲线就表现为低密度值。

二、函数型数据图形分析的主要方法

（一）彩虹图对于一些不依时间而按其他统计指标排序的数据，基于特定的排序指标或辅助标志也可以用彩虹图来表示，例如上面定义的数据深度或数据密度排序指标，绘图时根据o｛i｝的排序即可选择对应的线条颜色。为验证上述方法，Hyndman根据国际气象组织公布的厄尔尼诺现象的测量数据进行模拟计算，选择了1951年1月至2007年12月南太平洋赤道附近的厄尔尼诺浮标点（南纬0～10度，西经90～80度）海平面的月平均温度时间序列数据进行了分析［11］。本文更新了上述数据，采用1982年1月至2013年12月的数据进行了再次模拟，并绘制了一组彩虹图，这些数据没有显著的时间趋势，因此基于时序的彩虹图意义不大。图2显示的是基于深度函数和密度函数指标顺序所表示的彩虹图，颜色按照彩虹颜色顺序，最接近中心数据的曲线设为红色，远离中心的设为紫色。按照深度和密度描绘曲线，因此红色曲线是最模糊的，而紫色最清楚，即使曲线部分大多数数据重合。图2中横坐标表示数据测量的月度顺序，纵坐标表示海平面的温度；图2（a）中黑色实线表示中线，图2（b）中黑色实线表示参照曲线。

（二）打包图函数型数据的打包图是建立在二元主成分打包图基础上，并应用前两个主成分得分而得到。图3使用了Tukey所定义的位置深度函数的概念，这个深度区域Dk是所有θ的集合，d（θ，k）≥k，因此深度区域形成一个凸面体，对于任意k2＞k1有Dk1Dk2［8］。Tukey的二元深度中位线被定义为：给定一个独立的θ，使d（θ，k）达到最小的θ值，并被定义为最深区域的重心。类似于二维箱线图，二维打包图有一个中心点（即Tukey中位数）、一个内部区域（包）和一个外部区域，除此之外异常值以个别点显示。包被定义为至少包括总数50％的观测值的最小深度区域，打包图的外部区域（或称“围栏”），是包括由显著性水平决定的概率度因子ρ得到的包（与Tukey中位数有关）所组成域内的点组成的凸面体，当ρ＝2．58时，假设投影的二维得分服从标准正态分布，则该区域允许外部区域的围栏部分包含99％的观测值；同理，当ρ＝1．96时，围栏部分包含的数据观测值可以达到95％；当ρ＝1．64时，围栏部分包含的数据观测值可以达到90％。图3显示的是关于厄尔尼诺现象数据的描述图形。在图3（a）中，深灰色区域表示的是50％的包，浅灰色区域显示的是99％的边界，这些凸面体对应于图3（b）中函数型打包图的相同阴影部分。位于围栏区域之外的点一般可视之为异常值，不同颜色的异常值表明右边的个体函数曲线同左边的二元主成分得分相匹配。图3（a）中，红色星号表明了二元主成分得分的Tukey中位数，图3（b）每个平面的黑色实线表明的是中位数曲线，蓝色点线是中位线上对应点的95％的置信区间，深灰和浅灰域标明了包和边界。红色星号是Tukey深度平均数，右面的黑色实线是平均线，上下虚线是对应点的95％置信区间，而外域之外的曲线用不同颜色标明的是异常值。在厄尔尼诺数据中检测出的异常值出现在1982—1983年和1997—1998年。1982—1983年间厄尔尼诺指数在1982年6月升高较慢，在1982年9月至1983年6间海洋表面温度有一个极不正常的升高；1997—1998年间厄尔尼诺指数也不正常，尤其是在3月和5月。国际气象组织和有关国家公布的资料表明，1982年4月至1983年7月的ENSO现象是几个世纪来最严重的一次，太平洋东部至中部水面温度比正常高出约4～5℃，造成了全世界1300～1500人丧生，经济损失近百亿美元。2009年12月在丹麦举行的哥本哈根联合国气候变化大会上，2009年便被定为厄尔尼诺年。当异常值远离平均数时函数型打包图可能是一种较好的异常值检测方法，然而当异常值接近于平均数时，这种深度测量异常值的方法可能误判异值，在这种情况下函数型HDR箱线图则更适用。

（三）HDR箱线图HDR箱线图即高密度区域箱线图，是Hyndman最先提出并进行实证分析的。函数型HDR箱线图是建立在二维HDR箱线图基础上的，由前两个主成分得分并进行分析后绘制。二维箱线图HDR是由一个二维核密度估计＾f（z）构造，定义如下：区域，在这个区域内所有点要比区域外任何点有更高的密度估计，因此可以把这个区域命名为“高密度区域”。对于二维密度估计，HDR可看作是一种类似地理学上的等高线，随着α的减小，区域中包含的范围在不断扩大。二维HDR箱线图表示的是众数，被定义为sup＾f（z），即核心部分包含50％和外层包含99％的点的最高密度区域，分布在99％区域之外的点即可视为异常值。函数型HDR箱线图是一个由前两个稳健主成分得分的二元HDR箱线图到函数的映射，包括参照曲线（有最高密度的曲线）、内部和外部区域。内部区域是由50％二元HDR内的点所形成的曲线而组成的域，因此有50％的曲线在内部区域。类似地，外部区域被定义为二元HDR图上的点所形成的曲线而组成的区域（见图4）。图4显示了用厄尔尼诺统计数据进行实证分析的例子。在图4（a）中，深灰色和浅灰色区域分别表示50％的内部区域和99％的外部区域，直接对应于图4（b）中函数型HDR箱线图中的同色阴影区域，而外部区域之外的点被认为是异常值，不同颜色的异常值与图4（b）中个体曲线颜色及图4（a）中二维得分HDR相对应。图4（a）中的红色圆点表示二元主成分得分的众数，对应于图4（b）平面图中的黑色实线。从图4不难看出，1982－1983年和1997－1998年作为异常值被排除在外部区域之外，可以显著地定义为异常值，这个结论和前文的讨论是一致的，符合实际情况。包括打包图和HDR箱线图在内的任何一种异常值检测方法，都需要提前确定外部区域的覆盖率。在99％的概率保证程度下，厄尔尼诺数据集中检测出的异常值出现在1997—1998年。如果假定厄尔尼诺数据中的覆盖率分别为92％和93％，那么在每个例子中检测出的异常值将同打包图得到的结果相对应，这表明相对于其他数据，这些异常值有不同的数量大小和分布形状。

三、借助图形分析方法进行异常值检测

在函数型数据分析中，一般利用函数型数据打包图和函数型深度方法即可完成异常值的检测。

（一）函数型深度方法Febrero提出了一个对每一条曲线yi（x｛｝）计算其似然比检验统计量的异常值检测方法［6］。如果检验统计量的最大值比给定的置信值C大很多，那么可以断定该数据点为异常值。去掉该点，对余下的数据用该方法继续检测其他异常值，一直重复这个过程，直到不再有其他异常值出现。这种检测方法是建立在式（1）给定的函数型数据的深度测量基础上，故对形状异常值并不敏感。

（二）误差平方积分方法Hyndman等人提出了一种基于稳健型函数主成分分析的异常值检测方法［4］。设对每一个随机观测样本点i的误差平方项积分如下：

（三）稳健的马氏距离方法稳健的马氏距离方法是公认的可以用在离散型曲线yi（xj）；j＝1，2，…，｛p｝上的一种多元异常值检测方法。假定在相同的空间密度坐标x1，x2，…，x｛p｝上观测到函数型数据，那么平方的稳健的马氏距离可以定义如下：

四、异常值检测方法的实证分析

根据上述讨论，针对函数型数据的异常值数据，再次采用公开的厄尔尼诺现象在1982年1月到2014年1月的同步统计数据进行实证分析（图5）。

（一）正态分布假定下的异常值检测由图5可以直观地看出，在1982—2014年的33个年度数据中，除了4月份以外，其他月份的箱线图都有异常值出现，而且这些异常值主要集中在1982—1983年和1997—1998年两个厄尔尼诺现象严重的年份，这与前文所述的情形基本类似。为了进一步比较本文所述的异常值检测方法，对于同一组数据，首先使用传统的异常值检测方法进行初步分析，即考虑大样本情形下以均值为核心构建一个半径为3s的置信区间，如果数据溢出这个区间，则判定为异常值，结果见表1。

（二）函数型数据异常值检测对于厄尔尼诺数据，同样利用前文涉及到的三种方法和HDR箱线图与打包图进行异常值检测。R语言程序准确地记录到了每种方法的计算时间，考虑到全部计算程序在同一台计算机上完成，硬件的影响可以忽略不计。比较一种异常值检测方法的优劣，主要考察该方法的精确度和计算速度。以传统的箱线图方法为基准计算出的异常值检出率，即检测出的异常点占全部异常点的比例，除了函数型深度方法以外其他方法的检出率均为100％。检测方法的响应速度则直接采用统计软件显示的有效计算时间，计量单位是秒。计算表明，图形方法是异常值检测的较好方法，不但保持了较高的检出率，而且计算时间也较短，比较适合大数据场合下持续进行的计算和分析，其中打包图比HDR箱线图的计算效率更优。如果考虑到未来在线大数据的计算，动态的函数型HDR箱线图和打包图都是一种较好的选择方法，甚至可考虑生成类似股价图一类的动态图示方法，从而及时地检出异常值，排除异常原因，最终保证数据过程的正常运行（见表2）。

五、结论与进展

本文针对函数型数据的可视化分析，提出了三种描述工具和三种识别函数型数据异常值的方法。在所熟悉的二维空间中，对主成分得分按照数据深度和数据密度排列后，异常值和正常点自然就分开了。本文所提及的方法，其优点是可以对复杂的函数型数据进行快速分析，并通过图形直观地表示出来，有利于批量地处理大数据，也有利于在更广义的领域推广使用，有较高的应用价值。根据笔者对更新的厄尔尼诺统计数据的实证研究，所介绍的三种异常值检验方法，无论是检测速度、效率和直观性还是检测的精度，都比以往的传统方法更优。国外近期文献显示，在复杂数据背景下无论是自然科学还是社会科学，都对函数型数据的方法开展了研究，其研究主要包括三方面的进展：一是对传统主成分分析方法的扩展和进一步探索，主要目的是针对函数型数据的降维分析，在尽可能保持数据信息的前提下，降低数据的复杂性，以保留最多的因子贡献及方差贡献；二是探索更多的检测异常值的方法，以Tukey的深度函数和密度函数为基础，逐步发展了一些快速捕捉数据特性和检测函数型数据中异常点的方法，能辅助计算方法的改进，以替代传统的异常点检测方法；三是发展以稳健方法为主要趋势的基本统计方法，以函数型数据的众数和中位数为参照，对各类数据集进行排序和分割，以得到更直观有效的结果。

作者：米子川单位：山西财经大学统计学院