数据科学中级机器学习项目Python时间序列

使用ARIMA模型进行单变量时间序列异常检测

Srivignesh_R ,2021年8月9日

文章视频书

这篇文章是作为数据科学博客马拉松的一部分发表的

异常检测

异常检测是在数据中寻找异常的过程。异常点是指那些明显偏离数据一般行为的数据点。在训练之前，检测数据中的异常情况是非常有用的。异常检测的一些例子是欺诈检测、垃圾邮件过滤、CPU使用异常检测、检测服务器使用的异常情况等等。

时间序列数据

在时间序列数据中，数据中的观察值与时间戳有关。数据是按顺序放置的，所以我们不能在训练中洗数据。数据中的观测值是自相关的，这意味着这些观测值与它们之前的观测值高度相关。由于这些限制，时间序列数据必须以一种特殊的方式来处理。

时间序列异常检测

为了检测时间序列数据中的异常情况，我们不能使用传统的异常检测算法，如IQR、Isolation Forest、COPOD等。我们需要以单独的方式来处理时间序列异常检测的任务。

在时间序列数据中寻找异常情况需要遵循的步骤。

1. 检查数据是否是静止的。如果数据不是静止的，将数据转换为静止的。
2. 将时间序列模型与预处理过的数据相匹配。
3. 为数据中的每一个观测值找到平方误差。
4. 找到数据中误差的阈值
5. 如果误差超过该阈值，我们就可以将该观察值标记为异常值。

为什么我们要使用这种方法？

如前所述，时间序列数据是严格按顺序排列的，极易产生自相关。时间序列模型将使用数据进行训练，找到数据的一般行为，并试图对数据进行预测。如果一个观察是正常的，那么预测将尽可能地接近实际值，如果一个观察是异常的，那么预测将尽可能地接近实际值，所以如果我们检查预测的误差，我们可以找到数据中的异常。

图片来源

单变量与多变量的时间序列数据

• 单变量时间序列数据将只包含一个特征（或列）和一个与之相关的时间戳列。
• 多变量时间序列数据将包含一个以上的特征和一个与之相关的时间戳列。

在这篇文章中，我们将看到关于单变量时间序列异常检测。

单变量时间序列异常检测

我们将使用Kaggle的航空乘客数据。你可以在这里找到这些数据。该数据包含每月登上飞机的乘客数量。该数据包含两列，月份和乘客人数。

