Boom

Description

The Boom package provides access to the C++ BOOM library for Bayesian computation.

Details

Installation note for Linux users

If you are installing Boom using install.packages on a Linux machine (and thus compiling yourself) you will almost certainly want to set the Ncpus argument to a large number. Windows and Mac users can ignore this advice.

The main purpose of the Boom package is not to be used directly, but to provide the BOOM C++ library for other packages to link against. The Boom package provides additional utility code for C++ authors to use when writing R packages with C++ internals. These are described in .../inst/include/r_interface/boom_r_tools.hpp among the package's include files.

Boom provides a collection of R functions and objects to help users format data in the manner expected by the underlying C++ code. Standard distributions that are commonly used as Bayesian priors can be specified using BetaPrior, GammaPrior, etc.

Boom provides a set of utilities helpful when writing unit tests for Bayesian models. See CheckMcmcMatrix and CheckMcmcVector for MCMC output, and functions like check.probability.distribution for checking function inputs

Boom provides a collection of useful plots (using base R graphics) that have proven useful for summarizing MCMC output. See PlotDynamicDistribution, PlotManyTs, BoxplotTrue, and other code in the index with Plot in the title.

See Also

Please see the following pacakges

• bsts

• CausalImpact

Generate a data frame of all factor data

Description

This function is mainly intended for example code and unit testing. It generates a data.frame containing all factor data.

Usage

  GenerateFactorData(factor.levels.list, sample.size)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  foo <- GenerateFactorData(list(a = c("foo", "bar", "baz"),
                                 b = c("larry", "moe", "curly", "shemp")),
                            50)

  head(foo)
#     a     b
# 1 bar curly
# 2 foo curly
# 3 bar   moe
# 4 bar   moe
# 5 baz curly
# 6 bar curly

Conditional Multivaraite Normal Prior Given Variance

Description

A multivaraite normal prior distribution, typically used as the prior distribution for the mean of multivaraite normal data. The variance of this distribution is proportional to another parameter "Sigma" that exists elsewhere. Usually "Sigma" is the variance of the data.

This distribution is the "normal" part of a "normal-inverse Wishart" distribution.

Usage

MvnGivenSigmaMatrixPrior(mean, sample.size)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Time Series Boxplots

Description

Creates a series of boxplots showing the evolution of a distribution over time.

Usage

  TimeSeriesBoxplot(x, time, ylim = NULL, add = FALSE, ...)

Arguments

Value

Called for its side effect, which is to produce a plot on the current graphics device.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  x <- t(matrix(rnorm(1000 * 100, 1:100, 1:100), nrow=100))
  ## x has 1000 rows, and 100 columns.  Column i is N(i, i^2) noise.
  time <- as.Date("2010-01-01", format = "%Y-%m-%d") + (0:99 - 50)*7
  TimeSeriesBoxplot(x, time)

Convert to Character String

Description

Convert an object to a character string, suitable for including in error messages.

Usage

  ToString(object, ...)

Arguments

Value

A string (a character vector of length 1) containing the printed value of the object.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  m <- matrix(1:6, ncol = 2)
  printed.matrix <- ToString(m)

  y <- c(1, 2, 3, 3, 3, 3, 3, 3)
  tab <- table(y)
  printed.table <- ToString(tab)

Function to add horizontal line segments to an existing plot

Description

Adds horizontal line segments to an existing plot. The segments are centered at x with height y. The x values are assumed to be equally spaced, so that diff(x) is a constant 'dx'. The line segments go from x +/- half.width.factor *dx, so if half.width.factor=.5 there will be no gaps between segments. The default is to leave a small gap.

This function was originally used to add reference lines to side-by-side boxplots.

Usage

AddSegments(x, y, half.width.factor = 0.45, ...)

Arguments

Value

Called for its side effect.

Author(s)

Steven L. Scott

See Also

boxplot.true

Examples

x <- rnorm(100)
y <- rnorm(100, 1)
boxplot(list(x=x,y=y))
AddSegments(1:2, c(0, 1))  ## add segments to the boxplot

Normal prior for an AR1 coefficient

Description

A (possibly truncated) Gaussian prior on the autoregression coefficient in an AR1 model.

Usage

Ar1CoefficientPrior(mu = 0, sigma = 1, force.stationary = TRUE,
    force.positive = FALSE, initial.value = mu)

Arguments

Details

The Ar1CoefficientPrior() syntax is preferred, as it more closely matches R's syntax for other constructors.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Beta prior for a binomial proportion

Description

Specifies beta prior distribution for a binomial probability parameter.

Usage

  BetaPrior(a = 1, b = 1, mean = NULL, sample.size = NULL,
            initial.value = NULL)

Arguments

Details

The distribution should be specified either with a and b, or with mean and sample.size.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Plot the distribution of a matrix

Description

Plot the marginal distribution of each element in the Monte Carlo distribution of a matrix (e.g. a variance matrix or transition probability matrix). Rows and columns in the boxplots correspond to rows and columns in the matrix being plotted.

