diff --git a/DESCRIPTION b/DESCRIPTION index 000b937711d1e789ebe74a93193b623af5a9fbc1..3e9a394a2917c44f1b4f2dcc11c51667c392016a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,19 +14,20 @@ Authors@R: c( person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"), person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"), person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"), - person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "ctb")) + person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut")) Description: Exploits dynamical seasonal forecasts in order to provide information relevant to stakeholders at the seasonal timescale. The package contains process-based methods for forecast calibration, bias correction, statistical and stochastic downscaling, optimal forecast combination and multivariate verification, as well as basic and advanced tools to obtain - tailored products. This package was developed in the context of the + tailored products. This package was developed in the context of the ERA4CS project MEDSCOPE and the H2020 S2S4E project. Doblas-Reyes et al. (2005) . Mishra et al. (2018) . Terzago et al. (2018) . Torralba et al. (2017) . D'Onofrio et al. (2014) . + Van Schaeybroeck et al. (2015) . Depends: R (>= 3.2.0), maps diff --git a/NAMESPACE b/NAMESPACE index 50c19d9c68bca03a3886d7e0c6ae9655d12b7ee3..f57c5ac4083998f071d9b2d6f96da1c9df30d00e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ export(RFSlope) export(RainFARM) export(as.s2dv_cube) export(s2dv_cube) +import(abind) import(ggplot2) import(multiApply) import(ncdf4) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index bd1b6d366089317ce796114d278d45faf5fa6f75..89821dca448b11bda9c4d5c0f832c83f9d14fb1f 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -31,7 +31,7 @@ #'str(a) #'@export CST_BiasCorrection <- function(exp, obs) { - if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) { + if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R index 5115f1e831c16ce688adb2a9d9c723dc81348b67..c63dbb837cd573cb3236ef7d406296f9586ea15c 100644 --- a/R/CST_Calibration.R +++ b/R/CST_Calibration.R @@ -1,55 +1,84 @@ -#'Forecast Calibration based on the ensemble inflation -#' -#'@author Verónica Torralba, \email{veronica.torralba@bsc.es} -#'@description This function applies a variance inflation technique described in Doblas-Reyes et al. (2005) in leave-one-out cross-validation. This bias adjustment method produces calibrated forecasts with equivalent mean and variance to that of the reference dataset, but at the same time preserve reliability. +#'Forecast Calibration #' +#'@author Verónica Torralba, \email{veronica.torralba@bsc.es} +#'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be} +#'@description Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). +#'@description Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. #'@references Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x +#'@references Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 +#'@references Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post‐processing using member‐by‐member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 #' #'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}. #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}. -#' +#'@param cal.method is the calibration method used, can be either \code{bias}, \code{evmos}, \code{mse_min} or \code{crps_min}. Default value is \code{mse_min}. +#'@param eval.method is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation. +#'@param multi.model is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result. #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data. #' #'@import s2dverification -#'@import multiApply +#'@import abind #' #'@seealso \code{\link{CST_Load}} #' -#'@examples #'# Example -#'# Load data using CST_Load or use the sample data provided: -#'library(zeallot) -#'c(exp, obs) %<-% areave_data -#'exp_calibrated <- CST_Calibration(exp = exp, obs = obs) -#'str(exp_calibrated) +#'# Creation of sample s2dverification objects. These are not complete +#'# s2dverification objects though. The Load function returns complete objects. +#'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) +#'dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) +#'obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) +#'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) +#'lon <- seq(0, 30, 5) +#'lat <- seq(0, 25, 5) +#'exp <- list(data = mod1, lat = lat, lon = lon) +#'obs <- list(data = obs1, lat = lat, lon = lon) +#'attr(exp, 'class') <- 's2dv_cube' +#'attr(obs, 'class') <- 's2dv_cube' +#'a <- CST_Calibration(exp = exp, obs = obs, cal.method = "mse_min", eval.method = "in-sample") +#'str(a) #'@export -CST_Calibration <- function(exp, obs) { - if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) { + + +CST_Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-one-out", multi.model = F) { + if (!inherits(exp, "s2dv_cube") || !inherits(obs, "s2dv_cube")) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } - if (dim(obs$data)['member'] != 1) { + if (dim(obs$data)["member"] != 1) { stop("The length of the dimension 'member' in the component 'data' ", "of the parameter 'obs' must be equal to 1.") } - dimnames <- names(dim(exp$data)) - Calibrated <- Calibration(exp = exp$data, obs = obs$data) - pos <- match(dimnames, names(dim(Calibrated))) - Calibrated <- aperm(Calibrated, pos) - names(dim(Calibrated)) <- dimnames - exp$data <- Calibrated + + if(!missing(multi.model) & !(cal.method == "mse_min")){ + warning(paste0("The multi.model parameter is ignored when using the calibration method ", cal.method)) + } + exp$data <- .calibration.wrap(exp = exp$data, obs = obs$data, + cal.method = cal.method, eval.method = eval.method, + multi.model = multi.model) + exp$Datasets <- c(exp$Datasets, obs$Datasets) exp$source_files <- c(exp$source_files, obs$source_files) return(exp) + } -Calibration <- function(exp, obs) { - if (!all(c('member', 'sdate') %in% names(dim(exp)))) { +.calibration.wrap <- function(exp, obs, cal.method, eval.method, multi.model) { + + dim.exp <- dim(exp) + amt.dims.exp <- length(dim.exp) + dim.obs <- dim(obs) + amt.dims.obs <- length(dim.obs) + dim.names.exp <- names(dim.exp) + dim.names.obs <- names(dim.obs) + + target.dim.names.exp <- c("member", "sdate") + target.dim.names.obs <- c("sdate") + + if (!all(target.dim.names.exp %in% dim.names.exp)) { stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.") } - if (!all(c('sdate') %in% names(dim(obs)))) { + if (!all(c("sdate") %in% dim.names.obs)) { stop("Parameter 'obs' must have the dimension 'sdate'.") } @@ -60,65 +89,221 @@ Calibration <- function(exp, obs) { if (any(is.na(obs))) { warning("Parameter 'obs' contains NA values.") } - - target_dims_obs <- 'sdate' - if ('member' %in% names(dim(obs))) { - target_dims_obs <- c('member', target_dims_obs) + + if (dim(obs)['member']!=1){ + stop("Parameter 'obs' must have a member dimension with length=1") } + + obs <- Subset(obs, along = "member", indices = 1, drop = "selected") + + calibrated <- Apply(data = list(mod = exp, obs = obs), + cal.method = cal.method, + eval.method = eval.method, + multi.model = multi.model, + target_dims = list(mod = target.dim.names.exp, obs = target.dim.names.obs), + fun = .cal)$output1 - Calibrated <- Apply(data = list(var_obs = obs, var_exp = exp), - target_dims = list(target_dims_obs, c('member', 'sdate')), - fun = .cal)$output1 + dexes <- match(names(dim(exp)), names(dim(calibrated))) + calibrated <- aperm(calibrated, dexes) + dimnames(calibrated) <- dimnames(exp)[dexes] + - return(Calibrated) + return(calibrated) } -.cal <- function(var_obs, var_exp) { - ntime <- dim(var_exp)[which(names(dim(var_exp)) == 'sdate')][] - nmembers <- dim(var_exp)[which(names(dim(var_exp)) == 'member')][] - - if (all(names(dim(var_exp)) != c('member','sdate'))) { - var_exp <- t(var_exp) - } - climObs <- NA * var_obs - climPred <- NA * var_obs - for (t in 1 : length(var_obs)) - { - climObs[t] <- mean(var_obs[ , -t]) - climPred[t] <- mean(var_exp[ , -t]) +.make.eval.train.dexes <- function(eval.method, amt.points){ + if(eval.method == "leave-one-out"){ + dexes.lst <- lapply(seq(1, amt.points), function(x) return(list(eval.dexes = x, train.dexes = seq(1, amt.points)[-x]))) + } else if (eval.method == "in-sample"){ + dexes.lst <- list(list(eval.dexes = seq(1, amt.points), train.dexes = seq(1, amt.points))) + } else { + stop(paste0("unknown sampling method: ",eval.method)) } - var_obs <- var_obs - climObs - - calibrated <- NA * var_exp - for (t in 1 : ntime) { - # defining forecast,hindcast and observation in cross-validation - fcst <- NA * var_exp[ , t] - hcst <- NA * var_exp[ , -t] - for (i in 1 : nmembers) { - fcst[i] <- var_exp[i, t] - climPred[t] - hcst[i, ] <- var_exp[i, -t]- climPred[t] - } - obs <- var_obs[-t] - #coefficients - em_fcst <- mean(fcst) - em_hcst <- apply(hcst, c(2), mean) - corr <- cor(em_hcst, obs) - sd_obs <- sd(obs) - sd_em_hcst <- sd(em_hcst) - - fcst_diff <- fcst - em_fcst - hcst_diff <- NA * hcst - for (n in 1 : nmembers) { - hcst_diff[n,] <- hcst[n,] - em_hcst - } - sd_hcst_diff <- sd(hcst_diff) + return(dexes.lst) +} + +.cal <- function(mod, obs, cal.method, eval.method, multi.model) { + + var.obs <- as.vector(obs) + var.fc <- mod + dims.fc <- dim(var.fc) + amt.mbr <- dims.fc[1] + amt.sdate <- dims.fc[2] + var.cor.fc <- NA * var.fc + + eval.train.dexeses <- .make.eval.train.dexes(eval.method, amt.points = amt.sdate) + amt.resamples <- length(eval.train.dexeses) + for (i.sample in seq(1, amt.resamples)) { + # defining training (tr) and evaluation (ev) subsets + eval.dexes <- eval.train.dexeses[[i.sample]]$eval.dexes + train.dexes <- eval.train.dexeses[[i.sample]]$train.dexes - a <- corr * (sd_obs / sd_em_hcst) - b <- (sd_obs / sd_hcst_diff) * sqrt(1 - (corr ^ 2)) + fc.ev <- var.fc[ , eval.dexes, drop = FALSE] + fc.tr <- var.fc[ , train.dexes] + obs.tr <- var.obs[train.dexes , drop = FALSE] - calibrated[, t] <- (a * em_fcst) + (b * fcst_diff) + climObs[t] + if(cal.method == "bias"){ + var.cor.fc[ , eval.dexes] <- fc.ev + mean(obs.tr, na.rm = TRUE) - mean(fc.tr, na.rm = TRUE) + } else if(cal.method == "evmos"){ + #calculate ensemble and observational characteristics + quant.obs.fc.tr <- .calc.obs.fc.quant(obs = obs.tr, fc = fc.tr) + #calculate value for regression parameters + init.par <- c(.calc.evmos.par(quant.obs.fc.tr)) + #correct evaluation subset + var.cor.fc[ , eval.dexes] <- .correct.evmos.fc(fc.ev , init.par) + } else if (cal.method == "mse_min"){ + #calculate ensemble and observational characteristics + quant.obs.fc.tr <- .calc.obs.fc.quant(obs = obs.tr, fc = fc.tr) + #calculate value for regression parameters + init.par <- .calc.mse.min.par(quant.obs.fc.tr, multi.model) + #correct evaluation subset + var.cor.fc[ , eval.dexes] <- .correct.mse.min.fc(fc.ev , init.par) + } else if (cal.method == "crps_min"){ + #calculate ensemble and observational characteristics + quant.obs.fc.tr <- .calc.obs.fc.quant.ext(obs = obs.tr, fc = fc.tr) + #calculate initial value for regression parameters + init.par <- c(.calc.mse.min.par(quant.obs.fc.tr), 0.001) + init.par[3] <- sqrt(init.par[3]) + #calculate regression parameters on training dataset + optim.tmp <- optim(par = init.par, fn = .calc.crps.opt, + quant.obs.fc = quant.obs.fc.tr, method = "Nelder-Mead") + + mbm.par <- optim.tmp$par + #correct evaluation subset + var.cor.fc[ , eval.dexes] <- .correct.crps.min.fc(fc.ev , mbm.par) + } else { + stop("unknown calibration method: ",cal.method) + } } - names(dim(calibrated)) <- c('member', 'sdate') - return(calibrated) + names(dim(var.cor.fc)) <- c("member", "sdate") + return(var.cor.fc) +} + +.calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations + amt.mbr <- dim(fc)[1] + obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr) + fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) + cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") + obs.av <- mean(obs, na.rm = TRUE) + obs.sd <- sd(obs, na.rm = TRUE) + return( + append( + .calc.fc.quant(fc = fc), + list( + obs.per.ens = obs.per.ens, + cor.obs.fc = cor.obs.fc, + obs.av = obs.av, + obs.sd = obs.sd + ) + ) + ) +} + +.calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations + amt.mbr <- dim(fc)[1] + obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr) + fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) + cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs") + obs.av <- mean(obs, na.rm = TRUE) + obs.sd <- sd(obs, na.rm = TRUE) + + return( + append( + .calc.fc.quant.ext(fc = fc), + list( + obs.per.ens = obs.per.ens, + cor.obs.fc = cor.obs.fc, + obs.av = obs.av, + obs.sd = obs.sd + ) + ) + ) +} + + +.calc.fc.quant <- function(fc){ #function to calculate different quantities of a series of ensemble forecasts + amt.mbr <- dim(fc)[1] + fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE) + fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE) + fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE) + fc.ens.av.per.ens <- InsertDim(var = fc.ens.av, posdim = 1, lendim = amt.mbr) + fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE) + fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE)) + fc.dev <- fc - fc.ens.av.per.ens + fc.dev.sd <- sd(fc.dev, na.rm = TRUE) + fc.av <- mean(fc, na.rm = TRUE) + fc.sd <- sd(fc, na.rm = TRUE) + return( + list( + fc.ens.av = fc.ens.av, + fc.ens.av.av = fc.ens.av.av, + fc.ens.av.sd = fc.ens.av.sd, + fc.ens.av.per.ens = fc.ens.av.per.ens, + fc.ens.sd = fc.ens.sd, + fc.ens.var.av.sqrt = fc.ens.var.av.sqrt, + fc.dev = fc.dev, + fc.dev.sd = fc.dev.sd, + fc.av = fc.av, + fc.sd = fc.sd + ) + ) +} + +.calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts + + amt.mbr <- dim(fc)[1] + repmat1.tmp <- InsertDim(var = fc, posdim = 1, lendim = amt.mbr) + repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3)) + spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE) + spr.abs.per.ens <- InsertDim(var = spr.abs, posdim = 1, lendim = amt.mbr) + + return( + append(.calc.fc.quant(fc), + list(spr.abs = spr.abs, spr.abs.per.ens = spr.abs.per.ens)) + ) +} + +#Below are the core or elementary functions to calculate the regression parameters for the different methods +.calc.mse.min.par <- function(quant.obs.fc, multi.model = F){ + par.out <- rep(NA, 3) + + if(multi.model){ + par.out[3] <- with(quant.obs.fc, obs.sd * sqrt(1. - cor.obs.fc^2) / fc.ens.var.av.sqrt) + } else { + par.out[3] <- with(quant.obs.fc, obs.sd * sqrt(1. - cor.obs.fc^2) / fc.dev.sd) + } + par.out[2] <- with(quant.obs.fc, cor.obs.fc * obs.sd / fc.ens.av.sd) + par.out[1] <- with(quant.obs.fc, obs.av - par.out[2] * fc.ens.av.av, na.rm = TRUE) + + return(par.out) +} +.calc.evmos.par <- function(quant.obs.fc){ + par.out <- rep(NA, 2) + par.out[2] <- with(quant.obs.fc, obs.sd / fc.sd) + par.out[1] <- with(quant.obs.fc, obs.av - par.out[2] * fc.ens.av.av, na.rm = TRUE) + return(par.out) +} +#Below are the core or elementary functions to calculate the functions necessary for the minimization of crps +.calc.crps.opt <- function(par, quant.obs.fc){ + return( + with(quant.obs.fc, + mean(abs(obs.per.ens - (par[1] + par[2] * fc.ens.av.per.ens + + ((par[3])^2 + par[4] / spr.abs.per.ens) * fc.dev)), na.rm = TRUE) - + mean(abs((par[3])^2 * spr.abs + par[4]) / 2., na.rm = TRUE) + ) + ) } +#Below are the core or elementary functions to correct the evaluation set based on the regression parameters +.correct.evmos.fc <- function(fc, par){ + quant.fc.mp <- .calc.fc.quant(fc = fc) + return(with(quant.fc.mp, par[1] + par[2] * fc)) +} +.correct.mse.min.fc <- function(fc, par){ + quant.fc.mp <- .calc.fc.quant(fc = fc) + return(with(quant.fc.mp, par[1] + par[2] * fc.ens.av.per.ens + fc.dev * par[3])) +} +.correct.crps.min.fc <- function(fc, par){ + quant.fc.mp <- .calc.fc.quant.ext(fc = fc) + return(with(quant.fc.mp, par[1] + par[2] * fc.ens.av.per.ens + fc.dev * abs((par[3])^2 + par[4] / spr.abs))) +} diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index 485199eae252039266185188df10ee21cf1af52b..75207e2e6e4db8e934dc8a51ba589e90fd64d81d 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -42,3 +42,4 @@ Verónica Torralba, \email{veronica.torralba@bsc.es} Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS) } \encoding{UTF-8} + diff --git a/man/CST_Calibration.Rd b/man/CST_Calibration.Rd index c841907683cdafaa996291ad634b20a7ad0464e8..0e6a6a5423b292e0a62a4f360db81c6583681c11 100644 --- a/man/CST_Calibration.Rd +++ b/man/CST_Calibration.Rd @@ -2,36 +2,59 @@ % Please edit documentation in R/CST_Calibration.R \name{CST_Calibration} \alias{CST_Calibration} -\title{Forecast Calibration based on the ensemble inflation} +\title{Forecast Calibration} \usage{ -CST_Calibration(exp, obs) +CST_Calibration(exp, obs, cal.method = "mse_min", + eval.method = "leave-one-out", multi.model = F) } \arguments{ \item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.} \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.} + +\item{cal.method}{is the calibration method used, can be either \code{bias}, \code{evmos}, \code{mse_min} or \code{crps_min}. Default value is \code{mse_min}.} + +\item{eval.method}{is the sampling method used, can be either \code{in-sample} or \code{leave-one-out}. Default value is the \code{leave-one-out} cross validation.} + +\item{multi.model}{is a boolean that is used only for the \code{mse_min} method. If multi-model ensembles or ensembles of different sizes are used, it must be set to \code{TRUE}. By default it is \code{FALSE}. Differences between the two approaches are generally small but may become large when using small ensemble sizes. Using multi.model when the calibration method is \code{bias}, \code{evmos} or \code{crps_min} will not affect the result.} } \value{ an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data. } \description{ -This function applies a variance inflation technique described in Doblas-Reyes et al. (2005) in leave-one-out cross-validation. This bias adjustment method produces calibrated forecasts with equivalent mean and variance to that of the reference dataset, but at the same time preserve reliability. -} -\examples{ -# Example -# Load data using CST_Load or use the sample data provided: -library(zeallot) -c(exp, obs) \%<-\% areave_data -exp_calibrated <- CST_Calibration(exp = exp, obs = obs) -str(exp_calibrated) +Four types of member-by-member bias correction can be performed. The \code{bias} method corrects the bias only, the \code{evmos} method applies a variance inflation technique to ensure the correction of the bias and the correspondence of variance between forecast and observation (Van Schaeybroeck and Vannitsem, 2011). The ensemble calibration methods \code{"mse_min"} and \code{"crps_min"} correct the bias, the overall forecast variance and the ensemble spread as described in Doblas-Reyes et al. (2005) and Van Schaeybroeck and Vannitsem (2015), respectively. While the \code{"mse_min"} method minimizes a constrained mean-squared error using three parameters, the \code{"crps_min"} method features four parameters and minimizes the Continuous Ranked Probability Score (CRPS). + +Both in-sample or our out-of-sample (leave-one-out cross validation) calibration are possible. } \author{ Verónica Torralba, \email{veronica.torralba@bsc.es} + +Bert Van Schaeybroeck, \email{bertvs@meteo.be} } \references{ Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting-II calibration and combination. Tellus A. 2005;57:234-252. doi:10.1111/j.1600-0870.2005.00104.x + +Van Schaeybroeck, B., & Vannitsem, S. (2011). Post-processing through linear regression. Nonlinear Processes in Geophysics, 18(2), 147. doi:10.5194/npg-18-147-2011 + +Van Schaeybroeck, B., & Vannitsem, S. (2015). Ensemble post‐processing using member‐by‐member approaches: theoretical aspects. Quarterly Journal of the Royal Meteorological Society, 141(688), 807-818. doi:10.1002/qj.2397 } \seealso{ \code{\link{CST_Load}} + +# Example +# Creation of sample s2dverification objects. These are not complete +# s2dverification objects though. The Load function returns complete objects. +mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7) +dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7) +obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7) +dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7) +lon <- seq(0, 30, 5) +lat <- seq(0, 25, 5) +exp <- list(data = mod1, lat = lat, lon = lon) +obs <- list(data = obs1, lat = lat, lon = lon) +attr(exp, 'class') <- 's2dv_cube' +attr(obs, 'class') <- 's2dv_cube' +a <- CST_Calibration(exp = exp, obs = obs, cal.method = "mse_min", eval.method = "in-sample") +str(a) } \encoding{UTF-8} diff --git a/man/manual.pdf b/man/manual.pdf new file mode 100644 index 0000000000000000000000000000000000000000..4f5512918fead76fa41a65f8aed983f1f9a74e1c Binary files /dev/null and b/man/manual.pdf differ diff --git a/tests/testthat/test-CST_Calibration.R b/tests/testthat/test-CST_Calibration.R index 917b6ac5bd31c565709881b7ee7e1918c29d9d9f..a4491719ce63e224965a60000926c18799bffe7f 100644 --- a/tests/testthat/test-CST_Calibration.R +++ b/tests/testthat/test-CST_Calibration.R @@ -8,11 +8,11 @@ test_that("Sanity checks", { CST_Calibration(obs = 1), c("argument \"exp\" is missing, with no default") ) - library(zeallot) + library(zeallot) c(exp, obs) %<-% lonlat_data cal <- CST_Calibration(exp = exp, obs = obs) expect_equal(length(cal), 9) - expect_equal(dim(cal$data), dim(exp$data)) + expect_equal(as.numeric(dim(cal$data)), as.numeric(dim(exp$data))) expect_equal(cal$lat, exp$lat) expect_equal(cal$lat, obs$lat) expect_equal(cal$lon, exp$lon) diff --git a/tests/testthat/test-PlotForecastPDF.R b/tests/testthat/test-PlotForecastPDF.R index 1515a13fe084134aa5f44a13b6cde4624be8c77c..dfc35c4d7a82308a8a4dde8a618ee1f8664ec2f1 100644 --- a/tests/testthat/test-PlotForecastPDF.R +++ b/tests/testthat/test-PlotForecastPDF.R @@ -20,4 +20,3 @@ test_that("Sanity checks", { PlotForecastPDF(fcst = fcsts2, tercile.limits = c(-0.5, 0.5), extreme.limits = NA), "Parameter 'extreme.limits' should be an array with two limits for delimiting extreme categories") }) -