| Re: code for Kendall trend test [message #58133] |
Tue, 22 January 2008 10:10  |
robinson.inj
Messages: 32
Registered: August 2007
|
Member |
|
|
it can help you
=========
function mkstrend,x,undef=undef
;It uses the non parametric mann-kendal (mk) statistics to determine
if there
;is (or not) a signifcant trend 95%. The slope is calcu-
;lated by using the Sen's method(s). S is robust and less affected by
;outliers (sen p.k. 1968, ASAJ).
;ex. sl=mkstrend(x,undef=-999.).
; where x is your ts: x=[0.5, 0.1, 0.2, 0.1, 0.3]
; no slope then sl=undef. if the ts has one nan: sl=undef
if total(finite(x,/nan)) eq 0 then begin
x=float(x)
nx=n_elements(x)
nx1=nx-1.
n=nx*(nx-1)/2. ; The number of elements in d
d=fltarr(n)
m=0.
for i=0,nx1-1 do begin
for j=i+1,nx-1 do begin
d(m)=x(j)-x(i)
m=m+1
endfor
endfor
for i=0L,n-1 do begin
if d(i) lt 0 then d(i)=-1.
if d(i) eq 0 then d(i)= 0.
if d(i) gt 0 then d(i)= 1.
endfor
s=total(d)
U=x(uniq(x(sort(x))))
Corr=0 ;Correction for tied observations (equal value)
for y=0,n_elements(U)-1 do begin
find=where(x eq U(y))
uj=n_elements(find)
Corr=Corr+uj*(uj-1)*(2*uj+5)
endfor
Vs=(nx*(nx-1.)*(2*nx+5.)-Corr)/18. ;For long series it is
necessary to use the whole
;eq. 2.6 (Corr) (Sen p.k. 1968,
ASAJ)
if s gt 0. then z=(s-1)/sqrt(Vs)
if s lt 0. then z=(s+1)/sqrt(Vs)
if s eq 0. then z=0.
nor=gauss_cvf(0.025) ; Prob at 95% (two-side)
;The slope
Sn=fltarr(n)
m=0.
for i=0,nx1-1 do begin
for j=i+1,nx-1 do begin
Sn(m)=(x(i)-x(j))/(i-j)
m=m+1
endfor
endfor
Snsorted=Sn(sort(Sn))
m=float(fix(n/2.))
if abs(z) lt nor then begin
slope=undef
endif else begin
if 2*m eq n then slope=0.5*(Snsorted(m)+Snsorted(m+1))
if 2*m+1. eq n then slope=Snsorted(m+1)
endelse
endif else begin
slope=undef
endelse
return, slope
end
========
Robinson
On Jan 21, 4:46 pm, txoming...@gmail.com wrote:
> Hi folks,
>
> Does any one have or know a place to find free code in IDL to compute
> the Mann-Kendall trend test and the seasonal Kendall trend test?
>
> Thanks a lot,
>
> Domingo
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
