diff --git a/DESCRIPTION b/DESCRIPTION
index f2577bfeebb720f2d76e89bde01f267e54ade2ae..9690f44d6f6c97aa416a016cc1a46a9a0803f84d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -37,9 +37,9 @@ Imports:
     NbClust,
     multiApply (>= 2.1.1),
     SpecsVerification (>= 0.5.0),
-    easyNCDF
+    easyNCDF,
+    easyVerification
 Suggests:
-    easyVerification,
     testthat
 License: Apache License 2.0
 URL: https://earth.bsc.es/gitlab/es/s2dv/
diff --git a/NAMESPACE b/NAMESPACE
index 841814e89a1726a34ad3ad755b8fc24191eebfa8..1449ad4d9722d001583f7337f77d47ff90c91d1f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -54,6 +54,8 @@ export(ProjectField)
 export(REOF)
 export(RMS)
 export(RMSSS)
+export(RPS)
+export(RPSS)
 export(RandomWalkTest)
 export(RatioPredictableComponents)
 export(RatioRMS)
@@ -92,6 +94,7 @@ importFrom(ClimProjDiags,WeightedMean)
 importFrom(abind,abind)
 importFrom(abind,adrop)
 importFrom(easyNCDF,ArrayToNc)
+importFrom(easyVerification,convert2prob)
 importFrom(grDevices,bmp)
 importFrom(grDevices,col2rgb)
 importFrom(grDevices,colorRampPalette)
