查找时间序列数据中异常值的终极指南（第 2 部分）

数云界

于 2024-07-17 17:43:28 发布

阅读量848

点赞数 9

分类专栏：人工智能大语言模型多模大模型文章标签：人工智能大数据

本文链接：https://blog.csdn.net/2401_85233349/article/details/140495250

版权

人工智能同时被 3 个专栏收录

16 篇文章 0 订阅

订阅专栏

大语言模型

12 篇文章 0 订阅

订阅专栏

多模大模型

9 篇文章 0 订阅

订阅专栏

用于时间序列分析中异常值检测的有效机器学习方法和工具

异常值：那些可能颠覆统计模型、误导预测并扰乱决策过程的数据点。

本文是专门讨论时间序列数据中异常值的识别和管理的四部分系列文章的第二部分。

本系列的第一篇文章是关于探索视觉和统计方法，以有效识别时间序列数据中的异常值：

第二篇文章将专门介绍异常值检测的机器学习方法。

鉴于其重要性和复杂性，它们值得专门讨论！

第三篇文章探讨了如何管理这些异常值的各种策略，包括移除、保留和封顶技术，提供了一些处理异常值的实用方法。

在第四篇也是最后一篇文章中，我将继续探讨管理异常值的方法，重点关注归因和转换方法，以及评估异常值处理的影响。

本文的原始数据及代码可问博主索取

了解单变量和多变量数据

在第一篇文章中，我提到了两种不同类型的时间序列数据以及每种数据最合适的异常值检测方法。

由于我们现在正在探索机器学习方法，因此值得回顾一下哪些模型适用于单变量和多变量数据。

记住：

单变量数据：这涉及随时间变化的单个变量或特征。重点是检测单个时间序列中的异常。典型示例包括每日股票价格、每月销售数据或年度天气数据。

适合单变量数据的方法往往更简单、更直接。

多变量数据：这同时涉及多个变量或特征。在检测过程中考虑这些变量之间的关系和相互作用。

例如：多元时间序列可以包括温度、压力、风速的每日测量值，所有这些都同时记录。

对于具有多个变量的数据，隔离森林、LOF 和自动编码器自然适合处理高维数据。我们将在本文中探讨这些方法。

这些机器学习模型（例如自动编码器）是多功能工具，不仅限于分析多个变量，因为它们可以适应单变量数据。怎么做到的？

想象一下教一个模型解释单个数据序列的流动。通过将其压缩为更简单的格式，然后尝试重新创建它，模型可以学习典型的模式。

当呈现新数据时，原始版本和重建版本之间的显著差异（称为重建误差）会起到警报的作用，预示着潜在的异常。

数据说明

本文使用的所有数据集都是我专门为本次分析而生成的合成数据集。没有使用任何外部数据集。

使用机器学习方法检测异常值

用于异常值检测的孤立森林

孤立森林是一种广泛使用、功能强大的无监督机器学习算法，用于大数据集中的异常检测。

它之所以脱颖而出，是因为它采用独特的方法来隔离异常，而不是识别正常数据模式。它特别适用于检测大型数据集中的异常值。

假设：假设异常值比常规数据点出现的频率更低且更为孤立。

孤立森林构建了一组决策树，并通过识别需要较少分割的数据点来隔离树结构中的异常值。异常值在树中的路径预计较短，如下图所示。

该算法适用于具有低维特征空间（与观测值数量相比特征或维度较少）的时间序列数据，以及异常值的分布与常规数据点有显著差异的情况。

隔离森林是数据科学家的重要工具，特别是在异常检测领域。

想知道数据科学家还需要哪些其他技能？找出哪些技能对你的成功至关重要：

让我们逐步了解隔离森林的实施过程。

首先，让我们可视化数据：

