- 积分
- 1115
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-10-2
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
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年。
|
-
魏凤英老师书上的编程
-
MK检验结果
-
按照魏奉英老师的书写的编程算出来的结果
-
改完循环之后的结果
|