Usage

  BoxplotMcmcMatrix(X, ylim = range(X), col.names,
                    row.names, truth, colors = NULL,
                    las = 0, ...)

Arguments

Value

Called for its side effect, which is to draw a set of side-by-side boxplots on the current graphics device.

Author(s)

Steven L. Scott

See Also

boxplot.true, boxplot

Examples

  X <- array(rnorm(1000 * 3 * 4), dim=c(1000, 3, 4))
  dimnames(X)[[2]] <- paste("row", 1:3)
  dimnames(X)[[3]] <- paste("col", 1:4)
  BoxplotMcmcMatrix(X)

  truth <- 0
  BoxplotMcmcMatrix(X, truth=truth)

  truth <- matrix(rnorm(12), ncol=4)
  BoxplotMcmcMatrix(X, truth=truth)

Compare Boxplots to True Values

Description

Plots side-by-side boxplots of the columns of the matrix x. Each boxplot can have its own reference line (truth) and standard error lines se.truth, if desired. This function was originally written to display MCMC output, where the reference lines were true values used to test an MCMC simulation.

Usage

BoxplotTrue(x, truth = NULL, vnames = NULL, center = FALSE,
            se.truth = NULL, color = "white", truth.color = "black",
            ylim = NULL, ...)

Arguments

Value

called for its side effect

Author(s)

Steven L. Scott

See Also

boxplot.matrix, boxplot,

Examples

x <- t(matrix(rnorm(5000, 1:5, 1:5), nrow=5))
BoxplotTrue(x, truth=1:5, se.truth=1:5, col=rainbow(5), vnames =
  c("EJ", "TK", "JT", "OtherEJ", "TJ") )

Check MCMC Output

Description

Verify that MCMC output covers expected values.

Usage

CheckMcmcMatrix(draws, truth, confidence = .95,
                control.multiple.comparisons = TRUE,
                burn = 0)

CheckMcmcVector(draws, truth, confidence = .95, burn = 0)

McmcMatrixReport(draws, truth, confidence = .95, burn = 0)

Arguments

Details

CheckMcmcVector checks a vector of draws corresponding to a scalar random variable. CheckMcmcMatrix checks a matrix of draws corresponding to a vector of random variables. In either case the check is made by constructing a central confidence interval (obtained by removing half of 1 - confidence from the upper and lower tails of the distribution).

If a single variable is being checked with CheckMcmcVector then the check passes if and only if the interval covers the true value.

If multiple values are being checked with CheckMcmcMatrix then the user has control over how strict to make the check. If control.multiple.comparisons is FALSE then the check passes if and only if all intervals cover true values. Otherwise a fraction of intervals must cover. The fraction is the lower bound of the binomial confidence interval for the coverage rate under the hypothesis that the true coverage rate is confidence.

Value

CheckMcmcVector and CheckMcmcMatrix return TRUE if the check passes, and FALSE if it does not.

McmcMatrixReport returns a string that can be put in the info field of an expect_true expression, to give useful information about a failed test case. The return value is a textual representation of a three column matrix. Each row matches a variable in draws, and gives the lower and upper bounds for the credible interval used to check the values. The final column lists the true values that are supposed to be inside the credible intervals. The value is returned as a character string that is expected to be fed to cat() or print() so that it will render correctly in R CMD CHECK output.

Author(s)

Steven L. Scott

Examples

ndraws <- 100
draws <- rnorm(ndraws, 0, 1)
CheckMcmcVector(draws, 0)  ## Returns TRUE
CheckMcmcVector(draws, 17)  ## Returns FALSE

draws <- matrix(nrow = ndraws, ncol = 5)
for (i in 1:5) {
  draws[, i] <- rnorm(ndraws, i, 1)
}  
CheckMcmcMatrix(draws, truth = 1:5)  ## Returns TRUE
CheckMcmcMatrix(draws, truth = 5:1)  ## Returns FALSE

Checking data formats

Description

Checks that data matches a concept

Usage

check.scalar.probability(x)
check.positive.scalar(x)
check.nonnegative.scalar(x)
check.probability.distribution(x)
check.scalar.integer(x)
check.scalar.boolean(x)

Arguments

Details

If the object does not match the concept being checked, stop is called. Otherwise TRUE is returned.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Draw Circles

Description

Draw circles on the current graphics device.

Usage

circles(center, radius, ...) 

Arguments

Details

Draws circles on the current graphics device. This is a low-level plotting function similar to points, lines, segments, etc.

Value

Returns invisible NULL.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  plot(1:10, type = "n")
  circles(cbind(c(2, 3, 4), c(4, 5,6 )), radius = c(.3, .4, .5))

Compare several density estimates.

Description

Produces multiple density plots on a single axis, to compare the columns of a matrix or the elements of a list.

Usage

CompareDensities(x,
                 legend.text = NULL,
                 legend.location = "topright",
                 legend.title = NULL,
                 xlim = NULL,
                 ylim = NULL,
                 xlab = "parameter",
                 ylab = "density",
                 main = "",
                 lty = NULL,
                 col = "black",
                 axes = TRUE,
                 na.rm = TRUE,
                 ...)

