Matrix-Variate & Covariance Distributions¶
Statistical Distributions — priors on covariance and correlation matrices¶
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 distributions 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.
Because a distribution over matrices is hard to picture, we visualise it geometrically: a sampled covariance matrix is drawn as its confidence ellipse, so a distribution over covariance matrices becomes a cloud of ellipses; a correlation matrix is drawn as a heatmap; and the space of $2\times2$ covariances is a 3-D positive-definite cone we can scatter draws into.
The organising ideas¶
- Sampling a scatter matrix. 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. We sample it by the Bartlett decomposition, never summing outer products.
- A prior on covariance. Inverting a Wishart gives the Inverse-Wishart, the classical conjugate prior for a covariance $\Sigma$; pairing it with a Gaussian for the mean gives the Normal-Inverse-Wishart, the full conjugate prior for $(\boldsymbol\mu,\Sigma)$.
- A prior on correlation. The LKJ($\eta$) distribution puts a prior directly on correlation matrices, with a single shape $\eta$ controlling how strongly it favours the identity — the standard modern choice (Stan, PyMC), usually combined with separate scale priors.
Families covered¶
- Wishart — the scatter-matrix law (Bartlett sampling, ellipse clouds, the PD cone)
- Inverse-Wishart — the classical covariance prior
- Normal-Inverse-Wishart — the joint $(\boldsymbol\mu,\Sigma)$ prior
- LKJ — the correlation-matrix prior
- The modern covariance prior — LKJ + scale priors
- Timing
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa
import matrixvar as MV
from scipy import stats as S
plt.rcParams.update({"figure.figsize": (9,4), "axes.grid": True, "grid.alpha": .22,
"axes.spines.top": False, "axes.spines.right": False, "font.size": 10.5})
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#718096","#6b46c1"
rng = np.random.default_rng(0)
def ellipse_cloud(ax, mats, color, nsd=1.0, alpha=.10, lw=.7):
'''Draw each 2x2 covariance matrix as its nsd-sigma ellipse -- a distribution over
covariance matrices becomes a cloud of ellipses.'''
for Smat in mats:
E = MV.cov_to_ellipse(Smat, nsd)
ax.plot(E[0], E[1], color=color, alpha=alpha, lw=lw)
print("helpers ready")
helpers ready
1. Wishart — the scatter-matrix distribution¶
If $\mathbf z_1,\dots,\mathbf z_\nu\sim N_p(0,V)$, then $W=\sum_i\mathbf z_i\mathbf z_i^\top\sim\text{Wishart}_p(\nu,V)$: the distribution of a sum of outer products, i.e. an unnormalised sample covariance / scatter matrix. Its density on the cone of positive-definite matrices is $$f(W)\propto |W|^{(\nu-p-1)/2}\exp\!\big[-\tfrac12\text{tr}(V^{-1}W)\big],$$ its mean is $\nu V$, and each diagonal entry $W_{ii}\sim V_{ii}\,\chi^2_\nu$. As a prior it is conjugate for a Gaussian precision matrix $\Sigma^{-1}$. We sample it — without ever summing $\nu$ outer products — by the Bartlett decomposition: with $V=LL^\top$ and a lower-triangular $A$ whose diagonal is $\sqrt{\chi^2_{\nu-i}}$ and whose sub-diagonal is standard normal, $W=(LA)(LA)^\top$.
V = np.array([[2.0, 0.6],[0.6, 1.0]]); nu = 8
W = MV.wishart_rng(nu, V, 200000, rng)
print("pdf vs scipy (5 draws):", max(abs(MV.wishart_pdf(w, nu, V) - S.wishart(df=nu, scale=V).pdf(w)) for w in W[:5]))
print("sample mean\n", np.round(W.mean(0),3), "\nvs ν·V\n", np.round(nu*V,3))
print(f"W[0,0]: mean {W[:,0,0].mean():.2f} vs ν·V₀₀={nu*V[0,0]}, var {W[:,0,0].var():.1f} vs 2ν·V₀₀²={2*nu*V[0,0]**2}")
pdf vs scipy (5 draws): 1.1926223897340549e-18 sample mean [[16.003 4.807] [ 4.807 8. ]] vs ν·V [[16. 4.8] [ 4.8 8. ]] W[0,0]: mean 16.00 vs ν·V₀₀=16.0, var 64.1 vs 2ν·V₀₀²=64.0
# The Wishart as a distribution over covariance matrices: a cloud of ellipses
V = np.array([[2.0, 0.6],[0.6, 1.0]]); nu = 8
Wc = MV.wishart_rng(nu, V, 250, rng) / nu # divide by nu so ellipses centre on V
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ellipse_cloud(ax[0], Wc, BLUE, nsd=1.0, alpha=.10)
Emean = MV.cov_to_ellipse(V, 1.0); ax[0].plot(Emean[0], Emean[1], color=RED, lw=2.5, label="mean covariance V")
ax[0].set_title("Wishart(ν=8) — 250 sampled covariance ellipses (÷ν)"); ax[0].legend(fontsize=8); ax[0].set_aspect("equal"); ax[0].set_xlabel("x₁"); ax[0].set_ylabel("x₂")
# concentration: larger nu -> tighter cloud
for nu_,c,al in [(4,RED,.5),(20,GREEN,.7),(80,BLUE,.9)]:
Wc = MV.wishart_rng(nu_, V, 120, rng)/nu_
for Smat in Wc[:60]:
E = MV.cov_to_ellipse(Smat, 1.0); ax[1].plot(E[0], E[1], color=c, alpha=.10, lw=.6)
ax[1].plot([], [], color=c, label=f"ν={nu_}")
ax[1].set_title("Larger ν concentrates the Wishart around its mean"); ax[1].legend(fontsize=8); ax[1].set_aspect("equal"); ax[1].set_xlim(-3,3); ax[1].set_ylim(-2.5,2.5)
plt.tight_layout(); plt.show()
print("Each ellipse is one sampled covariance matrix; the red is the mean V. As ν grows the sampling spread shrinks (∝1/ν).")
Each ellipse is one sampled covariance matrix; the red is the mean V. As ν grows the sampling spread shrinks (∝1/ν).
# The space of 2x2 covariances is the positive-definite CONE: a 3-D scatter of draws
V = np.array([[1.5, 0.3],[0.3, 1.0]]); nu = 6
W = MV.wishart_rng(nu, V, 4000, rng)
fig = plt.figure(figsize=(12.5, 5))
ax = fig.add_subplot(1, 2, 1, projection="3d")
ax.scatter(W[:,0,0], W[:,1,1], W[:,0,1], s=3, alpha=.12, color=BLUE)
# the PD boundary: det = w11 w22 - w12^2 = 0 -> w12 = sqrt(w11 w22)
g = np.linspace(0.01, W[:,0,0].max(), 30); G1, G2 = np.meshgrid(g, g)
ax.plot_surface(G1, G2, np.sqrt(G1*G2), color=RED, alpha=.12, linewidth=0)
ax.plot_surface(G1, G2, -np.sqrt(G1*G2), color=RED, alpha=.12, linewidth=0)
ax.set_xlabel("W₁₁"); ax.set_ylabel("W₂₂"); ax.set_zlabel("W₁₂"); ax.set_title("Wishart draws in the positive-definite cone", fontsize=10); ax.view_init(22, -60)
ax2 = fig.add_subplot(1, 2, 2)
tt = np.linspace(0, W[:,0,0].max(), 300)
ax2.hist(W[:,0,0], bins=90, density=True, color="#cfe3f6", edgecolor="white", lw=.2, label="W₁₁ samples")
ax2.plot(tt, S.chi2.pdf(tt/V[0,0], nu)/V[0,0], color=RED, lw=2, label=f"V₁₁·χ²_{nu} density")
ax2.set_title("Diagonal marginal W₁₁ is a scaled χ²"); ax2.legend(fontsize=8); ax2.set_xlabel("W₁₁")
plt.tight_layout(); plt.show()
print("Every draw lands inside the PD cone |W₁₂|<√(W₁₁W₂₂) (red boundary); the diagonal marginal is a scaled χ².")
Every draw lands inside the PD cone |W₁₂|<√(W₁₁W₂₂) (red boundary); the diagonal marginal is a scaled χ².
2. Inverse-Wishart — the classical covariance prior¶
The Inverse-Wishart(ν, Ψ) is the distribution of $W^{-1}$ when $W\sim\text{Wishart}(\nu,\Psi^{-1})$ — a distribution over covariance matrices, and the classical conjugate prior for the covariance $\Sigma$ of a multivariate Normal. Its mean is $\Psi/(\nu-p-1)$ (for $\nu>p+1$), and larger $\nu$ concentrates it. It is the multivariate generalisation of the Inverse-Gamma variance prior from Folder 5. (Its one drawback — a single $\nu$ ties together the variability of all variances and correlations — is exactly what the LKJ construction in §5 fixes.)
Psi = np.array([[3.0, 0.6],[0.6, 2.0]]); nu = 8
Sg = MV.invwishart_rng(nu, Psi, 200000, rng)
print("pdf vs scipy (5 draws):", max(abs(MV.invwishart_pdf(s, nu, Psi) - S.invwishart(df=nu, scale=Psi).pdf(s)) for s in Sg[:5]))
print("sample mean\n", np.round(Sg.mean(0),3), "\nvs Ψ/(ν−p−1)\n", np.round(Psi/(nu-2-1),3))
# ellipse clouds at two concentrations
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
for k,nu_ in enumerate([6, 30]):
Sc = MV.invwishart_rng(nu_, Psi*(nu_-3), 150, rng) # scale Psi so the mean stays ~fixed
for Smat in Sc:
E = MV.cov_to_ellipse(Smat, 1.0); ax[k].plot(E[0], E[1], color=BLUE, alpha=.12, lw=.7)
Em = MV.cov_to_ellipse(Psi, 1.0); ax[k].plot(Em[0], Em[1], color=RED, lw=2.5)
ax[k].set_title(f"Inverse-Wishart covariance prior (ν={nu_})"); ax[k].set_aspect("equal"); ax[k].set_xlim(-4,4); ax[k].set_ylim(-3.5,3.5)
plt.tight_layout(); plt.show()
print("Each ellipse is a candidate covariance drawn from the prior; larger ν makes the prior more informative (tighter).")
pdf vs scipy (5 draws): 2.1316282072803006e-14 sample mean [[0.602 0.12 ] [0.12 0.4 ]] vs Ψ/(ν−p−1) [[0.6 0.12] [0.12 0.4 ]]
Each ellipse is a candidate covariance drawn from the prior; larger ν makes the prior more informative (tighter).
3. Normal-Inverse-Wishart — the joint prior for $(\boldsymbol\mu,\Sigma)$¶
When both the mean and covariance of a Gaussian are unknown, the conjugate prior is the Normal-Inverse-Wishart: draw the covariance from an Inverse-Wishart, then the mean from a Gaussian whose spread is that same covariance shrunk by $\kappa$, $$\Sigma\sim\text{IW}(\nu,\Psi),\qquad \boldsymbol\mu\mid\Sigma\sim N\!\big(\boldsymbol\mu_0,\ \Sigma/\kappa\big).$$ The coupling — $\boldsymbol\mu$'s uncertainty scales with the drawn $\Sigma$ — is what makes it conjugate: after seeing Gaussian data, the posterior is another NIW with updated parameters. Below, each draw is a $(\boldsymbol\mu,\Sigma)$ pair, shown as a centre point with its covariance ellipse.
mu0 = np.array([1.0, -0.5]); kappa = 1.0; Psi = np.array([[1.0, 0.3],[0.3, 1.0]])*4; nu = 8
mus, Sigs = MV.niw_rng(mu0, kappa, Psi, nu, 60, rng)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for m, Smat in zip(mus, Sigs):
E = MV.cov_to_ellipse(Smat, 1.0) + m[:, None]
ax[0].plot(E[0], E[1], color=BLUE, alpha=.25, lw=.7); ax[0].plot(*m, ".", color=RED, ms=4)
ax[0].plot(*mu0, "*", color="k", ms=14, label="μ₀"); ax[0].set_title("NIW: 60 draws of (μ, Σ)"); ax[0].legend(fontsize=8); ax[0].set_aspect("equal")
musB, SigsB = MV.niw_rng(mu0, kappa, Psi, nu, 100000, rng)
ax[1].scatter(musB[:,0], musB[:,1], s=2, alpha=.05, color=BLUE)
ax[1].plot(*mu0, "*", color=RED, ms=14); ax[1].set_title("Marginal spread of μ (integrating out Σ) — a multivariate t"); ax[1].set_aspect("equal"); ax[1].set_xlabel("μ₁"); ax[1].set_ylabel("μ₂")
plt.tight_layout(); plt.show()
print("E[μ] =", np.round(musB.mean(0),3), "= μ₀; E[Σ] =\n", np.round(SigsB.mean(0),3), "= IW mean.")
print("Marginally (Σ integrated out) μ follows a multivariate t — heavier-tailed than a Gaussian, the estimation-risk premium.")
E[μ] = [ 0.996 -0.504] = μ₀; E[Σ] = [[0.798 0.238] [0.238 0.798]] = IW mean. Marginally (Σ integrated out) μ follows a multivariate t — heavier-tailed than a Gaussian, the estimation-risk premium.
4. LKJ — the correlation-matrix prior¶
The LKJ(η) distribution (Lewandowski–Kurowicka–Joe) puts a prior directly on correlation matrices $R$ (unit diagonal, positive-definite), with density $f(R)\propto |R|^{\eta-1}$. The single shape $\eta$ tunes the strength of dependence it expects:
- $\eta=1$ — uniform over all valid correlation matrices;
- $\eta>1$ — concentrates near the identity (weak correlations), the usual weakly-informative choice;
- $\eta<1$ — favours strong correlations (mass near the boundary of the PD cone).
Its defining marginal is clean: any single off-diagonal correlation satisfies $(\rho+1)/2\sim\text{Beta}(a,a)$ with $a=\eta+(p-2)/2$, so $f(\rho)\propto(1-\rho^2)^{a-1}$. We sample by the C-vine method (partial correlations drawn from shifted Betas), and validate against that marginal.
# The signature LKJ figure: the 2x2 correlation density for varying eta
rr = np.linspace(-0.999, 0.999, 400)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
for eta,c in [(0.5,RED),(1.0,GREY),(2.0,GREEN),(5.0,BLUE)]:
ax[0].plot(rr, MV.lkj_marginal_rho_pdf(rr, eta, 2), color=c, lw=2.2, label=f"η={eta}")
if eta >= 1: # overlay a sample histogram where the density is bounded
R2 = MV.lkj_rng(eta, 2, 60000, rng)[:,0,1]
ax[0].hist(R2, bins=60, range=(-1,1), density=True, color=c, alpha=.12)
ax[0].set_title("LKJ correlation density on (−1,1) [p=2]"); ax[0].set_xlabel("ρ"); ax[0].legend(fontsize=8.5); ax[0].set_ylim(0,2.6)
# higher dimension: off-diagonal correlations shrink toward 0 as eta grows (p=6)
for eta,c in [(1.0,GREY),(2.0,GREEN),(10.0,BLUE)]:
R = MV.lkj_rng(eta, 6, 6000, rng)
off = np.array([Rm[np.triu_indices(6,1)] for Rm in R]).ravel()
ax[1].hist(off, bins=70, range=(-1,1), density=True, histtype="step", lw=2, color=c, label=f"η={eta}")
ax[1].set_title("LKJ off-diagonal correlations (p=6): larger η ⇒ nearer identity"); ax[1].set_xlabel("ρᵢⱼ"); ax[1].legend(fontsize=8.5)
plt.tight_layout(); plt.show()
print("η=1 is uniform (flat for p=2); η>1 pulls correlations toward 0; η<1 pushes them toward ±1. C-vine samples match the analytic marginal.")
η=1 is uniform (flat for p=2); η>1 pulls correlations toward 0; η<1 pushes them toward ±1. C-vine samples match the analytic marginal.
# Sampled correlation matrices as heatmaps (p=6), across eta
fig, ax = plt.subplots(1, 3, figsize=(14, 4.2))
for k,eta in enumerate([0.5, 1.0, 10.0]):
R = MV.lkj_rng(eta, 6, 1, rng)[0]
im = ax[k].imshow(R, cmap="RdBu_r", vmin=-1, vmax=1)
for i in range(6):
for j in range(6):
ax[k].text(j, i, f"{R[i,j]:.1f}", ha="center", va="center", fontsize=7, color="k")
ax[k].set_title(f"one LKJ(η={eta}) draw (p=6)"); ax[k].set_xticks([]); ax[k].set_yticks([]); plt.colorbar(im, ax=ax[k], shrink=.8)
plt.tight_layout(); plt.show()
print("A single draw at η=0.5 shows strong (dark) off-diagonals; at η=10 the matrix is nearly the identity.")
A single draw at η=0.5 shows strong (dark) off-diagonals; at η=10 the matrix is nearly the identity.
5. The modern covariance prior — LKJ + separate scales¶
The Inverse-Wishart's weakness is that one degrees-of-freedom parameter controls the uncertainty of every variance and correlation at once. The separation strategy (Barnard–McCulloch–Meng; the Stan/PyMC default) fixes this by splitting a covariance into scales × correlation: $$\Sigma=\text{diag}(\boldsymbol\tau)\,R\,\text{diag}(\boldsymbol\tau),\qquad \boldsymbol\tau\sim\text{(a scale prior)},\quad R\sim\text{LKJ}(\eta).$$ Now the standard deviations $\boldsymbol\tau$ get their own weakly-informative prior (a Half-Cauchy or Half-Normal from Folder 5) and the correlation structure gets an interpretable LKJ — independently. This is the recommended way to put a prior on a covariance matrix in modern Bayesian practice. Below we sample from it and show the induced cloud of covariance ellipses.
# Sigma = diag(tau) R diag(tau), tau ~ Half-Cauchy(1), R ~ LKJ(eta)
def modern_cov_rng(eta, p, scale, n, rng):
R = MV.lkj_rng(eta, p, n, rng)
tau = np.abs(scale * np.tan(np.pi*(rng.random((n,p))-0.5))) # Half-Cauchy scales
return np.array([np.diag(tau[i]) @ R[i] @ np.diag(tau[i]) for i in range(n)]), R, tau
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for k,eta in enumerate([1.0, 6.0]):
Sig, R, tau = modern_cov_rng(eta, 2, 1.0, 200, rng)
for Smat in Sig:
if np.all(np.linalg.eigvalsh(Smat) > 0) and np.trace(Smat) < 40:
E = MV.cov_to_ellipse(Smat, 1.0); ax[k].plot(E[0], E[1], color=BLUE, alpha=.12, lw=.7)
ax[k].set_title(f"LKJ(η={eta}) + Half-Cauchy scales — covariance prior"); ax[k].set_aspect("equal"); ax[k].set_xlim(-6,6); ax[k].set_ylim(-6,6)
plt.tight_layout(); plt.show()
print("η=1 admits any orientation/elongation; η=6 keeps the ellipses closer to axis-aligned (weak correlation).")
print("The Half-Cauchy scales let a few ellipses be large — the heavy-tailed scale prior in action (Folder 5).")
η=1 admits any orientation/elongation; η=6 keeps the ellipses closer to axis-aligned (weak correlation). The Half-Cauchy scales let a few ellipses be large — the heavy-tailed scale prior in action (Folder 5).
6. Timing¶
The Wishart (vectorised Bartlett) and Inverse-Wishart are fast batch operations. The LKJ C-vine and the NIW run a per-draw loop (a small Cholesky each), the honest cost of building a structured matrix one partial correlation, or one $(\boldsymbol\mu\mid\Sigma)$ pair, at a time.
import time
def best_time(fn, reps=3):
t = np.inf
for _ in range(reps):
s = time.perf_counter(); fn(); t = min(t, time.perf_counter()-s)
return t
V = np.array([[2.0,0.6],[0.6,1.0]]); Psi = np.array([[3.0,0.6],[0.6,2.0]])
benches = [
("Wishart(p=2)", lambda: MV.wishart_rng(8, V, 20000, rng), lambda: S.wishart(df=8, scale=V).rvs(20000, random_state=rng), 20000, "Bartlett (vectorised)"),
("Inv-Wishart(p=2)",lambda: MV.invwishart_rng(8, Psi, 20000, rng), lambda: S.invwishart(df=8, scale=Psi).rvs(20000, random_state=rng), 20000, "invert a Wishart"),
("LKJ(η=2, p=4)", lambda: MV.lkj_rng(2.0, 4, 4000, rng), None, 4000, "C-vine loop"),
("NIW(p=2)", lambda: MV.niw_rng(np.zeros(2),1,Psi,8, 4000, rng), None, 4000, "IW + per-draw μ"),
]
rows = []
for name, fs, sp, nn, kind in benches:
tf = best_time(fs); ts = best_time(sp) if sp else np.nan
rows.append((name, nn/tf/1e6, (nn/ts/1e6 if sp else np.nan), (ts/tf if sp else np.nan), kind))
df = pd.DataFrame(rows, columns=["dist","from-scratch (M/s)","scipy (M/s)","ratio","method"]).set_index("dist")
print(df.round(3).to_string())
# summary text computed from the table above -- never hardcode a timing claim
cmp_ = df.dropna(subset=["ratio"])
fmt = lambda sub: ", ".join("%s (%.1fx)" % (d, r) for d, r in sub["ratio"].items())
ahead = fmt(cmp_[cmp_["ratio"] >= 1]) or "none"
behind = fmt(cmp_[cmp_["ratio"] < 1]) or "none"
slowest = df["from-scratch (M/s)"].idxmin()
print("\nThe Bartlett decomposition vectorises across draws — one batched triangular construction rather than")
print("one per sample — so it runs well ahead of scipy's per-draw .rvs: " + ahead + ".")
print("Behind scipy: " + behind + ".")
print("The LKJ C-vine and the NIW cannot vectorise that way: each draw needs its own Cholesky, so both run")
print("as per-draw loops — the slowest here is %s at %.3f M draws/s." % (slowest, df["from-scratch (M/s)"].min()))
print("That is the natural cost of assembling a structured matrix one at a time, not a flaw in the algorithm.")
from-scratch (M/s) scipy (M/s) ratio method dist Wishart(p=2) 3.972 0.638 6.220 Bartlett (vectorised) Inv-Wishart(p=2) 3.277 0.532 6.163 invert a Wishart LKJ(η=2, p=4) 0.003 NaN NaN C-vine loop NIW(p=2) 0.074 NaN NaN IW + per-draw μ The Bartlett decomposition vectorises across draws — one batched triangular construction rather than one per sample — so it runs well ahead of scipy's per-draw .rvs: Wishart(p=2) (6.2x), Inv-Wishart(p=2) (6.2x). Behind scipy: none. The LKJ C-vine and the NIW cannot vectorise that way: each draw needs its own Cholesky, so both run as per-draw loops — the slowest here is LKJ(η=2, p=4) at 0.003 M draws/s. That is the natural cost of assembling a structured matrix one at a time, not a flaw in the algorithm.
7. Summary¶
Four matrix-variate laws — the distributions you place a prior on when the dependence structure itself is unknown — each built from scratch and validated against scipy.stats, and each made visible through covariance ellipses, the positive-definite cone, and correlation heatmaps.
The story. The Wishart is the scatter-matrix / precision law, sampled by the Bartlett decomposition and living in the PD cone with scaled-$\chi^2$ diagonals. Inverting it gives the Inverse-Wishart, the classical conjugate covariance prior, generalising the Inverse-Gamma; pairing that with a Gaussian mean gives the Normal-Inverse-Wishart, the joint conjugate prior for $(\boldsymbol\mu,\Sigma)$ whose marginal for $\boldsymbol\mu$ is a heavy-tailed multivariate $t$. The LKJ puts an interpretable prior on correlation matrices through a single shape $\eta$, sampled by the C-vine and validated on its exact off-diagonal marginal. And the separation strategy — LKJ correlations $\times$ independent scale priors — is the modern, recommended way to build a covariance prior, combining this folder with the Half-Cauchy scale priors of Folder 5.
Next in the catalog: compound & nonparametric distributions — the Tweedie family, the Dirichlet process and stick-breaking (GEM), and the Gaussian process, where the "parameter" becomes an entire function or an unbounded set of components.