plt.figure(figsize=( 10 , 6 )) 
plt.plot(summer_data.index, summer_data[ 'Mn_integ' ], marker= 'o' , linestyle= '-' ) 
plt.title( '夏季的alc_TA_integ' ) 
plt.xlabel( '日期' ) 
plt.ylabel( 'alc_TA_integ' ) 
plt.xticks(rotation= 45 ) 
plt.grid( True ) 
plt.tight_layout() 
plt.show()

该图显示了几年中几个夏季测量的大坝化学参数的数据点。

现在，让我们应用孤立森林：

从sklearn.preprocessing导入StandardScaler
从sklearn.ensemble导入IsolationForest 


scaler = StandardScaler() 
np_scaled = scaler.fit_transform(values.values.reshape(- 1 , 1 )) 
data = pd.DataFrame(np_scaled) 


outliers_fraction = float ( .2 ) 
model = IsolationForest(contamination=outliers_fraction) 
model.fit(data)

关于孤立森林的代码：

Scikit StandardScaler -learn 用于标准化数据集，通过赋予每个特征零平均值和一标准差来确保每个特征的贡献均等。

这种规范化对于时间序列数据尤其重要，因为值的范围可能变化很大。该fit_transform方法应用于数据，根据算法的需要对其进行重塑。

IsolationForest然后创建该模型的一个实例，并将contamination参数设置为 0.2。

该参数是对数据集中异常值比例的估计，指导模型预期 20% 的数据点是异常。

summer_data[ 'anomaly' ] = model.predict( data ) 
fig, ax = plt.subplots(figsize=( 10 , 6 )) 
a = summer_data.loc[summer_data[ 'anomaly' ] == - 1 , [ 'Mn_integ' ]] 
ax.plot(summer_data.index, summer_data[ 'Mn_integ' ], color= 'black' , label = 'Normal' ) 
ax.scatter(a.index,a[ 'Mn_integ' ], color= 'red' , label = 'Anomaly' ) 
plt.legend() 
plt.show()

该模型在识别异常数据点方面做得很好！

调整污染参数

那么参数呢contamination？我们如何调整它？

如果观察到大量的真阴性（即正常数据点被错误地归类为异常），这表明您的污染参数可能设置得太高。

该参数反映了数据中异常值的预期比例，将其设置得过高可能会导致许多数据点被错误地归类为正常。

另一方面，如果您发现在预期的位置缺少红点（表示异常），则可能意味着污染参数设置得太低。

调整此参数对于识别数据集中的真正异常值至关重要！

Prophet 用于异常值检测

您可能听说过这个，因为 Prophet 是一个著名的时间序列预测模型。

它旨在处理时间序列数据并对未来趋势做出预测，但也可以用于检测异常值。

它的工作原理是将时间序列分解为三个主要部分：

趋势：捕捉数据随时间变化的总体方向，解释长期的增加或减少，
季节性：模拟周期性波动，例如在固定时间段内重复的每日、每周或每年模式。
假期/活动：包括可能影响数据的已知事件或假期的影响。

我们首先准备一个关键指标的数据，在本例中，我们在大坝中测量了相同的化学参数，名为“Mn_integ”，捕获夏季数据：

# Prepare the data
df_mn = summer_data[['Date', 'Mn_integ']].rename(columns={'Date': 'ds', 'Mn_integ': 'y'})
# Train the Prophet model for 'Mn_integ'
model_mn = Prophet()
model_mn.fit(df_mn)
# Make predictions for both columns
future_mn = model_mn.make_future_dataframe(periods=0)
forecast_mn = model_mn.predict(future_mn)

在这里，我们分离出我们希望预测的特征和时间列“日期”。日期列重命名为“ds”，特征列重命名为“y”，这是Prophet 模型的要求。

什么是 periods 参数？

该periods参数至关重要，因为它指定了超出模型训练期间提供的历史数据范围的未来数据点的数量。

设置此参数有助于定义应生成预测的时间范围。

示例：设置periods为 0时，表示除了历史数据集中的最后一个日期之外，不会创建其他未来日期。本质上，该模型仅专注于为数据中的现有日期生成预测。

