异常点检测-CSDN博客

# 图像化正太检验
# QQ-正太检验
# PP图的思想是对比正态分布的累计概率值和实际分布的累计概率值
# QQ图比对正态分布的分位数和实际分布的分位数
# 分析散点是否均匀地分布在倾斜线附近，如果偏离倾斜线则不是正太分布
# ===================================================================================
# ===================================================================================
# ===================================================================================
# 定量正太检验
# Shapiro检验法和K-S检验法均属于非参数的统计方法，
# ***原假设：变量服从正态分布***
# 样本量低于5000时，Shapiro检验法比较合理，否则应使用k-S检验法
# 统计量值，概率值P ***概率值P<0.05 则拒绝原假设，即不满足正太分布

pp_qq_plot=sm.ProbPlot(fhData)
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q')

plt.rcParams['figure.figsize'] = (20.0, 12.0) # 单位是inches

# 测试线性数据
testLineData = np.array([i for i in range(100000)])
pp_qq_plot=sm.ProbPlot(testLineData)
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q')

KstestLineDataTest=stats.kstest(rvs=testLineData, args=(testLineData.mean(),testLineData.std()), cdf='norm')
ShapiroLineDataTest=stats.shapiro(x=testLineData)
print("LineData正太检验量化值(k-s)：", KstestLineDataTest)
print("LineData正太检验量化值(shapiro)：", ShapiroLineDataTest)

# 测试正太分布数据
testNormalData = np.random.normal(0.5, 1, 100000)
pp_qq_plot=sm.ProbPlot(testNormalData)
pp_qq_plot.qqplot(line='q')
plt.title('Q-Q')

KstestNormalDataTest=stats.kstest(rvs=testNormalData, args=(testNormalData.mean(),testNormalData.std()), cdf='norm')
ShapiroNormalDataTest=stats.shapiro(x=testNormalData)
print("NormalData正太检验量化值(k-s)：", KstestNormalDataTest)
print("NormalData正太检验量化值(shapiro)：", ShapiroNormalDataTest)

# 正太分布添加噪声测试
Outliers=np.array([1000, 10000, 100000])
testNormalDataOutliers = np.append(testNormalData, Outliers)
KstestNormalDataTestOutliers=stats.kstest(rvs=testNormalDataOutliers, args=(testNormalDataOutliers.mean(),testNormalDataOutliers.std()), cdf='norm')
ShapiroNormalDataTestOutliers=stats.shapiro(x=testNormalDataOutliers)
print("NormalData+Outliers正太检验量化值(k-s)：", KstestNormalDataTestOutliers)
print("NormalData+Outliers正太检验量化值(shapiro)：", ShapiroNormalDataTestOutliers)
print(threeSigmaCheck(pd.Series(testNormalDataOutliers))) # 3倍标准差检验离群值

NormalData正太检验量化值(k-s)： KstestResult(statistic=0.002135691438279741, pvalue=0.7517562982596927)
NormalData正太检验量化值(shapiro)： (0.999978244304657, 0.8738780617713928)

NormalData+Outliers正太检验量化值(k-s)： KstestResult(statistic=0.4967646694998119, pvalue=0.0) NormalData+Outliers正太检验量化值(shapiro)： (0.00036269426345825195, 0.0)

从上面数据的正太检验看出，十万正太分布的数据在插入3个异常值后就不满足正太分布的假设了。因此，包含异常值的数据最好还是不要使用基于正太分布的模型。

异常值检测算法

Grubbs检验

Grubbs' Test为一种假设检验的方法，常被用来检验服从正太分布的单变量数据集中的单个异常值。若有异常值，则其必为数据集中的最大值或最小值。原假设与备择假设如下：

H0H0: 数据集中没有异常值
H1H1: 数据集中有一个异常值

Grubbs' Test检验假设的所用到的检验统计量（test statistic）为:

$G=\frac{\max \left|Y_{i}-\bar{Y}\right|}{s}$

