Matrix-Variate & Covariance Distributions

Python · R · Download distributions module

The Catalog

Folder 7 of the Statistical Distributions catalog. Folder 6 treated random vectors; this one treats random matrices — specifically the positive-definite covariance and correlation matrices that are the dependence structure. These are the laws you put a prior on when the covariance itself is unknown: the Wishart and Inverse-Wishart (conjugate priors for a Gaussian's precision and covariance), the Normal-Inverse-Wishart (their joint), and the LKJ (the modern prior on correlation matrices). Each is built from scratch — density, sampling algorithm, mean — and validated against scipy.stats.

Distribution	Prior on	Sampled by
Wishart	a precision matrix — the scatter-matrix law	Bartlett decomposition (vectorised)
Inverse-Wishart	a covariance matrix — the matrix Inverse-Gamma	inverting a Wishart draw
Normal-Inverse-Wishart	the pair $(\boldsymbol\mu,\Sigma)$ jointly	Inverse-Wishart, then $\boldsymbol\mu$ given $\Sigma$
LKJ( $\eta$ )	a correlation matrix	C-vine construction

Making a distribution over matrices visible

A distribution over matrices is hard to picture, so the notebooks make it geometric. A sampled covariance matrix is drawn as its confidence ellipse, which turns a distribution over covariances into a cloud of ellipses — and makes the concentration visible, since the spread shrinks as $1/\nu$ . A correlation matrix becomes a heatmap. And the space of $2\times2$ covariances is scattered directly into its positive-definite cone, with every draw confirmed to land inside the boundary $|W_{12}|<\sqrt{W_{11}W_{22}}$ — the constraint that makes this a matrix problem rather than three scalar ones.

Wishart, Inverse-Wishart and their joint

The Wishart is the distribution of $\sum_i \mathbf{z}_i\mathbf{z}_i^\top$ for Gaussian $\mathbf{z}_i$ — the sampling law of an unnormalised sample covariance, and hence the conjugate prior for a precision matrix. It is sampled here by the Bartlett decomposition, which never forms a single outer product: it builds a lower-triangular factor with $\sqrt{\chi^2}$ on the diagonal and standard normals below, then multiplies through the Cholesky factor of $V$ . Inverting a Wishart draw gives the Inverse-Wishart, the classical conjugate prior for a covariance and the matrix generalisation of the Inverse-Gamma; pairing that with a Gaussian for the mean gives the Normal-Inverse-Wishart, whose marginal for $\boldsymbol\mu$ — once $\Sigma$ is integrated out — is a heavy-tailed multivariate $t$ rather than a Gaussian. That extra weight is the estimation-risk premium for not knowing the covariance.

W\sim\mathcal{W}_p(\nu,V):\quad \mathbb{E}[W]=\nu V,\qquad \operatorname{Var}(W_{ii})=2\nu V_{ii}^{2},\qquad \mathbb{E}[\Sigma]=\frac{\Psi}{\nu-p-1}\ \ (\text{inverse})

The LKJ prior and the separation strategy

The LKJ distribution puts a prior directly on correlation matrices, with a single shape $\eta$ controlling how strongly it favours the identity: $\eta=1$ is uniform over correlation matrices, $\eta>1$ pulls the off-diagonals toward zero, and $\eta<1$ pushes them toward $\pm1$ . It is sampled by the C-vine construction and validated against the exact off-diagonal marginal, which for $p=2$ is available in closed form. The folder closes on the separation strategy that is now the recommended way to build a covariance prior — $\Sigma=\operatorname{diag}(\boldsymbol\tau)\,R\,\operatorname{diag}(\boldsymbol\tau)$ with $R\sim\mathrm{LKJ}(\eta)$ and independent scale priors on $\boldsymbol\tau$ — drawn as ellipse clouds that show the two knobs acting separately: $\eta$ controlling orientation and elongation, the Half-Cauchy scales letting a few ellipses grow large.

Notebooks

The Python notebook checks the Wishart density against scipy at $1.2\times10^{-18}$ and the Inverse-Wishart at $5.3\times10^{-14}$ , and confirms the moment structure directly — sample mean against $\nu V$ , and $\operatorname{Var}(W_{ii})$ against $2\nu V_{ii}^2$ . The R notebook rebuilds everything in base R and adds a check the Python side cannot: at $p=1$ the Wishart is a scaled chi-square, verified against dchisq to $9\times10^{-17}$ , plus an independent cross-check of the sample mean against base R's compiled rWishart. The timing section carries the folder's one genuinely surprising result: because the Bartlett decomposition vectorises across draws — one batched triangular construction instead of one per sample — the from-scratch Wishart and Inverse-Wishart run roughly six times faster than scipy's per-draw .rvs. The LKJ C-vine and the NIW cannot vectorise the same way, since each draw needs its own Cholesky, and they run as per-draw loops accordingly.

Python Four matrix-variate laws from scratch — Bartlett-sampled Wishart, Inverse-Wishart, Normal-Inverse-Wishart and the LKJ C-vine, visualised as ellipse clouds, the PD cone and correlation heatmaps (Python) Open → R Independent R re-implementation — base R from scratch, with the p=1 Wishart verified against dchisq and an independent cross-check against compiled rWishart (R) Open →

Downloads

matrixvar.py All four laws — densities, means, Bartlett Wishart sampling, NIW and the LKJ C-vine with its exact off-diagonal marginal, NumPy + scipy.special only

Distributions Module — Source Code

"""
matrixvar.py -- From-scratch MATRIX-VARIATE & COVARIANCE distributions.

Backs the notebooks in  Z-Statistical Distributions-MatrixVariate-Covariance
(Folder 7 of the Statistical Distributions catalog).

These are distributions over MATRICES -- specifically over covariance and
correlation matrices -- the priors and sampling laws for the dependence structure
itself. Everything is implemented from the definitions using only linear algebra
and the (multivariate) log-gamma as a primitive; `scipy.stats` (Python) and R
(R notebook) are the INDEPENDENT references.

Distributions:
    Wishart(nu, V)            -- over p x p positive-definite matrices (scatter / precision)
    Inverse-Wishart(nu, Psi)  -- over covariance matrices (conjugate prior for a MVN Sigma)
    Normal-Inverse-Wishart    -- joint conjugate prior for (mu, Sigma)
    LKJ(eta)                  -- over correlation matrices

Highlights (the from-scratch sampling algorithms):
    Wishart           : Bartlett decomposition  W = (L A)(L A)^T, A lower-triangular
                        with chi diagonal and N(0,1) below,  V = L L^T
    Inverse-Wishart   : invert a Wishart draw
    LKJ               : the C-vine partial-correlation method

Parameterisations match scipy where a match exists:
    wishart(nu, V)        scipy.stats.wishart(df=nu, scale=V)
    invwishart(nu, Psi)   scipy.stats.invwishart(df=nu, scale=Psi)
"""

import numpy as np
from scipy.special import gammaln


# --------------------------------------------------------------------------- #
#  Generic helpers                                                              #
# --------------------------------------------------------------------------- #

def _mvgammaln(a, p):
    """log of the multivariate gamma  Gamma_p(a)."""
    j = np.arange(1, p + 1)
    return p * (p - 1) / 4 * np.log(np.pi) + np.sum(gammaln(a + (1 - j) / 2))

def cov_to_ellipse(S, nsd=1.0, npts=100):
    """Return the (x,y) points of the nsd-standard-deviation ellipse of a 2x2 cov S."""
    vals, vecs = np.linalg.eigh(S)
    t = np.linspace(0, 2 * np.pi, npts)
    circ = np.array([np.cos(t), np.sin(t)])
    return (vecs @ np.diag(nsd * np.sqrt(vals)) @ circ)


# --------------------------------------------------------------------------- #
#  Wishart(nu, V)                                                               #
# --------------------------------------------------------------------------- #

def wishart_logpdf(W, nu, V):
    W = np.asarray(W, float); V = np.asarray(V, float); p = V.shape[0]
    _, ldW = np.linalg.slogdet(W); _, ldV = np.linalg.slogdet(V)
    tr = np.trace(np.linalg.solve(V, W))
    return ((nu - p - 1) / 2 * ldW - tr / 2 - nu * p / 2 * np.log(2)
            - nu / 2 * ldV - _mvgammaln(nu / 2, p))

def wishart_pdf(W, nu, V):
    return np.exp(wishart_logpdf(W, nu, V))

def wishart_rng(nu, V, n, rng):
    """Bartlett decomposition: with V = L L^T and a lower-triangular A whose diagonal
    is A_ii = sqrt(chi^2_{nu-i}) and whose below-diagonal entries are N(0,1),
    W = (L A)(L A)^T ~ Wishart_p(nu, V). Vectorised over n draws."""
    V = np.asarray(V, float); p = V.shape[0]
    L = np.linalg.cholesky(V)
    A = np.zeros((n, p, p))
    for i in range(p):
        A[:, i, i] = np.sqrt(2.0 * _gamma_shape((nu - i) / 2.0, n, rng))   # chi^2_{nu-i}
        for j in range(i):
            A[:, i, j] = _std_normal(n, rng)
    M = np.einsum("ab,nbc->nac", L, A)
    return np.einsum("nac,nbc->nab", M, M)                                  # M M^T

def wishart_mean(nu, V):
    return nu * np.asarray(V, float)


# --------------------------------------------------------------------------- #
#  Inverse-Wishart(nu, Psi)                                                     #
# --------------------------------------------------------------------------- #

def invwishart_logpdf(S, nu, Psi):
    S = np.asarray(S, float); Psi = np.asarray(Psi, float); p = Psi.shape[0]
    _, ldS = np.linalg.slogdet(S); _, ldP = np.linalg.slogdet(Psi)
    tr = np.trace(np.linalg.solve(S, Psi))
    return (nu / 2 * ldP - nu * p / 2 * np.log(2) - _mvgammaln(nu / 2, p)
            - (nu + p + 1) / 2 * ldS - tr / 2)

def invwishart_pdf(S, nu, Psi):
    return np.exp(invwishart_logpdf(S, nu, Psi))

def invwishart_rng(nu, Psi, n, rng):
    """If W ~ Wishart(nu, Psi^{-1}) then W^{-1} ~ Inverse-Wishart(nu, Psi)."""
    Psi = np.asarray(Psi, float)
    W = wishart_rng(nu, np.linalg.inv(Psi), n, rng)
    return np.linalg.inv(W)

def invwishart_mean(nu, Psi):
    Psi = np.asarray(Psi, float); p = Psi.shape[0]
    return Psi / (nu - p - 1) if nu > p + 1 else np.full_like(Psi, np.nan)


# --------------------------------------------------------------------------- #
#  Normal-Inverse-Wishart(mu0, kappa, Psi, nu)  --  joint prior for (mu, Sigma)  #
# --------------------------------------------------------------------------- #

def niw_rng(mu0, kappa, Psi, nu, n, rng):
    """Sigma ~ Inverse-Wishart(nu, Psi);  mu | Sigma ~ Normal(mu0, Sigma/kappa).
    Returns (mus (n,p), Sigmas (n,p,p))."""
    mu0 = np.asarray(mu0, float); p = len(mu0)
    Sig = invwishart_rng(nu, Psi, n, rng)
    mus = np.empty((n, p))
    for r in range(n):
        L = np.linalg.cholesky(Sig[r] / kappa)
        mus[r] = mu0 + L @ _std_normal(p, rng)
    return mus, Sig


# --------------------------------------------------------------------------- #
#  LKJ(eta)  --  distribution over correlation matrices                          #
# --------------------------------------------------------------------------- #

def lkj_rng(eta, p, n, rng):
    """C-vine method: draw partial correlations from shifted Betas and convert to a
    correlation matrix. eta=1 is uniform over correlation matrices; eta>1 concentrates
    near the identity; eta<1 favours strong correlations. Returns (n, p, p)."""
    out = np.empty((n, p, p))
    for r in range(n):
        out[r] = _lkj_one(eta, p, rng)
    return out

def _lkj_one(eta, p, rng):
    P = np.zeros((p, p)); R = np.eye(p)
    beta = eta + (p - 1) / 2.0
    for k in range(p - 1):
        beta -= 0.5
        for i in range(k + 1, p):
            P[k, i] = 2.0 * _beta_rng_scalar(beta, beta, rng) - 1.0
            pcor = P[k, i]
            for l in range(k - 1, -1, -1):
                pcor = pcor * np.sqrt((1 - P[l, i] ** 2) * (1 - P[l, k] ** 2)) + P[l, i] * P[l, k]
            R[k, i] = pcor; R[i, k] = pcor
    return R

def lkj_logpdf_unnormalized(R, eta):
    """log density up to the (eta-independent-in-shape) constant:  (eta-1) log det R."""
    _, ld = np.linalg.slogdet(np.asarray(R, float))
    return (eta - 1) * ld

def lkj_marginal_rho_pdf(rho, eta, p):
    """The marginal density of any single off-diagonal correlation of an LKJ(eta)
    correlation matrix in dimension p:  (rho+1)/2 ~ Beta(a, a) with a = eta + (p-2)/2,
    so f(rho) proportional to (1-rho^2)^{a-1}."""
    a = eta + (p - 2) / 2.0
    rho = np.asarray(rho, float)
    logB = gammaln(a) + gammaln(a) - gammaln(2 * a)
    # density of rho on (-1,1): 2^{1-2a}/B(a,a) * (1-rho^2)^{a-1}
    return np.where(np.abs(rho) < 1,
                    (1 - rho ** 2) ** (a - 1) * 2.0 ** (1 - 2 * a) / np.exp(logB), 0.0)


# --------------------------------------------------------------------------- #
#  Private RNG primitives                                                        #
# --------------------------------------------------------------------------- #

def _std_normal(n, rng):
    m = (n + 1) // 2
    u1 = rng.random(m); u2 = rng.random(m)
    r = np.sqrt(-2 * np.log(u1)); z = np.empty(2 * m)
    z[:m] = r * np.cos(2 * np.pi * u2); z[m:] = r * np.sin(2 * np.pi * u2)
    return z[:n]

def _gamma_shape(k, n, rng):
    n = int(n)
    if k < 1:
        g = _gamma_shape(k + 1.0, n, rng); u = rng.random(n)
        return g * u ** (1.0 / k)
    d = k - 1.0 / 3.0; c = 1.0 / np.sqrt(9.0 * d)
    res = np.empty(n); todo = np.ones(n, bool)
    while todo.any():
        m = int(todo.sum()); x = _std_normal(m, rng)
        v = (1.0 + c * x) ** 3; u = rng.random(m)
        ok = (v > 0) & (np.log(u) < 0.5 * x ** 2 + d - d * v + d * np.log(np.where(v > 0, v, 1.0)))
        cur = np.where(todo)[0]; acc = cur[ok]
        res[acc] = d * v[ok]; todo[acc] = False
    return res

def _beta_rng_scalar(a, b, rng):
    x = _gamma_shape(a, 1, rng)[0]; y = _gamma_shape(b, 1, rng)[0]
    return x / (x + y)

References

Muirhead, R. J. (2005). Aspects of Multivariate Statistical Theory. Wiley. — the Wishart and Inverse-Wishart densities, moments and distributional theory
Bartlett, M. S. (1933). On the theory of statistical regression. Proceedings of the Royal Society of Edinburgh 53, 260–283. — the triangular decomposition used to sample the Wishart without summing outer products
Lewandowski, D., Kurowicka, D. & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis 100(9), 1989–2001. — the LKJ distribution and the C-vine sampler implemented here
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A. & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman & Hall/CRC. — the Normal-Inverse-Wishart conjugate analysis and its multivariate-t marginal for the mean
Barnard, J., McCulloch, R. & Meng, X.-L. (2000). Modeling covariance matrices in terms of standard deviations and correlations. Statistica Sinica 10(4), 1281–1311. — the separation strategy: scales and correlations given independent priors
Stan Development Team (2024). Stan User's Guide — Regression Models: priors on covariance. — the LKJ-plus-scales recommendation this folder reproduces