MannKendall突变检测算法MATLAB源码

MannKendall突变检测算法MATLAB源码 %% 以下是单边MannKendall突变检测算法 % GreenSim团队原创作品，转载请注明 % Email:greensim@163.com % GreenSim团队主页:[color=red]欢迎访问GreenSim——算法仿真团队?[url=][/url][/color] clear load data X=MJL; N=length(X); U=zeros(N-1,1); for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n S=S+sign(x(j)-x(k)); end end VarS=n*(n-1)*(2*n+5)/18; if S>0 Z=(S+1)/sqrt(VarS); elseif S==0 Z=0; else Z=(S-1)/sqrt(VarS); end U(t-1)=Z; end figure(1) plot(1:(N-1),U,'linewidth',1.5); hold on plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1); legend('统计量','0.05显著水平'); hold on plot(1:(N-1),0*ones(N-1,1),'-.','linewidth',1); plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1); plot(1:(N-1),-1.96*ones(N-1,1),':','linewidth',1); axis([1,N-1,-5,5]); xlabel('t (month)','FontName','TimesNewRoman','FontSize',12); ylabel('U统计量','FontName','TimesNewRoman','Fontsize',12); %grid on %计算趋势显著水平 Alpha=1-normcdf(U(end),0,1); disp('趋势的显著水平为'); disp(Alpha); %计算整体趋势的变化速率 Qi=zeros(N*(N-1)/2,1); counter=1; for k=1:(N-1) for j=(k+1):N Qi(counter)=(X(j)-X(k))/(j-k); counter=counter+1; end end Q=median(Qi); disp('整体趋势的变化速率为'); disp(Q); %% 以下是双边MannKendall突变检测算法 % GreenSim团队原创作品，转载请注明 % Email:greensim@163.com % GreenSim团队主页:[color=red]欢迎访问GreenSim——算法仿真团队 ?[url=][/url][/color] clear load data X=YJL; %计算UF统计量 N=length(X); UF=zeros(N-1,1); for t=2:N x=X(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)>x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UF(t-1)=Z; end %计算UB统计量 Y=flipud(X); UB=zeros(N-1,1); for t=2:N x=Y(1:t); S=0; n=length(x); for k=1:(n-1) for j=(k+1):n if x(j)>x(k) S=S+1; else S=S+0; end end end ES=n*(n+1)/4; VarS=n*(n-1)*(2*n+5)/72; Z=(S-ES)/sqrt(VarS); UB(t-1)=-Z; end %绘图 figure(2) plot(1:(N-1),UF,'r-','linewidth',1.5); hold on plot(1:(N-1),UB,'b-.','linewidth',1.5); plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1); axis([1,N-1,-4,4]); legend('UF统计量','UB统计量','0.05显著水平'); xlabel('t (year)','FontName','TimesNewRoman','FontSize',12); ylabel('统计量','FontName','TimesNewRoman','Fontsize',12); %grid on hold on plot(1:(N-1),0*ones(N-1,1),'-.','linewidth',1); plot(1:(N-1),1.96*ones(N-1,1),':','linewidth',1); plot(1:(N-1),-1.96*ones(N-1,1),':','linewidth',1);