add mktrend

src/Ipaper.jl

@@ -46,6 +46,7 @@ include("stringr.jl")
 include("missing.jl")
 include("Statistics/Statistics.jl")
 include("Climate/Climate.jl")
+include("Slope/Slope.jl")
 
 include("precompile.jl")
 include("tools_plot.jl")

src/Slope/Slope.jl

@@ -0,0 +1,3 @@
+include("base.jl")
+include("linreg.jl")
+include("mktrend.jl")

src/Slope/base.jl

@@ -0,0 +1,12 @@
+using LinearAlgebra
+# https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33
+
+function valid_input(y::AbstractVector, x::AbstractVector) 
+  inds = @.(!isnan(y) && !isnan(x))
+
+  y = @view y[inds]
+  x = @view x[inds]
+  x, y
+end
+
+

src/Statistics/linreg.jl → src/Slope/linreg.jl

@@ -1,14 +1,3 @@
-using LinearAlgebra
-# https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33
-
-function valid_input(y::AbstractVector{T}, x::AbstractVector{T}) where {T<:Real}
-  inds = @.(!isnan(y) && !isnan(x))
-
-  y = @view y[inds]
-  x = @view x[inds]
-  x, y
-end
-
 # function linreg1(y, X)
 #   β_hat = (X' * X) \ X' * y
 #   return (β_hat)
@@ -20,9 +9,9 @@ end
 # end
 
 """
-$(TYPEDSIGNATURES)
+    linreg_simple(y::AbstractVector, x::AbstractVector; na_rm=false) 
 """
-function linreg_simple(y::AbstractVector{T}, x::AbstractVector{T}; na_rm=false) where {T<:Real}
+function linreg_simple(y::AbstractVector, x::AbstractVector; na_rm=false) 
   (length(x)) == length(y) || throw(DimensionMismatch())
 
   if na_rm
@@ -32,12 +21,14 @@ function linreg_simple(y::AbstractVector{T}, x::AbstractVector{T}; na_rm=false)
   x̄ = mean(x)
   ȳ = mean(y)
 
-  sum_a = T(0.0)
-  sum_b = T(0.0)
-  @inbounds for i = 1:length(y)
+  sum_a = 0.0
+  sum_b = 0.0
+
+  @inbounds for i = eachindex(y)
     sum_a += (x[i] - x̄) * (y[i] - ȳ)
     sum_b += (x[i] - x̄) * (x[i] - x̄)
   end
+
   β1 = sum_a / sum_b
   β0 = ȳ - β1 * x̄
   # β1 = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄) .^ 2)
