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

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

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

%% MK突变检验
%% 修改日期 2022/7/29

function [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:n
    r1(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:n
   s(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:n
    time_series2(i)=time_series(n-i+1);
end
% 定义逆序统计量UB，长度=n，初始值=0
UB = zeros(size(time_series2));
for i= 1:n
    r2(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:n
   s2(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:n
   UB2(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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

注：代码并非全部原创，是修改自Mann-Kendall突变检测(mk突变检测)，优化了部分语法，修正了部分错误，更正了部分描述。

突变检验结果如下：

图像解读方法如下，来自mk突变点检测_（案例）时间序列分析之MK突变检测

Hope this may bring you some inspirations!