###  
###  S-plus Code for Developing Acceptance Sampling Plans
###  S P Millard and A Krause (ed), Applied Statistics in the 
###  Pharmaceutical Industry with Case Studies Using S-plus,
###  Springer 2001, Section 19.A.2.
###  
     
read.numeric <- function(text) { 
### A function to read numeric values 
  cat(text) 
  x <- as.numeric(readline()) 
  if(is.na(x)) 
    x <- read.numeric(text) 
  else return(x) 
}
     
AccS <- function() {
### A function for developing acceptance sampling plans 
  spacer <- "\n======================================\n"
  m1 <- 999 
  while(m1 != 6 & m1 != 0) {
## Generate a menu 
    cat(spacer) 
    choices <- c("Single Sampling Plan", 
      "Double Sampling Plan", 
      "Sampling Plan Based on Risk Analysis (SPBORA)", 
      "OC Curves of Single/Double Sampling Plan", 
      "Risk Plot of SPBORA", "Done") 
## Select an item from the menu 
    m1 <- menu(choices, title = "CHOOSE AN ACTIVITY") 
    if (m1 == 1) { 
## Create a single sampling plan 
    cat(spacer) 
    n <- read.numeric("ENTER SAMPLE SIZE N: ")
    c2 <- read.numeric("ENTER ACCEPTANCE NUMBER C: ")   #
## p - possible quality levels 
    p <- c(0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 0.5, 0.7, 
      c(1:15), 20)/100 
    CR <- pbinom(c2, n, p) 
    PR <- 1 - CR 
    DefectRate <- paste(p * 100, 
    rep("%", times = length(p)), sep = "") 
    sp1 <- data.frame("Defect Rate" = DefectRate, 
      PR = round(PR, 3), CR = round(CR, 3)) 
    cat(" \n") 
    cat("Single Sampling Plan\n") 
    cat(" \n") 
    print(sp1) 
  }
  if(m1 == 2) { 
## Create a double sampling plan 
      cat(spacer) 
      n1 <- read.numeric("ENTER INITIAL SAMPLE SIZE N1: ")
      n2 <- read.numeric("ENTER SECOND SAMPLE SIZE N2: ")
      c1 <- read.numeric( 
        "ENTER INITIAL ACCEPTANCE NUMBER C1: ")
      r1 <- read.numeric( 
        "ENTER INITIAL REJECTION NUMBER RI: ") 
      c2 <- read.numeric( 
        "ENTER SECOND ACCEPTANCE NUMBER C2: ")
      p <- c(0.001, 0.005, 0.01, 0.05, 0.1, 0.2, 0.5, 0.7, 
        c(1:15), 20)/100 
      DefectRate <- paste(p * 100, 
        rep("%", times = length(p)), sep ="")
      CR <- pbinom(c1, n1, p) 
      for(i in r1 + 1:n1) { 
        CR <- CR + dbinom(i, n1, p) * pbinom(c2 - i, n2, p) 
      }
      PR <- 1 - CR 
      sp1 <- data.frame("Defect Rate" = DefectRate, 
        PR = round(PR, 3), CR = round(CR, 3)) 
      cat(" \n") 
      cat("Double Sampling Plan\n") 
      cat(" \n") 
      print(sp1)
    } 
    if(m1 == 3) {
      cat(spacer) 
      n.up <- read.numeric( 
        "ENTER THE UPPER LIMIT FOR SAMPLE SIZE N: ")
      L0 <- read.numeric( 
        "ENTER COST OF ACCEPTING A GOOD LOT: ")
      L1 <- read.numeric( 
        "ENTER COST OF ACCEPTING A BAD LOT: ")
      L2 <- read.numeric( 
        "ENTER COST OF REJECTING A GOOD LOT: ")
      L3 <- read.numeric( 
        "ENTER COST OF REJECTING A BAD LOT: ")
      p <- read.numeric("ENTER PREVALENCE OF GOOD LOT: ")
      PQL <- read.numeric( 
        "ENTER PRODUCER'S QUALITY LEVEL: ")
      CQL <- read.numeric( 
        "ENTER CONSUMER'S QUALITY LEVEL: ")
      R <- matrix(c(rep(100000, times = n.up^2  + n.up)), 
        ncol = n.up + 1, byrow = T) 
      for(i in 1:n.up) {
        j <- 1 
        while(j <= i + 1) { 
          CR <- pbinom(j - 1, i, CQL) 
          PR <- 1 - pbinom(j - 1, i, PQL) 
          R[i, j] <- p * L0 * (1 - PR) + 
            (1 - p) * L1 * CR + p * L2 * PR + 
            (1 - p) * L3 * (1 - CR) 
          j <- j + 1 
        }
      }
      j0 <- c(rep(1, times = n.up)) 
      for(i in 1:n.up) {
        j <- 2 
        while(j <= i + 1) {
          if(R[i, j] < min(R[i, i:j - 1]))
            j0[i] <- j
          j <- j + 1 
        }
      }
      r <- c(rep(0, times = n.up)) 
      for(i in 1:n.up) {
        r[i] <- R[i, j0[i]]
      }
      i0 <- 1 
      i <- 2 
      while(i <= n.up) {
        if(r[i] < min(r[1:i - 1])) 
          i0 <- i 
        i <- i + 1
      }
### j0[i0]-1 represents the true acceptance number 
      sp1 <- data.frame("Sample Size" = i0, 
        "Acceptance Number" = j0[i0] - 1, 
        "Minimum Risk" = R[i0, j0[i0]]) 
      cat(" \n") 
      cat("Sampling Plan Based on Risk Analysis\n") 
      cat(" \n") 
      print(sp1) 
      R[R > 99999] <- max(R[R < 99999]) + 1 
      x <- c(1:n.up) 
      y <- c(0:n.up) 
      R <- round(R, dig = 5) 
    }
    if(m1 == 4) {
      xl <- read.numeric( 
        "ENTER LOWER PLOTTING LIMIT FOR X-AXIS: ")
      xu <- read.numeric( 
        "ENTER UPPER PLOTTING LIMIT FOR X~AXIS: ")
      yl <- read.numeric( 
        "ENTER LOWER PLOTTING LIMIT FOR Y-AXIS: ")
      yu <- read.numeric( 
        "ENTER UPPER PLOTTING LIMIT FOR Y-AXIS: ")
      cat("ENTER PLOT TITLE: ")
      pt.title <- readline() 
      plot(p * 100, CR, type = "n", 
        xlab = "Defect Rate (%)", 
        ylab = "Probability of Acceptance", 
        xlim = c(xl, xu), ylim = c(yl,yu), 
        main = pt.title) 
      lines(p * 100, CR) 
      }
    if(m1 == 5) {
      cat("ENTER PLOT TITLE: ") 
      pt.rtitle <- readline() 
      persp(x, y, R, xlab = "Sample Size", 
        ylab = "Acceptance Number", zlab = "Risk Index") 
      title(pt.rtitle) 
    } 
  }
}

accs <- function() { AccS() }