m-K检验 matlab代码

fisheryl

function [UFk,UBk2]= CalMK(y,x)
%Synax  [UFk,UBk2]= CalMK(y,x)
N=length(y);
n=length(y);
if nargin==1
    x=[1:1:n];
end

% 正序列计算---------------------------------
% 定义累计量序列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
if nargout==0
    figure
    ax1=gca;
    plot(x,UFk,'r-','linewidth',1.5); hold on
    plot(x,UBk2,'b-.','linewidth',1.5); hold on
    legend('UF统计量','UB统计量','0.05显著水平');
    xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
    ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
    hline([-1.96 0 1.96],'k--')
    set(ax1,'xlim',[x(1) x(end)])

    ax2 = axes('Position',get(ax1,'Position'),...
    'YAxisLocation','right','Color','none','XAxisLocation','top','xticklabel',[' ']);
   set(ax2,'xlim',get(ax1,'xlim'))

   line(x,y,'color',[0 0 0])


end

465582659

好人，收藏了下次试试

xiao豆豆

{:eb511:}{:eb511:}好东西

zhchlue

谢谢楼主的分享，有机会可以了解一下谢谢楼主的分享，有机会可以了解一下

757084134

谢谢了，不错的分享

胆小车迷

楼主的代码我用了，但是出来的图好像与魏凤英老师的例子不一样哦。
我看了半天也不知道到底是哪里出了问题。
http://bbs.06climate.com/forum.php?mod=viewthread&tid=34869，链接里的代码我检验过了，和书上是一样的。

taeyeon

谢谢了，不错的分享

何佳田力

楼主的代码我用了，但是出来的图好像与魏凤英老师的例子不一样哦