发一个Pettitt的程序

蓬蒿人

load('数据');
y=数据;
year=year;
N=length(y);
n=length(y);
% 正序列计算,定义累计量序列Sk，长度=y，初始值=0
Sk=zeros(n,1);
for i=2:n
   for j=1:n
         if y(i)>y(j) %能否取到等于号看到了好几种说法
           Sk(i)=Sk(i)+1;
         elseif y(i)==y(j)
               Sk(i)=Sk(i);
         else
           Sk(i)=Sk(i)-1;
         end;
   end;
end;
%输出Max的最大值及最大值所在的位置。
m=max(abs(Sk(2:n)));
p=find(abs(Sk)==m);
P_t=2*exp(-6*(m^2)*(n^3+n^2));
year_need=zeros(length(p),1);
% 获取一个和p一样的长度数列
if P_t<=0.5
    year_need=year(p);
else
    print('该序列趋势不显著');
end
%在图片上绘制出最大的m；
cankaoxian=zeros(n,1);
cankaoxian(:)=m;
plot(year,Sk,'k-','linewidth',2);
hold on;
plot(year,cankaoxian,'b--','linewidth',2);
plot(year,-cankaoxian,'b--','linewidth',2);
plot([year_need(1),year_need(1)],[-cankaoxian,cankaoxian],'r--','linewidth',1);
plot([year_need(2),year_need(2)],[-cankaoxian,cankaoxian],'r--','linewidth',1);
set(gca,'xlim',[1982 2015]);


跟一些程序相比，这个帖子强调的点就这一行代码：   for j=1:n
魏奉英老师的书上的内容如图所示。
但是按这张图写的程序算出来的突变点和MK检验的结果相差甚远，突变点显示在08年和13年。
MK检验代码移步：http://bbs.06climate.com/forum.p ... ight=MK%BC%EC%D1%E9，我也是参考这个画的。


这个用MK检验画出来的（因为是师兄论文要用到的图，所以我抹去了一些细节，就剩下交点，可以看到大概是在97年附近）：
改成for j=1:n以后就顺眼多了，突变点在97年。

蓬蒿人

我天，语文学得不好的我，发完才发现都是病句。。。反正意思就这个意思了。。。

puck66

Pettitt 。。。。。是啥来着？

蓬蒿人

puck66 发表于 2023-10-22 11:47
Pettitt 。。。。。是啥来着？

突变点检验呢。