时间序列异常检测算法总结

https://www.biaodianfu.com/timeseries-anomaly-detection.html

时间序列异常检测算法总结

异常的分类

时间序列的异常检测问题通常表示为相对于某些标准信号或常见信号的离群点。虽然有很多的异常类型，但是我们只关注业务角度中最重要的类型，比如意外的峰值、下降、趋势变化以及等级转换（level shifts）。

常见的异常有如下几种：

革新性异常：innovational outlier (IO)，造成离群点的干扰不仅作用于$X_T$，而且影响T时刻以后序列的所有观察值。

附加性异常：additive outlier (AO)，造成这种离群点的干扰，只影响该干扰发生的那一个时刻T上的序列值，而不影响该时刻以后的序列值。

水平移位异常：level shift (LS)，造成这种离群点的干扰是在某一时刻T，系统的结构发生了变化，并持续影响T时刻以后的所有行为，在数列上往往表现出T时刻前后的序列均值发生水平位移。

暂时变更异常temporary change (TC)：造成这种离群点的干扰是在T时刻干扰发生时具有一定初始效应，以后随时间根据衰减因子的大小呈指数衰减。

上面的解释可能不太容易理解，我们结合图片来看一下：

通常，异常检测算法应该将每个时间点标记为异常/非异常，或者预测某个点的信号，并衡量这个点的真实值与预测值的差值是否足够大，从而将其视为异常。使用后面的方法，你将能够得到一个可视化的置信区间，这有助于理解为什么会出现异常并进行验证。

常见异常检测方法

从分类看，当前发展阶段的时序异常检测算法和模型可以分为一下几类：

统计模型：优点是复杂度低，计算速度快，泛化能力强悍。因为没有训练过程，即使没有前期的数据积累，也可以快速的投入生产使用。缺点是准确率一般。但是这个其实是看场景的，并且也有简单的方法来提高业务层面的准确率。这个后面会提到。

机器学习模型：鲁棒性较好，准确率较高。需要训练模型，泛化能力一般。

深度学习模型：普遍需要喂大量的数据，计算复杂度高。整体看，准确性高，尤其是近段时间，强化学习的引入，进一步巩固其准确性方面的领先优势。

3-Sigma

3-Sigma原则又称为拉依达准则，该准则定义如下：假设一组检测数据只含有随机误差，对原始数据进行计算处理得到标准差，然后按一定的概率确定一个区间，认为误差超过这个区间的就属于异常值。

使用3-Sigma的前提是数据服从正态分布，满足这个条件之后，在3-Sigma范围$(\mu –3\sigma ,\mu +3\sigma )$内99.73%的为正常数据，其中$\sigma$代表标准差,$\mu$代表均值，$x=\mu$为图形的对称轴。下面是3-Sigma的Python实现：

import numpy as np


def three_sigma(df_col):
    """
    df_col：DataFrame数据的某一列
    """
    rule = (df_col.mean() - 3 * df_col.std() > df_col) | (
        df_col.mean() + 3 * df_col.std() < df_col
    )
    index = np.arange(df_col.shape[0])[rule]
    out_range = df_col.iloc[index]
    return out_range

对于异常值检测出来的结果，有多种处理方式，如果是时间序列中的值，那么我们可以认为这个时刻的操作属于异常的；如果是将异常值检测用于数据预处理阶段，处理方法有以下四种：

删除带有异常值的数据；
将异常值视为缺失值，交给缺失值处理方法来处理；
用平均值进行修正；
当然我们也可以选择不处理。

Grubbs测试

Grubbs Test为一种假设检验的方法，常被用来检验服从正态分布的单变量数据集（univariate data set）$Y$ 中的单个异常值。若有异常值，则其必为数据集中的最大值或最小值。原假设与备择假设如下：

H0: 数据集中没有异常值
H1: 数据集中有一个异常值

Grubbs Test检验假设的所用到的检验统计量（test statistic）为：

$$G=\frac{max|Y_i-\overline{Y}|}{s}$$

其中，$\overline{Y}$为均值，$s$为标准差。原假设H0被拒绝，当检验统计量满足以下条件：

$$G>\frac{(N-1)}{\sqrt{N}}\sqrt{\frac{(t_{\alpha /(2N),N-2}{)}^2}{N-2+(t_{\alpha /(2N),N-2}{)}^2}}$$

其中，$N$为数据集的样本数，$t_{\alpha /(2N),N-2}$为显著度(significance level)等于$\alpha /(2N)$、自由度（degrees of freedom）等于$N-2$的$t$分布临界值。实际上，Grubbs Test可理解为检验最大值、最小值偏离均值的程度是否为异常。

使用Grubbs测试需要总体是正态分布的。算法流程：

样本从小到大排序
求样本的mean和dev
计算min/max与mean的差距，更大的那个为可疑值
求可疑值的z-score (standard score)，如果大于Grubbs临界值，那么就是outlier

Grubbs临界值可以查表得到，它由两个值决定：检出水平$\alpha$（越严格越小），样本数量$n$，排除outlier，对剩余序列循环做1-4步骤。由于这里需要的是异常判定，只需要判断tail_avg是否outlier即可。

