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

1、精品非平稳时间序列突变检测的启发式分割算法( BG 算法) MATLAB源代码本源码的算法主要参考了下面参考文献：封国林 ,龚志强 ,董文杰等 .基于启发式分割算法的气候突变检测研究J物理学报，2005,54(11):5494-5499。function FLAG,AllT,AllTmax,AllPTmax=BGA(X,P0,L0)% 非平稳时间序列突变检测的启发式分割算法% 输入参数列表% X待检测的数据，列向量存储% P0显著性水平门限值，低于此值的不再分割% L0最小分割尺度，子段长度小于此值的不再分割% 输出参数列表% FLAG分割点标记，列向量存储，长度与 X 相同% AllT与分割

2、点对应的全部 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 PTmax0);M=length(pos)-1;% 当前子段数目for m=1:

3、Ms=pos(m);t=pos(m+1);L=length(SubX);if L=L0T,Tmax,p,PTmax=Tseries(SubX);T=s;T;if PTmax=P0TC=TC+1;FLAG(s+p-1)=counter;AllT=AllT;T;AllPTmax=AllPTmax;PTmax;counter=counter+1;endendend% 第四步：返回输出数据if TC=0flag=FLAG;flag(1)=0;flag(N)=0;pos3=flag(pos2);FLAG=pos2,pos3;returnend%TC=0;%endfunction T,Tmax,p,PTm

4、ax=Tseries(x)% 计算 t 检验统计序列的子函数% 参数列表% x时间序列，N X1列向量% Tt检验序列，N X1列向量% Tmax t 检验序列的最大值% pt 检验序列最大值对应的下标% PTmax Tmax 对应的统计显著性% 参数初始化N=length(x);T=zeros(N,1);% 以下是主循环，用于创建 t 检验序列for i=3:(N-2)% 最左边以及最右边的两个点没有对应的 t 检验值(或者说，x1=x(1:i);% 序列左边部分N1=length(x1);% 左边序列的长度x2=x(i:N);% 序列右边部分N2=length(x2);% 右边序列的长度m

5、ean_x2=mean(x2);% 右边部分的均值std_x2=std(x2);% 右边部分的标准差% 下面是计算合并偏差的公式， 中英文文献里的这个公式略有不同，准SD=sqrt(1/N1+1/N2)*sqrt(N1-1)*std_x1A2+(N2-1)*std_x2A2)/(N1+N2-2);T(i)=abs(mean_x1-mean_x2)/SD);end% 计算其它三个输出参数Tmax=max(T);%t 检验序列的最大值pos=find(T=Tmax);Eta=4.19*log(N)-11.54;% 计算 PTmax 用的参数Delta=0.40;% 计算 PTmax 用的参数v=N-2;% 计算 PTmax 用的参数c=v/(v+TmaxA2);%不完全B函数的下标PTmax=(1-betainc(c,Delta*Eta,Delta);% 调用不完全 beta 函数其值