R/igp_distributions.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Inverse Gaussian Distribution Functions
#'
#' Density, distribution function, quantile function, and random generation
#' for the Inverse Gaussian (Wald) distribution with mean \code{mu} and shape/scale
#' parameter \code{lambda}.
#'
#' @param x,q Numeric vector of quantiles. Must be positive.
#' @param p Numeric vector of probabilities (between 0 and 1).
#' @param n Number of observations to generate. Must be a single positive integer.
#' @param mu Mean parameter. Must be strictly positive. Default is \code{1}.
#' @param lambda Shape parameter. Must be strictly positive. Default is \code{1}.
#' @param log,log.p Logical; if \code{TRUE}, probabilities/densities are returned on the log scale. Default is \code{FALSE}.
#' @param lower.tail Logical; if \code{TRUE} (default), probabilities are \eqn{P[X \le x]}, otherwise \eqn{P[X > x]}.
#'
#' @details
#' The probability density function of the Inverse Gaussian distribution is:
#' \deqn{f(x; \mu, \lambda) = \sqrt{\frac{\lambda}{2\pi x^3}} \exp\left(-\frac{\lambda (x - \mu)^2}{2\mu^2 x}\right), \quad x > 0}
#' The cumulative distribution function is evaluated using:
#' \deqn{F(x; \mu, \lambda) = \Phi\left(\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu} - 1\right)\right) + \exp\left(\frac{2\lambda}{\mu}\right) \Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu} + 1\right)\right)}
#' where \eqn{\Phi(\cdot)} is the standard normal cumulative distribution function.
#'
#' @return
#' \code{d_ig} returns the density, \code{p_ig} returns the distribution function,
#' \code{q_ig} returns the quantile function, and \code{r_ig} returns a vector of
#' random variates of length \code{n}.
#'
#' @references
#' Wasan, M. T. (1968). On an inverse Gaussian process. \emph{Scandinavian Actuarial Journal}, 1968(1-2), 69-96.
#'
#' @seealso \code{\link{igp_fit}}, \code{\link{sim_igp}}
#'
#' @examples
#' # Density and CDF
#' d_ig(1.5, mu = 2, lambda = 5)
#' p_ig(1.5, mu = 2, lambda = 5)
#'
#' # Quantiles and random generation
#' q_ig(0.5, mu = 2, lambda = 5)
#' set.seed(123)
#' r_ig(5, mu = 2, lambda = 5)
#'
#' @export
d_ig <- function(x, mu = 1, lambda = 1, log = FALSE) {
  if (length(x) == 0) return(numeric(0))
  idx <- (!is.na(x) & !is.na(mu) & !is.na(lambda) & x > 0 & mu > 0 & lambda > 0)
  out <- if (log) rep(-Inf, length(x)) else numeric(length(x))
  if (any(idx)) {
    xi <- x[idx]
    mui <- if (length(mu) == 1) mu else mu[idx]
    li <- if (length(lambda) == 1) lambda else lambda[idx]
    
    log_dens <- 0.5 * log(li) - 0.5 * log(2 * pi) - 1.5 * log(xi) - (li * (xi - mui)^2) / (2 * mui^2 * xi)
    out[idx] <- if (log) log_dens else exp(log_dens)
  }
  out
}

#' @rdname d_ig
#' @export
p_ig <- function(q, mu = 1, lambda = 1, lower.tail = TRUE, log.p = FALSE) {
  if (length(q) == 0) return(numeric(0))
  idx <- (!is.na(q) & !is.na(mu) & !is.na(lambda) & q > 0 & mu > 0 & lambda > 0)
  out <- numeric(length(q))
  if (any(idx)) {
    xi <- q[idx]
    mui <- if (length(mu) == 1) mu else mu[idx]
    li <- if (length(lambda) == 1) lambda else lambda[idx]
    
    log_factor <- 2 * li / mui
    z1 <- -sqrt(li / xi) * (1 + xi / mui)
    log_p1 <- pnorm(z1, log.p = TRUE)
    log_tot <- log_factor + log_p1
    term1 <- ifelse(log_tot > -700, exp(log_tot), 0)
    
    z2 <- sqrt(li / xi) * (xi / mui - 1)
    term2 <- pnorm(z2)
    
    prob <- pmin(pmax(term1 + term2, 0), 1)
    out[idx] <- prob
  }
  
  if (!lower.tail) {
    out <- 1 - out
  }
  if (log.p) {
    out <- log(pmax(out, 1e-300))
  }
  out
}

