1、非平穩(wěn)時間序列突變檢測的啟發(fā)式分割算法（BG算法）MATLAB源代碼本源碼的算法主要參考了下面參考文獻：封國林,龔志強,董文杰等.基于啟發(fā)式分割算法的氣候突變檢測研究J.物理學(xué)報,2005,54(11):5494-5499。function FLAG,AllT,AllTmax,AllPTmax=BGA(X,P0,L0)% 非平穩(wěn)時間序列突變檢測的啟發(fā)式分割算法% 輸入?yún)?shù)列表% X 待檢測的數(shù)據(jù)，列向量存儲% P0 顯著性水平門限值，低于此值的不再分割% L0 最小分割尺度，子段長度小于此值的不再分割% 輸出參數(shù)列表% FLAG 分割點標記，列向量存儲，長度與X相同% AllT 與分割點對應(yīng)的

2、全部t檢驗序列，其首位數(shù)字為起點坐標% AllTmax 與分割點對應(yīng)的全部t檢驗序列的最大值% AllPTmax 與分割點對應(yīng)的全部t檢驗序列對應(yīng)的統(tǒng)計顯著性% 第一步：變量初始化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);% 第二步：產(chǎn)生第一個突變點，并對序列進行分割T,Tmax,p,PTmax=Tseries(X);T=1;T;if PTmax0); M=length(pos)-1;%當前子段數(shù)目 for m=1:M s=pos

3、(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% 第四步：返回輸出數(shù)據(jù) if TC=0 flag=FLAG; flag(1)=0; flag(N)=0; pos3=flag(pos2); FLAG=pos2,pos3; return end% TC=0;%endfunc

4、tion T,Tmax,p,PTmax=Tseries(x)% 計算t檢驗統(tǒng)計序列的子函數(shù)% 參數(shù)列表% x 時間序列，N1列向量% T t檢驗序列，N1列向量% Tmax t檢驗序列的最大值% p t檢驗序列最大值對應(yīng)的下標% PTmax Tmax對應(yīng)的統(tǒng)計顯著性% 參數(shù)初始化N=length(x);T=zeros(N,1);% 以下是主循環(huán)，用于創(chuàng)建t檢驗序列for i=3:(N-2)%最左邊以及最右邊的兩個點沒有對應(yīng)的t檢驗值（或者說，其值初始化為0） x1=x(1:i);%序列左邊部分 N1=length(x1);%左邊序列的長度 x2=x(i:N);%序列右邊部分 N2=length

5、(x2);%右邊序列的長度 mean_x2=mean(x2);%右邊部分的均值 std_x2=std(x2);%右邊部分的標準差 %下面是計算合并偏差的公式，中英文文獻里的這個公式略有不同，此處以英文文獻為準 SD=sqrt(1/N1+1/N2)*sqrt(N1-1)*std_x12+(N2-1)*std_x22)/(N1+N2-2); T(i)=abs(mean_x1-mean_x2)/SD);end% 計算其它三個輸出參數(shù)Tmax=max(T);%t檢驗序列的最大值pos=find(T=Tmax);Eta=4.19*log(N)-11.54;%計算PTmax用的參數(shù)Delta=0.40;%計算PTmax用的參數(shù)v=N-2;%計算PTmax用的參數(shù)c