Appendix B: Root-Finding Algorithm for Confidence Intervals

confint <- function(FUN,param,stat,lpb=-1.e+4,upb=1.e+4,
                    bound="two-sided",alpha=0.05,tolb=1.e-6,tol=1.e-6,...)
{
  f <- function(x,FUN,param,stat,c,...)
  {
    pars <- function(...) { list(...) }
    args <- pars(...)
    args <- append(args,stat)
    names(args)[length(args)] <- "q"
    if ( is.null(param) == FALSE ) {
      args <- append(args,x)
      names(args)[length(args)] <- param
    }
    do.call(FUN,args) - c
  }
  
  findRoot <- function(f,lpb,upb,tolb,FUN,param,stat,q,tol,...)
  {
    tryCatch(
      {uniroot(f,interval=c(lpb+tolb,upb-tolb),tol=tol,FUN,param,stat,q,...)$root},
      warning = function(w) { cat("The interval bounds are invalid.\n"); return(NULL); },
      error = function(e) { cat("Root not computable.\n"); return(NULL) }
    )
  }

  if ( bound == "two-sided" ) {
    lo <- findRoot(f,lpb,upb,tolb,FUN,param,stat,1-alpha/2,tol,...)
    if ( is.null(lo) ) { return(NULL) }
    hi <- findRoot(f,lpb,upb,tolb,FUN,param,stat,alpha/2,tol,...)
    if ( hi < lo ) { tmp = hi; hi = lo; lo = tmp; }
    return(c(lo,hi))
  } else if ( bound == "lower" ) {
    lo.inc <- findRoot(f,lpb,upb,tolb,FUN,param,stat,1-alpha,tol,...)
    if ( is.null(lo.inc) ) { return(NULL) }
    lo.dec <- findRoot(f,lpb,upb,tolb,FUN,param,stat,alpha,tol,...)
    lo <- min(c(lo.inc,lo.dec))
    return(lo)
  } else if ( bound == "upper" ) {
    hi.inc <- findRoot(f,lpb,upb,tolb,FUN,param,stat,alpha,tol,...)
    if ( is.null(hi.inc) ) { return(NULL) }
    hi.dec <- findRoot(f,lpb,upb,tolb,FUN,param,stat,1-alpha,tol,...)
    hi <- max(c(hi.inc,hi.dec))
    return(hi)
  } else {
    stop("invalid choice for bound")
  }
  return(NULL)
}