Arguments

Value

Called for its side effect, which is to produce multiple density plots on the current graphics device.

Author(s)

Steven L. Scott

See Also

density

Examples

x <- t(matrix(rnorm(5000, 1:5, 1:5), nrow=5))
CompareDensities(x, legend.text=c("EJ", "TK", "JT", "OtherEJ", "TJ"),
                 col=rainbow(5), lwd=2)

Compare Dynamic Distributions

Description

Produce a plot showing several stacked dynamic distributions over the same horizontal axis.

Usage

CompareDynamicDistributions(
    list.of.curves,
    timestamps,
    style = c("dynamic", "boxplot"),
    xlab = "Time",
    ylab = "",
    frame.labels = rep("", length(list.of.curves)),
    main = "",
    actuals = NULL,
    col.actuals = "blue",
    pch.actuals = 1,
    cex.actuals = 1,
    vertical.cuts = NULL,
    ...)

Arguments

Author(s)

Steven L. Scott

Compare several density estimates.

Description

Produce a plot that compares the kernel density estimates for each element in a series of Monte Carlo draws of a vector or matrix.

Usage

CompareManyDensities(list.of.arrays,
                     style = c("density", "box"),
                     main = "",
                     color = NULL,
                     gap = 0,
                     burn = 0,
                     suppress.labels = FALSE,
                     x.same.scale = TRUE,
                     y.same.scale = FALSE,
                     xlim = NULL,
                     ylim = NULL,
                     legend.location = c("top", "right"),
                     legend.cex = 1,
                     reflines = NULL,
                     ...)

Arguments

Author(s)

Steven L. Scott

See Also

density, CompareManyTs

Examples

x <- array(rnorm(9000), dim = c(1000, 3, 3))
dimnames(x) <- list(NULL, c("Larry", "Moe", "Curly"), c("Larry", "Eric", "Sergey"))
y <- array(rnorm(9000), dim = c(1000, 3, 3))
z <- array(rnorm(9000), dim = c(1000, 3, 3))
data <- list(x = x, y = y, z = z)
CompareManyDensities(data, color = c("red", "blue", "green"))
CompareManyDensities(data, style = "box")

x <- matrix(rnorm(5000), nrow = 1000)
colnames(x) <- c("Larry", "Moe", "Curly", "Shemp", "???")
y <- matrix(rnorm(5000), nrow = 1000)
z <- matrix(rnorm(5000), nrow = 1000)
data <- list(x = x, y = y, z = z)
CompareManyDensities(data, color = c("red", "blue", "green"))
CompareManyDensities(data, style = "box")

Compares several density estimates.

Description

Produce a plot that compares the kernel density estimates for each element in a series of Monte Carlo draws of a vector or matrix.

Usage

CompareManyTs(list.of.ts, burn = 0, type = "l", gap = 0,
              boxes = TRUE, thin = 1, labels = NULL,
              same.scale = TRUE, ylim = NULL, refline = NULL,
              color = NULL, ...)

Arguments

Author(s)

Steven L. Scott

See Also

PlotManyTs, CompareManyDensities

Examples

x <- array(rnorm(9000), dim = c(1000, 3, 3))
dimnames(x) <- list(NULL, c("Larry", "Moe", "Curly"), c("Larry", "Eric", "Sergey"))
y <- array(rnorm(9000), dim = c(1000, 3, 3))
z <- array(rnorm(9000), dim = c(1000, 3, 3))
data <- list(x = x, y = y, z = z)
CompareManyTs(data, color = c("red", "blue", "green"))

x <- matrix(rnorm(5000), nrow = 1000)
colnames(x) <- c("Larry", "Moe", "Curly", "Shemp", "???")
y <- matrix(rnorm(5000), nrow = 1000)
z <- matrix(rnorm(5000), nrow = 1000)
data <- list(x = x, y = y, z = z)
CompareManyTs(data, color = c("red", "blue", "green"))

Boxplots to compare distributions of vectors

Description

Uses boxplots to compare distributions of vectors.

Usage

CompareVectorBoxplots(draws, main = NULL, colors = NULL, burn = 0, ...)

Arguments

Details

Creates side-by-side boxplots with the dimensions of each vector gropued together.

Examples

x <- matrix(rnorm(300, mean = 1:3, sd = .4), ncol = 3, byrow = TRUE)
y <- matrix(rnorm(600, mean = 3:1, sd = .2), ncol = 3, byrow = TRUE)
CompareVectorBoxplots(list(x = x, y = y), colors = c("red", "blue"))

DiffDoubleModel

Description

A 'DiffDoubleModel' is tag given to a probability distribution that measures a real-valued scalar outcome, and whose log density is twice differentiable with respect to the random variable. The tag is a signal to underlying C++ code that the object being passed is one of a subset of understood distributions. Presently that subset includes the following distributions.

• SdPrior

• NormalPrior

• BetaPrior

• UniformPrior

• GammaPrior

• TruncatedGammaPrior

• LognormalPrior

Clearly this list is non-exhaustive, and other distributions may be added in the future.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

