文档库

最新最全的文档下载
当前位置:文档库 > mk突变检测

mk突变检测

寄信人: fig (【晋】『康杰中学』)
标 题: 没主题
发信站: 南京大学小百合站 (Sun May 23 19:18:25 2010)
来 源: 172.16.74.204


Fortran:

program main
implicit none
c This is a program for testing climate jumpi by use of
c 'Mann-Kendall test'.
c-------------------------------------------------
integer,parameter :: iy=100
integer i
real x(iy),u1(iy),u2(iy)
c-------Read Data

c-------
open(31,file='d:mk.txt',form='formatted')
do i=1,50
x(i)=0.5
x(i+50)=-0.5
enddo
c-------Mann-Kendall test method
call MKtest(iy,x,u1,u2)
c------
do i=1,iy
write(31,10)i,x(i),u1(i),u2(i)
end do
10 format(1x,i4,4f8.2)
stop
end
!c------------------------------------------------------------------
subroutine MKtest(m,x,u1,u2)
dimension x(m),u1(m),u2(m)
dimension xr(m),ur(m) !work array
do i=1,m
xr(i)=x(m+1-i)
end do
call MKrank(m,x,u1)
call MKrank(m,xr,ur)
do i=1,m
u2(i)=-ur(m+1-i)
enddo
return
end
!c------------------------------------------------------------------
!c-------Mann-Kendall Rank Statistic
subroutine MKrank(m,x,u)
dimension x(m),u(m)
u(1)=0.0
d=0.0
do i=2,m
num=0
do j=1,i
if(x(i).gt.x(j))then
num=num+1
endif
enddo
d=d+float(num)
e=i*(i-1.)/4.
var=i*(i-1.)*(2.*i+5.)/72.
u(i)=(d-e)/sqrt(var)
enddo
return
end

Matlab:

%
% 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
%
% Bob Newell, February 1996
%
%--------------------------------------------------
%
function [ z, sl, lcl, ucl ] = trendMK( 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 statis

tic = ' 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;
%----------------------------------------------------
%
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 ) ] );
% variance ( assuming no tied groups )
v = ( n * ( n - 1 ) * ( 2 * n + 5 ) ) / 18;
m1 = fix( ( ndash - nor * sqrt( v ) ) / 2 );
m2 = fix( ( ndash + nor * sqrt( v ) ) / 2 );
s = sort( s );
lcl = s( m1 );
ucl = s( m2 + 1 );
disp( [ ' Lower Confidence Limit = ' ...
num2str( lcl ) ] );
disp( [ ' Upper Confidence Limit = ' ...
num2str( ucl ) ] );
%----------------------------------------------------
% the end
--


※ 来源:.南京大学小百合站 http://www.wendangku.net/doc/c6b0eada7f1922791688e8aa.html [FROM: 172.16.74.204]