GESD 离群值检验——理论与 Python 实现

Posted 2023-03-09 zhuo木鸟

文章目录

• GESD 离群值检验
• • 步骤
  • Python 代码

本文主要根据标准 ISO 16269-4: 2010 的 4.3.2 条款

设样本为 $x_1,x_2,\cdots ,x_n$

GESD 离群值检验

GESD 的适用条件是：当样本来源于正态分布时，方可使用。

步骤

1. 首先，根据正态分布和样本的 Q-Q 图，找出 $x_1,x_2,\cdots ,x_n$ 的所有可能的离群值的个数。
2. 定义一个最大离群值数 $m$，$m$ 可以大于等于步骤 1 中得到的可能离群值的个数。
3. 定义 $l=0$，开始迭代过程，记 $X=x_1,x_2,\cdots ,x_n$。
4. 计算检验统计量：
   $$R_l=\frac{{\max }_{x\in X}|x-\bar{X}|}{s(X)}$$
   其中，s 为样本 $X$ 的标准差。
5. 找到 $R_l$ 对应的样本 $x$，将 $x$ 记为 $x^{(l)}$ 并从 $X$ 中删除。
6. 计算临界值：
   $${\lambda }_l=\frac{(n-l-1)\cdot t(p,n-l-2)}{\sqrt{(n-l-2+t^2(p,n-l-2)\cdot (n-l)}}$$
   其中 $p$ 为 $p=(1-\alpha /2{)}^{1/(n-l)}$。$t(p,v)$ 为学生分布 t 的 probability percentage function（百占比点），$v$ 为自由度。也可以说是学生分布 t 的分布函数的反函数吧。
7. 另 $l=l+1$，重复 4-6 直到 $l=m$ 为止。
8. 于是，对于 $x^{(1)},x^{(2)},\cdots ,x^{(m)}$，以及对应的 $R_l$，找到最大的 $l={\max }_{1\le l\le m}l:R_l>{\lambda }_l$，则 $x^{(1)},x^{(2)},\cdots ,x^{(l)}$ 都是离群值

Python 代码

import numpy as np
import pandas as pd
from scipy.stats import t

def perform_GESD(x_original, m=1, alpha=0.05):
    x = x_original.copy()
    l = 1
    # 样本个数
    n = len(x)
    # 定义临时变量
    x_l_list = []
    # flags = True 意味着 Rl > lambda_l
    flags_list = []
    # 迭代直到 l == m
    while l <= m:
        diff = np.abs(x-np.mean(x))/np.std(x, ddof=1)
        # 求检验统计量
        R_idx = np.argmax(diff)
        R_l = diff[R_idx]
        print(R_l)
        # 计算 p
        p = (1-alpha/2)**(1/(n-l))
        # 根据公式 3 计算临界值
        t_percentage = t.ppf(q=p, df=(n-l-2))
        numerator = (n-l-1)*t.ppf(q=p, df=(n-l-2))
        denominator = np.sqrt((n-l-2+t_percentage**2)*(n-l))
        lambda_l = numerator/denominator
        if R_l > lambda_l:
            flags_list.append(True)
        else:
            flags_list.append(False)

        # 保存那些极大值对应的 x
        x_l = x[R_idx]
        x_l_list.append(x_l)
        # 更新迭代 l
        l = l + 1
        x = np.delete(x, R_idx)
    
    # 根据标准定义找出离群值
    idx_flags = list(enumerate(flags_list))
    print(idx_flags)
    True_list = [item[0] for item in idx_flags if item[1]]
    extreme_idx = max(True_list)
    outliers = x_l_list[:(extreme_idx+1)]
    
    outlier_indexs = []
    for outlier in outliers:
        outlier_index = np.where(x_original == outlier)
        outlier_indexs.append(outlier_index)
    # 删除离群值
    x_without_outliers = np.delete(x_original, outlier_indexs)
    
    return outliers, x_without_outliers