Random draws from a Normal-Inverse-Wishart distribution.
Generates random draws from a Normal-Inverse-Wishart (NIW) distribution. Can be used to compare prior to posterior parameter distributions.
rniw(n, lambda, kappa, Psi, nu)
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. |
The NIW distribution p(μ, Σ | λ, κ, Ψ, ν) is defined as
Σ \sim W^{-1}(Ψ, ν), \quad μ | Σ \sim N(λ, Σ/κ).
Returns a list with elements mu and Sigma of sizes c(n,length(lambda)) and c(nrow(Psi),ncol(Psi),n).
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")Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.