python包 outlier_utils 0.0.3

from outliers import smirnov_grubbs as grubbs

print(grubbs.test([8, 9, 10, 1, 9], alpha=0.05))
print(grubbs.min_test_outliers([8, 9, 10, 1, 9], alpha=0.05))
print(grubbs.max_test_outliers([8, 9, 10, 1, 9], alpha=0.05))
print(grubbs.max_test_indices([8, 9, 10, 50, 9], alpha=0.05))

[ 8  9 10  9]
[1]
[]
[3]

ESD（generalized ESD test）

在现实数据集中，异常值往往是多个而非单个。为了将Grubbs Test扩展到k个异常值检测，则需要在数据集中逐步删除与均值偏离最大的值（为最大值或最小值），同步更新对应的t分布临界值，检验原假设是否成立。基于此，Rosner提出了Grubbs’Test的泛化版ESD（Extreme Studentized Deviate test）。算法流程如下：

1、计算与均值偏离最远的残差，注意计算均值时的数据序列应是删除上一轮最大残差样本数据后。
$$R_j=\frac{ma{x}_i|Y_i-\overline{Y^{\prime }}|}{s},1\le j\le k$$

2、计算临界值（critical value）
$${\lambda }_j=\frac{(n-j)t_{p,n-j-1}}{\sqrt{(n-j-1+t_{p,n-j-1}^2)(n-j+1)}},1\le j\le k$$

3、检验原假设，比较检验统计量与临界值；若$R_i>{\lambda }_j$，则原假设H0不成立，该样本点为异常点

4、重复以上步骤k次至算法结束。

python包 pip install PyAstronomy

import matplotlib.pylab as plt
import numpy as np
from PyAstronomy import pyasl

# Convert data given at:
# http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm
# to array.
x = np.array(
    [
        float(x)
        for x in "-0.25 0.68 0.94 1.15 1.20 1.26 1.26 1.34 1.38 1.43 1.49 1.49 \
          1.55 1.56 1.58 1.65 1.69 1.70 1.76 1.77 1.81 1.91 1.94 1.96 \
          1.99 2.06 2.09 2.10 2.14 2.15 2.23 2.24 2.26 2.35 2.37 2.40 \
          2.47 2.54 2.62 2.64 2.90 2.92 2.92 2.93 3.21 3.26 3.30 3.59 \
          3.68 4.30 4.64 5.34 5.42 6.01".split()
    ]
)

# Apply the generalized ESD
r = pyasl.generalizedESD(x, 10, 0.05, fullOutput=True)

print("Number of outliers: ", r[0])
print("Indices of outliers: ", r[1])
print("        R      Lambda")
for i in range(len(r[2])):
    print("%2d  %8.5f  %8.5f" % ((i + 1), r[2][i], r[3][i]))

# Plot the "data"
plt.plot(x, "b.")
# and mark the outliers.
for i in range(r[0]):
    plt.plot(r[1][i], x[r[1][i]], "rp")
plt.show()

Number of outliers:  3
Indices of outliers:  [53, 52, 51]
        R      Lambda
 1   3.14819   3.15879
 2   2.97114   3.15143
 3   3.21044   3.14389
 4   2.83814   3.13616
 5   2.84416   3.12825
 6   2.87769   3.12013
 7   2.30345   3.11180
 8   2.33534   3.10324
 9   2.12480   3.09446
10   2.09054   3.08542

S-ESD与S-H-ESD

鉴于时间序列数据具有周期性（seasonal）、趋势性（trend），异常检测时不能作为孤立的样本点处理；故而Twitter的工程师提出了S- ESD (Seasonal ESD)与S-H-ESD (Seasonal Hybrid ESD)算法，将ESD扩展到时间序列数据。

STL分解

具体查看时间序列分解总结

STL (Seasonal-Trend decomposition procedure based on Loess) 为时序分解中一种常见的算法，基于LOESS将某时刻的数据$Y_v$分解为趋势分量（trend component）、季节性分量（seasonal component）和残差（remainder component）:

$$Y_v=T_v+S_v+R_{v}v=1,\dots ,N$$

由上到下依次为：原始时间序列和使用 STL 分解得到的季节变化部分、趋势变化部分以及残差部分。

STL分为内循环（inner loop）与外循环（outer loop），其中内循环主要做了趋势拟合与周期分量的计算。假定$T_v^{(k)}$、$Sv(k)$为内循环中第$k-1$次pass结束时的趋势分量、周期分量，初始时$T(k)v=0$；并有以下参数：

$n(i)$内层循环数
$n(o)$外层循环数
$n(p)$为一个周期的样本数
$n(s)$为Step 2中LOESS平滑参数
$n(l)$为Step 3中LOESS平滑参数
$n(t)$为Step 6中LOESS平滑参数

每个周期相同位置的样本点组成一个子序列（subseries），容易知道这样的子序列共有共有$n(p)$个，我们称其为cycle-subseries。内循环主要分为以下6个步骤：

