niiw.post: Random draws from the posterior distribution with Normal-Independent-Inverse-Wishart (NIIW) prior.

Description

Given iid $d$-dimensional niche indicators $X = (X_1,\ldots,X_N)$ with $X_i \sim N(\mu, \Sigma)$, this function generates random draws from $p(\mu,\Sigma | X)$ for the Normal-Independent-Inverse-Wishart (NIIW) prior.

Usage

niiw.post(nsamples, X, lambda, Omega, Psi, nu, mu0 = lambda, burn)

Arguments

nsamples

the number of posterior draws.

a data matrix with observations along the rows.

lambda

mean of mu. See Details.

Omega

variance of mu. Defaults to Omega = 0. See Details.

Psi

scale matrix of Sigma. Defaults to Psi = 0. See Details.

degrees of freedom of Sigma. Defaults to nu = ncol(X)+1. See Details.

mu0

initial value of mu to start the Gibbs sampler. See Details.

burn

burn-in for the MCMC sampling algorithm. Either an integer giving the number of initial samples to discard, or a fraction with 0 < burn < 1. Defaults to burn = floor(nsamples/10).

Value

Returns a list with elements mu and Sigma of sizes c(nsamples, length(lambda)) and c(dim(Psi), nsamples).

Details

The NIIW distribution $p(\mu, \Sigma | \lambda, \kappa, \Psi, \nu)$ is defined as $$\Sigma \sim W^{-1}(\Psi, \nu), \quad \mu | \Sigma \sim N(\lambda, \Omega).$$ The default value Omega = 0 uses the Lebesque prior on $\mu$: $p(\mu) \propto 1$. In this case the NIW and NIIW priors produce identical resuls, but niw.post is faster. The default value Psi = 0 uses the scale-invariant prior on $\Sigma$: $p(\Sigma) \propto |\Sigma|^{-(\nu+d+1)/2}$. The default value nu = ncol(X)+1 for Omega = 0 and Psi = 0 makes $E[\mu|X]=\code{colMeans(X)}$ and $E[\Sigma | X]=\code{var(X)}$. Random draws are obtained by a Markov chain Monte Carlo (MCMC) algorithm; specifically, a Gibbs sampler alternates between draws from $p(\mu | \Sigma, X)$ and $p(\Sigma | \mu, X)$, which are Normal and Inverse-Wishart distributions respectively.

Examples

Run this code

# NOT RUN {
# simulate data
d <- 4
mu0 <- rnorm(d)
Sigma0 <- matrix(rnorm(d^2), d, d)
Sigma0 <- Sigma0 %*% t(Sigma0)
N <- 100
X <- rmvnorm(N, mean = mu0, sigma = Sigma0)

# prior parameters
# flat prior on mu
lambda <- 0
Omega <- 0
# informative prior on Sigma
Psi <- crossprod(matrix(rnorm(d^2), d, d))
nu <- 5

# sample from NIIW posterior
nsamples <- 2e3
system.time({
 siiw <- niiw.post(nsamples, X, lambda, Omega, Psi, nu, burn = 100)
})

# sample from NIW posterior
kappa <- 0
system.time({
 siw <- niw.post(nsamples, X, lambda, kappa, Psi, nu)
})

# check that posteriors are the same

# p(mu | X)
clrs <- c("black", "red")
par(mar = c(4.2, 4.2, 2, 1)+.1)
niche.par.plot(list(siiw, siw), col = clrs, plot.mu = TRUE, plot.Sigma = FALSE)
legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs)

# p(Sigma | X)
par(mar = c(4.2, 4.2, 2, 1)+.1)
niche.par.plot(list(siiw, siw), col = clrs, plot.mu = FALSE, plot.Sigma = TRUE)
legend(x = "topright", legend = c("NIIW Prior", "NIW Prior"), fill = clrs)
# }