非平穩時間序列突變檢測 -- Bernaola Galvan分割算法

• • 引言
  • 原理
  • 實現
  • 結果

引言

非平穩序列是指包含趨勢、季節性或周期性的序列，它可能只含有其中的一種成分， 也可能是幾種成分的組合，例如溫度、降雨等數據。在一些研究中，如氣候突變檢測中，經常需要對氣候數據進行突變檢測。常用的突變檢測方法有滑動t-檢驗、Cramer’s方法、Yamamoto方法、M-K突變檢測方法、Pettitt方法、Bernaola Galvan分割算法等。今天主要介紹一種比較適用于非平穩時間序列的啟發式突變檢測方法–Bernaola Galvan分割算法（簡稱BG分割算法）。下圖是BG算法在檢測到的4個突變點。

去看原文

原理

對于包含N個樣本的一時間序列數據X，通過第$i(2\le i\le N?1)$個數據將X分割為左右兩個子序列X1（包含N1個樣本）和X2（包含N2個樣本），分別計算兩個子序列的均值$m_{i1}$、$m_{i2}$和標準差$st{d}_{i1}$、$st{d}_{i2}$，以及t-檢驗統計值$T(i)$，其計算公式為
$T(i)=abs(m_{i1}?m_{i2})/SD(i)$
其中$SD(i)$第$i$個點的合并偏差，
$SD(i)=sqrt(((N1?1)?st{d}_{i1}^2+(N2?1)?st{d}_{i2}^2)/(N1+N2?2))?sqrt(\frac{1}{N1}+\frac{1}{N2})$。
這樣從左到右依次計算每個數據的t-檢驗統計值$T(i)$，可以得到一個T序列，從中找到最大值$T_{max}$及其所對應的索引$j$，如果$T_{max}$的統計顯著性$P(T_{max})>=P0$（P0為給定的參數），那么便可對序列X在第$j$個樣本處進行分割，也就是突變點。$P(T_{max})$的計算可近似為

同樣的，當分割完后，可以對分割后的兩個子序列重復進行上述操作，直到不可分割為止，便可以得到所有的突變點。

實現

在這里，采用Python實現上述原理，并通過測試用例數據（用例數據隨機生成）進行驗證，具體如下：

# -*- coding: utf-8 -*- # @Author: 武辛 # @Email: geo_data_analysis@163.com # @Note: 如有疑問，可加微信"wxid-3ccc" # @All Rights Reserved! ? import pandas as pd import numpy as np import scipy.special import math, sys import random import matplotlib.pyplot as plt np.random.seed(0) ? def main(): ?data = pd.read_csv("data.csv")X = list(data["X"]) ?P0 = 0.95MIN_SUB_LENGTH = 100N = len(X)flag = [0] * Nresults = {} ?T, Tmax, Tmax_idx, PTmax = FindTmax(X)print(PTmax)if PTmax < P0:print("No find!")sys.exit() ?flag[Tmax_idx] = 1break_idx = [0, Tmax_idx, len(X)-1]results[Tmax_idx] = {"Tmax":Tmax, "T":T, "PTmax":PTmax, "start_idx":0, "end_idx":len(X)-1, "break_order":1}...

結果

第一欄為原始序列數據，可以看出，存在4個比較明顯的突變。第2-5欄為BG算法計算t-檢驗統計值T的結果，在第2欄圖中，紅色線部分找到最大Tmax，即在第799個樣本數據處存在第一個突變點，第3-5欄圖中分別在7000、4999、2799處找到突變點。

原文有源碼，更多內容，請關注地學分析與算法。

總結

以上是生活随笔為你收集整理的非平稳时间序列突变检测 -- Bernaola Galvan分割算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。