@@ -51,9 +42,9 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+    linreg_fast(y::AbstractVector, x::AbstractVector; na_rm=false)
 """
-function linreg_fast(y::AbstractVector{T}, x::AbstractVector{T}; na_rm=false) where {T<:Real}
+function linreg_fast(y::AbstractVector, x::AbstractVector; na_rm=false) 
   (length(x)) == length(y) || throw(DimensionMismatch())
 
   if na_rm
@@ -62,7 +53,7 @@ function linreg_fast(y::AbstractVector{T}, x::AbstractVector{T}; na_rm=false) wh
 
   N = length(y)
   ldiv!(
-    cholesky!(Symmetric([T(N) sum(x); zero(T) sum(abs2, x)], :U)),
+    cholesky!(Symmetric([N sum(x); 0 sum(abs2, x)], :U)),
     [sum(y), dot(x, y)])
 end
 
@@ -73,4 +64,13 @@ function linreg_fast(y::AbstractVector{T}; kw...) where {T<:Real}
 end
 
 
+function lm_resid(y::AbstractVector, x::AbstractVector)
+  β0, β1 = linreg_simple(y, x)
+  ysim = β0 .+ β1 .* x
+  y .- ysim
+end
+
 lm = linreg = linreg_fast
+
+export lm, linreg, linreg_fast, linreg_simple
+export lm_resid

src/Slope/mktrend.jl

@@ -0,0 +1,124 @@
+using Statistics: mean, median, quantile
+using StatsBase: autocor, tiedrank
+using Distributions: ccdf, Normal
+
+
+"""
+    trend_mk(y::AbstractVector{T}; ci=0.95) where {T<:Real}
+
+# Arguments
+
+- `y`: numeric vector
+- `x`: (optional) numeric vector
+- `ci`: critical value of autocorrelation
+ 
+# Return
+
+- `Z0`    : The original (non corrected) Mann-Kendall test Z statistic.
+
+- `pval0` : The original (non corrected) Mann-Kendall test p-value
+
+- `Z`     : The new Z statistic after applying the correction
+
+- `pval`  : Corrected p-value after accounting for serial autocorrelation
+  `N/n*s` Value of the correction factor, representing the quotient of the
+  number of samples N divided by the effective sample size `n*s`
+
+- `slp`   : Sen slope, The slope of the (linear) trend according to Sen test.
+  slp is significant, if pval < alpha.
+
+# References
+
+1. Hipel, K.W. and McLeod, A.I. (1994), \emph{Time Series Modelling of Water
+   Resources and Environmental Systems}. New York: Elsevier Science.
+
+2. Libiseller, C. and Grimvall, A., (2002), Performance of partial Mann-Kendall
+   tests for trend detection in the presence of covariates.
+   \emph{Environmetrics} 13, 71--84, \doi{10.1002/env.507}.
+
+# Example
+```julia
+julia> x = [4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69];
+julia> trend_mk(x)
+(z0 = 0.35777087639996635, pval0 = 0.7205147871362552, z = 0.35777087639996635, pval = 0.7205147871362552, slope = 0.040000000000000036, intercept = 4.441)
+```
+"""
+function trend_mk(y::AbstractVector{T}; ci=0.95) where {T<:Real}
+  z0 = z = pval0 = pval = slp = intercept = NaN
+
+  # y = dropmissing(y)
+  n = length(y)
+  if n < 5
+    return (; z0, pval0, z, pval, slp, intercept)
+  end
+
+  S = 0
+  for i in 1:(n-1)
+    for j in (i+1):n
+      S += sign(y[j] - y[i])
+    end
+  end
+
+  sig = quantile(Normal(), (1 + ci) / 2) / sqrt(n) # qnorm((1 + ci)/2)/sqrt(n)
+
+  rank = tiedrank(lm_resid(y, 1:n))
+  ro = autocor(rank, 1:n-1)
+  ro[abs.(ro) .<= sig] .= 0.0 # modified by dongdong Kong, 2017-04-03
+
+  cte = 2 / (n * (n - 1) * (n - 2))
+
+  ess = 0.0
+  for i in 1:n-1
+    ess += (n - i) * (n - i - 1) * (n - i - 2) * ro[i]
+  end
+
+  essf = 1 + ess * cte
+  var_S = n * (n - 1) * (2n + 5) * (1.0 / 18)
+
+  if length(unique(y)) < n
+    aux = unique(y)
+    for i in eachindex(aux)
+      tie = count(y .== aux[i])
+      if tie > 1
+        var_S -= tie * (tie - 1) * (2tie + 5) * (1 / 18)
+      end
+    end
+  end
+
+  VS = var_S * essf
+  if S == 0
+    z = 0.0
+    z0 = 0.0
+  elseif S > 0
+    z = (S - 1) / sqrt(VS)
+    z0 = (S - 1) / sqrt(var_S)
+  else
+    z = (S + 1) / sqrt(VS)
+    z0 = (S + 1) / sqrt(var_S)
+  end
+
+  pval = 2 * ccdf(Normal(), abs(z))
+  pval0 = 2 * ccdf(Normal(), abs(z0))
+  # Tau = S / (0.5 * n * (n - 1))
+  slope = slope_sen(y)
+  intercept = mean(y .- slope .* (1:n))
+
+  (; z0, pval0, z, pval, slope, intercept)
+end
+
+
+function slope_sen(y::AbstractVector{T}, x=1:length(y)) where {T<:Real}
+  n = length(x)
+  V = fill(NaN, Int((n^2 - n) / 2))
+  k = 0
+  for i in 2:n
+    for j in 1:(i-1)
+      k += 1
+      V[k] = (y[i] - y[j]) / (x[i] - x[j])
+    end
+  end
+  median(V)
+end
+
+
+export trend_mk

src/Statistics/Statistics.jl

@@ -1,7 +1,6 @@
 import Statistics
 using Statistics: mean, median
 
-include("linreg.jl")
 include("movmean.jl")
 include("apply.jl")
 include("NanQuantile.jl")
@@ -10,4 +9,3 @@ include("NanQuantile.jl")
 export mean, median, quantile, movmean, weighted_mean, weighted_movmean
 export apply
 export _nanquantile!, nanquantile, NanQuantile
-export lm, linreg, linreg_fast, linreg_simple