Nothing
#' @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
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.