R/preddist.R

In RBesT: R Bayesian Evidence Synthesis Tools

#' Predictive Distributions for Mixture Distributions
#'
#' @description
#' Predictive distribution for mixture of conjugate distributions
#' (beta, normal, gamma).
#'
#' @param mix mixture distribution
#' @param n predictive sample size, set by default to 1
#' @param ... includes arguments which depend on the specific prior-likelihood pair, see description below.
#'
#' @details Given a mixture density (either a posterior or a prior)
#'
#' \deqn{p(\theta,\mathbf{w},\mathbf{a},\mathbf{b})}{p(\theta,w,a,b)}
#'
#' and a data likelihood of
#'
#' \deqn{y|\theta \sim f(y|\theta),}{y|\theta ~ f(y|\theta),}
#'
#' the predictive distribution of a one-dimensional summary \eqn{y_n}
#' of $n$ future observations is distributed as
#'
#' \deqn{y_n \sim \int p(\theta,\mathbf{w},\mathbf{a},\mathbf{b}) \, f(y_n|\theta) \, d\theta .}{y_n ~ \int p(u,w,a,b) \, f(y_n|u) du .}
#'
#' This distribution is the marginal distribution of the data under
#' the mixture density. For binary and Poisson data \eqn{y_n =
#' \sum_{i=1}^n y_i} is the sum over future events. For normal data,
#' it is the mean\eqn{\bar{y}_n = 1/n \sum_{i=1}^n y_i}.
#'
#' @return The function returns for a normal, beta or gamma mixture
#' the matching predictive distribution for \eqn{y_n}.
#'
#' @template conjugate_pairs
#'
#' @examples
#'
#' # Example 1: predictive distribution from uniform prior.
#' bm <- mixbeta(c(1,1,1))
#' bmPred <- preddist(bm, n=10)
#' # predictive proabilities and cumulative predictive probabilities
#' x <- 0:10
#' d <- dmix(bmPred, x)
#' names(d) <- x
#' barplot(d)
#' cd <- pmix(bmPred, x)
#' names(cd) <- x
#' barplot(cd)
#' # median
#' mdn <- qmix(bmPred,0.5)
#' mdn
#'
#' # Example 2: 2-comp Beta mixture
#'
#' bm <- mixbeta( inf=c(0.8,15,50),rob=c(0.2,1,1))
#' plot(bm)
#' bmPred <- preddist(bm,n=10)
#' plot(bmPred)
#' mdn <- qmix(bmPred,0.5)
#' mdn
#' d <- dmix(bmPred,x=0:10)
#' \donttest{
#' n.sim <- 100000
#' r <-  rmix(bmPred,n.sim)
#' d
#' table(r)/n.sim
#' }
#'
#' # Example 3: 3-comp Normal mixture
#'
#' m3 <- mixnorm( c(0.50,-0.2,0.1),c(0.25,0,0.2), c(0.25,0,0.5), sigma=10)
#' print(m3)
#' summary(m3)
#' plot(m3)
#' predm3 <- preddist(m3,n=2)
#' plot(predm3)
#' print(predm3)
#' summary(predm3)
#'
#' @export
preddist <- function(mix, ...) UseMethod("preddist")
#' @export
preddist.default <- function(mix, ...) stop("Unknown distribution")

#' @describeIn preddist Obtain the matching predictive distribution
#'     for a beta distribution, the BetaBinomial.
#' @export
preddist.betaMix <- function(mix, n=1, ...) {
    attr(mix, "n") <- n
    class(mix) <- c("betaBinomialMix", "mix")
    mix
}

#' @describeIn preddist Obtain the matching predictive distribution
#'     for a Normal with constant standard deviation. Note that the
#'     reference scale of the returned Normal mixture is meaningless
#'     as the individual components are updated appropriatley.
#' @param sigma The fixed reference scale of a normal mixture. If left
#'     unspecified, the default reference scale of the mixture is
#'     assumed.
#' @export
preddist.normMix <- function(mix, n=1, sigma, ...) {
    if(missing(sigma)) {
        sigma <- RBesT::sigma(mix)
        message("Using default mixture reference scale ", sigma)
    }
    assert_number(sigma, lower=0)
    sigma_ref <- sigma
    ## note: this is effectivley a hierarchical model as we give the
    ## distribution of the sum of n variables which have exactly the
    ## same mean (which is sampled from the normal)
    ## old: sum over y_i
    ##mix[2,] <- mix[2,] * n
    ##mix[3,] <- sqrt(mix[3,]^2 * n^2 + tau^2 * n )
    ## now: predictive for \bar{y}_n, the mean
    mix[3,] <- sqrt(mix[3,]^2  + sigma_ref^2/n )
    class(mix) <- c("normMix", "mix")
    mix
}

#' @describeIn preddist Obtain the matching predictive distribution
#'     for a Gamma. Only Poisson likelihoods are supported.
#' @export
preddist.gammaMix <- function(mix, n=1, ...) {
    assert_set_equal(likelihood(mix), "poisson")
    attr(mix, "n") <- n
    class(mix) <- c(switch(likelihood(mix),
                           poisson="gammaPoissonMix",
                           exp="gammaExpMix"), "mix")
    mix
}

#' @describeIn preddist Multivariate normal mixtures predictive
#'     densities are not (yet) supported.
#' @export
preddist.mvnormMix <- function(mix, ...) stop("Multivariate normal mixture predictive density not supported.")