Nothing
#' @aliases summary summary.imprecise
#'
#' @title
#' Summarizing Imprecise Class Objcects
#'
#' @description
#' The function \code{summary.imprecise} is the \code{imprecise} method of
#' generic function \code{summary} which is used for producing summaries of
#' the results from various model fitting functions performed on the imprecise
#' inferential framework.
#' See \sQuote{Details}.
#'
#' @param object
#' the object of class \code{imprecise}.
#'
# @param prob
# a numeric vector of probabilities with values in [0,1].
#'
#' @param HT.est
#' a logical value; Defaults to \code{FALSE}.
#' Would like to see Horvitz-Thompson estimator?
#'
#' @param silent
#' a logical value; Would like to see the report of progress messages
#' (including warnings and error messages) generated from numerical methods
#' that are used for a numerical approximation?
#' Defaults to \code{FALSE}.
#'
#' @details
#' Object summaries depend on the option chosen in the \code{method}.
#' The commonly returned object is \code{est} which is a collection of
#' all Bayes estimates from the imprecise posterior.
#'
#' The Horvitz-Thompson estimator is the implementation of that in the study
#' of \cite{Heijden 2003}.
#' \code{HT.est=TRUE} is not allowed when \code{ztrunc=TRUE}.
#'
#' @return
#' A list with the following compoentns:
#' \item{est}{a list of vectors or matrices of the estimates of canonical
#' parameters.}
#' \item{N}{Horvitz-Thompson estimates if \code{HT.est=TRUE}}
#' \item{mui}{Estimate of the mean parameter if \code{xi} is supplied.}
#' \item{inf}{a numeric vector of length \code{p}, having an infimum of each
#' canonical parameter.}
#' \item{sup}{a numeric vector of length \code{p}, having a supremum of each
#' canonical parameter.}
#' \item{inf.idx}{an index vector corresponding to \code{inf} in \code{est}.}
#' \item{sup.idx}{an index vector corresponding to \code{sup} in \code{est}.}
#' \item{delta}{the vector of imprecision}
#'
#' @references
#' Peter GM van der Heijden, Rami Bustami, Maarten JLF
#' Cruyff, Godfried Engbersen and Hans C van Houwelingen (2003), Point
#' and interval estimation of the population size using the truncated
#' Poisson regression model, \emph{Statistical Modelling}, \bold{3},
#' pp. 305-322.
#'
#' Lee (2013) ``Imprecise inferential framework'', PhD thesis.
#'
#' Walley P (1991). _Statistical reasoning with imprecise probabilities_,
#' series Monographs on statistics and applied probability. Chapman and
#' Hall. ISBN 9780412286605, <URL:
#' http://books.google.ca/books?id=Nk9Qons1kHsC>.
#'
#' @keywords
#' imprecision, Horvitz-Thompson estimator
#'
#' @concept
#' imprecision, prior ignorance, prior modelling, imprecise prior,
#' quantification of prior ignorance, indeterminate decision
#'
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @method summary imprecise
#' @S3method summary imprecise
summary.imprecise <- function(object, HT.est=FALSE, silent=FALSE, ...){
mc <- match.call()
# check and update the stage on the inferential framework
if (object$stage != "update") {
stop("Not correct sequence of imprecise inferential framework.\n",
sQuote(mc[[1]]), " is followed by ",
sQoute("update()"), ".\n")
}
object$stage <- "summary"
# sanity check
stopifnot(!missing(object), is.logical(HT.est))
# stopifnot(is.numeric(prob), is.vector(prob))
# stopifnot(is.numeric(alpha), is.vector(alpha))
stopifnot(is.logical(HT.est), length(HT.est)==1)
stopifnot(is.logical(silent), length(silent)==1)
# naming conventions
xtms <- object$xtms
xreg <- object$xreg
ztrunc <- object$ztrunc
method <- object$method
apriori <- object$apriori
y <- object$y
X <- object$X
m1 <- object$m1
xi <- object$xi
m0shape <- object$m0shape
formula <- object$formula
object$HT.est <- HT.est
if (!ztrunc & HT.est) {
stop("HT.est=", sQuote(HT.est), "cannot be used with ztrunc=",
sQuote(ztrunc))
}
if (!silent) {
message("Input sanity check .................... PASS! \n",
"in ", mc[[1]], ".\n")
}
# definitions for estimating 'mui' and 'N'
fnxregmui <- function(x, ...){
# Computing the value of a linear predictor 'mu_i'
#
# Args:
# x: The estimate of 'beta'
#
# Return:
# The evaluated scalar 'mu_i'
beta <- as.vector(x)
val <- exp(xi%*%beta)
return(val)
}
fnxregN <- function(x, ...){
# Computing the Horvitz-Thompson estimate of 'N' in the context of
# a generalized linear model
#
# Args:
# x: The estimate of 'beta'
#
# Return:
# The evaluated scalar 'N'
beta <- as.vector(x)
mu <- exp(X%*%beta)
val <- sum(1/ppois(q=0, lambda=mu, lower.tail=FALSE, log.p=FALSE))
return(val)
}
fnN <- function(t, ...) {
# Computing the Horvtiz-Thompson estimate of 'N' in the context of
# models that do not involve with covariates
#
# Args:
# t: The estimate of a canonical parameter 'theta'
#
# Return:
# The evaluated scalar 'N'
val <- length(y)/ppois(q=0, lambda=exp(t), lower.tail=FALSE, log.p=FALSE)
return(val)
}
# tracking the correct extreme point
xid <- do.call(c, lapply(m1, "[[", "xid"))
# nxid <- length(xid)
# General summarization for all numerical methods
if (xreg) { # 'nrow:xtms' x 'p:cfs'
est <- as.data.frame(do.call(rbind, lapply(m1, "[[", "cfs")))
inf <- do.call(c, lapply(est, min))
sup <- do.call(c, lapply(est, max))
inf.idx <- do.call(c, lapply(est, which.min))
sup.idx <- do.call(c, lapply(est, which.max))
tslpars <- est
} else { # 'length:xtms' x '1:theta'
est <- do.call(c, lapply(m1, "[[", "theta"))
names(est) <- xid
inf <- min(est)
sup <- max(est)
inf.idx <- which.min(est)
sup.idx <- which.max(est)
tslpars <- do.call(rbind, lapply(m1, "[[", "eta"))
}
delta <- abs(sup-inf)
object <- c(object, list(est=est, inf=inf, sup=sup, inf.idx=inf.idx,
sup.idx=sup.idx, delta=delta, tslpars=tslpars))
if (!is.null(object$qcLA)) {
message("Minimal summaries are kept for imprecise inferential framework.")
invisible(object)
}
# additional summaries for 'MH'
if (method == "MH" & is.null(object$qcLA)) {
if (xreg) {
MHchain <- lapply(lapply(m1, "[[", "MHchain"), as.data.frame)
} else {
MHchain <- as.data.frame(do.call(cbind, lapply(m1, "[[", "MHchain")))
}
accept.rate <- do.call(c, lapply(m1, "[[", "accept.rate"))
object$MHchain <- MHchain
object$accept.rate <- accept.rate
if (ztrunc) {
Nchain <- muichain <- list()
for (i in xid) {
if (!silent) {
message(xid[i], " is processing now ...")
}
mci <- as.matrix(MHchain[[i]])
if (xreg) {
Nchain[[i]] <- apply(mci, 1, fnxregN)
muichain[[i]] <- apply(mci, 1, fnxregmui)
} else {
Nchain[[i]] <- fnN(mci)
muichain[[i]] <- exp(mci)
}
}
N <- do.call(c, lapply(Nchain, mean))
mui <- do.call(c, lapply(muichain, mean))
object$N <- N
object$mui <- mui
object$muichain <- muichain
object$Nchain <- Nchain
}
}
# additional summaries for 'IS'
if(method == "IS" & is.null(object$qcLA)){
if (xreg) {
isample <- lapply(m1, "[[", "isample")
lweight <- lapply(m1, "[[", "lweight")
} else {
isample <- as.data.frame(do.call(cbind, lapply(m1, "[[", "isample")))
lweight <- as.data.frame(do.call(cbind, lapply(m1, "[[", "lweight")))
}
ess <- do.call(c, lapply(m1, "[[", "ess"))
object$isample <- isample
object$lweight <- lweight
object$ess <- ess
if (ztrunc) {
Nsmp <- muismp <- list()
for(i in xid) {
if (!silent) {
message(xid[i], "is processing now ...")
}
isi <- as.matrix(isample[[i]])
lwi <- as.vector(lweight[[i]])
if (xreg) {
Nsmp[[i]] <- apply(isi, 1, fnxregN)
muismp[[i]] <- apply(isi, 1, fnxregmui)
} else {
Nsmp[[i]] <- fnN(isi)
muismp[[i]] <- isi*exp(lwi)/sum(exp(lwi))
}
}
N <- do.call(c, lapply(Nsmp, function(x){
val <- sum(x*exp(lwi)/sum(exp(lwi)))
return(val)
}))
mui <- do.call(c, lapply(muismp, function(x){
val <- sum(x*exp(lwi)/sum(exp(lwi)))
return(val)
}))
object$N <- N
object$mui <- mui
object$Nsmp <- Nsmp
object$muismp <- muismp
}
}
# additional summaries for 'LA'
if (method=="LA" & is.null(object$qcLA)) {
sratio <- do.call(rbind, lapply(m1, "[[", "sratio"))
fdiff <- do.call(rbind, lapply(m1, "[[", "fdiff"))
niter <- do.call(rbind, lapply(m1, "[[", "niter"))
colnames(niter) <- colnames(fdiff) <- colnames(sratio) <- colnames(X)
object <- c(object, list(sratio=sratio, fdiff=fdiff, niter=niter))
if (!is.null(xi)) {
if (xreg) {
object$mui <- do.call(c, lapply(m1, "[[", "mui"))
} else {
object$mui <- NA # TODO
}
}
if (ztrunc & HT.est) {
if (xreg) {
object$N <- do.call(c, lapply(m1, "[[", "N"))
} else {
object$N <- fn(est)
}
}
}
class(object) <- c("summary.imprecise")
return(object)
}
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.