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 covariance and correlation matrices that are the dependence structure. Such a matrix cannot be arbitrary: it must be positive-definite, meaning every portfolio built from the variables has positive variance, which is the matrix version of “a variance cannot be negative” and is what makes this a matrix problem rather than a handful of separate scalar ones. 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 — the inverse of the covariance, which is the more convenient object because a zero in it means two variables are independent once the others are accounted for, and conjugate meaning the posterior stays in the same family as the prior, so the update is arithmetic rather than integration. 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  "Matrix-Variate & Covariance Distributions"
(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