MK趋势检验和MK突变检验（代码分享及结果分析）

2023-09-10 08:38

短信预约 -IT技能 免费直播动态提醒

MK趋势检验

在时间序列趋势分析中，Mann-Kendall检验是世界气象组织推荐并已被广泛使用的非参数检验方法，最初由Mann和Kendall提出，现已被很多学者用来分析降雨、气温、径流和水质等要素时间序列的趋势变化。Mann-Kendall检验不需要样本遵从一定的分布，也不受少数异常值的干扰，适用于水文、气象等非正态分布的数据，计算简便。
代码如下：
这是代码1

% Mann-Kendall趋势检测 % Time Series Trend Detection Tests% [ z, sl, lcl, ucl ] = trend( y, dt )% where z = Mann-Kendall Statistic% sl = Sen's Slope Estimate% lcl = Lower Confidence Limit of sl% ucl = Upper Confidence Limit of sl% y = Time Series of Data% dt = Time Interval of Data% nor给定的显著性水平 当Z的绝对值大于等于1.64、 1.96、 2.58时则说明该时间序列分别通过了置信水平90%、95%、99%的显著性检验%y是待检测数据序列function [ z, sl, lcl, ucl ] = mk( y )n = length( y );dt=1;% Mann-Kendall Test for N > 40    % disp( 'Mann-Kendall Test;' );    % if n < 41,    % disp( 'WARNING - sould be more than 40 points' );    % end;% calculate statistics = 0;for k = 1:n-1    for j = k+1:n        s = s + sign( y(j) - y(k) );    endend% variance ( assuming no tied groups )v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;% test statisticif s == 0z = 0;elseif s > 0z = ( s - 1 ) / sqrt( v );elsez = ( s + 1 ) / sqrt( v );end% should calculate Normal value herenor = 1.96;%当Z的绝对值大于等于1.64、 1.96、 2.58时则说明该时间序列分别通过了置信水平90%、95%、99%的显著性检验% resultsdisp( [ ' n = ' num2str( n ) ] );disp( [ ' Mean Value = ' num2str( mean( y ) ) ] );disp( [ ' Z statistic = ' num2str( z ) ] );if abs( z ) < nor    disp( ' No significant trend' );    z = 0;elseif z > 0    disp( ' Upward trend detected' );else    disp( ' Downward trend detected' );end    disp( 'Sens Nonparametric Estimator:' );% calculate slopesndash = n * ( n - 1 ) / 2;s = zeros( ndash, 1 );i = 1;for k = 1:n-1    for j = k+1:n    s(i) = ( y(j) - y(k) ) / ( j - k ) / dt;    i = i + 1;    endend% the estimatesl = median( s );%M = median(A) 返回 A 的中位数值。disp( [ ' Slope Estimate = ' num2str( sl ) ] );% variance ( assuming no tied groups )v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;m1 = fix( ( ndash - nor * sqrt( v ) ) / 2 );m2 = fix( ( ndash + nor * sqrt( v ) ) / 2 );s = sort( s );lcl = s( m1 );ucl = s( m2 + 1 );disp( [ ' Lower Confidence Limit = ' ...num2str( lcl ) ] );disp( [ ' Upper Confidence Limit = ' ...num2str( ucl ) ] );

对于**标准值Z（Z statistic），其大于0，则序列呈上升趋势；若其小于0，则序列呈下降趋势。**当Z的绝对值大于等于1.64、 1.96、 2.58时则说明该时间序列分别通过了置信水平90%、95%、99%的显著性检验
斜率（Slope Estimate）用于衡量趋势大小，当斜率大于0 反应上升趋势，反之亦然。

还有代码2（计算结果一致，制图结果不一样）：

