时间序列异常值检验替换—

时间序列异常值检验替换——基于Hampel滤波器

news2024/9/23 9:26:10

Hampel滤波器作为一种强大的时间序列异常值处理工具，在数据清洗和分析中发挥着重要作用。本文将深入研究Hampel滤波器的原理和数学推导，并通过实际代码演示其在异常值处理中的应用

1. Hampel滤波器简介

1.1 什么是Hampel滤波器？

Hampel滤波器是一种基于中值和中值绝对偏差（MAD）的滤波器，旨在识别和去除时间序列数据中的异常值。相对于传统均值和标准差方法，Hampel滤波器对异常值更具鲁棒性

1.2 为什么选择Hampel滤波器？

在真实世界的数据中，异常值往往会对统计分析产生较大影响。Hampel滤波器通过使用中值和MAD，适应异常值的存在，提高异常值检测的准确性

2. Hampel滤波器原理

2.1 中值和MAD

Hampel滤波器的核心在于中值的计算和MAD的求解。中值代表数据的中间值，而MAD度量了数据点与中值之间的离散程度

2.2 数学推导

3. 代码实例

3.1 数据生成与可视化

首先，我们生成一个包含异常值的正弦波数据，并通过可视化展示原始数据，感受一下异常值在数据集的呈现

import matplotlib.pyplot as plt
import numpy as np

# 设置中文字体
plt.rcParams['font.sans-serif'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False

# 生成包含异常值的正弦波数据
def generate_data_with_outliers():
    time = np.linspace(0, 10, 100)
    signal = np.sin(time) + np.random.normal(0, 0.1, 100)

    # 添加异常值
    outliers_indices = [20, 40, 60, 80]
    outliers_values = [2.0, -1.9, 2.1, -0.5]
    
    for index, value in zip(outliers_indices, outliers_values):
        signal[index] = value

    return time, signal

# 生成包含异常值的数据
time, data = generate_data_with_outliers()

# 获取异常值的索引和数值
outliers_indices = [20, 40, 60, 80]
outliers_values = [data[i] for i in outliers_indices]

# 绘制原始数据
plt.figure(figsize=(10, 6))
plt.plot(time, data, label='原始数据')

# 仅将异常值标红
plt.scatter(time[outliers_indices], outliers_values, c='red', marker='o', label='异常值', s=50)

plt.title('带有异常值的原始数据')
plt.legend()
plt.show()

3.2 Hampel滤波器定义

接下来，定义Hampel滤波器函数

def hampel(vals_orig, k=7, t0=3):
    """
    使用Hampel滤波器去除时间序列中的异常值。

    参数：
        - vals_orig: numpy数组，原始时间序列数据
        - k: 整数，滤波器窗口的大小（半窗口大小为k/2）
        - t0: 浮点数，用于异常值检测的阈值

    返回：
        - vals_filt: numpy数组，经过滤波的时间序列数据
        - outliers_indices: list，异常值的索引列表
    """
    # 创建输入数据的副本
    vals_filt = np.copy(vals_orig)
    outliers_indices = []

    # 定义Hampel滤波器函数
    n = len(vals_orig)

    for i in range(k, n - k):
        # 提取窗口
        window = vals_orig[i - k:i + k + 1]

        # 计算中值和中值绝对偏差（MAD）
        median = np.median(window)
        mad = np.median(np.abs(window - median))

        # 检查当前值是否为异常值
        if np.abs(vals_orig[i] - median) > t0 * mad:
            # 用中值替换异常值
            vals_filt[i] = median
            # 记录异常值的索引
            outliers_indices.append(i)

    return vals_filt, outliers_indices

3.3 Hampel滤波器实现

接下来，我们使用Hampel滤波器对数据进行处理，去除异常值

# 对数据应用Hampel滤波器
filtered_data, outliers_indices = hampel(data)

3.4 滤波效果展示

最后，我们展示Hampel滤波器处理后的数据，并通过可视化对比原始数据，特别突出异常值的识别和去除效果


# 绘制原始数据和滤波后的数据
plt.figure(figsize=(10, 6))
plt.plot(time, data, label='原始数据')
plt.scatter(time[outliers_indices], data[outliers_indices], c='red', marker='o', label='异常值', s=50)
plt.plot(time, filtered_data, label='滤波后数据', linestyle='--', color='red')
plt.fill_between(time, filtered_data, data, color='red', alpha=0.2, label='异常值区域')
plt.title('使用Hampel滤波器去除异常值')
plt.grid(True)
plt.legend()
plt.show()

# 打印异常值的索引
print("异常值的索引：", outliers_indices)