Package {wishartinference}

wishartinference: Bayesian Inference for the Wishart Distribution Shape Parameter

Description

Posterior inference for the shape parameter alpha and mean matrix mu in the model X_i ~ Wishart_p(2*alpha, Sigma).

Author(s)

Maintainer: Hanqi Shi hshi1@swarthmore.edu

Authors:

• Hanqi Shi hshi1@swarthmore.edu

• Phil Everson peverso1@swarthmore.edu

Bayesian inference for the Wishart shape parameter, gamma (p = 1) case

Description

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.

Usage

gamma_inference(
  x,
  mu0 = 1,
  beta = 0,
  eta = 1,
  kappa = 0,
  nsamp = 10000L,
  tol = 1e-06,
  prnt = FALSE,
  max_em_iter = 1000L
)

Arguments

Value

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.

Examples

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$ahat

Log determinant of a positive definite matrix

Description

Uses Armadillo's log_det for numerical stability.

Usage

ldet(X)

Arguments

Value

log|X|

Examples

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))

Log unnormalized posterior of alpha, gamma (p = 1) case

Description

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.

Usage

lfafun_gamma(a, n, xbar, ldetxbarg, beta = 0, eta = 1, kappa = 0, mu0 = 1)

Arguments

Value

log f*(alpha), or -Inf if a <= 0

Examples

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)

Log unnormalized marginal posterior of alpha under the improper prior

Description

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.

Usage

lfafun_improper(a, n, p, xbar, ldetxbarg)

Arguments

Value

log f*(alpha), or -Inf if alpha <= (p-1)/2

Examples

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)

Log unnormalized marginal posterior of alpha under the proper prior

Description

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.

Usage

lfafun_proper(a, n, p, ldet_muhat, ldetxbarg, ldet_mu0, beta, eta, kappa)

Arguments

Value

log f*(alpha), or -Inf if outside domain

Examples

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)

Log of the multivariate Gamma function

Description

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.

Usage

lgammap_export(a, p)

Arguments

Value

log Gamma_p(a)

Examples

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))

Find the posterior mode of alpha

Description

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.

Usage

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
)

Arguments

Value

A numeric vector c(ahat, log f*(ahat))

Examples

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)

EM algorithm to find the posterior mode of alpha, gamma (p = 1) case

Description

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.

Usage

mode_alphaEM_gamma(
  n,
  xbar,
  ldetxbarg,
  beta = 0,
  eta = 1,
  kappa = 0,
  mu0 = 1,
  tol = 1e-06,
  prnt = FALSE,
  max_em_iter = 1000L
)

Arguments

Value

A numeric vector c(ahat, log f*(ahat))

Examples

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)

Rejection sampler for the joint posterior of (alpha, mu)

Description

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().

Usage

rejection_sampler(
  ahat,
  mxlfa,
  lambda,
  nu_star,
  p,
  n,
  xbar,
  ldetxbarg,
  mu0 = NULL,
  beta = 0,
  eta = 1,
  kappa = 0,
  nsamp = 10000L
)

Arguments

Value

A list with alpha_sample, mu_sample, empirical_acpt_rate, theoretical_acpt_rate

Examples

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_rate

Rejection sampler for the joint posterior of (alpha, mu), gamma (p = 1) case

Description

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.

Usage

rejection_sampler_gamma(
  n,
  xbar,
  ldetxbarg,
  ahat,
  mxlfa,
  beta = 0,
  eta = 1,
  kappa = 0,
  mu0 = 1,
  nsamp = 10000L
)

Arguments

Value

A list with alpha_sample, mu_sample, empirical_acpt_rate, theoretical_acpt_rate

Examples

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_rate

Generate draws from a Wishart distribution

Description

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.

Usage

rwishart(n, p, nu, Sigma)

Arguments

Value

A cube of size (p, p, n) containing the Wishart draws

Examples

set.seed(1)
X <- rwishart(5, 2, 6, diag(2))
dim(X)
X[, , 1]

Verify that the compiled library loaded correctly

Description

A minimal smoke test, useful immediately after installation to confirm that the package's compiled code is available.

Usage

wishart_hello()

Value

Returns 0. Called for the confirmation message it prints.

Examples

wishart_hello()

Bayesian inference for the Wishart shape parameter

Description

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.

Usage

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
)

Arguments

Value

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.

Examples

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$ahat

Compute sufficient statistics for Wishart observations

Description

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).

Usage

wishart_stats(X)

Arguments

Value

A list with xbar (p x p sample mean matrix) and ldetxbarg (log geometric mean of determinants)

Examples

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)