格鲁布斯算法在stm32的代码实现

只看该作者 · 2019-1-2 09:45

格拉布斯法—异常值判断
▲概述：一组测量数据中，如果个别数据偏离平均值很远，那么这个(这些)数据称作“可疑值”。如果用统计方法—例如格拉布斯(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个数据中没有异常值。

只看该作者 · 2019-1-2 09:46

只看该作者 · 2019-1-2 09:53

对异常值及统计检验法的解释
■测量过程是对一个无限大总体的抽样：对固定条件下的一种测量，理论上可以无限次测量下去，可以得到无穷多的测量数据，这些测量数据构成一个容量为无限大的总体；或者换一个角度看，本来就存在一个包含无穷多测量数据的总体。实际的测量只不过是从该无限大总体中随机抽取一个容量为n(例如n＝10)的样本。这种样本也可以有无数个，每个样本相当于总体所含测量数据的不同随机组合。样本中的正常值应当来自该总体。通常的目的是用样本的统计量来估计总体参量。总体一般假设为正态分布。
■异常值区分：样本中的正常值应当属于同一总体；而异常值有两种情况：第一种情况异常值不属于该总体，抽样抽错了，从另外一个总体抽出一个(一些)数据，其值与总体平均值相差较大；第二种情况异常值虽属于该总体，但可能是该总体固有随机变异性的极端表现，比如说超过3σ的数据，出现的概率很小。用统计判断方法就是将异常值找出来，舍去。
■犯错误1：将本来不属于该总体的、第一种情况的异常值判断出来舍去，不会犯错误；将本来属于该总体的、出现的概率小的、第二种情况的异常值判断出来舍去，就会犯错误。
■犯错误2：还有一种情况，不属于该总体但数值又和该总体平均值接近的数据被抽样抽出来，统计检验方法判断不出它是异常值，就会犯另外一种错误。
■异常值检验法：判断异常值的统计检验法有很多种，例如格拉布斯法、狄克逊法、偏度-峰度法、拉依达法、奈尔法等等。每种方法都有其适用范围和优缺点。
■格拉布斯法最佳：每种统计检验法都会犯犯错误1和错误2。但是有人做过统计，在所有方法中，格拉布斯法犯这两种错误的概率最小，所以推荐使用格拉布斯法。
■多种方法结合使用：为了减少犯错误的概率，可以将3种以上统计检验法结合使用，根据多数方法的判断结果，确定可疑值是否为异常值。
■异常值来源：测量仪器不正常，测量环境偏离正常值较大，计算机出错，看错，读错，抄错，算错，转移错误。

只看该作者 · 2019-1-2 09:54

代码实现：
#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;
}