diff --git a/R/RPS.R b/R/RPS.R
new file mode 100644
index 0000000000000000000000000000000000000000..84956ac85a2ee39250e793ec80dab9f26df75eaa
--- /dev/null
+++ b/R/RPS.R
@@ -0,0 +1,193 @@
+#'Compute the Ranked Probability Score
+#'
+#'The Ranked Probability Score (RPS; Wilks, 2011) is defined as the sum of the
+#'squared differences between the cumulative forecast probabilities (computed 
+#'from the ensemble members) and the observations (defined as 0% if the category
+#'did not happen and 100% if it happened). It can be used to evaluate the skill 
+#'of multi-categorical probabilistic forecasts. The RPS ranges between 0 
+#'(perfect forecast) and n-1 (worst possible forecast), where n is the number of
+#'categories. In the case of a forecast divided into two categories (the lowest
+#'number of categories that a probabilistic forecast can have), the RPS 
+#'corresponds to the Brier Score (BS; Wilks, 2011), therefore, ranges between 0 
+#'and 1.
+#'
+#'@param exp A named numerical array of the forecast with at least time 
+#'  dimension.  
+#'@param obs A named numerical array of the observation with at least time 
+#'  dimension. The dimensions must be the same as 'exp' except 'memb_dim'.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'sdate'.
+#'@param memb_dim A character string indicating the name of the member dimension
+#'  to compute the probabilities of the forecast. The default value is 'member'.
+#'@param prob_thresholds A numeric vector of the relative thresholds (from 0 to
+#'  1) between the categories. The default value is c(1/3, 2/3), which 
+#'  corresponds to tercile equiprobable categories.
+#'@param indices_for_clim A vector of the indices to be taken along 'time_dim' 
+#'  for computing the thresholds between the probabilistic categories. If NULL,
+#'  the whole period is used. The default value is NULL.
+#'@param Fair A logical indicating whether to compute the FairRPSS (the 
+#'  potential RPSS that the forecast would have with an infinite ensemble size).
+#'  The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return 
+#'A numerical array of RPS with the same dimensions as "exp" except the 
+#''time_dim' and 'memb_dim' dimensions.
+#'
+#'@references 
+#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+#'
+#'@examples
+#'exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+#'obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
+#'res <- RPS(exp = exp, obs = obs)
+#'
+#'@import multiApply
+#'@importFrom easyVerification convert2prob
+#'@export
+RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', 
+                prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, Fair = FALSE, 
+                ncores = NULL) {
+  
+  # Check inputs
+  ## exp and obs (1)
+  if (!is.array(exp) | !is.numeric(exp))
+    stop('Parameter "exp" must be a numeric array.')
+  if (!is.array(obs) | !is.numeric(obs))
+    stop('Parameter "obs" must be a numeric array.')
+  if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+     any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1)
+    stop('Parameter "time_dim" must be a character string.')
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) > 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!memb_dim %in% names(dim(exp))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  name_exp <- name_exp[-which(name_exp == memb_dim)]
+  if (memb_dim %in% name_obs) {
+    name_obs <- name_obs[-which(name_obs == memb_dim)]
+  }
+  if (!identical(length(name_exp), length(name_obs)) |
+      !identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimensions expect 'memb_dim'."))
+  }
+  ## prob_thresholds
+  if (!is.numeric(prob_thresholds) | !is.vector(prob_thresholds) |
+      any(prob_thresholds <= 0) | any(prob_thresholds >= 1)) {
+    stop("Parameter 'prob_thresholds' must be a numeric vector between 0 and 1.")
+  }
+  ## indices_for_clim
+  if (is.null(indices_for_clim)) {
+    indices_for_clim <- 1:dim(obs)[time_dim]
+  } else {
+    if (!is.numeric(indices_for_clim) | !is.vector(indices_for_clim)) {
+      stop("Parameter 'indices_for_clim' must be NULL or a numeric vector.")
+    } else if (length(indices_for_clim) > dim(obs)[time_dim] |
+               max(indices_for_clim) > dim(obs)[time_dim] |
+               any(indices_for_clim) < 1) {
+      stop("Parameter 'indices_for_clim' should be the indices of 'time_dim'.")
+    }
+  }
+  ## Fair
+  if (!is.logical(Fair)  | length(Fair) > 1) {
+    stop("Parameter 'Fair' must be either TRUE or FALSE.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be either NULL or a positive integer.")
+    }
+  }
+
+  ###############################
+ 
+  # Compute RPS
+  if (!memb_dim %in% names(dim(obs))) {
+    target_dims_obs <- time_dim
+  } else {
+    target_dims_obs <- c(time_dim, memb_dim)
+  }
+
+  rps <- Apply(data = list(exp = exp, obs = obs), 
+               target_dims = list(exp = c(time_dim, memb_dim), 
+                                  obs = target_dims_obs),
+               output_dims = time_dim,
+               fun = .RPS,
+               prob_thresholds = prob_thresholds, 
+               indices_for_clim = indices_for_clim, Fair = Fair,
+               ncores = ncores)$output1
+
+  # Return only the mean RPS
+  rps <- MeanDims(rps, time_dim, na.rm = FALSE)
+ 
+  return(rps)
+}
+
+
+.RPS <- function(exp, obs, prob_thresholds = c(1/3, 2/3), 
+                 indices_for_clim = NULL, Fair = FALSE) {
+  # exp: [sdate, memb]
+  # obs: [sdate, (memb)]
+
+  exp_probs <- .get_probs(data = exp, indices_for_quantiles = indices_for_clim, 
+                          prob_thresholds = prob_thresholds)
+  # exp_probs: [bin, sdate] 
+  obs_probs <- .get_probs(data = obs, indices_for_quantiles = indices_for_clim, 
+                          prob_thresholds = prob_thresholds)
+  # obs_probs: [bin, sdate]
+
+  probs_exp_cumsum <- apply(exp_probs, 2, cumsum)
+  probs_obs_cumsum <- apply(obs_probs, 2, cumsum)
+  rps <- apply((probs_exp_cumsum - probs_obs_cumsum)^2, 2, sum)
+  # rps: [sdate]
+
+  if (Fair) { # FairRPS
+    ## adjustment <- rowSums(-1 * (1/R - 1/R.new) * ens.cum * (R - ens.cum)/R/(R - 1)) [formula taken from SpecsVerification::EnsRps]
+    R <- dim(exp)[2]  #memb
+    R_new <- Inf
+    adjustment <- (-1) / (R - 1) * probs_exp_cumsum * (1 - probs_exp_cumsum)
+    adjustment <- apply(adjustment, 2, sum)
+    rps <- rps + adjustment
+  }
+
+  return(rps)
+}
+
+.get_probs <- function(data, indices_for_quantiles, prob_thresholds) {
+  # if exp: [sdate, memb]
+  # if obs: [sdate, (memb)]
+
+  # Add dim [memb = 1] to obs if it doesn't have memb_dim
+  if (length(dim(data)) == 1) dim(data) <- c(dim(data), 1)
+
+  # Absolute thresholds
+  quantiles <- quantile(x = as.vector(data[indices_for_quantiles, ]), probs = prob_thresholds, type = 8, na.rm = TRUE)
+
+  # Probabilities
+  probs <- array(dim = c(bin = length(quantiles) + 1, dim(data)[1])) # [bin, sdate]
+  for (i_time in 1:dim(data)[1]) {
+    if (anyNA(data[i_time, ])) {
+      probs[, i_time] <- rep(NA, dim = length(quantiles) + 1)
+    } else {
+      probs[, i_time] <- colMeans(easyVerification::convert2prob(data[i_time, ], threshold = quantiles))
+    }
+  }
+
+  return(probs)
+}
+
diff --git a/R/RPSS.R b/R/RPSS.R
new file mode 100644
index 0000000000000000000000000000000000000000..04c31378f6e1d733d83a5fb1ea1e5d297612973b
--- /dev/null
+++ b/R/RPSS.R
@@ -0,0 +1,223 @@
+#'Compute the Ranked Probability Skill Score
+#'
+#'The Ranked Probability Skill Score (RPSS; Wilks, 2011) is the skill score 
+#'based on the Ranked Probability Score (RPS; Wilks, 2011). It can be used to 
+#'assess whether a forecast presents an improvement or worsening with respect to
+#'a reference forecast. The RPSS ranges between minus infinite and 1. If the 
+#'RPSS is positive, it indicates that the forecast has higher skill than the 
+#'reference forecast, while a negative value means that it has a lower skill. 
+#'Examples of reference forecasts are the climatological forecast (same 
+#'probabilities for all categories for all time steps), persistence, a previous
+#'model version, and another model. It is computed as RPSS = 1 - RPS_exp / RPS_ref.
+#'The statistical significance is obtained based on a Random Walk test at the 
+#'95% confidence level (DelSole and Tippett, 2016).
+#'
+#'@param exp A named numerical array of the forecast with at least time 
+#'  dimension.
+#'@param obs A named numerical array of the observation with at least time 
+#'  dimension. The dimensions must be the same as 'exp' except 'memb_dim'. 
+#'@param ref A named numerical array of the reference forecast data with at 
+#'  least time dimension. The dimensions must be the same as 'exp' except 
+#'  'memb_dim'. If it is NULL, the climatological forecast is used as reference
+#'  forecast. The default value is NULL.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'sdate'.
+#'@param memb_dim A character string indicating the name of the member dimension
+#'  to compute the probabilities of the forecast and the reference forecast. The
+#'  default value is 'member'.
+#'@param prob_thresholds A numeric vector of the relative thresholds (from 0 to
+#'  1) between the categories. The default value is c(1/3, 2/3), which 
+#'  corresponds to tercile equiprobable categories.
+#'@param indices_for_clim A vector of the indices to be taken along 'time_dim' 
+#'  for computing the thresholds between the probabilistic categories. If NULL,
+#'  the whole period is used. The default value is NULL.
+#'@param Fair A logical indicating whether to compute the FairRPSS (the 
+#'  potential RPSS that the forecast would have with an infinite ensemble size).
+#'  The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return
+#'\item{$rpss}{
+#'  A numerical array of the RPSS with the same dimensions as "exp" except the
+#'  'time_dim' and 'memb_dim' dimensions.
+#'}
+#'\item{$sign}{
+#'  A logical array of the statistical significance of the RPSS with the same
+#'  dimensions as 'exp' except the 'time_dim' and 'memb_dim' dimensions.
+#'}
+#'
+#'@references 
+#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+#'DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1
+#'
+#'@examples
+#'exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+#'obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
+#'ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+#'res <- RPSS(exp = exp, obs = obs) ## climatology as reference forecast
+#'res <- RPSS(exp = exp, obs = obs, ref = ref) ## ref as reference forecast
+#'
+#'@import multiApply
+#'@export
+RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member',
+                 prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, Fair = FALSE, 
+                 ncores = NULL) {
+ 
+  # Check inputs
+  ## exp, obs, and ref (1)
+  if (!is.array(exp) | !is.numeric(exp))
+    stop('Parameter "exp" must be a numeric array.')
+  if (!is.array(obs) | !is.numeric(obs))
+    stop('Parameter "obs" must be a numeric array.')
+  if (!is.null(ref)) {
+    if (!is.array(ref) | !is.numeric(ref))
+      stop('Parameter "ref" must be a numeric array.')
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1)
+    stop('Parameter "time_dim" must be a character string.')
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  if (!is.null(ref) & !time_dim %in% names(dim(ref))) {
+    stop("Parameter 'time_dim' is not found in 'ref' dimension.")
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) > 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!memb_dim %in% names(dim(exp))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+  }
+  if (!is.null(ref) & !memb_dim %in% names(dim(ref))) {
+    stop("Parameter 'memb_dim' is not found in 'ref' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  name_exp <- name_exp[-which(name_exp == memb_dim)]
+  if (memb_dim %in% name_obs) {
+    name_obs <- name_obs[-which(name_obs == memb_dim)]
+  }
+  if (!identical(length(name_exp), length(name_obs)) |
+      !identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimensions expect 'memb_dim'."))
+  }
+  if (!is.null(ref)) {
+    name_ref <- sort(names(dim(ref)))
+    name_ref <- name_ref[-which(name_ref == memb_dim)]
+    if (!identical(length(name_exp), length(name_ref)) |
+        !identical(dim(exp)[name_exp], dim(ref)[name_ref])) {
+      stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                  "all dimensions expect 'memb_dim'."))
+    }
+  }
+  ## prob_thresholds
+  if (!is.numeric(prob_thresholds) | !is.vector(prob_thresholds) |
+      any(prob_thresholds <= 0) | any(prob_thresholds >= 1)) {
+    stop("Parameter 'prob_thresholds' must be a numeric vector between 0 and 1.")
+  }
+  ## indices_for_clim
+  if (is.null(indices_for_clim)) {
+    indices_for_clim <- 1:dim(obs)[time_dim]
+  } else {
+    if (!is.numeric(indices_for_clim) | !is.vector(indices_for_clim)) {
+      stop("Parameter 'indices_for_clim' must be NULL or a numeric vector.")
+    } else if (length(indices_for_clim) > dim(obs)[time_dim] |
+               max(indices_for_clim) > dim(obs)[time_dim] |
+               any(indices_for_clim) < 1) {
+      stop("Parameter 'indices_for_clim' should be the indices of 'time_dim'.")
+    }
+  }
+  ## Fair
+  if (!is.logical(Fair)  | length(Fair) > 1) {
+    stop("Parameter 'Fair' must be either TRUE or FALSE.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be either NULL or a positive integer.")
+    }
+  }
+
+  ###############################
+
+  # Compute RPSS
+  if (!memb_dim %in% names(dim(obs))) {
+    target_dims_obs <- time_dim
+  } else {
+    target_dims_obs <- c(time_dim, memb_dim)
+  }
+
+  if (!is.null(ref)) { # use "ref" as reference forecast
+    data <- list(exp = exp, obs = obs, ref = ref)
+    target_dims = list(exp = c(time_dim, memb_dim),
+                       obs = target_dims_obs,
+                       ref = c(time_dim, memb_dim))
+  } else {
+    data <- list(exp = exp, obs = obs)
+    target_dims = list(exp = c(time_dim, memb_dim),
+                       obs = target_dims_obs)
+  }
+  output <- Apply(data,
+                  target_dims = target_dims,
+                  fun = .RPSS,
+                  prob_thresholds = prob_thresholds,
+                  indices_for_clim = indices_for_clim, Fair = Fair,
+                  ncores = ncores)
+  
+  return(output)
+}
+
+.RPSS <- function(exp, obs, ref = NULL, prob_thresholds = c(1/3, 2/3), 
+                  indices_for_clim = NULL, Fair = FALSE) {
+  # exp: [sdate, memb]
+  # obs: [sdate, (memb)]
+  # ref: [sdate, memb] or NULL
+
+  # RPS of the forecast
+  rps_exp <- .RPS(exp = exp, obs = obs, prob_thresholds = prob_thresholds, 
+                  indices_for_clim = indices_for_clim, Fair = Fair)
+  
+  # RPS of the reference forecast
+  if (is.null(ref)) { ## using climatology as reference forecast
+    obs_probs <- .get_probs(data = obs, indices_for_quantiles = indices_for_clim, 
+                            prob_thresholds = prob_thresholds)
+    # obs_probs: [bin, sdate] 
+
+    clim_probs <- c(prob_thresholds[1], diff(prob_thresholds), 1 - prob_thresholds[length(prob_thresholds)])
+    clim_probs <- array(clim_probs, dim = dim(obs_probs))
+    # clim_probs: [bin, sdate]  
+
+    # Calculate RPS for each time step
+    probs_clim_cumsum <- apply(clim_probs, 2, cumsum)
+    probs_obs_cumsum <- apply(obs_probs, 2, cumsum)
+    rps_ref <- apply((probs_clim_cumsum - probs_obs_cumsum)^2, 2, sum)
+    # rps_ref: [sdate]
+
+#    if (Fair) { # FairRPS
+#      ## adjustment <- rowSums(-1 * (1/R - 1/R.new) * ens.cum * (R - ens.cum)/R/(R - 1)) [formula taken from SpecsVerification::EnsRps]
+#      R <- dim(exp)[2]  #memb
+#      R_new <- Inf
+#      adjustment <- (-1) / (R - 1) * probs_clim_cumsum * (1 - probs_clim_cumsum)
+#      adjustment <- apply(adjustment, 2, sum)
+#      rps_ref <- rps_ref + adjustment
+#    }
+
+  } else { # use "ref" as reference forecast
+    rps_ref <- .RPS(exp = ref, obs = obs, prob_thresholds = prob_thresholds, 
+                    indices_for_clim = indices_for_clim, Fair = Fair)
+  }
+  
+  # RPSS
+  rpss <- 1 - mean(rps_exp) / mean(rps_ref)
+  
+  # Significance
+  sign <- .RandomWalkTest(skill_A = rps_exp, skill_B = rps_ref)$signif
+  
+  return(list(rpss = rpss, sign = sign))
+}
+
diff --git a/R/RandomWalkTest.R b/R/RandomWalkTest.R
index 494be6520dba4d6aedff0449b23e59a321e00732..adeadc1ec94b62920c885640938f966c91e75ddc 100644
--- a/R/RandomWalkTest.R
+++ b/R/RandomWalkTest.R
@@ -1,8 +1,8 @@
 #'Random walk test for skill differences
 #'
 #'Forecast comparison of the skill obtained with 2 forecasts (with respect to a 
