MK检验程序_mk检验程序
MK检验程序由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“mk检验程序”。
%The script was written by Dr.Qiang Zhang from
%Nanjing Institute of Geography and Limnology,%Chinese Academy of Sciences on March 26, 2006.%Financially supported by Alexander von Humboldt Stiftung, Germany
% if you use this script for research, please cite the references below:
%%Reference:
%Gerstengarbe, F.W.and Werner, P.C., 1999.Estimation of the beginning and end of%%recurrent events within a climate regime.Climate Research 11: 97-107.%Qiang Zhang, Chunling Liu, Chong-yu Xu, Youpeng Xu, Tong Jiang.Observed trends of annual maximum
%water level and streamflow during past 130 years in the Yangtze River basin, China.Journal of Hydrology, 2006, 324, 255-265.function [z]=rain(x)
% x为一时间序列(列向量)
x=[7.1
7.3
6.2
6.8
6.3
7.1
6.8
5.6
6.3
5.8
6.1
6.2
6.8
7.0
6.3
7.5
6.0
6.8
7.0
6.8
6.4
6.3
7.1
7.0
5.9
6.2
6.2
7.7
6.6
7.3
7.4
7.0
6.9
6.6
7.6
6.9
6.8
7.8
8.8
8.8
7.4
7.9
8.1
7.2
7.8
7.2
8.3
8.4
7.4
8.0
7.7
];
% define the time
time=1960+(0:length(x)-1);% please alter the starting year of this series.n=length(x);%the size of the series
[rows,cols]=size(x);
%---------JUDGEMENT OF TIME SERIES------------
if rows==1
error('x should be a column vector')
end
if cols~=1
error('x should be a column vector')
end
%-----------------COMPUTATION in normal sequence--------------------
for k=2:n
count1=0;
for j=1:k-1
if x(k)>x(j)
count1=count1+1;
end
%junzhi=k*(k-1)/4;
% fangcha=k*(k-1)*(2*k+5)/72;
%z(1,1)=0;
% z(k,1)=(count1-junzhi)/sqrt(fangcha);
end
junzhi1(k,1)=k*(k-1)/4;
fangcha1(k,1)=k*(k-1)*(2*k+5)/72;
zonghe1(1,1)=0;
zonghe1(k,1)=count1;
end
%-------------------COMPUTATION IN ADVERSE SEQUENCY-----------------x1=flipud(x);
for m=2:n
count2=0;
for i=1:m-1
%bijiao=x1(m)-x1(i);
if x1(m)>x1(i)
count2=count2+1;
end
%junzhi=m*(m-1)/4;
%fangcha=m*(m-1)*(2*m+5)/72;
%z(1,2)=0;
%z(m,2)=(count-junzhi)/sqrt(fangcha);
end
junzhi2(m,1)=m*(m-1)/4;
fangcha2(m,1)=m*(m-1)*(2*m+5)/72;
zonghe2(1,1)=0;
zonghe2(m,1)=count2;
end
%------------------------autocorrelation analysis------------------------
kkk = fix(n/3);
for ii = 0:kkk;
rr(ii+1)= sum(((x(1:n-ii,1)-mean(x))/std(x)).*((x(1+ii:n,1)-mean(x))/std(x)))/(n-ii);end
%---------------confidence level for autocorrelation coefficient--------
for jj = 0:kkk
pk(jj+1)= length(x(1:n-jj));
end
lower =(-1-1.96*(pk-2).^0.5)./(pk-1);
upper =(-1+1.96*(pk-2).^0.5)./(pk-1);
%----------计算调整系数----------------
adjust=1+2*(rr(2)^(n+1)-n*rr(2)^2+(n-1)*rr(2))/(n*(rr(2)-1)^2);
%-----------调整系数计算完毕-----------
%对原MK值根据自相关系数计算情况进行调值-------------
if rr(2)>upper(2)
%--------------------Z2 and Z1 value computation----------------
junzhi1(1,1)=0;
fangcha1(1,1)=0;
z(1,1)=0;
zonghe1=cumsum(zonghe1);
z(2:n,1)=(zonghe1(2:end)-junzhi1(2:end))./sqrt(fangcha1(2:end)*adjust);%%Z2 value%%
z(1,2)=0;
junzhi2(1,1)=0;
fangcha2(1,1)=0;
zonghe2=cumsum(zonghe2);
z(2:n,2)=(zonghe2(2:end)-junzhi2(2:end))./sqrt(fangcha2(2:end)*adjust);else
junzhi1(1,1)=0;
fangcha1(1,1)=0;
z(1,1)=0;
zonghe1=cumsum(zonghe1);
z(2:n,1)=(zonghe1(2:end)-junzhi1(2:end))./sqrt(fangcha1(2:end));%%Z2 value%%
z(1,2)=0;
junzhi2(1,1)=0;
fangcha2(1,1)=0;
zonghe2=cumsum(zonghe2);
z(2:n,2)=(zonghe2(2:end)-junzhi2(2:end))./sqrt(fangcha2(2:end));end
z(:,2)=-z(:,2);
xx=flipud(z(:,2));
z(:,2)=xx;
figure(1)
plot(time,z(:,1),'r-','LineWidth',1.5)
plot(time,z(:,2),'b','LineWidth',1.5)
box off
hold on
plot(time,z(:,1),'r-','LineWidth',1.5)
plot(time,ones(size(time))*1.96, ':k',time,ones(size(time))*-1.96, ':k')plot(time,zeros(size(time)),'-k')
hold on
xlabel('Time(year)','FontSize',12)
ylabel('统计量','FontSize',12)
