# How find controls limits for CUSUM chart

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.

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)

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.

• 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)

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)
col=c("Black","Blue","Red"),
pch=c(1,18,18),
lty=c(1,1,1))


The result is this plot: 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:

## First to attempt to answer the question given the criteria set

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, 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;
}
}

// [[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.