[R-sig-Geo] Mann-kendall

[Nuno Sá]

Hello

Using the integration by part summation was quite fast (I owe a beer to the
ones who've developed the raster package!!)

So, let variable Z.combi be the Z value of the Mann-kendall.

And PHI(x) be the CDF of the Standar normal distribution w/ mean 0 and var
1.

According to Neeti the p-value is given by:

p = 2 * (1 - PHI(|Z|)

Now, on the code, using the integral expansion:

-----------------------

library(phangorn) #to have double factorials

phi.aprox <- 0.5+exp((-1/2)*abs(Z.combi)^2)/sqrt(2*pi)
phi.integ <- 0

for (n in 1:100){

  phi.integ <- phi.integ + Z.combi^(2*n+1)/dfactorial(2*n+1)
}

phi.aprox <- phi-aprox*phi.integ

p.value.rst.neeti <- 2*(1-phi.aprox)

-------------------

The result falls within 0 to 1 which is a good indication. But all depends
on everything I've asked before.

Thank you!



On 8 October 2014 22:23, Nuno Sá <nunocesardesa at gmail.com> wrote:

> Hello!
>
> I've been trying for two days to implement the Mann-Kendall test. It is
> available in IDRISI and it would be cool that it would be available here.
>
> I believe I made part of it correctly but I am stuck in two steps.
>
> 1. Calculating the ties/NA coefficients. (Algorithm problem..)
>
> 2. Estimating the p-value on the last step (Lack of statistical knowledge
> problem maybe)
>
> So here it goes what I have so far. Bear with me that I have tried to keep
> my variables names with a certain "instant" meaning.
>
> ---------------------- Problem here, is identifying the coefficients for
> the ties/NA -------------------
>
> dummy.stack <- raster(stack.ndvi,layer=1) #required to initiate stack
> within loop
> mtx <- matrix(0,ncol=nlayers(stack.ndvi),nrow=(nlayers(stack.ndvi)-1))
> #this allows the construction of a matrix with tie values signalled to 1
>
> for (i in 1:(nlayers(stack.ndvi)-1)){
>
>   print("----- i ----")
>   print(i)
>   for (j in (i+1):nlayers(stack.ndvi)){
>
>     print(j)
>
>     dummy.stack <-
> stack(dummy.stack,sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)))
>     tie.break <-
> sign(raster(stack.ndvi,layer=i)-raster(stack.ndvi,layer=j)) == 0
>     cc <- maxValue(tie.break)
>     if( cc == 1){
>       mtx[i,j-1]=1
>     }
>   }
> }
>
> S.stack <- dropLayer(dummy.stack,1) #removing dummy
>
> S.stati <- raster(S.stack,layer=1) #for the somatory - initiating in the
> second term (can be substituted by in-built overlay function)
>
> for (i in 2:nlayers(S.stack)){
>
>   S.stati <- S.stati + raster(S.stack,layer=i)
> }
>
> -----------------------------------------------------------------------
>
> #Ignoring the existence of ties/na (which I've read is acceptable for
> large enough datasets)
> #then:
>
>   var.mk.sqr <- sqrt(n*(n-1)*(2*n+5)/18)
>
>   S.c1 <- ((S.stati > 0)*S.stati - 1 )/ var.mk.sqr
>   S.c2 <- (S.stati ==0)*S.stati* 0
>   S.c3 <- ((S.stati < 0)*S.stati + 1 )/ var.mk.sqr
>
>   Z.combi <- overlay(S.c1,S.c2,S.c3,fun="sum")
>
> #"The P value (probability value p) of the MK statistic S of sample data
> can
> be estimated using the normal cumulative distribution function:" (Yue,2002
> and Neeti, 2011)
> # So, how to implement the CDF using raster data? Am I forced to do series
> expansions? would this not be to slow with raster? Or is there a quicker
> way to iterate pixel by pixel?
>
> Sorry for the long email and thank you in advance for any
> comments/corrections/recommendations!
>
> Best regards to all,
> Nuno Sá
>
>
> --
>
> Nuno César de Sá
> +351 91 961 90 37
>
>


-- 

Nuno César de Sá
+351 91 961 90 37

	[[alternative HTML version deleted]]