其中，Y¯为均值，s为标准差。

当检验统计量满足以下条件:

$G>\frac{(N-1)}{\sqrt{N}} \sqrt{\frac{\left(t_{\alpha /(2 N), N-2}\right)^{2}}{N-2+\left(t_{\alpha /(2 N), N-2}\right)^{2}}}$

原假设H0被拒绝，其中，N为数据集的样本数，tα/(2N),N−2为显著度(significance level)等于α/(2N)、自由度（degrees of freedom）等于N−2的t分布临界值。实际上，Grubbs' Test可理解为检验最大值、最小值偏离均值的程度是否为异常。

# pip3 install outlier_utils==0.0.3
from outlier import smirnov_grubbs as grubbs
# 默认双边检测，alpha是显著性，输出的结果会删除不满足要求的数据
data = pd.Series([1, 8, 9, 10, 9])
grubbs.test(data, alpha=0.05)
# 1     8
# 2     9
# 3    10
# 4     9
# dtype: int64

data = np.array([1, 8, 9, 10, 9])
grubbs.test(data, alpha=0.05) # array([8, 9, 10, 9])
# 返回离群指数的单侧（最小）检验
grubbs.min_test_indices([8, 9, 10, 1, 9], alpha=0.05) #[3]
#返回异常值的单侧（max）测试
grubbs.max_test_outliers([8, 9, 10, 1, 9], alpha=0.05) #[]
grubbs.max_test_outliers([8, 9, 10, 50, 9], alpha=0.05) #[50]

ESD检验（Extreme Studentized Deviate test）

在现实数据集中，异常值往往是多个而非单个。为了将Grubbs' Test扩展到k个异常值检测，则需要在数据集中逐步删除与均值偏离最大的值（为最大值或最小值），同步更新对应的t分布临界值，检验原假设是否成立。基于此，Rosner提出了Grubbs' Test的泛化版ESD。算法流程如下：

计算与均值偏离最远的残差，注意计算均值时的数据序列应是删除上一轮最大残差样本数据后；

$R_{j}=\frac{\max _{i}\left|Y_{i}-\overline{Y^{\prime}}\right|}{s}, \quad 1 \leq j \leq k$

计算临界值（critical value）

$\lambda_{j}=\frac{(n-j) * t_{p, n-j-1}}{\sqrt{\left(n-j-1+t_{p, n-j-1}^{2}\right)(n-j+1)}}, \quad 1 \leq j \leq k$

检验原假设，比较检验统计量与临界值；若Ri>λj，则原假设H0不成立，该样本点为异常点；
重复以上步骤k次至算法结束。

import math
import numpy as np
from scipy import stats

def get_r(arr):
    m_arr = np.mean(arr)
    d_arr = abs(arr - m_arr)
    s = np.std(arr, ddof=1)  # 样本标准差，注意分母n-1
    out_ind = np.argmax(d_arr)
    return np.max(d_arr) / s, out_ind

def esd(data, alpha=0.05, max_anoms=0.10):
    n = len(data)
    if isinstance(max_anoms, float):
        r = math.ceil(max_anoms * n)
    else:
        r = max_anoms
    outliers = []
    for i in range(1, r + 1):
        p = 1 - alpha / (n - i + 1) / 2
        t = stats.t.ppf(p, n - i - 1)  # p分位点
        _lambda = (n - i) * t / math.sqrt((n - i - 1 + t ** 2) * (n - i + 1))
        arr = np.delete(data, outliers)
        _r, out_ind = get_r(arr)
        if _r > _lambda:  # 超出临界值，视为异常点
            outliers.append(out_ind)
        else:
            break
    return np.delete(data, outliers), data[outliers]

S-ESD(Seasonal ESD)

时间序列数据具有周期性（seasonal）、趋势性（trend），时间序列数据之间是存在时序关联的，因此时间序列数据异常检测时不能作为孤立的样本点处理。