这在目标不是预测未来而是评估模型在已知数据上的性能或检测现有数据集中的异常（就像我们的情况一样！）的情况下特别有用。

plt.figure(figsize=(10, 6))
model_mn.plot(forecast_mn, xlabel='Date', ylabel='Mn_integ', ax=plt.gca())
plt.title('Prophet Forecast for Mn_integ')
plt.show()

上图中，该图展示了预测的“Mn_integ”值。该图将实际数据点显示为黑点，将预测值显示为蓝线，将预测区间显示为阴影蓝色区域，突出显示预测准确性和潜在异常。如您所见，初始结果非常糟糕。但我们的目标是检测异常，而不是预测。让我们看看即使预测不佳，模型也会检测到哪些异常值。

我们应该怎么做？

请注意，forecast_mn 第一个代码片段包含“Mn_integ”特征的未来预测。但不仅仅是这样。它还有更多组件，主要是：

ds：进行预测的日期。
yhat：每个日期的“Mn_integ”预测值。
yhat_lower：预测区间的下限，表示最小预期值。
yhat_upper：预测区间的上限，表示最大预期值。

我们将利用这些列来检测异常值。

# Merging forecasted data with your original data
forecasting_final_mn = pd.merge(forecast_mn[['ds', 'yhat', 'yhat_lower', 
                       'yhat_upper']], df_mn, how='inner', on='ds')

# Calculate the prediction error and uncertainty
forecasting_final_mn['error'] = forecasting_final_mn['y'] - forecasting_final_mn['yhat']
forecasting_final_mn['uncertainty'] = forecasting_final_mn['yhat_upper'] - forecasting_final_mn['yhat_lower']

# Anomaly detection
factor = 1.5
forecasting_final_mn['anomaly'] = forecasting_final_mn.apply(
    lambda x: 'Yes' if (np.abs(x['error']) > factor * x['uncertainty']) else 'No', axis=1
)

为了识别异常，我们可以计算预测误差（y - yhat）和不确定区间（yhat_upper - yhat_lower）。

然后根据定义的阈值检测异常值：如果绝对误差超过不确定区间的 1.5 倍，则该数据点被视为异常。

这个 1.5 的系数可以根据异常检测过程所需的灵敏度进行调整。

import plotly.express as px
# Visualization with Plotly
color_discrete_map = {'Yes': 'rgb(255,12,0)', 'No': 'blue'}
fig = px.scatter(forecasting_final_mn, x='ds', y='y', color='anomaly', title='Anomaly Detection in Mn_integ',
                 color_discrete_map=color_discrete_map)
fig.show()

很好，该方法检测到了一些异常值。请记住，您可以根据特定数据和异常值检测需求，通过调整变化点、假期、季节性、趋势灵活性和不确定性间隔等参数来优化Prophet 模型！

局部异常因子 (LOF)

局部异常因子 (LOF) 因其识别异常的巧妙方法而脱颖而出。

该方法对于簇密度差异很大的数据集尤其有效。

#Apply LOF Algorithm
lof = LocalOutlierFactor(n_neighbors=20, contamination=0.08)
y_pred = lof.fit_predict(y.reshape(-1, 1))
outlier_scores = lof.negative_outlier_factor_
is_outlier = y_pred == -11</span></span></span></span></span>

LOF 的运行方式是首先检查每个数据点的局部邻域，由最近的邻居定义。

这第一步至关重要，因为它为评估一个集群的典型“人群”密度奠定了基础。

n_neighbors指定测量局部密度时要考虑的邻居数。它是 k-最近邻居中的“k”。典型的起点是 10 或 20。

如果你增加这个数字密度估计变得更加稳定，对单个数据点的敏感度降低，这有助于带有噪音的数据集。

如果减少它，它可以使算法对局部数据结构更加敏感，这可能有助于检测密集集群中的异常。

