| Title: | Bayesian Inference for the Wishart Distribution Parameters |
|---|---|
| Description: | Posterior inference for the shape parameter alpha and mean matrix mu in the model X_i ~ Wishart_p(2*alpha, Sigma), under both an improper prior and a proper Gamma/inverse-Wishart prior. The posterior mode is found via a Newton-within-EM algorithm and joint samples are drawn via rejection sampling. |
| Authors: | Phil Everson [aut], Hanqi Shi [aut, cre] |
| Maintainer: | Hanqi Shi <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2026-07-17 04:51:02 UTC |
| Source: | https://github.com/sunnyhq-shi/wishart-rejection-sampling |
Uses Armadillo's log_det for numerical stability.
ldet(X)ldet(X)
X |
Square positive definite matrix |
log|X|
The improper prior p(alpha,mu) propto (alpha-(p-1)/2)^(-1) *
|mu|^(-(p+1)/2), after marginalizing over mu, gives a closed-form
log posterior for alpha. Returns -Inf when a <= (p-1)/2.
lfafun_improper(a, n, p, xbar, ldetxbarg)lfafun_improper(a, n, p, xbar, ldetxbarg)
a |
Shape parameter alpha, must be > (p-1)/2 |
n |
Number of Wishart observations |
p |
Dimension |
xbar |
Sample mean matrix (p x p) |
ldetxbarg |
Log geometric mean of determinants (1/n)*sum log|X_i| |
log f*(alpha), or -Inf if alpha <= (p-1)/2
The proper prior alpha ~ Gamma(beta, beta*eta) (truncated at
(p-1)/2) and mu|alpha ~ inv-Wishart(2*kappa*alpha, 2*kappa*alpha*mu0),
after marginalizing over mu, gives a closed-form log posterior for
alpha. Returns -Inf outside the domain alpha > (p-1)/2.
lfafun_proper(a, n, p, ldet_muhat, ldetxbarg, ldet_mu0, beta, eta, kappa)lfafun_proper(a, n, p, ldet_muhat, ldetxbarg, ldet_mu0, beta, eta, kappa)
a |
Shape parameter alpha, must be > (p-1)/2 |
n |
Number of Wishart observations |
p |
Dimension |
ldet_muhat |
log|muhat| where muhat = (n*xbar + kappa*mu0)/(n+kappa) |
ldetxbarg |
Log geometric mean of determinants (1/n)*sum log|X_i| |
ldet_mu0 |
log|mu0| |
beta |
Gamma prior shape, must be > 1 |
eta |
Gamma prior rate parameter |
kappa |
inv-Wishart prior strength, must be >= 1 |
log f*(alpha), or -Inf if outside domain
Defined as log Gamma_p(a) = p*(p-1)/4 * log(pi) plus the sum over i = 1, ..., p of log Gamma(a - (i-1)/2). Required for the normalizing constant of the Wishart and inverse-Wishart distributions.
lgammap_export(a, p)lgammap_export(a, p)
a |
Argument, must be > (p-1)/2 |
p |
Dimension |
log Gamma_p(a)
Dispatches to the improper- or proper-prior EM mode finder depending
on which parameters are supplied. The prior type is auto-detected:
omit mu0, beta, eta, kappa for the
improper prior; supply mu0 and beta > 1, kappa >= 1
for the proper prior. Uses a Newton-within-EM algorithm.
mode_alphaEM( n, p, xbar, ldetxbarg, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L, max_nr_iter = 100L )mode_alphaEM( n, p, xbar, ldetxbarg, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L, max_nr_iter = 100L )
n |
Number of Wishart observations |
p |
Dimension |
xbar |
Sample mean matrix (p x p) |
ldetxbarg |
Log geometric mean of determinants |
mu0_ |
Prior center matrix (p x p); omit for improper prior |
beta |
Gamma prior shape, must be > 1 (ignored if improper) |
eta |
Gamma prior rate (ignored if improper) |
kappa |
inv-Wishart prior strength, must be >= 1 (ignored if improper) |
tol |
Convergence tolerance |
prnt |
If TRUE, print EM iterates and bounds |
max_em_iter |
Maximum number of EM iterations |
max_nr_iter |
Maximum number of inner Newton-Raphson iterations |
A numeric vector c(ahat, log f*(ahat))
Requires the posterior mode ahat and mxlfa = log f*(ahat)
from mode_alphaEM(), plus the covering Gamma parameters
lambda and nu_star. Prior type is auto-detected from
parameters, exactly as in mode_alphaEM().
rejection_sampler( ahat, mxlfa, lambda, nu_star, p, n, xbar, ldetxbarg, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, nsamp = 10000L )rejection_sampler( ahat, mxlfa, lambda, nu_star, p, n, xbar, ldetxbarg, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, nsamp = 10000L )
ahat |
Posterior mode from |
mxlfa |
log f*(ahat) from |
lambda |
Covering Gamma rate |
nu_star |
Covering Gamma shape (= ahat*lambda + 1) |
p |
Dimension |
n |
Number of Wishart observations |
xbar |
Sample mean matrix (p x p) |
ldetxbarg |
Log geometric mean of determinants |
mu0_ |
Prior center matrix; omit for improper prior |
beta |
Gamma prior shape; 0 = improper, must be > 1 if proper |
eta |
Gamma prior rate; default 1.0 |
kappa |
inv-Wishart prior strength; 0 = improper, must be >= 1 if proper |
nsamp |
Number of posterior samples (default 10000) |
A list with alpha_sample, mu_sample,
empirical_acpt_rate, theoretical_acpt_rate
Uses the Bartlett decomposition: if L is the lower Cholesky factor of Sigma, then A = L * T satisfies W = A * A^T ~ Wishart_p(nu, Sigma), where T is lower triangular with T(i,i) ~ sqrt(Chi^2(nu - i)) and T(i,j) ~ N(0,1) for i > j.
rwishart(n, p, nu, Sigma)rwishart(n, p, nu, Sigma)
n |
Number of draws |
p |
Dimension of the matrix, must be >= 2 |
nu |
Degrees of freedom, must be > p - 1 |
Sigma |
Scale matrix (p x p), must be positive definite |
A cube of size (p, p, n) containing the Wishart draws
Hello-world test to verify the library loaded correctly
test()test()
Always returns 0. Also prints "Hello from C++" as a side effect.
Given X_1,...,X_n ~ iid Wishart_p(2*alpha, Sigma), computes the
posterior of (alpha, mu) where mu = alpha*Sigma, under either an
improper or proper prior. Runs the full pipeline: sufficient
statistics, EM mode finding, and rejection sampling.
wishart_inference( X, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, nsamp = 10000L, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L, max_nr_iter = 100L )wishart_inference( X, mu0_ = NULL, beta = 0, eta = 1, kappa = 0, nsamp = 10000L, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L, max_nr_iter = 100L )
X |
Cube of n Wishart observations, dimensions (p, p, n) |
mu0_ |
Prior center matrix (p x p); omit for improper prior |
beta |
Gamma prior shape, must be > 1 (ignored if improper) |
eta |
Gamma prior rate parameter (ignored if improper) |
kappa |
inv-Wishart prior strength, must be >= 1 (ignored if improper) |
nsamp |
Number of posterior samples |
tol |
Convergence tolerance for the EM mode finder |
prnt |
If TRUE, print EM iterates |
max_em_iter |
Maximum number of EM iterations |
max_nr_iter |
Maximum number of inner Newton-Raphson iterations |
A list with results (alpha_samples, mu_samples, ahat,
theoretical_acpt_rate, empirical_acpt_rate) and statistics
(xbar, muhat, log_det_geometric_mean, cover_shape, cover_rate,
elapsed_seconds). If a convergence failure is caught internally,
returns list(error = "...") instead.
Returns the sample mean matrix xbar and the log geometric mean of determinants ldetxbarg, the sufficient statistics for (alpha, Sigma) under X_i ~ Wishart_p(2*alpha, Sigma).
wishart_stats(X)wishart_stats(X)
X |
Cube of n Wishart draws, dimensions (p, p, n) |
A list with xbar (p x p sample mean matrix) and
ldetxbarg (log geometric mean of determinants)