将时间序列数据分解为趋势分量、周期分量和余项分量（残差分量）：

$Y_{v}=T_{v}+S_{v}+R_{v} \quad v=1, \cdots, N$

余项分量（残差分量）在时间序列数据消除趋势分量、周期分量后的结果，个人理解可以看作是近似孤立的样本点。想当然的，将ESD运用于STL分解后的余项分量中，即可得到时间序列上的异常点。但含异常点的数据直接使用基于正太分布的ESD会造成很多误检，因为异常值会改变真实分布的方差和均值，后面会解释。

时间序列数据分解为趋势分量、周期分量和余项分量（残差分量）可以使用STL(Seasonal-Trend decomposition procedure based on Loess)

R提供STL函数，底层为作者Cleveland的Fortran实现，封装为python的包为rstl。

Python的statsmodels实现了一个简单版的时序分解，通过加权滑动平均提取趋势分量，然后对cycle-subseries每个时间点数据求平均组成周期分量：

S-ESD算法用中位数（median）替换掉趋势分量；余项计算公式如下：

$R_{X}=X-S_{X}-\tilde{X}$

其中，X为原时间序列数据，SX为STL分解后的周期分量，X~为X的中位数。

data.value - decomp.seasonal - data.value.median()

S-H-ESD (Seasonal Hybrid ESD)

在得到余项后，由于个别异常值会极大地拉伸均值和方差，为了解决这个问题，S-H-ESD采用了更具鲁棒性的中位数与绝对中位差（Median Absolute Deviation, MAD）替换ESD中的均值与标准差。MAD的计算公式如下：

$M A D=\operatorname{median}\left(\left|X_{i}-\operatorname{median}(X)\right|\right)$

可以使用from statsmodels.robust.scale import mad ，注意区别pandas的mad函数【平均绝对偏差】

实现：

def detect_anoms(data, max_outliers=10, alpha=0.05, one_tail=False, upper_tail=False):
    """
    H-ESD算法
    max_outliers H-ESD算法检测的异常值最大数量
    alpha 接受或拒绝异常的统计显著性水平
    one_tail = False # 只有最大或最小值检测异常值， False时为最大、最小值都检测异常值
    upper_tail = False # 最大值是异常值，配合one_tail使用
    """
    # Maximum number of outliers that S-H-ESD can detect (e.g. 49% of data)
    n = len(data)
    data = fhDataDiffClear
    R_idx = list(range(max_outliers))
    for i in range(1, max_outliers + 1):
        if one_tail:
            if upper_tail:
                ares = data - data.median()
            else:
                ares = data.median() - data
        else:
            ares = (data - data.median()).abs()

        # protect against constant time series
        data_sigma = mad(data)
        if data_sigma == 0:
            break

        ares = ares / float(data_sigma)
        R = ares.max()
        temp_max_idx = ares[ares == R].index.tolist()[0]
        R_idx[i - 1] = temp_max_idx
        data = data[data.index != R_idx[i - 1]]
        if one_tail:
            p = 1 - alpha / float(n - i + 1)
        else:
            p = 1 - alpha / float(2 * (n - i + 1))

        t = student_t.ppf(p, (n - i - 1))
        lam = t * (n - i) / float(sqrt((n - i - 1 + t**2) * (n - i + 1)))

        if R > lam:
            num_anoms = i

    if num_anoms > 0:
        R_idx = R_idx[:num_anoms]
    else:
        R_idx = None
    return R_idx

# 测试
R_idx = detect_anoms(fhDataDiffClear)
print(R_idx)

普通时序数据异常值检测

普通时序信号有时周期性并不平稳，例如变电站的负荷数据，有时会并入其他线路的负荷，导致数据更为复杂，因此直接使用提取周期和趋势，进而获取余项分量（残差分量）的S-H-ESD算法会导致异常检测性能下降（误检和漏检）。而周期和趋势近似于时序数据的平滑，而平滑并不会因为周期性被破坏导致性能下降，因此算法流程为：