contamination 表示预计为异常值的数据点占总数据点的比例。如果你认为你的数据集包含更多异常值，或者初始结果缺少一些潜在的异常值，则增加此值将使算法在将点标记为异常值时更具包容性。

plt.figure(figsize=(12, 6))
plt.plot(x, y, 'o', label='Data points')
plt.plot(x[is_outlier], y[is_outlier], 'ro', label='Outliers')
plt.title("Time-Series Data with Outliers Detected by LOF")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()
# Output the indices of outliers and their scores for review
outliers_detected = {
    "indices": np.where(is_outlier)[0],
    "scores": outlier_scores[is_outlier]
}
outliers_detected

图 7：使用 LOF（局部异常值因子）检测异常值 | 图片来自作者。

一旦建立了邻域， LOF 就会计算每个点相对于该本地人群的密度。

该方法的本质在于将一个点的密度与其相邻点的密度进行比较。

密度明显低于其邻居的点被标记为异常值，表明它是一个异常点，而不是集群的常规成员。

LOF 算法在具有多样化聚类密度的数据集中表现出色。它旨在检测与其特定聚类相关的异常值，使其成为数据集未标记的无监督学习场景的可靠选择。

然而，使用 LOF 并非没有挑战。它对设置的参数（例如邻居数量）非常敏感；这里的失误可能会导致检测精度不理想。

此外，由于对点之间的多次距离计算需要大量计算，LOF 在较大的数据集中的可扩展性会受到影响。

基于聚类的异常检测

我可以就这些方法写一整篇文章。在这里，我将简要介绍一些常用方法以及如何将它们应用于时间序列异常检测。

每种方法都有其独特的特点，适用于不同类型的数据和分析目标。

分区方法通过优化标准（例如最小化簇内距离）将数据划分为特定数量的簇，K 均值就是一个典型的例子。

对于时间序列数据，这通常涉及将数据转换为合适的格式，例如表示时间序列段的特征向量或使用降维技术，如主成分分析（PCA）。

它最适合于具有球形聚类且聚类数量已知的大型数据集，并且这些方法计算效率高且易于实现。

然而，它们对质心的初始位置很敏感，并且难以处理大小和密度各异的聚类。在时间序列中，定义适当的特征和处理时间依赖性可能具有挑战性。

# Transform time-series data into feature vectors using a sliding window approach
window_size = 10
X = np.array([amplitude[i:i+window_size] for i in range(len(amplitude) - window_size)])

# Standardize the data
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Apply K-means clustering
n_clusters = 3
kmeans = KMeans(n_clusters=n_clusters, random_state=42)
labels = kmeans.fit_predict(X_scaled)

# Prepare data for visualization
# We'll plot the first point in each window
time_plot = time[:len(labels)]
amplitude_plot = amplitude[:len(labels)]

# Visualize the clusters
plt.figure(figsize=(12, 6))
colors = ['blue', 'green', 'red']
for i in range(n_clusters):
    plt.scatter(time_plot[labels == i], amplitude_plot[labels == i], color=colors[i], label=f'Cluster {i}')
plt.scatter(time[outlier_indices], amplitude[outlier_indices], color='gold', edgecolor='black', s=100, label='Outliers')
plt.title("Time-Series Data Clustering with K-means")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()

# Check cluster centers
cluster_centers = kmeans.cluster_centers_
print("Cluster Centers (Scaled):", cluster_centers)

分层方法可以构建聚类树而不需要预先指定聚类数量。

当集群结构未知且可能聚集（建立）或分裂（分解）时，此方法很有用。

这些方法通过树状图提供直观的可视化，但计算量大，这会限制它们在大型数据集中的使用。

示例：凝聚聚类

凝聚聚类是一种层次聚类方法。它是最常用的层次聚类技术之一，在将数据转换为合适的格式后，特别适用于识别时间序列数据中的嵌套聚类。

import numpy as np
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.cluster import AgglomerativeClustering
from sklearn.preprocessing import StandardScaler

