Wishart {mniw}

Wishart and Inverse-Wishart distributions.

Description

Densities and random sampling for the Wishart and Inverse-Wishart distributions.

Usage

dwish(X, Psi, nu, log = FALSE)

rwish(n, Psi, nu)

diwish(X, Psi, nu, log = FALSE)

riwish(n, Psi, nu)

dwishart(X, Psi, nu, inverse = FALSE, log = FALSE)

rwishart(n, Psi, nu, inverse = FALSE)

Arguments

X

Argument to the density function. Either a q x q matrix or a q x q x n array.

Psi

Scale parameter. Either a q x q matrix or a q x q x n array.

nu

Degrees-of-freedom parameter. A scalar or vector of length n.

log

Logical; whether or not to compute the log-density.

n

Integer number of random samples to generate.

inverse

Logical; whether or not to use the Inverse-Wishart distribution.

Details

The Wishart distribution \boldsymbol{X} \sim \textrm{Wishart}(\boldsymbol{\Psi}, \nu) on a symmetric positive-definite random matrix \boldsymbol{X} of size q \times q has PDF

f(\boldsymbol{X} \mid \boldsymbol{\Psi}, \nu) = \frac{|\boldsymbol{X}|^{(\nu-q-1)/2}\exp\big\{-\textrm{tr}(\boldsymbol{\Psi}^{-1}\boldsymbol{X})/2\big\}}{2^{q\nu/2}|\boldsymbol{\Psi}|^{\nu/2} \Gamma_q(\frac \nu 2)},

where \Gamma_q(\alpha) is the multivariate gamma function,

\Gamma_q(\alpha) = \pi^{q(q-1)/4} \prod_{i=1}^q \Gamma(\alpha + (1-i)/2).

The Inverse-Wishart distribution \boldsymbol{X} \sim \textrm{Inverse-Wishart}(\boldsymbol{\Psi}, \nu) is defined as \boldsymbol{X}^{-1} \sim \textrm{Wishart}(\boldsymbol{\Psi}^{-1}, \nu).

dwish and diwish are convenience wrappers for dwishart, and similarly rwish and riwish are wrappers for rwishart.

Value

A vector length n for density evaluation, or an array of size q x q x n for random sampling.

Examples

# Random sampling

n <- 1e5
q <- 3
Psi1 <- crossprod(matrix(rnorm(q^2),q,q))
nu <- q + runif(1, 0, 5)

X1 <- rwish(n,Psi1,nu) # Wishart

# plot it
plot_fun <- function(X) {
  q <- dim(X)[1]
  par(mfrow = c(q,q))
  for(ii in 1:q) {
    for(jj in 1:q) {
      hist(X[ii,jj,], breaks = 100, freq = FALSE,
           xlab = "", main = parse(text = paste0("X[", ii, jj, "]")))
    }
  }
}

plot_fun(X1)

# "vectorized" scale parameeter
Psi2 <- 5 * Psi1
vPsi <- array(c(Psi1, Psi2), dim = c(q, q, n))
X2 <- rwish(n, Psi = vPsi, nu = nu)
plot_fun(X2)

# Inverse-Wishart
X3 <- riwish(n, Psi2, nu)
plot_fun(X3)

# log-density calculation for sampled values
par(mfrow = c(1,1))
hist(dwish(X2, vPsi, nu, log = TRUE),
     breaks = 100, freq = FALSE, xlab = "",
     main = expression("log-p"*(X[2]*" | "*list(Psi,nu))))