R/rSSP.R

In rSSP: Single Acceptance Sampling Plans for Time-Truncated Life Test

# ============================================================
# SINGLE ACCEPTANCE SAMPLING INSPECTION PLAN
# ============================================================


#' Single Acceptance Sampling Inspection Plan
#'
#' Calculates the minimum required sample size for a single
#' acceptance sampling inspection plan under a time-truncated
#' life test using the producer's risk condition Pa <= be.
#'
#' The probability p is user-defined. Therefore, no particular
#' lifetime distribution is embedded in this function.
#'
#' @param p User-defined probability of failure before the
#'   termination time. It must lie between 0 and 1.
#' @param a Termination ratio, defined as a = t/theta0.
#'   It must be a positive numeric value.
#' @param b Quality ratio, defined as b = theta/theta0.
#'   It must be a positive numeric value.
#' @param be Producer's risk. It must be a value between 0 and 1.
#' @param c Acceptance number. It must be a non-negative integer.
#'
#' @return A data frame containing a, b, p, be, c, required
#'   sample size n, and probability of acceptance Pa.
#'
#' @examples
#'
#' # Example 1: User-defined failure probabilities
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' single_asip(
#'   p = p,
#'   a = c(0.5, 1, 1.5, 2),
#'   b = 1,
#'   be = 0.25,
#'   c = 0
#' )
#'
#' # Example 2: Failure probabilities obtained from
#' # a Burr Type X lifetime distribution
#'
#' al <- 0.5
#' b <- 1
#' a <- c(0.4, 0.6, 0.8, 1, 1.5, 2, 2.5, 3)
#'
#' h <- sqrt(
#'   -log(
#'     1 - (1/2)^(1/al)
#'   )
#' )
#'
#' p <- (1 - exp(-(a * h / b)^2))^al
#'
#' single_asip(
#'   p = p,
#'   a = a,
#'   b = b,
#'   be = 0.25,
#'   c = 6
#' )
#' @importFrom stats pbinom
#' @importFrom graphics grid
#' @export


single_asip <- function(
    p,
    a,
    b,
    be = 0.25,
    c = 0) {

  # ----------------------------------------------------------
  # Input checking
  # ----------------------------------------------------------

  if (!is.numeric(p) || any(!is.finite(p))) {
    stop("'p' must be a numeric value or numeric vector.")
  }

  if (any(p < 0 | p > 1)) {
    stop("'p' must lie between 0 and 1.")
  }


  if (!is.numeric(a) ||
      any(!is.finite(a)) ||
      any(a <= 0)) {

    stop("'a' must contain positive finite values.")
  }


  if (!is.numeric(b) ||
      any(!is.finite(b)) ||
      any(b <= 0)) {

    stop("'b' must contain positive finite values.")
  }


  if (!is.numeric(be) ||
      length(be) != 1 ||
      !is.finite(be) ||
      be <= 0 ||
      be >= 1) {

    stop("'be' must be between 0 and 1.")
  }


  if (!is.numeric(c) ||
      length(c) != 1 ||
      !is.finite(c) ||
      c < 0 ||
      c != floor(c)) {

    stop("'c' must be a non-negative integer.")
  }


  # ----------------------------------------------------------
  # Make a and b compatible with p
  # ----------------------------------------------------------

  if (length(a) == 1) {
    a <- rep(a, length(p))
  }

  if (length(b) == 1) {
    b <- rep(b, length(p))
  }


  if (length(a) != length(p) ||
      length(b) != length(p)) {

    stop("'a', 'b', and 'p' must have compatible lengths.")
  }


  # ----------------------------------------------------------
  # Find minimum sample size
  # ----------------------------------------------------------

  find_n <- function(p_value) {

    n <- 0

    repeat {

      n <- n + 1

      # Probability of acceptance
      Pa <- stats::pbinom(
        q = c,
        size = n,
        prob = p_value
      )


      # Producer's risk condition
      if (Pa <= be) {

        return(
          c(
            n = n,
            Pa = Pa
          )
        )
      }
    }
  }


  # ----------------------------------------------------------
  # Calculate n for each p
  # ----------------------------------------------------------

  result <- t(
    sapply(
      p,
      find_n
    )
  )


  # ----------------------------------------------------------
  # Final output
  # ----------------------------------------------------------

  output <- data.frame(
    a = a,
    b = b,
    p = p,
    be = be,
    c = c,
    n = result[, "n"],
    Pa = result[, "Pa"]
  )


  return(output)
}



# ============================================================
# PLOT SAMPLE SIZE AGAINST TERMINATION RATIO
# ============================================================


#' Plot Sample Size Against Termination Ratio
#'
#' Calculates the minimum sample size for different termination
#' ratios and plots the resulting sample size against the
#' termination ratio.
#'
#' The user supplies the corresponding values of p calculated
#' from their own lifetime distribution.
#'
#' @param p User-defined probability of failure corresponding
#'   to each termination ratio.
#' @param a Numeric vector of termination ratios.
#' @param b Quality ratio. It may be a single value or a vector
#'   having the same length as a.
#' @param be Producer's risk.
#' @param c Acceptance number.
#' @param ylim Numeric vector of length two specifying the limits
#'   of the y-axis. Default is c(0, 50).
#'
#' @return A data frame containing a, b, p and the required
#'   sample size n.
#'
#' @examples
#'
#' # Example 1: User-defined failure probabilities
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' plot_single_asip(
#'   p = p,
#'   a = c(0.5, 1, 1.5, 2),
#'   b = 1,
#'   be = 0.25,
#'   c = 0,
#'   ylim=c(0,50)
#' )
#'
#' # Example 2: Failure probabilities obtained from
#' # a Burr Type X lifetime distribution
#'
#' al <- 0.5
#' b <- 1
#' a <- c(0.4, 0.6, 0.8, 1, 1.5, 2, 2.5, 3)
#'
#' h <- sqrt(
#'   -log(
#'     1 - (1/2)^(1/al)
#'   )
#' )
#'
#' p <- (1 - exp(-(a * h / b)^2))^al
#'
#' plot_single_asip(
#'   p = p,
#'   a = a,
#'   b = b,
#'   be = 0.25,
#'   c = 6,
#'   ylim = c(0, 50)
#' )
#'
#' @export


plot_single_asip <- function(
    p,
    a,
    b = 1,
    be = 0.25,
    c = 0,
    ylim=c(0,50)) {


  # Calculate sample sizes
  result <- single_asip(
    p = p,
    a = a,
    b = b,
    be = be,
    c = c
  )


  # ----------------------------------------------------------
  # Plot a against n
  # ----------------------------------------------------------

  plot(
    result$a,
    result$n,
    type = "b",
    pch = 19,
    ylim=ylim,
    xlab = "Termination Ratio (a)",
    ylab = "Required Sample Size (n)",
    main = "Termination Ratio vs Required Sample Size"
  )

  graphics::grid()


  return(result)
}