# niiw.post: Random draws from the posterior distribution with... In nicheROVER: Niche Region and Niche Overlap Metrics for Multidimensional Ecological Niches

## Description

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

## Usage

 1 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 μ. See 'Details'. Omega Variance of μ. Defaults to Omega = 0. See 'Details'. Psi Scale matrix of Σ. Defaults to Psi = 0. See 'Details'. nu Degrees of freedom of Σ. Defaults to nu = ncol(X)+1. See 'Details'. mu0 Initial value of μ 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).

## Details

The NIIW distribution p(μ, Σ | λ, κ, Ψ, ν) is defined as

Σ \sim W^{-1}(Ψ, ν), \quad μ | Σ \sim N(λ, Ω).

The default value Omega = 0 uses the Lebesque prior on μ: p(μ) \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 Σ: p(Σ) \propto |Σ|^{-(ν+d+1)/2}.

The default value nu = ncol(X)+1 for Omega = 0 and Psi = 0 makes E[μ|X]=colMeans(X) and E[Σ | X]=var(X).

Random draws are obtained by a Markov chain Monte Carlo (MCMC) algorithm; specifically, a Gibbs sampler alternates between draws from p(μ | Σ, X) and p(Σ | μ, X), which are Normal and Inverse-Wishart distributions respectively.

## Value

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

niw.post(), rwish().
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 # 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)