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; however,
#'it gives no information about the typical magnitude of individual forecast errors.
#'
#'@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 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 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 the dimensions of
#''exp' except 'time_dim' and 'memb_dim' dimensions.
#'
#'@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, year = 50))
#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50))
#'bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1)
#'
#'@import multiApply
#'@export
Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) {
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
# 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.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.")
}
## dat_dim
# 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)))
name_exp <- name_exp[-which(name_exp == memb_dim)]
if (memb_dim %in% name_obs) {
name_obs <- name_obs[-which(name_obs == 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.")
}
## 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
bias <- multiApply::Apply(data = list(exp, obs),
target_dims = time_dim,
fun = .Bias,
ncores = ncores)$output1
## Return the mean bias
bias <- MeanDims(bias, time_dim, na.rm = na.rm)
return(bias)
}
.Bias <- function(exp, obs) {
bias <- exp - obs
return(bias)
}