R/cdfinvsim.R

In cdfinv: Confidence Interval Estimation via CDF Inversion

#' Computation of confidence intervals via simulations and CDF inversion
#'
#' cdfinv.sim() returns one- or two-sided confidence interval estimates.
#'
#' @param DISTR name of distribution (in R) from which each datum is sampled
#' @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 seed random number generator seed
#' @param numsim number of simulated datasets
#' @param nsamp sample size for each simulated dataset
#' @param stat.func pointer to function computing the statistic for each dataset
#' @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}
#'
cdfinv.sim <- 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,
                       seed=1,numsim=100000,nsamp=1,
                       stat.func=mean,...)
{
  RFUN = paste("r",DISTR,sep="")
  if ( is.character(RFUN) ) {
    RFUN <- get(RFUN,mode="function",envir=parent.frame())
  }
  if ( is.function(RFUN) == 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.")
  }
  if ( is.null(stat.func) ) {
    stop("Please provide a statistic-computing function.")
  }

  f <- function(x,RFUN,PARAM,STAT,q,seed,numsim,nsamp,stat.func,...)
  {
    pars <- function(...) { list(...) }
    args <- pars(...)
    args <- append(args,numsim*nsamp)
    names(args)[length(args)] <- "n"
    if ( is.null(PARAM) == FALSE ) {
      args <- append(args,x)
      names(args)[length(args)] <- PARAM
    }

    set.seed(seed)
    r <- do.call(RFUN,args)
    Y <- apply(matrix(r,nrow=numsim),1,stat.func)
    sum(Y <= STAT)/numsim - q
  }

  findRoot <- function(f,lpb,upb,tolb,RFUN,PARAM,STAT,q,tol,
                       seed,numsim,nsamp,stat.func,...)
  {
    tryCatch(
      {
        uniroot(f,interval=c(lpb+tolb,upb-tolb),tol=tol,RFUN,PARAM,STAT,
                q,seed,numsim,nsamp,stat.func,...)$root
      },
      warning = function(w) {
        warning("The interval bounds are invalid."); return(NULL);
      },
      error = function(e) {
        stop("The root is not derivable.")
      }
    )
  }

  if ( bound == "two-sided" ) {
    lo <- findRoot(f,lpb,upb,tolb,RFUN,PARAM,STAT,1-alpha/2,tol,
                   seed,numsim,nsamp,stat.func,...)
    if ( is.null(lo) ) { return(NULL) }
    hi <- findRoot(f,lpb,upb,tolb,RFUN,PARAM,STAT,alpha/2,tol,
                   seed,numsim,nsamp,stat.func,...)
    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,RFUN,PARAM,STAT,1-alpha,tol,
                       seed,numsim,nsamp,stat.func,...)
    if ( is.null(lo.inc) ) { return(NULL) }
    lo.dec <- findRoot(f,lpb,upb,tolb,RFUN,PARAM,STAT,alpha,tol,
                       seed,numsim,nsamp,stat.func,...)
    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,RFUN,PARAM,STAT,alpha,tol,
                       seed,numsim,nsamp,stat.func,...)
    if ( is.null(hi.inc) ) { return(NULL) }
    hi.dec <- findRoot(f,lpb,upb,tolb,RFUN,PARAM,STAT,1-alpha,tol,
                       seed,numsim,nsamp,stat.func,...)
    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)
}