本帖最后由 vb1987 于 2013-6-12 23:27 编辑
最近写论文需要%MK检验法在网上收集了大量检验法matlab但是没有代码可以
%完全正确运行或分析信息不完整，结合多位网友编写MK经过我的改编，检验法顺利得到
%正确的操作结果，感谢所有网友，希望对有需要的盆友有所帮助
% Mann-Kendall突变检测
% 数据序列y
% 结果序列UFk，UBk2
%--------------------------------------------
%读取excel中的数据，赋给矩阵y
y的样本数为%
%A时间和径流数据列
A=xlswrite('数据.xls')
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(y,1);
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)的方差
% 由于对逆序列的累积量Sk在2的构造中，仍采用累加法，即后者大于前者时s加1，
% S的大小表示上升趋势的大小，在序列逆序后，应该表现出与原序列相反
% 因此，采用累加法统计趋势表现Sk2序列，统计公式(S(i)-E(i))/sqrt(Var(i))
% 也不应该改变，但统计量UBk应用相反的数字来表示正确的逆序序列趋势
UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end;
% ------------------------------逆序列计算end
% 上一步到了UBk逆序列在逆序时间上的趋势统计
% 与UFk在寻找突变点时，两条曲线应具有相同的时间轴，因此
% 逆转结果按时间顺序统计UBk，得到时间正序UBk2，做图用
UBk2=zeros(size(y));
% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
for i=1:n
UBk2(i)=UBk(n-i 1);
end;
% 使用突变检测图时UFk和UBk2
% 写入目标xls文件：f:\test2.xls
% 目标表单：Sheet1
% 目标区域：UFk从A1开始，UBk2从B1开始
xlswrite('f:\test2.xls',UFk,'Sheet1','A1');
xlswrite('f:\test2.xls',UBk2,'Sheet1','B1');
figure(3)%画图
plot(x,UFk,'r-','linewidth',1.5);
hold on
plot(x,UBk2,'b-.','linewidth',1.5);
plot(x,1.96*ones(N,一、linewidth',1);
axis([min(x),max(x),-5,5]);
legend('UF统计量，UB统计量，0.05显著水平)；
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel统计量，FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(x,0*ones(N,1),'-.','linewidth',1);
plot(x,1.96*ones(N,一、linewidth',1);
plot(x,-1.96*ones(N,一、linewidth',1);