| 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
|
|
|
|
|
|
|
|
| Re: code for Kendall trend test [message #58315 is a reply to message #58133] |
Wed, 23 January 2008 16:01  |
txomingo99
Messages: 3
Registered: January 2008
|
Junior Member |
|
|
Thank you very much, Robinson. This is really helpful!!!.
Anyway, I've been comparing the results with those from the equivalent
test in S-Plus and they are slightly different. I've two questions. I
didn't manage to know where the p-value is calculated/stored and how
to calculate the Confidence interval.
Thanks again! I'm very grateful.
Domingo.
On Jan 22, 1:10 pm, robinson....@gmail.com wrote:
> 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- Hide quoted text -
>
> - Show quoted text -
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
