tests/testthat/test-MDSP.R

In rMDSPTT: Multiple Dependent State Sampling Inspection Plan for Time Truncated Life Test

# ============================================================
# TESTS FOR rMDSP
# Multiple Dependent State Sampling Inspection Plan
# New method: Pa = A + B A^i
# ============================================================

library(testthat)


# ------------------------------------------------------------
# 1. Internal MDS probability calculation
# ------------------------------------------------------------

test_that("MDS acceptance probability uses Pa = A + B A^i", {

  A <- 0.60
  B <- 0.25
  i <- 3

  expected <- A + B * A^i

  expect_equal(
    .mds_pa(A, B, i),
    expected,
    tolerance = 1e-12
  )
})


# ------------------------------------------------------------
# 2. MDS minimum sample size
# ------------------------------------------------------------

test_that("mds_asip satisfies Pa <= beta", {

  result <- mds_asip(
    p = c(0.05, 0.10, 0.15, 0.20),
    a = c(0.5, 1, 1.5, 2),
    b = 1,
    i = 3,
    beta = 0.25,
    c1 = 0,
    c2 = 1,
    n_max = 1000
  )

  expect_s3_class(result, "data.frame")
  expect_equal(nrow(result), 4)
  expect_true(all(result$Pa <= 0.25 + 1e-12))
  expect_true(all(result$n >= 2))
  expect_equal(result$ASN, result$n)
})


# ------------------------------------------------------------
# 3. Exact verification of A, B and Pa
# ------------------------------------------------------------

test_that("mds_asip calculates A, B and Pa correctly", {

  p <- 0.10
  a <- 1
  b <- 1
  i <- 3
  c1 <- 0
  c2 <- 1
  beta <- 0.25

  result <- mds_asip(
    p = p,
    a = a,
    b = b,
    i = i,
    beta = beta,
    c1 = c1,
    c2 = c2,
    n_max = 1000
  )

  n <- result$n

  A <- stats::pbinom(c1, size = n, prob = p)
  B <- stats::pbinom(c2, size = n, prob = p) - A
  Pa <- A + B * A^i

  expect_equal(result$A, A, tolerance = 1e-12)
  expect_equal(result$B, B, tolerance = 1e-12)
  expect_equal(result$Pa, Pa, tolerance = 1e-12)
})


# ------------------------------------------------------------
# 4. Verify minimum n
# ------------------------------------------------------------

test_that("mds_asip returns the first n satisfying Pa <= beta", {

  p <- 0.10
  beta <- 0.25
  i <- 3
  c1 <- 0
  c2 <- 1

  result <- mds_asip(
    p = p,
    a = 1,
    b = 1,
    i = i,
    beta = beta,
    c1 = c1,
    c2 = c2,
    n_max = 1000
  )

  n <- result$n

  expect_lte(result$Pa, beta + 1e-12)

  if (n > c2) {
    A_prev <- stats::pbinom(c1, size = n - 1, prob = p)
    B_prev <- stats::pbinom(c2, size = n - 1, prob = p) - A_prev
    Pa_prev <- A_prev + B_prev * A_prev^i

    expect_gt(Pa_prev, beta)
  }
})


# ------------------------------------------------------------
# 5. Vectorized p and a
# ------------------------------------------------------------

test_that("mds_asip accepts vector p and a", {

  p <- c(0.05, 0.10, 0.15, 0.20)
  a <- c(0.5, 1, 1.5, 2)

  result <- mds_asip(
    p = p,
    a = a,
    b = 1,
    i = 3,
    beta = 0.25,
    c1 = 0,
    c2 = 1,
    n_max = 1000
  )

  expect_equal(result$p, p)
  expect_equal(result$a, a)
  expect_equal(nrow(result), length(p))
  expect_true(all(result$Pa <= 0.25 + 1e-12))
})


# ------------------------------------------------------------
# 6. General c1 and c2
# ------------------------------------------------------------

test_that("mds_asip works for general c1 and c2", {

  result <- mds_asip(
    p = 0.10,
    a = 1,
    b = 1,
    i = 2,
    beta = 0.25,
    c1 = 1,
    c2 = 2,
    n_max = 1000
  )

  n <- result$n
  A <- stats::pbinom(1, n, 0.10)
  B <- stats::pbinom(2, n, 0.10) - A

  expect_equal(result$Pa, A + B * A^2, tolerance = 1e-12)
  expect_lte(result$Pa, 0.25 + 1e-12)
})


