非平稳时间序列突变检测 -- Bernaola Galvan分割算法

地学分析与算法

已于 2022-04-10 17:15:18 修改

阅读量9k

点赞数 9

分类专栏： Python 突变检测算法文章标签：突变检测 Bernaola Galvan 非平稳时间序列

于 2021-05-22 18:30:55 首次发布

本文链接：https://blog.csdn.net/qq_25102303/article/details/117168905

版权

Python 同时被 3 个专栏收录

11 篇文章 8 订阅

订阅专栏

算法

10 篇文章 0 订阅

订阅专栏

突变检测

2 篇文章 6 订阅

订阅专栏

非平稳时间序列突变检测 -- Bernaola Galvan分割算法

引言

非平稳序列是指包含趋势、季节性或周期性的序列，它可能只含有其中的一种成分，也可能是几种成分的组合，例如温度、降雨等数据。在一些研究中，如气候突变检测中，经常需要对气候数据进行突变检测。常用的突变检测方法有滑动t-检验、Cramer’s方法、Yamamoto方法、M-K突变检测方法、Pettitt方法、Bernaola Galvan分割算法等。今天主要介绍一种比较适用于非平稳时间序列的启发式突变检测方法–Bernaola Galvan分割算法（简称BG分割算法）。下图是BG算法在检测到的4个突变点。
在这里插入图片描述

去看原文

原理

对于包含N个样本的一时间序列数据X，通过第 $\leq i \leq N-1)$ 个数据将X分割为左右两个子序列X1（包含N1个样本）和X2（包含N2个样本），分别计算两个子序列的均值 $m_{i1}$ 、 $m_{i2}$ 和标准差 $std_{i1}$ 、 $std_{i2}$ ，以及t-检验统计值 $T (i)$ ，其计算公式为
$T(i)=abs(m_{i1}-m_{i2})/SD(i)$
其中 $S D (i)$ 第 $i$ 个点的合并偏差，
$SD(i)=sqrt(((N1-1)*std_{i1}^2 + (N2-1)*std_{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.95
  MIN_SUB_LENGTH = 100
  N = len(X)
  flag = [0] * N
  results = {}

  T, Tmax, Tmax_idx, PTmax = FindTmax(X)
  print(PTmax)
  if PTmax < P0:
    print("No find!")
    sys.exit()

  flag[Tmax_idx] = 1
  break_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}
  ...