《MK突变检测程序.docx》由会员分享，可在线阅读，更多相关《MK突变检测程序.docx（7页珍藏版）》请在三一办公上搜索。
1、MK突变检测程序MannKendall突变检测算法MATLAB源码 题目：MannKendall突变检测算法MATLAB源码 % 以下是单边MannKendall突变检测算法 % GreenSim团队原创作品，转载请注明 % Email:greensim % GreenSim团队主页： % color=red欢迎访问GreenSim算法仿真团队url= 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+
2、sign(x(j)-x(k); end end VarS=n*(n-1)*(2*n+5)/18; if S0 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
3、); 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); %计算整体趋
4、势的变化速率 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 % GreenSim团队主页： % color=red欢迎访问GreenSim算法仿真团队url= clear load data X=
5、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
6、(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-
7、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); %最近写论文
8、需要用到MK检验法，网上收集到大量的matlab代码，但是没有一个代码能够 %完全正确运行或者分析信息不全，结合多位网友编写的MK检验法，经过我的改编，顺利得到 %正确的运行结果，谢谢各位网友，希望对有需要的盆友有帮助 % Mann-Kendall突变检测 % 数据序列y % 结果序列UFk，UBk2 %- %读取excel中的数据，赋给矩阵y %获取y的样本数 %A为时间和径流数据列 A=xlswrite(数据.xls) x=A(:,1);%时间序列 y=A(:,2);%径流数据列 N=length(y); n=length(y); % 正序列计算- % 定义累计量序列Sk，长度=y，初始值
9、=0 Sk=zeros(size(y); % 定义统计量UFk，长度=y，初始值=0 UFk=zeros(size(y); % 定义Sk序列元素s s = 0; % i从2开始，因为根据统计量UFk公式，i=1时，Sk(1)、E(1)、Var(1)均为0 % 此时UFk无意义，因此公式中，令UFk(1)=0 for i=2:n for j=1:i if y(i)y(j) s=s+1; else s=s+0; end; end; Sk(i)=s; E=i*(i-1)/4; % Sk(i)的均值 Var=i*(i-1)*(2*i+5)/72; % Sk(i)的方差 UFk(i)=(Sk(i)-E)
10、/sqrt(Var); end; % -正序列计算end % 逆序列计算- % 构造逆序列y2，长度=y，初始值=0 y2=zeros(size(y); % 定义逆序累计量序列Sk2，长度=y，初始值=0 Sk2=zeros(size(y); % 定义逆序统计量UBk，长度=y，初始值=0 UBk=zeros(size(y); % s归0 s=0; % 按时间序列逆转样本y % 也可以使用y2=flipud(y);或者y2=flipdim(y,1); for i=1:n y2(i)=y(n-i+1); end; % i从2开始，因为根据统计量UBk公式，i=1时，Sk2(1)、E(1)、Var
11、(1)均为0 % 此时UBk无意义，因此公式中，令UBk(1)=0 for i=2:n for j=1:i if y2(i)y2(j) s=s+1; else s=s+0; end; end; Sk2(i)=s; E=i*(i-1)/4; % Sk2(i)的均值 Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差 % 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1， % 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反 % 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i)/sqrt(Var(
12、i) % 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势 UBk(i)=0-(Sk2(i)-E)/sqrt(Var); end; % -逆序列计算end % 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量 % 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此 % 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，做图用 UBk2=zeros(size(y); % 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1); for i=1:n UBk2(i)=UBk(n-i+1); end; % 做突变检测图时，使用UF
13、k和UBk2 % 写入目标xls文件：f:test2.xls % 目标表单：Sheet1 % 目标区域：UFk从A1开始，UBk2从B1开始 xlswrite(f:test2.xls,UFk,Sheet1,A1); xlswrite(f:test2.xls,UBk2,Sheet1,B1); figure(3)%画图 plot(x,UFk,r-,linewidth,1.5); hold on plot(x,UBk2,b-.,linewidth,1.5); plot(x,1.96*ones(N,1),:,linewidth,1); axis(min(x),max(x),-5,5); legend(UF统计量,UB统计量,0.05显著水平); xlabel(t (year),FontName,TimesNewRoman,FontSize,12); ylabel(统计量,FontName,TimesNewRoman,Fontsize,12); %grid on hold on plot(x,0*ones(N,1),-.,linewidth,1); plot(x,1.96*ones(N,1),:,linewidth,1); plot(x,-1.96*ones(N,1),:,linewidth,1);