Step 1: 去趋势（Detrending），减去上一轮结果的趋势分量，$Y_v-T_v^{(k)}$
Step 2: 周期子序列平滑（Cycle-subseries smoothing），用$LOESS(q=n_{n(s)},d=1)$对每个子序列做回归，并向前向后各延展一个周期；平滑结果组成temporary seasonal series，记为$C_v^{(k+1)},v=-n_{(p)}+1,\dots ,-N+n_{(p)}$
Step 3: 周期子序列的低通量过滤（Low-Pass Filtering），对上一个步骤的结果序列$C_v^{(k+1)}$依次做长度为$n(p)$、$n(p)$、3的滑动平均（moving average），然后做$LOESS(q=n_{n(l)},d=1)$回归，得到结果序列$L_v^{(k+1)},v=1,\dots ,N$；相当于提取周期子序列的低通量
Step 4: 去除平滑周期子序列趋势（Detrending of Smoothed Cycle-subseries），$S_v^{(k+1)}=C_v^{(k+1)}-L_v^{(k+1)}$
Step 5: 去周期（Deseasonalizing），减去周期分量，$Y_v–S_v^{(k+1)}$
Step 6: 趋势平滑（Trend Smoothing），对于去除周期之后的序列做$LOESS(q=n_{n(t)},d=1)$回归，得到趋势分量$T_v^{(k+1)}$。

外层循环主要用于调节robustness weight。如果数据序列中有outlier，则余项会较大。定义：

$$h=6\ast median(|R_v|)$$
对于位置为$v$的数据点，其robustness weight为：

$${\rho }_v=B(|R_v|/h)$$
其中B函数为bisquare函数：

