%?Mann-Kendall突變檢測?
%?數據序列y
%?結果序列UFk，UBk2
%--------------------------------------------
%讀取excel中的數據，賦給矩陣y
%獲取y的樣本數
%A為時間和徑流數據列
A=xlsread（‘kk.xls‘‘Sheet1‘）
x=A（:1）;%時間序列
y=A（:2）;%徑流數據列
N=length（y）;
n=length（y）;
%?正序列計算---------------------------------
%?定義累計量序列Sk，長度=y，初始值=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）/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（y1）;
for?i=1:n
????y2（i）=y（n-i+1）;
end;
%?i從2開始，因為根據統計量UBk公式，i=1時，Sk2（1）、E（1）、Var（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（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（UBk1）;
for?i=1:n
???UBk2（i）=UBk（n-i+1）;
end;
%?做突變檢測圖時，使用UFk和UBk2
%?寫入目標xls文件：f:\test2.xls
%?目標表單：Sheet1
%?目標區域：UFk從A1開始，UBk2從B1開始
xlswrite（‘E:\matlab\manasMK\manasoutput.xls‘UFk‘Sheet1‘‘A1‘）;
xlswrite（‘E:\matlab\manasMK\manasoutput.xls‘UBk2‘Sheet2‘‘B1‘）;
figure（3）%畫圖
plot（xUFk‘r-‘‘linewidth‘1.5）;
hold?on
plot（xUBk2‘b-.‘‘linewidth‘1.5）;
plot（x1.96*ones（N1）‘:‘‘linewidth‘1）;
axis（[min（x）max（x）-55]）;
legend（‘UF統計量‘‘UB統計量‘‘0.05顯著水平‘）;
xlabel（‘t?（year）‘‘FontName‘‘TimesNewRoman‘‘FontSize‘12）;
ylabel（‘統計量‘‘FontName‘‘TimesNewRoman‘‘Fontsize‘12）;
%grid?on
hold?on
plot（x0*ones（N1）‘-.‘‘linewidth‘1）;
plot（x1.96*ones（N1）‘:‘‘linewidth‘1）;
plot（x-1.96*ones（N1）‘:‘‘linewidth‘1）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件??????19968??2015-03-06?10:24??MKtest\kk.xls

?????文件???????2960??2015-03-06?11:39??MKtest\mk.m

?????目錄??????????0??2015-03-06?11:10??MKtest

-----------?---------??----------?-----??----

????????????????22928????????????????????3