利用MATLAB进行MK检验
作者：第八星系-石头人
邮箱：2205455617@qq.com
读取数据
clc;clear;close all;
load data1.mat
year = data1(:,1);
data = data1(:,2);趋势检验
s = 0;
len=size(data,1);
for m=1:len-1
for n=m+1:len
if data(n) > data(m)
s = s+1;
elseif data(n) == data(m)
s = s+0;
elseif data(n) < data(m)
s = s-1;
end
end
end
vars = len*(len-1)*(2*len+5)/18;
% 计算方差
if s > 0
zc = (s-1)/sqrt(vars);
elseif s < 0
zc = (s+1)/sqrt(vars);
end
% zc>0趋势为增，zc<0趋势为减
% 当zc的绝对值大于等于1.64、1.96、2.58时则说明该时间序列分别通过了置信水平90%、95%、99%的显著性检验突变检验
% 计算UF
Sk=zeros(size(data)); UF=zeros(size(data)); s1=0;
for i=2:n
for j=1:i
if data(i)>data(j)
s1=s1+1;
end
end
Sk(i)=s1;
E=i*(i-1)/4;
Var=i*(i-1)*(2*i+5)/72;
UF(i)=(Sk(i)-E)/sqrt(Var);
end
% 计算UB
data2=zeros(size(data));
Sk2=zeros(size(data)); UBk=zeros(size(data)); s2=0;
for i=1:n
data2(i)=data(n-i+1);
end
for i=2:n
for j=1:i
if data2(i)>data2(j)
s2=s2+1;
end
end
Sk2(i)=s2;
E=i*(i-1)/4;
Var=i*(i-1)*(2*i+5)/72;
UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end
UB=zeros(size(data));
for i=1:n
UB(i)=UBk(n-i+1);
end画图
figure(i)
plot(year,UF,'bd-','LineWidth',1.2);
hold on
plot(year,UB,'ro-','LineWidth',1.2);
hold on
plot(xlim,,'g-.','LineWidth',0.8) % 画显著性线
hold on
plot(xlim,[-1.96,-1.96],'g-.','LineWidth',0.8)% 画显著性线
hold on
plot(xlim,,'k-','LineWidth',0.8)
hold on
axis()
xlabel('年份','FontSize',13,'FontWeight','bold','FontName','宋体');
ylabel('统计量','FontSize',13,'FontWeight','bold','FontName','宋体');
legend('UF','UB','0.05显著性水平','Location','best', ...
'FontSize',10,'FontWeight','bold','FontName','宋体');
% UF值>0,说明持续增长趋势，值在0.05显著性水平线上，说明通过0.05显著性检验
% UF和UB曲线的交点在置信水平区间[-1.96，1.96]内，并且确定交点具体年份，说明该年份呈现突变性增长状态；
% UF和UB曲线的交点不位于检验范围内，说明交点没有通过0.05 的检验，所以该年份参数突变性上升不具有突变性
完整代码
clc;clear;close all;
load data1.mat
year = data1(:,1);
data = data1(:,2);
s = 0;
len=size(data,1);
for m=1:len-1
for n=m+1:len
if data(n) > data(m)
s = s+1;
elseif data(n) == data(m)
s = s+0;
elseif data(n) < data(m)
s = s-1;
end
end
end
vars = len*(len-1)*(2*len+5)/18;
if s > 0
zc = (s-1)/sqrt(vars);
elseif s < 0
zc = (s+1)/sqrt(vars);
end
Sk=zeros(size(data)); UF=zeros(size(data)); s1=0;
for i=2:n
for j=1:i
if data(i)>data(j)
s1=s1+1;
end
end
Sk(i)=s1;
E=i*(i-1)/4;
Var=i*(i-1)*(2*i+5)/72;
UF(i)=(Sk(i)-E)/sqrt(Var);
end
data2=zeros(size(data));
Sk2=zeros(size(data)); UBk=zeros(size(data)); s2=0;
for i=1:n
data2(i)=data(n-i+1);
end
for i=2:n
for j=1:i
if data2(i)>data2(j)
s2=s2+1;
end
end
Sk2(i)=s2;
E=i*(i-1)/4;
Var=i*(i-1)*(2*i+5)/72;
UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end
UB=zeros(size(data));
for i=1:n
UB(i)=UBk(n-i+1);
end
year=1991:2022;
figure(i)
plot(year,UF,'bd-','LineWidth',1.2);
hold on
plot(year,UB,'ro-','LineWidth',1.2);
hold on
plot(xlim,,'g-.','LineWidth',0.8)
hold on
plot(xlim,[-1.96,-1.96],'g-.','LineWidth',0.8)
hold on
plot(xlim,,'k-','LineWidth',0.8)
hold on
axis()
xlabel('年份','FontSize',13,'FontWeight','bold','FontName','宋体');
ylabel('统计量','FontSize',13,'FontWeight','bold','FontName','宋体');
legend('UF','UB','0.05显著性水平','Location','best', ...
'FontSize',10,'FontWeight','bold','FontName','宋体')
微信搜索“第八星系人造大气理论爱好者”公众号，关注获取文章数据
页:
[1]