MATLAB Mann-Kendall突变检验 (mk突变检验)

任务描述：对时间序列进行MK突变检验：

将MK突变检验的代码封装为函数，直接调用即可，代码如下：

%% MK突变检验
%% 修改日期 2022/7/29function [UF,UB] = MKbreak(time_series)n = length(time_series);%%  ---------------------------------正序列计算--------------------------------
% 定义统计量UF，长度=n，初始值=0；
UF=zeros(size(time_series));E = n*(n-1)/4;                         % E是累计数s(k)的均值
Var = n*(n-1)*(2*n+5)/72;              % Var是累计数s(k)的方差% 定义秩序列，r(i)记录的是第i个时刻，其数值大于j时刻(其中j<=i)数值的个数的
for i= 1:nr1(i) = sum(time_series(i)>time_series(1:i));
end% 定义累计量序列s，s(k)记录的是第i个时刻(其中i<=k)，其数值大于j时刻(其中j<=i)数值个数的累计数
s = zeros(size(time_series));
% k从2开始，因为根据统计量UF(k)公式，k=1时，s(1)、E(1)、Var(1)均为0，此时UF无意义，因此公式中，令UFk(1)=0
for k = 2:ns(k) = sum(r1(1:k));E = k*(k-1)/4;                     % s(k)的均值Var = k*(k-1)*(2*k+5)/72;          % s(k)的方差UF(k) = (s(k)-E)/sqrt(Var);
end%% ---------------------------------逆序列计算--------------------------------
% 按时间序列逆转样本,构造逆序列time_series2和
for i=1:ntime_series2(i)=time_series(n-i+1);
end
% 定义逆序统计量UB，长度=n，初始值=0
UB = zeros(size(time_series2));
for i= 1:nr2(i) = sum(time_series2(i)>time_series2(1:i));
end% 定义逆序累计量序列s2，长度=n，初始值=0
s2 = zeros(size(time_series2));
% i从2开始，因为根据统计量UB(k)公式，i=1时，s2(1)、E(1)、Var(1)均为0，此时UB无意义，因此公式中，令UB(1)=0
for k = 2:ns2(k) = sum(r2(1:k));E = k*(k-1)/4;                     % s2(k)的均值Var = k*(k-1)*(2*k+5)/72;          % s2(k)的方差
% 由于对逆序序列的累计量S2的构建中，依然用的是累加法，即后者大于前者时r加1，
% 则r的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
% 的趋势表现，因此，用累加法统计S2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势UB(k) = -(s2(k)-E)/sqrt(Var);
end%% ---------------------------------绘图
% 此时上一步的到UB表现的是逆序列在逆序时间上的趋势统计量
% 与UF做图寻找突变点时，2条曲线应具有同样的时间轴，因此
% 再按时间序列逆转结果统计量UB，得到时间正序的UB，做图用
for i=1:nUB2(i)=UB(n-i+1);
end
x = 1:n;plot(x,UF,'g-','linewidth',1.5);
hold on
plot(x,UB2,'m-','linewidth',1.5);
% plot(x,UB,'m-','linewidth',1.5);
plot(x,1.96*ones(n,1),'-.r','linewidth',1);
plot(x,0*ones(n,1),'-','color',[0.2,0.2,0.2],'linewidth',1);
plot(x,-1.96*ones(n,1),'-.r','linewidth',1);
grid(gca,'minor')% 画出次格网axis([min(x),max(x),-5,5]);
legend('UF统计量','UB统计量','0.05显著水平');
set(gca,'Fontsize',12)
set(gca,'ytick',-5:1:5)
xlabel('{\itt} (year)','FontSize',15);
ylabel('统计量','Fontsize',15);end