Description Usage Arguments Details Value See Also Examples
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.
1 | niiw.post(nsamples, X, lambda, Omega, Psi, nu, mu0 = lambda, burn)
|
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 |
Psi |
Scale matrix of Σ. Defaults to |
nu |
Degrees of freedom of Σ. Defaults to |
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 |
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.
Returns a list with elements mu
and Sigma
of sizes c(nsamples, length(lambda))
and c(dim(Psi), nsamples)
.
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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.