R/densities.R In CholWishart: Cholesky Decomposition of the Wishart Distribution

Documented in dInvWishartdWishart

#   densities.R
#   CholWishart: Sample the Cholesky Factor of the Wishart and Other Functions
#   Copyright (C) 2018  GZ Thompson <gzthompson@gmail.com>
#
#   This program is free software; you can redistribute it and/or modify
#   the Free Software Foundation; either version 3 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#     You should have received a copy of the GNU General Public License
#   along with this program; if not, a copy is available at

#' Density for Random Wishart Distributed Matrices
#'
#' Compute the density of an observation of a random Wishart distributed matrix
#' (\code{dWishart}) or an observation
#' from the inverse Wishart distribution (\code{dInvWishart}).
#'
#'    Note there are different ways of parameterizing the Inverse
#'    Wishart distribution, so check which one you need.
#'     Here,  If \eqn{X \sim IW_p(\Sigma, \nu)}{X ~ IW_p(Sigma, df)} then
#'     \eqn{X^{-1} \sim W_p(\Sigma^{-1}, \nu)}{X^{-1} ~ W_p(Sigma^{-1}, df)}.
#'     Dawid (1981) has a different definition: if
#'     \eqn{X \sim W_p(\Sigma^{-1}, \nu)}{X ~ W_p(Sigma^{-1}, df)} and
#'     \eqn{\nu > p - 1}{df > p - 1}, then
#'     \eqn{X^{-1} = Y \sim IW(\Sigma, \delta)}{X^{-1} = Y ~ IW(Sigma, delta)},
#'     where
#'     \eqn{\delta = \nu - p + 1}{delta = df - p + 1}.
#'
#' @param x positive definite \eqn{p \times p}{p * p} observations for density
#'          estimation - either one matrix or a 3-D array.
#' @param df numeric parameter, "degrees of freedom".
#' @param Sigma positive definite \eqn{p \times p}{p * p} "scale" matrix,
#'              the matrix parameter of the distribution.
#' @param log logical, whether to return value on the log scale.
#'
#' @return Density or log of density
#'
#' @references
#' Dawid, A. (1981). Some Matrix-Variate Distribution Theory:
#' Notational Considerations and a Bayesian Application.
#' \emph{Biometrika}, 68(1), 265-274. \doi{10.2307/2335827}
#'
#' Gupta, A. K.  and D. K. Nagar (1999). \emph{Matrix variate distributions}.
#' Chapman and Hall.
#'
#' Mardia, K. V., J. T. Kent, and J. M. Bibby (1979)
#' \emph{Multivariate Analysis},
#' @export
#'
#' @examples
#' set.seed(20180222)
#' A <- rWishart(1, 10, diag(4))[, , 1]
#' A
#' dWishart(x = A, df = 10, Sigma = diag(4L), log = TRUE)
#' dInvWishart(x = solve(A), df = 10, Sigma = diag(4L), log = TRUE)
dWishart <- function(x, df, Sigma, log = TRUE) {
if (!is.numeric(Sigma)) {
stop("'Sigma' must be numeric")
}
sigma <- as.matrix(Sigma)
dims <- dim(sigma)
if (!is.numeric(x)) {
stop("'x' must be numeric")
}
if (dim(x)[1L] != dim(x)[2L] ||
dim(x)[1L] != dims[1L] || dims[1L] != dims[2L]) {
stop("non-conformable dimensions")
}
if (!isSymmetric(sigma)) {
stop("non-symmetric input")
}
if (length(dim(x)) < 3L) x <- array(x, dim = c(dim(x), 1L))
dimx <- dim(x)
ldensity <- rep(0, dimx[3L])
chol_s <- chol(sigma)
ldet_s <- sum(log(diag(chol_s))) * 2
for (i in seq(dimx)) {
if (!isSymmetric(x[, , i])) {
stop("non-symmetric input")
}
chol_x <- chol(x[, , i])
ldet_x <- sum(log(diag(chol_x))) * 2
ldensity[i] <-
.5 * (df - dims[1L] - 1) * ldet_x +
-.5 * sum(diag(chol2inv(chol_s) %*% x[, , i])) +
-(df * dims[1L] / 2 * log(2)) - .5 * df * ldet_s +
-lmvgamma(df / 2, dims[1L])
}
if (log) {
return(ldensity)
} else {
return(exp(ldensity))
}
}

#' @describeIn dWishart density for the inverse Wishart distribution.
#' @export
dInvWishart <- function(x, df, Sigma, log = TRUE) {
if (!is.numeric(Sigma)) {
stop("'Sigma' must be numeric")
}
sigma <- as.matrix(Sigma)
dims <- dim(sigma)
if (!is.numeric(x)) {
stop("x must be numeric.")
}
if (dim(x)[1L] != dim(x)[2L] ||
dim(x)[1L] != dims[1L] || dims[1L] != dims[2L]) {
stop("non-conformable dimensions")
}
if (!isSymmetric(sigma)) {
stop("non-symmetric input")
}
if (length(dim(x)) < 3L) x <- array(x, dim = c(dim(x), 1L))
dimx <- dim(x)
ldensity <- rep(0, dimx[3L])
chol_s <- chol(sigma)
ldet_s <- sum(log(diag(chol_s))) * 2
for (i in seq(dimx[3L])) {
if (!isSymmetric(x[, , i])) {
stop("non-symmetric input")
}
chol_x <- chol(x[, , i])
ldet_x <- sum(log(diag(chol_x))) * 2
ldensity[i] <-
-.5 * (df + dims[1L] + 1) * ldet_x + .5 * df * ldet_s +
-.5 * sum(diag(chol2inv(chol_x) %*% sigma)) -
(df * dims[1L] / 2 * log(2)) - lmvgamma(df / 2, dims[1L])
}
if (log) {
return(ldensity)
} else {
return(exp(ldensity))
}
}

Try the CholWishart package in your browser

Any scripts or data that you put into this service are public.

CholWishart documentation built on Oct. 8, 2021, 9:09 a.m.