The Dirichlet Distribution

Description

Density and random generation for the Dirichlet distribution.

Usage

ddirichlet(probabilities, nu, logscale = FALSE)
rdirichlet(n, nu)

Arguments

Details

The Dirichlet distribution is a generalization of the beta distribution. Whereas beta distribution is a model for probabilities, the Dirichlet distribution is a model for discrete distributions with several possible outcome values.

Let \pi denote a discrete probability distribution (a vector of positive numbers summing to 1), and let \nu be a vector of positive numbers (the parameters of the Dirichlet distribution), which can be thought of as prior counts. Then the density of the Dirichlet distribution can be written

f(\pi) = \frac{\Gamma(\sum_i\nu_i)}{\prod_i \Gamma(\nu_i)} \prod_i \pi_i^{\nu_i-1}.

Value

ddirichlet returns a vector of density values, with one entry per row in probabilities. rdirichlet returns a matrix (if n > 1) or a vector (if n==1) containing the draws from the Dirichlet distribution with the specified parameters.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Dirichlet prior for a multinomial distribution

Description

Specifies Dirichlet prior for a discrete probability distribution.

Usage

 DirichletPrior(prior.counts, initial.value = NULL)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Discrete prior distributions

Description

Prior distributions over a discrete quantities.

Usage

PointMassPrior(location)
PoissonPrior(mean, lower.limit = 0, upper.limit = Inf)
DiscreteUniformPrior(lower.limit, upper.limit)

Arguments

Value

Each function returns a prior object whose class is the same as the function name. All of these inherit from "DiscreteUniformPrior" and from "Prior".

The PoissonPrior assumes a potentially truncated Poisson distribution with the given mean.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

## Specify an exact number of trees in a Bart model (see the BoomBart
## package).

ntrees <- PointMassPrior(200)

## Uniform prior between 50 and 100 trees, including the endpoints.
ntrees <- DiscreteUniformPrior(50, 100)

## Truncated Poisson prior, with a mean of 20, a lower endpoint of 1,
## and an upper endpoint of 50.
ntrees <- PoissonPrior(20, 1, 50)

Multivariate Normal Density

Description

Evaluate the multivariate normal density.

Usage

dmvn(y, mu, sigma, siginv = NULL, ldsi = NULL, logscale = FALSE)

Arguments

Value

A vector containing the density of each row of y.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Prior distributions for a real valued scalar

Description

A DoubleModel is a class of prior distributions for real valued scalar parameters. A DoubleModel is sometimes used by a probability model that does not have a conjugate prior.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

See Also

BetaPrior, GammaPrior, LognormalPrior, NormalPrior, SdPrior, TruncatedGammaPrior, and UniformPrior.

Add an external legend to an array of plots.

Description

ExternalLegendLayout sets up a plotting region to plot a regular grid of points, with an optional legend on the top or right of the grid.

AddExternalLegend adds a legend to a grid of plots that was set up using ExternalLegendLayout.

Usage

ExternalLegendLayout(nrow,
                     ncol,
                     legend.labels,
                     legend.location = c("top", "right"),
                     outer.margin.lines = rep(4, 4),
                     gap.between.plots = rep(0, 4),
                     legend.cex = 1,
                     x.axis = TRUE,
                     y.axis = TRUE)

AddExternalLegend(legend.labels,
                  legend.location = c("top", "right"),
                  legend.cex =1,
                  bty = "n",
                  ...)

Arguments

Value

ExternalLegendLayout returns the original graphical parameters, intended for use with on.exit.

AddExternalLegend returns invisible NULL.

Author(s)

Steven L. Scott

See Also

legend layout

Examples

example.plot <- function() {
  x <- rnorm(100)
  y <- rnorm(100)
  scale <- range(x, y)
  opar <- ExternalLegendLayout(nrow = 2,
                               ncol = 2,
                               legend.labels = c("foo", "bar"))
  on.exit({par(opar); layout(1)})
  hist(x, xlim = scale, axes = FALSE, main = "")
  mtext("X", side = 3, line = 1)
  box()
  plot(x, y, xlim = scale, ylim = scale, axes = FALSE)
  box()
  axis(3)
  axis(4)
  plot(y, x, xlim = scale, ylim = scale, axes = FALSE, pch = 2, col = 2)
  box()
  axis(1)
  axis(2)
  hist(y, xlim = scale, axes = FALSE, main = "")
  mtext("Y", side = 1, line = 1)
  box()
  AddExternalLegend(legend.labels = c("foo", "bar"),
                    pch = 1:2,
                    col = 1:2,
                    legend.cex = 1.5)
}

## Now call example.plot().
example.plot()

## Because of the call to on.exit(), in example.plot,
## the original plot layout is restored.
hist(1:10)

Gamma prior distribution

Description

Specifies gamma prior distribution.

Usage

  GammaPrior(a = NULL, b = NULL, prior.mean = NULL, initial.value = NULL)
  TruncatedGammaPrior(a = NULL, b = NULL, prior.mean = NULL,
                      initial.value = NULL,
                      lower.truncation.point = 0,
                      upper.truncation.point = Inf)

