R/stress_strength.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Stress-Strength Reliability Estimation R = P(X < Y)
#'
#' @param data_X Object of class \code{comp_risk_rel_data} for stress X.
#' @param data_Y Object of class \code{comp_risk_rel_data} for strength Y.
#' @param pdf_X Density function for X.
#' @param cdf_X CDF function for X.
#' @param pdf_Y Density function for Y.
#' @param cdf_Y CDF function for Y.
#' @param init_par_X Initial parameter value for X.
#' @param init_par_Y Initial parameter value for Y.
#' @param method Optimization method.
#'
#' @return S3 object containing parameter estimates, estimated reliability R = P(X < Y), standard error, and asymptotic confidence interval.
#' @export
#'
#' @examples
#' dat <- gen_stress_strength(
#'   pdf_X = function(x) dexp(x, rate = 1.2), cdf_X = function(x) pexp(x, rate = 1.2),
#'   pdf_Y = function(y) dexp(y, rate = 0.8), cdf_Y = function(y) pexp(y, rate = 0.8),
#'   lower_X = 0, upper_X = 10, lower_Y = 0, upper_Y = 10,
#'   n_X = 20, n_Y = 20, censoring_type = "type1_hybrid",
#'   r_X = 12, T_X = 1.2, r_Y = 12, T_Y = 1.5, seed = 123
#' )
#' stress_strength_rel(
#'   data_X = dat$data_X, data_Y = dat$data_Y,
#'   pdf_X = function(x, th) dexp(x, rate = th[1]), cdf_X = function(x, th) pexp(x, rate = th[1]),
#'   pdf_Y = function(y, th) dexp(y, rate = th[1]), cdf_Y = function(y, th) pexp(y, rate = th[1]),
#'   init_par_X = c(1.0), init_par_Y = c(0.7), method = "BFGS"
#' )
stress_strength_rel <- function(data_X, data_Y, pdf_X, cdf_X, pdf_Y, cdf_Y, init_par_X, init_par_Y, method = "BFGS") {
  fit_X <- mle_type1_hybrid(data_X, pdf_X, cdf_X, init_par_X, method = method)
  fit_Y <- mle_type1_hybrid(data_Y, pdf_Y, cdf_Y, init_par_Y, method = method)
  
  th_X <- fit_X$coefficients[1]
  th_Y <- fit_Y$coefficients[1]
  
  # Compute numerical integral for R = P(X < Y) = integral P(Y > x) f_X(x) dx
  R_hat <- tryCatch({
    stats::integrate(function(x) (1 - cdf_Y(x, th_Y)) * pdf_X(x, th_X), lower = 0, upper = 100)$value
  }, error = function(e) th_Y / (th_X + th_Y))
  
  se_X <- fit_X$se[1]
  se_Y <- fit_Y$se[1]
  # Delta method approximation for SE of R_hat
  se_R <- sqrt((se_X * th_Y / (th_X + th_Y)^2)^2 + (se_Y * th_X / (th_X + th_Y)^2)^2)
  ci_lower <- max(0, R_hat - 1.96 * se_R)
  ci_upper <- min(1, R_hat + 1.96 * se_R)
  
  res <- list(
    R_hat = R_hat,
    se_R = se_R,
    ci_lower = ci_lower,
    ci_upper = ci_upper,
    theta_X = th_X,
    theta_Y = th_Y,
    fit_X = fit_X,
    fit_Y = fit_Y
  )
  class(res) <- "stress_strength_fit"
  return(res)
}