# Transform time-series data into feature vectors using a sliding window approach
window_size = 10
X = np.array([amplitude[i:i+window_size] for i in range(len(amplitude) - window_size)])

# Standardize the data
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Step 1: Use Agglomerative Clustering
linked = linkage(X_scaled, method='ward')

# Step 2: Plot Dendrogram to determine the optimal number of clusters
plt.figure(figsize=(12, 6))
dendrogram(linked, orientation='top', distance_sort='descending', show_leaf_counts=True)
plt.title('Dendrogram for Hierarchical Clustering')
plt.xlabel('Time Series Segments')
plt.ylabel('Euclidean Distance')
plt.show()

上图的树状图表示了层次聚类过程。它显示了各个数据点如何根据其欧几里得距离逐渐合并到聚类中。分析树状图以选择截止距离，从而产生合理的聚类数量（在本例中为 3）。

树状图的纵轴表示簇合并的距离。距离越大，簇之间的相似性越低。

# Choosing a cutoff (distance) that results in a sensible number of clusters
#From the dendogram, 3 clusters
agg_clustering = AgglomerativeClustering(n_clusters=3, affinity='euclidean', linkage='ward')
labels_agg = agg_clustering.fit_predict(X_scaled)

# Prepare data for visualization
# We'll plot the first point in each window
time_plot = time[:len(labels_agg)]
amplitude_plot = amplitude[:len(labels_agg)]

# Step 3: Visualization of Clusters
plt.figure(figsize=(12, 6))
colors = ['blue', 'green', 'red']
for i in range(3):
    plt.scatter(time_plot[labels_agg == i], amplitude_plot[labels_agg == i], color=colors[i], label=f'Cluster {i}')