Arguments

Details

The mean of the Gamma(a, b) distribution is a/b and the variance is a/b^2. If prior.mean is not NULL, then one of either a or b must be non-NULL as well.

GammaPrior is the conjugate prior for a Poisson mean or an exponential rate. For a Poisson mean a corresponds to a prior sum of observations and b to a prior number of observations. For an exponential rate the roles are reversed a represents a number of observations and b the sum of the observed durations. The gamma distribution is a generally useful for parameters that must be positive.

The gamma distribution is the conjugate prior for the reciprocal of a Guassian variance, but SdPrior should usually be used in that case.

A TruncatedGammaPrior is a GammaPrior with support truncated to the interval (lower.truncation.point, upper.truncation.point).

If an object specifically needs a GammaPrior you typically cannot pass a TruncatedGammaPrior.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

A Bunch of Histograms

Description

Plot a bunch of histograms describing the marginal distributions the columns in a data frame.

Usage

histabunch(x, gap = 1, same.scale = FALSE, boxes = FALSE,
           min.continuous = 12, max.factor = 40,
           vertical.axes = FALSE, ...)

Arguments

Value

Called for its side effect, which is to produce multiple histograms on the current graphics device.

Author(s)

Steven L. Scott

See Also

hist barplot

Examples

  data(airquality)
  histabunch(airquality)

Inverse Wishart Distribution

Description

Density for the inverse Wishart distribution.

Usage

dInverseWishart(Sigma, sum.of.squares, nu, logscale = FALSE,
                log.det.sumsq = log(det(sum.of.squares)))

InverseWishartPrior(variance.guess, variance.guess.weight)

Arguments

Details

The inverse Wishart distribution has density function

\frac{|Sigma|^{-\frac{\nu + p + 1}{2}} \exp(-tr(\Sigma^{-1}S) / 2)}{ 2^{\frac{\nu p}{2}}|\Sigma|^{\frac{\nu}{2}}\Gamma_p(\nu / 2)}%

Value

dInverseWishart returns the scalar density (or log density) at the specified value. This function is not vectorized, so only one random variable (matrix) can be evaluated at a time.

InverseWishartPrior returns a list that encodes the parameters of the distribution in a format expected by underlying C++ code.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

See Also

dWishart, rWishart, NormalInverseWishartPrior

Inverse Gamma Distribution

Description

Density, distribution function, quantile function, and random draws from the inverse gamma distribution.

Usage

dinvgamma(x, shape, rate, logscale = FALSE)
pinvgamma(x, shape, rate, lower.tail = TRUE, logscale = FALSE)
qinvgamma(p, shape, rate, lower.tail = TRUE, logscale = FALSE)
rinvgamma(n, shape, rate)

Arguments

Value

• rinvgamma returns draws from the distribution.

• dinvgamma returns the density function.

• pinvgamma returns the cumulative distribution function (or survivor function, if lower.tail == FALSE).

• qinvgamma returns quantiles from the distribution. qinvgamma and pinvgamma are inverse functions.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Check whether a number is even or odd.

Description

Check whether a number is even or odd.

Usage

IsEven(x)
IsOdd(x)

Arguments

Value

Logical indicating whether the argument is even (or odd).

Author(s)

Steven L. Scott

Examples

IsEven(2) ## TRUE
IsOdd(2)  ## FALSE

Log Multivariate Gamma Function

Description

Returns the log of the multivariate gamma function.

Usage

lmgamma(y, dimension)

Arguments

Details

The multivariate gamma function is

\Gamma_p(y) = \pi^{p(p-1)/4} \prod_{j =1}^p \Gamma(y + (1-j)/2).

The multivariate gamma function shows up as part of the normalizing constant for the Wishart and inverse Wishart distributions.

Value

Returns the log of the multivariate gamma function. Note that this function is not vectorized. Both y and dimension must be scalars, and the return value is a scalar.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Log Integrated Gaussian Likelihood

Description

Compute the log of the integrated Gaussian likelihood, where the model paramaters are integrated out with respect to a normal-inverse gamma prior.

Usage

LogIntegratedGaussianLikelihood(suf, prior)

Arguments

Value

Returns a scalar giving the log integrated likelihood.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

prior <- NormalInverseGammaPrior(10, 2, sqrt(2), 1)
y <- c(7.8949, 9.17438, 8.38808, 5.52521)
suf <- GaussianSuf(y)
loglike <- LogIntegratedGaussianLikelihood(suf, prior)
# -9.73975

Lognormal Prior Distribution

Description

Specifies a lognormal prior distribution.

Usage

  LognormalPrior(mu = 0.0, sigma = 1.0, initial.value = NULL)

Arguments

Details

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.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

https://en.wikipedia.org/wiki/Log-normal_distribution

Prior for a Markov chain

Description

The conjugate prior distribution for the parameters of a homogeneous Markov chain. The rows in the transition probability matrix modeled with independent Dirichlet priors. The distribution of the initial state is modeled with its own independent Dirichlet prior.

Usage

