用于时间序列中的变点检测算法

图 (1)

变点检测是指在时间序列中发生了重大结构性断裂或者转变的点，这些变化可能是由于数据生成、技术或消费者行为等外部因素造成的。检测这些变点非常重要，因为它有助于我们理解和量化变化。我们需要及时准确地检测这些变化并立即发出警报。

Change point detection (CPD) 被称为变点检测，其基本定义是在一个序列或过程中，当某个统计特性（分布类型、分布参数）在某时间点受系统性因素而非偶然因素影响发生变化，我们就称该时间点为变点。变点识别即利用统计量或统计方法或机器学习方法将该变点位置估计出来。

CPD在金融、医疗保健和环境监测等诸多领域都有着广泛的应用。其中，它在质量控制过程中可以帮助识别产品或服务质量的变化，也可以应用于医疗诊断，帮助确定病人的健康状况或疾病的变化。举个例子，重症监护室（ICU）中的病人需要持续进行心率监测，而及时发现心率的任何变化并迅速提醒医生做出反应对于病人的生命至关重要。

在CPD中，我们主要寻找时间序列中基本统计属性（比如均值、方差或自相关性）发生明显变化的点。虽然有多种算法可以检测这些变化点，但一个重要的方面是要明确数据类型（即实时数据流还是离线数据），因为这将决定算法的选择和发展。

算法取决于实时数据还是离线数据

CPD算法的运行方式取决于数据的类型，即实时数据或离线数据。对于离线数据，我们可以利用历史数据来分析整个序列，这种情况下适用的是离线CPD。离线CPD涉及分析已经收集的数据集，适用于历史数据分析或检测数据集中的异常情况。

然而，在实时环境中，我们需要快速检测变点，而此时并没有历史数据可用。举例来说，无人机每秒都会向家庭设备发送位置信息流，我们无法等待并长时间收集数据来检测变化。在线视频流也是一个例子。实时CPD需要在数据流到达时对其进行监控，并立即做出响应。这种情况常出现在许多实时应用中，比如金融市场监控、欺诈检测或关键基础设施监控。

对于离线 CPD，我们将引入Ruptures Python 模块。对于实时 CPD，我们将演示changefinder模块。

生成数据

生成两个时间序列来测试算法。

• 恒定方差 (ts1)：有十个数据段，方差恒定。
• 变化方差 (ts2)：有十个数据段，但方差变化。

import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt

# Example 1: 恒定方差
ts1 = []
mu, sigma, seg = 0.0, 1.0, 1000
for i in range(10):
    ts = np.random.normal(mu, sigma, seg) + np.random.randint(low=-10, high=10)
    ts1 = np.append(ts1,ts, axis=0)

plt.figure(figsize=(16,4))
plt.plot(ts1)


# Example 2: 变化方差
ts2 = []
mu, sigma, seg = 0.0, 1.0, 1000
for i in range(10):
    sig = np.random.randint(low=1, high=50)
    ts = np.random.normal(mu, sigma * sig, seg) 
    ts2 = np.append(ts2,ts, axis=0)
    
plt.figure(figsize=(16,4))
plt.plot(ts2)

图（2）显示了时间序列。第一个时间序列中的变点比较容易发现，而第二个时间序列中的变点就比较难发现了。

图(2)：恒定方差时间序列和变化方差时间序列

离线 - ruptures模块

在离线分析中，我们能够利用时间序列的历史数据。对于 CPD，我们可以应用线性回归的概念。然而，如果存在变点，直线就无法很好地拟合数据，这时候分段线能够更好地适应数据。建立分段线的一种直观算法是确定变点作为断点。这种方法被称为 精确线性时间（PELT）。

图 (3.A) 和图 (3.B) 解释了PELT。在时间序列中（蓝色显示）存在一个变点和两个分段。橙色线代表了回归线，而橙色的垂直线表示了各点（用白色圆圈表示）到回归线的距离。通过最小化所有数据点的距离之和来确定回归线。

图 (3)：剪枝后的精确线性时间（PELT）

在图（3.B）中，分段线更适合数据。实际点到线条的距离和小于图（3.A）中的距离之和。该算法通过从时间序列的左侧滑动到右侧来找到合适的变点，使得距离或误差之和最小。

下面是用于搜索变点数量和位置的算法。C(.)代表距离或成本函数。我们还需要控制不要创建过多的线段，以防止对时间序列进行过度拟合。因此，b（β）项作为惩罚线段数量的参数，以防止搜索生成过多的线段。

该算法在Python 模块ruptures中编码。

首次使用，需要用 pip install ruptures 安装。

(1)恒定方差

# !pip install ruptures
import ruptures as rpt

# Detect the change points
algo1 = rpt.Pelt(model="rbf").fit(ts1)
change_location1 = algo1.predict(pen=10)

# Point the change points:
def plot_change_points(ts,ts_change_loc):
    plt.figure(figsize=(16,4))
    plt.plot(ts)
    for x in ts_change_loc:
        plt.axvline(x,lw=2, color='red')

plot_change_points(ts1,change_location1)

图 (4) 显示，算法检测到了所有十个变化点。

图 (4)：检测到恒定方差时间序列的所有十个变点

当方差随时间变化时，CPD 是否仍然有效。

（2）变化方差

图 (5) 显示 PELT 算法在[3000, 4000, 6005, 7000, 8000, 8995, 10000]处发现了变点。

# detect the change points #
algo2 = rpt.Pelt(model="rbf").fit(ts2)
change_location2 = algo2.predict(pen=10)
change_location2

# Plot the change points #
plot_change_points(ts2,change_location2)

我们知道应该有两个变点 [1000, 2000, 5000]。但由于这些点前后的数据过于相似，PELT 无法发现这些差异。

图 (5)：PELT 检测到变化方差时间序列的一些变点

