Nothing
SdPrior <- function(sigma.guess,
sample.size = .01,
initial.value = sigma.guess,
fixed = FALSE,
upper.limit = Inf) {
## Generates an object of class SdPrior that can be used as an input
## to a Bayesian model for a standard deviation paramter.
##
## Args:
## sigma.guess: An a priori guess at the value of the standard
## deviation modeled by this distribution.
## sample.size: The number of observations worth of weight to give
## to sigma.guess.
## initial.value: The initial value to use (in an MCMC algorithm)
## for a deviate modeled by this distribution.
## fixed: Should the deviate modeled by this distribution be fixed
## at its initial value? (Used for debugging by some code. Not
## universal.)
## upper.limit: The largest acceptable value for the standard
## deviation modeled by this distribution. May be Inf.
stopifnot(is.numeric(sigma.guess), length(sigma.guess) == 1, sigma.guess > 0)
stopifnot(is.numeric(sample.size), length(sample.size) == 1, sample.size > 0)
stopifnot(is.numeric(initial.value), length(initial.value) == 1,
initial.value > 0)
stopifnot(is.logical(fixed), length(fixed) == 1)
stopifnot(is.numeric(upper.limit), length(upper.limit) == 1)
ans <- list(prior.guess = sigma.guess,
prior.df = sample.size,
initial.value = initial.value,
fixed = fixed,
upper.limit = upper.limit)
class(ans) <- c("SdPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
NormalPrior <- function(mu, sigma, initial.value = mu, fixed = FALSE) {
## Returns a list with the information needed to specify a Gaussian
## prior on a scalar parameter.
##
## Args:
## mu: Mean of the distribution.
## sigma: Standard deviation.
## initial.value: The initial value of the deviate to be modeled
## by this distribution in an MCMC algorithm.
## fixed: Should the deviate modeled by this distribution be fixed
## at its initial value? (Used for debugging by some code. Not
## universal.)
stopifnot(is.numeric(mu), length(mu) == 1)
stopifnot(is.numeric(sigma), length(sigma) == 1, sigma >= 0)
stopifnot(is.numeric(initial.value), length(initial.value) == 1)
stopifnot(is.logical(fixed), length(fixed) == 1)
ans <- list(mu = mu,
sigma = sigma,
initial.value = initial.value,
fixed = fixed)
class(ans) <- c("NormalPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
Ar1CoefficientPrior <- function(mu = 0,
sigma = 1,
force.stationary = TRUE,
force.positive = FALSE,
initial.value = mu) {
## Returns a list with the information needed to supply a prior
## distribution on an AR1 coefficient.
ans <- NormalPrior(mu, sigma, initial.value)
ans$force.stationary <- force.stationary
ans$force.positive <- force.positive
class(ans) <- c("Ar1CoefficientPrior", class(ans))
return(ans)
}
BetaPrior <- function(a = 1, b = 1, mean = NULL, sample.size = NULL,
initial.value = NULL) {
## Returns an object of class "BetaPrior", which is a list
## containing the parameters of a beta distribution. The prior can
## either be given in terms of 'a' and 'b', or it can be given in
## terms of mean and sample.size, where mean = a/a+b and sample.size
## = a+b.
if (!is.null(sample.size) && !is.null(mean)) {
stopifnot(is.numeric(mean) &&
is.numeric(sample.size) &&
length(mean) == 1 &&
length(sample.size) == 1)
a <- mean * sample.size
b <- (1 - mean) * sample.size
}
stopifnot(is.numeric(a),
length(a) == 1,
a > 0)
stopifnot(is.numeric(b),
length(b) == 1,
b > 0)
if (is.null(initial.value)) {
initial.value <- a / (a + b)
}
stopifnot(is.numeric(initial.value),
length(initial.value) == 1,
initial.value > 0,
initial.value < 1);
ans <- list(a = a, b = b, initial.value = initial.value)
class(ans) <- c("BetaPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
UniformPrior <- function(lo = 0, hi = 1, initial.value = NULL) {
## A uniform prior distribution on the interval [lo, hi].
## Args:
## lo: The lower limit of prior support.
## hi: The upper limit of prior support.
## Returns:
## An object of class UniformPrior.
stopifnot(is.numeric(lo),
length(lo) == 1)
stopifnot(is.numeric(hi),
length(lo) == 1,
lo <= hi)
if (is.null(initial.value)) {
initial.value <- .5 * (lo + hi)
}
stopifnot(is.numeric(initial.value),
length(initial.value) == 1,
initial.value >= lo,
initial.value <= hi)
ans <- list(lo = lo, hi = hi, initial.value = initial.value)
class(ans) <- c("UniformPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
GammaPrior <- function(a = NULL, b = NULL, prior.mean = NULL,
initial.value = NULL) {
## Gamma distribution with parameters (a, b), where the mean is a/b
## and variance is a/b^2.
## Args:
## a: The shape parameter 'a'.
## b: The scale parameter 'b'.
## prior.mean: This specifies a/b. If non-NULL at least one of
## 'a' or 'b' must be specified.
## initial.value: The initial.value of the variable to be modeled
## in the MCMC algorithm.
## Returns:
## An object of class GammaPrior.
if (is.null(prior.mean)) {
stopifnot(!is.null(a))
stopifnot(is.numeric(a))
} else {
stopifnot(is.numeric(prior.mean))
if (!is.null(a) && is.numeric(a)) {
b <- a / prior.mean
} else if (!is.null(b) && is.numeric(b)) {
a <- b * prior.mean
}
}
stopifnot(is.numeric(a))
stopifnot(is.numeric(b))
stopifnot(length(a) == 1)
stopifnot(length(b) == 1)
stopifnot(all(a > 0))
stopifnot(all(b > 0))
if (is.null(initial.value)) {
initial.value <- a / b
}
ans <- list(a = a, b = b, initial.value = initial.value)
class(ans) <- c("GammaPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
TruncatedGammaPrior <- function(a = NULL, b = NULL, prior.mean = NULL,
initial.value = NULL,
lower.truncation.point = 0,
upper.truncation.point = Inf) {
## A Gamma distribution with parameter matching 'GammaPrior' but
## with support truncated to between lower.truncation.point and
## upper.truncation.point.
ans <- GammaPrior(a, b, prior.mean, initial.value)
stopifnot(is.numeric(lower.truncation.point),
length(lower.truncation.point) == 1)
stopifnot(is.numeric(upper.truncation.point),
length(upper.truncation.point) == 1,
upper.truncation.point > lower.truncation.point)
ans$lower.truncation.point <- lower.truncation.point
ans$upper.truncation.point <- upper.truncation.point
if (!is.null(initial.value) && (
initial.value < lower.truncation.point ||
initial.value > upper.truncation.point)) {
stop("Initial value must be between the lower and upper truncation points.")
} else if (ans$initial.value < lower.truncation.point
|| ans$initial.value > upper.truncation.point) {
if (is.finite(upper.truncation.point)) {
warning("Initial value for TruncatedGammaPrior is outside the ",
"supported range. Changing it to the midpoint of the ",
"support region.")
ans$initial.value <- .5 * (lower.truncation.point +
upper.truncation.point)
} else {
warning("Initial value for TruncatedGammaPrior is outside the ",
"supported range. Placing it above the lower bound.")
ans$initial.value <- lower.truncation.point + 1
}
}
class(ans) <- c("TruncatedGammaPrior", "DiffDoubleModel",
"DoubleModel", "Prior")
return(ans)
}
LognormalPrior <- function(mu = 0.0, sigma = 1.0, initial.value = NULL) {
## A lognormal distribution, where log(y) ~ N(mu, sigma). The mean
## of this distribution is exp(mu + 0.5 * sigma^2), so don't only
## focus on the mean parameter here.
##
## Args:
## mu: mean of the corresponding normal distribution.
## sigma: standard deviation of the corresponding normal
## distribution. WARNING: If something looks strange in your
## program, look out for SD != Variance errors.
stopifnot(is.numeric(mu),
length(mu) == 1,
is.finite(mu))
stopifnot(is.numeric(sigma),
length(sigma) == 1,
sigma > 0,
is.finite(sigma))
if (is.null(initial.value)) {
initial.value <- exp(mu + .5 * sigma^2)
}
stopifnot(is.numeric(initial.value),
length(initial.value) == 1,
initial.value > 0,
is.finite(initial.value))
ans <- list(mu = mu, sigma = sigma, initial.value = initial.value)
class(ans) <- c("LognormalPrior", "DiffDoubleModel", "DoubleModel", "Prior")
return(ans)
}
MarkovPrior <- function(prior.transition.counts = NULL,
prior.initial.state.counts = NULL,
state.space.size = NULL,
uniform.prior.value = 1) {
## Prior distribution for a MarkovModel.
## Args:
## prior.transition.counts: A matrix of non-negative numbers
## interpretable as prior counts. Transitions are from rows to
## columns.
## prior.initial.state.counts: A vector of non-negative numbers
## interpretable as prior counts.
## state.space.size: If both prior.transition.counts and
## prior.initial.state.counts are missing then they will be filled
## with an object of dimension state.space.size where all
## entries are set to uniform.prior.value.
## uniform.prior.value: The default value to use for entries of
## prior.transition.counts and prior.initial.state.counts, when
## they are not supplied by the user.
if (is.null(state.space.size)) {
if (is.null(prior.transition.counts) && is.null(prior.transition.counts)) {
stop("Either 'state.space.size' or one of 'prior.transition.counts' or",
"'prior.initial.state.counts' must be supplied to MarkovPrior")
}
if (!is.null(prior.transition.counts)) {
state.space.size <- nrow(prior.transition.counts)
} else {
state.space.size <- length(prior.initial.state.counts)
}
}
stopifnot(state.space.size > 0)
stopifnot(uniform.prior.value > 0)
if (is.null(prior.transition.counts)) {
prior.transition.counts <- matrix(uniform.prior.value, nrow =
state.space.size, ncol = state.space.size)
} else {
stopifnot(is.matrix(prior.transition.counts))
stopifnot(nrow(prior.transition.counts) == ncol(prior.transition.counts))
stopifnot(all(prior.transition.counts >= 0))
.CheckForPositiveValue <- function(x) { return(any(x > 0)) }
## Check that at least one positive value is present in each row.
stopifnot(all(apply(prior.transition.counts, 1, .CheckForPositiveValue)))
## If state.space.size and prior.transition.counts are both
## present then issue a warning if they don't match, and use
## prior.transition.counts.
if (state.space.size != nrow(prior.transition.counts)) {
warning("state.space.size is ", state.space.size,
", but nrow(prior.transition.counts) is",
nrow(prior.transition.counts),
". Changing state.space.size to nrow(prior.transition.counts).")
state.space.size <- nrow(prior.transition.counts)
}
}
if (is.null(prior.initial.state.counts)) {
prior.initial.state.counts <- rep(uniform.prior.value, state.space.size)
} else {
stopifnot(is.numeric(prior.initial.state.counts))
## Allow a prior that places all its mass on a single point.
stopifnot(all(prior.initial.state.counts >= 0))
stopifnot(any(prior.initial.state.counts > 0))
}
prior <- list(prior.transition.counts = prior.transition.counts,
prior.initial.state.counts = prior.initial.state.counts)
class(prior) <- c("MarkovPrior", "Prior")
return(prior)
}
DirichletPrior <- function(prior.counts, initial.value = NULL) {
## Encodes a Dirichlet priors, the conjugate prior for the
## multinomial distribution.
## Args:
## prior.counts: A vector of positive numbers with dimension
## matching the probability distribution it is modeling.
## Returns:
## An object of class DirichletPrior, which is a list containing
## the prior.counts argument, after some type and sanity checking.
stopifnot(is.numeric(prior.counts))
stopifnot(length(prior.counts) > 0)
stopifnot(all(prior.counts > 0))
if (is.null(initial.value)) {
initial.value <- prior.counts / sum(prior.counts)
}
ans <- list(prior.counts = prior.counts)
class(ans) <- c("DirichletPrior", "Prior")
return(ans)
}
NormalInverseGammaPrior <- function(mu.guess,
mu.guess.weight = .01,
sigma.guess,
sigma.guess.weight = 1,
...) {
## A conjugate prior for the mean and variance of a Gaussian
## distribution.
## Args:
## mu.guess: Prior guess at the normal mean parameter.
## mu.guess.weight: Number of observations worth of weight
## assigned to mu.guess.
## sigma.guess: A prior guess at the value of the normal standard
## deviation parameter.
## sigma.guess.weight: Number of observations worth of weight
## assigned to sigma.guess.
## ...: extra parameters passed to SdPrior
## Returns:
## An object of class NormalInverseGammaPrior, which contains the
## 'mu' arguments and an element of type SdPrior
stopifnot(is.numeric(mu.guess))
stopifnot(mu.guess.weight > 0)
ans <- list(mu.guess = mu.guess,
mu.guess.weight = mu.guess.weight,
sigma.prior = SdPrior(
sigma.guess = sigma.guess,
sample.size = sigma.guess.weight,
...))
class(ans) <- c("NormalInverseGammaPrior", "Prior")
return(ans)
}
MvnPrior <- function(mean, variance) {
## Encodes the mean and variance of a multivariate normal
## distribution for use as a prior distribution.
## Args:
## mean: A numeric vector
## variance: A symmetric positive definite matrix.
## Returns:
## An object of class MvnPrior containing the mean and variance
## arguments, after some sanity checking.
if (length(variance) == 1) {
variance <- diag(as.numeric(variance), length(mean), length(mean))
}
stopifnot(is.matrix(variance))
stopifnot(nrow(variance) == ncol(variance))
stopifnot(nrow(variance) == length(mean))
ans <- list(mean = mean,
variance = variance)
class(ans) <- c("MvnPrior", "Prior")
return(ans)
}
MvnGivenSigmaMatrixPrior <- function(mean, sample.size) {
stopifnot(is.numeric(mean))
stopifnot(is.numeric(sample.size),
length(sample.size) == 1,
sample.size > 0)
ans <- list("mean" = mean, "sample.size" = sample.size)
class(ans) <- c("MvnGivenSigmaMatrixPrior", "Prior")
return(ans)
}
InverseWishartPrior <- function(variance.guess,
variance.guess.weight) {
## Conjugate prior distribution for the variance matrix in a multivariate
## normal wtih known mean.
stopifnot(is.matrix(variance.guess))
stopifnot(is.numeric(variance.guess.weight),
length(variance.guess.weight) == 1,
variance.guess.weight > ncol(variance.guess))
ans <- list(variance.guess = variance.guess,
variance.guess.weight = variance.guess.weight)
class(ans) <- c("InverseWishartPrior", "Prior")
return(ans)
}
NormalInverseWishartPrior <- function(
mean.guess,
mean.guess.weight = .01,
variance.guess,
variance.guess.weight = nrow(variance.guess) + 1) {
## The conjugate prior distribution for the multivariate normal. A
## multivariate generalization of the NormalInverseGammaPrior
## Args:
## mean.guess: A numeric vector that is a guess at the value of
## the multivariate normal mean parameter.
## mean.guess.weight: The number of observations worth of weight
## assigned to mean.guess.
## variance.guess: A symmetric positive definite matrix that is a
## guess at the value of the multivariate normal variance
## parameter.
## variance.guess.weight: The number of observations worth of
## weight assigned to variance.guess.
stopifnot(is.vector(mean.guess))
stopifnot(is.matrix(variance.guess))
stopifnot(all(dim(variance.guess) == length(mean.guess)))
stopifnot(is.numeric(mean.guess.weight),
length(mean.guess.weight) == 1,
mean.guess.weight > 0)
stopifnot(is.numeric(variance.guess.weight),
length(variance.guess.weight) == 1,
variance.guess.weight > 0)
ans <- list(mean.guess = mean.guess,
mean.guess.weight = mean.guess.weight,
variance.guess = variance.guess,
variance.guess.weight = variance.guess.weight)
class(ans) <- c("NormalInverseWishartPrior", "Prior")
return(ans)
}
MvnDiagonalPrior <- function(mean.vector, sd.vector) {
## A multivariate normal distribution with a diagonal variance
## matrix (i.e. the product of several independent normals).
##
## Args:
## mean.vector: The mean of the multivariate normal.
## sd.vector: The standard deviations of the elements of the
## multivariate normal.
##
## Returns:
## An object of class MvnDiagonalPrior, which is a list containing
## mean.vector and sd.vector.
stopifnot(is.numeric(mean.vector))
stopifnot(is.numeric(sd.vector))
stopifnot(all(sd.vector > 0))
stopifnot(length(mean.vector) == length(sd.vector))
ans <- list(mean = mean.vector, sd = sd.vector)
class(ans) <- c("MvnDiagonalPrior", "Prior")
return(ans)
}
MvnIndependentSigmaPrior <- function(mvn.prior, sd.prior.list) {
## A prior for the parameters of the multivariate normal
## distribution that assumes Sigma to be a diagonal matrix with
## elements modeled by independent inverse Gamma priors.
##
## Args:
## mvn.prior: An object of class MvnPrior that is the prior
## distribution of the multivariate normal mean parameter.
## sd.prior.list: A list of SdPrior object modeling the diagonal
## elements of the multivariate normal variance matrix. The
## off-diagonal elements are assumed to be zero.
##
## Returns:
## An object of clas MvnIndependentSigmaPrior, which is a list
## containing the function arguments.
stopifnot(inherits(mvn.prior, "MvnPrior"))
stopifnot(is.list(sd.prior.list))
stopifnot(length(sd.prior.list) == length(mvn.prior$mean))
stopifnot(all(sapply(sd.prior.list, inherits, "SdPrior")))
ans <- list(mu.prior = mvn.prior,
sigma.prior = sd.prior.list)
class(ans) <- c("MvnIndependentSigmaPrior", "Prior")
return(ans)
}
ScaledMatrixNormalPrior <- function(mean, nu) {
# A prior for multivariate regression coefficients Beta, such that
# Beta | Sigma ~ MatrixNormal(mu, X'X / nu, Sigma)
#
# Args:
# mean: A matrix of the same size as Beta (ydim * xdim).
# nu: A scalar controlling the amount of prior correlation attributable to
# rows of beta. In the context of multivariate regression, rows of beta
# correspond ot different Y variables. A useful value for nu is kappa /
# sample.size, where kappa is the number of prior observations worth of
# weight you wish to give to 'mean'.
#
# Details:
# The Sigma parameter is obtained elsewhere.
stopifnot(is.matrix(mean))
stopifnot(is.numeric(nu), length(nu) == 1, nu > 0)
ans <- list("mean" = mean, "nu" = nu)
class(ans) <- c("ScaledMatrixNormalPrior", "Prior")
return(ans)
}
MultivariateRegressionConjugatePrior <-function(coefficient.prior,
variance.prior) {
## A conjguate prior for a multivariate regression model.
##
## Args:
## coefficient.prior: An object of class ScaledMatrixNormalPrior
## representing the prior on the regression coefficients.
## variance.prior: An object of class InverseWishartPrior representing the
## prior on the variance matrix.
##
## Returns:
## A list containing the two arguments, after checking that they are the
## right type. The return value is marked with a class attribute
## documenting its role.
stopifnot(inherits(coefficient.prior, "ScaledMatrixNormalPrior"),
inherits(variance.prior, "InverseWishartPrior"))
ans <- list(coefficient = coefficient.prior, variance = variance.prior)
class(ans) <- c("MultivariateRegressionConjugatePrior", "Prior")
return(ans)
}
RegressionCoefficientConjugatePrior <- function(
mean,
sample.size,
additional.prior.precision = numeric(0),
diagonal.weight = 0) {
## A conditional prior for the coefficients (beta) in a linear regression
## model. The prior is conditional on the residual variance sigma^2, the
## sample size n, and the design matrix X. The prior is
##
## beta | sigsq, X ~ N(b, sigsq * (Lambda^{-1} + V))
##
## where V^{-1} = ((1 - w) * XTX + w * Diag(XTX)) * kappa / n.
##
## Args:
## mean: The mean of the prior distribution, denoted 'b' above.
## sample.size: The value denoted 'kappa' above. This can be
## interpreted as a number of observations worth of weight to be
## assigned to 'b' in the posterior distribution.
## additional.prior.precision: A vector of non-negative numbers
## representing the diagonal matrix Lambda^{-1} above. Positive
## values for additional.prior.precision will ensure the
## distribution is proper even if the regression model has no
## data. If all columns of the design matrix have positive
## variance then additional.prior.precision can safely be set to
## zero. A zero-length numeric vector is a slightly more
## efficient equivalent to a vector of all zeros.
## diagonal.weight: The weight given to the diagonal when XTX is
## averaged with its diagonal. The purpose of diagonal.weight
## is to keep the prior distribution proper even if X is less
## than full rank. If the design matrix is full rank then
## diagonal.weight can be set to zero.
##
## Details:
## The prior distribution also depends on the cross product matrix
## XTX and the sample size n, which are not arguments to this
## function. It is expected that the underlying C++ code will get
## those quantities elsewhere (presumably from the regression
## modeled by this prior).
stopifnot(is.numeric(mean))
stopifnot(is.numeric(sample.size),
length(sample.size) == 1,
sample.size > 0)
stopifnot(is.numeric(additional.prior.precision),
length(additional.prior.precision) == length(mean)
|| length(additional.prior.precision) == 0)
if (length(additional.prior.precision) > 0) {
stopifnot(all(additional.prior.precision >= 0))
}
stopifnot(is.numeric(diagonal.weight),
length(diagonal.weight) == 1,
diagonal.weight >= 0,
diagonal.weight <= 1)
ans <- list(mean = mean,
sample.size = sample.size,
additional.prior.precision = additional.prior.precision,
diagonal.weight = diagonal.weight)
class(ans) <- c("RegressionCoefficientConjugatePrior")
return(ans)
}
PointMassPrior <- function(location) {
## A prior putting a point mass (probability = 1) on a single scalar
## location.
## Args:
## location: The location of the point mass on the real line.
ans <- list(location = location)
class(ans) <- c("PointMassPrior", "DiscretePrior", "Prior")
return(ans)
}
PoissonPrior <- function(mean, lower.limit = 0, upper.limit = Inf) {
## Prior over the positive integers based on a Poisson distribution.
## Can optionally be truncated to {lower.limit, ..., upper.limit},
## including the endpoints.
ans <- list(mean = mean,
lower.limit = lower.limit,
upper.limit = upper.limit)
class(ans) <- c("PoissonPrior", "DiscretePrior", "Prior")
return(ans)
}
DiscreteUniformPrior <- function(lower.limit, upper.limit) {
## A discrete uniform distribution over the integers in
## {lower.limit, ..., upper.limit}. The end points are included.
## Args:
## lower.limit: The lower limit of the support.
## upper.limit: The upper limit of the support.
ans <- list(lower.limit = lower.limit,
upper.limit = upper.limit)
class(ans) <- c("DiscreteUniformPrior", "DiscretePrior", "Prior")
return(ans)
}
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.