MarkovPrior(prior.transition.counts = NULL,
            prior.initial.state.counts = NULL,
            state.space.size = NULL,
            uniform.prior.value = 1)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

MatchDataFrame

Description

Given two data frames with the same data, but with rows and columns in potentially different orders, produce a pair of permutations such that data2[row.permutation, column.permutation] matches data1.

Usage

MatchDataFrame(data.to.match, data.to.permute)

Arguments

Value

Returns a list with two elements.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

 x1 <- data.frame(larry = rnorm(10), moe = 1:10, curly = rpois(10, 2))
 x2 <- x1[c(1:5, 10:6), c(3, 1, 2)]

 m <- MatchDataFrame(x1, x2)
 x2[m$row.permutation, m$column.permutation] == x1  ## all TRUE

Scan a Matrix

Description

Quickly scan a matrix from a file.

Usage

  mscan(fname, nc = 0, header = FALSE, burn = 0, thin = 0, nlines = 0L,
        sep = "", ...)

Arguments

Details

This function is similar to read.table, but scanning a matrix of homogeneous data is much faster because there is much less format deduction.

Value

The matrix stored in the data file.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

filename <- file.path(tempdir(), "example.data")
cat("foo bar baz", "1 2 3", "4 5 6", file = filename, sep = "\n")
m <- mscan(filename, header = TRUE)
m
##      foo bar baz
## [1,]   1   2   3
## [2,]   4   5   6

diagonal MVN prior

Description

A multivariate normal prior distribution formed by the product of independent normal margins.

Usage

MvnDiagonalPrior(mean.vector, sd.vector)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Independence prior for the MVN

Description

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.

Usage

MvnIndependentSigmaPrior(mvn.prior, sd.prior.list)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Multivariate normal prior

Description

A multivariate normal prior distribution.

Usage

MvnPrior(mean, variance)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Normal inverse gamma prior

Description

The NormalInverseGammaPrior is the conjugate prior for the mean and variance of the scalar normal distribution. The model says that

\frac{1}{\sigma^2} \sim Gamma(df / 2, ss/2) \mu|\sigma \sim N(\mu_0, \sigma^2/\kappa)

Usage

NormalInverseGammaPrior(mu.guess, mu.guess.weight = .01,
       sigma.guess, sigma.guess.weight = 1, ...)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Normal inverse Wishart prior

Description

The NormalInverseWishartPrior is the conjugate prior for the mean and variance of the multivariate normal distribution. The model says that

\Sigma^{-1} \sim Wishart(\nu, S) \mu|\sigma \sim N(\mu_0, \Sigma/\kappa)

The Wishart(S, \nu) distribution is parameterized by S, the inverse of the sum of squares matrix, and the scalar degrees of freedom parameter nu.

The distribution is improper if \nu < dim(S).

Usage

NormalInverseWishartPrior(mean.guess,
                          mean.guess.weight = .01,
                          variance.guess,
                          variance.guess.weight = nrow(variance.guess) + 1)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Normal (scalar Gaussian) prior distribution

Description

Specifies a scalar Gaussian prior distribution.

Usage

NormalPrior(mu, sigma, initial.value = mu, fixed = FALSE)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Pairs plot for posterior distributions.

Description

A pairs plot showing the posterior distribution of the given list of Monte Carlo draws. Plots above the diagonal show the posterior distribution on a scale just wide enough to fit the plots. The diagonal shows a marginal density plot, and the subdiagonal shows the distribution with all plots on a common scale.

Usage

PairsDensity(draws,
             nlevels = 20,
             lty = NULL,
             color = NULL,
             subset = NULL,
             labels,
             legend.location = "top",
             legend.cex = 1,
             label.cex = 1,
             ...)

Arguments

Author(s)

Steven L. Scott

See Also

pairs, CompareDensities, CompareManyDensities

Examples

## You can see the pairs plot for a single set of draws.
y <- matrix(rnorm(5000, mean = 1:5), ncol = 5, byrow = TRUE)
PairsDensity(y)

## You can also compare two or more sets of draws.
z <- matrix(rnorm(2500, mean = 2:6), ncol = 5, byrow = TRUE)
PairsDensity(list("first set" = y, "second set" = z))

Contour plot of a bivariate density.

Description

Contour plot of one ore more bivariate densities. This function was originally created to implement PairsDensity, but might be useful on its own.

Usage

PlotDensityContours(draws,
                    x.index = 1,
                    y.index = 2,
                    xlim = NULL,
                    ylim = NULL,
                    nlevels = 20,
                    subset = NULL,
                    color = NULL,
                    lty = NULL,
                    axes = TRUE,
                    ...)

Arguments

Author(s)

Steven L. Scott

See Also

contour, kde2d

Examples

## You can see the pairs plot for a single set of draws.
y <- matrix(rnorm(5000, mean = 1:5), ncol = 5, byrow = TRUE)
PlotDensityContours(y, 3, 1)

## You can also compare two or more sets of draws.
z <- matrix(rnorm(2500, mean = 2:6), ncol = 5, byrow = TRUE)
PlotDensityContours(list("first set" = y, "second set" = z), 3, 1)

