Nothing
#' Generic function for predictive intervals
#'
#' See [predictive_interval.stanreg()](https://mc-stan.org/rstanarm/reference/predictive_interval.stanreg.html)
#' in the \pkg{rstanarm} package for an example.
#'
#' @export
#' @template args-object
#' @template args-dots
#' @param prob A number \eqn{p \in (0,1)}{p (0 < p < 1)} indicating the desired
#' probability mass to include in the intervals.
#'
#' @return `predictive_interval()` methods should return a matrix with two
#' columns and as many rows as data points being predicted. For a given value
#' of `prob`, \eqn{p}, the columns correspond to the lower and upper
#' \eqn{100p}\% interval limits and have the names \eqn{100\alpha/2}\% and
#' \eqn{100(1 - \alpha/2)}\%, where \eqn{\alpha = 1-p}. For example, if
#' `prob=0.9` is specified (a \eqn{90}\% interval), then the column names
#' would be `"5%"` and `"95%"`, respectively.
#'
#' The default method just takes `object` to be a matrix and computes
#' quantiles, with `prob` defaulting to `0.9`.
#'
#' @template seealso-rstanarm-pkg
#' @template seealso-vignettes
#'
#' @examples
#' # Default method takes a numeric matrix (of draws from posterior
#' # predictive distribution)
#' ytilde <- matrix(rnorm(100 * 5, sd = 2), 100, 5) # fake draws
#' predictive_interval(ytilde, prob = 0.8)
#'
#' # Also see help("predictive_interval", package = "rstanarm")
#'
predictive_interval <- function(object, ...) {
UseMethod("predictive_interval")
}
#' @rdname predictive_interval
#' @export
predictive_interval.default <- function(object, prob = 0.9, ...) {
if (!is.matrix(object))
stop("For the default method 'object' should be a matrix.")
.central_intervals(object, prob)
}
# internal ----------------------------------------------------------------
#' Compute central intervals
#'
#' @noRd
#' @param object A numeric matrix
#' @param prob Probability mass to include in intervals (in (0,1))
#' @return See above.
#'
#' @importFrom stats quantile
.central_intervals <- function(object, prob) {
stopifnot(is.matrix(object))
if (length(prob) != 1L || prob <= 0 || prob >= 1)
stop("'prob' should be a single number greater than 0 and less than 1.",
call. = FALSE)
alpha <- (1 - prob) / 2
probs <- c(alpha, 1 - alpha)
labs <- paste0(100 * probs, "%")
out <- t(apply(object, 2, quantile, probs = probs))
structure(out, dimnames = list(colnames(object), labs))
}
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.