MK检验法

dwhuang

A=xlsread('test-mk.xlsx');

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)的方差

% 由于对逆序序列的累计量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(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,1),':','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,1),':','linewidth',1);

plot(x,-1.96*ones(N,1),':','linewidth',1);

dwhuang

function [ UF ] = MannKendall( Runoff )
%UNTITLED 此处显示有关此函数的摘要
%输入数据格式，第一列为runoff
%a=size(Data);
%b=max(a);
%Year=Data(:,1);
%Runoff=Data(:,2);
%构建r序列
a=size(Runoff);
b=max(a);
r=[];
for i=1:b
    xi=Runoff(i);
    xj=Runoff(1:i);
    R=xj(xj<xi);
    m=size(R);
    if min(m)==0
        r_v=0;
    else
        r_v=max(m);
    end
    r=[r;r_v];
end
%构建s序列
s=zeros(b,1);
for j=1:b
    s(j)=sum(r(1:j));
end
%构建UF序列
UF=[];
for k=2:b
    E_sk=k*(k-1)/4;
    Var_sk=k*(k-1)*(2*k+5)/72;
    UF_k=(s(k)-E_sk)./(Var_sk.^0.5);
    UF=[UF;UF_k];
end
UF=[0;UF];
end

manrabbit

赞一个，最近也学到了这个，参考一下您的代码

缪大肉喵喵喵

赞赞赞，之前看了好几个fortran和matlab脚本都没有注释，理解起来非常费劲，看了这个感觉好多了！非常感谢！

缪大肉喵喵喵

我现在有一个问题，我按照这个程序改编成了ncl 脚本，但是画出来的两根折线是完全对称的，这应该不对吧？这是为什么呢？?

syzeos

缪大肉喵喵喵 发表于 2018-12-20 22:23
我现在有一个问题，我按照这个程序改编成了ncl 脚本，但是画出来的两根折线是完全对称的，这应该不对吧？这 ...

您好，问题解决了吗？我也遇到了这个问题！让人头大

帅蛋fly

syzeos 发表于 2019-5-29 13:20
您好，问题解决了吗？我也遇到了这个问题！让人头大

您好，请问可以问一下您的问题解决了嘛？

1271756664

楼主您好，请问您这个是MATLAB脚本吗？初学，请赐教！

风倾云歌

楼主你好，请问为什么我的图最后出来3条DATA123，明明是在hold on后面，怎么会这样呢？还有就是你发的第二大段是干什么用的呢，谢谢

潘W七七七

这个太好了，解决了我的问题，谢谢了！