R/cdfinv.R

In cdfinv: Confidence Interval Estimation via CDF Inversion

#' Computation of confidence intervals via CDF inversion
#'
#' cdfinv() returns one- or two-sided confidence interval estimates.
#'
#' @param DISTR name of sampling distribution in R
#' @param PARAM name of distribution parameter for which we are computing an interval estimate
#' @param STAT observed value of the chosen statistic
#' @param lpb lower bound of search interval
#' @param upb upper bound of search interval
#' @param bound one of "two-sided", "lower", or "upper"
#' @param alpha the confidence coefficient is 1 - alpha
#' @param tolb search interval bound offset value
#' @param tol convergence tolerance for uniroot function
#' @param ... additional arguments for DISTR's cdf function
#'
#' @returns A list with interval bounds and associated cdf values.
#' \itemize{
#'   \item DISTR - The distribution name (as given in R)
#'   \item PARAM - The parameter name (as given in R)
#'   \item STAT - The observed statistic value
#'   \item bound - The interval bound(s)
#'   \item q - The cdf quantile(s) associated with the interval bound(s)
#' }
#'
#' @author Peter E. Freeman, \email{pfreeman@@cmu.edu}
#'
#' @examples
#' cdfinv("norm","mean",3.45,sd=2) ## returns -0.4699279 and 7.3699277
#' cdfinv("gamma","rate",12.25,lpb=0,bound="upper",shape=10) ## returns 1.282058
#' cdfinv("nbinom","prob",22,lpb=0,upb=1,bound="lower",size=10) ## returns 0.1803843
cdfinv <- function(DISTR,PARAM,STAT,lpb=-1.e+4,upb=1.e+4,
                   bound="two-sided",alpha=0.05,tolb=1.e-6,tol=1.e-6,...)
{
  PFUN = paste("p",DISTR,sep="")
  if ( is.character(PFUN) ) {
    PFUN <- get(PFUN,mode="function",envir=parent.frame())
  }
  if ( is.function(PFUN) == FALSE ) {
    stop("The distribution name ",DISTR," is not recognized.")
  }
  if ( is.character(PARAM) == FALSE ) {
    stop("The name of a param dictating the cdf shape is required.")
  }
  if ( is.numeric(STAT) == FALSE ) {
    stop("An observed statistic value is required.")
  }

  f <- function(x,PFUN,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(PFUN,args) - c
  }

  findRoot <- function(f,lpb,upb,tolb,PFUN,PARAM,STAT,q,tol,...)
  {
    tryCatch(
      {uniroot(f,interval=c(lpb+tolb,upb-tolb),tol=tol,PFUN,PARAM,STAT,
               q,...)$root},
      warning = function(w) { warning("The interval bounds are invalid."); return(NULL); },
      error = function(e) { stop("The root cannot be derived.") }
    )
  }

  if ( bound == "two-sided" ) {
    lo <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,1-alpha/2,tol,...)
    if ( is.null(lo) ) { return(NULL) }
    hi <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,alpha/2,tol,...)
    o <- order(c(lo,hi))
    return(list(DISTR=DISTR,PARAM=PARAM,STAT=STAT,bound=c(lo,hi)[o],q=c(1-alpha/2,alpha/2)[o]))
  } else if ( bound == "lower" ) {
    lo.inc <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,1-alpha,tol,...)
    if ( is.null(lo.inc) ) { return(NULL) }
    lo.dec <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,alpha,tol,...)
    o <- order(c(lo.inc,lo.dec))
    return(list(DISTR=DISTR,PARAM=PARAM,STAT=STAT,bound=c(lo.inc,lo.dec)[o[1]],q=c(1-alpha,alpha)[o[1]]))
  } else if ( bound == "upper" ) {
    hi.inc <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,alpha,tol,...)
    if ( is.null(hi.inc) ) { return(NULL) }
    hi.dec <- findRoot(f,lpb,upb,tolb,PFUN,PARAM,STAT,1-alpha,tol,...)
    o <- order(c(hi.inc,hi.dec))
    return(list(DISTR=DISTR,PARAM=PARAM,STAT=STAT,bound=c(hi.inc,hi.dec)[o[2]],q=c(alpha,1-alpha)[o[2]]))
  } else {
    stop("Invalid choice for bound.")
  }
  return(NULL)
}