#' @rdname d_ig
#' @export
q_ig <- function(p, mu = 1, lambda = 1, lower.tail = TRUE, log.p = FALSE) {
  if (log.p) p <- exp(p)
  if (!lower.tail) p <- 1 - p
  
  out <- numeric(length(p))
  for (i in seq_along(p)) {
    pi <- p[i]
    if (is.na(pi) || pi < 0 || pi > 1) {
      out[i] <- NA
      next
    }
    if (pi == 0) {
      out[i] <- 0
      next
    }
    if (pi == 1) {
      out[i] <- Inf
      next
    }
    
    mui <- if (length(mu) == 1) mu else mu[i]
    li <- if (length(lambda) == 1) lambda else lambda[i]
    
    f_root <- function(x) p_ig(x, mu = mui, lambda = li) - pi
    lower_b <- 1e-8
    upper_b <- max(mui * 50, 100)
    while (f_root(upper_b) < 0) {
      upper_b <- upper_b * 2
    }
    res <- tryCatch(
      uniroot(f_root, lower = lower_b, upper = upper_b, tol = 1e-9)$root,
      error = function(e) NA
    )
    out[i] <- res
  }
  out
}

#' @rdname d_ig
#' @export
r_ig <- function(n, mu = 1, lambda = 1) {
  if (n <= 0) return(numeric(0))
  v <- rnorm(n)
  y <- v^2
  x <- mu + (mu^2 * y) / (2 * lambda) - (mu / (2 * lambda)) * sqrt(4 * mu * lambda * y + mu^2 * y^2)
  u <- runif(n)
  res <- ifelse(u <= mu / (mu + x), x, mu^2 / x)
  res
}

# Internal increment density
.d_igp_inc <- function(y, dg, eta) {
  idx <- (!is.na(y) & !is.na(dg) & !is.na(eta) & y > 0 & dg > 0 & eta > 0)
  out <- numeric(length(y))
  if (any(idx)) {
    yi <- y[idx]
    dgi <- if (length(dg) == 1) dg else dg[idx]
    etai <- if (length(eta) == 1) eta else eta[idx]
    out[idx] <- sqrt(etai / (2 * pi * yi^3)) * dgi * exp(-etai * (yi - dgi)^2 / (2 * yi))
  }
  out
}

# Internal increment CDF with precision buffer to avoid zero survival
.p_igp_inc <- function(y, dg, eta) {
  idx <- (!is.na(y) & !is.na(dg) & !is.na(eta) & y > 0 & dg > 0 & eta > 0)
  out <- numeric(length(y))
  if (any(idx)) {
    yi <- y[idx]
    dgi <- if (length(dg) == 1) dg else dg[idx]
    etai <- if (length(eta) == 1) eta else eta[idx]
    
    log_factor <- 2 * etai * dgi
    z1 <- -sqrt(etai / yi) * (yi + dgi)
    log_p1 <- pnorm(z1, log.p = TRUE)
    log_tot <- log_factor + log_p1
    term1 <- ifelse(log_tot > -700, exp(log_tot), 0)
    
    z2 <- sqrt(etai / yi) * (yi - dgi)
    term2 <- pnorm(z2)
    
    out[idx] <- pmin(pmax(term1 + term2, 0), 1 - 1e-15)
  }
  out
}

# Internal increment reliability
.r_igp_inc <- function(y, dg, eta) {
  pmax(1 - .p_igp_inc(y, dg, eta), 1e-300)
}

# Internal increment hazard
.h_igp_inc <- function(y, dg, eta) {
  d <- .d_igp_inc(y, dg, eta)
  r <- .r_igp_inc(y, dg, eta)
  d / r
}

# Internal increment cumulative hazard
.H_igp_inc <- function(y, dg, eta) {
  -log(.r_igp_inc(y, dg, eta))
}