R/UnivariateObtainAverageRunLength.R

Defines functions getARL getRL

Documented in getARL getRL

#' @title Run Length
#' @description Get the run length
#' @inheritParams getDist
#' @inheritParams NS
#' @param replica scalar. It is used for the parallel version of the function (\code{parallel=TRUE}). Default \code{1}.
#' @param n scalar. Subroup size
#' @param m scalar. Reference sample size
#' @param mu vector. Two elements, the first one is the mean of the reference sample and the second one is the mean of the monitoring sample.
#' @param sigma vector. Two elements, the first one is the sd of the reference sample and the second one is the sd of the monitoring sample.
#' @param dist.par vector. Distribution parameters. \code{c(par.a, par.b)}. Default \code{c(0,1)}.
#' @param chart character string. Selected type of chart. Three options are available: Shewhart, CUSUM, EWMA
#' @param chart.par vector. The size depends on the selected chart:
#' \describe{
#'   \item{Shewhart scheme: }{is \code{c(k)}, where \code{k} comes from \eqn{UCL = mu + k\sigma, LCL = mu - k\sigma.}}
#'   \item{CUSUM scheme: }{is \code{c(k, h, t)} where \code{k} is the reference value and \code{h} is the control limit,
#'   and \code{t} is the type of the chart (1:positive, 2:negative, 3:two sides)}
#'   \item{EWMA scheme: }{is \code{c(lambda, L)}, where \code{lambda} is the smoothing constant
#'   and \code{L} multiplies standard deviation to get the control limit}
#' }
#' @param calibrate logical. If \code{TRUE} the RL is limit to 10 times the target ARL.
#' @param arl0 scalar. Expected value of the RL. Default \code{370}.
#' @param isFixed logical. If \code{TRUE} the reference sample does not update, otherwise the reference sample is updated whenever the batch is in control.
#' @return \code{RL} vector. The run length of the chart for the parameter setting.
#' @export
#' @import stats
#' @examples
#' n <- 5 # subgroup size
#' m <- 100 # reference-sample size
#' dist <- "Normal"
#' mu <- c(0, 0) # c(reference sample mean, monitoring sample mean)
#' sigma <- c(1, 1) # c(reference sample sd, monitoring sample sd)
#'
#' #### Distribution parameters
#' dist.par <- c(0, 1, 1) # c(location, scale, shape)
#'
#' #### Other Parameters
#' replicates <- 2
#' print.RL <- TRUE
#' calibrate <- FALSE
#' progress <- TRUE
#' arl0 <- 370
#'
#' #### Control chart parameters
#' chart <- "Shewhart"
#' chart.par <- c(3)
#' shewhart <- getRL(1, n, m,
#'   theta = NULL, Ftheta = NULL,dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, calibrate = calibrate, arl0 = arl0
#' )
#'
#' chart <- "CUSUM"
#' chart.par <- c(0.25, 4.4181, 3)
#' cusum <- getRL(1, n, m,
#'   theta = NULL, Ftheta = NULL, dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, calibrate = calibrate, arl0 = arl0
#' )
#'
#' chart <- "EWMA"
#' chart.par <- c(0.2, 2.962)
#' shewhart <- getRL(1, n, m,
#'   theta = NULL, Ftheta = NULL,dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, calibrate = calibrate, arl0 = arl0
#' )
getRL <- function(replica = 1, n, m, theta = NULL, Ftheta = NULL,
                  dist, mu, sigma, dist.par = c(0,1,1), scoring = "Z",
                  chart, chart.par, calibrate = FALSE, arl0 = 370,
                  alignment = "unadjusted", constant = NULL, absolute=FALSE,
                  isFixed=FALSE,Chi2corrector="None", rounding.factor = NULL) {
  # initilize the reference sample
  Y <- NULL
  if (m > 0) { # if there are reference sample
    # generate the reference sample
    Y <- SNSchart::getDist(n = m, dist = dist, mu = mu[1], sigma = sigma[1], dist.par = dist.par, rounding.factor = rounding.factor)

  }

  RL <- 0
  in.Control <- TRUE

  switch(chart,
    Shewhart = {
      k <- chart.par[1]
    },
    CUSUM = {
      #type is always the last value in vector
      type = chart.par[length(chart.par)]

      if(type == 3 && scoring == "Z-SQ"){
        kp <- chart.par[1]
        km <- chart.par[2]
        h <- chart.par[3]
      }else{
        k <- chart.par[1]
        h <- chart.par[2]

        kp <- k
        km <- k
      }

      Cplus <- 0
      Cminus <- 0
    },
    EWMA = {
      lambda <- chart.par[1]
      L <- chart.par[2]
      E <- 0
    }
  )

  while (in.Control) {
    # add one iteration to run length
    RL <- RL + 1

    # generate the subgroup to monitor
    X <- SNSchart::getDist(n = n, dist = dist, mu = mu[2], sigma = sigma[2], dist.par = dist.par, rounding.factor = rounding.factor)

    # get the normal scores
    ns <- SNSchart::NS(X = X, Y = Y, theta = theta, Ftheta = Ftheta, alignment = alignment, constant = constant, scoring = scoring, Chi2corrector=Chi2corrector)
    Z <- ns$Z

    switch(scoring,
           "Z" = {# it is a vector with a subgroup size so it is needed to average them
             Z = mean(Z)
           },
           "Z-SQ" = {# it is a vector with a subgroup size so it is needed to sum them
             Z = sum(Z)
           }
    )

    # if the subgroup is out of the limits
    # an alarm is detected
    switch(chart,
      Shewhart = {
        # if the subgroup is out of the limits an alarm is detected
        ucl = k
        if (scoring == "Z"){
          ucl = ucl / sqrt(n) #k / sqrt(n)
        }
        if (abs(Z) >= ucl) in.Control <- FALSE
      },
      CUSUM = {
        if(scoring == "Z-SQ"){# for obtain variance Z/n
          Z = Z /(n * sqrt(n))
        }
        switch(type,
          "1" = {
            Cplus <- max(c(0, Cplus + Z * sqrt(n) - kp))
          },
          "2" = {
            Cminus <- min(c(0, Cminus + Z * sqrt(n) + km))
          },
          "3" = {
            Cplus <- max(c(0, Cplus + Z * sqrt(n) - kp))
            Cminus <- min(c(0, Cminus + Z * sqrt(n) + km))
          }
        )
        if (Cplus >= h || Cminus <= -h) in.Control <- FALSE
      },
      EWMA = {
        E <- lambda * Z + (1 - lambda) * E

        UCL <- L / sqrt(n) * sqrt(lambda / (2 - lambda) * (1 - (1 - lambda)^(2 * RL)))
        # LCL = - UCL

        if (abs(E) >= UCL) in.Control <- FALSE
      }
    )

    if (calibrate){
      if (RL >= arl0){
        in.Control <- FALSE
      }
    }
    if (RL >= arl0 * 1000){
      in.Control <- FALSE
    }

    # update the reference sample
    if(!isFixed){#if the reference sample is updated (not fixed)
      Y <- c(Y, X)
    }
  }
  return(RL)
}

