Matlab完成MK检验

MK突变检测需要，查阅资料，发现口径不一，传说中的DPS软件不能用（正版的要钱）。大多数代码只能自己看图说话，找UB、UF交点作为突变点，但好多没有格网，放大都看不了，还有，作图上没有加For循环，一次画一个图，倒不是缺点，只是处理的数据多，懒，所以，根据网上说明自己修改了代码，感谢相关博主（见附录）。另，这个检测还可以用一个神奇额R包：Trend，进行相应检验。这里不详细论述。若有其他好方法请指点。

代码：

clear;clc;

Data=xlsread('AnnualIntensity_16.csv');

Year=Data(:,1);

[a,b]=size(Year);

Sig1=1.96*ones(a,b);

Sig2=-1.96*ones(a,b);

for i =1:17

Runoff=Data(:,i);

UF=MannKendall(Runoff);

Runoff2=flipud(Runoff);

UF2=MannKendall(Runoff2);

UB=-flipud(UF2);

figure(i)

plot(Year,UF,'k');

hold on

plot(Year,UB,'r');

plot(Year,Sig1,'--');

plot(Year,Sig2,'--');

xlabel('Year');

ylabel('统计量');

legend('UF','UB','α<0.05');

hold off

title={'Year','UF','UB','α','α'};

shuju=[Year UF UB Sig1 Sig2];

hold on

%hold on

grid(gca,'minor')% 画出次格网

Runoff=[] %清空数组

end;

%xlswrite('test01.xlsx',title,'Sheet1','A1');

%xlswrite('test01.xlsx',shuju,'Sheet1','A2');

            ------------------------------分-----------------------割----------------------线---------------------------------

增加保存图像功能：

clear;clc

Data=xlsread('numofExtreme_16.csv');

StationALL=xlsread('StationID.xlsx');

Year=Data(:,1);

[a,b]=size(Year);

Sig1=1.96*ones(a,b);

Sig2=-1.96*ones(a,b);

for i =2:17

Runoff=Data(:,i);

StationID=StationALL(i-1);

UF=MannKendall(Runoff);

Runoff2=flipud(Runoff);

UF2=MannKendall(Runoff2);

UB=-flipud(UF2);

figure(i)

plot(Year,UF,'k');

hold on

plot(Year,UB,'r');

plot(Year,Sig1,'--');

plot(Year,Sig2,'--');

xlabel('Year');

ylabel('统计量');

legend('UF','UB','α<0.05');

hold off

title={'Year','UF','UB','α','α'};

%shuju=[Year UF UB Sig1 Sig2];



hold on;%hold on

grid(gca,'minor');% 画出次格网

Runoff=[]; %清空数组

savepath=['D:\Figures\EPF',num2str(StationID),'.png'];%保存路径

saveas(i,savepath);%保存为png文件

end;

%xlswrite('test01.xlsx',title,'Sheet1','A1');

%xlswrite('test01.xlsx',shuju,'Sheet1','A2');

---------------------------------------------------------------------------------------------------------------------------------------------------------------------

应邀加个R的trend包趋势分析代码，需要先安装trend包。为防止乱码，在这里写一下，仅供参考

#EPI
rm(list=ls());
library(trend);
data <- ts(data=c(64.9434375,68.8955625,75.9810625,59.627625,53.6349375,58.2556875,53.319375,53.4730625,59.842,61.718125,56.6345,55.6553125,52.0409375,58.372875,56.81325,60.070625,58.4079375,52.10657143,54.649125,53.1128125,63.21625,57.8413125,55.9164375,56.5718,57.909375,54.8268125,55.968625,55.69225,60.05925,61.7584375,63.6616875,59.565375,57.103625,61.4910625,60.4355,58.321875,63.4930625,54.5916,63.313,56.3408,61.45193333,56.2614,56.25506667,58.99106667,71.18806667,54.0644,72.23966667,58.05253333,66.1172,54.56493333,62.973,67.73566667),start=1961);
res <- pettitt.test(data);res

参考：

Mann-Kendall突变检测(mk突变检测) - CSDN博客  http://blog.csdn.net/liyanzhong/article/details/41867859

科学网—matlab版本的MK程序，包含突变和趋势检验 - 张乐乐的博文  http://blog.sciencenet.cn/blog-1103122-851049.html

Mann-Kendall突变检测(mk突变检测) – MATLAB中文论坛  http://www.ilovematlab.cn/thread-246993-1-1.html