假設(shè)你正在進行運動時，使用數(shù)字設(shè)備監(jiān)測心率。你先跑了四分之一公里，然后走了十分鐘，接著又跑了四分之一公里里。可能你的心率情況與圖 (1) 中的時間序列相似。圖中展示了一段高心率、一段低心率，然后又回到高心率。時間序列的突然變化提示我們，你的活動狀態(tài)發(fā)生了重大變化。

圖 (1)

變點檢測是指在時間序列中發(fā)生了重大結(jié)構(gòu)性斷裂或者轉(zhuǎn)變的點，這些變化可能是由于數(shù)據(jù)生成、技術(shù)或消費者行為等外部因素造成的。檢測這些變點非常重要，因為它有助于我們理解和量化變化。我們需要及時準確地檢測這些變化并立即發(fā)出警報。

Change point detection (CPD) 被稱為變點檢測，其基本定義是在一個序列或過程中，當某個統(tǒng)計特性（分布類型、分布參數(shù)）在某時間點受系統(tǒng)性因素而非偶然因素影響發(fā)生變化，我們就稱該時間點為變點。變點識別即利用統(tǒng)計量或統(tǒng)計方法或機器學(xué)習(xí)方法將該變點位置估計出來。

CPD在金融、醫(yī)療保健和環(huán)境監(jiān)測等諸多領(lǐng)域都有著廣泛的應(yīng)用。其中，它在質(zhì)量控制過程中可以幫助識別產(chǎn)品或服務(wù)質(zhì)量的變化，也可以應(yīng)用于醫(yī)療診斷，幫助確定病人的健康狀況或疾病的變化。舉個例子，重癥監(jiān)護室（ICU）中的病人需要持續(xù)進行心率監(jiān)測，而及時發(fā)現(xiàn)心率的任何變化并迅速提醒醫(yī)生做出反應(yīng)對于病人的生命至關(guān)重要。

在CPD中，我們主要尋找時間序列中基本統(tǒng)計屬性（比如均值、方差或自相關(guān)性）發(fā)生明顯變化的點。雖然有多種算法可以檢測這些變化點，但一個重要的方面是要明確數(shù)據(jù)類型（即實時數(shù)據(jù)流還是離線數(shù)據(jù)），因為這將決定算法的選擇和發(fā)展。

算法取決于實時數(shù)據(jù)還是離線數(shù)據(jù)

CPD算法的運行方式取決于數(shù)據(jù)的類型，即實時數(shù)據(jù)或離線數(shù)據(jù)。對于離線數(shù)據(jù)，我們可以利用歷史數(shù)據(jù)來分析整個序列，這種情況下適用的是離線CPD。離線CPD涉及分析已經(jīng)收集的數(shù)據(jù)集，適用于歷史數(shù)據(jù)分析或檢測數(shù)據(jù)集中的異常情況。

然而，在實時環(huán)境中，我們需要快速檢測變點，而此時并沒有歷史數(shù)據(jù)可用。舉例來說，無人機每秒都會向家庭設(shè)備發(fā)送位置信息流，我們無法等待并長時間收集數(shù)據(jù)來檢測變化。在線視頻流也是一個例子。實時CPD需要在數(shù)據(jù)流到達時對其進行監(jiān)控，并立即做出響應(yīng)。這種情況常出現(xiàn)在許多實時應(yīng)用中，比如金融市場監(jiān)控、欺詐檢測或關(guān)鍵基礎(chǔ)設(shè)施監(jiān)控。

對于離線 CPD，我們將引入Ruptures Python 模塊。對于實時 CPD，我們將演示changefinder模塊。

生成數(shù)據(jù)

生成兩個時間序列來測試算法。

• 恒定方差 (ts1)：有十個數(shù)據(jù)段，方差恒定。
• 變化方差 (ts2)：有十個數(shù)據(jù)段，但方差變化。

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）顯示了時間序列。第一個時間序列中的變點比較容易發(fā)現(xiàn)，而第二個時間序列中的變點就比較難發(fā)現(xiàn)了。

圖(2)：恒定方差時間序列和變化方差時間序列

離線 - ruptures模塊

在離線分析中，我們能夠利用時間序列的歷史數(shù)據(jù)。對于 CPD，我們可以應(yīng)用線性回歸的概念。然而，如果存在變點，直線就無法很好地擬合數(shù)據(jù)，這時候分段線能夠更好地適應(yīng)數(shù)據(jù)。建立分段線的一種直觀算法是確定變點作為斷點。這種方法被稱為 精確線性時間（PELT）。

圖 (3.A) 和圖 (3.B) 解釋了PELT。在時間序列中（藍色顯示）存在一個變點和兩個分段。橙色線代表了回歸線，而橙色的垂直線表示了各點（用白色圓圈表示）到回歸線的距離。通過最小化所有數(shù)據(jù)點的距離之和來確定回歸線。

圖 (3)：剪枝后的精確線性時間（PELT）

