add slope_mk and slope_p

ext/IpaperSlopeExt/IpaperSlopeExt.jl

@@ -7,9 +7,9 @@ using StatsBase: autocor, tiedrank
 using Distributions: ccdf, Normal
 
 using Ipaper
-using Ipaper: lm_resid
-# , trend_mk
+using Ipaper: lm_resid, slope_mk, slope_p, mkTrend
 
-include("mktrend.jl")
+include("mkTrend.jl")
+include("slope_fun.jl")
 
 end

ext/IpaperSlopeExt/mktrend.jl

@@ -1,3 +1,17 @@
+function _slope_sen(y::AbstractVector, x::AbstractVector=1:length(y))
+  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
+
+
 """
     trend_mk(y::AbstractVector{T}; ci=0.95) where {T<:Real}
 
@@ -6,7 +20,7 @@
 - `y`: numeric vector
 - `x`: (optional) numeric vector
 - `ci`: critical value of autocorrelation
- 
+
 # Return
 
 - `Z0`    : The original (non corrected) Mann-Kendall test Z statistic.
@@ -33,12 +47,11 @@
 
 # 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)
+julia> mkTrend([4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69])
 (z0 = 0.35777087639996635, pval0 = 0.7205147871362552, z = 0.35777087639996635, pval = 0.7205147871362552, slope = 0.040000000000000036, intercept = 4.441)
 ```
 """
-function Ipaper.trend_mk(y::AbstractVector{T}; ci=0.95) where {T<:Real}
+function Ipaper.mkTrend(y::AbstractVector, x::AbstractVector; ci=0.95)
   z0 = z = pval0 = pval = slp = intercept = NaN
 
   # y = dropmissing(y)
@@ -92,27 +105,13 @@ function Ipaper.trend_mk(y::AbstractVector{T}; ci=0.95) where {T<:Real}
     z0 = (S + 1) / sqrt(var_S)
   end
 
-  pval = 2 * ccdf(Normal(), abs(z))
-  pval0 = 2 * ccdf(Normal(), abs(z0))
+  pvalue = 2 * ccdf(Normal(), abs(z))
+  pvalue0 = 2 * ccdf(Normal(), abs(z0))
   # Tau = S / (0.5 * n * (n - 1))
-  slope = slope_sen(y)
+  slope = _slope_sen(y, x)
   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)
+  (; z0, pvalue0, z, pvalue, slope, intercept)
 end
 
 # export trend_mk

ext/IpaperSlopeExt/slope_fun.jl

@@ -0,0 +1,48 @@
+using Distributions: TDist, cdf
+using StatsBase: cov, var
+
+
+function pt(x::Float64, df::Int)
+  t = TDist(df)
+  cdf(t, x)
+end
+
+
+"""
+    slope_p(y::AbstractVector, x::AbstractVector=1:length(y))
+
+# Reference
+1. https://zhuanlan.zhihu.com/p/642186978
+
+# Example
+```julia
+x = [1, 2, 3, 4, 5];
+y = [2, 4, 5, 4, 6];
+slope_p(y)
+```
+"""
+function Ipaper.slope_p(y::AbstractVector, x::AbstractVector=1:length(y); ignored...)
+  n = length(x)
+  df = n - 2
+
+  covxy = cov(x, y)
+  varx = var(x) # (x - mean(x))' * (x - mean(x)) / (n - 1)
+  slope = covxy / varx
+  resid = lm_resid(y, x)
+
+  _var = sum((x .- mean(x)) .^ 2)
+  resvar = sum(resid .^ 2) / df
+
+  se = sqrt(resvar / _var)
+  tval = slope / se
+  pvalue = 2 * (1 - pt(abs(tval), df))
+
+  # (; slope, se, tval, pvalue)
+  (; slope, pvalue) #, sd = se
+end
+
+
+function Ipaper.slope_mk(y::AbstractVector, x::AbstractVector=1:length(y); ci=0.95, ignored...)
+  r = mkTrend(y, x; ci)
+  (;slope = r.slope, pvalue = r.pvalue)
+end

src/Slope/Slope.jl

@@ -2,7 +2,10 @@ include("base.jl")
 include("linreg.jl")
 # include("mktrend.jl")
 
-function trend_mk end
+function mkTrend end
+function slope_mk end
+function slope_p end
 
 
-export trend_mk
+export mkTrend,
+  slope_mk, slope_p

test/runtests.jl

@@ -24,4 +24,5 @@ include("test-stat_movmean.jl")
 include("test-stat_anomaly.jl")
 include("test-stat_threshold.jl")
 include("test-stat_warmingLevel.jl")
-# end
+
+include("test-slope.jl")

test/test-slope.jl

@@ -0,0 +1,12 @@
+using Distributions
+using Ipaper
+using Test
+
+@testset "slope_fun" begin
+  x = [4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69]
+  t_mk = slope_mk(x)
+  t_lm = slope_p(x)
+
+  @test t_mk == (slope=0.040000000000000036, pvalue=0.7205147871362552)
+  @test t_lm == (slope=0.04636363636363642, pvalue=0.6249862523021623)
+end