#' @title Average Run Length (ARL)
#' @description Get the ARL \code{\link{getRL}}
#' @inheritParams getRL
#' @param print.RL logical. If \code{TRUE} return the vectors of RL for each iteration.
#' @param replicates scalar. Number of replicates to get the ARL
#' @param progress logical. If \code{TRUE} it shows the progress in the console.
#' @param isParallel logical. If \code{TRUE} the code runs in parallel according to the
#' number of cores in the computer,otherwise the code runs sequentially. Default \code{TRUE}.
#' @return Multiple output. Select by \code{output$}
#' \itemize{
#'   \item \code{ARL}: scalar. Average Run Length for the \code{RL}s of all the \code{replicates}.
#'   \item \code{SDRL}: scalar. Standard Deviation Run Length for the \code{RL} in all the \code{replicates}.
#'   \item \code{MRL}: bolean. Median Run Length for the \code{RL}s of all the \code{replicates}.
#'   \item \code{QRL}: vector. It retrieve the quantiles (0.05, 0.1, 0.2, 0.25, 0.5, 0.75, 0.8, 0.9, 0.95) for all the \code{RL}s.
#' }
#' @export
#' @import parallel
#' @import stats
#' @examples
#' n <- 5 # subgroup size
#' m <- 100 # reference-sample size
#' dist <- "Normal"
#' mu <- c(0, 0) # c(reference sample mean, monitoring sample mean)
#' sigma <- c(1, 1) # c(reference sample sd, monitoring sample sd)
#'
#' #### Normal distribution parameters
#' dist.par <- c(0, 1) # c(location, scale)
#'
#' #### Other Parameters
#' replicates <- 2
#' print.RL <- TRUE
#' isParallel <- FALSE
#' calibrate <- FALSE
#' progress <- TRUE
#' arl0 <- 370
#'
#' #### Control chart parameters
#' chart <- "Shewhart"
#' chart.par <- c(3)
#' shewhart <- getARL(n, m,
#'   theta = NULL, Ftheta = NULL, dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, print.RL = print.RL,
#'   replicates = replicates, isParallel = isParallel,
#'   calibrate = calibrate, arl0 = arl0
#' )
#'
#' chart <- "CUSUM"
#' chart.par <- c(0.25, 4.4181, 3)
#' cusum <- getARL(n, m,
#'   theta = NULL, Ftheta = NULL, dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, print.RL = print.RL,
#'   replicates = replicates, isParallel = isParallel,
#'   calibrate = calibrate, arl0 = arl0
#' )
#'
#' chart <- "EWMA"
#' chart.par <- c(0.2, 2.962)
#' shewhart <- getARL(n, m,
#'   theta = NULL, Ftheta = NULL, dist, mu, sigma, dist.par = dist.par,
#'   chart = chart, chart.par = chart.par, print.RL = print.RL,
#'   replicates = replicates, isParallel = isParallel,
#'   calibrate = calibrate, arl0 = arl0
#' )
getARL <- function(n, m, theta = NULL, Ftheta = NULL,
                   dist, mu, sigma, dist.par = c(0, 1, 1),
                   chart, chart.par, scoring = "Z",Chi2corrector="None",
                   replicates = 10000, isParallel = TRUE,
                   print.RL = FALSE, progress = FALSE,
                   calibrate = FALSE, arl0 = 370,
                   alignment = "unadjusted", constant = NULL, absolute=FALSE,
                   isFixed=FALSE, rounding.factor = NULL) {

  type = chart.par[length(chart.par)]
  if(type == 3 && scoring == "Z-SQ"){
    if(length(chart.par) < 4){
      stop("Missing argument in chart.par. Four arguments needed.")
    }
  }
  RLs <- NULL
  if (isParallel) {
    cluster <- parallel::makeCluster(parallel::detectCores() - 1)
    parallel::clusterExport(cluster, "NS")
    parallel::clusterExport(cluster, "getDist")
    parallel::clusterExport(cluster, "getRL")
    RLs <- parallel::parSapply(cluster, 1:replicates, getRL, n = n, m = m, theta = theta, Ftheta = Ftheta, dist = dist, mu = mu, sigma = sigma, dist.par = dist.par, chart = chart, chart.par = chart.par, calibrate = calibrate, arl0 = arl0, alignment=alignment, constant=constant,absolute=absolute,isFixed=isFixed,scoring=scoring,Chi2corrector=Chi2corrector, rounding.factor = rounding.factor)
    parallel::stopCluster(cluster)
  } else {
    t0 <- Sys.time()
    for (r in 1:replicates) {
      RL <- SNSchart::getRL(1, n = n, m = m, theta = theta, Ftheta = Ftheta, dist = dist, mu = mu, sigma = sigma, dist.par = dist.par, chart = chart, chart.par = chart.par, calibrate = calibrate, arl0 = arl0, alignment=alignment, constant=constant,absolute=absolute,isFixed=isFixed,scoring=scoring,Chi2corrector=Chi2corrector, rounding.factor = rounding.factor)

      RLs <- c(RLs, RL)

      # print out progress
      if (progress) { # if is TRUE
        if (r %% 10 == 0) { # every 10 replicates
          t1 <- Sys.time()
          remaining.iterations <- replicates - r
          remaining.time <- remaining.iterations * difftime(t1, t0, units = "min") / r
          message("ARL", round(mean(RLs), digits = 1), "-- SDRL", round(sd(RLs), digits = 1), "--> Time remaining", remaining.time, "in minutes to complete", remaining.iterations, "iterations", "\n", sep = " ")
        }
      }
    }
  }

  output <- list(
    ARL = mean(RLs),
    SDRL = sd(RLs),
    MRL = median(RLs),
    QRL = quantile(x = RLs, probs = c(0.05, 0.1, 0.2, 0.25, 0.5, 0.75, 0.8, 0.9, 0.95), names = TRUE, type = 3)
  )
  if (print.RL) output$RL <- RLs

  if (progress) message("Final ARL", round(mean(RLs), digits = 1), "-- SDRL", round(sd(RLs), digits = 1), "\n", "See output variable for more.\n\n", sep = " ")

  return(output)
}

Try the SNSchart package in your browser

Any scripts or data that you put into this service are public.

SNSchart documentation built on April 7, 2021, 9:07 a.m.