关于机器学习:用于时间序列异常检测的学生化残差-studentized-residual的理论和代码实现

69次阅读

共计 3763 个字符,预计需要花费 10 分钟才能阅读完成。

异样检测是指数据迷信中可帮忙发现数据集中的异样值有用的技术。异样检测在解决工夫序列数据时特地有用。例如工夫序列数据来自传感器测量后果(例如压力和温度),因为设施故障和瞬态景象等问题蕴含许多异样点,异样检测有助于打消这些点异样值,以优化工夫序列数据中的信号。对于销量预测等需要异样点也能够示意为流动或者营销的记录,能够进行重点剖析。

在这篇文章中,将介绍一个可用于检测异样值的简略但高效的算法,该算法来自论文(https://www.researchgate.net/…)

工夫序列异样检测算法

下图阐明了能够在测量传感器的日常操作中察看到的工夫序列数据的典型示例。橙色线示意根底信号,而蓝色峰示意可能因为测量读数中的尖峰而呈现的异样点。在这种状况下,咱们所需的异样检测工具的目标是通过删除那些异样点来简略地细化信号。

咱们将点异样定义为与其预期值齐全不同的任何点。在这篇文章中展现的算法是通过应用多项式回归和学生化残差(studentized residual 也叫学生化删除的残差)来辨认这些异样。

第一步是定义一条多项式曲线,为数据集的根底信号提供预计。

为了将这条曲线拟合到数据中,必须通过最小化某个损失函数来确定系数(直到 N 级)。通常损失函数能够定义为一般残差的最小化,其计算为理论值与其预测值之间的差别。

然而应用这种形式辨认异样值存在一些局限性。异样的存在可能会导致回归系数呈现偏差,从而无奈标记异样值。这个限度能够通过移除评估为残差的数据点并在数据上从新拟合多项式回归来解决,并且这个操作能够反复屡次。

上述办法对于每个数据点,必须从新拟合回归模型。然而有一个数学技巧能够通过仅在整个数据集上计算一次回归拟合来确定删除的残差并将它们标准化。最终残差被称为学生化删除残差(行将残差除以其标准差),所以能够按如下形式计算:

数学技巧是应用 Hat 矩阵的对角线来调整每个观测值 i 的 SSE(误差平方和)。这个 Hat 矩阵计算为:

而后,学生化删除的残差可用于通过查找异样大的偏差来查找异样点。这些残差遵循具备 n-1-p 自由度的 T 散布,因而能够通过计算定义为的 Bonferroni 临界值来建设适合的阈值:

α 是显著性值(通常设置为 0.05),能够辨认咱们冀望在预期置信区间内的值。

而后能够应用此阈值来辨认和删除数据集中的任何点异样。此外还能够对 BC 值利用一个校对因子以取得更好的后果(在论文中发现 1/6 的值能够提供最佳性能)。

Python 实现

为了生成简略的试验数据集,咱们应用增加了高斯噪声的基线多项式曲线。而后,咱们将 20 个随机点增加到咱们认为是异样的数据中。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
rng = np.random.default_rng(12345)

# Define index
nx = 1000
index = pd.date_range(start="1970", periods=nx, freq="1T")

# Define signal and noise
x = np.linspace(0, 10, nx)
signal =  2*x**2- 10* x + 2
noise = np.random.normal(loc=100, size=nx, scale=2)
y = noise + signal

# Add anomalies
anom_num = rng.integers(low=0, high=200, size=20)
anom_ids = rng.integers(low=0, high=nx, size=20)
y[anom_ids] = anom_num
is_anom = [item in anom_ids for item in range(nx)]

# Pandas DataFrame and plot
raw_data = pd.Series(y, index = index)
clean_data = raw_data[np.invert(is_anom)]
raw_data.plot(figsize=(15,5))
clean_data.plot()
raw_data[anom_ids].plot(style='o')
plt.legend(['Raw Data', 'Clean Data','Anomalous Points'])

在将日期工夫索引转换为整数后,能够应用 numpy 对该数据集执行多项式回归(在这种状况下,它转换为从 1970-01-01 开始的以毫秒为单位的工夫)。

# Transform variables to lists
x = (np.array(raw_data.index, dtype=np.int64) - raw_data.index[0].value) / 1e9
y = raw_data.to_numpy()

# Create a polynomial fit and apply the fit to data
poly_order = 2
coefs = np.polyfit(x, y, poly_order)
y_pred = np.polyval(coefs, x)

Hat 矩阵的计算能够如下进行:

# Calculate hat matrix
X_mat = np.vstack((np.ones_like(x), x)).T
X_hat = X_mat @ np.linalg.inv(X_mat.T @ X_mat) @ X_mat.T
hat_diagonal = X_hat.diagonal()

依据 T 散布计算学生化删除残差及其对应的 p 值能够如下执行:

from scipy.stats import t as student_dist

# Calculate degrees of freedom
n = len(y)
dof = n - 3  # Using p = 2 from paper

# Calculate standardised residuals 
res = y - y_pred
sse = np.sum(res ** 2)
t_res = res * np.sqrt(dof / (sse * (1 - hat_diagonal) - res**2))

最初,用 Bonferroni 临界值过滤掉异样,阈值和绘制后果如下。

# Return filtered dataframe with the anomalies removed using BC value
alpha=0.05
bc_relaxation = 1/6
bc = student_dist.ppf(1 - alpha / (2 * n), df=dof) * bc_relaxation
mask = np.logical_and(t_res < bc, t_res > - bc)

# Plot anomalous and filtered data
ax=raw_data.plot(figsize=(15,5))

raw_data[mask].plot(ax=ax)
raw_data[anom_ids].plot(style='o')

raw_data[np.invert(mask)].plot(style='ok')
plt.legend(['Raw Data','Filtered with Anomaly Detector', 'Missed Anomalies', 'Anomalous Points Found'],bbox_to_anchor=(0.5, 1.05))

后果

将生成数据应用典型的准确率和召回率分类指标来确定咱们的模型的工作状况。在咱们的合成数据集上应用上述办法,失去了 95% 的召回率和 86% 的准确率。这意味着咱们只错过了 20 个异样中的一个,并且只有 3 个误报。

到目前为止所有都很好,然而这个算法将如何解决理论的实在数据呢?为了测试这一点,咱们应用凋谢工业数据,这是一个可供公众应用的公开数据集(https://openindustrialdata.com/)。那数据量比拟大,但让咱们抉择一个对工业过程十分重要的传感器数据。在此示例中,将应用一个压力变送器来测量第一级压缩机的冲击压力(标签的内部 ID 为 pi:160696)并查看过来 50 天的每小时值。

人肉的查看表明异样的确已胜利打消,并且信号已被细化以供进一步剖析。

总结

最初大家可能对这个术语听起来感到十分的奇怪,这里做一个简略的解释

残差(residual)= 察看值 – 预期值

一个好的线性回归残差应该是合乎正态分布的,然而能够通过变换使得残差合乎自由度为 n -1-p 的 t 散布。t 散布因为叫做 student‘s t distribution,所以这个变换后的残差值就是 studentized residual。通过查看这个值能够晓得察看值的散布状况,能够用来寻找 outlier 及确定其 p value。

还记得 t 散布的背景吧:过后的业余统计学研究者 W. Gosset 以笔名 Student 发表对于 T 散布的统计学史地标性文献,所以 T 散布又被称作学生 t - 散布(Student’s t-distribution)

学生化这个词其实就是 studentized 的中文直译,因为约定俗成了所以也没什方法,studentized 就是把其余散布转换成 t 散布,所以其实 studentized residual 翻译为 T 化残差,要比 学生化残差 更天然,也更好了解。

https://www.overfit.cn/post/7dba63d4464c4e8f8881331457541e29

作者:Andris Piebalgs

正文完
 0