Pettitt突变点检测

Pettitt突变点检测

• • 引言
  • 原理
  • 实现
  • 结果

引言

前期介绍了一些常用的突变检测方法，如MK突变检测法、Bernaola Galvan分割算法：Manner-Kendall(M-K)—突变检验；非平稳时间序列突变检测–BG算法。今天介绍一种非参数突变检测方法–Pettitt突变检测方法。下图是本文实现的Pettitt突变检测方法检测得到的突变点以及满足的置信水平，第一栏为原始序列数据及突变点，数据是两段均匀采样，第一段的均值0.9，第二段的均值为0.6，从图中可以看出明显的突变位置，第2栏为统计量Ut及显著性水平，可以看出该突变满足显著性水平p=0.001。

去看原文

原理

实现

在这里，采用Python实现上述原理，并通过测试用例数据进行验证，具体如下：

# -*- coding: utf-8 -*-
# @Author: 武辛
# @Email: geo_data_analysis@163.com
# @Note: 如有疑问，可加微信"wxid-3ccc"
# @All Rights Reserved!
​
import numpy as np
import pandas as pd
import sys, math
import matplotlib.pyplot as plt

def main():
  alpha = 0.05
  data = pd.read_csv("data_random.csv")
  Year = list(data["Year"])
  X = list(data["X"])
​
  Ut, max_idx, Kt = pettitt(X)
  print("The %dth sample is change point. Kt=%.2f, Year:%d, p=%.5f" % (max_idx + 1, Kt, Year[max_idx], CalP(Kt, len(X))))
  # print(Ut, max_idx, Kt)
  if CalP(Kt, len(X)) < alpha:
    print("The %dth sample is change point. Year:%d, p=%.9f" % (max_idx + 1, Year[max_idx], CalP(Kt, len(X))))
  else:
    print("Does not satisfy the significance level! P(Kt) = %.5f" % CalP(Kt, len(X)))
  plot(X, Ut, Kt, max_idx)
  
if __name__ == '__main__':
  main()

结果

第一栏为原始序列数据及突变点。第2栏为统计量Ut及显著性水平，可以看出该突变满足显著性水平p=0.001。

原文有源码，更多内容，请关注地学分析与算法。