| 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: | Philip Everson [aut], Hanqi Shi [aut, cre] |
| Maintainer: | Hanqi Shi <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2026-10-04 10:03:51 UTC |
| Source: | https://github.com/cran/wishartinference |
Given x_1,...,x_n ~ iid Gamma(alpha, alpha/mu) – the p = 1
reduction of X_i ~ Wishart_p(2*alpha, Sigma) – computes the
posterior of (alpha, mu) under either an improper or proper prior.
Scalar (p = 1) analogue of wishart_inference(): takes a plain
numeric vector instead of a p x p x n array, since each observation
is already a positive scalar. There is no unifying dispatch with
wishart_inference() – call this directly for scalar data.
gamma_inference( x, mu0 = 1, beta = 0, eta = 1, kappa = 0, nsamp = 10000L, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L )gamma_inference( x, mu0 = 1, beta = 0, eta = 1, kappa = 0, nsamp = 10000L, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L )
x |
Numeric vector of n scalar observations, each > 0 |
mu0 |
Prior center for mu, must be > 0 (ignored if improper) |
beta |
Gamma prior shape, must be >= 0 (ignored if improper) |
eta |
Gamma prior rate, must be >= 0 (ignored if improper) |
kappa |
Prior strength on mu, must be >= 0 (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 |
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.
set.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3, mu = 1 # improper prior (the default) res <- gamma_inference(x, nsamp = 200) res$results$ahat quantile(res$results$alpha_samples, c(0.025, 0.975)) # proper prior res2 <- gamma_inference(x, mu0 = 1, beta = 5, eta = 0.5, kappa = 2, nsamp = 200) res2$results$ahatset.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3, mu = 1 # improper prior (the default) res <- gamma_inference(x, nsamp = 200) res$results$ahat quantile(res$results$alpha_samples, c(0.025, 0.975)) # proper prior res2 <- gamma_inference(x, mu0 = 1, beta = 5, eta = 0.5, kappa = 2, nsamp = 200) res2$results$ahat
Uses Armadillo's log_det for numerical stability.
ldet(X)ldet(X)
X |
Square positive definite matrix |
log|X|
ldet(diag(c(1, 2, 3))) # log(6) set.seed(1) A <- crossprod(matrix(rnorm(9), 3, 3)) + diag(3) ldet(A) log(det(A))ldet(diag(c(1, 2, 3))) # log(6) set.seed(1) A <- crossprod(matrix(rnorm(9), 3, 3)) + diag(3) ldet(A) log(det(A))
The p = 1 reduction of the Wishart model: X_i ~ Gamma(alpha, alpha/mu). Auto-detects improper vs. proper exactly as the p >= 2 functions do (beta = 0, eta = 1, kappa = 0 = improper), except kappa = 0 alone (with beta > 0) is also valid here – a proper prior on alpha with a flat/improper treatment of mu – unlike the p >= 2 backend, which requires kappa >= 1 whenever the prior isn't fully improper.
lfafun_gamma(a, n, xbar, ldetxbarg, beta = 0, eta = 1, kappa = 0, mu0 = 1)lfafun_gamma(a, n, xbar, ldetxbarg, beta = 0, eta = 1, kappa = 0, mu0 = 1)
a |
Shape parameter alpha, must be > 0 |
n |
Number of scalar observations |
xbar |
Sample mean of the observations |
ldetxbarg |
mean(log(x)) – scalar analogue of ldetxbarg |
beta |
Gamma prior shape, must be >= 0 (0 = improper) |
eta |
Gamma prior rate, must be >= 0 |
kappa |
Prior strength on mu, must be >= 0 (0 = improper on mu) |
mu0 |
Prior center for mu, must be > 0 |
log f*(alpha), or -Inf if a <= 0
set.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3, mu = 1 lfafun_gamma(3, 20, mean(x), mean(log(x))) # proper prior on alpha lfafun_gamma(3, 20, mean(x), mean(log(x)), beta = 5, eta = 0.5, kappa = 2, mu0 = 1)set.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3, mu = 1 lfafun_gamma(3, 20, mean(x), mean(log(x))) # proper prior on alpha lfafun_gamma(3, 20, mean(x), mean(log(x)), beta = 5, eta = 0.5, kappa = 2, mu0 = 1)
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
set.seed(1) X <- rwishart(10, 2, 6, diag(2)) st <- wishart_stats(X) lfafun_improper(3, 10, 2, st$xbar, st$ldetxbarg) # -Inf outside the domain alpha > (p-1)/2 lfafun_improper(0.2, 10, 2, st$xbar, st$ldetxbarg)set.seed(1) X <- rwishart(10, 2, 6, diag(2)) st <- wishart_stats(X) lfafun_improper(3, 10, 2, st$xbar, st$ldetxbarg) # -Inf outside the domain alpha > (p-1)/2 lfafun_improper(0.2, 10, 2, st$xbar, st$ldetxbarg)
The proper prior (alpha - (p-1)/2) ~ Gamma(beta, rate = beta*eta) 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 >= 0 |
eta |
Gamma prior rate parameter, must be >= 0 |
kappa |
inv-Wishart prior strength, must be >= 1 |
log f*(alpha), or -Inf if outside domain
set.seed(1) n <- 10; p <- 2; kappa <- 1 X <- rwishart(n, p, 6, diag(p)) st <- wishart_stats(X) mu0 <- diag(2, p) muhat <- (n * st$xbar + kappa * mu0) / (n + kappa) lfafun_proper(3, n, p, ldet(muhat), st$ldetxbarg, ldet(mu0), beta = 5, eta = 0.5, kappa = kappa)set.seed(1) n <- 10; p <- 2; kappa <- 1 X <- rwishart(n, p, 6, diag(p)) st <- wishart_stats(X) mu0 <- diag(2, p) muhat <- (n * st$xbar + kappa * mu0) / (n + kappa) lfafun_proper(3, n, p, ldet(muhat), st$ldetxbarg, ldet(mu0), beta = 5, eta = 0.5, kappa = kappa)
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)
lgammap_export(5, 3) # equals the defining sum a <- 5; p <- 3 p * (p - 1) / 4 * log(pi) + sum(lgamma(a - (0:(p - 1)) / 2))lgammap_export(5, 3) # equals the defining sum a <- 5; p <- 3 p * (p - 1) / 4 * log(pi) + sum(lgamma(a - (0:(p - 1)) / 2))
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 >= 0, 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 >= 0 (ignored if improper) |
eta |
Gamma prior rate, must be >= 0 (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))
set.seed(1) n <- 10; p <- 2 X <- rwishart(n, p, 6, diag(p)) # true alpha = 3 st <- wishart_stats(X) # improper prior: omit mu0, beta, eta, kappa em <- mode_alphaEM(n, p, st$xbar, st$ldetxbarg) em[1] # posterior mode of alpha em[2] # log f*(ahat) # proper prior: (alpha - (p-1)/2) ~ Gamma(beta, rate = beta*eta) mode_alphaEM(n, p, st$xbar, st$ldetxbarg, mu0 = diag(2, p), beta = 5, eta = 0.5, kappa = 1)set.seed(1) n <- 10; p <- 2 X <- rwishart(n, p, 6, diag(p)) # true alpha = 3 st <- wishart_stats(X) # improper prior: omit mu0, beta, eta, kappa em <- mode_alphaEM(n, p, st$xbar, st$ldetxbarg) em[1] # posterior mode of alpha em[2] # log f*(ahat) # proper prior: (alpha - (p-1)/2) ~ Gamma(beta, rate = beta*eta) mode_alphaEM(n, p, st$xbar, st$ldetxbarg, mu0 = diag(2, p), beta = 5, eta = 0.5, kappa = 1)
Scalar (p = 1) analogue of mode_alphaEM(). See the C++ doc
comment for details on the bisection-style ascent used here and the
n = 2, kappa = 0, beta = 0 structural boundary case.
mode_alphaEM_gamma( n, xbar, ldetxbarg, beta = 0, eta = 1, kappa = 0, mu0 = 1, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L )mode_alphaEM_gamma( n, xbar, ldetxbarg, beta = 0, eta = 1, kappa = 0, mu0 = 1, tol = 1e-06, prnt = FALSE, max_em_iter = 1000L )
n |
Number of scalar observations |
xbar |
Sample mean of the observations |
ldetxbarg |
mean(log(x)) |
beta |
Gamma prior shape, must be >= 0 (0 = improper) |
eta |
Gamma prior rate, must be >= 0 |
kappa |
Prior strength on mu, must be >= 0 (0 = improper on mu) |
mu0 |
Prior center for mu, must be > 0 |
tol |
Convergence tolerance |
prnt |
If TRUE, print EM iterates and bounds |
max_em_iter |
Maximum number of EM iterations |
A numeric vector c(ahat, log f*(ahat))
set.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3 em <- mode_alphaEM_gamma(20, mean(x), mean(log(x))) em[1] # posterior mode of alpha # proper prior mode_alphaEM_gamma(20, mean(x), mean(log(x)), beta = 5, eta = 0.5, kappa = 2, mu0 = 1)set.seed(1) x <- rgamma(20, shape = 3, rate = 3) # true alpha = 3 em <- mode_alphaEM_gamma(20, mean(x), mean(log(x))) em[1] # posterior mode of alpha # proper prior mode_alphaEM_gamma(20, mean(x), mean(log(x)), beta = 5, eta = 0.5, kappa = 2, mu0 = 1)
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 (with defaults below), must be >= 0 if proper |
eta |
Gamma prior rate; default 1.0, must be >= 0 if proper |
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
set.seed(1) n <- 10; p <- 2 X <- rwishart(n, p, 6, diag(p)) st <- wishart_stats(X) em <- mode_alphaEM(n, p, st$xbar, st$ldetxbarg) # covering Gamma parameters under the improper prior lambda <- max(n * (ldet(st$xbar) - st$ldetxbarg), 1e-10) nu_star <- em[1] * lambda + 1 samp <- rejection_sampler(em[1], em[2], lambda, nu_star, p, n, st$xbar, st$ldetxbarg, nsamp = 100) quantile(samp$alpha_sample, c(0.025, 0.5, 0.975)) samp$empirical_acpt_rateset.seed(1) n <- 10; p <- 2 X <- rwishart(n, p, 6, diag(p)) st <- wishart_stats(X) em <- mode_alphaEM(n, p, st$xbar, st$ldetxbarg) # covering Gamma parameters under the improper prior lambda <- max(n * (ldet(st$xbar) - st$ldetxbarg), 1e-10) nu_star <- em[1] * lambda + 1 samp <- rejection_sampler(em[1], em[2], lambda, nu_star, p, n, st$xbar, st$ldetxbarg, nsamp = 100) quantile(samp$alpha_sample, c(0.025, 0.5, 0.975)) samp$empirical_acpt_rate
Scalar (p = 1) analogue of rejection_sampler(). mu | alpha, x is
drawn from reciprocal-Gamma(nk_eff*alpha, alpha*(n*xbar+kappa*mu0)) –
see the C++ doc comment for the derivation.
rejection_sampler_gamma( n, xbar, ldetxbarg, ahat, mxlfa, beta = 0, eta = 1, kappa = 0, mu0 = 1, nsamp = 10000L )rejection_sampler_gamma( n, xbar, ldetxbarg, ahat, mxlfa, beta = 0, eta = 1, kappa = 0, mu0 = 1, nsamp = 10000L )
n |
Number of scalar observations |
xbar |
Sample mean of the observations |
ldetxbarg |
mean(log(x)) |
ahat |
Posterior mode from |
mxlfa |
log f*(ahat) from |
beta |
Gamma prior shape, must be >= 0 (0 = improper) |
eta |
Gamma prior rate, must be >= 0 |
kappa |
Prior strength on mu, must be >= 0 (0 = improper on mu) |
mu0 |
Prior center for mu, must be > 0 |
nsamp |
Number of posterior samples (default 10000) |
A list with alpha_sample, mu_sample,
empirical_acpt_rate, theoretical_acpt_rate
set.seed(1) x <- rgamma(20, shape = 3, rate = 3) em <- mode_alphaEM_gamma(20, mean(x), mean(log(x))) samp <- rejection_sampler_gamma(20, mean(x), mean(log(x)), em[1], em[2], nsamp = 100) quantile(samp$alpha_sample, c(0.025, 0.5, 0.975)) samp$empirical_acpt_rateset.seed(1) x <- rgamma(20, shape = 3, rate = 3) em <- mode_alphaEM_gamma(20, mean(x), mean(log(x))) samp <- rejection_sampler_gamma(20, mean(x), mean(log(x)), em[1], em[2], nsamp = 100) quantile(samp$alpha_sample, c(0.025, 0.5, 0.975)) samp$empirical_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
set.seed(1) X <- rwishart(5, 2, 6, diag(2)) dim(X) X[, , 1]set.seed(1) X <- rwishart(5, 2, 6, diag(2)) dim(X) X[, , 1]
A minimal smoke test, useful immediately after installation to confirm that the package's compiled code is available.
wishart_hello()wishart_hello()
Returns 0. Called for the confirmation message it prints.
wishart_hello()wishart_hello()
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 >= 0 (ignored if improper) |
eta |
Gamma prior rate parameter, must be >= 0 (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.
set.seed(1) X <- rwishart(10, 2, 6, diag(2)) # true alpha = 3 # improper prior (the default) res <- wishart_inference(X, nsamp = 200) res$results$ahat quantile(res$results$alpha_samples, c(0.025, 0.975)) # proper prior res2 <- wishart_inference(X, mu0 = diag(2, 2), beta = 5, eta = 0.5, kappa = 1, nsamp = 200) res2$results$ahatset.seed(1) X <- rwishart(10, 2, 6, diag(2)) # true alpha = 3 # improper prior (the default) res <- wishart_inference(X, nsamp = 200) res$results$ahat quantile(res$results$alpha_samples, c(0.025, 0.975)) # proper prior res2 <- wishart_inference(X, mu0 = diag(2, 2), beta = 5, eta = 0.5, kappa = 1, nsamp = 200) res2$results$ahat
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)
set.seed(1) X <- rwishart(10, 2, 6, diag(2)) st <- wishart_stats(X) st$xbar st$ldetxbarg # ldetxbarg is a log geometric mean, so it sits below log|xbar| st$ldetxbarg < ldet(st$xbar)set.seed(1) X <- rwishart(10, 2, 6, diag(2)) st <- wishart_stats(X) st$xbar st$ldetxbarg # ldetxbarg is a log geometric mean, so it sits below log|xbar| st$ldetxbarg < ldet(st$xbar)