R/utils-numeric.R

In MultiFrailty: Shared Frailty Regression Models with Inverse Gaussian, Generalized Lindley, and Gamma Frailty Distributions

#' Numerically stable log(1 - exp(-x))
#'
#' Accurately evaluates \code{log(1 - exp(-x))} for any \code{x >= 0}
#' avoiding catastrophic loss of precision for small or large x.
#'
#' @param x Numeric vector of non-negative values.
#' @return Numeric vector of \code{log(1 - exp(-x))}.
#' @keywords internal
log1mexp <- function(x) {
  res <- numeric(length(x))
  idx_na <- is.na(x)
  idx_le0 <- !idx_na & (x <= 0)
  idx_high <- !idx_na & (x > log(2))
  idx_low <- !idx_na & (x > 0 & x <= log(2))

  res[idx_na] <- NA_real_
  res[idx_le0] <- -Inf
  res[idx_high] <- log1p(-exp(-x[idx_high]))
  res[idx_low] <- log(-expm1(-x[idx_low]))
  res
}

#' Numerically stable log(1 + exp(x))
#'
#' Accurately evaluates \code{log(1 + exp(x))} for all real x.
#'
#' @param x Numeric vector.
#' @return Numeric vector of \code{log(1 + exp(x))}.
#' @keywords internal
log1pexp <- function(x) {
  res <- numeric(length(x))
  idx_na <- is.na(x)
  idx_low <- !idx_na & (x <= -37)
  idx_mid <- !idx_na & (x > -37 & x <= 18)
  idx_high <- !idx_na & (x > 18 & x <= 33.3)
  idx_very_high <- !idx_na & (x > 33.3)

  res[idx_na] <- NA_real_
  res[idx_low] <- exp(x[idx_low])
  res[idx_mid] <- log1p(exp(x[idx_mid]))
  res[idx_high] <- x[idx_high] + log1p(exp(-x[idx_high]))
  res[idx_very_high] <- x[idx_very_high]
  res
}

#' Log-sum-exp helper for two log-scale components
#'
#' Computes \code{log(exp(a) + exp(b))} in log-space.
#'
#' @param a Numeric vector of log-scale component 1.
#' @param b Numeric vector of log-scale component 2.
#' @return Numeric vector of \code{log(exp(a) + exp(b))}.
#' @keywords internal
log_sum_exp <- function(a, b) {
  m <- pmax(a, b)
  m + log(exp(a - m) + exp(b - m))
}

#' Safe exponential function with overflow capping
#'
#' @param x Numeric vector.
#' @return Numeric vector capped to prevent overflow/underflow.
#' @keywords internal
safe_exp <- function(x) {
  exp(pmin(pmax(x, -700), 700))
}

#' Safe power evaluation
#'
#' @param base Positive numeric base.
#' @param exponent Numeric exponent.
#' @return Numeric result of base^exponent.
#' @keywords internal
safe_pow <- function(base, exponent) {
  exp(exponent * log(pmax(base, 1e-300)))
}