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