Nothing
#' 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.