用於時間序列中的變點偵測演算法，你學會了嗎？

假設你正在進行運動時，使用數位設備監測心率。你先跑了四分之一公里，然後走了十分鐘，接著又跑了四分之一公里裡。可能你的心率狀況與圖(1) 中的時間序列相似。圖中展示了一段高心率、一段低心率，然後又回到高心率。時間序列的突然變化提示我們，你的活動狀態有了重大變化。

圖(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)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.

圖（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)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.

圖(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)
1.
2.
3.
4.
5.
6.
7.

我們知道應該有兩個變點[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)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.

圖(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)
1.
2.
3.

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

總的來說，考慮到這是即時操作，演算法做得還不錯。它能偵測到許多變點，但卻漏掉了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)
1.
2.
3.
4.
5.
6.
7.

圖(8)：SDAR 演算法偵測到了恆定方差時間序列10 個變點中的大多數變點

（2）變化方差

變化方差時間序列中的變點很難找到。如圖（9）所示，ts2 看起來像三到四個簇，而不是我們設計的十個簇。

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

圖(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)
1.
2.
3.

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

該方法確定了1、4000、8000，但錯過了1000、2000、3000、5000、6000、7000 和9000 的變點。

plot_change_points(ts2,ts_change_loc2)
1.

如果我們直觀地觀察圖（10），就會發現有三、四個聚類。 SDAR 演算法可以偵測到這些主要變點。

圖(10)：SDAR 演算法偵測變化方差時間序列的主要變點