tests passed

Author: kongdd

DESCRIPTION

@@ -25,3 +25,6 @@ Config/testthat/edition: 3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.1
+Depends: 
+    R (>= 4.0)
+LazyData: true

NAMESPACE

@@ -1,25 +1,37 @@
 # Generated by roxygen2: do not edit by hand
 
+export(GOF)
+export(KGE)
 export(Ksat)
+export(NSE)
+export(R2_sign)
 export(Wcr_FT)
 export(Wpwp_FT)
 export(XAJ_calib)
 export(XAJ_predict)
 export(bubble)
 export(calib_routing_muskingum_nL)
+export(cv_coef)
 export(expt)
 export(field_capacity)
 export(grid_date)
+export(kurtosis)
 export(plot_calib)
 export(plot_gof)
 export(plot_runoff)
 export(routing_muskingum)
 export(routing_muskingum_nL)
 export(saturated_mois)
+export(skewness)
 export(soil_class)
 export(wilting_point)
 import(magrittr)
+importFrom(data.table,data.table)
+importFrom(dplyr,tibble)
+importFrom(lubridate,day)
 importFrom(lubridate,make_date)
 importFrom(lubridate,make_datetime)
+importFrom(lubridate,month)
+importFrom(lubridate,year)
 importFrom(purrr,map_int)
 importFrom(stats,uniroot)

R/GOF.R

