%?Mann-Kendall突變檢測(cè)?
%?數(shù)據(jù)序列y
%?結(jié)果序列UFk，UBk2
%--------------------------------------------
%讀取excel中的數(shù)據(jù)，賦給矩陣y
%獲取y的樣本數(shù)
%A為時(shí)間和徑流數(shù)據(jù)列
A=xlsread（‘kk.xls‘‘Sheet1‘）
x=A（:1）;%時(shí)間序列
y=A（:2）;%徑流數(shù)據(jù)列
N=length（y）;
n=length（y）;
%?正序列計(jì)算---------------------------------
%?定義累計(jì)量序列Sk，長(zhǎng)度=y，初始值=0
Sk=zeros（size（y））;
%?定義統(tǒng)計(jì)量UFk，長(zhǎng)度=y，初始值=0
UFk=zeros（size（y））;
%?定義Sk序列元素s
s?=?0;
%?i從2開始，因?yàn)楦鶕?jù)統(tǒng)計(jì)量UFk公式，i=1時(shí)，Sk（1）、E（1）、Var（1）均為0
%?此時(shí)UFk無(wú)意義，因此公式中，令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;
%?------------------------------正序列計(jì)算end
%?逆序列計(jì)算---------------------------------
%?構(gòu)造逆序列y2，長(zhǎng)度=y，初始值=0
y2=zeros（size（y））;
%?定義逆序累計(jì)量序列Sk2，長(zhǎng)度=y，初始值=0
Sk2=zeros（size（y））;
%?定義逆序統(tǒng)計(jì)量UBk，長(zhǎng)度=y，初始值=0
UBk=zeros（size（y））;
%?s歸0
s=0;
%?按時(shí)間序列逆轉(zhuǎn)樣本y
%?也可以使用y2=flipud（y）;或者y2=flipdim（y1）;
for?i=1:n
????y2（i）=y（n-i+1）;
end;
%?i從2開始，因?yàn)楦鶕?jù)統(tǒng)計(jì)量UBk公式，i=1時(shí)，Sk2（1）、E（1）、Var（1）均為0
%?此時(shí)UBk無(wú)意義，因此公式中，令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）的方差
%?由于對(duì)逆序序列的累計(jì)量Sk2的構(gòu)建中，依然用的是累加法，即后者大于前者時(shí)s加1，
%?則s的大小表征了一種上升的趨勢(shì)的大小，而序列逆序以后，應(yīng)當(dāng)表現(xiàn)出與原序列相反
%?的趨勢(shì)表現(xiàn)，因此，用累加法統(tǒng)計(jì)Sk2序列，統(tǒng)計(jì)量公式（S（i）-E（i））/sqrt（Var（i））
%?也不應(yīng)改變，但統(tǒng)計(jì)量UBk應(yīng)取相反數(shù)以表征正確的逆序序列的趨勢(shì)
??UBk（i）=0-（Sk2（i）-E）/sqrt（Var）;
end;
%?------------------------------逆序列計(jì)算end
%?此時(shí)上一步的到UBk表現(xiàn)的是逆序列在逆序時(shí)間上的趨勢(shì)統(tǒng)計(jì)量
%?與UFk做圖尋找突變點(diǎn)時(shí)，2條曲線應(yīng)具有同樣的時(shí)間軸，因此
%?再按時(shí)間序列逆轉(zhuǎn)結(jié)果統(tǒng)計(jì)量UBk，得到時(shí)間正序的UBk2，做圖用
UBk2=zeros（size（y））;
%?也可以使用UBk2=flipud（UBk）;或者UBk2=flipdim（UBk1）;
for?i=1:n
???UBk2（i）=UBk（n-i+1）;
end;
%?做突變檢測(cè)圖時(shí)，使用UFk和UBk2
%?寫入目標(biāo)xls文件：f:\test2.xls
%?目標(biāo)表單：Sheet1
%?目標(biāo)區(qū)域：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統(tǒng)計(jì)量‘‘UB統(tǒng)計(jì)量‘‘0.05顯著水平‘）;
xlabel（‘t?（year）‘‘FontName‘‘TimesNewRoman‘‘FontSize‘12）;
ylabel（‘統(tǒng)計(jì)量‘‘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）;

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件??????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