1
$\begingroup$

I am trying to re-do a simulation work in quality control. In this simulation they use bisection search algorithm to find the control limits h (the threshold) that gives the desire average run length ARL = 370. This is the information they give

  1. they simulate a random data with n=20
  2. the desired shift .15
  3. the desired ARL0= or close to 370
  4. they calculate the CUSUM statistics c_i= min(0, x_i - mu_0+ k + c_(i-1))
  5. the range for choosing h is between 0 and 5 and h1=0 and h2=5
  6. use the bisection search algorithm (in this step they just give the number of steps which the algorithm should run =30 and the replication =15000 and the midpoint: h0=(h1+h2)/2 ).

I never use the bisection algorithm and if there is another way to find h that give the desire ARL i guess it should be fine.

$\endgroup$
5
$\begingroup$

First thoughts:

It looks like there are several packages for CUSUM charts in R, inculding but not limited to the following:

  • 'qcc' (v2.6) - quality control charts (2014)
  • 'spc' (v0.5.1) - statistical process control (2015)
  • 'strucchange' (v1.5-1) - testing structural change in linear regression models (2015)
  • 'surveillance' (v1.9-1) - spatio-temporal modeling of epidemic phenomena (2015)
  • 'spcadjust'(v0.1-1) - early revision of bootstrap based adjust control charts

The 'qcc' function form is:

cusum(data, sizes, center, std.dev, decision.interval = 5, se.shift = 1,
      data.name, labels, newdata, newsizes, newlabels, plot = TRUE, ...)

It is not the library that you are looking for.

The 'spc' function has "xcusum.crit.L0L1' with functional form:

xcusum.crit.L0L1(L0, L1, hs=0, sided="one", r=30, L1.eps=1e-6, k.eps=1e-8)

You likely will find this (or a nearby function from the same library) can do what you are looking for.

Note: So I have a problem with your statement: "they simulate a random data with n=20".

There are an infinite number of breeds/species of random data, within which there are infinite number of distributions, and within which there are somewhere between many and infinite number of samples that could be drawn. This is amazingly general. It is like describing the creation of yogurt as requiring 'some organism' to transform the milk. In reality there is only one distribution in mind, just like there is only one type of organism that transforms milk into yogurt - and it isn't an elephant.

Comments:

  • I like how a lag-plot shows transitions. The standard 4-plot is just very useful when you are first looking at your data. It is called a fundamental, a "Gross Reality Check" (aka GRC) for a good reason.

Exercise:

Here is the basic code (no library) for the NIST example.

#library
library(stats)

#data
mydata <- c(324.925, 324.675, 324.725, 324.350, 325.350, 325.225, 324.125, 324.525, 325.225, 324.600, 324.625, 325.150, 328.325, 327.250, 327.825, 328.500, 326.675, 327.775, 326.875, 328.350)

#hand-made model
n <- length(mydata)
n

#parameters
mu <- 325
h <- 4.1959
k <- 0.3175

#predeclarations
group <- numeric(length = n)
x <- numeric(length = n)

x_less_mu <- numeric(length = n)
x_less_mu_less_k <- numeric(length = n)
S_hi <- numeric(length = n)

mu_less_k_less_x <- numeric(length = n)
S_lo <- numeric(length = n)

my_cusum <- numeric(length = n)

for (i in 1:n){

     group[i]  <- i
     x[i]         <- mydata[i]

     x_less_mu[i] <- x[i] - mu
     x_less_mu_less_k[i] <- x[i]-mu-k
     mu_less_k_less_x[i] <- mu-k-x[i]

     if (i==1){
          S_hi[i]        <- max(c(0,x_less_mu_less_k[i]))
          S_lo[i]        <- max(c(0,mu_less_k_less_x[i]))
          my_cusum[i]    <- x_less_mu[i]
     } else {
          S_hi[i]        <- max(c(0,S_hi[i-1]+x_less_mu_less_k[i]))
          S_lo[i]        <- max(c(0,S_lo[i-1]+mu_less_k_less_x[i]))
          my_cusum[i]    <- x_less_mu[i] + my_cusum[i-1]
     }

}

plot(my_cusum,ylim=c(-5,20),xlab="Index",ylab="Cumulative Sum",type="b")
grid()
lines(S_hi,col="Blue")
points(S_hi,col="Blue",pch=18)
lines(S_lo,col="Red")
points(S_lo,col="Red",pch=18)
legend(x=1,y=20,c("CUSUM","Upper_alert","Lower_alert"),
       col=c("Black","Blue","Red"),
       pch=c(1,18,18),
       lty=c(1,1,1))

The result is this plot:
enter image description here

