exp_fixed_duration_tests: Exponential Fixed Duration Test Plans

In jjw3952/mcotear: Contains helpful functions for use within MCOTEA

Exponential Fixed Duration Test Plans

Description

exp_fixed_duration_tests See pg. 67 of MIL-HDBK-781A. Fixed duration tests for exponentially distributed MTBF with desired alpha and beta.

Usage

exp_fixed_duration_tests(mtbf0, mtbfa, alpha = 0.2, beta = 0.2)

Arguments

mtbf0	The MTBF below which to reject the system with probability alpha.
alpha	The type I error rate, producer's risk - the probability we reject a system with MTBF > MTBFA.
beta	The type II error rate, consumer's risk - the probability we accept a system with MTBF < MTBF0.
mtbfA	The MTBF above which to accept the system with probability 1-beta.

Value

The output is a data.frame containing the allowable number of failures at which the system would pass the test (Accept), the number of failures at which the system would fail the test (reject), the test duration (Duration), and the actual beta value (the probability we accept a system with MTBF = MTBF0) given the proposed test.

References

Mil-Hdbk-781A

See Also

exp_mtbf_req, exp_reliability_req, exp_test_duration, exp_mean_lcb, exp_mean_ci, exp_oc

Examples

(fd <- exp_fixed_duration_tests(mtbf0 = 500, mtbfa = 1000, alpha = .2, beta = .2))
theta <- seq(250, 2000, 50)
prob_pass <- exp_oc(accept = fd$Accept[6], duration = fd$Duration[6], mtbf = theta)

df <- data.frame(x = c(500,1000), y = c(0.20, 1-fd$beta[6]), label = c("alpha", "1-beta"))

ggplot(data.frame(prob_pass = prob_pass, theta = theta)) +
  geom_point(ggplot2::aes(x = theta, y = prob_pass)) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0,1,.2)) +
  geom_point(
    data = df,
    mapping = aes(x, y), colour = "red"
  ) +
  geom_text(
    data = df,
    mapping = aes(x, y, label = label), parse = TRUE, colour = "red", hjust = c(-1,1.2)
  )

# Recipe for plotting a single OC Curve
alpha <- .2
beta <- .2
mtbf0 <- 500   # threshold requirement
mtbfa <- 1000  # setting this closer to the requirement will result in
               # requiring a longer test and hence the table returned
               # from exp_fixed_duration_tests will be loner

(fd <- exp_fixed_duration_tests(mtbf0 = mtbf0, mtbfa = mtbfa, alpha = alpha, beta = beta))
test <- 6 # select row of test from fd
theta <- seq(250, 2000, 50)

# The probability of passing the test of a given duration with an allowable number of failures, given the true MTBF
prob_pass <- exp_oc(accept = fd$Accept[test], duration = fd$Duration[test], mtbf = theta)

df <- data.frame(x = c(mtbf0, mtbfa), y = c(alpha, 1-fd$beta[test]), label = c("alpha", "1-beta"))

# The multiples of mtbf0 I want plotted on the top x-axis
 m <- seq(0, 9, 1)
 labels <- parse(text = paste0(m, "*theta"))

ggplot(data.frame(prob_pass = prob_pass, theta = theta)) +
  geom_point(aes(x = theta, y = prob_pass)) +
  geom_path(aes(x = theta, y = prob_pass)) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0,1,.2)) +
  geom_point(
    data = df,
    mapping = aes(x, y), colour = "red"
  ) +
  geom_text(
    data = df,
    mapping = aes(x, y, label = label), parse = TRUE, colour = "red", hjust = c(-1,1.2)
  ) +
   labs(
     x = "MTBF", y = "Prob of Acceptance",
     title = "Exponential OC Curve",
     subtitle = bquote(atop(H[0]*": MTBF"<=.(mtbf0)*phantom(0),H[a]*": MTBF"==.(mtbfa))),
     caption = bquote(
                 atop(
                   "Test Duration = "*.(ceiling(fd$Duration[test]))*". Allowable Failures = "*.(fd$Accept[test])*".",
                   theta==MTBF~Threshold==.(mtbf0)
                 )
               )
   ) +
   scale_x_continuous(
     breaks = seq(0, 3000, 200),
     sec.axis = sec_axis(~.,
       breaks = m*mtbf0,
       labels = labels
     )
   )

