- 文章概述
- 程序设计
- 学习总结
文章概述
本文提出了一种用于在线推理的贝叶斯变点检测算法。与回顾性分割不同，我们关注的是因果预测滤波；只考虑已经观察到的数据，生成序列中下一个不可见数据的精确分布。
程序设计
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
import numpy as np
import math
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
config={"font.family":'Times New Roman',"font.size":24,"mathtext.fontset":'stix'}
rcParams.update(config)
# input data
# 读取数据
filename=r'F:/Rpython/lp37/testdata3.xlsx'
df2=pd.read_excel(filename)#读取文件
years=df2['data1'].values.ravel()
d=df2['data2'].values.ravel()
d_series=pd.Series(d,index=years)
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Plot the input data. Note this works best in iPython notebook
# d_series.plot(title='Number of tropical storms in the North Atlantic per year')
# main changepoint detection algorithm
# single changepoint model
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
def step4(d):
n = len(d)
# dbar = sum(d)/float(n)
dbar = np.mean(d)
# dsbar = sum (d*d)/float(n)
dsbar = np.mean(np.multiply(d,d))
fac = dsbar-np.square(dbar)
summ = 0
summup = []
for z in range(n):
summ+=d[z]
summup.append(summ)
y = []
for m in range(n-1):
pos=m+1
mscale = 4*(pos)*(n-pos)
Q = summup[m]-(summ-summup[m])
U = -np.square(dbar*(n-2*pos) + Q)/float(mscale) + fac
y.append(-(n/float(2)-1)*math.log(n*U/2) - 0.5*math.log((pos*(n-pos))))
z, zz = np.max(y), np.argmax(y)
mean1 = sum(d[:zz+1])/float(len(d[:zz+1]))
mean2=sum(d[(zz+1):n])/float(n-1-zz)
return y, zz, mean1, mean2
# calling of main function
step_like = step4(d)
step_series = pd.Series(step_like[0],index=years[1:])
print(step_series)
print(step_series.shape)
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Plot the result. Note this works best in iPython notebook
# plt.figure();
# step_series.plot(title='Log likelihood of step change location in number of tropical storms per year')
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#开始绘图1
fig,ax=plt.subplots(figsize=(12,9),dpi=600)
ax.plot(years,d,linewidth=2,color='r')
# ax.grid(True)
ax.set_xlim((1958,2018))
ax.set_ylim((1200,2100))
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ax.set_xticks(np.arange(1958,2018.1,step=10))
ax.set_yticks(np.arange(1200,2100.1,step=100))
plt.savefig('F:/Rpython/lp37/plot94.png',dpi=800,bbox_inches='tight',pad_inches=0)
plt.show()
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#开始绘图2
fig,ax=plt.subplots(figsize=(12,9),dpi=600)
years2=range(1958,2016)
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ax.plot(years2,step_series,linewidth=2,color='r')
# ax.grid(True)
ax.set_xlim((1958,2018))
ax.set_ylim((-380,-370))
ax.set_xticks(np.arange(1958,2018.1,step=10))
ax.set_yticks(np.arange(-380,-369,step=2))
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plt.savefig('F:/Rpython/lp37/plot94.2.png',dpi=800,bbox_inches='tight',pad_inches=0)
plt.show()
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
变点检测是识别序列数据的生成参数的突变。作为一种在线或者离线信号处理工具，它已被证明在过程控制、脑电图分析、DNA分割、计量经济学和疾病人口统计学等应用中非常有用。
主要有几种较好的变点检测方法:
-
1：贝叶斯在线变点检测-Bayesian Online Changepoint Detection
本文提出了一种用于在线推理的贝叶斯变点检测算法。与回顾性分割不同，我们关注的是因果预测滤波；只考虑已经观察到的数据，生成序列中下一个不可见数据的精确分布。 -
2：findchangepts
变化点是信号的某些统计特性突然改变的采样点或时间瞬间。所讨论的属性可以是信号的平均值、方差或频谱特性等。为了找到信号变化点，findchangepts使用参数化全局方法。功能：
选择一个点并将信号分成两个部分。
计算每个部分所需统计特性的经验估计值。
在一个部分中的每一点上，测量属性偏离经验估计的程度。添加所有点的偏差。
将偏差部分添加到节中，以查找总残余误差。
改变分割点的位置，直到总残差达到最小值。
- 3：Recursive mean difference maximization
该算法原理很简单，我们要解决的问题是非平稳时间序列的划分，它由许多段组成用不同的平均值，以最大化的方式对于mleft和mright之间相邻seg差的平均值差。
接下来我们确定指针的位置t达到其最大值tmax，并计算tmax的统计显著性。
我们检查这个显著性是否超过了所选的阈值P0，通常取95%。如果是，则信号被切断此时分成两个子序列；否则信号不可分割。如果顺序被切断，程序对于每个切割创建的两个结果子序列中的每一个，递归地继续。
在接受新的切割点之前，我们还计算右侧新线段之间的t以及它的右邻居（通过前一次切割获得）和t在左边的新段和它的左邻居之间（也可通过先前的切割获得）并检查t的两个值是否具有超过P0的统计显著性。如果是的话，我们继续新的削减，否则我们不会削减。
因此，所有产生的片段都具有统计学意义不同。如果没有，则进程停止
