Nothing
R/rMDSPTT.R
In rMDSPTT: Multiple Dependent State Sampling Inspection Plan for Time Truncated Life Test
Defines functions
plot_compare_mds_ssp
compare_mds_ssp
plot_mds_oc
plot_mds_n
mds_oc
mds_asip
.mds_pa
Documented in
compare_mds_ssp
mds_asip
mds_oc
plot_compare_mds_ssp
plot_mds_n
plot_mds_oc
# ============================================================
# rMDSP
# Multiple Dependent State Sampling Inspection Plan
# ============================================================
# ============================================================
# Internal function: MDS probability of acceptance
# ============================================================
.mds_pa <- function(A, B, i) {
# ----------------------------------------------------------
# New MDS acceptance probability from the proposed method:
#
# Pa = A + B A^i
#
# where
# A = P(D <= c1)
# B = P(c1 < D <= c2)
# ----------------------------------------------------------
if (!is.finite(A) || !is.finite(B) ||
A < 0 || A > 1 || B < 0 || B > 1) {
return(NA_real_)
}
Pa <- A + B * A^i
# Protect against numerical round-off outside [0,1].
Pa <- min(max(Pa, 0), 1)
Pa
}
# ============================================================
# MDS: Minimum sample size
# ============================================================
#' Multiple Dependent State Sampling Plan (MDS)
#'
#' Calculates the minimum sample size for a Multiple Dependent
#' State Sampling (MDS) plan for inspection by attributes under
#' a time-truncated life test.
#'
#' The failure probability before the termination time is
#' supplied directly by the user through `p`. Thus, the function
#' is distribution-free.
#'
#' The MDS plan is specified by `(n, c1, c2, i)`.
#'
#' Let D denote the number of defective units in a sample.
#' Under the binomial model,
#'
#' \deqn{
#' D\sim Binomial(n,p)
#' }{
#' D ~ Binomial(n,p)
#' }
#'
#' Define
#'
#' \deqn{
#' A=P(D\leq c_1)
#' }{
#' A=P(D<=c1)
#' }
#'
#' and
#'
#' \deqn{
#' B=P(c_1<D\leq c_2).
#' }{
#' B=P(c1<D<=c2)
#' }
#'
#' The probability of acceptance of the MDS plan is obtained
#' from
#'
#' \deqn{
#' P_a=A+B A^i.
#' }{
#' Pa=A+B A^i
#' }
#'
#' The minimum sample size is the smallest `n` satisfying
#'
#' \deqn{
#' P_a\leq\beta.
#' }{
#' Pa<=beta
#' }
#'
#' For this MDS plan, the sample size of every inspected lot is
#' fixed at `n`. Therefore, the average sample number (ASN) is
#' exactly equal to `n`.
#'
#' @param p User-defined probability of failure before the
#' termination time. It must lie strictly between 0 and 1.
#' @param a Termination ratio, defined as `a = t/theta0`.
#' It must contain positive values.
#' @param b Quality ratio, defined as `b = theta/theta0`.
#' It must contain positive values. It is included to identify
#' the design condition associated with `p`.
#' @param i Number of preceding lots considered in the MDS
#' decision. It must be a positive integer.
#' @param beta Consumer's risk. It must be between 0 and 1.
#' @param c1 First acceptance number. It must be a non-negative
#' integer.
#' @param c2 Second acceptance number. It must be greater than
#' or equal to `c1`.
#' @param n_max Maximum sample size to be searched.
#'
#' @return A data frame containing the design values, minimum
#' sample size `n`, ASN, probabilities `A` and `B`, and the
#' probability of acceptance `Pa`.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Weibull distribution
#' # ----------------------------------------------------------
#'
#' shape <- 2
#' b <- 1
#'
#' a <- c(
#' 0.5, 0.75, 1, 1.25,
#' 1.5, 1.75, 2
#' )
#'
#' p <- 1 - exp(
#' -((a / b)^shape)
#' )
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = b,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 3: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#' b <- 1
#'
#' a <- c(
#' 0.5, 0.75, 1, 1.25,
#' 1.5, 1.75, 2
#' )
#'
#' p <- (
#' 1 - exp(-a / b)
#' )^alpha
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = b,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' @importFrom stats pbinom
#' @export
mds_asip <- function(
p,
a,
b,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p) || any(!is.finite(p))) {
stop("'p' must be numeric and finite.")
}
if (any(p <= 0 | p >= 1)) {
stop("'p' must lie strictly 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(i) ||
length(i) != 1 ||
!is.finite(i) ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
!is.finite(beta) ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
!is.finite(c1) ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
!is.finite(c2) ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
!is.finite(n_max) ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive 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."
)
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find minimum sample size for one p
# ----------------------------------------------------------
find_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(
c(
A = A,
B = B,
n = n,
Pa = Pa
)
)
}
}
stop(
"No sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Calculate design results
# ----------------------------------------------------------
result <- t(
sapply(
p,
find_n
)
)
# ----------------------------------------------------------
# Final output
# ----------------------------------------------------------
output <- data.frame(
a = a,
b = b,
p = p,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
A = result[, "A"],
B = result[, "B"],
n = as.integer(result[, "n"]),
ASN = as.integer(result[, "n"]),
Pa = result[, "Pa"]
)
output
}
# ============================================================
# MDS: OC VALUES
# ============================================================
#' Operating Characteristic Values for an MDS Plan
#'
#' Calculates the operating characteristic (OC) values for a
#' Multiple Dependent State Sampling (MDS) plan.
#'
#' The sample size is first determined at the design condition
#' using `p_design` and the constraint `Pa <= beta`. This sample
#' size is then kept fixed while the failure probabilities in
#' `p_oc` are used to calculate the probability of acceptance
#' for the quality ratios in `b_oc`.
#'
#' Let
#'
#' \deqn{
#' A=P(D\leq c_1)
#' }{
#' A=P(D<=c1)
#' }
#'
#' and
#'
#' \deqn{
#' B=P(c_1<D\leq c_2).
#' }{
#' B=P(c1<D<=c2)
#' }
#'
#' The MDS probability of acceptance satisfies
#'
#' \deqn{
#' P_a=A+B A^i.
#' }{
#' Pa=A+B A^i
#' }
#'
#' The resulting `Pa` values are the OC values for the specified
#' quality ratios in `b_oc`, using `Pa = A + B A^i`.
#'
#' @param p_design Failure probability at the design condition.
#' It must lie strictly between 0 and 1.
#' @param a Termination ratio, defined as `a = t/theta0`.
#' @param p_oc Matrix or data frame of failure probabilities
#' used for OC calculation. Rows correspond to `a` and columns
#' correspond to `b_oc`.
#' @param b_oc Quality ratios used for OC calculation.
#' @param i Number of preceding lots considered in the MDS plan.
#' @param beta Consumer's risk.
#' @param c1 First acceptance number.
#' @param c2 Second acceptance number.
#' @param n_max Maximum sample size searched at the design
#' condition.
#'
#' @return A data frame containing `a`, `b`, `p`, fixed `n`,
#' `ASN`, `A`, `B`, and `Pa`. The `Pa` values are the OC values.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' p_design <- c(
#' 0.05, 0.10, 0.15, 0.20
#' )
#'
#' b_oc <- 2:12
#'
#' # User-defined probabilities for OC calculation.
#' # Rows correspond to a and columns correspond to b.
#' p_oc <- outer(
#' a,
#' b_oc,
#' function(a, b) pmin(0.95, a / (10 * b))
#' )
#'
#' # The n values are determined using p_design and then
#' # kept fixed for all b values.
#' #
#' # The Pa values represent the OC values.
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Weibull distribution
#' # ----------------------------------------------------------
#'
#' shape <- 2
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' # Design quality ratio
#' b_design <- 1
#'
#' # Failure probabilities at b = 1
#' p_design <- 1 - exp(
#' -((a / b_design)^shape)
#' )
#'
#' # Quality ratios for OC calculation
#' b_oc <- 2:12
#'
#' # Failure probabilities for each a and b
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' 1 - exp(-((a / b)^shape))
#' )
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 3: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' # Design quality ratio
#' b_design <- 1
#'
#' # Failure probabilities at b = 1
#' p_design <- (
#' 1 - exp(-a / b_design)
#' )^alpha
#'
#' # Quality ratios for OC calculation
#' b_oc <- 2:12
#'
#' # Failure probabilities for each a and b
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' (1 - exp(-a / b))^alpha
#' )
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' @importFrom stats pbinom
#' @export
mds_oc <- function(
p_design,
a,
p_oc,
b_oc,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p_design) ||
any(!is.finite(p_design))) {
stop("'p_design' must be numeric and finite.")
}
if (any(p_design <= 0 | p_design >= 1)) {
stop("'p_design' must lie strictly 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(p_oc) ||
any(!is.finite(p_oc))) {
stop("'p_oc' must be numeric and finite.")
}
if (any(p_oc <= 0 | p_oc >= 1)) {
stop(
"All values of 'p_oc' must lie strictly between 0 and 1."
)
}
if (!is.numeric(b_oc) ||
any(!is.finite(b_oc)) ||
any(b_oc <= 0)) {
stop("'b_oc' must contain positive finite values.")
}
if (!is.numeric(i) ||
length(i) != 1 ||
!is.finite(i) ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
!is.finite(beta) ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
!is.finite(c1) ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
!is.finite(c2) ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
!is.finite(n_max) ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive integer.")
}
# ----------------------------------------------------------
# Convert p_oc to matrix
# ----------------------------------------------------------
p_oc <- as.matrix(p_oc)
if (length(p_design) != length(a)) {
stop(
"'p_design' and 'a' must have the same length."
)
}
if (nrow(p_oc) != length(a)) {
stop(
"The number of rows of 'p_oc' must equal the length of 'a'."
)
}
if (ncol(p_oc) != length(b_oc)) {
stop(
"The number of columns of 'p_oc' must equal the length of 'b_oc'."
)
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find design sample size
# ----------------------------------------------------------
find_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No sample size satisfies Pa <= beta within n_max."
)
}
n_design <- sapply(
p_design,
find_n
)
# ----------------------------------------------------------
# Calculate OC values using fixed design sample size
# ----------------------------------------------------------
output_list <- vector(
"list",
length(a) * length(b_oc)
)
counter <- 1
for (j in seq_along(a)) {
n <- n_design[j]
for (k in seq_along(b_oc)) {
p_value <- p_oc[j, k]
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
output_list[[counter]] <- data.frame(
a = a[j],
b = b_oc[k],
p = p_value,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
n = as.integer(n),
ASN = as.integer(n),
A = A,
B = B,
Pa = Pa
)
counter <- counter + 1
}
}
result <- do.call(
rbind,
output_list
)
rownames(result) <- NULL
result
}
# ============================================================
# Plot minimum sample size / ASN against termination ratio
# ============================================================
#' Plot MDS Sample Size Against Termination Ratio
#'
#' Produces a base R plot of the minimum sample size against the
#' termination ratio `a`.
#'
#' For the MDS plan considered here, ASN is exactly equal to the
#' sample size `n`. Therefore, the plot can also be interpreted
#' as ASN against the termination ratio.
#'
#' @param x Output from `mds_asip()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#' a <- c(0.5, 1, 1.5, 2)
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' result <- mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_mds_n(result)
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- ( 1- exp(-a / b))^alpha
#'
#' result <- mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'plot_mds_n(result)
#'
#' @export
plot_mds_n <- function(x, ...) {
if (!is.data.frame(x)) {
stop("'x' must be a data frame returned by 'mds_asip()'.")
}
required <- c("a", "n")
if (!all(required %in% names(x))) {
stop(
"'x' must contain columns 'a' and 'n'."
)
}
b_values <- unique(x$b)
if (length(b_values) == 1) {
graphics::plot(
x$a,
x$n,
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Minimum sample size (n)",
main = "MDS Sample Size vs. Termination Ratio",
...
)
} else {
graphics::plot(
x$a[x$b == b_values[1]],
x$n[x$b == b_values[1]],
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Minimum sample size (n)",
main = "MDS Sample Size vs. Termination Ratio",
...
)
if (length(b_values) > 1) {
for (j in 2:length(b_values)) {
graphics::lines(
x$a[x$b == b_values[j]],
x$n[x$b == b_values[j]],
type = "b",
pch = 19
)
}
graphics::legend(
"topright",
legend = paste0("b = ", b_values),
lty = 1,
pch = 19
)
}
}
invisible(x)
}
# ============================================================
# Plot MDS OC curve
# ============================================================
#'
#' Plot OC Values for an MDS Plan
#'
#' Produces a base R OC curve using the `Pa` values returned by
#' `mds_oc()`.
#'
#' Different termination ratios `a` are distinguished using
#' different plotting symbols (`pch`) and line types (`lty`).
#'
#' @param x Output from `mds_oc()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' shape <- 2
#' a <- c(0.5, 1, 1.5, 2)
#'
#' p_design <- 1 - exp(-(a / 1)^shape)
#'
#' b_oc <- 2:12
#'
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' 1 - exp(-((a / b)^shape))
#' )
#'
#' result <- mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_mds_oc(result)
#'
#' @export
plot_mds_oc <- function(x, ...) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.data.frame(x)) {
stop("'x' must be a data frame returned by 'mds_oc()'.")
}
required <- c("a", "b", "Pa")
if (!all(required %in% names(x))) {
stop(
"'x' must contain columns 'a', 'b', and 'Pa'."
)
}
# ----------------------------------------------------------
# Values of a
# ----------------------------------------------------------
a_values <- unique(x$a)
n_a <- length(a_values)
# ----------------------------------------------------------
# Plotting symbols and line types
#
# These are deliberately different so that the curves
# remain distinguishable even in black-and-white printing.
# ----------------------------------------------------------
pch_values <- c(
1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18
)
lty_values <- c(
1, 2, 3, 4, 5, 6
)
# ----------------------------------------------------------
# If there are more curves than available symbols,
# recycle the symbols and line types.
# ----------------------------------------------------------
pch_values <- rep(
pch_values,
length.out = n_a
)
lty_values <- rep(
lty_values,
length.out = n_a
)
# ----------------------------------------------------------
# Plot first curve
# ----------------------------------------------------------
first <- a_values[1]
current <- x[x$a == first, ]
plot(
current$b,
current$Pa,
type = "b",
pch = pch_values[1],
lty = lty_values[1],
lwd = 1.2,
ylim = c(0, 1),
xlab = "Quality ratio (b)",
ylab = "Probability of acceptance (Pa)",
main = "OC Curve of MDS Plan",
...
)
# ----------------------------------------------------------
# Add remaining curves
# ----------------------------------------------------------
if (n_a > 1) {
for (j in 2:n_a) {
current <- x[x$a == a_values[j], ]
graphics::lines(
current$b,
current$Pa,
type = "b",
pch = pch_values[j],
lty = lty_values[j],
lwd = 1.2
)
}
}
# ----------------------------------------------------------
# Legend
# ----------------------------------------------------------
graphics::legend(
"bottomright",
legend = paste0("a = ", a_values),
pch = pch_values,
lty = lty_values,
lwd = 1.2,
bty = "n"
)
# ----------------------------------------------------------
# Return original data invisibly
# ----------------------------------------------------------
invisible(x)
}
# ============================================================
# Compare MDS and SSP
# ============================================================
#' Compare MDS and Single Sampling Plans
#'
#' Compares the minimum sample size required by an MDS plan and
#' a corresponding single sampling plan (SSP).
#'
#' The SSP uses acceptance number `c1`, while the MDS plan uses
#' `(c1, c2, i)`. Both plans use the same failure probability
#' `p` and consumer's risk `beta`.
#'
#' The MDS sample size is determined from
#'
#' \deqn{
#' P_a=A+B A^i\leq\beta.
#' }{
#' Pa=A+B A^i<=beta.
#' }
#'
#' The SSP sample size is determined from
#'
#' \deqn{
#' P(D\leq c_1)\leq\beta.
#' }{
#' P(D<=c1)<=beta.
#' }
#'
#' Since the sample size is fixed for each inspected lot under
#' both plans, ASN is equal to sample size for both plans.
#'
#' @param p User-defined failure probability.
#' @param a Termination ratio.
#' @param b Quality ratio.
#' @param i Number of preceding lots in the MDS plan.
#' @param beta Consumer's risk.
#' @param c1 MDS first acceptance number and SSP acceptance
#' number.
#' @param c2 MDS second acceptance number.
#' @param n_max Maximum sample size searched.
#'
#' @return A data frame containing the sample sizes required by
#' the MDS and SSP plans.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#' a <- c(0.5, 1, 1.5, 2)
#'
#' compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#'alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- ( 1- exp(-a / b))^alpha
#'
#' compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#' @importFrom stats pbinom
#' @export
compare_mds_ssp <- function(
p,
a,
b,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p) ||
any(!is.finite(p)) ||
any(p <= 0 | p >= 1)) {
stop(
"'p' must contain values strictly 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 (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."
)
}
if (!is.numeric(i) ||
length(i) != 1 ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive integer.")
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find MDS sample size
# ----------------------------------------------------------
find_mds_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A,
B,
i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No MDS sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Find SSP sample size
# ----------------------------------------------------------
find_ssp_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c1 > n) {
next
}
Pa <- stats::pbinom(
c1,
size = n,
prob = p_value
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No SSP sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Calculate both plans
# ----------------------------------------------------------
mds_n <- sapply(
p,
find_mds_n
)
ssp_n <- sapply(
p,
find_ssp_n
)
# ----------------------------------------------------------
# Final output
# ----------------------------------------------------------
data.frame(
a = a,
b = b,
p = p,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
MDS_n = as.integer(mds_n),
MDS_ASN = as.integer(mds_n),
SSP_n = as.integer(ssp_n),
SSP_ASN = as.integer(ssp_n)
)
}
# ============================================================
# Plot MDS versus SSP sample size
# ============================================================
#' Plot MDS and SSP Sample Size Comparison
#'
#' Produces a base R plot comparing MDS and SSP sample sizes
#' against the termination ratio.
#'
#' Since ASN equals sample size for both plans, the same plot
#' also represents the ASN comparison.
#'
#' @param x Output from `compare_mds_ssp()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#' a <- c(0.5, 1, 1.5, 2)
#'
#' result <- compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_compare_mds_ssp(result)
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#'alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- (1- exp(-a / b))^alpha
#' result <- compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#' plot_compare_mds_ssp(result)
#' @export
plot_compare_mds_ssp <- function(x, ...) {
if (!is.data.frame(x)) {
stop(
"'x' must be a data frame returned by 'compare_mds_ssp()'."
)
}
required <- c(
"a",
"MDS_n",
"SSP_n"
)
if (!all(required %in% names(x))) {
stop(
"'x' must contain 'a', 'MDS_n', and 'SSP_n'."
)
}
plot(
x$a,
x$MDS_n,
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Sample size (n)",
main = "MDS vs. SSP Sample Size",
...
)
graphics::lines(
x$a,
x$SSP_n,
type = "b",
pch = 17
)
graphics::legend(
"topright",
legend = c(
"MDS",
"SSP"
),
lty = 1,
pch = c(19, 17)
)
invisible(x)
}
Try the rMDSPTT package in your browser
Any scripts or data that you put into this service are public.
rMDSPTT documentation built on Oct. 2, 2026, 5:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot_compare_mds_ssp compare_mds_ssp plot_mds_oc plot_mds_n mds_oc mds_asip .mds_pa
Documented in compare_mds_ssp mds_asip mds_oc plot_compare_mds_ssp plot_mds_n plot_mds_oc
# ============================================================
# rMDSP
# Multiple Dependent State Sampling Inspection Plan
# ============================================================
# ============================================================
# Internal function: MDS probability of acceptance
# ============================================================
.mds_pa <- function(A, B, i) {
# ----------------------------------------------------------
# New MDS acceptance probability from the proposed method:
#
# Pa = A + B A^i
#
# where
# A = P(D <= c1)
# B = P(c1 < D <= c2)
# ----------------------------------------------------------
if (!is.finite(A) || !is.finite(B) ||
A < 0 || A > 1 || B < 0 || B > 1) {
return(NA_real_)
}
Pa <- A + B * A^i
# Protect against numerical round-off outside [0,1].
Pa <- min(max(Pa, 0), 1)
Pa
}
# ============================================================
# MDS: Minimum sample size
# ============================================================
#' Multiple Dependent State Sampling Plan (MDS)
#'
#' Calculates the minimum sample size for a Multiple Dependent
#' State Sampling (MDS) plan for inspection by attributes under
#' a time-truncated life test.
#'
#' The failure probability before the termination time is
#' supplied directly by the user through `p`. Thus, the function
#' is distribution-free.
#'
#' The MDS plan is specified by `(n, c1, c2, i)`.
#'
#' Let D denote the number of defective units in a sample.
#' Under the binomial model,
#'
#' \deqn{
#' D\sim Binomial(n,p)
#' }{
#' D ~ Binomial(n,p)
#' }
#'
#' Define
#'
#' \deqn{
#' A=P(D\leq c_1)
#' }{
#' A=P(D<=c1)
#' }
#'
#' and
#'
#' \deqn{
#' B=P(c_1<D\leq c_2).
#' }{
#' B=P(c1<D<=c2)
#' }
#'
#' The probability of acceptance of the MDS plan is obtained
#' from
#'
#' \deqn{
#' P_a=A+B A^i.
#' }{
#' Pa=A+B A^i
#' }
#'
#' The minimum sample size is the smallest `n` satisfying
#'
#' \deqn{
#' P_a\leq\beta.
#' }{
#' Pa<=beta
#' }
#'
#' For this MDS plan, the sample size of every inspected lot is
#' fixed at `n`. Therefore, the average sample number (ASN) is
#' exactly equal to `n`.
#'
#' @param p User-defined probability of failure before the
#' termination time. It must lie strictly between 0 and 1.
#' @param a Termination ratio, defined as `a = t/theta0`.
#' It must contain positive values.
#' @param b Quality ratio, defined as `b = theta/theta0`.
#' It must contain positive values. It is included to identify
#' the design condition associated with `p`.
#' @param i Number of preceding lots considered in the MDS
#' decision. It must be a positive integer.
#' @param beta Consumer's risk. It must be between 0 and 1.
#' @param c1 First acceptance number. It must be a non-negative
#' integer.
#' @param c2 Second acceptance number. It must be greater than
#' or equal to `c1`.
#' @param n_max Maximum sample size to be searched.
#'
#' @return A data frame containing the design values, minimum
#' sample size `n`, ASN, probabilities `A` and `B`, and the
#' probability of acceptance `Pa`.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Weibull distribution
#' # ----------------------------------------------------------
#'
#' shape <- 2
#' b <- 1
#'
#' a <- c(
#' 0.5, 0.75, 1, 1.25,
#' 1.5, 1.75, 2
#' )
#'
#' p <- 1 - exp(
#' -((a / b)^shape)
#' )
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = b,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 3: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#' b <- 1
#'
#' a <- c(
#' 0.5, 0.75, 1, 1.25,
#' 1.5, 1.75, 2
#' )
#'
#' p <- (
#' 1 - exp(-a / b)
#' )^alpha
#'
#' mds_asip(
#' p = p,
#' a = a,
#' b = b,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' @importFrom stats pbinom
#' @export
mds_asip <- function(
p,
a,
b,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p) || any(!is.finite(p))) {
stop("'p' must be numeric and finite.")
}
if (any(p <= 0 | p >= 1)) {
stop("'p' must lie strictly 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(i) ||
length(i) != 1 ||
!is.finite(i) ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
!is.finite(beta) ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
!is.finite(c1) ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
!is.finite(c2) ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
!is.finite(n_max) ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive 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."
)
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find minimum sample size for one p
# ----------------------------------------------------------
find_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(
c(
A = A,
B = B,
n = n,
Pa = Pa
)
)
}
}
stop(
"No sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Calculate design results
# ----------------------------------------------------------
result <- t(
sapply(
p,
find_n
)
)
# ----------------------------------------------------------
# Final output
# ----------------------------------------------------------
output <- data.frame(
a = a,
b = b,
p = p,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
A = result[, "A"],
B = result[, "B"],
n = as.integer(result[, "n"]),
ASN = as.integer(result[, "n"]),
Pa = result[, "Pa"]
)
output
}
# ============================================================
# MDS: OC VALUES
# ============================================================
#' Operating Characteristic Values for an MDS Plan
#'
#' Calculates the operating characteristic (OC) values for a
#' Multiple Dependent State Sampling (MDS) plan.
#'
#' The sample size is first determined at the design condition
#' using `p_design` and the constraint `Pa <= beta`. This sample
#' size is then kept fixed while the failure probabilities in
#' `p_oc` are used to calculate the probability of acceptance
#' for the quality ratios in `b_oc`.
#'
#' Let
#'
#' \deqn{
#' A=P(D\leq c_1)
#' }{
#' A=P(D<=c1)
#' }
#'
#' and
#'
#' \deqn{
#' B=P(c_1<D\leq c_2).
#' }{
#' B=P(c1<D<=c2)
#' }
#'
#' The MDS probability of acceptance satisfies
#'
#' \deqn{
#' P_a=A+B A^i.
#' }{
#' Pa=A+B A^i
#' }
#'
#' The resulting `Pa` values are the OC values for the specified
#' quality ratios in `b_oc`, using `Pa = A + B A^i`.
#'
#' @param p_design Failure probability at the design condition.
#' It must lie strictly between 0 and 1.
#' @param a Termination ratio, defined as `a = t/theta0`.
#' @param p_oc Matrix or data frame of failure probabilities
#' used for OC calculation. Rows correspond to `a` and columns
#' correspond to `b_oc`.
#' @param b_oc Quality ratios used for OC calculation.
#' @param i Number of preceding lots considered in the MDS plan.
#' @param beta Consumer's risk.
#' @param c1 First acceptance number.
#' @param c2 Second acceptance number.
#' @param n_max Maximum sample size searched at the design
#' condition.
#'
#' @return A data frame containing `a`, `b`, `p`, fixed `n`,
#' `ASN`, `A`, `B`, and `Pa`. The `Pa` values are the OC values.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' p_design <- c(
#' 0.05, 0.10, 0.15, 0.20
#' )
#'
#' b_oc <- 2:12
#'
#' # User-defined probabilities for OC calculation.
#' # Rows correspond to a and columns correspond to b.
#' p_oc <- outer(
#' a,
#' b_oc,
#' function(a, b) pmin(0.95, a / (10 * b))
#' )
#'
#' # The n values are determined using p_design and then
#' # kept fixed for all b values.
#' #
#' # The Pa values represent the OC values.
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Weibull distribution
#' # ----------------------------------------------------------
#'
#' shape <- 2
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' # Design quality ratio
#' b_design <- 1
#'
#' # Failure probabilities at b = 1
#' p_design <- 1 - exp(
#' -((a / b_design)^shape)
#' )
#'
#' # Quality ratios for OC calculation
#' b_oc <- 2:12
#'
#' # Failure probabilities for each a and b
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' 1 - exp(-((a / b)^shape))
#' )
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#'
#' # ----------------------------------------------------------
#' # Example 3: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#'
#' a <- c(0.5, 1, 1.5, 2)
#'
#' # Design quality ratio
#' b_design <- 1
#'
#' # Failure probabilities at b = 1
#' p_design <- (
#' 1 - exp(-a / b_design)
#' )^alpha
#'
#' # Quality ratios for OC calculation
#' b_oc <- 2:12
#'
#' # Failure probabilities for each a and b
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' (1 - exp(-a / b))^alpha
#' )
#'
#' mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' @importFrom stats pbinom
#' @export
mds_oc <- function(
p_design,
a,
p_oc,
b_oc,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p_design) ||
any(!is.finite(p_design))) {
stop("'p_design' must be numeric and finite.")
}
if (any(p_design <= 0 | p_design >= 1)) {
stop("'p_design' must lie strictly 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(p_oc) ||
any(!is.finite(p_oc))) {
stop("'p_oc' must be numeric and finite.")
}
if (any(p_oc <= 0 | p_oc >= 1)) {
stop(
"All values of 'p_oc' must lie strictly between 0 and 1."
)
}
if (!is.numeric(b_oc) ||
any(!is.finite(b_oc)) ||
any(b_oc <= 0)) {
stop("'b_oc' must contain positive finite values.")
}
if (!is.numeric(i) ||
length(i) != 1 ||
!is.finite(i) ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
!is.finite(beta) ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
!is.finite(c1) ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
!is.finite(c2) ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
!is.finite(n_max) ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive integer.")
}
# ----------------------------------------------------------
# Convert p_oc to matrix
# ----------------------------------------------------------
p_oc <- as.matrix(p_oc)
if (length(p_design) != length(a)) {
stop(
"'p_design' and 'a' must have the same length."
)
}
if (nrow(p_oc) != length(a)) {
stop(
"The number of rows of 'p_oc' must equal the length of 'a'."
)
}
if (ncol(p_oc) != length(b_oc)) {
stop(
"The number of columns of 'p_oc' must equal the length of 'b_oc'."
)
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find design sample size
# ----------------------------------------------------------
find_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No sample size satisfies Pa <= beta within n_max."
)
}
n_design <- sapply(
p_design,
find_n
)
# ----------------------------------------------------------
# Calculate OC values using fixed design sample size
# ----------------------------------------------------------
output_list <- vector(
"list",
length(a) * length(b_oc)
)
counter <- 1
for (j in seq_along(a)) {
n <- n_design[j]
for (k in seq_along(b_oc)) {
p_value <- p_oc[j, k]
A <- stats::pbinom(
q = c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
q = c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A = A,
B = B,
i = i
)
output_list[[counter]] <- data.frame(
a = a[j],
b = b_oc[k],
p = p_value,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
n = as.integer(n),
ASN = as.integer(n),
A = A,
B = B,
Pa = Pa
)
counter <- counter + 1
}
}
result <- do.call(
rbind,
output_list
)
rownames(result) <- NULL
result
}
# ============================================================
# Plot minimum sample size / ASN against termination ratio
# ============================================================
#' Plot MDS Sample Size Against Termination Ratio
#'
#' Produces a base R plot of the minimum sample size against the
#' termination ratio `a`.
#'
#' For the MDS plan considered here, ASN is exactly equal to the
#' sample size `n`. Therefore, the plot can also be interpreted
#' as ASN against the termination ratio.
#'
#' @param x Output from `mds_asip()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#' a <- c(0.5, 1, 1.5, 2)
#' p <- c(0.05, 0.10, 0.15, 0.20)
#'
#' result <- mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_mds_n(result)
#'
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#' alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- ( 1- exp(-a / b))^alpha
#'
#' result <- mds_asip(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'plot_mds_n(result)
#'
#' @export
plot_mds_n <- function(x, ...) {
if (!is.data.frame(x)) {
stop("'x' must be a data frame returned by 'mds_asip()'.")
}
required <- c("a", "n")
if (!all(required %in% names(x))) {
stop(
"'x' must contain columns 'a' and 'n'."
)
}
b_values <- unique(x$b)
if (length(b_values) == 1) {
graphics::plot(
x$a,
x$n,
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Minimum sample size (n)",
main = "MDS Sample Size vs. Termination Ratio",
...
)
} else {
graphics::plot(
x$a[x$b == b_values[1]],
x$n[x$b == b_values[1]],
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Minimum sample size (n)",
main = "MDS Sample Size vs. Termination Ratio",
...
)
if (length(b_values) > 1) {
for (j in 2:length(b_values)) {
graphics::lines(
x$a[x$b == b_values[j]],
x$n[x$b == b_values[j]],
type = "b",
pch = 19
)
}
graphics::legend(
"topright",
legend = paste0("b = ", b_values),
lty = 1,
pch = 19
)
}
}
invisible(x)
}
# ============================================================
# Plot MDS OC curve
# ============================================================
#'
#' Plot OC Values for an MDS Plan
#'
#' Produces a base R OC curve using the `Pa` values returned by
#' `mds_oc()`.
#'
#' Different termination ratios `a` are distinguished using
#' different plotting symbols (`pch`) and line types (`lty`).
#'
#' @param x Output from `mds_oc()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' shape <- 2
#' a <- c(0.5, 1, 1.5, 2)
#'
#' p_design <- 1 - exp(-(a / 1)^shape)
#'
#' b_oc <- 2:12
#'
#' p_oc <- sapply(
#' b_oc,
#' function(b)
#' 1 - exp(-((a / b)^shape))
#' )
#'
#' result <- mds_oc(
#' p_design = p_design,
#' a = a,
#' p_oc = p_oc,
#' b_oc = b_oc,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_mds_oc(result)
#'
#' @export
plot_mds_oc <- function(x, ...) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.data.frame(x)) {
stop("'x' must be a data frame returned by 'mds_oc()'.")
}
required <- c("a", "b", "Pa")
if (!all(required %in% names(x))) {
stop(
"'x' must contain columns 'a', 'b', and 'Pa'."
)
}
# ----------------------------------------------------------
# Values of a
# ----------------------------------------------------------
a_values <- unique(x$a)
n_a <- length(a_values)
# ----------------------------------------------------------
# Plotting symbols and line types
#
# These are deliberately different so that the curves
# remain distinguishable even in black-and-white printing.
# ----------------------------------------------------------
pch_values <- c(
1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18
)
lty_values <- c(
1, 2, 3, 4, 5, 6
)
# ----------------------------------------------------------
# If there are more curves than available symbols,
# recycle the symbols and line types.
# ----------------------------------------------------------
pch_values <- rep(
pch_values,
length.out = n_a
)
lty_values <- rep(
lty_values,
length.out = n_a
)
# ----------------------------------------------------------
# Plot first curve
# ----------------------------------------------------------
first <- a_values[1]
current <- x[x$a == first, ]
plot(
current$b,
current$Pa,
type = "b",
pch = pch_values[1],
lty = lty_values[1],
lwd = 1.2,
ylim = c(0, 1),
xlab = "Quality ratio (b)",
ylab = "Probability of acceptance (Pa)",
main = "OC Curve of MDS Plan",
...
)
# ----------------------------------------------------------
# Add remaining curves
# ----------------------------------------------------------
if (n_a > 1) {
for (j in 2:n_a) {
current <- x[x$a == a_values[j], ]
graphics::lines(
current$b,
current$Pa,
type = "b",
pch = pch_values[j],
lty = lty_values[j],
lwd = 1.2
)
}
}
# ----------------------------------------------------------
# Legend
# ----------------------------------------------------------
graphics::legend(
"bottomright",
legend = paste0("a = ", a_values),
pch = pch_values,
lty = lty_values,
lwd = 1.2,
bty = "n"
)
# ----------------------------------------------------------
# Return original data invisibly
# ----------------------------------------------------------
invisible(x)
}
# ============================================================
# Compare MDS and SSP
# ============================================================
#' Compare MDS and Single Sampling Plans
#'
#' Compares the minimum sample size required by an MDS plan and
#' a corresponding single sampling plan (SSP).
#'
#' The SSP uses acceptance number `c1`, while the MDS plan uses
#' `(c1, c2, i)`. Both plans use the same failure probability
#' `p` and consumer's risk `beta`.
#'
#' The MDS sample size is determined from
#'
#' \deqn{
#' P_a=A+B A^i\leq\beta.
#' }{
#' Pa=A+B A^i<=beta.
#' }
#'
#' The SSP sample size is determined from
#'
#' \deqn{
#' P(D\leq c_1)\leq\beta.
#' }{
#' P(D<=c1)<=beta.
#' }
#'
#' Since the sample size is fixed for each inspected lot under
#' both plans, ASN is equal to sample size for both plans.
#'
#' @param p User-defined failure probability.
#' @param a Termination ratio.
#' @param b Quality ratio.
#' @param i Number of preceding lots in the MDS plan.
#' @param beta Consumer's risk.
#' @param c1 MDS first acceptance number and SSP acceptance
#' number.
#' @param c2 MDS second acceptance number.
#' @param n_max Maximum sample size searched.
#'
#' @return A data frame containing the sample sizes required by
#' the MDS and SSP plans.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#' a <- c(0.5, 1, 1.5, 2)
#'
#' compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#'alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- ( 1- exp(-a / b))^alpha
#'
#' compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#' @importFrom stats pbinom
#' @export
compare_mds_ssp <- function(
p,
a,
b,
i = 1,
beta = 0.25,
c1 = 0,
c2 = 1,
n_max = 10000) {
# ----------------------------------------------------------
# Input checking
# ----------------------------------------------------------
if (!is.numeric(p) ||
any(!is.finite(p)) ||
any(p <= 0 | p >= 1)) {
stop(
"'p' must contain values strictly 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 (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."
)
}
if (!is.numeric(i) ||
length(i) != 1 ||
i <= 0 ||
i != floor(i)) {
stop("'i' must be a positive integer.")
}
if (!is.numeric(beta) ||
length(beta) != 1 ||
beta <= 0 ||
beta >= 1) {
stop("'beta' must be between 0 and 1.")
}
if (!is.numeric(c1) ||
length(c1) != 1 ||
c1 < 0 ||
c1 != floor(c1)) {
stop("'c1' must be a non-negative integer.")
}
if (!is.numeric(c2) ||
length(c2) != 1 ||
c2 < c1 ||
c2 != floor(c2)) {
stop(
"'c2' must be an integer greater than or equal to 'c1'."
)
}
if (!is.numeric(n_max) ||
length(n_max) != 1 ||
n_max <= 0 ||
n_max != floor(n_max)) {
stop("'n_max' must be a positive integer.")
}
i <- as.integer(i)
c1 <- as.integer(c1)
c2 <- as.integer(c2)
n_max <- as.integer(n_max)
# ----------------------------------------------------------
# Find MDS sample size
# ----------------------------------------------------------
find_mds_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c2 > n) {
next
}
A <- stats::pbinom(
c1,
size = n,
prob = p_value
)
B <- stats::pbinom(
c2,
size = n,
prob = p_value
) - A
Pa <- .mds_pa(
A,
B,
i
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No MDS sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Find SSP sample size
# ----------------------------------------------------------
find_ssp_n <- function(p_value) {
for (n in seq_len(n_max)) {
if (c1 > n) {
next
}
Pa <- stats::pbinom(
c1,
size = n,
prob = p_value
)
if (is.finite(Pa) &&
Pa <= beta) {
return(n)
}
}
stop(
"No SSP sample size satisfies Pa <= beta within n_max."
)
}
# ----------------------------------------------------------
# Calculate both plans
# ----------------------------------------------------------
mds_n <- sapply(
p,
find_mds_n
)
ssp_n <- sapply(
p,
find_ssp_n
)
# ----------------------------------------------------------
# Final output
# ----------------------------------------------------------
data.frame(
a = a,
b = b,
p = p,
i = i,
beta = beta,
c1 = c1,
c2 = c2,
MDS_n = as.integer(mds_n),
MDS_ASN = as.integer(mds_n),
SSP_n = as.integer(ssp_n),
SSP_ASN = as.integer(ssp_n)
)
}
# ============================================================
# Plot MDS versus SSP sample size
# ============================================================
#' Plot MDS and SSP Sample Size Comparison
#'
#' Produces a base R plot comparing MDS and SSP sample sizes
#' against the termination ratio.
#'
#' Since ASN equals sample size for both plans, the same plot
#' also represents the ASN comparison.
#'
#' @param x Output from `compare_mds_ssp()`.
#' @param ... Additional graphical arguments passed to `plot()`.
#'
#' @return Invisibly returns the supplied data frame.
#'
#' @examples
#'
#' # ----------------------------------------------------------
#' # Example 1: User-defined failure probabilities
#' # ----------------------------------------------------------
#'
#' p <- c(0.05, 0.10, 0.15, 0.20)
#' a <- c(0.5, 1, 1.5, 2)
#'
#' result <- compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#'
#' plot_compare_mds_ssp(result)
#'
#' # ----------------------------------------------------------
#' # Example 2: Generalized Exponential distribution
#' # ----------------------------------------------------------
#'
#'alpha <- 2
#' b <- 1
#' a <- c(0.5, 0.75, 1, 1.25, 1.5, 1.75, 2)
#' p <- (1- exp(-a / b))^alpha
#' result <- compare_mds_ssp(
#' p = p,
#' a = a,
#' b = 1,
#' i = 3,
#' beta = 0.25,
#' c1 = 0,
#' c2 = 1
#' )
#' plot_compare_mds_ssp(result)
#' @export
plot_compare_mds_ssp <- function(x, ...) {
if (!is.data.frame(x)) {
stop(
"'x' must be a data frame returned by 'compare_mds_ssp()'."
)
}
required <- c(
"a",
"MDS_n",
"SSP_n"
)
if (!all(required %in% names(x))) {
stop(
"'x' must contain 'a', 'MDS_n', and 'SSP_n'."
)
}
plot(
x$a,
x$MDS_n,
type = "b",
pch = 19,
xlab = "Termination ratio (a)",
ylab = "Sample size (n)",
main = "MDS vs. SSP Sample Size",
...
)
graphics::lines(
x$a,
x$SSP_n,
type = "b",
pch = 17
)
graphics::legend(
"topright",
legend = c(
"MDS",
"SSP"
),
lty = 1,
pch = c(19, 17)
)
invisible(x)
}
Try the rMDSPTT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.