def test_stationarity(ts_data, column='', signif=0.05, series=False):
    if series:
        adf_test = adfuller(ts_data, autolag='AIC')
    else:
        adf_test = adfuller(ts_data

"""Univariate Time Series"""
data = pd.read_csv('airline-passengers.csv')
data = data.set_index('Month')

图片来源。作者的Jupyter笔记本

首先，我们要用Augmented Dickey-Fuller（ADF）测试来检查数据是否是静止的。该测试将告诉我们数据是否是静止的。

ADF检验会返回统计结果，其中包含p值。p值用于确定数据是否是静止的。如果p值小于0.05，那么数据是静止的；如果数据大于0.05，那么数据是非静止的。0.05被称为显著性水平，对应于95%。你也可以尝试不同的显著性水平来测试静止性。

test_stationarity(data, 'Passengers') # Output: Non-Stationary

检验结果告诉我们，数据是非平稳的，这意味着数据没有恒定的平均数、方差和自相关。为了将数据转换为静止状态，我们可以对数据进行差分。

def differencing(data, column, order):
    differenced_data = data

上面的代码将对数据进行差分，从而使数据成为静止的。现在让我们用ADF检验来检查被差分的数据是否是静止的。

test_stationarity(preprocessed_data, series=True) # Output: Stationary

经过差分后，数据已经成为静止的。现在我们可以使用这些数据来训练时间序列模型。

在这里，我们将使用ARMA（自回归移动平均）模型来预测数据。ARMA模型有两个参数，即p和q。p是指AR（自动回归），q是指MA（移动平均）。自动回归使用以前的滞后期来建立数据模型，而移动平均则使用以前的预测误差来建立数据模型。

自动回归将数据的滞后期作为特征，将提供的数据作为目标，并使用最小二乘法（或线性回归）对数据之间的关系进行建模。如前所述，我们用它来为时间序列数据建模，这些数据是严格意义上的顺序，并且具有自相关。由于残差的自相关性，我们在这里不能直接使用最小二乘法或线性回归。残差的自相关会导致虚假的预测。这就是为什么我们使用一个专门为时间序列数据构建的特殊模型。

为了找到模型的最佳阶数（p, q ），我们必须选择能降低模型AIC的阶数（p, q）。我们可以用下面的代码找到。

max_p, max_q = 5, 5 
def get_order_sets(n, n_per_set) -> list:
    n_sets = [i for i in range(n)]
    order_sets = [
        n_sets[i:i + n_per_set]
        for i in range(0, n, n_per_set)
    ]
    return order_sets
def find_aic_for_model(data, p, q, model, model_name):
    try:
        msg = f"Fitting {model_name} with order p, q = {p, q}n"
        print(msg)
        if p == 0 and q == 1:
            # since p=0 and q=1 is already
            # calculated
            return None, (p, q)
        ts_results = model(data, order=(p, q)).fit(disp=False)
        curr_aic = ts_results.aic
        return curr_aic, (p, q)
    except Exception as e:
        f"""Exception occurred continuing {e}"""
        return None, (p, q)
def find_best_order_for_model(data, model, model_name):
    p_ar, q_ma = max_p, max_q
    final_results = []
    ts_results = model(data, order=(0, 1)).fit(disp=False)
    min_aic = ts_results.aic
    final_results.append((min_aic, (0, 1)))
    # example if q_ma is 6
    # order_sets would be [[0, 1, 2, 3, 4], [5]]
    order_sets = get_order_sets(q_ma, 5)
    for p in range(0, p_ar):
        for order_set in order_sets:
            # fit the model and find the aic
            results = Parallel(n_jobs=len(order_set), prefer='threads')(
                delayed(find_aic_for_model)(data, p, q, model, model_name)
                for q in order_set
            )
            final_results.extend(results)
    results_df = pd.DataFrame(
        final_results,
        columns=['aic', 'order']
    )
    min_df = results_df[
        results_df['aic'] == results_df['aic'].min()
    ]
    min_aic = min_df['aic'].iloc[0]
    min_order = min_df['order'].iloc[0]
    return min_aic, min_order, results_df
min_aic, min_order, results_df = find_best_order_for_model(
    preprocessed_data, ARMA, "ARMA"
)
print(min_aic, min_order)# Output: 1341.1035677173795, (4, 4)

上述代码是如何工作的？

• 我们应该提供p和q的最大值来探索。
• 上面的代码将拟合不同p和q值的模型，并记录每个模型的AIC。
  • 在这里，我使用了多线程来快速执行代码。
• 完成后，使用模型的最小AIC值找到最佳顺序。

现在我们已经发现（4，4）是最佳顺序，我们可以用它来拟合ARMA模型。

def find_anomalies(squared_errors):
    threshold = np.mean(squared_errors) + np.std(squared_errors)
    predictions = (squared_errors >= threshold).astype(int)
    return predictions, threshold

arma = ARMA(preprocessed_data, order=min_order)
arma_fit = arma.fit()
squared_errors = arma_fit.resid ** 2
predictions, threshold = find_anomalies(squared_errors) # threshold = 1428.0094261298002

我们正在用（4，4）的顺序来拟合模型。拟合结果的剩余属性将有每个观察和预测的残差。

我们用下面的公式来寻找误差平方的阈值，超过这个阈值的数据就被认为是异常的。

threshold = mean(squared_errors) + (z * standard_deviation(squared_errors))

z是一个整数。这里我们使用了z=1。你可以根据你的需要调整z的值以获得不同的阈值。还有其他方法可以找到阈值，比如找到误差的第90个百分点。这种方法也可以用来寻找阈值。

预测结果将包含整数0和1，分别表示正常和异常的观察。

data['predictions'] = predictions
print(data)

图片来源。作者的Jupyter笔记本

摘要

• 找出数据是否是静止的。
• 如果数据不是静止的，用差分法将其转换为静止的。
• 找到ARMA模型的最佳顺序
• 拟合ARMA模型并找出误差平方
• 利用误差平方找到阈值
• 超过阈值的误差被认为是异常数据。

本文所展示的媒体不属于Analytics Vidhya所有，由作者自行决定使用。

你也可以在我们的手机APP上阅读这篇文章

标签 :异常检测, ARIMA, 博客马拉松, 时间序列

[

上一篇文章

探索Python中的神奇方法

](www.analyticsvidhya.com/blog/2021/0…)

Srivignesh_R