function  [Zs ,beta ,UFk ,UBk2 ]= MKtrend2(Data,n)% MKTest函数用于趋势和突变检验% 输入参数% Data    序列数据% n         序列长度% 输出参数% Zs       统计量% beta    斜率% UFk     统计量UFk% UBk2   逆序统计量%% 趋势分析线性:Mann-Kendall检验Sgn=zeros(n-1,n-1);        %初始化分配内存for i=1:n-1    for j=i+1:n        if((Data(j)-Data(i))>0)            Sgn(i,j)=1;        else            if((Data(j)-Data(i))==0)                Sgn(i,j)=0;            else                if((Data(j)-Data(i))<0)                    Sgn(i,j)=-1;                end            end        end    endendSmk=sum(sum(Sgn));VarS=n*(n-1)*(2*n+5)/18;if n>10    if Smk>0        Zs=(Smk-1)/sqrt(VarS);    else        if Smk==0            Zs=0;        else            if  Smk<0                Zs=(Smk+1)/sqrt(VarS);            end        end    endend%% beta 斜率 描述单调趋势t=1;for i=2:n   for j=1:(i-1)       temp(t)=( Data(i)-Data(j) )/( i-j );       t=t+1;   endendbeta=median( temp );%% 突变检验Sk=zeros(n,1);          % 定义累计量序列SkUFk=zeros(n,1);        % 定义统计量UFks = 0;% 正序列计算start---------------------------------for i=2:n   for j=1:i         if Data(i)>Data(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---------------------------------% 逆序列计算start---------------------------------Sk2=zeros(n);  % 定义逆序累计量序列Sk2UBk=zeros(n,1); s=0;Data2=flipud(Data);        % 按时间序列逆转样本yfor i=2:n   for j=1:i         if Data2(i)>Data2(j)           s=s+1;         else           s=s+0;         end   end   Sk2(i)=s;   E=i*(i-1)/4;   Var=i*(i-1)*(2*i+5)/72;   UBk(i,1)=0-(Sk2(i)-E)/sqrt(Var);end% 逆序列计算end------------------------------UBk2=flipud(UBk);UFk=UFk';UBk2=UBk2';%{figure(3)%画图plot(1:n,UFk,'r-','linewidth',1.5);hold onplot(1:n,UBk2,'b-.','linewidth',1.5);plot(1:n,1.96*ones(n,1),':','linewidth',1);% axis([1,n,-5,8]);legend('UF统计量','UB统计量','0.05显著水平');xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);%grid onhold onplot(1:n,0*ones(n,1),'-.','linewidth',1);plot(1:n,1.96*ones(n,1),':','linewidth',1);plot(1:n,-1.96*ones(n,1),':','linewidth',1);%}end

MK突变检验

MK突变趋势检验可以寻找到时间序列的突变发生点。
代码：

clcclear all% A=xlsread('C:\Users\DI cOS\Desktop\tt.xlsx','Sheet1','A2:B22');load('D:\00.研究生学习\08.EEMD\eemd_TOT_40_all.mat');B=allresidual(1732,:);A=B';% x=A(:,1); % 时间列aaa=1981:1:2020;x=aaa'; % 时间列% y=A(:,2); % 数据列y=A; % 数据列N=length(y);n=length(y);Sk=zeros(size(y));UFk=zeros(size(y));s=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;   Var=i*(i-1)*(2*i+5)/72;   UFk(i)=(Sk(i)-E)/sqrt(Var);endy2=zeros(size(y));Sk2=zeros(size(y));UBk=zeros(size(y));s=0;for i=1:n    y2(i)=y(n-i+1);endfor 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;   Var=i*(i-1)*(2*i+5)/72;  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);endUBk2=zeros(size(y));for i=1:n   UBk2(i)=UBk(n-i+1);endh1=plot(x,UFk,'r-','linewidth',1.5);hold onh2=plot(x,UBk2,'b-.','linewidth',1.5);axis([min(x),max(x),-5,6]);xlabel('年份','FontName','TimesNewRoman','FontSize',12);ylabel('时间序列数据','FontName','TimesNewRoman','Fontsize',12);hold onplot(x,0*ones(N,1),'-.','linewidth',1); ylim([-3 7]);h3=plot(x,2.58*ones(N,1),':','linewidth',1);plot(x,-2.58*ones(N,1),':','linewidth',1);legend([h1 h2 h3],'UF统计量','UB统计量','0.01显著水平','Location','NorthEast');f1=UFk;f2=UBk2;

如何解读结果，如何判断突变点？
得到的结果如下：
在这里插入图片描述
检验曲线图中若UF线在临界线内（两条显著性检验线之间，置信区间内）变动，表明变化曲线趋势和突变不明显;UF的值大于零，表明序列呈上升趋势，反之呈下降趋势;当其超过临界线时表明上升或下降趋势显著。如果UF和UB 2条曲线在临界线之间出现交点，则交点对应的时刻即为突变开始的时间;若交点出现在临界线外或出现多个交点，可结合其他检验方法进-步判定是否为突变点。