R/calculate_threshold_probs.R

In DUToolkit: Visualizing and Quantifying Decision Uncertainty

#' Calculate threshold exceedance probabilities
#'
#' @description For each policy alternative, this function calculates the
#' probability that the peak (or minimum) value exceeds (or is below)
#' a specified threshold(s) using a Riemann sum approach.
#'
#' @param max_min_values_list A list generated by [get_max_min_values()]
#' @param Dp A vector of threshold values to calculate the probabilities for.
#' @param Dt_max A logical value indicating whether the decision threshold
#' is a maximum (`TRUE`) or a minimum (`FALSE`). The default is `TRUE`.
#'
#' @return  A list of vectors of threshold values and corresponding
#' probabilities.
#' @export
#'
#' @examples
#' tmin <- "2021-01-01"
#' tmax <- "2021-04-10"
#' D <- 750
#' Dp <- c(750, 1000, 2000)
#'
#' peak_values_list <- get_max_min_values(
#'   psa_data,
#'   tmin = tmin,
#'   tmax = tmax,
#'   Dt_max = TRUE
#' )
#'
#' calculate_threshold_probs(
#'   peak_values_list,
#'   Dp = Dp,
#'   Dt_max = TRUE
#' )
calculate_threshold_probs <- function(
    max_min_values_list,
    Dp,
    Dt_max = TRUE) {
  if (inherits(max_min_values_list, "list")) {
    values <- lapply(max_min_values_list, calculate_threshold_probs_1, Dp, Dt_max)
  } else if (inherits(max_min_values_list, "data.frame")) {
    values <- calculate_threshold_probs_1(max_min_values_list, Dp, Dt_max)
  } else {
    rlang::abort("The first argument is not a data.frame or list of data.frames",
      class = "data_type"
    )
  }
  return(values)
}

#' Calculates the probability that the peak (or minimum) value exceeds (or is below)
#' a specified threshold for a single set of simulations
#'
#' @inheritParams calculate_threshold_probs
#' @noRd
#' @return A vector of threshold values and corresponding probabilities
calculate_threshold_probs_1 <- function(
    max_min_values_list,
    Dp,
    Dt_max = TRUE) {
  d <- stats::density(max_min_values_list$outcome)
  xx <- d$x # 512 evenly spaced points on [min(x) - 3 * d$bw, max(x) + 3 * d$bw]
  dx <- xx[2] - xx[1] # spacing between x
  yy <- d$y # 512 density values for xx
  C <- sum(yy) * dx # Riemann Sum normalizing constant
  v_p <- c()
  if (Dt_max == TRUE) { # if threshold is a maximum
    for (j in 1:length(Dp)) {
      p.unscaled <- sum(yy[xx > Dp[j]]) * dx
      p.scaled <- round(p.unscaled / C, 4)
      v_p <- append(v_p, p.scaled)
    }
  } else {
    for (j in 1:length(Dp)) {
      p.unscaled <- sum(yy[xx < Dp[j]]) * dx
      p.scaled <- round(p.unscaled / C, 4)
      v_p <- append(v_p, p.scaled)
    }
  }
  names(v_p) <- Dp
  return(v_p)
}