# Random draws from the posterior distribution with Normal-Independent-Inverse-Wishart (NIIW) prior. 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. ```r niiw.post(nsamples, X, lambda, Omega, Psi, nu, mu0 = lambda, burn) ``` ## Arguments - `nsamples`: The number of posterior draws. - `X`: 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'. - `nu`: 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)`. ## Returns 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]=`colMeans(X)`$ and $E[\Sigma | X]=`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 ```r # simulate normal data with mean and variance (mu0, Sigma0) d <- 4 mu0 <- rnorm(d) Sigma0 <- matrix(rnorm(d^2), d, d) Sigma0 <- Sigma0 %*% t(Sigma0) N <- 1e2 X <- matrix(rnorm(N*d), N, d) # iid N(0,1) X <- t(t(X %*% chol(Sigma0)) + mu0) # each row is N(mu0, Sigma) # 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) ``` ## See Also `niw.post()`, `rwish()`.

niiw.post function