mk检验matlab程序

% Mann-Kendall突变检测 % 数据序列y % 结果序列UFk，UBk2 %-------------------------------------------- %读取excel中的数据，赋给矩阵y %获取y的样本数 %A为时间和径流数据列 A=xlsread('kk.xls','Sheet1') x=A(:,1);%时间序列 y=A(:,2);%径流数据列

% 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