-#'common reference) based on Random Walks, with significance estimate at the 5%
-#'confidence level, as in DelSole and Tippett (2015).
+#'common reference) based on Random Walks, with significance estimate at the 95%
+#'confidence level, as in DelSole and Tippett (2016).
 #'
 #'@param skill_A A numerical array of the time series of the skill with the 
 #'  forecaster A's.
diff --git a/R/ResidualCorr.R b/R/ResidualCorr.R
index dae0cb861ebcccb9132af3c629d7c987611cb416..7428d1b08be1ab462cbd37750e82835a2c91f797 100644
--- a/R/ResidualCorr.R
+++ b/R/ResidualCorr.R
@@ -17,11 +17,11 @@
 #'
 #'@param exp A named numerical array of the forecast with at least time 
 #'  dimension. 
+#'@param obs A named numerical array of the observations with at least time 
+#'  dimension. The dimensions must be the same as "exp" except 'memb_dim'.
 #'@param ref A named numerical array of the reference forecast data with at 
 #'  least time dimension. The dimensions must be the same as "exp" except 
 #'  'memb_dim'.
-#'@param obs A named numerical array of the observations with at least time 
-#'  dimension. The dimensions must be the same as "exp" except 'memb_dim'.
 #'@param N.eff Effective sample size to be used in the statistical significance
 #'  test. It can be NA (and it will be computed with the s2dv:::.Eno), a 
 #'  numeric (which is used for all cases), or an array with the same dimensions
