R/ci_plot.R

In cdfinv: Confidence Interval Estimation via CDF Inversion

#' Plotting of confidence intervals derived via CDF inversion
#'
#' ci_plot() creates a plot showing the cumulative distribution function(s) for the
#' sampling distribution(s) of the chosen statistic, evaluated at the interval bound(s).
#' The horizontal dashed line(s) show(s) the assumed quantile(s)
#' (e.g., 0.95 for a 95\% lower-bound), while the vertical dashed line shows the statistic
#' value.
#'
#' @param cdfinv.out output list from cdfinv()
#' @param ... those additional arguments that were passed to cdfinv()
#'
#' @returns None
#'
#' @author Peter E. Freeman, \email{pfreeman@@cmu.edu}
#'
#' @examples
#' ci_plot(cdfinv("norm","mean",2.5,sd=3),sd=3)
ci_plot <- function(cdfinv.out,...)
{
  PFUN = paste("p",cdfinv.out$DISTR,sep="")
  QFUN = paste("q",cdfinv.out$DISTR,sep="")
  if ( is.character(PFUN) ) {
    PFUN <- get(PFUN,mode="function",envir=parent.frame())
  }
  if ( is.function(PFUN) == FALSE ) {
    stop("The distribution name ",cdfinv.out$DISTR," is not recognized.")
  }
  if ( is.character(cdfinv.out$PARAM) == FALSE ) {
    stop("The name of a PARAM dictating the cdf shape is required.")
  }
  if ( is.numeric(cdfinv.out$STAT) == FALSE ) {
    stop("An observed STAT value is required.")
  }

  pf <- function(FUN,param,bound,q,...)
  {
    pars <- function(...) { list(...) }
    args <- pars(...)
    args <- append(args,q)
    names(args)[length(args)] <- "q"
    if ( is.null(param) == FALSE ) {
      args <- append(args,bound)
      names(args)[length(args)] <- param
    }
    do.call(FUN,args)
  }

  qf <- function(FUN,param,bound,p,...)
  {
    pars <- function(...) { list(...) }
    args <- pars(...)
    args <- append(args,p)
    names(args)[length(args)] <- "p"
    if ( is.null(param) == FALSE ) {
      args <- append(args,bound)
      names(args)[length(args)] <- param
    }
    do.call(FUN,args)
  }

  if ( length(cdfinv.out$q) == 2 ) {
    x.lo       <- min(c(qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],0.001,...),
                        qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[2],0.001,...)))
    x.hi       <- max(c(qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],0.999,...),
                        qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[2],0.999,...)))
    x.plot     <- seq(x.lo,x.hi,by=(x.hi-x.lo)/1000)
    F.x.plot.1 <- rep(NA,length(x.plot))
    F.x.plot.2 <- rep(NA,length(x.plot))
    for ( ii in 1:length(x.plot) ) {
      F.x.plot.1[ii] <- pf(PFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],x.plot[ii],...)
      F.x.plot.2[ii] <- pf(PFUN,cdfinv.out$PARAM,cdfinv.out$bound[2],x.plot[ii],...)
    }
    plot(x.plot,F.x.plot.1,ylim=c(0,1),lwd=1.5,xlab="x",
         ylab=expression(F[X]*"(x)"),typ="l")
    lines(x.plot,F.x.plot.2,lwd=1.5)
    abline(v=cdfinv.out$STAT,lwd=1.5,lty=2)
    abline(h=cdfinv.out$q[1],lwd=1.5,lty=3)
    abline(h=cdfinv.out$q[2],lwd=1.5,lty=3)
  } else {
    x.lo       <- qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],0.001,...)
    x.hi       <- qf(QFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],0.999,...)
    x.plot     <- seq(x.lo,x.hi,by=(x.hi-x.lo)/1000)
    F.x.plot   <- rep(NA,length(x.plot))
    for ( ii in 1:length(x.plot) ) {
      F.x.plot[ii] <- pf(PFUN,cdfinv.out$PARAM,cdfinv.out$bound[1],x.plot[ii],...)
    }
    plot(x.plot,F.x.plot,ylim=c(0,1),lwd=1.5,xlab="x",
         ylab=expression(F[X]*"(x)"),typ="l")
    abline(v=cdfinv.out$STAT,lwd=1.5,lty=2)
    abline(h=cdfinv.out$q[1],lwd=1.5,lty=3)
  }
}