From efc482d81936edaa8d59d84eab72e5ee93f8f143 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 22 May 2024 16:39:27 +0200 Subject: [PATCH 1/4] changes to fair rps rpss when probs are provided --- R/RPS.R | 37 ++++++++++++++++++++++++++++--------- R/RPSS.R | 32 ++++++++++++++++++++------------ 2 files changed, 48 insertions(+), 21 deletions(-) diff --git a/R/RPS.R b/R/RPS.R index 59b2d01..0b06b56 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,14 @@ 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 (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 +263,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 +280,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 +360,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 +372,16 @@ 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 + } adjustment <- (-1) / (R - 1) * probs_exp_cumsum * (1 - probs_exp_cumsum) adjustment <- colSums(adjustment) rps[, i, j] <- rps[, i, j] + adjustment diff --git a/R/RPSS.R b/R/RPSS.R index de4e257..51033a4 100644 --- a/R/RPSS.R +++ b/R/RPSS.R @@ -127,7 +127,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, + weights_exp = NULL, weights_ref = NULL, cross.val = FALSE, na.rm = FALSE, sig_method.type = 'two.sided.approx', alpha = 0.05, ncores = NULL) { @@ -403,6 +404,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, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -415,7 +417,8 @@ 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, + Fair = FALSE, nmemb = NULL, + weights_exp = NULL, weights_ref = NULL, cross.val = FALSE, na.rm = FALSE, sig_method.type = 'two.sided.approx', alpha = 0.05) { #--- if memb_dim: # exp: [sdate, memb, (dat)] @@ -493,7 +496,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], @@ -501,17 +505,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, + 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 @@ -560,15 +566,16 @@ 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 + 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 <- dim(probs_clim_cumsum)[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 - # } + adjustment <- (-1) / (R - 1) * probs_clim_cumsum * + (1 - probs_clim_cumsum) + adjustment <- apply(adjustment, 2, sum) + rps_ref <- rps_ref + adjustment + } } } @@ -626,3 +633,4 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', return(list(rpss = rpss, sign = sign)) } + -- GitLab From 4dea52bbb9a6bf59516e2cfb12eb2cba0c680ecc Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 22 May 2024 17:32:25 +0200 Subject: [PATCH 2/4] reference nmembers rpss --- R/RPSS.R | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/R/RPSS.R b/R/RPSS.R index 51033a4..043bc84 100644 --- a/R/RPSS.R +++ b/R/RPSS.R @@ -127,7 +127,7 @@ #'@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, nmemb = 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, ncores = NULL) { @@ -404,7 +404,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 = nmemb, nmemb_ref = nmemb_ref, weights_exp = weights_exp, weights_ref = weights_ref, cross.val = cross.val, @@ -417,7 +417,7 @@ 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, nmemb = 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) { #--- if memb_dim: @@ -505,7 +505,7 @@ 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, nmemb = nmemb, + Fair = Fair, nmemb = nmemb_ref, weights = weights_ref[good_values, , k], na.rm = na.rm, cross.val = cross.val) } @@ -567,16 +567,19 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', 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(probs_clim_cumsum)[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 (!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_exp_cumsum * (1 - probs_exp_cumsum) + adjustment <- colSums(adjustment) + rps[, i, j] <- rps[, i, j] + adjustment } - } } } -- GitLab From 57c13661bc8142685b05a6e8da16197aaa8a874a Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 23 May 2024 12:12:49 +0200 Subject: [PATCH 3/4] fix probs clim in FairRPSS --- R/RPSS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/RPSS.R b/R/RPSS.R index 043bc84..978bc47 100644 --- a/R/RPSS.R +++ b/R/RPSS.R @@ -576,7 +576,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', } else { R <- nmemb_ref } - adjustment <- (-1) / (R - 1) * probs_exp_cumsum * (1 - probs_exp_cumsum) + adjustment <- (-1) / (R - 1) * probs_clim_cumsum * (1 - probs_clim_cumsum) adjustment <- colSums(adjustment) rps[, i, j] <- rps[, i, j] + adjustment } -- GitLab From 7a7bf9e7b0104dbdbf6ed1bc69bc2cca6fc5aaf3 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 23 May 2024 18:18:50 +0200 Subject: [PATCH 4/4] fix typo rps ref --- R/RPSS.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/RPSS.R b/R/RPSS.R index 978bc47..dd47442 100644 --- a/R/RPSS.R +++ b/R/RPSS.R @@ -578,7 +578,7 @@ RPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', } adjustment <- (-1) / (R - 1) * probs_clim_cumsum * (1 - probs_clim_cumsum) adjustment <- colSums(adjustment) - rps[, i, j] <- rps[, i, j] + adjustment + rps_ref[, i, j] <- rps_ref[, i, j] + adjustment } } } -- GitLab