| code for Kendall trend test [message #58163] |
Mon, 21 January 2008 14:46  |
txomingo99
Messages: 3
Registered: January 2008
|
Junior Member |
|
|
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 #58313 is a reply to message #58163] |
Wed, 23 January 2008 23:21  |
robinson.inj
Messages: 32
Registered: August 2007
|
Member |
|
|
On Jan 23, 6:01 pm, txoming...@gmail.com wrote:
> 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 -- Hide quoted text -
>
> - Show quoted text -
nor=gauss_cvf(0.025) ; Prob at 95% (two-side)
use the "help" to know more about gauss_cvf
Robinson
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
