powerprior (version 1.0.0)

sample_posterior_multivariate: Sample from Posterior Distribution (Multivariate)

Description

Generates samples from the multivariate posterior distribution using exact closed-form expressions from the Normal-Inverse-Wishart conjugate family.

Usage

sample_posterior_multivariate(posterior, n_samples = 1000, marginal = FALSE)

Value

If marginal=FALSE, a list with components:

mu

n_samples x p matrix of mean samples

Sigma

p x p x n_samples array of covariance samples

If marginal=TRUE, an n_samples x p matrix of mean samples.

Arguments

posterior

Object of class "posterior_multivariate" from posterior_multivariate()

n_samples

Number of samples to generate (default: 1000). For large n_samples or high dimensions, computation time increases.

marginal

Logical. If TRUE, samples only \(\mu\) from the multivariate t-distribution. If FALSE (default), samples the joint (\(\mu\), \(\Sigma\)) from the NIW distribution, which is more computationally intensive but provides uncertainty in the covariance structure.

Details

Sampling Algorithms

Joint Sampling (marginal=FALSE):

Implements the standard hierarchical sampling algorithm for the NIW distribution:

1. Draw \(\Sigma\) ~ Inverse-Wishart(\(\nu^*\), \(\Lambda^*\))

2. Draw \(\mu\) | \(\Sigma\) ~ N_p(\(\mu^*\), \(\Sigma\) / \(\kappa^*\))

This produces samples from the joint distribution P(\(\mu\), \(\Sigma\) | Y, X, \(a_0\)) and captures both uncertainty in the mean and uncertainty in the covariance structure, including their dependence.

Marginal Sampling (marginal=TRUE):

Uses the marginal t-distribution of the mean:

$$\mu | Y, X, a_0 \sim t_{\nu^*-p+1}(\mu^*, \Lambda^* / (\kappa^*(\nu^*-p+1)))$$

This is computationally faster and useful when you primarily care about inference on the mean vector, marginalizing over uncertainty in the covariance.

Examples

library(MASS)
Sigma_true <- matrix(c(4, 1, 1, 2), 2, 2)
historical <- mvrnorm(50, mu = c(10, 5), Sigma = Sigma_true)
current <- mvrnorm(30, mu = c(10.5, 5.2), Sigma = Sigma_true)

pp <- powerprior_multivariate(historical, a0 = 0.5)
posterior <- posterior_multivariate(pp, current)

# Sample from joint distribution (covariance included)
samples_joint <- sample_posterior_multivariate(posterior, n_samples = 100)
cat("Mean samples shape:", dim(samples_joint$mu), "\n")
cat("Covariance samples shape:", dim(samples_joint$Sigma), "\n")

# Sample only means (faster)
samples_marginal <- sample_posterior_multivariate(posterior, n_samples = 100,
                                                   marginal = TRUE)