@@ -66,13 +66,13 @@
 #'
 #'@examples
 #' exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
-#' ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 5, sdate = 50))
 #' obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
-#' res <- ResidualCorr(exp, ref, obs, memb_dim = 'member')
+#' ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 5, sdate = 50))
+#' res <- ResidualCorr(exp = exp, obs = obs, ref = ref, memb_dim = 'member')
 #'
 #'@import multiApply
 #'@export
-ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
+ResidualCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate',
                          memb_dim = NULL, method = 'pearson', alpha = NULL, 
                          handle.na = 'return.na', ncores = NULL) {
 
@@ -80,10 +80,10 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
   ## exp, ref, and obs (1)
   if (!is.array(exp) | !is.numeric(exp))
     stop('Parameter "exp" must be a numeric array.')
-  if (!is.array(ref) | !is.numeric(ref))
-    stop('Parameter "ref" must be a numeric array.')
   if (!is.array(obs) | !is.numeric(obs))
     stop('Parameter "obs" must be a numeric array.')
+  if (!is.array(ref) | !is.numeric(ref))
+    stop('Parameter "ref" must be a numeric array.')
 
   ## N.eff
   if (is.array(N.eff)) {
@@ -101,9 +101,9 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) != 1)
     stop('Parameter "time_dim" must be a character string.')
-  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(ref)) | 
-      !time_dim %in% names(dim(obs))) {
-    stop("Parameter 'time_dim' is not found in 'exp', 'ref', or 'obs' dimension.")
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs)) | 
+      !time_dim %in% names(dim(ref))) {
+    stop("Parameter 'time_dim' is not found in 'exp', 'obs', or 'ref' dimension.")
   }
   ## memb_dim
   if (!is.null(memb_dim)) {
@@ -116,15 +116,15 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
   }
   ## exp, ref, and obs (2)
   name_exp <- sort(names(dim(exp)))
-  name_ref <- sort(names(dim(ref)))
   name_obs <- sort(names(dim(obs)))
+  name_ref <- sort(names(dim(ref)))
   if (!is.null(memb_dim)) {
     name_exp <- name_exp[-which(name_exp == memb_dim)]
     name_ref <- name_ref[-which(name_ref == memb_dim)]
   }
