用于时间序列异常检测的学生化残差的理论和代码实现

异常检测是指数据科学中的技术，可以帮助找到数据集中的异常值。异常检测在处理时间序列数据时特别有用。例如，时间序列数据来自压力和温度等传感器测量值。由于设备故障和瞬态现象等问题，它包含许多异常值。异常检测有助于消除这些点异常值以优化时间序列数据中的信号。销售预测等需求异常点也可以表现为活动或营销记录，可以在关键点进行分析。在本文中，将介绍一种简单但高效的可用于检测异常值的算法，该算法来自论文（https://www.researchgate.net/...）时间序列异常检测算法下图说明了可以在测量传感器的日常操作中观察到的时间序列数据的典型示例。橙色线代表基础信号，而蓝色峰代表可能由于测量读数中的尖峰而出现的异常值。在这种情况下，我们所需的异常检测工具的目的是通过删除那些异常值来简单地优化信号。我们将点异常定义为与预期值完全不同的任何点。本文中介绍的算法通过使用多项式回归和学生化残差（也称为学生化残差）来识别这些异常。第一步是定义多项式曲线，为数据集的基础信号提供估计。为了使这条曲线适合数据，必须通过最小化一些损失函数来确定系数（直到N级）。通常损失函数可以定义为普通残差的最小化，计算为实际值与其预测值之间的差值。但是以这种方式识别异常值存在一些局限性。异常的存在会使回归系数产生偏差，从而无法标记异常值。这个限制可以通过移除评估为残差的数据点并重新拟合数据的多项式回归来解决，并且这个操作可以重复多次。上述方法必须为每个数据点重新拟合回归模型。但是有一个数学技巧可以确定移除的残差，并通过在整个数据集上仅计算一次回归拟合来对它们进行归一化。最后的残差称为学生化删除残差（即残差除以其标准差），因此可以计算如下：数学技巧是使用Hat矩阵的对角线来调整SSE（误差平方和）.这个Hat矩阵被计算为：学生化的移除残差然后可以用于通过发现异常大的偏差来找到异常值。这些残差服从具有n-1-p自由度的T分布，因此可以通过计算Bonferroni临界值来建立合适的阈值，定义为：α是识别我们期望的显着性值（通常设置为0.05）预期置信区间内的值。然后可以使用此阈值来识别和删除数据集中的任何点异常。也可以对BC值应用校正因子以获得更好的结果（在论文中发现1/6的值可提供最佳性能）。Python实现为了生成一个简单的实验数据集，我们使用添加了高斯噪声的基线多项式曲线。然后我们将20个随机点添加到我们认为异常的数据中。importnumpyasnpimportpandasaspdimportmatplotlib.pyplotaspltrng=np.random.default_rng(12345)#定义indexnx=1000index=pd.date_range(start="1970",periods=nx,freq="1T")#定义信号和noisex=np.linspace(0,10,nx)signal=2*x**2-10*x+2noise=np.random.normal(loc=100,size=nx,scale=2)y=noise+signal#添加anomaliesanom_num=rng.integers(low=0,high=200,size=20)anom_ids=rng.integers(low=0,high=nx,size=20)y[anom_ids]=anom_numis_anom=[iteminanom_idsforiteminrange(nx)]#PandasDataFrameandplotraw_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(['RawData','CleanData','AnomalousPoints'])将日期时间索引转换为整数后，您可以使用numpy进行此操作数据集执行多项式回归（在本例中，它转换为自1970-01-01以来的毫秒数）。#将变量转换为列表x=(np.array(raw_data.index,dtype=np.int64)-raw_data.index[0].value)/1e9y=raw_data.to_numpy()#创建多项式拟合并将拟合应用于datapoly_order=2coefs=np.polyfit(x,y,poly_order)y_pred=np.polyval(coefs,x)帽子矩阵可以这样计算：#计算帽子矩阵X_mat=np.vstack((np.ones_like(x),x)).TX_hat=X_mat@np.linalg.inv(X_mat.T@X_mat)@X_mat.That_diagonal=X_hat.diagonal()从T分布计算学生化删除残差及其相应的p值可以按如下方式执行：来自scipy。statsimporttasstudent_dist#Calculatedegreesoffreedom=len(y)dof=n-3#Usingp=2frompaper#计算标准化残差res=y-y_predsse=np.sum(res**2)t_res=res*np.sqrt(dof/(sse*(1-hat_diagonal)-res**2))最后用Bonferroni临界值过滤异常。阈值和绘图结果如下。#使用BCvaluealpha=0.05bc_relaxation=1/6bc=student_dist.ppf(1-alpha/(2*n),df=dof)*bc_relaxationmask=np.logical_and(t_res-bc)#绘制异常和过滤数据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(['RawData','FilteredwithAnomalyDetector','MissedAnomalies','AnomalousPointsFound'],bbox_to_anchor=(0.5,1.05))结果将使用典型的精度和召回分类指标生成数据，以确定我们的模型的工作情况。在我们的合成数据集上使用上述方法，我们得到了95%的召回率和86%的准确率。这意味着我们只遗漏了20个异常中的1个，并且只有3个误报。到目前为止一切顺利，但该算法将如何处理实际的真实数据？为了对此进行测试，我们使用开放工业数据，这是一个向公众开放的开放数据集(https://openindustrialdata.com/)。那个数据量比较大，但是让我们选择一个对工业过程非常重要的传感器数据。在这个例子中，压力变送器将用于测量第一级压缩机的冲击压力（外部IDpi:160696的标签）并检查最近50天的每小时值。对人肉的检查显示，异常确实已经成功消除，信号也经过提炼以供进一步分析。总结说到最后，大家可能会对这个名词感到很陌生。这是一个简单的解释。残差=观测值-期望值。一个好的线性回归残差应该符合正态分布，但是可以通过变换使残差拟合一个自由度为n-1-p的t分布。t分布称为学生t分布，因此变换后的残差值是学生化残差。通过检查此值，您可以了解观察值的分布，这可用于查找异常值并确定其p值。还记得t分布的背景吗：当时的业余统计研究员W.Gosset以笔名Student发表了一篇关于T分布的统计史上里程碑式的文献，所以T分布也被称为Student'st-distribution(Student'st-distribution)studentized这个词其实就是studentized的中文直译。因为约定俗成了，就没有办法了。studentized就是将其他分布转化为t分布，所以实际上就是将studentized残差转化为T残差，比studentized残差更自然。也比较好理解。https://www.overfit.cn/post/7dba63d4464c4e8f8881331457541e29作者：AndrisPiebalgs