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


*

*they simulate a random data with n=20

*the desired shift .15

*the desired ARL0= or close to 370

*they calculate the CUSUM statistics c_i= min(0, x_i - mu_0+ k + c_(i-1))

*the range for choosing h is between 0 and 5 and h1=0 and h2=5

*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. 
 A: 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.