-  if (!identical(length(name_exp), length(name_ref), length(name_obs)) |
-      !identical(dim(exp)[name_exp], dim(ref)[name_ref], dim(obs)[name_obs])) {
-    stop(paste0("Parameter 'exp', 'ref', and 'obs' must have same length of ",
+  if (length(name_exp) != length(name_obs) | length(name_exp) != length(name_ref) |
+      any(dim(exp)[name_exp] != dim(obs)[name_obs]) | any(dim(exp)[name_exp] != dim(ref)[name_ref])) {
+    stop(paste0("Parameter 'exp', 'obs', and 'ref' must have same length of ",
                 "all dimensions expect 'memb_dim'."))
   }
   ## method
@@ -153,13 +153,13 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
   ###############################
 
   # NA check: na.fail
-  if (handle.na == "na.fail" & (anyNA(exp) | anyNA(ref) | anyNA(obs)))
+  if (handle.na == "na.fail" & (anyNA(exp) | anyNA(obs) | anyNA(ref)))
     stop('The data contain NAs.')
 
   # Calculate ensemble mean
   dim_exp <- dim(exp)
-  dim_ref <- dim(ref)
   dim_obs <- dim(obs)
+  dim_ref <- dim(ref)
 
   if (!is.null(memb_dim)) {
     exp_memb_dim_ind <- which(names(dim_exp) == memb_dim)
@@ -180,18 +180,17 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
 
   # Residual correlation
   if (is.array(N.eff)) {
-    output <- Apply(data = list(exp = exp, ref = ref,
-                                obs = obs, N.eff = N.eff),
-                    target_dims = list(exp = time_dim, ref = time_dim,
-                                       obs = time_dim, N.eff = NULL),
+    output <- Apply(data = list(exp = exp, obs = obs, ref = ref,
+                                N.eff = N.eff),
+                    target_dims = list(exp = time_dim, obs = time_dim,
+                                       ref = time_dim, N.eff = NULL),
                     output_dims = output_dims,
                     fun = .ResidualCorr, method = method,
                     alpha = alpha, handle.na = handle.na, ncores = ncores)
   } else { 
-    output <- Apply(data = list(exp = exp, ref = ref, 
-                                obs = obs),
-                    target_dims = list(exp = time_dim, ref = time_dim, 
-                                       obs = time_dim), 
+    output <- Apply(data = list(exp = exp, obs = obs, ref = ref),
+                    target_dims = list(exp = time_dim, obs = time_dim, 
+                                       ref = time_dim), 
                     output_dims = output_dims, N.eff = N.eff,
                     fun = .ResidualCorr, method = method, 
                     alpha = alpha, handle.na = handle.na, ncores = ncores)
@@ -200,9 +199,9 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
   return(output)
 }
 
-.ResidualCorr <- function(exp, ref, obs, N.eff, method, alpha, handle.na) {
+.ResidualCorr <- function(exp, obs, ref, N.eff, method, alpha, handle.na) {
   # exp and ref and obs: [time]
-  .residual.corr <- function(exp, ref, obs, method, N.eff, alpha) {
+  .residual.corr <- function(exp, obs, ref, method, N.eff, alpha) {
 
     # Residuals of 'exp' and 'obs' (regressing 'ref' out in both 'exp' and 'obs')
     exp_res <- lm(formula = y ~ x, data = list(y = exp, x = ref), na.action = NULL)$residuals
@@ -235,16 +234,16 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
 
   #==================================================
   
-  if (anyNA(exp) | anyNA(ref) | anyNA(obs)) { ## There are NAs 
+  if (anyNA(exp) | anyNA(obs) | anyNA(ref)) { ## There are NAs 
     if (handle.na == 'only.complete.triplets') {
-      nna <- is.na(exp) | is.na(ref) | is.na(obs) # A vector of T/F
+      nna <- is.na(exp) | is.na(obs) | is.na(ref) # A vector of T/F
       if (all(nna)) stop("There is no complete set of forecasts and observations.")
       # Remove the incomplete set
       exp <- exp[!nna]
-      ref <- ref[!nna]
       obs <- obs[!nna]
+      ref <- ref[!nna]
 
-      output <- .residual.corr(exp = exp, ref = ref, obs = obs, method = method, 
+      output <- .residual.corr(exp = exp, obs = obs, ref = ref, method = method, 
                                N.eff = N.eff, alpha = alpha)
       
     } else if (handle.na == 'return.na') {
@@ -257,7 +256,7 @@ ResidualCorr <- function(exp, ref, obs, N.eff = NA, time_dim = 'sdate',
     }
     
   } else { ## There is no NA  
-    output <- .residual.corr(exp = exp, ref = ref, obs = obs, method = method, 
+    output <- .residual.corr(exp = exp, obs = obs, ref = ref, method = method, 
                          N.eff = N.eff, alpha = alpha)
   }
 
diff --git a/man/RPS.Rd b/man/RPS.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d1acc0282797fe813543f0ba47298bf0e1801148
--- /dev/null
+++ b/man/RPS.Rd
@@ -0,0 +1,70 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RPS.R
+\name{RPS}
+\alias{RPS}
+\title{Compute the Ranked Probability Score}
+\usage{
+RPS(
+  exp,
+  obs,
+  time_dim = "sdate",
+  memb_dim = "member",
+  prob_thresholds = c(1/3, 2/3),
+  indices_for_clim = NULL,
+  Fair = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numerical array of the forecast with at least time 
+dimension.}
+
+\item{obs}{A named numerical array of the observation with at least time 
+dimension. The dimensions must be the same as 'exp' except 'memb_dim'.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'sdate'.}
+
+\item{memb_dim}{A character string indicating the name of the member dimension
+to compute the probabilities of the forecast. The default value is 'member'.}
+
+\item{prob_thresholds}{A numeric vector of the relative thresholds (from 0 to
+1) between the categories. The default value is c(1/3, 2/3), which 
+corresponds to tercile equiprobable categories.}
+
+\item{indices_for_clim}{A vector of the indices to be taken along 'time_dim' 
+for computing the thresholds between the probabilistic categories. If NULL,
+the whole period is used. The default value is NULL.}
+
+\item{Fair}{A logical indicating whether to compute the FairRPSS (the 
+potential RPSS that the forecast would have with an infinite ensemble size).
+The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A numerical array of RPS with the same dimensions as "exp" except the 
+'time_dim' and 'memb_dim' dimensions.
+}
+\description{
+The Ranked Probability Score (RPS; Wilks, 2011) is defined as the sum of the
+squared differences between the cumulative forecast probabilities (computed 
+from the ensemble members) and the observations (defined as 0% if the category
+did not happen and 100% if it happened). It can be used to evaluate the skill 
+of multi-categorical probabilistic forecasts. The RPS ranges between 0 
+(perfect forecast) and n-1 (worst possible forecast), where n is the number of
+categories. In the case of a forecast divided into two categories (the lowest
+number of categories that a probabilistic forecast can have), the RPS 
+corresponds to the Brier Score (BS; Wilks, 2011), therefore, ranges between 0 
+and 1.
+}
+\examples{
+exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
+res <- RPS(exp = exp, obs = obs)
+
+}
+\references{
+Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+}
diff --git a/man/RPSS.Rd b/man/RPSS.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..6324ea2e1bf8537d5dff33ea137ea1a8b41e1ab1
--- /dev/null
+++ b/man/RPSS.Rd
@@ -0,0 +1,87 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RPSS.R
+\name{RPSS}
+\alias{RPSS}
+\title{Compute the Ranked Probability Skill Score}
+\usage{
+RPSS(
+  exp,
+  obs,
+  ref = NULL,
+  time_dim = "sdate",
+  memb_dim = "member",
+  prob_thresholds = c(1/3, 2/3),
+  indices_for_clim = NULL,
+  Fair = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numerical array of the forecast with at least time 
+dimension.}
+
+\item{obs}{A named numerical array of the observation with at least time 
+dimension. The dimensions must be the same as 'exp' except 'memb_dim'.}
+
+\item{ref}{A named numerical array of the reference forecast data with at 
+least time dimension. The dimensions must be the same as 'exp' except 
+'memb_dim'. If it is NULL, the climatological forecast is used as reference
+forecast. The default value is NULL.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'sdate'.}
+
+\item{memb_dim}{A character string indicating the name of the member dimension
+to compute the probabilities of the forecast and the reference forecast. The
+default value is 'member'.}
+
+\item{prob_thresholds}{A numeric vector of the relative thresholds (from 0 to
+1) between the categories. The default value is c(1/3, 2/3), which 
+corresponds to tercile equiprobable categories.}
+
+\item{indices_for_clim}{A vector of the indices to be taken along 'time_dim' 
+for computing the thresholds between the probabilistic categories. If NULL,
+the whole period is used. The default value is NULL.}
+
+\item{Fair}{A logical indicating whether to compute the FairRPSS (the 
+potential RPSS that the forecast would have with an infinite ensemble size).
+The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+\item{$rpss}{
+ A numerical array of the RPSS with the same dimensions as "exp" except the
+ 'time_dim' and 'memb_dim' dimensions.
+}
+\item{$sign}{
+ A logical array of the statistical significance of the RPSS with the same
+ dimensions as 'exp' except the 'time_dim' and 'memb_dim' dimensions.
+}
+}
+\description{
+The Ranked Probability Skill Score (RPSS; Wilks, 2011) is the skill score 
+based on the Ranked Probability Score (RPS; Wilks, 2011). It can be used to 
+assess whether a forecast presents an improvement or worsening with respect to
+a reference forecast. The RPSS ranges between minus infinite and 1. If the 
+RPSS is positive, it indicates that the forecast has higher skill than the 
+reference forecast, while a negative value means that it has a lower skill. 
+Examples of reference forecasts are the climatological forecast (same 
+probabilities for all categories for all time steps), persistence, a previous
+model version, and another model. It is computed as RPSS = 1 - RPS_exp / RPS_ref.
+The statistical significance is obtained based on a Random Walk test at the 
+95% confidence level (DelSole and Tippett, 2016).
+}
+\examples{
+exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
+ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
+res <- RPSS(exp = exp, obs = obs) ## climatology as reference forecast
+res <- RPSS(exp = exp, obs = obs, ref = ref) ## ref as reference forecast
+
+}
+\references{
+Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1
+}
diff --git a/man/RandomWalkTest.Rd b/man/RandomWalkTest.Rd
index 11106487e6ab0e82b7299c0f378ebbb7247d3a55..bd1460a74fc8c9cbe69f810ec6ef6c23497dcbd6 100644
--- a/man/RandomWalkTest.Rd
+++ b/man/RandomWalkTest.Rd
@@ -37,8 +37,8 @@ A list of 2:
 }
 \description{
 Forecast comparison of the skill obtained with 2 forecasts (with respect to a 
-common reference) based on Random Walks, with significance estimate at the 5%
-confidence level, as in DelSole and Tippett (2015).
+common reference) based on Random Walks, with significance estimate at the 95%
+confidence level, as in DelSole and Tippett (2016).
 }
 \examples{
 fcst_A <- array(c(11:50), dim = c(sdate = 10, lat = 2, lon = 2))
diff --git a/man/ResidualCorr.Rd b/man/ResidualCorr.Rd
index e98ea7a439d56191eeca8b2d135924c401c23ad9..fe7dd1012f5f81a4feaf574891b2468b79e88572 100644
--- a/man/ResidualCorr.Rd
+++ b/man/ResidualCorr.Rd
@@ -6,8 +6,8 @@
 \usage{
 ResidualCorr(
   exp,
-  ref,
   obs,
+  ref,
   N.eff = NA,
   time_dim = "sdate",
   memb_dim = NULL,
@@ -21,13 +21,13 @@ ResidualCorr(
 \item{exp}{A named numerical array of the forecast with at least time 
 dimension.}
 
+\item{obs}{A named numerical array of the observations with at least time 
+dimension. The dimensions must be the same as "exp" except 'memb_dim'.}
+
 \item{ref}{A named numerical array of the reference forecast data with at 
 least time dimension. The dimensions must be the same as "exp" except 
 'memb_dim'.}
 
-\item{obs}{A named numerical array of the observations with at least time 
-dimension. The dimensions must be the same as "exp" except 'memb_dim'.}
-
 \item{N.eff}{Effective sample size to be used in the statistical significance
 test. It can be NA (and it will be computed with the s2dv:::.Eno), a 
 numeric (which is used for all cases), or an array with the same dimensions
@@ -95,8 +95,8 @@ effective degrees of freedom to account for the time series' autocorrelation
 }
 \examples{
 exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
-ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 5, sdate = 50))
 obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
-res <- ResidualCorr(exp, ref, obs, memb_dim = 'member')
+ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 5, sdate = 50))
+res <- ResidualCorr(exp = exp, obs = obs, ref = ref, memb_dim = 'member')
 
 }
diff --git a/tests/testthat/test-RPS.R b/tests/testthat/test-RPS.R
new file mode 100644
index 0000000000000000000000000000000000000000..e5a5f0c9f3027cbdb5a3d2cb7183688bf5e3e32f
--- /dev/null
+++ b/tests/testthat/test-RPS.R
@@ -0,0 +1,138 @@
+context("s2dv::RPS tests")
+
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(60), dim = c(member = 3, sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+set.seed(2)
+obs2 <- array(rnorm(10), dim = c(sdate = 10))
+
+# dat2_1
+set.seed(1)
+exp2_1 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+set.seed(2)
+obs2_1 <- array(rnorm(10), dim = c(member = 1, sdate = 10))
+
+
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs (1)
+  expect_error(
+  RPS(c()),
+  'Parameter "exp" must be a numeric array.'
+  )
+  expect_error(
+  RPS(exp1, c()),
+  'Parameter "obs" must be a numeric array.'
+  )
+
+  # time_dim
+  expect_error(
+  RPS(exp1, obs1, time_dim = 1),
+  'Parameter "time_dim" must be a character string.'
+  )
+  expect_error(
+  RPS(exp1, obs1, time_dim = 'time'),
+  "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # memb_dim
+  expect_error(
+  RPS(exp1, obs1, memb_dim = TRUE),
+  "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+  RPS(exp1, obs1, memb_dim = 'memb'),
+  "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  # exp, ref, and obs (2)
+  expect_error(
+  RPS(exp1, array(1:9, dim = c(sdate = 9))),
+  "Parameter 'exp' and 'obs' must have same length of all dimensions expect 'memb_dim'."
+  )
+  # prob_thresholds
+  expect_error(
+  RPS(exp1, obs1, prob_thresholds = 1),
+  "Parameter 'prob_thresholds' must be a numeric vector between 0 and 1."
+  )
+  # indices_for_clim
+  expect_error(
+  RPS(exp1, obs1, indices_for_clim = array(1:6, dim = c(2, 3))),
+  "Parameter 'indices_for_clim' must be NULL or a numeric vector."
+  )
+  expect_error(
+  RPS(exp1, obs1, indices_for_clim = 3:11),
+  "Parameter 'indices_for_clim' should be the indices of 'time_dim'."
+  )
+  # Fair
+  expect_error(
+  RPS(exp1, obs1, Fair = 1),
+  "Parameter 'Fair' must be either TRUE or FALSE."
+  )
+  # ncores
+  expect_error(
+  RPS(exp2, obs2, ncores = 1.5),
+  "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+dim(RPS(exp1, obs1)),
+c(lat = 2)
+)
+expect_equal(
+as.vector(RPS(exp1, obs1)),
+c(0.3555556, 0.4444444),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPS(exp1, obs1, Fair = T)),
+c(0.2000000, 0.2666667),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPS(exp1, obs1, indices_for_clim = 3:5)),
+c(0.5000000, 0.4666667),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPS(exp1, obs1, prob_thresholds = seq(0.1, 0.9, 0.1))),
+c(1.600000, 1.888889),
+tolerance = 0.0001
+)
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+expect_equal(
+dim(RPS(exp2, obs2)),
+NULL
+)
+expect_equal(
+RPS(exp2, obs2),
+0.75,
+tolerance = 0.0001
+)
+expect_equal(
+RPS(exp2, obs2, indices_for_clim = 2:5, prob_thresholds = seq(0.1, 0.9, 0.1)),
+2.75,
+tolerance = 0.0001
+)
+expect_equal(
+RPS(exp2, obs2),
+RPS(exp2_1, obs2_1)
+)
+})
diff --git a/tests/testthat/test-RPSS.R b/tests/testthat/test-RPSS.R
new file mode 100644
index 0000000000000000000000000000000000000000..b4c0b7a4fa5edb8804e08348ce4ec042de9adc4f
--- /dev/null
+++ b/tests/testthat/test-RPSS.R
@@ -0,0 +1,276 @@
+context("s2dv::RPSS tests")
+
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(60), dim = c(member = 3, sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+set.seed(3)
+ref1 <- array(rnorm(40), dim = c(member = 2, sdate = 10, lat = 2))
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+set.seed(2)
+obs2 <- array(rnorm(10), dim = c(sdate = 10))
+set.seed(2)
+obs2_1 <- array(rnorm(10), dim = c(member = 1, sdate = 10))
+set.seed(3)
+ref2 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs (1)
+  expect_error(
+  RPSS(c()),
+  'Parameter "exp" must be a numeric array.'
+  )
+  expect_error(
+  RPSS(exp1, c()),
+  'Parameter "obs" must be a numeric array.'
+  )
+
+  # time_dim
+  expect_error(
+  RPSS(exp1, obs1, time_dim = 1),
+  'Parameter "time_dim" must be a character string.'
+  )
+  expect_error(
+  RPSS(exp1, obs1, time_dim = 'time'),
+  "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+  RPSS(exp2, obs2, array(rnorm(20), dim = c(member = 2, time = 10))),
+  "Parameter 'time_dim' is not found in 'ref' dimension."
+  )
+  # memb_dim
+  expect_error(
+  RPSS(exp1, obs1, memb_dim = TRUE),
+  "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+  RPSS(exp1, obs1, memb_dim = 'memb'),
+  "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  expect_error(
+  RPSS(exp2, obs2, array(rnorm(20), dim = c(memb = 2, sdate = 10))),
+  "Parameter 'memb_dim' is not found in 'ref' dimension."
+  )
+  # exp, ref, and obs (2)
+  expect_error(
+  RPSS(exp1, array(1:9, dim = c(sdate = 9))),
+  "Parameter 'exp' and 'obs' must have same length of all dimensions expect 'memb_dim'."
+  )
+  expect_error(
+  RPSS(exp1, obs1, ref2),
+  "Parameter 'exp' and 'obs' must have same length of all dimensions expect 'memb_dim'."
+  )
+  # prob_thresholds
+  expect_error(
+  RPSS(exp1, obs1, ref1, prob_thresholds = 1),
+  "Parameter 'prob_thresholds' must be a numeric vector between 0 and 1."
+  )
+  # indices_for_clim
+  expect_error(
+  RPSS(exp1, obs1, ref1, indices_for_clim = array(1:6, dim = c(2, 3))),
+  "Parameter 'indices_for_clim' must be NULL or a numeric vector."
+  )
+  expect_error(
+  RPSS(exp1, obs1, indices_for_clim = 3:11),
+  "Parameter 'indices_for_clim' should be the indices of 'time_dim'."
+  )
+  # Fair
+  expect_error(
+  RPSS(exp1, obs1, Fair = 1),
+  "Parameter 'Fair' must be either TRUE or FALSE."
+  )
+  # ncores
+  expect_error(
+  RPSS(exp2, obs2, ncores = 1.5),
+  "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+names(RPSS(exp1, obs1)),
+c("rpss", "sign")
+)
+expect_equal(
+names(RPSS(exp1, obs1, ref1)),
+c("rpss", "sign")
+)
+expect_equal(
+dim(RPSS(exp1, obs1)$rpss),
+c(lat = 2)
+)
+expect_equal(
+dim(RPSS(exp1, obs1)$sign),
+c(lat = 2)
+)
+expect_equal(
+dim(RPSS(exp1, obs1, ref1)$rpss),
+c(lat = 2)
+)
+expect_equal(
+dim(RPSS(exp1, obs1, ref1)$sign),
+c(lat = 2)
+)
+# ref = NULL
+expect_equal(
+as.vector(RPSS(exp1, obs1)$rpss),
+c(0.15789474, -0.05263158),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1)$sign),
+c(FALSE, FALSE),
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, Fair = T)$rpss),
+c(0.5263158, 0.3684211),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, indices_for_clim = 3:5)$rpss),
+c(-0.4062500, -0.1052632),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, prob_thresholds = seq(0.1, 0.9, 0.1))$rpss),
+c(0.03030303, -0.14478114),
+tolerance = 0.0001
+)
+# ref = ref
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1)$rpss),
+c(0.5259259, 0.4771242),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1)$sign),
+c(FALSE, FALSE)
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1, Fair = T)$rpss),
+c(0.6000000, 0.6190476),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1, indices_for_clim = 3:5)$rpss),
+c(0.2857143, 0.4166667),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1, indices_for_clim = 3:5)$sign),
+c(FALSE, FALSE)
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1, prob_thresholds = seq(0.1, 0.9, 0.1))$rpss),
+c(0.4754098, 0.3703704),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp1, obs1, ref1, prob_thresholds = seq(0.1, 0.9, 0.1))$sign),
+c(T, F)
+)
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+expect_equal(
+names(RPSS(exp2, obs2)),
+c("rpss", "sign")
+)
+expect_equal(
+names(RPSS(exp2, obs2, ref2)),
+c("rpss", "sign")
+)
+expect_equal(
+dim(RPSS(exp2, obs2)$rpss),
+NULL
+)
+expect_equal(
+dim(RPSS(exp2, obs2)$sign),
+NULL
+)
+expect_equal(
+dim(RPSS(exp2, obs2, ref2)$rpss),
+NULL
+)
+expect_equal(
+dim(RPSS(exp2, obs2, ref2)$sign),
+NULL
+)
+# ref = NULL
+expect_equal(
+as.vector(RPSS(exp2, obs2)$rpss),
+c(-0.7763158),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2)$sign),
+FALSE,
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, Fair = T)$rpss),
+c(-0.1842105),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, indices_for_clim = 3:5)$rpss),
+-0.8984375,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, prob_thresholds = seq(0.1, 0.9, 0.1))$rpss),
+c(-0.7272727),
+tolerance = 0.0001
+)
+# ref = ref
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2)$rpss),
+0,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2)$sign),
+FALSE
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2, Fair = T)$rpss),
+0,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2, indices_for_clim = 3:5)$rpss),
+0.03571429,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2, indices_for_clim = 3:5)$sign),
+FALSE
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2, prob_thresholds = seq(0.1, 0.9, 0.1))$rpss),
+0.06557377,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(RPSS(exp2, obs2, ref2, prob_thresholds = seq(0.1, 0.9, 0.1))$sign),
+FALSE
+)
+expect_equal(
+RPSS(exp2, obs2),
+RPSS(exp2, obs2_1)
+)
+
+})
diff --git a/tests/testthat/test-ResidualCorr.R b/tests/testthat/test-ResidualCorr.R
index c9ac4b9cee0bcd9b02de29622e124a6709966b99..c15276c93368322eec48e8b33657474fb8f4a1c4 100644
--- a/tests/testthat/test-ResidualCorr.R
+++ b/tests/testthat/test-ResidualCorr.R
@@ -29,63 +29,63 @@ test_that("1. Input checks", {
   'Parameter "exp" must be a numeric array.'
   )
   expect_error(
-  ResidualCorr(exp1, c('a')),
+  ResidualCorr(exp1, obs1, c('a')),
   'Parameter "ref" must be a numeric array.'
   )
   expect_error(
-  ResidualCorr(exp1, ref1, c()),
+  ResidualCorr(exp1, c()),
   'Parameter "obs" must be a numeric array.'
   )
   # N.eff
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, N.eff = array(1, dim = c(lat = 2, lon = 2))),
+  ResidualCorr(exp1, obs1, ref1, N.eff = array(1, dim = c(lat = 2, lon = 2))),
 'If parameter "N.eff" is provided with an array, it must have the same dimensions as "obs" except "time_dim".'
   )
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, N.eff = 1:3),
+  ResidualCorr(exp1, obs1, ref1, N.eff = 1:3),
   'Parameter "N.eff" must be NA, a numeric, or an array with the same dimensions as "obs" except "time_dim".'
   )
   # time_dim
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, time_dim = 1),
+  ResidualCorr(exp1, obs1, ref1, time_dim = 1),
   'Parameter "time_dim" must be a character string.'
   )
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, time_dim = 'time'),
-  "Parameter 'time_dim' is not found in 'exp', 'ref', or 'obs' dimension."
+  ResidualCorr(exp1, obs1, ref1, time_dim = 'time'),
+  "Parameter 'time_dim' is not found in 'exp', 'obs', or 'ref' dimension."
   )
   # memb_dim
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, memb_dim = TRUE),
+  ResidualCorr(exp1, obs1, ref1, memb_dim = TRUE),
   "Parameter 'memb_dim' must be a character string."
   )
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, memb_dim = 'member'),
+  ResidualCorr(exp1, obs1, ref1, memb_dim = 'member'),
   "Parameter 'memb_dim' is not found in 'exp' or 'ref' dimension."
   )
   # exp, ref, and obs (2)
   expect_error(
-  ResidualCorr(exp1, array(1:10, dim = c(sdate = 10, memb = 1)), obs1, memb_dim = 'memb'),
-  "Parameter 'exp', 'ref', and 'obs' must have same length of all dimensions expect 'memb_dim'."
+  ResidualCorr(exp1, obs1, array(1:10, dim = c(sdate = 10, memb = 1)), memb_dim = 'memb'),
+  "Parameter 'exp', 'obs', and 'ref' must have same length of all dimensions expect 'memb_dim'."
   )
   # method
   expect_error(
-  ResidualCorr(exp1, ref1, obs1, method = 'asd', memb_dim = 'memb'),
+  ResidualCorr(exp1, obs1, ref1, method = 'asd', memb_dim = 'memb'),
   'Parameter "method" must be "pearson", "kendall", or "spearman".'
   )
   # alpha
   expect_error(
-  ResidualCorr(exp2, ref2, obs2, alpha = 1),
+  ResidualCorr(exp2, obs2, ref2, alpha = 1),
   'Parameter "alpha" must be NULL or a number between 0 and 1.'
   )
   # handle.na
   expect_error(
-  ResidualCorr(exp2, ref2, obs2, handle.na = TRUE),
+  ResidualCorr(exp2, obs2, ref2, handle.na = TRUE),
   'Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".'
   )
   # ncores
   expect_error(
-  ResidualCorr(exp2, ref2, obs2, ncores = 1.5),
+  ResidualCorr(exp2, obs2, ref2, ncores = 1.5),
   'Parameter "ncores" must be either NULL or a positive integer.'
   )
 
