R/sampling_plans.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Reliability Acceptance Sampling Plan for Exponential Distribution (10.5.1)
#'
#' @param n Sample size placed on test.
#' @param T_star Fixed test termination time.
#' @param r Required failure limit.
#' @param alpha Producer's risk.
#' @param beta Consumer's risk.
#' @param theta0 Acceptable quality level (AQL) mean lifetime.
#' @param theta1 Rejectable quality level (RQL) mean lifetime.
#'
#' @return S3 object of class \code{sampling_plan} with acceptance constant c, operating characteristic (OC) values, and decision rule.
#' @export
#'
#' @examples
#' sampling_plan_exponential(
#'   n = 20, T_star = 1.0, r = 10,
#'   alpha = 0.05, beta = 0.10,
#'   theta0 = 2.0, theta1 = 0.8
#' )
sampling_plan_exponential <- function(n, T_star, r, alpha = 0.05, beta = 0.10, theta0, theta1) {
  # Acceptance decision rule: Accept H0 if MLE hat_theta > c
  # Critical value c determined via Chi-Square pivoting under Exponential censoring
  df <- 2 * r
  chi_alpha <- stats::qchisq(alpha, df = df)
  c_crit <- (2 * r * theta0) / chi_alpha
  
  # Operating characteristic (OC) function L(theta) = P_theta(hat_theta > c)
  oc_theta0 <- 1 - stats::pchisq(2 * r * theta0 / c_crit, df = df)
  oc_theta1 <- 1 - stats::pchisq(2 * r * theta1 / c_crit, df = df)
  
  res <- list(
    critical_value = c_crit,
    n = n,
    r = r,
    T_star = T_star,
    oc_AQL = oc_theta0,
    oc_RQL = oc_theta1,
    producer_risk_achieved = 1 - oc_theta0,
    consumer_risk_achieved = oc_theta1,
    distribution = "Exponential"
  )
  class(res) <- "sampling_plan"
  return(res)
}

#' Reliability Acceptance Sampling Plan for Weibull Distribution (10.5.2)
#'
#' @param n Sample size.
#' @param T_star Termination time.
#' @param r Minimum failure limit.
#' @param beta_shape Known shape parameter of Weibull distribution.
#' @param alpha Producer's risk.
#' @param beta Consumer's risk.
#' @param theta0 AQL scale parameter.
#' @param theta1 RQL scale parameter.
#'
#' @return S3 object of class \code{sampling_plan}.
#' @export
#'
#' @examples
#' sampling_plan_weibull(
#'   n = 25, T_star = 1.5, r = 12, beta_shape = 1.5,
#'   alpha = 0.05, beta = 0.10, theta0 = 3.0, theta1 = 1.0
#' )
sampling_plan_weibull <- function(n, T_star, r, beta_shape, alpha = 0.05, beta = 0.10, theta0, theta1) {
  # Transformation of Weibull to Exponential via Y = X^beta_shape
  T_star_trans <- T_star^beta_shape
  theta0_trans <- theta0^beta_shape
  theta1_trans <- theta1^beta_shape
  
  plan_exp <- sampling_plan_exponential(
    n = n, T_star = T_star_trans, r = r,
    alpha = alpha, beta = beta,
    theta0 = theta0_trans, theta1 = theta1_trans
  )
  
  plan_exp$critical_value <- (plan_exp$critical_value)^(1 / beta_shape)
  plan_exp$distribution <- "Weibull"
  plan_exp$weibull_shape <- beta_shape
  return(plan_exp)
}

#' Bayesian Reliability Acceptance Sampling Plan for Exponential lifetimes (10.5.3)
#'
#' @param n Sample size.
#' @param T_star Termination time limit.
#' @param r Failure count limit.
#' @param prior_shape Shape parameter of Inverse-Gamma prior.
#' @param prior_rate Rate parameter of Inverse-Gamma prior.
#' @param c_accept Acceptance threshold constant.
#'
#' @return S3 object of class \code{sampling_plan} with posterior risk and decision criteria.
#' @export
#'
#' @examples
#' bayesian_sampling_plan_exp(
#'   n = 20, T_star = 1.0, r = 10,
#'   prior_shape = 2, prior_rate = 1, c_accept = 1.2
#' )
bayesian_sampling_plan_exp <- function(n, T_star, r, prior_shape = 2, prior_rate = 1, c_accept = 1.2) {
  # Posterior distribution of theta given r failures and total time on test TTT
  # Inverse-Gamma(prior_shape + r, prior_rate + TTT)
  post_shape <- prior_shape + r
  # Expected total time on test under T_star cutoff
  expected_ttt <- r * T_star * 0.8
  post_rate <- prior_rate + expected_ttt
  
  bayesian_accept_prob <- 1 - stats::pgamma(1 / c_accept, shape = post_shape, rate = post_rate)
  
  res <- list(
    critical_value = c_accept,
    post_shape = post_shape,
    post_rate = post_rate,
    posterior_acceptance_prob = bayesian_accept_prob,
    distribution = "Bayesian Exponential"
  )
  class(res) <- "sampling_plan"
  return(res)
}