matlab程序突变点 检验

fangcha1(k,1)=k*(k-1)*(2*k+5)/72;
zonghe1(1,1)=0;
zonghe1(k,1)=count1;
end
%-------------------COMPUTATION IN ADVERSE SEQUENCY-----------------
% 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.
% 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.
ylabel('统计量','FontSize',12)
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 %%
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;
% function [z]=rain (x)
% x为一时间序列(列向量)
x=[42.4
40.3
42.8
42.6
46.6
43.9
44.0
46.0
46.1
46.4
42.0
43.8
46.1
42.1
45.9
43.9
45.1
46.8
44.0
47.3
44.6
45.4
[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')
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);
40.3
42.9
39.8
44.8
41.3
42.6
];
% define the time
time=1955+(0:length(x)-1); % please alter the starting year of this series.
n=length(x); %the size of the series
%-------------------- 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);
zonghe2(1,1)=0;
zonghe2(m,1)=count2;
end
%------------------------autocorrelation analysis------------------------
kkk = fix(n/3);
for ii = 0:kkk;
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)
% The script was written by Dr. Qiang Zhang from
% Nanjing Institute of Geography and Limnology,
% Chinese Academy of Sciences on March 26, 2006.
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;
end
% -----------------COMPUTATION in normal sequence--------------------
for k=2:n
count1=0;
for j=1:k-1
if x(k)>x(j)
count1=count1+1;
adjust=1+2*(rr(2)^(n+1)-n*rr(2)^2+(n-1)*rr(2))/(n*(rr(2)-1)^2);
%-----------调整系数计算完毕-----------
%对原MK值根据自相关系数计算情况进行调值---pper(2)
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--------
%% 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);
%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;
46.9
47.0
43.9
44.1
45.4
46.4
45.4
43.0
42.3
43.1
47.3
47.3
46.9
46.1
48.5
49.6
49.5
46.0
49.3
47.8
43.3
48.5
46.0
45.3
44.9
48.0
49.5
45.5
44.0
45.6
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);
%----------计算调整系数----------------
plot(time,ones(size(time))*1.96, ':k',time,ones(size(time))*-1.96, ':k')
plot(time,zeros(size(time)),'-k')