R/inverse.wishart.R

In Boom: Bayesian Object Oriented Modeling

dInverseWishart <- function(Sigma, sum.of.squares, nu, logscale = FALSE,
                            log.det.sumsq = log(det(sum.of.squares))) {
  ## Inverse Wishart density function
  ## Args:
  ##   Sigma:  The argument at which to evaluate the density.  A variance matrix.
  ##   sum.of.squares: A positive definite matrix.  Typically this is
  ##     the sum of squares that is the sufficient statistic for the
  ##     inverse Wishart distribution.
  ##   nu: A positive scalar giving the "degrees of freedom" parameter
  ##     for the distribution.  This must satisfy nu > nrow(Sigma) - 1.
  ##   logscale: Logical.  If TRUE then the log of the density is
  ##     returned.  Otherwise the density is returned on the absolute
  ##     scale.
  ##   log.det.sumsq: If this function is to be called many times with
  ##     the same value of sum.of.squares, then precomputing
  ##     log(det(sum.of.squares)) can save considerable time.
  ## Returns:
  ##   The scalar value of the density function (or its log) evaluated
  ##   at the argument.
  dimension <- ncol(Sigma)
  stopifnot(nu > dimension - 1)
  stopifnot(nrow(Sigma) == dimension,
            dim(sum.of.squares) == dim(Sigma))
  Sigma.chol <- chol(Sigma)
  logdet.Sigma <- 2 * sum(log(diag(Sigma.chol)))
  if (is.null(log.det.sumsq)) {
    log.det.sumsq <- log(det(sum.of.squares))
  }
  siginv <- chol2inv(Sigma.chol)
  ans <- .5 * nu * log.det.sumsq -
      .5 * (nu + dimension + 1) * logdet.Sigma -
      .5 * TraceProduct(sum.of.squares, siginv) -
      .5 * nu * dimension * log(2) -
      lmgamma(.5 * nu, dimension)
  if (logscale) return(ans)
  return(exp(ans))
}

dWishart <- function(W, Sigma, nu, logscale = FALSE) {
  ## Evaluates the density of the Wishart distribution.
  ##
  ## Args:
  ##   W: The argument of the distribution (i.e. the random variable).
  ##     A positive definite matrix.
  ##   Sigma: The variance of the generating normal distribution.  See
  ##     the Details section below.
  ##   nu: A positive scalar giving the 'degrees.of.freedom' parameter
  ##     of the Wishart distribution.  The distribution is only
  ##     defined for nu >= nrow(W) - 1.
  ##   logscale: Logical.  If TRUE then the log of the density is
  ##     returned.  Otherwise the density is returned on the absolute
  ##     scale.
  ##
  ## Returns:
  ##   The scalar value of the density function (or its log) evaluated
  ##   at the argument.
  ##
  ## Details:
  ##
  ##   If 'nu' is a scalar then a Wishart random variable can be
  ##   thought of as the sum of 'nu' outer products of zero-mean
  ##   multivariate normal deviates, each with the same variance
  ##   matrix Sigma.  The mean of the distribution is 'nu * Sigma'.
  sigma.chol <- chol(Sigma);
  siginv <- chol2inv(sigma.chol)
  sigma.logdet <- 2 * sum(log(diag(sigma.chol)))
  W.logdet <- 2 * sum(log(diag(chol(W))))
  dimension <- nrow(W)
  log2 <- 0.69314718055994528623
  ans <- .5 * (nu - dimension - 1) * W.logdet -
      TraceProduct(siginv, W) / 2 -
      .5 * nu * dimension * log2 -
      .5 * nu * sigma.logdet -
      lmgamma(nu/2, dimension)
  if (logscale) return(ans)
  return(exp(ans))
}

rWishart <- function(nu,  scale.matrix, inverse = FALSE) {
  ## Generates a single draw from the Wishart or inverse Wishart
  ## distribution.  The Wishart distribution is generated by taking
  ## 'nu' draws from the MVN(0, Sigma) and adding their outer
  ## products.
  ##
  ## The inverse Wishart distribution is the conjugate prior for Sigma
  ## in a multivariate normal model.
  ##
  ## Args:
  ##   nu: The degrees of freedom parameter in the Wishart or inverse
  ##     Wishart distribution.  The distribution requires nu > dimension - 1.
  ##   scale.matrix: For the Wishart, 'scale.matrix' is the 'Sigma'
  ##     parameter of the generating multivariate normal distribution.
  ##     For the inverse Wishart, scale.matrix is the INVERSE of the
  ##     "sum of squares" matrix portion of the Mvn sufficient
  ##     statistics.
  ##   inverse: Logical.  If TRUE then simulate from the Wishart
  ##     distribution.  Otherwise simulate from the inverse: Wishart
  ##     distribution.
  dimension = nrow(scale.matrix)
  R <- matrix(0, nrow = dimension, ncol = dimension)
  R[upper.tri(R)] <- rnorm(dimension * (dimension - 1) / 2)
  diag(R) <- sqrt(rchisq(dimension, df = (nu + 1) - (1:dimension)))
  R <- R %*% chol(scale.matrix)
  if (!inverse) return(t(R) %*% R)  # draw of siginv
  return(chol2inv(R))
}

TraceProduct <- function(A, B, b.is.symmetric = FALSE) {
  ## Returns the trace of the product of the two matrices A and B.
  ## Args:
  ##   A:  The matrix on the left hand side of the product.
  ##   B:  The matrix on the right hand side of the product.
  ##   b.is.symmetric: Logical.  A TRUE value indicates that B is a
  ##     symmetric matrix.  A slight computational savings is possible
  ##     if B is symmetric.
  ## Returns:
  ##   A scalar value giving trace(A %*% B).
  stopifnot(is.matrix(A))
  stopifnot(is.matrix(B))

  ## The trace of AB is sum_i AB(i, i).
  ## AB(i, i) is sum_j A(i, j) * B(j, i).
  ## Therefore, trace(AB) = sum_i sum_j A(i, j) * B(j, i).  This is
  ## the elementwise product of A times B.transpose.
  if (b.is.symmetric) {
    return(sum(A * B))
  } else {
    return(sum(A * t(B)))
  }
}