该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
分享一段用matlab进行MK趋势的代码给你，希望你能用上
%α=0.05
% Time Series Trend Detection Tests
%
% [ z, sl, lcl, ucl ] = trend( y, dt )
%
% where z = 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
%
%
%--------------------------------------------------
%
function [ z, sl, lcl, ucl ] = trend( y, dt )
%
n = length( y );
%--------------------------------------------------
% Mann-Kendall Test for N > 40
%
disp( 'Mann-Kendall Test:' );
if n < 41,
disp( 'WARNING - sould be more than 40 points' );
end;
% calculate statistic
s = 0;
for k = 1:n-1,
for j = k+1:n,
s = s + sign( y(j) - y(k) );
end;
end;
% variance ( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
% test statistic
if s == 0,
z = 0;
elseif s > 0,
z = ( s - 1 ) / sqrt( v );
else,
z = ( s + 1 ) / sqrt( v );
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;
%----------------------------------------------------
%
转播