While it is a pretty picture and shows that a CUSUM might work, it also has "prepackaged" parameters. It is not (yet) reproducible. The question being asked (how to find ARL, h, k, ...) must be answered.

Some useful references:

  1. http://www.itl.nist.gov/div898/handbook/pmc/section3/pmc323.htm
  2. https://stats.stackexchange.com/questions/107562/multivariate-control-chart?rq=1
$\endgroup$
0
$\begingroup$

First to attempt to answer the question given the criteria set

I am not entirely sure about the comments

the desired shift .15

the number of steps which the algorithm should run =30

However if the CUSUM is set for an odds ratio of 1.5, the baseline rate of an in-contol process is 0.15 and the random data drawn from a binomial distribution, the R code could look something like this:

bisect.search <- function(delegate, target, low, high, tolerance=0.005, max.iterations=25){
  for(i in 1:max.iterations){
    half <- (low + high)/2
    val <- delegate(half)
    if (!is.numeric(val)) {
      warning("non numeric value returned - aborting")
      return(NA)
    }
    if ((abs(val - target) <= tolerance)){
      return(half)
    }
    else if (val < target) {
      low <- half
    } else {
      high <- half
    }
  } 
  warning(paste("iteration limit exceeded - current low", low, ", high", high, ", current f(x)", val))
  return(NA)
}

cusum <- function(observed, predicted, odds.ratio, H) {
  # note predicted can be a scalar or vector of the same length as observed 
  si <- log(ifelse(observed, odds.ratio, 1)/(1 - predicted + (odds.ratio * predicted)))
  # sprt <- cumsum(si)
  prev <- 0
  result <- vapply(si, function(j) {
    prev <<- prev + j
    if (prev < 0) {
      prev <<- 0
    } else if (prev >= H) {
      current <- prev
      prev <<- 0 # reset
      return(current)
    }
    return(prev)
  },numeric(1))
  return(result)
}

ARL0 <- function(H, in.control.rate = 0.15, odds.ratio = 1.5, n.trials = 15000 * 370) {
  rnd.results <- rbinom(n.trials, 1, in.control.rate)
  result <- cusum(rnd.results, in.control.rate, odds.ratio, H)
  exceed.h <- which(result >= H)
  if (length(exceed.h) == 0) {
    warning("Never exceeded threshold, try again with lower threshold (H) or higher number of trials")
    return(Inf)
  }
  run.lengths <- c(exceed.h[1], diff(exceed.h))
  return(mean(run.lengths))
}

print(bisect.search(ARL0, target = 370, low = 0, high = 5, tolerance = 2))

In truth most of the time this would be written in RCPP for performance (note the 10 times increase in the number of replications in the example below):

#include <Rcpp.h>
using namespace Rcpp;

double ARL0(double H, double p, int replications, double OR) {
  double si0(log(1 / (1 - p + (OR * p))));
  double si1(log(OR / (1 - p + (OR * p))));
  // we could fill a vector and figure out inter-quartile range etc with std::vector<int> run(replication);
  // and then figure out each run length with int drawNo(0);
  // however adding together the total draws to reach replication will give us the same answer as adding together run lengths
  unsigned long long total(0); // at very large replications and very tolerant H values, this may overflow
  double current(0.0);
  for(int i = 0; i < replications; ++i) {
    while (current < H) {
      ++total;
      current += R::runif(0,1) < p ? si1 : si0;
      if (current < 0.0) {
        current = 0.0;
      }
    }
    current = 0.0;
  }
  return total / replications;
}

// [[Rcpp::export]]
double findH(double p, double targetRL, int replications, double OR = 1.5) {
  double low(0.0);
  double high(5.0);
  const double tolerance(2.0);
  for(int i = 0; i < 25; ++i){
    double half((low + high) / 2.0);
    double val(ARL0(half, p, replications, OR));
    if ((abs(val - targetRL) <= tolerance)){
      return half;
    }
    else if (val < targetRL) {
      low = half;
    }
    else {
      high = half;
    }
  }
  return -1.0;
}

/*** R
findH(0.15, 370, 150000)
*/

if you wanted to be really clever, you would alter the above code to increase the number of replications as the bisection search algorithm honed in closer to the desired run length. Both the above should be giving an H (threshold) value around 1.82

Then to give a more real world approach

The VLAD package for R includes the ability to perform CUSUM control limit calculations based on fast and accurate Markov chain approximations, and alternatively to run Monte-Carlo trials with random numbers generated from beta distributions etc. Different control limit generation techniques can be applied and multi-threading is used to speed up processing time. The peer reviewed references given in the package read-me (co-authored by the R package author) explore the use of the different techniques on offer.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.