Mann-kendall

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]]