# Recipe for plotting multiple OC curves
alpha <- .2
beta <- .2
mtbf0 <- 500   # threshold requirement
mtbfa <- 1000  # setting this closer to the requirement will result in
               # requiring a longer test and hence the table returned
               # from exp_fixed_duration_tests will be loner

(fd <- exp_fixed_duration_tests(mtbf0 = mtbf0, mtbfa = mtbfa, alpha = alpha, beta = beta))
test <- 1:7 # select row(s) of test from fd
theta <- seq(mtbf0/2, mtbf0*8.8, 50) # the range of the true MTBF over which I want to plot

#library(purrr)
#library(reshape2)

# exp_oc gives the probability of passing the test of a given duration with an allowable number of failures, given the true MTBF
prob_pass <- purrr::map2(
  .x = fd$Accept[test],
  .y = fd$Duration[test],
  .f = ~exp_oc(accept = .x, duration = .y, mtbf = theta))
names(prob_pass) <- paste0("OC", 1:length(prob_pass))
prob_pass <- as.data.frame(prob_pass)
prob_pass <- reshape2::melt(prob_pass, id.vars = NULL, value.name = "prob_pass", variable.name = "OC")
prob_pass$Accept <- rep(ceiling(fd$Accept[test]), each = length(theta))
prob_pass$Duration <- rep(ceiling(fd$Duration[test]), each = length(theta))
prob_pass$MTBF <- theta
head(prob_pass)


df <- unique(data.frame(
  x = rep(c(mtbf0, mtbfa), each = length(test)),
  y = c(rep(alpha, length(test)), 1-fd$beta[test]),
  label = rep(c("alpha", "1-beta"), each = length(test))
  ))

exp_test_duration(0, 500)

# The multiples of mtbf0 I want plotted on the top x-axis
m <- seq(0, 9, 1)
labels <- parse(text = paste0(m, "*theta"))

ggplot(
  prob_pass,
  aes(
    x = MTBF,
    y = prob_pass,
    colour = interaction(Accept, Duration, sep = ", ")
    )
  ) +
  geom_hline(yintercept = 1-beta, colour = "red", linetype = "longdash") +
  geom_point(size = .75) +
  geom_path() +
  scale_y_continuous(limits = c(0,1), breaks = seq(0,1,.2)) +
  scale_x_continuous(
    breaks = seq(0, 5000, 500),
    limits = c(0, max(prob_pass$MTBF)),
    sec.axis = sec_axis(~.,
      breaks = m*mtbf0,
      labels = labels
    )
  ) +
  geom_point(
    inherit.aes = FALSE,
    data = df,
    mapping = aes(x, y, fill = label), shape = 21, colour = "transparent"
  ) +
  scale_fill_manual(
    "Measure of Merit",
    values = c("alpha" = "blue", "1-beta" = "red"),
    labels = scales::parse_format()
  ) +
  labs(
    colour = "Allowable Failures,\nTest Duration",
    x = "MTBF", y = "Prob of Acceptance",
    title = "Exponential OC Curve",
    subtitle = bquote(atop(H[0]*": MTBF"<=.(mtbf0)*phantom(0),H[a]*": MTBF"==.(mtbfa))),
    caption = bquote(theta==MTBF~Threshold==.(mtbf0))
  ) +
  guides(
    fill = guide_legend(
      override.aes = list(
        colour = "transparent" #c("blue", "red")
        )
      )
  ) +
  theme(
    # plot.title = element_text(vjust = 1, lineheight = 1, margin = margin(0,0,-35.5,0)),
    # plot.subtitle = element_text(hjust = 1, vjust = 0, lineheight = 1, margin = margin(0,0,5.5,0) ),
    # plot.margin = margin(30,0,0,0)
  )
# getwd()
# ggsave("Exp OC.png", height = 5, width = 6.5)
# grid::grid.ls(grid::grid.force())
# grid::grid.gedit("key-3-1-1.4-2-4-2", size = grid::unit(7, "points"))
# grid::grid.gedit("key-4-1-1.5-2-5-2", size = grid::unit(7, "points"))

# Clean up workspace
rm(list =
  c("alpha", "beta", "df", "fd", "mtbf0", "m", "labels",
    "mtbfa", "prob_pass","test", "theta")
)

jjw3952/mcotear documentation built on Sept. 2, 2023, 10:30 a.m.