@@ -0,0 +1,286 @@
+#' @param ... ignored
+#' 
+#' @rdname GOF
+#' @export
+NSE <- function(yobs, ysim, w, ...) {
+    if (missing(w)) w <- rep(1, length(yobs))
+
+    ind <- valindex(yobs, ysim)
+    w <- w[ind]
+
+    y_mean <- sum(yobs[ind] * w) / sum(w)
+    # R2: the portion of regression explained variance, also known as
+    # coefficient of determination
+
+    # SSR <- sum((ysim - y_mean)^2 * w)
+    SST <- sum((yobs[ind] - y_mean)^2 * w)
+    # R2     <- SSR / SST
+
+    RE <- ysim[ind] - yobs[ind]
+    # Bias <- sum(w * RE) / sum(w) # bias
+    # Bias_perc <- Bias / y_mean # bias percentage
+    # MAE <- sum(w * abs(RE)) / sum(w) # mean absolute error
+    RMSE <- sqrt(sum(w * (RE)^2) / sum(w)) # root mean sqrt error
+
+    NSE <- 1 - sum((RE)^2 * w) / SST # NSE coefficient
+    NSE
+}
+
+#' @rdname GOF
+#' @export
+KGE <- function(yobs, ysim, ...) {
+  ind <- valindex(yobs, ysim)
+  yobs <- yobs[ind]
+  ysim <- ysim[ind]
+  
+  c1 = cor(yobs, ysim)
+  c2 = sd(ysim) / sd(yobs)
+  c3 = mean(ysim) / mean(yobs)  
+  
+  1 - sqrt((c1 - 1)^2 + (c2 - 1)^2 + (c3 - 1)^2)
+}
+
+#' GOF
+#'
+#' Good of fitting
+#'
+#' @param yobs Numeric vector, observations
+#' @param ysim Numeric vector, corresponding simulated values
+#' @param w Numeric vector, weights of every points. If w included, when
+#' calculating mean, Bias, MAE, RMSE and NSE, w will be taken into considered.
+#' @param include.cv If true, cv will be included.
+#' @param include.r If true, r and R2 will be included.
+#' 
+#' @return
+#' * `RMSE` root mean square error
+#' * `NSE` NASH coefficient
+#' * `MAE` mean absolute error
+#' * `AI` Agreement index (only good points (w == 1)) participate to
+#' calculate. See details in Zhang et al., (2015).
+#' * `Bias` bias
+#' * `Bias_perc` bias percentage
+#' * `n_sim` number of valid obs
+#' * `cv` Coefficient of variation
+#' * `R2` correlation of determination
+#' * `R` pearson correlation
+#' * `pvalue` pvalue of `R`
+#'
+#' @references
+#' 1. https://en.wikipedia.org/wiki/Coefficient_of_determination
+#' 2. https://en.wikipedia.org/wiki/Explained_sum_of_squares
+#' 3. https://en.wikipedia.org/wiki/Nash%E2%80%93Sutcliffe_model_efficiency_coefficient
+#' 4. Zhang Xiaoyang (2015), http://dx.doi.org/10.1016/j.rse.2014.10.012
+#'
+#' @examples
+#' yobs = rnorm(100)
+#' ysim = yobs + rnorm(100)/4
+#' GOF(yobs, ysim)
+#' 
+#' @importFrom dplyr tibble
+#' @export
+GOF <- function(yobs, ysim, w, include.cv = FALSE, include.r = TRUE){
+    if (missing(w)) w <- rep(1, length(yobs))
+
+    # remove NA_real_ and Inf values in ysim, yobs and w
+    valid <- function(x) !is.na(x) & is.finite(x)
+
+    I <- which(valid(ysim) & valid(yobs) & valid(w))
+    # n_obs <- length(yobs)
+    n_sim <- length(I)
+
+    ysim <- ysim[I]
+    yobs <- yobs[I]
+    w     <- w[I]
+
+    if (include.cv) {
+        CV_obs <- cv_coef(yobs, w)
+        CV_sim <- cv_coef(ysim, w)
+    }
+    if (is_empty(yobs)){
+        out <- c(RMSE = NA_real_, 
+            KGE = NA_real_,
+            NSE = NA_real_, MAE = NA_real_, AI = NA_real_,
+            Bias = NA_real_, Bias_perc = NA_real_, n_sim = NA_real_)
+
+        if (include.r) out <- c(out, R2 = NA_real_, R = NA_real_, pvalue = NA_real_)
+        if (include.cv) out <- c(out, obs = CV_obs, sim = CV_sim)
+        return(out)
+    }
+
+    # R2: the portion of regression explained variance, also known as
+    # coefficient of determination
+    KGE = KGE(ysim, yobs)
+    # https://en.wikipedia.org/wiki/Coefficient_of_determination
+    # https://en.wikipedia.org/wiki/Explained_sum_of_squares
+    y_mean <- sum(yobs * w) / sum(w)
+
+    SSR    <- sum( (ysim - y_mean)^2 * w)
+    SST    <- sum( (yobs - y_mean)^2 * w)
+    # R2     <- SSR / SST
+
+    RE     <- ysim - yobs
+    Bias   <- sum ( w*RE)     /sum(w)                     # bias
+    Bias_perc <- Bias/y_mean                              # bias percentage
+    MAE    <- sum ( w*abs(RE))/sum(w)                     # mean absolute error
+    RMSE   <- sqrt( sum(w*(RE)^2)/sum(w) )                # root mean sqrt error
+
+    # https://en.wikipedia.org/wiki/Nash%E2%80%93Sutcliffe_model_efficiency_coefficient
+    NSE  <- 1  - sum( (RE)^2 * w) / SST # NSE coefficient
+
+    # Observations number are not same, so comparing correlation coefficient
+    # was meaningless.
+    # In the current, I have no idea how to add weights `R`.
+    if (include.r){
+        R      <- NA_real_
+        pvalue <- NA_real_
+        
+        tryCatch({
+            cor.obj <- cor.test(yobs, ysim, use = "complete.obs")
+            R       <- cor.obj$estimate[[1]]
+            pvalue  <- cor.obj$p.value
+        }, error = function(e){
+            message(e$message)
+        })
+        R2 = R^2
+    }
+    # In Linear regression, R2 = R^2 (R is pearson cor)
+    # R2     <- summary(lm(ysim ~ yobs))$r.squared # low efficient
+
+    # AI: Agreement Index (only good values(w==1) calculate AI)
+    AI <- NA_real_
+    I2 <- which(w == 1)
+    if (length(I2) >= 2) {
+        yobs = yobs[I2]
+        ysim = ysim[I2]
+        y_mean = mean(yobs)
+        AI = 1 - sum( (ysim - yobs)^2 ) / sum( (abs(ysim - y_mean) + abs(yobs - y_mean))^2 )
+    }
+
+    out <- tibble(R, pvalue, R2, NSE, KGE, RMSE, MAE, 
+             Bias, Bias_perc, AI = AI, n_sim = n_sim)
+    if (include.cv) out <- cbind(out, CV_obs, CV_sim)
+    return(out)
+}
+
+#' skewness and kurtosis
+#'
+#' Inherit from package `e1071`
+#' @param x a numeric vector containing the values whose skewness is to be
+#' computed.
+#' @param na.rm a logical value indicating whether NA values should be stripped
+#' before the computation proceeds.
+#' @param type an integer between 1 and 3 selecting one of the algorithms for
+#' computing skewness.
+#'
+#' @examples
+#' x <- rnorm(100)
+#' coef_kurtosis <- kurtosis(x)
+#' coef_skewness <- skewness(x)
+#' @keywords internal
+#' @export
+kurtosis <- function(x, na.rm = FALSE, type = 3) {
+    if (any(ina <- is.na(x))) {
+        if (na.rm) {
+            x <- x[!ina]
+        } else {
+            return(NA_real_)
+        }
+    }
+    if (!(type %in% (1:3))) stop("Invalid 'type' argument.")
+    n <- length(x)
+    x <- x - mean(x)
+    r <- n * sum(x^4) / (sum(x^2)^2)
+
+    y <- if (type == 1) {
+        r - 3
+    } else if (type == 2) {
+        if (n < 4) {
+            warning("Need at least 4 complete observations.")
+            return(NA_real_)
+        }
+        ((n + 1) * (r - 3) + 6) * (n - 1) / ((n - 2) * (n - 3))
+    } else {
+        r * (1 - 1 / n)^2 - 3
+    }
+    y
+}
+
+#' @rdname kurtosis
+#' @export
+skewness <- function(x, na.rm = FALSE, type = 3) {
+    if (any(ina <- is.na(x))) {
+        if (na.rm) {
+            x <- x[!ina]
+        } else {
+            return(NA_real_)
+        }
+    }
+    if (!(type %in% (1:3))) stop("Invalid 'type' argument.")
+    n <- length(x)
+    x <- x - mean(x)
+    y <- sqrt(n) * sum(x^3) / (sum(x^2)^(3 / 2))
+
+    if (type == 2) {
+        if (n < 3) {
+            warning("Need at least 3 complete observations.")
+            return(NA_real_)
+        }
+        y <- y * sqrt(n * (n - 1)) / (n - 2)
+    } else if (type == 3) {
+        y <- y * ((1 - 1 / n))^(3 / 2)
+    }
+    y
+}
+
+#' weighted CV
+#' @param x Numeric vector
+#' @param w weights of different point
+#'
+#' @return Named numeric vector, (mean, sd, cv).
+#' @examples
+#' x <- rnorm(100)
+#' coefs <- cv_coef(x)
+#' @keywords internal
+#' @export
+cv_coef <- function(x, w) {
+    if (missing(w)) w <- rep(1, length(x))
+    if (length(x) == 0) {
+        return(c(mean = NA_real_, sd = NA_real_, cv = NA_real_))
+    }
+    # rm NA_real_
+    I <- is.finite(x)
+    x <- x[I]
+    w <- w[I]
+
+    mean <- sum(x * w) / sum(w)
+    sd <- sqrt(sum((x - mean)^2 * w) / sum(w))
+    cv <- sd / mean
+    c(mean = mean, sd = sd, cv = cv) # quickly return
+}
+
+#' Critical value of determined correlation
+#'
+#' @param n length of observation.
+#' @param NumberOfPredictor Number of predictor, including constant.
+#' @param alpha significant level.
+#'
+#' @return `F` statistic and `R2` at significant level.
+#'
+#' @keywords internal
+#' @references
+#' Chen Yanguang (2012), Geographical Data analysis with MATLAB.
+#' @examples
+#' R2_critical <- R2_sign(30, NumberOfPredictor = 2, alpha = 0.05)
+#' @export
+R2_sign <- function(n, NumberOfPredictor = 2, alpha = 0.05) {
+  freedom_r <- NumberOfPredictor - 1 # regression
+  freedom_e <- n - NumberOfPredictor # error
+
+  F <- qf(1 - alpha, freedom_r, freedom_e)
+  R2 <- 1 - 1 / (1 + F * freedom_r / freedom_e)
+
+  # F = 485.1
+  # F = R2/freedom_r/((1-R2)/freedom_e)
+  # Rc = sqrt(/(qf(1 - alpha, 1, freedom) + freedom)) %TRUE>% print  # 0.11215
+  return(list(F = F, R2 = R2))
+}

R/HydroALL.R

@@ -1,4 +1,5 @@
 #' @import magrittr
+#' @importFrom data.table data.table
 #' @keywords internal
 "_PACKAGE"