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

《16非平稳时间序列突变检测的启发式分割算法(bg算法)matlab源代码》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、欢迎访问GreenSim团队主页→http://blog.sina.com.cn/greensim邮箱：greensim@163.com非平稳时间序列突变检测的启发式分割算法（BG算法）MATLAB源代码本源码的算法主要参考了下面参考文献：封国林,龚志强,董文杰等.基于启发式分割算法的气候突变检测研究[J].物理学报,2005,54(11):5494-5499。function[FLAG,AllT,AllTmax,AllPTmax]=BGA(X,P0,L0)%%非平稳时间序列突变检测的启发式分割算法%%输入参数列表%X待检测的数据，列向量存储%P0显著性水平门限值，低于此值的

2、不再分割%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,

3、p,PTmax]=Tseries(X);T=[1;T];ifPTmax

4、le1%设置死循环%%第三步：对每个段进行突变检测，能分割则分割，直到不能分割为止pos=find(FLAG>0);M=length(pos)-1;%当前子段数目form=1:Ms=pos(m);t=pos(m+1);L=length(SubX);ifL>=L0[T,Tmax,p,PTmax]=Tseries(SubX);T=[s;T];ifPTmax>=P0TC=TC+1;FLAG(s+p-1)=counter;AllT=[AllT;T];AllPTmax=[AllPTmax;PTmax];counter=counter+1;endendend%%第四步：返回输出数据ifT

5、C==0flag=FLAG;flag(1)=0;flag(N)=0;pos3=flag(pos2);FLAG=[pos2,pos3];returnend%%TC=0;%%endfunction[T,Tmax,p,PTmax]=Tseries(x)%%计算t检验统计序列的子函数%%参数列表%x时间序列，N×1列向量%Tt检验序列，N×1列向量%Tmaxt检验序列的最大值%pt检验序列最大值对应的下标%PTmaxTmax对应的统计显著性第4页欢迎访问GreenSim团队主页→http://blog.sina.com.cn/greensim邮箱：greensim@163.com%%

6、参数初始化N=length(x);T=zeros(N,1);%%以下是主循环，用于创建t检验序列fori=3:(N-2)%最左边以及最右边的两个点没有对应的t检验值（或者说，其值初始化为0）x1=x(1:i);%序列左边部分N1=length(x1);%左边序列的长度x2=x(i:N);%序列右边部分N2=length(x2);%右边序列的长度mean_x2=mean(x2);%右边部分的均值std_x2=std(x2);%右边部分的标准差%下面是计算合并偏差的公式，中英文文献里的这个公式略有不同，此处以英文文献为准SD=sqrt(1/N1+1/N2)*sqrt(((N1-1

7、)*std_x1^2+(N2-1)*std_x2^2)/(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+Tmax^2);%不完全B函数的下标PTmax=(1-betainc(c,Delta*Eta,