$$B(u)=\{\begin{array}{c}(1-u^2{)}^{2}if0\le u<1\\ 0ifu\ge 0\end{array}$$
然后每一次迭代的内循环中，在Step 2与Step 6中做LOESS回归时，邻域权重（neighborhood weight）需要乘以${\rho }_v$，以减少outlier对回归的影响。STL的具体流程如下：

outer loop:
    计算robustness weight;
    inner loop:
        Step 1 去趋势;
        Step 2 周期子序列平滑;
        Step 3 周期子序列的低通量过滤;
        Step 4 去除平滑周期子序列趋势;
        Step 5 去周期;
        Step 6 趋势平滑;

为了使得算法具有足够的robustness，所以设计了内循环与外循环。特别地，当$n(i)$足够大时，内循环结束时趋势分量与周期分量已收敛；若时序数据中没有明显的outlier，可以将$n(o)$设为0。

Python的statsmodels实现了一个简单版的时序分解，通过加权滑动平均提取趋势分量，然后对cycle-subseries每个时间点数据求平均组成周期分量：

使用示例：

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose

# Generate some data
np.random.seed(0)
n = 1500
dates = np.array("2019-01-01", dtype=np.datetime64) + np.arange(n)
data = 12 * np.sin(2 * np.pi * np.arange(n) / 365) + np.random.normal(12, 2, 1500)
df = pd.DataFrame({"data": data}, index=dates)

# Reproduce the example in OP
seasonal_decompose(df, model="additive", period=1).plot()
plt.show()

S-ESD

STL将时间序列数据分解为趋势分量、周期分量和余项分量。想当然的解法——将ESD运用于STL分解后的余项分量中，即可得到时间序列上的异常点。但是，我们会发现在余项分量中存在着部分假异常点（spurious anomalies）。如下图所示：

在红色矩形方框中，向下突起点被误报为异常点。为了解决这种假阳性降低准确率的问题，S-ESD算法用中位数（median）替换掉趋势分量；余项计算公式如下：

$$R_X=X-S_X-\tilde{X}$$
其中，$X$为原时间序列数据，$S_X$为STL分解后的周期分量，$\tilde{X}$为$X$的中位数。

import numpy as np
import scipy.stats as stats
from statsmodels.tsa.seasonal import STL


def calculate_test_statistic(ts, hybrid=False):
    """
    Calculate the test statistic defined by being the top z-score in the time series.

    Args:
        ts (list or np.array): The time series to compute the test statistic.
        hybrid: A flag that determines the type of z-score. See the paper.

    Returns:
        tuple(int, float): The index of the top z-score and the value of the top z-score.

    """
    if hybrid:
        median = np.ma.median(ts)
        mad = np.ma.median(np.abs(ts - median))
        scores = np.abs((ts - median) / mad)
    else:
        scores = np.abs((ts - ts.mean()) / ts.std())
    max_idx = np.argmax(scores)
    return max_idx, scores[max_idx]


def calculate_critical_value(size, alpha):
    """
    Calculate the critical value with the formula given for example in
    https://en.wikipedia.org/wiki/Grubbs%27_test_for_outliers#Definition

    Args:
        size: The current size of the time series
        alpha (float): The significance level.

    Returns:
        float: The critical value for this test.

    """
    t_dist = stats.t.ppf(1 - alpha / (2 * size), size - 2)

    numerator = (size - 1) * t_dist
    denominator = np.sqrt(size**2 - size * 2 + size * t_dist**2)

    return numerator / denominator


def seasonal_esd(ts, periodicity=None, hybrid=False, max_anomalies=10, alpha=0.05):
    """
    Compute the Seasonal Extreme Studentized Deviate of a time series.
    The steps taken are first to to decompose the time series into STL
    decomposition (trend, seasonality, residual). Then, calculate
    the Median Absolute Deviate (MAD) if hybrid (otherwise the median)
    and perform a regular ESD test on the residual, which we calculate as:
                    R = ts - seasonality - MAD or median

    Note: The statsmodel library requires a seasonality to compute the STL
    decomposition, hence the parameter seasonality. If none is given,
    then it will automatically be calculated to be 20% of the total
    timeseries.

    Args:
    ts (list or np.array): The timeseries to compute the ESD.
    periodicity (int): Number of time points for a season.
    hybrid (bool): See Twitter's research paper for difference.
    max_anomalies (int): The number of times the Grubbs' Test will be applied to the ts.
    alpha (float): The significance level.

    Returns:
    list int: The indices of the anomalies in the timeseries.

    """
    if max_anomalies >= len(ts) / 2:
        raise ValueError(
            "The maximum number of anomalies must be less than half the size of the time series."
        )

    ts = np.array(ts)
    period = periodicity or int(
        0.2 * len(ts)
    )  # Seasonality is 20% of the ts if not given.
    stl = STL(ts, period=period, robust=True)
    decomposition = stl.fit()
    residual = ts - decomposition.seasonal - np.median(ts)
    outliers = generalized_esd(
        residual, max_anomalies=max_anomalies, alpha=alpha, hybrid=hybrid
    )
    return outliers


def generalized_esd(ts, max_anomalies=10, alpha=0.05, hybrid=False):
    """
    Compute the Extreme Studentized Deviate of a time series.
    A Grubbs Test is performed max_anomalies times with the caveat
       that each time the top value is removed. For more details visit
       http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm

    Args:
        ts (list or np.array): The time series to compute the ESD.
        max_anomalies (int): The number of times the Grubbs' Test will be applied to the ts.
        alpha (float): The significance level.
        hybrid: A flag that determines the type of z-score. See the paper.

    Returns:
        list int: The indices of the anomalies in the time series.

    """
    ts = np.ma.array(
        ts
    )  # A masked array needed to ignore outliers in subsequent ESD tests.
    test_statistics = []
    num_outliers = 0
    for curr in range(max_anomalies):
        test_idx, test_val = calculate_test_statistic(ts, hybrid=hybrid)
        critical_val = calculate_critical_value(len(ts) - curr, alpha)
        if test_val > critical_val:
            num_outliers = curr
        test_statistics.append(test_idx)
        ts[
            test_idx
        ] = (
            np.ma.masked
        )  # Mask this index so that we don't consider it in subsequent ESD tests.
    anomalous_indices = test_statistics[: num_outliers + 1] if num_outliers > 0 else []
    return anomalous_indices

import numpy as np

ts = np.random.random(100)
# Introduce artificial anomalies
ts[14] = 9
ts[83] = 10
outliers_indices = seasonal_esd(ts, periodicity=20, hybrid=True, max_anomalies=2)
for idx in outliers_indices:
    print(f"Anomaly index: {idx}, anomaly value: {ts[idx]}")

Anomaly index: 83, anomaly value: 10.0
Anomaly index: 14, anomaly value: 9.0

S-H-ESD

由于个别异常值会极大地拉伸均值和方差，从而导致S-ESD未能很好地捕获到部分异常点，召回率偏低。为了解决这个问题，S-H-ESD采用了更具鲁棒性的中位数与绝对中位差（Median Absolute Deviation, MAD）替换公式中的均值与标准差。MAD的计算公式如下：

$$MAD=median(|X_i-median(X)|)$$

三方包为pyculiarity

import matplotlib.pyplot as plt
import pandas as pd
from pyculiarity import detect_ts

# first run the models
twitter_example_data = pd.read_csv("raw_data.csv", usecols=["timestamp", "count"])
# 需要将源码中detect_anoms 67行添加sum()
results = detect_ts(
    twitter_example_data, max_anoms=0.05, alpha=0.001, direction="both", only_last=None
)

D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\pyculiarity\date_utils.py:37: FutureWarning: In a future version, `df.iloc[:, i] = newvals` will attempt to set the values inplace instead of always setting a new array. To retain the old behavior, use either `df[df.columns[i]] = newvals` or, if columns are non-unique, `df.isetitem(i, newvals)`
  indf.iloc[:,index] = column
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\pyculiarity\detect_anoms.py:58: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  if not isinstance(data.index, ps.Int64Index):
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\pyculiarity\detect_ts.py:256: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  all_anoms = all_anoms.append(anoms)
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\pyculiarity\detect_ts.py:257: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
  seasonal_plus_trend = seasonal_plus_trend.append(data_decomp)

# format the twitter data nicely
twitter_example_data["timestamp"] = pd.to_datetime(twitter_example_data["timestamp"])
twitter_example_data.set_index("timestamp", drop=True)

# make a nice plot
f, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(twitter_example_data["timestamp"], twitter_example_data["value"], "b")
ax[0].plot(results["anoms"].index, results["anoms"]["anoms"], "ro")
ax[0].set_title("Detected Anomalies")
ax[1].set_xlabel("Time Stamp")
ax[0].set_ylabel("Count")
ax[1].plot(results["anoms"].index, results["anoms"]["anoms"], "b")
ax[1].set_ylabel("Anomaly Magnitude")
plt.show()

移动平均/加权移动平均/指数加权移动平均

移动平均 moving average

给定一个时间序列和窗口长度$N$，moving average等于当前data point之前N个点（包括当前点）的平均值。不停地移动这个窗口，就得到移动平均曲线。

累加移动平均 cumulative moving average

设 { x i : i ≥ 1 } \{x_i:i\ge 1\} {xi​:i≥1}是观察到的数据序列。 累积移动平均线是所有数据的未加权平均值。 如果若干天的值是$x_1,\dots ,x_i$，那么$CM{A}_i=\frac{x_1+\cdots +x_i}{i}$

当有新的值$x_{i+1}$，那么累积移动平均为：

$$CM{A}_{i+1}=\frac{x_1+\cdots +x_i+x_{i+1}}{i+1}$$

加权移动平均 weighted moving average

加权移动平均值是先前$w$个数据的加权平均值。假设${\sum }_{j=0}^{w-1}weigh{t}_j=1$，其中$weigh{t}_j>0$，则加权移动平均值为$WM{A}_i={\sum }_{j=0}^{w-1}weigh{t}_j{x}_j$

指数加权移动平均 exponential weighted moving average

指数移动与移动平均有些不同：

并没有时间窗口，用的是从时间序列第一个data point到当前data point之间的所有点。

每个data point的权重不同，离当前时间点越近的点的权重越大，历史时间点的权重随着离当前时间点的距离呈指数衰减，从当前data point往前的data point，权重依次为$\alpha ,\alpha (1-\alpha ),\alpha (1-\alpha {)}^2,\dots ,\alpha (1-\alpha {)}^n$

该算法可以检测一个异常较短时间后发生另外一个异常的情况，异常持续一段时间后可能被判定为正常。

双指数平滑 double exponential smoothing (Holt-Linear)

假设 { Y t : t ≥ 1 } \{Y_t:t\ge 1\} {Yt​:t≥1}是观测数据序列，有两个与双指数平滑相关的方程：

$$S_t=\alpha Y_t+(1-\alpha )(S_t-1+b_t-1)$$
$$b_t=\beta (S_t-S_{t-1})+(1-\beta )b_t-1$$

其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子。

这里，初始值为$S_1=Y_1$，$b_1$有三种可能：

$$b_1=Y_2-Y_1$$
$$b_1=\frac{(Y_2-Y_1)+(Y_3-Y_2)+(Y_4-Y_3)}{3}=\frac{Y_4-Y_1}{3}$$
$$b_1=\frac{Y_n-Y_1}{n-1}$$

预测：

单个时刻：$F_{t+1}=S_t+b_t$

$m$个时刻： $F_{t+m}=S_t+m{b}_t$

三指数平滑 triple exponential smoothing （Holt-Winters）

multiplicative seasonality

设 { Y t : t ≥ 1 } \{Y_t:t\ge 1\} {Yt​:t≥1}是观察到的数据序列，则三指数平滑为：

$$S_t=\alpha \frac{Y_t}{c_{t-L}}+(1-\alpha )(S_{t-1}+b_{t-1}),overallsmoothing$$
$$b_t=\beta (S_t-S_{t-1})+(1-\beta )b_{t-1},trendsmoothing$$
$$c_t=\gamma \frac{Y_t}{S_t}+(1-\gamma )c_{t-L},seasonalsmoothing$$

$L$来表示季节频率，即一年中包含的季节数。例如，季度数据的$L=4$，月度数据的$L=12$

其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子，$\gamma \in [0,1]$是季节变化平滑因子。

预测$m$个时刻：$F_{t+m}=(S_t+m{b}_t)c_{(t-L+m)modL}$
初始值设定：

$$S_1=Y_1$$
$$b_0=\frac{(Y_{L+1}-Y_1)+(Y_{L+2}-Y_2)+\cdots +(Y_{L+L}-Y_L)}{L}$$
$$c_i=\frac{1}{N}\sum _{j=1}^N\frac{Y_{L(j-1)+i}}{A_j},\forall i\in 1,\dots ,L$$
$$A_j=\frac{{\sum }_{i=1}^L{Y}_{L(j-1)+i}}{L},\forall \in 1,\dots ,N$$

additive seasonality

设 { Y t : t ≥ 1 } \{Y_t:t\ge 1\} {Yt​:t≥1}是观察到的数据序列，则三指数平滑为：

$$S_t=\alpha (Y_t-c_{t-L})+(1-\alpha )(S_{t-1}+b_{t-1}),overallsmoothing$$
$$b_t=\beta (S_t-S_{t-1})+(1-\beta )b_{t-1},trendsmoothing$$
$$c_t=\gamma (Y_t-S_{t-1}-b_{t-1})+(1-\gamma )c_{t-L},seasonalsmoothing$$
$L$来表示季节频率，即一年中包含的季节数。例如，季度数据的$L=4$，月度数据的$L=12$

其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子，$\gamma \in [0,1]$是季节变化平滑因子。

预测$m$个时刻：$F_{t+m}=S_t+m{b}_t+c_{(t-L+m)modL}$

ARIMA 模型

具体查看arima
自回归移动平均模型（ARIMA）是一种设计上非常简单的方法，但其效果足够强大，可以预测信号并发现其中的异常。该方法的思路是从过去的几个数据点来生成下一个数据点的预测，在过程中添加一些随机变量（通常是添加白噪声）。以此类推，预测得到的数据点可以用来生成新的预测。很明显：它会使得后续预测信号数据更平滑。使用这种方法最困难的部分是选择差异数量、自回归数量和预测误差系数。另一个障碍是信号经过差分后应该是固定的。也就是说，这意味着信号不应该依赖于时间，这是一个比较显著的限制。

异常检测是利用离群点来建立一个经过调整的信号模型，然后利用t-统计量来检验该模型是否比原模型能更好的拟合数据。

该方法最受欢迎的实现是 R 语言中的 tsoutliers 包。在这种情况下，你可以找到适合信号的 ARIMA 模型，它可以检测出所有类型的异常。

神经网络

与CART方法一样，神经网络有两种应用方式：监督学习和无监督学习。我们处理的数据是时间序列，所以最适合的神经网络类型是 LSTM。如果构建得当，这种循环神经网络将可以建模实现时间序列中最复杂的依赖关系，包括高级的季节性依赖关系。如果存在多个时间序列相互耦合，该方法也非常有用。该领域还在研究中，可以参考这里，构建时序模型需要大量的工作。构建成功完成后，就可能在精确度方面取得优异的成绩。

时间序列异常检测方案

360时间序列检测机制

短期环比（SS）

对于时间序列来说，$T$ 时刻的数值对于 $T-1$ 时刻有很强的依赖性。比如流量在8:00很多，在8:01时刻的概率是很大的，但是07:01时刻对于8:01时刻影响不是很大。

首先，我们可以使用最近时间窗口（$T$）内的数据遵循某种趋势的现象来做文章。比如我们将 $T$ 设置为7，则我们取检测值（now_value）和过去7个（记为 $i$）点进行比较，如果大于阈值我们将 count 加1，如果 count 超过我们设置的 count_num，则认为该点是异常点。

$$count(\sum _{i=0}^T(nowvalue-i)>threshold)>coun{t}_{n}um$$
上面的公式涉及到$threshold$和$coun{t}_{n}um$两个参数， $coun{t}_{n}um$可以根据的需求进行设置，比如对异常敏感，可以设置$coun{t}_{n}um$小一些，而如果对异常不敏感，可以将$coun{t}_{n}um$设置的大一些。

$Threshold$可以采用动态阈值。业界关于动态阈值设置的方法有很多，通常阈值设置方法会参考过去一段时间内的均值、最大值以及最小值，我们也同样应用此方法。取过去一段时间（比如 $T$ 窗口）的平均值、最大值以及最小值，然后取 $max-avg$ 和 $avg-min$ 的最小值。之所以取最小值的原因是让筛选条件设置的宽松一些，让更多的值通过此条件，减少一些漏报的事件。

$$threshold=min(max-avg,avg-min)$$

长期环比（LS）

短期环比参考的是短期内的数据，而仅仅有短期内的数据是不够的，我们还需要参考更长时间内数据的总体走势。

通常使用一条曲线对该趋势进行拟合来反应曲线的走势，如果新的数据打破了这种趋势，使曲线变得不平滑，则该点就出现了异常。曲线拟合的方法有很多，比如回归、moving aver-age、…。在这我们使用 EWMA，即指数权重移动平均方法来拟合曲线。在 EWMA 中，下一点的平均值是由上一点的平均值，加上当前点的实际值修正而来。对于每一个 EWMA 值，每个数据的权重是不一样的，越近的数据将拥有越高的权重。有了平均值之后，我们就可以使用 3-sigma 理论来判断新的 input 是否超过了容忍范围。比较实际的值是否超出了这个范围就可以知道是否可以告警了。

同比（chain）

很多监控项都具有一定的周期性，其中以一天为周期的情况比较常见。为了将监控项的周期性考虑进去，我们选取了某个监控项过去14天的数据。对于某个时刻，将得到14个点可以作为参考值，我们记为 $x_i$，其中 $i=1,\dots ,14$。

我们先考虑静态阈值的方法来判断 input 是否异常（突增和突减）。如果 input 比过去14天同一时刻的最小值乘以一个阈值还小，就会认为该输入为异常点（突减）；而如果 input 比过去14天同一时刻的最大值乘以一个阈值还大，就会认为该输入为异常点（突增）。

静态阈值的方法是根据历史经验得出的值，实际中如何给 max_threshold 和 min_threshold 是一个需要讨论的话题。根据目前动态阈值的经验规则来说，取平均值是一个比较好的思路。

同比振幅（CA）

过去14天的历史曲线必然会比今天的曲线低很多。假设今天出了一个小故障，曲线下跌了，相对于过去14天的曲线仍然是高很多的。这样的故障，使用方法二就检测不出来，那么我们将如何改进我们的方法呢？

怎么计算 $t$ 时刻的振幅呢？我们使用$\frac{x(t)-x(t-1)}{x(t-1)}$来表示振幅。举个例子，例如 $t$ 时刻的流量为 900bit，$t-1$ 时刻的是 1000bit，那么可以计算出掉线人数是10%。如果参考过去的数据，我们会得到14个振幅值。使用14个振幅的绝对值作为标准，如果 amplitude threshold 小于 $m$ 时刻的振幅$\frac{m(t)-m(t-1)}{m(t-1)}$并且 $m$ 时刻的振幅大于0，则我们认为该时刻发生突增，而如果 $m$ 时刻的振幅大于 amplitude threshold 并且 $m$ 时刻的振幅小于0， 则认为该时刻发生突减。

$$amplitude=max[|\frac{x_i(t)-x_i(t-1)}{x_i(t-1)}|],x=1,2,\dots ,14$$

算法组合

上面介绍了四种方法，这四种方法里面，SS 和 LS 是针对非周期性数据的验证方法，而 chain 和 CA 是针对周期性数据的验证方法。那这四种方法应该如何选择和使用呢？下面我们介绍两种使用方法：

1、根据周期性的不同来选择合适的方法。这种方法需要首先验证序列是否具有周期性，如果具有周期性则进入左边分支的检测方法，如果没有周期性则选择进入右分支的检测方法。

上面涉及到了如何检测数据周期的问题，我们可以使用差分的方法来检测数据是否具有周期性。比如取最近两天的数据来做差分，如果是周期数据，差分后就可以消除波动，然后结合方差阈值判断的判断方法来确定数据的周期性。当然，如果数据波动范围比较大，可以在差分之前先对数据进行归一化（比如 z-score）。

2、不区分周期性，直接根据“少数服从多数”的方法来去检测，这种方法比较好理解，在此就不说明了。

美团外卖订单量预测异常报警模型

https://tech.meituan.com/2017/04/21/order-holtwinter.html

异常检测主要有两种策略：

异常驱动的异常检测（敏感性）：宁愿误报，也不能错过任何一个异常，这适用于非常重要的检测。简单概括，就是“宁可错杀一千，不能放过一个”。
预算驱动的异常检测（准确性）：这种策略的异常检测，从字面理解就是只有定量的一些预算去处理这些报警，那么只能当一定是某种问题时，才能将报警发送出来。

这两种策略不可兼容的。对于检测模型的改善，可以从两个方面入手，一是预测器的优化，二是比较器的优化。我们从这两个方面描述模型的改善。

预测器，就是用一批历史数据预测当前的数据。使用的历史数据集大小，以及使用的预测算法都会影响最终的预测效果。外卖订单量具有明显的周期性，同时相邻时刻的订单量数据也有很强的相关性，我们的目标，就是使用上面说的相关数据预测出当前的订单量。下面，我们分析几种常用的预测器实现。

同比环比预测器

同比环比是比较常用的异常检测方式，它是将当前时刻数据和前一时刻数据（环比）或者前一天同一时刻数据（同比）比较，超过一定阈值即认为该点异常。如果将不同日期、时刻的监控数据以矩阵方式存储，每一行表示一天内不同时刻的监控数据，每一列表示同一时刻不同日期的监控数据，那么存储矩阵如下图所示：

假如需要预测图中黄色数据，那么环比使用图中的蓝色数据作为预测黄点的源数据，同比使用图中红色数据作为预测黄点的源数据。

基线预测器

同比环比使用历史上的单点数据来预测当前数据，误差比较大。$t$时刻的监控数据，与$ t-1,t-2,\dots$时刻的监控数据存在相关性。同时，与$t-k,t-2k,\dots$时刻的数据也存在相关性（$k$为周期），如果能利用上这些相关数据对$t$时刻进行预测，预测结果的误差将会更小。

比较常用的方式是对历史数据求平均，然后过滤噪声，可以得到一个平滑的曲线（基线），使用基线数据来预测当前时刻的数据。该方法预测t时刻数据（图中黄色数据）使用到的历史数据如下图所示（图中红色数据）：

基线数据预测器广泛应用在业务大盘监控中，预测效果如图所示。从图中可以看出，基线比较平滑，在低峰期预测效果比较好，但是在外卖的午高峰和晚高峰预测误差比较大。

Holt-Winters预测器

同比环比预测到基线数据预测，使用的相关数据变多，预测的效果也较好。但是基线数据预测器只使用了周期相关的历史数据，没有使用上同周期相邻时刻的历史数据，相邻时刻的历史数据对于当前时刻的预测影响是比较大的。如外卖订单量，某天天气不好，很多用户不愿意出门，那么当天的外卖的订单量就会呈现整体的上涨，这种整体上涨趋势只能从同一周期相邻时刻的历史数据中预测出来。如图，预测图中黄色数据，如果使用上图中所有的红色数据，那么预测效果会更好。

本文使用了Holt-Winters来实现这一目标。Holt-Winters是三次指数滑动平均算法，它将时间序列数据分为三部分：残差数据$a(t)$，趋势性数据$b(t)$，季节性数据$s(t)$。使用Holt-Winters预测$t$时刻数据，需要$t$时刻前包含多个周期的历史数据。

外卖报警模型中的预测器

在外卖订单量异常检测中，使用Holt-Winters预测器实时预测下一分钟订单量，每次需要至少5天以上的订单量数据才能有较好的预测效果，数据量要求比较大。在实际的异常检测模型中，我们对Holt-Winters预测器进行了简化。预测器的趋势数据表示的是时间序列的总体变化趋势，如果以天为周期看待外卖的订单量时间序列，是没有明显的趋势性的，因此，我们可以去掉其中的趋势数据部分。

计算序列的周期性数据

时间序列的周期性数据不需要实时计算，按周期性更新即可，如外卖订单大盘监控，s(t)只需要每天更新一次即可。对于s(t)的计算，可以有多种方法，可以使用上面提到的Holt-Winters按公式计算出时间序列的周期性数据，或直接使用前一天的监控数据作为当天的周期数据（这两种方式都需要对输入序列进行预处理，保证算法的输入序列不含有异常数据）。也可以将历史数据做平均求出基线作为序列的周期性数据。

目前外卖订单中心报警模型采用的是Holt-Winters计算周期数据的方式。在将该模型推广到外卖其他业务线监控时，使用了计算基线数据作为周期数据的方式，这里简单对比一下两种方式的优劣。

1、使用Holt-Winters算法计算周期数据

优点：如果序列中含有周期性的陡增陡降点，Holt-Winters计算出的周期数据中会保留这些陡增陡降趋势，因此可以准确的预测出这些趋势，不会产生误报。比如外卖订单的提单数据，在每天的某个时刻都有一个定期陡降，使用该方式可以正确的预测出下降的趋势。如图所示，蓝色线是真实数据，棕色线是预测数据，在该时刻，棕色线准确的预测出了下降点。
缺点：需要对输入数据进行预处理，去除异常数据。如果输入序列中含有异常数据，使用Holt-Winters时可能会把这些异常数据计算到周期数据中，影响下一周期的预测从而产生误报（Holt-Winters理论上也只是滑动平均的过程，因此如果输入数据中含有比较大的异常数据时，存在这种可能性，实际应用中订单的报警模型也出现过这种误报）。

2、历史数据平均求基线

优点：计算出的周期数据比较平滑，不需要对输入序列进行预处理，计算过程中可以自动屏蔽掉异常数据的影响，计算过程简单。
缺点：周期数据比较平滑，不会出现陡增陡降点，因此对于周期性出现的陡增陡降不能很好的预测，出现误报。比如外卖活动的大盘（如图所示，红线是真实数据，黑线是预测数据），提前下单优惠在每天某个时刻会出现周期性陡降，使用该方式会出现误报。

残差数据实时预测

计算出周期数据后，下一个目标就是对残差数据的预测。实际监控数据与周期数据相减得到残差数据，对残差数据做一次滑动平均，预测出下一刻的残差，将该时刻的残差、周期数据相加即可得到该时刻的预测数据。残差序列的长度设为60，即可以得到比较准确的预测效果。

对于实时预测，使用的是当天的周期数据和前60分钟数据。最终的预测结果如图所示，其中蓝色线是真实数据，红色线是预测数据。

预测器预测出当前时刻订单量的预测值后，还需要与真实值比较来判断当前时刻订单量是否异常。一般的比较器都是通过阈值法，比如实际值超过预测值的一定比例就认为该点出现异常，进行报警。这种方式错误率比较大。在订单模型的报警检测中没有使用这种方式，而是使用了两个串联的Filter，只有当两个Fliter都认为该点异常时，才进行报警，下面简单介绍一下两个Filter的实现。

离散度Filter：根据预测误差曲线离散程度过滤出可能的异常点。一个序列的方差表示该序列离散的程度，方差越大，表明该序列波动越大。如果一个预测误差序列方差比较大，那么我们认为预测误差的报警阈值相对大一些才比较合理。离散度Filter利用了这一特性，取连续15分钟的预测误差序列，分为首尾两个序列（e1,e2），如果两个序列的均值差大于e1序列方差的某个倍数，我们就认为该点可能是异常点。
阈值Filter：根据误差绝对值是否超过某个阈值过滤出可能的异常点。利用离散度Filter进行过滤时，报警阈值随着误差序列波动程度变大而变大，但是在输入数据比较小时，误差序列方差比较小，报警阈值也很小，容易出现误报。所以设计了根据误差绝对值进行过滤的阈值Filter。阈值Filter设计了一个分段阈值函数y=f(x)，对于实际值x和预测值p，只有当|x-p|>f(x)时报警。实际使用中，可以寻找一个对数函数替换分段阈值函数，更易于参数调优。
最终的外卖订单异常报警模型结构图如图所示，每天会有定时Job从ETL中统计出最近10天的历史订单量，经过预处理模块，去除异常数据，经过周期数据计算模块得到周期性数据。对当前时刻预测时，取60分钟的真实数据和周期性数据，经过实时预测模块，预测出当前订单量。将连续15分钟的预测值和真实值通过比较器，判断当前时刻是否异常。