寄信人: 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 &gt 40

%

disp( 'Mann-Kendall Test:' );

if n &lt 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 &gt 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