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