1）数据平滑

2）使用H-SED算法检测异常值

数据平滑算法有：

一维曲线平滑

window = 7
# 简单移动平均
simp_moving_avg = fhData.rolling(window=window, min_periods=1).mean()
# 加权移动平均
weighted_moving_avg = fhData.rolling(window=window, min_periods=1, win_type="cosine").mean()
# 指数加权移动平均
ewma = fhData.ewm(alpha=alpha, min_periods=1).mean()

一阶指数平滑：

def exponential_smoothing(alpha, s):
    '''
    一次指数平滑
    :param alpha:  平滑系数
    :param s:      数据序列， list
    :return:       返回一次指数平滑模型参数， list
    '''
    s_temp = []
    s_temp.append(s[0])
    print(s_temp)
    for i in range(1, len(s), 1):
        s_temp.append(alpha * s[i-1] + (1 - alpha) * s_temp[i-1])
    return s_temp

二阶指数平滑：

def double_exponential_smoothing(series, alpha, beta):
    """
        series - dataset with timeseries
        alpha - float [0.0, 1.0], smoothing parameter for level
        beta - float [0.0, 1.0], smoothing parameter for trend
    """
    # first value is same as series
    result = [series[0]]
    for n in range(1, len(series)+1):
        if n == 1:
            level, trend = series[0], series[1] - series[0]
        if n >= len(series): # forecasting
            value = result[-1]
        else:
            value = series[n]
        last_level, level = level, alpha*value + (1-alpha)*(level+trend)
        trend = beta*(level-last_level) + (1-beta)*trend
        result.append(level+trend)

三阶指数平滑（Holt-Winters）：

def exponential_smoothing_3(alpha, s):
     '''
    三次指数平滑
    :param alpha:  平滑系数
    :param s:      数据序列， list
    :return:       返回三次指数平滑模型参数a, b, c， list
    '''
    s_single = exponential_smoothing_1(alpha, s)
    s_double = exponential_smoothing_1(alpha, s_single)
    s_triple = exponential_smoothing_1(alpha, s_double)
    
    a_triple = [0 for i in range(len(s))]
    b_triple = [0 for i in range(len(s))]
    c_triple = [0 for i in range(len(s))] 
    for i in range(len(s)):
        a_triple[i] = 3 * s_single[i] - 3 * s_double[i] + s_triple[i]
        b_triple[i] = (alpha / (2 * ((1 - alpha) ** 2))) * ((6 - 5 * alpha) * s_single[i] - 2 * ((5 - 4 * alpha) * s_double[i]) + (4 - 3 * alpha) * s_triple[i])
        c_triple[i] = ((alpha ** 2) / (2 * ((1 - alpha) ** 2))) * (s_single[i] - 2 * s_double[i] + s_triple[i])
    return a_triple, b_triple, c_triple
  

def predict_value_with_exp_smoothing_3(alpha,s):
      a,b,c=exponential_smoothing_3(alpha,s)
      s_temp=[]
      s_temp.append(a[0])
      for i in range(len(a)):
         s_temp.append(a[i]+b[i]+c[i])
     return s_temp
dataset=predict_value_with_exp_smoothing_3(alpha=0.2,s=dataset)

开尔曼滤波

def Kalman1D(x, damping=1):
    """
    卡尔曼滤波，缺点：耗时较长
    :param x 时间序列
    :damping 协方差，控制参数
    :return x_hat 平滑后的时间序列
    """
    # To return the smoothed time series data
    observation_covariance = damping
    initial_value_guess = x[0]
    transition_matrix = 1
    transition_covariance = 0.1
    kf = KalmanFilter(
        initial_state_mean=initial_value_guess,
        initial_state_covariance=observation_covariance,
        observation_covariance=observation_covariance,
        transition_covariance=transition_covariance,
        transition_matrices=transition_matrix
    )
    x_hat, state_cov = kf.smooth(x)
    return x_hat