@@ -94,69 +94,69 @@ test_that("1. Input checks", {
 ##############################################
 test_that("2. Output checks: dat1", {
 expect_equal(
-names(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')),
+names(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')),
 c("res.corr", "p.val")
 )
 expect_equal(
-dim(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')$res.corr),
+dim(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')$res.corr),
 c(lat = 3, lon = 2)
 )
 expect_equal(
-dim(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')$p.val),
+dim(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')$p.val),
 c(lat = 3, lon = 2)
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')$res.corr),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')$res.corr),
 c(0.002784318, 0.537697647, -0.240071018, 0.258706464, 0.338160748, 0.432107476),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')$p),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')$p),
 c(0.49695468, 0.05446055, 0.25203961, 0.23522967, 0.16960864, 0.10618145),
 tolerance = 0.0001
 )
 expect_equal(
-names(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', alpha = 0.05)),
+names(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)),
 c("res.corr", "sign")
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', alpha = 0.05)$res.corr),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$res.corr),
 c(0.002784318, 0.537697647, -0.240071018, 0.258706464, 0.338160748, 0.432107476),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', alpha = 0.05)$sign),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$sign),
 rep(FALSE, 6)
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', method = "kendall")$res),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "kendall")$res),
 c(0.11111111, 0.42222222, -0.20000000, 0.02222222, 0.20000000, 0.42222222),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', method = "kendall")$p),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "kendall")$p),
 c(0.3799615, 0.1120891, 0.2897920, 0.4757064, 0.2897920, 0.1120891),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', method = "spearman")$res),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$res),
 c(0.13939394, 0.61212121, -0.27272727, 0.07878788, 0.28484848, 0.47878788),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', method = "spearman")$p),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$p),
 c(0.35046594, 0.02998607, 0.22291917, 0.41435870, 0.21251908, 0.08076146),
 tolerance = 0.0001
 )
 
 
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', N.eff = Neff1)$res),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', N.eff = Neff1)$res),
 c(0.002784318, 0.537697647, -0.240071018, 0.258706464, 0.338160748, 0.432107476),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', N.eff = Neff1)$p),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', N.eff = Neff1)$p),
 c(0.49716394, 0.05446055, 0.26690777, 0.23522967, 0.18671025, 0.10618145),
 tolerance = 0.0001
 )
@@ -165,17 +165,17 @@ tolerance = 0.0001
 #---------------------------
 exp1[1] <- NA
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb')$res.corr),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb')$res.corr),
 c(NA, 0.537697647, -0.240071018, 0.258706464, 0.338160748, 0.432107476),
 tolerance = 0.0001
 )
 expect_equal(
-as.vector(ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', handle.na = 'only.complete.triplets')$res.corr),
+as.vector(ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', handle.na = 'only.complete.triplets')$res.corr),
 c(-0.1082588, 0.537697647, -0.240071018, 0.258706464, 0.338160748, 0.432107476),
 tolerance = 0.0001
 )
 expect_error(
-ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', handle.na = 'na.fail'),
+ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', handle.na = 'na.fail'),
 "The data contain NAs."
 )
 
@@ -183,7 +183,7 @@ exp1[1,1,1,1:3] <- NA
 ref1[1,1,1,4:8] <- NA
 obs1[1,1,9:10] <- NA
 expect_error(
-ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', handle.na = 'only.complete.triplets'),
+ResidualCorr(exp1, obs1, ref1, memb_dim = 'memb', handle.na = 'only.complete.triplets'),
 "There is no complete set of forecasts and observations."
 )
 
@@ -194,24 +194,24 @@ ResidualCorr(exp1, ref1, obs1, memb_dim = 'memb', handle.na = 'only.complete.tri
 
 test_that("3. Output checks: dat2", {
 expect_equal(
-names(ResidualCorr(exp2, ref2, obs2)),
+names(ResidualCorr(exp2, obs2, ref2)),
 c("res.corr", "p.val")
 )
 expect_equal(
-dim(ResidualCorr(exp2, ref2, obs2)$res.corr),
+dim(ResidualCorr(exp2, obs2, ref2)$res.corr),
 NULL
 )
 expect_equal(
-dim(ResidualCorr(exp2, ref2, obs2)$p),
+dim(ResidualCorr(exp2, obs2, ref2)$p),
 NULL
 )
 expect_equal(
-ResidualCorr(exp2, ref2, obs2)$p,
+ResidualCorr(exp2, obs2, ref2)$p,
 0.2209847,
 tolerance = 0.0001
 )
 expect_equal(
-ResidualCorr(exp2, ref2, obs2)$res,
+ResidualCorr(exp2, obs2, ref2)$res,
 -0.274962,
 tolerance = 0.0001
 )