R/igp_simulation.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Simulation of Inverse Gaussian Process Degradation Paths
#'
#' Simulates longitudinal degradation measurements for multiple units under
#' classical IGP, IGP-Gamma frailty, or IGP-IG frailty models.
#'
#' @param n Integer specifying the number of experimental units to simulate.
#' @param times Numeric vector of inspection times (e.g. \code{seq(0, 4, by = 0.25)}). Must start at 0.
#' @param theta Mean parameter \eqn{\theta > 0} (or numeric vector for non-linear mean functions).
#' @param eta Precision/scale parameter \eqn{\eta > 0}.
#' @param xi Frailty variance parameter \eqn{\xi > 0}. Ignored if \code{frailty = "none"}.
#' @param frailty Frailty specification: \code{"none"} (default), \code{"gamma"}, or \code{"ig"}.
#' @param mean_fun Mean degradation function \code{"linear"} (default), \code{"power"}, \code{"exponential"}, or a custom function.
#' @param seed Optional integer random seed for reproducibility.
#'
#' @return A data frame containing simulated degradation paths:
#' \describe{
#'   \item{unit}{Integer unit identifier (1 to \code{n}).}
#'   \item{t}{Inspection time.}
#'   \item{increment}{Simulated degradation increment \eqn{\Delta D(t)}.}
#'   \item{degradation}{Simulated cumulative degradation \eqn{D(t)}.}
#'   \item{frailty_z}{Realized individual frailty multiplier \eqn{z_i} for unit \eqn{i}.}
#' }
#'
#' @references
#' Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021).
#' Inverse Gaussian process model with frailty term in reliability analysis. \emph{Quality and Reliability Engineering International},
#' 37(2), 763-784. \doi{10.1002/qre.2762}.
#'
#' @seealso \code{\link{igp_fit}}, \code{\link{lifetime_dist}}
#'
#' @examples
#' set.seed(42)
#' sim_data <- sim_igp(n = 5, times = seq(0, 2, by = 0.5), theta = 1.5,
#'                     eta = 10, xi = 0.3, frailty = "gamma")
#' head(sim_data)
#'
#' @export
sim_igp <- function(n = 10, 
                    times = seq(0, 4, by = 0.25), 
                    theta = 2, 
                    eta = 15, 
                    xi = 0.2, 
                    frailty = c("none", "gamma", "ig"), 
                    mean_fun = "linear", 
                    seed = NULL) {
  frailty <- match.arg(frailty)
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  if (times[1] != 0) {
    times <- c(0, times)
  }
  times <- sort(unique(times))
  dt_vec <- diff(times)
  m <- length(dt_vec)
  
  # Mean function
  if (is.character(mean_fun)) {
    if (mean_fun == "linear") {
      g_fun <- function(t) theta[1] * t
    } else if (mean_fun == "power") {
      g_fun <- function(t) theta[1] * (t^theta[2])
    } else if (mean_fun == "exponential") {
      g_fun <- function(t) exp(theta[1] * t) - 1
    } else {
      stop("Unknown mean function.")
    }
  } else if (is.function(mean_fun)) {
    g_fun <- function(t) mean_fun(t, theta)
  } else {
    stop("Argument 'mean_fun' must be character or function.")
  }
  
  dg_vec <- g_fun(times[-1]) - g_fun(times[-length(times)])
  
  # Generate frailties
  if (frailty == "none") {
    z_vec <- rep(1, n)
  } else if (frailty == "gamma") {
    # Gamma(shape = 1/xi, scale = xi)
    z_vec <- rgamma(n, shape = 1 / xi, scale = xi)
  } else if (frailty == "ig") {
    # IG(mu = 1, lambda = 1/xi)
    z_vec <- r_ig(n, mu = 1, lambda = 1 / xi)
  }
  
  rows_list <- list()
  for (i in seq_len(n)) {
    zi <- z_vec[i]
    dy <- numeric(m)
    for (j in seq_len(m)) {
      mu_inc <- zi * dg_vec[j]
      lambda_inc <- (zi^2) * eta * (dg_vec[j]^2)
      dy[j] <- r_ig(1, mu = mu_inc, lambda = lambda_inc)
    }
    
    cum_y <- c(0, cumsum(dy))
    df_u <- data.frame(
      unit = i,
      t = times,
      increment = c(0, dy),
      degradation = cum_y,
      frailty_z = zi,
      stringsAsFactors = FALSE
    )
    rows_list[[i]] <- df_u
  }
  
  res <- do.call(rbind, rows_list)
  rownames(res) <- NULL
  res
}