rniw: Random draws from a Normal-Inverse-Wishart distribution.

Description

Generates random draws from a Normal-Inverse-Wishart (NIW) distribution. Can be used to compare prior to posterior parameter distributions.

Usage

rniw(n, lambda, kappa, Psi, nu)

Value

Returns a list with elements mu and Sigma of sizes c(n,length(lambda)) and c(nrow(Psi),ncol(Psi),n).

Arguments

n: Number of samples to draw.
lambda: Location parameter. See 'Details'.
kappa: Scale parameter. See 'Details'.
Psi: Scale matrix. See 'Details'.
nu: Degrees of freedom. See 'Details'.

Details

The NIW distribution $p(\mu, \Sigma | \lambda, \kappa, \Psi, \nu)$ is defined as $$ \Sigma \sim W^{-1}(\Psi, \nu), \quad \mu | \Sigma \sim N(\lambda, \Sigma/\kappa). $$

Examples

Run this code

d <- 4 # number of dimensions
nu <- 7 # degrees of freedom
Psi <- crossprod(matrix(rnorm(d^2), d, d)) # scale
lambda <- rnorm(d)
kappa <- 2
n <- 1e4

niw.sim <- rniw(n, lambda, kappa, Psi, nu)

# diagonal elements of Sigma^{-1} are const * chi^2
S <- apply(niw.sim$Sigma, 3, function(M) diag(solve(M)))

ii <- 2
const <- solve(Psi)[ii,ii]
hist(S[ii,], breaks = 100, freq = FALSE,
     main = parse(text = paste0("\"Histogram of \"*(Sigma^{-1})[", ii,ii,"]")),
     xlab = parse(text = paste0("(Sigma^{-1})[", ii,ii,"]")))
curve(dchisq(x/const, df = nu)/const,
      from = min(S[ii,]), to = max(S[ii,]), col = "red", add = TRUE)

# elements of mu have a t-distribution
mu <- niw.sim$mu

ii <- 4
const <- sqrt(Psi[ii,ii]/(kappa*(nu-d+1)))
hist(mu[,ii], breaks = 100, freq = FALSE,
     main = parse(text = paste0("\"Histogram of \"*mu[", ii, "]")),
     xlab = parse(text = paste0("mu[", ii, "]")))
curve(dt((x-lambda[ii])/const, df = nu-d+1)/const, add = TRUE, col = "red")