From 24e6a2d7055acb604de9302bb59d64cd903fc3ba Mon Sep 17 00:00:00 2001 From: Carlos Delgado Date: Wed, 11 Jan 2023 16:25:31 +0100 Subject: [PATCH 1/2] cross-validation for climatlogical forecast --- R/CRPSS.R | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/R/CRPSS.R b/R/CRPSS.R index a6b4a14..d31a264 100644 --- a/R/CRPSS.R +++ b/R/CRPSS.R @@ -222,11 +222,20 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', if (is.null(ref)) { ## using climatology as reference forecast ## all the time steps are used as if they were members - ## then, ref dimensions are [sdate, memb], both with length(sdate) + ## then, ref dimensions are [sdate, memb] + ## memb dimension has length(sdate)-1 due to cross-validation obs_time_len <- dim(obs)[time_dim] if (is.null(dat_dim)) { - ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + + ## Without cross-validation: + ## ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + ## With cross-validation (excluding the value of that year to create ref for that year): + ref <- array(data = NA, dim = c(obs_time_len, obs_time_len-1)) + for (i in 1:obs_time_len){ + ref[i,] <- obs[-i] + } + names(dim(ref)) <- c(time_dim, memb_dim) # ref: [sdate, memb]; obs: [sdate] crps_ref <- .CRPS(exp = ref, obs = obs, time_dim = time_dim, memb_dim = memb_dim, @@ -237,7 +246,15 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', crps_ref <- array(dim = c(obs_time_len, nobs)) names(dim(crps_ref)) <- c(time_dim, 'nobs') for (i_obs in 1:nobs) { - ref <- array(data = rep(obs[, i_obs], each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + + ## Without cross-validation: + ## ref <- array(data = rep(obs[, i_obs], each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + ## With cross-validation (excluding the value of that year to create ref for that year): + ref <- array(data = NA, dim = c(obs_time_len, obs_time_len-1)) + for (i in 1:obs_time_len){ + ref[i,] <- obs[-i, i_obs] + } + names(dim(ref)) <- c(time_dim, memb_dim) crps_ref[, i_obs] <- .CRPS(exp = ref, obs = ClimProjDiags::Subset(obs, dat_dim, i_obs, drop = 'selected'), time_dim = time_dim, memb_dim = memb_dim, dat_dim = NULL, Fair = Fair) -- GitLab From 27b8767d5ca5ac34c81c7a480b3c1cf07b2d2f35 Mon Sep 17 00:00:00 2001 From: aho Date: Tue, 17 Jan 2023 12:29:54 +0100 Subject: [PATCH 2/2] Correct ttests for cross-validation ref when it is NULL --- R/CRPSS.R | 14 +++++++------- tests/testthat/test-CRPSS.R | 24 ++++++++++++------------ 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/R/CRPSS.R b/R/CRPSS.R index d31a264..cfd8fdf 100644 --- a/R/CRPSS.R +++ b/R/CRPSS.R @@ -223,7 +223,7 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', ## using climatology as reference forecast ## all the time steps are used as if they were members ## then, ref dimensions are [sdate, memb] - ## memb dimension has length(sdate)-1 due to cross-validation + ## memb dimension has length(sdate) - 1 due to cross-validation obs_time_len <- dim(obs)[time_dim] if (is.null(dat_dim)) { @@ -231,9 +231,9 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', ## Without cross-validation: ## ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len)) ## With cross-validation (excluding the value of that year to create ref for that year): - ref <- array(data = NA, dim = c(obs_time_len, obs_time_len-1)) - for (i in 1:obs_time_len){ - ref[i,] <- obs[-i] + ref <- array(data = NA, dim = c(obs_time_len, obs_time_len - 1)) + for (i in 1:obs_time_len) { + ref[i, ] <- obs[-i] } names(dim(ref)) <- c(time_dim, memb_dim) @@ -250,9 +250,9 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member', ## Without cross-validation: ## ref <- array(data = rep(obs[, i_obs], each = obs_time_len), dim = c(obs_time_len, obs_time_len)) ## With cross-validation (excluding the value of that year to create ref for that year): - ref <- array(data = NA, dim = c(obs_time_len, obs_time_len-1)) - for (i in 1:obs_time_len){ - ref[i,] <- obs[-i, i_obs] + ref <- array(data = NA, dim = c(obs_time_len, obs_time_len - 1)) + for (i in 1:obs_time_len) { + ref[i, ] <- obs[-i, i_obs] } names(dim(ref)) <- c(time_dim, memb_dim) diff --git a/tests/testthat/test-CRPSS.R b/tests/testthat/test-CRPSS.R index 6937edb..db0eecd 100644 --- a/tests/testthat/test-CRPSS.R +++ b/tests/testthat/test-CRPSS.R @@ -141,7 +141,7 @@ test_that("2. Output checks: dat1", { # ref = NULL expect_equal( as.vector(CRPSS(exp1, obs1)$crpss), - c(-0.1582765, -0.2390707), + c(0.061796021, -0.003647287), tolerance = 0.0001 ) expect_equal( @@ -150,12 +150,12 @@ test_that("2. Output checks: dat1", { ) expect_equal( as.vector(CRPSS(exp1, obs1, Fair = T)$crpss), - c(0.07650872, -0.09326681), + c(0.2612070, 0.1253865), tolerance = 0.0001 ) expect_equal( as.vector(CRPSS(exp1, obs1)$crpss), - c(-0.1582765, -0.2390707), + c(0.061796021, -0.003647287), tolerance = 0.0001 ) # ref = ref @@ -175,7 +175,7 @@ test_that("2. Output checks: dat1", { ) expect_equal( as.vector(CRPSS(exp1, obs1, ref1)$crpss), - c( 0.3491793, 0.3379610), + c(0.3491793, 0.3379610), tolerance = 0.0001 ) @@ -211,7 +211,7 @@ test_that("3. Output checks: dat2", { # ref = NULL expect_equal( as.vector(CRPSS(exp2, obs2)$crpss), - c(-0.8209236), + -0.4749481, tolerance = 0.0001 ) expect_equal( @@ -222,11 +222,11 @@ test_that("3. Output checks: dat2", { expect_equal( as.vector(CRPSS(exp2, obs2)$sign), - TRUE, + FALSE, ) expect_equal( as.vector(CRPSS(exp2, obs2, Fair = T)$crpss), - c(-0.468189), + c(-0.1745512), tolerance = 0.0001 ) expect_equal( @@ -265,12 +265,12 @@ test_that("4. Output checks: dat3", { ) expect_equal( mean(CRPSS(exp3, obs3, dat_dim = 'dataset')$crpss), - c(-0.7390546), + c(-0.4086342), tolerance = 0.0001 ) expect_equal( as.vector(CRPSS(exp3, obs3, dat_dim = 'dataset')$crpss), - c(-0.8209236, -0.6270744, -1.0829403, -0.6574485, -0.5970569, -0.6488837), + c(-0.4749481, -0.3179303, -0.6871816, -0.3425333, -0.2936161, -0.3355958), tolerance = 0.0001 ) expect_equal( @@ -279,12 +279,12 @@ test_that("4. Output checks: dat3", { ) expect_equal( mean(CRPSS(exp3, obs3, dat_dim = 'dataset', Fair = T)$crpss), - c(-0.5302703), + -0.2242162, tolerance = 0.0001 ) expect_equal( as.vector(CRPSS(exp3, obs3, dat_dim = 'dataset', Fair = T)$crpss), - c(-0.4681890, -0.3580685, -1.0116628, -0.3730576, -0.3965618, -0.5740819), + c(-0.17455119, -0.08645477, -0.60933021, -0.09844606, -0.11724947, -0.25926556), tolerance = 0.0001 ) # ref = ref3 @@ -310,7 +310,7 @@ test_that("5. Output checks: dat4", { ) expect_equal( as.vector(CRPSS(exp4, obs4, dat_dim = 'dataset', Fair = T)$crpss)[1:6], - c(-0.4681890, -0.3580685, -1.0116628, -0.3730576, -0.3965618, -0.5740819), + c(-0.17455119, -0.08645477, -0.60933021, -0.09844606, -0.11724947, -0.25926556), tolerance = 0.0001 ) # ref = ref3 -- GitLab