Plots the pointwise evolution of a distribution over an index set.

Description

Produces an dynamic distribution plot where gray scale shading is used to show the evolution of a distribution over an index set. This function is particularly useful when the index set is too large to do side-by-side boxplots.

Usage

PlotDynamicDistribution(curves,
                        timestamps = NULL,
                        quantile.step=.01,
                        xlim = NULL,
                        xlab = "Time",
                        ylim = range(curves, na.rm = TRUE),
                        ylab = "distribution",
                        add = FALSE,
                        axes = TRUE,
                        ...)

Arguments

Details

The function works by passing many calls to polygon. Each polygon is associated with a quantile level, with darker shading near the median.

Value

This function is called for its side effect, which is to produce a plot on the current graphics device.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  x <- t(matrix(rnorm(1000 * 100, 1:100, 1:100), nrow=100))
  ## x has 1000 rows, and 100 columns.  Column i is N(i, i^2) noise.

  PlotDynamicDistribution(x)
  time <- as.Date("2010-01-01", format = "%Y-%m-%d") + (0:99 - 50)*7
  PlotDynamicDistribution(x, time)

Plots individual autocorrelation functions for many-valued time series

Description

Produces individual autocorrelation functions for many-valued time series such as those produced by highly multivariate MCMC output. Cross-correlations such as those produced by acf are not shown.

Usage

PlotMacf(x, lag.max = 40, gap = 0.5, main = NULL, boxes = TRUE,
         xlab = "lag", ylab = "ACF", type = "h")

Arguments

Value

Called for its side effect

Author(s)

Steven L. Scott

See Also

acf, plot.many.ts.

Examples

x <- matrix(rnorm(1000), ncol=10)
PlotMacf(x)

Multiple time series plots

Description

Plots many time series plots on the same graphical device. Each plot gets its own frame. Scales can be adjusted to see variation in each plot (each plot gets its own scale), or variation between plots (common scale).

Usage

PlotManyTs(x, type = "l", gap = 0, boxes = TRUE, truth = NULL,
           thin = 1, labs, same.scale = TRUE, ylim = NULL,
           refline = NULL, color = NULL, ...)

Arguments

Author(s)

Steven L. Scott

See Also

plot.ts (for plotting a small number of time series) plot.macf

Examples

x <- matrix(rnorm(1000), ncol = 10)
PlotManyTs(x)
PlotManyTs(x, same = FALSE)

Regression Coefficient Conjugate Prior

Description

A conjugate prior for regression coefficients, conditional on residual variance, and sample size.

Usage

  RegressionCoefficientConjugatePrior(
    mean,
    sample.size,
    additional.prior.precision = numeric(0),
    diagonal.weight = 0)

Arguments

