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