trendMMK_MK_MMK趋势检验_MK趋势分析_

% % Time Series Trend Detection Tests % % [ z, sl, lcl, ucl ] = trend( y, dt ) % % where z = Modified Mann-Kendall Statistic % sl = Sen's Slope Estimate % lcl = Lower Confidence Limit of sl % ucl = Upper Confidence Limit of sl % y = Time Series of Data % dt = Time Interval of Data % % Bob Newell, February 1996 % %-------------------------------------------------- % function [ z, sl, lcl, ucl ] = trendMMK( y, dt ) n = length( y ); %-------------------------------------------------- % Modified Mann-Kendall Test for N > 40 % disp( 'Modified Mann-Kendall Test:' ); if n < 41, disp( 'WARNING - sould be more than 40 points' ); end; % calculate statistic s = []; v = 1; x = zeros(n,1); for i = 1:n-1, for j = i+1:n, s(v) = ( y(j) - y(i) )/( j - i ); end; end; beta = median( s ); for t = 1:n x(t) = y(t) - beta * t; end R = FractionalRankings(x); Rmean = sum(R)/length(R); r=[]; for z=1:(n-1); yy=[]; zz=[]; for q=1:(n-z); yy(q)=(R(q)-Rmean)*(R(q+z)-Rmean); end for p=1:n; zz(p)=(R(p)-Rmean)*(R(p)-Rmean); end r(z)=sum(yy)/sum(zz); end % variance ( with tied groups ) xo=[]; for z=1:(n-1); xo(z)=(n-z)*(n-z-1)*(n-z-2)*r(z); end eta=1+(2/(n*(n-1)*(n-2)))*sum(xo); sm=tabulate(y); tindice=find(sm(:,2)>1); tnum=length(tindice); ox=[]; for xx=1:tnum; ox(xx)=sm(tindice(xx),2)*(sm(tindice(xx),2)-1)*(2*sm(tindice(xx),2)+5)/18; end sox=sum(ox); va=n*(n-1)*(2*n-5)/18-sox; va1=va*eta; ss = 0; for k = 1:n-1, for j = k+1:n, ss = ss + sign( y(j) - y(k) ); end; end; % test statistic if ss == 0, z = 0; elseif ss > 0, z = ( ss - 1 ) / sqrt( va1 ); else, z = ( ss + 1 ) / sqrt( va1 ); end; % should calculate Normal value here nor = 1.96; % results disp( [ ' n = ' num2str( n ) ] ); disp( [ ' Mean Value = ' num2str( mean( y ) ) ] ); disp( [ ' Z statistic = ' num2str( z ) ] ); if abs( z ) < nor, disp( ' No significant trend' ); z = 0; elseif z > 0, disp( ' Upward trend detected' ); else, disp( ' Downward trend detected' ); end; %---------------------------------------------------- % 这一下为sen斜率估计！！！！！！！！！！！！！！！！ disp( 'Sen''s Nonparametric Estimator:' ); % calculate slopes ndash = n * ( n - 1 ) / 2; s = zeros( ndash, 1 ); i = 1; for k = 1:n-1, for j = k+1:n, s(i) = ( y(j) - y(k) ) / ( j - k ) / dt; i = i + 1; end; end; % the estimate sl = median( s ); disp( [ ' Slope Estimate = ' num2str( sl ) ] ); m1 = fix( ( ndash - nor * sqrt( va1 ) ) / 2 ); %和原来的trendMK程序相比仅仅把此处的方差改为va1，即修正的方差，其实我不懂。。！！ m2 = fix( ( ndash + nor * sqrt( va1 ) ) / 2 ); %和原来的trendMK程序相比仅仅把此处的方差改为va1，即修正的方差，其实我不懂。。！！！ s = sort( s ); lcl = s( m1 ); ucl = s( m2 + 1 ); disp( [ ' Lower Confidence Limit = ' ... num2str( lcl ) ] ); disp( [ ' Upper Confidence Limit = ' ... num2str( ucl ) ] ); %---------------------------------------------------- % the end