Newer
Older
#'Compute the Mean Bias
#'
#'The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference
#'between the ensemble mean forecast and the observations. It is a deterministic
#'metric. Positive values indicate that the forecasts are on average too high
#'and negative values indicate that the forecasts are on average too low.
#'It also allows to compute the Absolute Mean Bias.
#'
#'@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' and
#' 'dat_dim'.
#'@param time_dim A character string indicating the name of the time dimension.
#' The default value is 'sdate'.
#'@param dat_dim A character string indicating the name of dataset dimension.
#' The length of this dimension can be different between 'exp' and 'obs'.
#' The default value is NULL.
#'@param memb_dim A character string indicating the name of the member dimension
#' to compute the ensemble mean; it should be set to NULL if the parameter
#' 'exp' is already the ensemble mean. The default value is NULL.
#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or
#' kept (FALSE) for computation. The default value is FALSE.
#'@param absolute A logical value indicating whether to compute the absolute
#' bias. The default value is FALSE.
#'@param time_dim A logical value indicating whether to compute the temporal
#' mean of the bias. The default value is TRUE.
#'@param ncores An integer indicating the number of cores to use for parallel
#' computation. The default value is NULL.
#'
#'@return
#'A numerical array of bias with dimensions nexp, nobs and the rest dimensions
#'of 'exp' except 'time_dim' (if time_mean = T). nexp is the number of
#'experiment (i.e., 'dat_dim' in exp), and nobs is the number of observation
#'(i.e., 'dat_dim' in obs). If dat_dim is NULL, nexp and nobs are omitted.
#'
#'@references
#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
#'
#'@examples
#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
#'bias <- Bias(exp = exp, obs = obs, memb_dim = 'member')
Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE,
absolute = FALSE, time_mean = TRUE, 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.null(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 (memb_dim %in% names(dim(obs))) {
if (identical(as.numeric(dim(obs)[memb_dim]), 1)) {
obs <- ClimProjDiags::Subset(x = obs, along = memb_dim, indices = 1, drop = 'selected')
} else {
stop("Not implemented for observations with members ('obs' can have 'memb_dim', ",
"but it should be of length = 1).")
}
}
if (!is.null(dat_dim)) {
if (!is.character(dat_dim) | length(dat_dim) > 1) {
stop("Parameter 'dat_dim' must be a character string.")
}
if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.",
" Set it as NULL if there is no dataset dimension.")
}
}
## exp and obs (2)
name_exp <- sort(names(dim(exp)))
name_obs <- sort(names(dim(obs)))
if (!is.null(memb_dim)) {
name_exp <- name_exp[-which(name_exp == memb_dim)]
if (!is.null(dat_dim)) {
name_exp <- name_exp[-which(name_exp == dat_dim)]
name_obs <- name_obs[-which(name_obs == dat_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 except 'memb_dim' and 'dat_dim'."))
}
## na.rm
if (!is.logical(na.rm) | length(na.rm) > 1) {
stop("Parameter 'na.rm' must be one logical value.")
}
## absolute
if (!is.logical(absolute) | length(absolute) > 1) {
stop("Parameter 'absolute' must be one logical value.")
}
## time_mean
if (!is.logical(time_mean) | length(time_mean) > 1) {
stop("Parameter 'time_mean' must be one logical value.")
}
## 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.")
}
}
## Ensemble mean
if (!is.null(memb_dim)) {
exp <- MeanDims(exp, memb_dim, na.rm = na.rm)
}
## (Mean) Bias
target_dims = c(time_dim, dat_dim),
fun = .Bias,
dat_dim = dat_dim,
na.rm = na.rm,
absolute = absolute,
time_mean = time_mean,
.Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, absolute = FALSE, time_mean = TRUE) {
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
if (is.null(dat_dim)) {
bias <- exp - obs
if (isTRUE(absolute)){
bias <- abs(bias)
}
if (isTRUE(time_mean)){
bias <- mean(bias, na.rm = na.rm)
}
} else {
nexp <- as.numeric(dim(exp)[dat_dim])
nobs <- as.numeric(dim(obs)[dat_dim])
bias <- array(dim = c(dim(exp)[1], nexp = nexp, nobs = nobs))
for (i in 1:nexp) {
for (j in 1:nobs) {
exp_data <- exp[ , i]
obs_data <- obs[ , j]
if (is.null(dim(exp_data))) dim(exp_data) <- c(dim(exp)[1])
if (is.null(dim(obs_data))) dim(obs_data) <- c(dim(obs)[1])
bias[ , i, j] <- exp_data - obs_data
}
}
if (isTRUE(absolute)){
bias <- abs(bias)
}
if (isTRUE(time_mean)){
bias <- MeanDims(bias, time_dim, na.rm = na.rm)