diff --git a/R/RPS.R b/R/RPS.R index 59b2d01a0d842967cdaa4d0351ca1321e19ccf8c..e1b41d88da974b5f8c28e5a9d055484c3031213c 100644 --- a/R/RPS.R +++ b/R/RPS.R @@ -43,6 +43,7 @@ #'@param Fair A logical indicating whether to compute the FairRPS (the #' potential RPS that the forecast would have with an infinite ensemble size). #' The default value is FALSE. +#'@param nmemb A numeric value indicating the number of members used to compute the probabilities. This parameter is necessary when calculating FairRPS from probabilities. Default is NULL. #'@param weights A named numerical array of the weights for 'exp' probability #' calculation. If 'dat_dim' is NULL, the dimensions should include 'memb_dim' #' and 'time_dim'. Else, the dimension should also include 'dat_dim'. The @@ -88,7 +89,8 @@ #'@export RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NULL, dat_dim = NULL, prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, - Fair = FALSE, weights = NULL, cross.val = FALSE, return_mean = TRUE, + Fair = FALSE, nmemb = NULL, weights = NULL, + cross.val = FALSE, return_mean = TRUE, na.rm = FALSE, ncores = NULL) { # Check inputs @@ -178,6 +180,16 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL if (!is.logical(Fair) | length(Fair) > 1) { stop("Parameter 'Fair' must be either TRUE or FALSE.") } + if (Fair) { + if (!is.null(cat_dim)) { + if (cat_dim %in% names(dim(exp))) { + if (is.null(nmemb)) { + stop("Parameter 'nmemb' necessary to compute Fair", + "score from probabilities") + } + } + } + } ## return_mean if (!is.logical(return_mean) | length(return_mean) > 1) { stop("Parameter 'return_mean' must be either TRUE or FALSE.") @@ -253,7 +265,7 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL fun = .RPS, dat_dim = dat_dim, time_dim = time_dim, memb_dim = memb_dim, cat_dim = cat_dim, - prob_thresholds = prob_thresholds, + prob_thresholds = prob_thresholds, nmemb = nmemb, indices_for_clim = indices_for_clim, Fair = Fair, weights = weights, cross.val = cross.val, na.rm = na.rm, ncores = ncores)$output1 @@ -270,7 +282,8 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL .RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NULL, dat_dim = NULL, prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, - Fair = FALSE, weights = NULL, cross.val = FALSE, na.rm = FALSE) { + Fair = FALSE, nmemb = NULL, weights = NULL, + cross.val = FALSE, na.rm = FALSE) { #--- if memb_dim: # exp: [sdate, memb, (dat)] # obs: [sdate, (memb), (dat)] @@ -349,8 +362,11 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL # obs_probs: [bin, sdate] } else { # inputs are probabilities already - exp_probs <- t(exp_data) - obs_probs <- t(obs_data) + if (all(names(dim(exp_data)) == c(time_dim, memb_dim)) || + all(names(dim(exp_data)) == c(time_dim, cat_dim))) { + exp_probs <- t(exp_data) + obs_probs <- t(obs_data) + } } probs_exp_cumsum <- apply(exp_probs, 2, cumsum) @@ -358,11 +374,17 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL # rps: [sdate, nexp, nobs] rps [good_values, i, j] <- colSums((probs_exp_cumsum - probs_obs_cumsum)^2) - 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 + if (!is.null(memb_dim)) { + if (memb_dim %in% names(dim(exp))) { + ## 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 + } + } else { + R <- nmemb + } + warning("Applying fair correction.") adjustment <- (-1) / (R - 1) * probs_exp_cumsum * (1 - probs_exp_cumsum) adjustment <- colSums(adjustment) rps[, i, j] <- rps[, i, j] + adjustment @@ -384,5 +406,3 @@ RPS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', cat_dim = NUL return(rps) } - - diff --git a/R/RPSS.R b/R/RPSS.R index db73a494ee8a1504cbe38a2cc21d8f9f2a22c75c..4c9becdd63a8d48e87017e730a843b12c0661264 100644 --- a/R/RPSS.R +++ b/R/RPSS.R @@ -133,7 +133,8 @@ #'@export RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = NULL, dat_dim = NULL, prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, - Fair = FALSE, weights_exp = NULL, weights_ref = NULL, + Fair = FALSE, nmemb = NULL, nmemb_ref = NULL, + weights_exp = NULL, weights_ref = NULL, cross.val = FALSE, na.rm = FALSE, sig_method.type = 'two.sided.approx', alpha = 0.05, N.eff = NA, ncores = NULL) { @@ -425,8 +426,8 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', if (is.array(N.eff)) { data$N.eff <- N.eff - target_dims[length(target_dims)+1] <- list(NULL) - if (!is.null(ref)){ + target_dims[length(target_dims) + 1] <- list(NULL) + if (!is.null(ref)) { output <- Apply(data, target_dims = target_dims, fun = .RPSS, @@ -434,6 +435,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = dat_dim, prob_thresholds = prob_thresholds, indices_for_clim = indices_for_clim, Fair = Fair, + nmemb = nmemb, nmemb_ref = nmemb_ref, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -448,6 +450,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = dat_dim, prob_thresholds = prob_thresholds, indices_for_clim = indices_for_clim, Fair = Fair, + nmemb = nmemb, nmemb_ref = nmemb_ref, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -455,7 +458,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', ncores = ncores) } } else { # N.eff not an array - if (!is.null(ref)){ + if (!is.null(ref)) { output <- Apply(data, target_dims = target_dims, fun = .RPSS, @@ -463,6 +466,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = dat_dim, prob_thresholds = prob_thresholds, indices_for_clim = indices_for_clim, Fair = Fair, + nmemb = nmemb, nmemb_ref = nmemb_ref, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -477,6 +481,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = dat_dim, prob_thresholds = prob_thresholds, indices_for_clim = indices_for_clim, Fair = Fair, + nmemb = nmemb, nmemb_ref = nmemb_ref, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -489,10 +494,13 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', } -.RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = NULL, - dat_dim = NULL, prob_thresholds = c(1/3, 2/3), indices_for_clim = NULL, - Fair = FALSE, weights_exp = NULL, weights_ref = NULL, cross.val = FALSE, - na.rm = FALSE, sig_method.type = 'two.sided.approx', alpha = 0.05, N.eff = NA) { +.RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', + cat_dim = NULL, dat_dim = NULL, prob_thresholds = c(1/3, 2/3), + indices_for_clim = NULL, Fair = FALSE, + nmemb = NULL, nmemb_ref = NULL, + weights_exp = NULL, weights_ref = NULL, cross.val = FALSE, + na.rm = FALSE, sig_method.type = 'two.sided.approx', + alpha = 0.05, N.eff = NA) { #--- if memb_dim: # exp: [sdate, memb, (dat)] # obs: [sdate, (memb), (dat)] @@ -569,7 +577,8 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = NULL, prob_thresholds = prob_thresholds, indices_for_clim = good_indices_for_clim, - Fair = Fair, weights = weights_exp[good_values, , i], + Fair = Fair, nmemb = nmemb, + weights = weights_exp[good_values, , i], cross.val = cross.val, na.rm = na.rm) rps_ref[good_values, i, j] <- .RPS(exp = ref[good_values, , k], obs = obs[good_values, , j], @@ -577,17 +586,19 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', cat_dim = cat_dim, dat_dim = NULL, prob_thresholds = prob_thresholds, indices_for_clim = good_indices_for_clim, - Fair = Fair, weights = weights_ref[good_values, , k], + Fair = Fair, nmemb = nmemb_ref, + weights = weights_ref[good_values, , k], na.rm = na.rm, cross.val = cross.val) } } } } - } else { # ref is NULL + } else { # ref is NULL rps_exp <- .RPS(exp = exp, obs = obs, time_dim = time_dim, memb_dim = memb_dim, cat_dim = cat_dim, dat_dim = dat_dim, prob_thresholds = prob_thresholds, - indices_for_clim = indices_for_clim, Fair = Fair, weights = weights_exp, + indices_for_clim = indices_for_clim, Fair = Fair, + nmemb = nmemb, weights = weights_exp, cross.val = cross.val, na.rm = na.rm) # RPS of the reference forecast @@ -636,16 +647,20 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', probs_obs_cumsum <- apply(obs_probs, 2, cumsum) rps_ref[good_values, i, j] <- colSums((probs_clim_cumsum - probs_obs_cumsum)^2) } - # 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 - # } - + if (Fair) { # FairRPS + if (!is.null(memb_dim)) { + if (memb_dim %in% names(dim(exp))) { + ## adjustment <- rowSums(-1 * (1/R - 1/R.new) * ens.cum * (R - ens.cum)/R/(R - 1)) + ## [formula taken from SpecsVerification::EnsRps] + R <- dim(obs)[1] #number of years + } + } else { + R <- nmemb_ref + } + adjustment <- (-1) / (R - 1) * probs_clim_cumsum * (1 - probs_clim_cumsum) + adjustment <- colSums(adjustment) + rps_ref[, i, j] <- rps_ref[, i, j] + adjustment + } } } } @@ -708,3 +723,4 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', return(list(rpss = rpss, sign = sign)) } +