R/misc.R In HuraultMisc: Guillem Hurault Functions' Library

Documented in compute_resolutioncompute_RPS

# Compute resolution ------------------------------------------------------

#' Compute resolution of forecasts, normalised by the uncertainty
#'
#' The resolution is computed as the mean squared distance to a base rate (reference forecast) and
#' is then normalised by the uncertainty (maximum resolution).
#' This means the output is between 0 and 1, 1 corresponding to the maximum resolution.
#'
#' @param f Vector of forecasts
#' @param p0 Vector of base rate.
#' In the case rate is usually the prevalence of a uniform forecast (e.g. 1 / number of categories)
#' but can depend on the observation (hence the vector).
#'
#' @return Vector of resolution values
#' @export
#'
#' @examples
#' compute_resolution(seq(0, 1, .1), 0.5)
compute_resolution <- function(f, p0) {

stopifnot(min(f) >= 0,
max(f) <= 1,
min(p0) >= 0,
max(p0) <= 1)

reso <- mean((f - p0)^2)
uncertainty <- mean(p0 * (1 - p0))
return(reso / uncertainty)
}

# Compute RPS -------------------------------------------------------------

#' Compute RPS for a single forecast
#'
#' @param forecast Vector of length N (forecast).
#' @param outcome Index of the true outcome (between 1 and N).
#'
#' @return RPS (numeric scalar)
#' @export
#'
#' @examples
#' compute_RPS(c(.2, .5, .3), 2)
compute_RPS <- function(forecast, outcome) {

stopifnot(is.vector(forecast, mode = "numeric"),
length(forecast) > 1)

if (any(is.na(c(forecast, outcome)))) {
RPS <- NA
} else {
stopifnot(all(between(forecast, 0, 1)),
round(sum(forecast), 2) == 1,
outcome %in% 1:length(forecast))
dummy_outcome <- 0 * forecast
dummy_outcome[outcome] <- 1
RPS <- sum((cumsum(forecast) - cumsum(dummy_outcome))^2) / (length(forecast) - 1)
}
return(RPS)
}

Try the HuraultMisc package in your browser

Any scripts or data that you put into this service are public.

HuraultMisc documentation built on Sept. 6, 2021, 9:09 a.m.