Details

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 | \sigma^2, X \sim % N(b, \sigma^2 (\Lambda^{-1} + V

where

V^{-1} = \frac{\kappa}{n} ((1 - w) X^TX + w diag(X^TX)) .

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).

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Repeated Lists of Objects

Description

Produces repeated copies of an object.

Usage

RepList(object, times)

Arguments

Value

Returns a list containing times copies of object.

Author(s)

Steven L. Scott

Examples

  alist <- list(x = "foo", y = 12, z = c(1:3))
  three.copies <- RepList(alist, 3)

Multivariate Normal Simulation

Description

Simulate draws from the multivariate normal distribution.

Usage

rmvn(n = 1, mu, sigma = diag(rep(1., length(mu))))

Arguments

Details

Note that mu and sigma are the same for all n draws. This function cannot handle separate parameters for each draw the way rnorm and similar functions for scalar random variables can.

Value

If n == 1 the return value is a vector. Otherwise it is a matrix with n rows and length(mu) columns.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

y1 <- rnorm(1, 1:3)
## y1 is a vector
y2 <- rnorm(10, 1:3)
## y2 is a matrix

RVectorFunction

Description

A wrapper for passing R functions to C++ code.

Usage

  RVectorFunction(f, ...)

Arguments

Details

The Boom library can handle the output of this function as a C++ function object. Note that evaluating R functions in C is no faster than evaluating them in R, but a wrapper like this is useful for interacting with C and C++ libraries that expect to operate on C and C++ functions.

Value

A list containing the information needed to evaluate the function f in C++ code.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Scaled Matrix-Normal Prior

Description

A matrix-normal prior distribution, intended as the conjugate prior for the regression coefficients in a multivariate linear regression.

Usage

ScaledMatrixNormalPrior(mean, nu)

Arguments

Details

The matrix normal distribution is a 3-parameter distribution MN(mu, Omega, V), where mu is the mean. A deviate from the distribution is a matrix B, where Cov(B[i, j], B[k, m]) = Omega[i, k] * Sigma[j, m]. If b = Vec(B) is the vector obtained by stacking columns of B, then b is multivariate normal with mean Vec(mu) and covariance matrix

\Sigma \otimes Omega

(the kronecker product).

This prior distribution assumes the underlying C++ code knows where to find the predictor (X) matrix in the regression, and the residual variance matrix Sigma. The assumed prior distribution is B ~ MN(mu, X'X / nu, Sigma).

Like most other priors in Boom, this function merely encodes information expected by the underlying C++ code, ensuring correct names and formatting.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Prior for a standard deviation or variance

Description

Specifies an inverse Gamma prior for a variance parameter, but inputs are defined in terms of a standard deviation.

Usage

  SdPrior(sigma.guess, sample.size = .01, initial.value = sigma.guess,
          fixed = FALSE, upper.limit = Inf)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Sufficient Statistics

Description

Sufficient statistics for various models.

Usage

  RegressionSuf(X = NULL,
                y = NULL,
                xtx = crossprod(X),
                xty = crossprod(X, y),
                yty = sum(y^2),
                n = length(y),
                xbar = colMeans(X),
                ybar = mean(y))

  GaussianSuf(y)

Arguments

Value

The returned value is a function containing the sufficient statistics for a regression model. Arguments are checked to ensure they have legal values. List names match the names expected by underlying C++ code.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

  X <- cbind(1, matrix(rnorm(3 * 100), ncol = 3))
  y <- rnorm(100)

  ## Sufficient statistics can be computed from raw data, if it is
  ## available.
  suf1 <- RegressionSuf(X, y)

  ## The individual components can also be computed elsewhere, and
  ## provided as arguments.  If n is very large, this can be a
  ## substantial coomputational savings.
  suf2 <- RegressionSuf(xtx = crossprod(X),
                        xty = crossprod(X, y),
                        yty = sum(y^2),
                        n = 100,
                        xbar = colMeans(X))

Suggest MCMC Burn-in from Log Likelihood

Description

Suggests a burn-in period for an MCMC chain based on the log likelihood values simulated in the final few iterations of the chain.

Usage

  SuggestBurnLogLikelihood(log.likelihood, fraction = .10, quantile = .9)

Arguments

Details

Looks at the last fraction of the log.likelihood sequence and finds a specified quantile to use as a threshold. Returns the first iteration where log.likelihood exceeds this threshold.

Value

Returns a suggested number of iterations to discard. This can be 0 if fraction == 0, which is viewed as a signal that no burn-in is desired.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Thin the rows of a matrix

Description

Systematic sampling of every thin'th row of a matrix or vector. Useful for culling MCMC output or denoising a plot.

Usage

thin(x, thin)

Arguments

Value

The thinned vector, matrix, or array is returned.

Author(s)

Steven L. Scott

Examples

x <- rnorm(100)
thin(x, 10)
# returns a 10 vector

y <- matrix(rnorm(200), ncol=2)
thin(y, 10)
# returns a 10 by 2 matrix

Thin a Matrix

Description

Return discard all but every k'th row of a matrix.

Usage

ThinMatrix(mat, thin)

Arguments

Details

The bigger the value of thin the more thinning that gets done. For example, thin = 10 will keep every 10 lines from mat.

Value

The matrix mat, after discarding all but every thin lines.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Examples

m <- matrix(1:100, ncol = 2)
ThinMatrix(m, thin = 10)
##      [,1] [,2]
## [1,]   10   60
## [2,]   20   70
## [3,]   30   80
## [4,]   40   90
## [5,]   50  100

Trace of the Product of Two Matrices

Description

Returns the trace of the product of two matrices.

Usage

TraceProduct(A, B, b.is.symmetric = FALSE)

Arguments

Value

Returns a number equivalent to sum(diag(A %*% B)).

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

Uniform prior distribution

Description

Specifies a uniform prior distribution on a real-valued scalar parameter.

Usage

  UniformPrior(lo = 0, hi = 1, initial.value = NULL)

Arguments

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com

References

Gelman, Carlin, Stern, Rubin (2003), "Bayesian Data Analysis", Chapman and Hall.

Wishart Distribution

Description

Density and random generation for the Wishart distribution.

Usage

dWishart(W, Sigma, nu, logscale = FALSE)
rWishart(nu, scale.matrix, inverse = FALSE)

Arguments

Details

If nu is an integer then a W(\Sigma, \nu) random variable can be thought of as the sum of nu outer products: y_iy_i^T, where y_i is a zero-mean multivariate normal with variance matrix Sigma.

The Wishart distribution is

\frac{|W|^{\frac{\nu - p - 1}{2}} \exp(-tr(\Sigma^{-1}W) / 2)}{ 2^{\frac{\nu p}{2}}|\Sigma|^{\frac{\nu}{2}}\Gamma_p(\nu / 2)}%

where p == nrow(W) and \Gamma_p(\nu) is the multivariate gamma function (see lmgamma).

Value

dWishart returns the density of the Wishart distribution. It is not vectorized, so only one random variable (matrix) can be evaluated at a time.

rWishart returns one or more draws from the Wishart or inverse Wishart distributions. If n > 0 the result is a 3-way array. Unlike the rWishart function from the stats package, the first index corresponds to draws. This is in keeping with the convention of other models from the Boom package.

Author(s)

Steven L. Scott steve.the.bayesian@gmail.com