|
|
本帖最后由 linccfzu 于 2017-12-20 11:29 编辑
格拉布斯法—异常值判断
▲概述：一组测量数据中，如果个别数据偏离平均值很远，那么这个(这些)数据称作“可疑值”。如果用统计方法—例如格拉布斯(Grubbs)法判断，能将“可疑值”从此组测量数据中剔除而不参与平均值的计算，那么该“可疑值”就称作“异常值(粗大误差)”。本文就是介绍如何用格拉布斯法判断“可疑值”是否为“异常值”。
▲测量数据：例如测量10次(n＝10)，获得以下数据：8.2、5.4、14.0、7.3、4.7、9.0、6.5、10.1、7.7、6.0。
▲排列数据：将上述测量数据按从小到大的顺序排列，得到4.7、5.4、6.0、6.5、7.3、7.7、8.2、9.0、10.1、14.0。可以肯定，可疑值不是最小值就是最大值。
▲计算平均值x-和标准差s：x-＝7.89；标准差s＝2.704。计算时，必须将所有10个数据全部包含在内。
▲计算偏离值：平均值与最小值之差为7.89－4.7＝3.19；最大值与平均值之差为14.0－7.89＝6.11。
▲确定一个可疑值：比较起来，最大值与平均值之差6.11大于平均值与最小值之差3.19，因此认为最大值14.0是可疑值。
▲计算Gi值：Gi＝(xi－x- )/s；其中i是可疑值的排列序号
——10号；因此G10＝( x10－x- )/s＝(14.0－7.89)/2.704＝2.260。由于 x10－x-是残差，而s是标准差，因而可认为G10是残差与标准差的比值。下面要把计算值Gi与格拉布斯表给出的临界值GP(n)比较，如果计算的Gi值大于表中的临界值GP(n)，则能判断该测量数据是异常值，可以剔除。但是要提醒，临界值GP(n)与两个参数有关：检出水平α (与置信概率P有关)和测量次数n (与自由度f有关)。
▲定检出水平α：如果要求严格，检出水平α可以定得小一些，例如定α＝0.01，那么置信概率P＝1－α＝0.99；如果要求不严格，α可以定得大一些，例如定α＝0.10，即P＝0.90；通常定α＝0.05，P＝0.95。
▲查格拉布斯表获得临界值：根据选定的P值(此处为0.95)和测量次数n(此处为10)，查格拉布斯表，横竖相交得临界值G95(10)＝2.176。
▲比较计算值Gi和临界值G95(10)：Gi＝2.260，G95(10)＝2.176，Gi＞G95(10)。
▲判断是否为异常值：因为Gi＞G95(10)，可以判断测量值14.0为异常值，将它从10个测量数据中剔除。
▲余下数据考虑：剩余的9个数据再按以上步骤计算，如果计算的Gi＞G95(9)，仍然是异常值，剔除；如果Gi＜G95(9)，不是异常值，则不剔除。本例余下的9个数据中没有异常值。
格拉布斯表——临界值GP(n)
P
n 0.95 0.99 P
n 0.95 0.99
3 1.135 1.155 17 2.475 2.785
4 1.463 1.492 18 2.504 2.821
5 1.672 1.749 19 2.532 2.854
6 1.822 1.944 20 2.557 2.884
7 1.938 2.097 21 2.580 2.912
8 2.032 2.231 22 2.603 2.939
9 2.110 2.323 23 2.624 2.963
10 2.176 2.410 24 2.644 2.987
11 2.234 2.485 25 2.663 3.009
12 2.285 2.550 30 2.745 3.103
13 2.331 2.607 35 2.811 3.178
14 2.371 2.659 40 2.866 3.240
15 2.409 2.705 45 2.914 3.292
16 2.443 2.747 50 2.956 3.336
对异常值及统计检验法的解释
■测量过程是对一个无限大总体的抽样：对固定条件下的一种测量，理论上可以无限次测量下去，可以得到无穷多的测量数据，这些测量数据构成一个容量为无限大的总体；或者换一个角度看，本来就存在一个包含无穷多测量数据的总体。实际的测量只不过是从该无限大总体中随机抽取一个容量为n(例如n＝10)的样本。这种样本也可以有无数个，每个样本相当于总体所含测量数据的不同随机组合。样本中的正常值应当来自该总体。通常的目的是用样本的统计量来估计总体参量。总体一般假设为正态分布。
■异常值区分：样本中的正常值应当属于同一总体；而异常值有两种情况：第一种情况异常值不属于该总体，抽样抽错了，从另外一个总体抽出一个(一些)数据，其值与总体平均值相差较大；第二种情况异常值虽属于该总体，但可能是该总体固有随机变异性的极端表现，比如说超过3σ的数据，出现的概率很小。用统计判断方法就是将异常值找出来，舍去。
■犯错误1：将本来不属于该总体的、第一种情况的异常值判断出来舍去，不会犯错误；将本来属于该总体的、出现的概率小的、第二种情况的异常值判断出来舍去，就会犯错误。
■犯错误2：还有一种情况，不属于该总体但数值又和该总体平均值接近的数据被抽样抽出来，统计检验方法判断不出它是异常值，就会犯另外一种错误。
■异常值检验法：判断异常值的统计检验法有很多种，例如格拉布斯法、狄克逊法、偏度-峰度法、拉依达法、奈尔法等等。每种方法都有其适用范围和优缺点。
■格拉布斯法最佳：每种统计检验法都会犯犯错误1和错误2。但是有人做过统计，在所有方法中，格拉布斯法犯这两种错误的概率最小，所以推荐使用格拉布斯法。
■多种方法结合使用：为了减少犯错误的概率，可以将3种以上统计检验法结合使用，根据多数方法的判断结果，确定可疑值是否为异常值。
■异常值来源：测量仪器不正常，测量环境偏离正常值较大，计算机出错，看错，读错，抄错，算错，转移错误。
代码实现：
#include "math_flash.h"
#define MIN_SMOOTHING_SAMPLE_AMOUNT 5 //最小平滑样本量
typedef enum //数据有效标志定义
{
VALID_FLAG_VALID = 0x00, //数据有效
VALID_FLAG_VALID_LESS = 0x01, //有效数据不足
VALID_FLAG_MORE_ERROR_DATA = 0x02 //粗大误差太多
}ValidFlagDef;
const float f32GrubbsTable[SMOOTHING_SAMPLE_SIZE] = { //格鲁布斯临界值
0.0f, 0.0f, 0.0f, 0.0f, 1.672f,
1.822f, 1.938f, 2.032f, 2.11f, 2.176f,
2.234f, 2.285f, 2.331f, 2.371f, 2.409f,
2.443f, 2.475f, 2.504f, 2.532f, 2.557f
};
/*******************************************************************
功能：格鲁布斯法求平均值和标准差
*******************************************************************/
uint8_t calc_grubbs_value(float32_t * pSrc, uint32_t blockSize, float32_t * mean, float32_t * stdErr)
{
volatile float32_t f32Min = 0.0f, f32Max = 0.0f, f32StdError = 0.0f, f32GMin = 0.0f, f32GMax = 0.0f, f32Mean = 0.0f;
uint32_t u32MinIndex = 0, u32MaxIndex = 0;
uint8_t result = VALID_FLAG_VALID;
uint8_t u8Exit = 0;
//剔除粗大误差
while(!u8Exit)
{
arm_min_f32(pSrc, blockSize, (float32_t*)&f32Min, &u32MinIndex); //最小值
arm_max_f32(pSrc, blockSize, (float32_t*)&f32Max, &u32MaxIndex); //最大值
arm_std_f32(pSrc, blockSize, (float32_t*)&f32StdError); //标准差
arm_mean_f32(pSrc, blockSize, (float32_t*)&f32Mean); //平均值
if(f32StdError <= 0.001f)
{
break;
}
f32GMin = (f32Mean - f32Min) / f32StdError;
f32GMax = (f32Max - f32Mean) / f32StdError;
if(f32GMin <= f32GMax && f32GMax > f32GrubbsTable[blockSize-1])
{
arm_copy_f32(&(pSrc[u32MaxIndex+1]), &(pSrc[u32MaxIndex]), blockSize-u32MaxIndex-1);
}
else if(f32GMin > f32GMax && f32GMin > f32GrubbsTable[blockSize-1])
{
arm_copy_f32(&(pSrc[u32MinIndex+1]), &(pSrc[u32MinIndex]), blockSize-u32MinIndex-1);
}
else
{
u8Exit = 1;
}
blockSize--;
if(blockSize < MIN_SMOOTHING_SAMPLE_AMOUNT)
{
result = VALID_FLAG_MORE_ERROR_DATA;
u8Exit = 1;
}
}
*mean = f32Mean;
*stdErr = f32StdError;
return result;
} |
阿莫论坛20周年了！感谢大家的支持与爱护！！
知道什么是神吗？其实神本来也是人，只不过神做了人做不到的事情 所以才成了神。 （头文字D, 杜汶泽）
|
[复制链接]
发表于 2017-12-20 11:27:51
