R/pbox_imprecise.R

Defines functions pbox pbox.default pbox.summary.imprecise

Documented in pbox pbox.default pbox.summary.imprecise

#' @title Construction of Probability Box
#'
#' @description The function \code{pbox} generates a probability box of 
#' imprecise posterior distributions using the Markov chain produced by 
#' \code{method=MH}.
#'
#' @param x An object of class \code{summary.imprecise} or 
#'          extraced by \code{extractMCHAIN}.
#' @param control control variables
#' @param pretty Logical; Defaults to TRUE; Would you like to smoothe 
#'        the empirical cumulative density curve?  See \emph{Details}.
#' @param ... other arguments to be passed to \code{pbox}.
#' 
#' @details
#' The function \code{ecdf} is applied to the Markov chain at each extrem point.
#' The function \code{density} is used to make the empirical cdf smooth.
#' 
#' If no object that contains a set of Markov chains is found,
#' the \code{update.imprecise} procedure with \code{method=MH} is automatically 
#' taken.
#'
#' @keywords probability box
#' @seealso 
#' \code{\link{density}}, \code{\link{ecdf}}, \code{\link{update.imprecise}}
#'
#' @references
#' Lee (2013) ``Imprecise inferential framework'', PhD thesis.
#'
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @export
pbox <- function(x, ...) UseMethod("pbox")
NULL

#' @rdname pbox
#' @method pbox default
#' @S3method pbox default
pbox.default <- function(x, ...){
  invisible(x)
}
NULL

#' @rdname pbox
#' @method pbox summary.imprecise
#' @S3method pbox summary.imprecise
pbox.summary.imprecise <- function(x, control=list(), 
                                   pretty=FALSE, verbose=FALSE, ...){
                                   
  # naming convention
  xtms <- x$xtms
  method <- x$method 
  xreg <- x$xreg
  y <- x$y
  X <- x$X
  ztrunc <- x$ztrunc
  B <- x$B
  init <- x$init
  inf.idx <- x$inf.idx
  sup.idx <- x$sup.idx
  mc <- match.call()
  xid <- names(x$m1)
  
  # control variable check
  ctl <- list(beta=1, xlim=NULL)
  optimidx <- c(inf.idx[ctl$beta], sup.idx[ctl$beta])
  ctl$xtms <- optimidx
  if (all(names(control) %in% names(ctl))) {
    ctl[names(control)] <- control
  } else {
    stop("Unknown names in ", sQuote("control"), ".")
  }
  
  stopifnot(is.numeric(ctl$xtms), is.vector(ctl$xtms), length(ctl$xtms)<=length(xid))
  ctl$xtms <- ctl$xtms[!duplicated(ctl$xtms)]
  
  if (xreg) {
    stopifnot(length(ctl$beta)==1)
  } else {
    ctl$beta <- NULL
  } 
  if (verbose) {
    message(mc[[1]], ":\nInput sanity check .................... PASS!")
  }
  
  # function definitions
  fn.chain <- function(...){
    xtms <- xtms[xid[ctl$xtms]]
    # op <- lapply(names(xtms), FUN=function(p){
      # xtms.i <- p
      # p <- as.vector(xtms[[p]])
    op <- lapply(xtms, FUN=function(p){
      if(xreg){
        rvals <- cpef2reg(b=p, B=B, y=y, X=X, ztrunc=ztrunc, method="MH", start=init, initrun=TRUE, verbose=FALSE, apriori="normal")$MHchain
        rvals <- as.data.frame(rvals)
      } else {
        rvals <- cpef(hparam=p, y=y, ztrunc=ztrunc, method="MH", apriori=apriori, initrun=TRUE, verbose=FALSE)$MHchain
      }
      return(rvals)
    })
    return(op)
  }
  
  if (all(method=="MH", exists(x="MHchain", where=x))) {
    PAR <- x$MHchain
  } else {
    PAR <- fn.chain(...)
  } 
  
  # print(names(PAR))
  
  if (xreg) {
    PAR1 <- lapply(PAR, "[[", ctl$beta)
    if (verbose) {
      message(sQuote(ctl$beta),"-th regression parameter is selected.")
    }
  } else {
    PAR1 <- PAR
  }
  
  # print(names(PAR1))
  chain1 <- PAR1[xid[ctl$xtms]]
  if (verbose) {
    message("Following extreme points are selected: ", sQuote(xid[ctl$xtms]))
  }
  
  # print(chain1)
  
  colorset <- rainbow(length(names(chain1)))
  
  if (is.null(ctl$xlim)) {
		xlim <- range(chain1)
  } else {
		xlim <- ctl$xlim
  }
  
  plot(0, type='n', xlim=xlim, ylim=c(0,1), ylab=expression(F~"("~theta~"|"~y~")"), xlab=expression(theta), ...)
  abline(h=c(0,1), lty="dashed")
  
  for (i in seq_len(length(chain1))) {
  
    vn <- names(chain1)[i]    
    if (pretty) {
      vec <- density(chain1[[vn]])
      mat <- cbind(vec$x, cumsum(vec$y)/sum(vec$y))
      lines(mat, type='l', lwd=2, col=colorset[i])
    } else {
      env <- environment(ecdf(chain1[[vn]]))
      lines(cbind(env$x, env$y), lwd=2, type='l', col=colorset[i])
    }
  }
  legend("topleft", legend=names(chain1), fill=colorset)
  
  if (verbose) {
    message("Probability-box is produced!")
  }
  invisible(x)
}
NULL

Try the ipeglim package in your browser

Any scripts or data that you put into this service are public.

ipeglim documentation built on May 2, 2019, 4:31 p.m.