在圖（3.B）中，分段線更適合數(shù)據(jù)。實際點到線條的距離和小于圖（3.A）中的距離之和。該算法通過從時間序列的左側(cè)滑動到右側(cè)來找到合適的變點，使得距離或誤差之和最小。

下面是用于搜索變點數(shù)量和位置的算法。C(.)代表距離或成本函數(shù)。我們還需要控制不要創(chuàng)建過多的線段，以防止對時間序列進行過度擬合。因此，b（β）項作為懲罰線段數(shù)量的參數(shù)，以防止搜索生成過多的線段。

圖片

該算法在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]處發(fā)現(xiàn)了變點。

# 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)

我們知道應(yīng)該有兩個變點 [1000, 2000, 5000]。但由于這些點前后的數(shù)據(jù)過于相似，PELT 無法發(fā)現(xiàn)這些差異。

圖 (5)：PELT 檢測到變化方差時間序列的一些變點

當使用 PELT 算法時，找到圖（4）以及圖（5）中的變化點可能需要相對較長的處理時間，特別是針對圖（5）。這樣可能無法滿足實時流數(shù)據(jù)的需求。因此，為了實時應(yīng)用，我們設(shè)計了名為changefinder 的 Python 模塊。

實時 CPD

時間序列可以用自回歸（AR）移動平均過程來描述。在AR模型中，下一個數(shù)據(jù)點是過去數(shù)據(jù)點的加權(quán)移動平均值，并且?guī)в须S機噪聲。具體而言，下式表示了AR模型，其中 θi 是過去 p 個數(shù)據(jù)點的權(quán)重。

圖片

如果存在一個變點，那么預(yù)計變點前后的自回歸過程將是不同的。 基于這種直覺，Yamanishi & Takeuchi 提出了順序貼現(xiàn)自回歸（SDAR）學(xué)習(xí)算法，這種直覺推動了他們的研究。SDAR 方法包括兩個學(xué)習(xí)階段。在第一個學(xué)習(xí)階段，它生成一個被稱為"異常分數(shù)"的中間分數(shù)。在第二個學(xué)習(xí)階段，它生成可以檢測變點的“變點分數(shù)”。以下是該算法的概述。

• 第 1 個 AR 模型 讀入一大塊數(shù)據(jù) t = 1, ..., N， 然后建立一個 AR（自動回歸）模型。然后產(chǎn)生一個 "異常得分"，即 AR 預(yù)測的 Xt 值與實際數(shù)據(jù) Xt 之間的差值。請注意，在這一步驟中只提取了 N 個數(shù)據(jù)點。由于它不使用整個歷史數(shù)據(jù)，因此是為在線數(shù)據(jù)流設(shè)置的。
• 第 1 次平滑處理上述異常得分將非常不穩(wěn)定，并產(chǎn)生錯誤信號。該算法為異常得分生成移動平均值 Yt，以平滑異常值。
• 第 2 個 AR 模型為 Yt 建立一個 AR 模型，并根據(jù)新建立的 AR 模型和 Yt 生成另一個 "異常得分"。
• 第 2 次平滑，同樣，新的異常得分也會出現(xiàn)波動。算法會生成移動平均值來平滑。如圖（6）所示，最終生成的分數(shù)稱為 "變點分數(shù)"。

這種算法不需要整個時間序列來檢測變點，因此大大減少了計算時間。

圖 (6)：順序貼現(xiàn)自動回歸（SDAR）學(xué)習(xí)算法

來研究兩種時間序列情況。

（1）恒定方差

適用于恒定方差時間序列 (ts1) 的前述代碼。Changefinder 需要三個參數(shù)：

• r：貼現(xiàn)率（0 至 1）。較高的貼現(xiàn)率會導(dǎo)致過去的時間序列迅速減少，意味著您可能不希望從過去的時間序列中學(xué)習(xí)。建議設(shè)置為 0.01。由于貼現(xiàn)率不是很敏感，設(shè)置為 0.01 或 0.05 都是可以的，可以自行嘗試。
• order：AR 模型階數(shù)
• smooth：用于計算平滑移動平均值的最近 N 個數(shù)據(jù)的大小。

在 changefinder 模塊中，我們對變點得分非常感興趣，它可以顯示時間序列是否突然偏離其常態(tài)。

# !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，第二張圖中的紅色時間序列是變點分數(shù)。發(fā)生變化的位置就是那些大的變點分數(shù)。

圖 (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 個點。這是因為這些點之前和之后的數(shù)據(jù)頻率過于相似，不符合差異條件。在運行代碼時，您可能會發(fā)現(xiàn)計算時間比破裂模塊的 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 個變點中的大多數(shù)變點

（2）變化方差

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

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

圖 (9) 顯示了第一張圖中的 ts2 和第二張圖中的變點分數(shù)。有三個高點。

圖(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），就會發(fā)現(xiàn)有三四個聚類。SDAR 算法可以檢測到這些主要變點。

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