# ------------------------------------------------------------
# 7. mds_oc uses the new method
# ------------------------------------------------------------

test_that("mds_oc uses Pa = A + B A^i", {

  a <- c(0.5, 1)
  p_design <- c(0.05, 0.10)
  b_oc <- c(1, 2, 3)

  p_oc <- sapply(
    b_oc,
    function(b) 1 - exp(-((a / b)^2))
  )

  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,
    n_max = 1000
  )

  expect_s3_class(result, "data.frame")
  expect_equal(nrow(result), length(a) * length(b_oc))
  expect_true(all(result$Pa >= 0 & result$Pa <= 1))

  for (j in seq_along(a)) {
    for (k in seq_along(b_oc)) {

      row <- result[result$a == a[j] & result$b == b_oc[k], ]
      n <- row$n
      p_value <- p_oc[j, k]

      A <- stats::pbinom(0, n, p_value)
      B <- stats::pbinom(1, n, p_value) - A
      expected <- A + B * A^3

      expect_equal(row$A, A, tolerance = 1e-12)
      expect_equal(row$B, B, tolerance = 1e-12)
      expect_equal(row$Pa, expected, tolerance = 1e-12)
    }
  }
})


# ------------------------------------------------------------
# 8. MDS versus SSP comparison uses new MDS method
# ------------------------------------------------------------

test_that("compare_mds_ssp uses the new MDS formulation", {

  p <- c(0.05, 0.10, 0.15)
  a <- c(0.5, 1, 1.5)

  result <- compare_mds_ssp(
    p = p,
    a = a,
    b = 1,
    i = 3,
    beta = 0.25,
    c1 = 0,
    c2 = 1,
    n_max = 1000
  )

  expect_s3_class(result, "data.frame")
  expect_true(all(result$MDS_n >= 2))
  expect_true(all(result$SSP_n >= 1))

  for (j in seq_along(p)) {

    n <- result$MDS_n[j]
    A <- stats::pbinom(0, n, p[j])
    B <- stats::pbinom(1, n, p[j]) - A
    Pa <- A + B * A^3

    expect_lte(Pa, 0.25 + 1e-12)
  }
})


# ------------------------------------------------------------
# 9. Input validation
# ------------------------------------------------------------

test_that("mds_asip validates inputs", {

  expect_error(
    mds_asip(
      p = 1.1, a = 1, b = 1, i = 3,
      beta = 0.25, c1 = 0, c2 = 1
    )
  )

  expect_error(
    mds_asip(
      p = 0.1, a = 1, b = 1, i = 0,
      beta = 0.25, c1 = 0, c2 = 1
    )
  )

  expect_error(
    mds_asip(
      p = 0.1, a = 1, b = 1, i = 3,
      beta = 0.25, c1 = 2, c2 = 1
    )
  )

  expect_error(
    mds_asip(
      p = 0.1, a = 1, b = 1, i = 3,
      beta = 0.25, c1 = -1, c2 = 1
    )
  )
})


# ------------------------------------------------------------
# 10. Plot functions run successfully
# ------------------------------------------------------------

test_that("MDS plot functions run without error", {

  result <- mds_asip(
    p = c(0.05, 0.10, 0.15),
    a = c(0.5, 1, 1.5),
    b = 1,
    i = 3,
    beta = 0.25,
    c1 = 0,
    c2 = 1,
    n_max = 1000
  )

  oc_result <- mds_oc(
    p_design = c(0.05, 0.10, 0.15),
    a = c(0.5, 1, 1.5),
    p_oc = sapply(
      1:3,
      function(b) 1 - exp(-((c(0.5, 1, 1.5) / b)^2))
    ),
    b_oc = 1:3,
    i = 3,
    beta = 0.25,
    c1 = 0,
    c2 = 1,
    n_max = 1000
  )

  f1 <- tempfile(fileext = ".pdf")
  f2 <- tempfile(fileext = ".pdf")

  grDevices::pdf(f1)
  expect_no_error(plot_mds_n(result))
  grDevices::dev.off()

  grDevices::pdf(f2)
  expect_no_error(plot_mds_oc(oc_result))
  grDevices::dev.off()

  expect_true(file.exists(f1))
  expect_true(file.exists(f2))

  unlink(c(f1, f2))
})


# ============================================================
# END OF TEST FILE
# ============================================================