- 积分
- 12226
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-12-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
请问大家,所用Mann-kendall突变检验的程序代码的数据该如何加载和修改?本人是新手,有一个57年的时间序列,改如何使用这个代码进行突变检验?需要修改哪些参数?
function Aburptchangetime=mk(array,time,alpha,drawing) %突变检验函数
%array为待检验的时间序列
%time为与array对应的时间,默认按时间升序排列
%array与time均为行向量
%alpha为显著水平
%drawing为判定条件,如果drawing=1,则绘图
%Abruptchangetime为返回的突变时间
%正序计算
n=length(array);%返回array元素个数
s=0;%初始化s
for i=2:n
for j=1:i
if array(i)>array(j)
s=s+1;
end
end
Sk(i)=s;
E=i*(i+1)/4;
Var=i*(i-1)*(2*i+5)/72;
UFk(i)=(Sk(i)-E)/Var^0.5;
end
%逆序计算
clear Sk;
s=0;%再次初始化s
for i=1:n
arrayinverting(i)=array(n-i+1);%反转array为arrayinverting
end
for i=2:n
for j=1:i
if arrayinverting(i)>arrayinverting(j)
s=s+1;
end
end
Sk(i)=s;
E=i*(i+1)/4;
Var=i*(i-1)*(2*i+5)/72;
UB(i)=-(Sk(i)-E)/Var^0.5;
end
for i=1:n
UBk(i)=UB(n-i+1);%反转UB得到UBk
end
U=norminv(1-alpha/2,0,1);%由显著水平alpha计算临界值
%计算UFk与UBk曲线的交点,并判断是否为突变点
m=0;%突变点数初始化为0
for i=2:n
x=time(i-1:i);
y1=UFk(i-1:i);
y2=UBk(i-1:i);
xishu1=polyfit(x,y1,1);
xishu2=polyfit(x,y2,1);
xishu=xishu1-xishu2;
changetime=roots(xishu);
if changetime>time(i-1)
if changetime<time(i)
if abs(xishu1(1)*changetime+xishu1(2))<U
m=m+1;
Aburptchangetime(m)=changetime;
end
end
end
end
if m==0
Aburptchangetime(1)=0;
end
%绘图
if drawing==1
figure(1);%准备绘图
plot(time,UFk,'r-','linewidth',2.0);
hold on
plot(time,UBk,'g-.','linewidth',2.0);
plot(time,U*ones(n,1),':','linewidth',1);
plot(time,-U*ones(n,1),':','linewidth',1);
axis([min(time),max(time),-10,10]);
legend('UF统计量','UB统计量',[num2str(alpha),'显著水平']);
xlabel('time','Fontsize',12);
ylabel('U','Fontsize',12);
end
end
|
|