MK突变分析UF统计量先超过显著线后又不超过显著线了，该如何分析？

jiangsaiping

jiangsaiping

代码是在家园里找的，如下。是不是没有问题啊？
%数据读取
%读取excel中的数据，赋给矩阵y
%获取y的样本数
%A为时间和数据列，x=时间列，y=检验数据列
A=xlsread('weifengyingtestdata.xls')%(魏凤英教材-上海年均气温1900-1990)
%A=xlsread('G:\LinshiData\weifengyingtestdata.xls')
x=A(:,1);%时间列
y=A(:,2);%数据列
N=length(y);
n=length(y);
%--------------------------------------------
% 正序列计算
Sk=zeros(size(y));% 定义累计量序列Sk，长度=y，初始值=0
UFk=zeros(size(y));% 定义统计量UFk，长度=y，初始值=0
s = 0;% 定义Sk序列元素s
% i从2开始。因i=1时，Sk(1)、E(1)、Var(1)均为0，此时UFk无意义，故令UFk(1)=0
for i=2:n
   for j=1:i-1
         if y(i)>=y(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;
%--------------------------------------------
% 逆序列计算
y2=zeros(size(y));% 构造逆序列y2，长度=y，初始值=0
Sk2=zeros(size(y));% 定义逆序累计量序列Sk2，长度=y，初始值=0
UBk=zeros(size(y));% 定义逆序统计量UBk，长度=y，初始值=0
s=0;% s归0
% 按时间序列逆转样本y。也可以使用y2=flipud(y);或者y2=flipdim(y,1);
for i=1:n
    y2(i)=y(n-i+1);
end;
% i从2开始。因i=1时，Sk2(1)、E(1)、Var(1)均为0，此时UBk无意义，故令UBk(1)=0
for i=2:n
   for j=1:i-1
         if y2(i)>=y2(j)
           s=s+1;
         else
           s=s+0;
         end;
   end;
   Sk2(i)=s;
   E=i*(i-1)/4; % Sk2(i)的均值
  Var=i*(i-1)*(2*i+5)/72; % Sk2(i)的方差
% 由于对逆序序列的累计量Sk2的构建中，依然用的是累加法，即后者大于前者时s加1，
% 则s的大小表征了一种上升的趋势的大小，而序列逆序以后，应当表现出与原序列相反
% 的趋势表现，因此，用累加法统计Sk2序列，统计量公式(S(i)-E(i))/sqrt(Var(i))
% 也不应改变，但统计量UBk应取相反数以表征正确的逆序序列的趋势，因此
  UBk(i)=0-(Sk2(i)-E)/sqrt(Var);
end;
%--------------------------------------------
% 此时上一步的到UBk表现的是逆序列在逆序时间上的趋势统计量
% 与UFk做图寻找突变点时，2条曲线应具有同样的时间轴，因此
% 再按时间序列逆转结果统计量UBk，得到时间正序的UBk2，
% 做突变检测图时，使用UFk和UBk2
UBk2=zeros(size(y));% 也可以使用UBk2=flipud(UBk);或者UBk2=flipdim(UBk,1);
for i=1:n
   UBk2(i)=UBk(n-i+1);
end;
%--------------------------------------------
%数据储存
%目标文件：G:\LinshiData\weifengyingtestout.xls；表单：Sheet1；区域：UFk从A1开始，UBk2从B1开始
xlswrite('weifengyingtestout.xls',UFk,'Sheet1','A1');
xlswrite('weifengyingtestout.xls',UBk2,'Sheet1','B1');
%--------------------------------------------
%绘图
figure(3)
plot(x,UFk,'r-','linewidth',1.5);
hold on
plot(x,UBk2,'b-.','linewidth',1.5);
plot(x,1.96*ones(N,1),':','linewidth',1);
T1=max(max(UFk));
T2=max(max(UBk));
T3=max(T1,T2);
axis([min(x),max(x),-ceil(T3)+3,ceil(T3)+1]);%根据实际情况调整
legend('UF统计量','UB统计量','0.05显著水平');
xlabel('t (year)','FontName','TimesNewRoman','FontSize',12);
ylabel('统计量','FontName','TimesNewRoman','Fontsize',12);
%grid on
hold on
plot(x,0*ones(N,1),'-.','linewidth',1);
plot(x,1.96*ones(N,1),':','linewidth',1);
plot(x,-1.96*ones(N,1),':','linewidth',1);

越饼干

我也想知道(●—●)