plt.scatter(time[outlier_indices], amplitude[outlier_indices], color='gold', edgecolor='black', s=100, label='Outliers')
plt.title("Time-Series Data Clustering with Agglomerative Hierarchical Clustering")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show(

让我们检查一下代码片段。

n_clusters指定要查找的聚类数。这决定了聚集聚类过程结束时预期的聚类数。

为了优化n_clusters，您可以使用树状图直观地识别自然的聚类划分，应用肘部法在聚类紧凑性和数量之间找到平衡，或者计算轮廓分数以确定将数据最好地分成不同组的聚类配置。

基于密度的方法将聚类识别为比其余数据密度更高的区域，将稀疏区域视为噪声或边界点，其中DBSCAN是一种流行的算法。

这些方法适用于不规则或相互交织的簇，并且对噪声和异常值具有很强的鲁棒性。然而，它们的有效性取决于准确设置密度定义参数，而这些参数会随着数据密度的变化而变化。

DBSCAN 的工作原理是将时间序列数据转换为合适的特征空间，例如使用滑动窗口、特征提取或嵌入技术。它可以识别相似模式的聚类，并将与这些聚类有显著偏差的点检测为异常值。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler


# Create a synthetic dataset with different parameters
data, _ = make_blobs(n_samples=350, centers=4, cluster_std=1.0, random_state=24)

# Introduce some outliers
outliers = np.random.uniform(low=-12, high=12, size=(25, 2))
data = np.vstack([data, outliers])

# Standardize the dataset
scaler = StandardScaler()
data = scaler.fit_transform(data)

# Apply DBSCAN with different parameters
epsilon = 0.35
min_pts = 8
dbscan = DBSCAN(eps=epsilon, min_samples=min_pts)
dbscan.fit(data)

该min_samples参数降低了聚类的密度要求，允许更多的点成为核心点，从而可以产生更多的聚类，包括较小或密度较低的聚类。
epsDBSCAN 中的(epsilon)参数定义了两个点之间的最大距离，使它们被视为同一邻域（聚类）的一部分。要调整它，通常从较小的值开始，然后逐渐增加它，检查它如何影响聚类的数量和质量，以在过多的小聚类和过少的大聚类之间取得平衡。

# Apply DBSCAN
epsilon = 2.8  # Increase epsilon to include more points in each cluster
min_samples = 3  # Decrease min_samples to allow smaller clusters

dbscan = DBSCAN(eps=epsilon, min_samples=min_samples)
labels = dbscan.fit_predict(windows_scaled)


# Identify outlier indices
outlier_indices = np.where(labels == -1)[0]

# Convert outlier indices to the corresponding time points
outlier_time_indices = np.array([i + window_size // 2 for i in outlier_indices])

# Plot the time series data with detected outliers
plt.figure(figsize=(12, 6))
plt.plot(time, data, label='Time Series Data')
plt.scatter(time[outlier_time_indices], data[outlier_time_indices], color='red', label='Detected Outliers', marker='x')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.show()

DBSCAN 成功识别了一些异常值，但从图 11 中我们可以发现，可能存在一些误报。

该模型的有效性在很大程度上取决于其参数（eps和min_samples）的正确设置。请确保根据您的需要调整这些参数！

最后，基于网格的方法可以通过将时间维度离散化为区间来适应时间序列数据，从而有效地沿时间轴创建网格。

这种方法有利于高维时间序列数据，可以实现高效处理和并行化。然而，这些方法在聚类时间序列中的有效性在很大程度上取决于所选时间间隔的粒度，并且可能难以应对不同密度的时间序列模式。

因此，虽然它们提供了高效的处理，但它们的实现可能很复杂，并且可能需要自定义才能在时间序列应用程序中有效处理各种模式。鉴于它们的复杂性，我不会在这种情况下提供示例。

自动编码器

自动编码器是一类用于无监督学习的神经网络，特别用于学习数据的有效表示或编码。

有几种使用自动编码器进行时间序列数据异常检测的方法。在这里，我将举例说明其中一种：

原始自动编码器

原始自动编码器由两个主要组件组成——编码器和解码器。编码器将时间序列数据压缩为潜在空间表示，而解码器则从这种压缩形式重建原始数据。

import numpy as np
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Dense

# Normalize data
amplitude = (amplitude - np.mean(amplitude)) / np.std(amplitude)

# Reshape for input to autoencoder (time steps, features)
amplitude = np.expand_dims(amplitude, axis=-1)

# Define autoencoder model
input_layer = Input(shape=(amplitude.shape[1],))
encoded = Dense(32, activation='relu')(input_layer)
decoded = Dense(amplitude.shape[1], activation='linear')(encoded)

autoencoder = Model(input_layer, decoded)
autoencoder.compile(optimizer='adam', loss='mse')

# Train the autoencoder
autoencoder.fit(amplitude, amplitude,
                epochs=50,
                batch_size=32,
                shuffle=True,
                validation_split=0.2)

# Predict on the entire dataset
reconstructed_data = autoencoder.predict(amplitude)

# Calculate Mean Squared Error (MSE) as reconstruction error
mse = np.mean(np.power(amplitude - reconstructed_data, 2), axis=1)

# Set a threshold for outlier detection
threshold = np.percentile(mse, 95)  # Adjust percentile as needed

# Identify outliers
outliers = np.where(mse > threshold)[0]

此代码首先通过减去平均值并除以标准差来对数据进行归一化，然后对其进行重塑以适应自动编码器所需的输入格式。它使用 TensorFlow/Keras 定义和训练Vanilla 自动编码器，计算原始数据和重建数据之间的均方误差 (MSE) 以识别异常值，最后根据 MSE 值第 95 个百分位数设置的阈值检测异常值。

import matplotlib.pyplot as plt

# Plot original and reconstructed data
plt.figure(figsize=(14, 7))
plt.plot(amplitude, label='Original Data', color='blue')
plt.plot(reconstructed_data, label='Reconstructed Data', color='red', linestyle='--')
plt.title('Original vs Reconstructed Data')
plt.xlabel('Time Steps')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()