16非平稳时间序列突变检测的启发式分割算法(BG算法)MATLAB源代码.docx

欢迎访问GreenSim团队主页→ HYPERLINK /greensim /greensim 邮箱： HYPERLINK greensim@163.com greensim@163.com 第 PAGE \* MERGEFORMAT 4页 非平稳时间序列突变检测的启发式分割算法（BG算法）MATLAB源代码 本源码的算法主要参考了下面参考文献：封国林,龚志强,董文杰等.基于启发式分割算法的气候突变检测研究[J].物理学报,2005,54(11):5494-5499。 function [FLAG,AllT,AllTmax,AllPTmax]=BGA(X,P0,L0) %% 非平稳时间序列突变检测的启发式分割算法 %% 输入参数列表 % X 待检测的数据，列向量存储 % P0 显著性水平门限值，低于此值的不再分割 % L0 最小分割尺度，子段长度小于此值的不再分割 %% 输出参数列表 % FLAG 分割点标记，列向量存储，长度与X相同 % AllT 与分割点对应的全部t检验序列，其首位数字为起点坐标 % AllTmax 与分割点对应的全部t检验序列的最大值 % AllPTmax 与分割点对应的全部t检验序列对应的统计显著性 %% 第一步：变量初始化 N=length(X); FLAG=zeros(N,1); FLAG(1)=0.1; FLAG(N)=0.1; AllT=cell(0,0); AllTmax=cell(0,0); AllPTmax=cell(0,0); %% 第二步：产生第一个突变点，并对序列进行分割 [T,Tmax,p,PTmax]=Tseries(X); T=[1;T]; if PTmaxP0 flag=FLAG; flag(1)=0; flag(N)=0; pos3=flag(pos2); return end %记录输出数据 FLAG(p)=1; AllT=[AllT;T]; AllTmax=[AllTmax;Tmax]; AllPTmax=[AllPTmax;PTmax]; %以下为两个控制计数器 counter=2;%下一个突变点的序号 TC=0;%临时计数器 %% while 1%设置死循环 %% 第三步：对每个段进行突变检测，能分割则分割，直到不能分割为止 pos=find(FLAG0); M=length(pos)-1;%当前子段数目 for m=1:M s=pos(m); t=pos(m+1); L=length(SubX); if L=L0 [T,Tmax,p,PTmax]=Tseries(SubX); T=[s;T]; if PTmax=P0 TC=TC+1; FLAG(s+p-1)=counter; AllT=[AllT;T]; AllPTmax=[AllPTmax;PTmax]; counter=counter+1; end end end %% 第四步：返回输出数据 if TC==0 flag=FLAG; flag(1)=0; flag(N)=0; pos3=flag(pos2); FLAG=[pos2,pos3]; return end %% TC=0; %% end function [T,Tmax,p,PTmax]=Tseries(x) %% 计算t检验统计序列的子函数 %% 参数列表 % x 时间序列，N×1列向量 % T t检验序列，N×1列向量 % Tmax t检验序列的最大值 % p t检验序列最大值对应的下标 % PTmax Tmax对应的统计显著性 %% 参数初始化 N=length(x); T=zeros(N,1); %% 以下是主循环，用于创建t检验序列 for i=3:(N-2)%最左边以及最右边的两个点没有对应的t检验值（或者说，其值初始化为0） x1=x(1:i);%序列左边部分 N1=length(x1);%左边序列的长度 x2=x(i: