Wishart: Wishart and Inverse-Wishart distributions.

In mniw: The Matrix-Normal Inverse-Wishart Distribution

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 X ~ Wishart(Ψ, ν) on a symmetric positive-definite random matrix X of size q x q has PDF

f(X | Ψ, ν) = [ |X|^((ν-q-1)/2) * e^(-trace(Ψ^{-1} X)/2) ] / [ 2^(qν/2) * |Ψ|^(ν/2) * Γ_q(ν/2) ],

where Γ_q(α) is the multivariate gamma function,

Γ_q(α) = π^(q(q-1)/4) ∏_{i=1}^q Γ(α + (1-i)/2).

The Inverse-Wishart distribution X ~ Inverse-Wishart(Ψ, ν) is defined as X^{-1} ~ Wishart(Ψ^{-1}, ν).

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

mniw documentation built on Aug. 22, 2022, 5:05 p.m.