当使用 PELT 算法时，找到图（4）以及图（5）中的变化点可能需要相对较长的处理时间，特别是针对图（5）。这样可能无法满足实时流数据的需求。因此，为了实时应用，我们设计了名为changefinder 的 Python 模块。

实时 CPD

时间序列可以用自回归（AR）移动平均过程来描述。在AR模型中，下一个数据点是过去数据点的加权移动平均值，并且带有随机噪声。具体而言，下式表示了AR模型，其中 θi 是过去 p 个数据点的权重。

如果存在一个变点，那么预计变点前后的自回归过程将是不同的。 基于这种直觉，Yamanishi & Takeuchi 提出了顺序贴现自回归（SDAR）学习算法，这种直觉推动了他们的研究。SDAR 方法包括两个学习阶段。在第一个学习阶段，它生成一个被称为"异常分数"的中间分数。在第二个学习阶段，它生成可以检测变点的“变点分数”。以下是该算法的概述。

• 第 1 个 AR 模型 读入一大块数据 t = 1, ..., N， 然后建立一个 AR（自动回归）模型。然后产生一个 "异常得分"，即 AR 预测的 Xt 值与实际数据 Xt 之间的差值。请注意，在这一步骤中只提取了 N 个数据点。由于它不使用整个历史数据，因此是为在线数据流设置的。
• 第 1 次平滑处理上述异常得分将非常不稳定，并产生错误信号。该算法为异常得分生成移动平均值 Yt，以平滑异常值。
• 第 2 个 AR 模型为 Yt 建立一个 AR 模型，并根据新建立的 AR 模型和 Yt 生成另一个 "异常得分"。
• 第 2 次平滑，同样，新的异常得分也会出现波动。算法会生成移动平均值来平滑。如图（6）所示，最终生成的分数称为 "变点分数"。

这种算法不需要整个时间序列来检测变点，因此大大减少了计算时间。

图 (6)：顺序贴现自动回归（SDAR）学习算法

来研究两种时间序列情况。

（1）恒定方差

适用于恒定方差时间序列 (ts1) 的前述代码。Changefinder 需要三个参数：

• r：贴现率（0 至 1）。较高的贴现率会导致过去的时间序列迅速减少，意味着您可能不希望从过去的时间序列中学习。建议设置为 0.01。由于贴现率不是很敏感，设置为 0.01 或 0.05 都是可以的，可以自行尝试。
• order：AR 模型阶数
• smooth：用于计算平滑移动平均值的最近 N 个数据的大小。

在 changefinder 模块中，我们对变点得分非常感兴趣，它可以显示时间序列是否突然偏离其常态。

# !pip install changefinder
import changefinder

def findChangePoints(ts, r, order, smooth):
    '''
       r: Discounting rate
       order: AR model order
       smooth: smoothing window size T
    '''
    cf = changefinder.ChangeFinder(r=r, order=order, smooth=smooth)
    ts_score = [cf.update(p) for p in ts]
    plt.figure(figsize=(16,4))
    plt.plot(ts)
    plt.figure(figsize=(16,4))
    plt.plot(ts_score, color='red')
    return(ts_score)
    
ts_score1 = findChangePoints(ts1, r = 0.01, order = 3, smooth = 5)

图 (7) 显示了第一张图中的 ts1，第二张图中的红色时间序列是变点分数。发生变化的位置就是那些大的变点分数。

图 (7)：针对恒定方差时间序列的 SDAR 算法的变点得分

在此，我打印出了前 20 名的位置（您可以选择更多）。

ts_change_loc1 = pd.Series(ts_score1).nlargest(20)
ts_change_loc1 = ts_change_loc1.index
np.sort(ts_change_loc1)

array([ 11, 12, 13, 42, 43, 159, 160, 4001, 4008,
      5001, 5007, 5008, 6007, 6008, 7000, 7001,
      9000, 9001, 9007, 9008])

总的来说，考虑到这是实时操作，该算法做得还不错。它能检测到许多变点，但却漏掉了 1000、2000 和 8000 个点。这是因为这些点之前和之后的数据频率过于相似，不符合差异条件。在运行代码时，您可能会发现计算时间比破裂模块的 PELT 方法要少得多。

在图 (8) 中绘制带有变点的时间序列。

def plot_change_points(ts,ts_change_loc):
    plt.figure(figsize=(16,4))
    plt.plot(ts)
    for x in ts_change_loc:
        plt.axvline(x,lw=2, color='red')
        
plot_change_points(ts1,ts_change_loc1)

图(8)：SDAR 算法检测到了恒定方差时间序列 10 个变点中的大多数变点

（2）变化方差

变化方差时间序列中的变点很难找到。如图（9）所示，ts2 看起来像三到四个簇，而不是我们设计的十个簇。

ts_score2 = findChangePoints(ts2, r = 0.01, order = 3, smooth = 5)

图 (9) 显示了第一张图中的 ts2 和第二张图中的变点分数。有三个高点。

图(9)：变化方差时间序列的 SDAR 算法变点得分

打印出前 20 名的位置。

ts_change_loc2 = pd.Series(ts_score2).nlargest(20)
ts_change_loc2 = ts_change_loc2.index
np.sort(ts_change_loc2)

array([ 11, 12, 13, 19, 20, 21, 22, 80, 81, 107, 
        108, 117, 4003, 4004, 4010, 4011, 8000,
        8001, 8007, 8008])*

该方法确定了 1、4000、8000，但错过了 1000、2000、3000、5000、6000、7000 和 9000 的变点。

plot_change_points(ts2,ts_change_loc2)

如果我们直观地观察图（10），就会发现有三四个聚类。SDAR 算法可以检测到这些主要变点。

图(10)：SDAR 算法检测变化方差时间序列的主要变点