Multivariate Distributions¶
Statistical Distributions — joint laws and the dependence structure¶
Folder 6 of the Statistical Distributions catalog. A multivariate distribution describes a vector $\mathbf X=(X_1,\dots,X_k)$ jointly. The new object — absent from every univariate folder — is the dependence structure: a covariance (or precision) matrix, a shape vector, a spatial graph. This notebook builds the core multivariate laws from scratch — densities, sampling algorithms, marginals and conditionals — and validates against scipy.stats.
Because these live in two or three dimensions, we can see them. Throughout, each continuous law is drawn two ways: a 3-D density surface $z=f(x_1,x_2)$, and its iso-density contours — the curves along which the density is constant (the level sets). For the Gaussian and $t$ these level sets are ellipses whose orientation and elongation are the covariance; for the skew-normal they bend.
The organising ideas¶
- Elliptical laws. The multivariate Normal and multivariate $t$ depend on $\mathbf x$ only through the Mahalanobis quadratic form $(\mathbf x-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf x-\boldsymbol\mu)$, so their contours are concentric ellipses — the $t$ merely adds heavier tails. The multivariate skew-normal tilts them.
- Count vectors. The Multinomial is the joint law of category counts (the multivariate Binomial), and the Dirichlet-Multinomial is its over-dispersed compound — the discrete companions of the Dirichlet from Folder 5.
- Structured Gaussians. A Gaussian Markov random field (CAR) is a multivariate Normal whose precision matrix encodes a graph: zeros in $\Sigma^{-1}$ are conditional independences. This is the engine of spatial statistics.
- Copulas separate the dependence from the margins entirely — summarised here, with a cross-reference to the dedicated Copulas project.
Families covered¶
- Multivariate Normal — surfaces, contour ellipses, correlation, marginals & conditionals
- Multivariate Student-$t$ — elliptical heavy tails
- Multivariate Skew-Normal — tilted contours
- Multinomial & Dirichlet-Multinomial — count vectors and over-dispersion
- Gaussian Markov random fields / CAR — structured spatial Gaussians
- Copulas — the dependence structure, separated (reference)
- Timing
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa
import multivariate as M
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 surf_contour(title, pdf_fn, xr, yr, samples=None, npts=160, levels=9, cmap="viridis", log_levels=False):
'''3-D density surface + iso-density contour curves (the level sets) for a bivariate law.
log_levels=True places the contour levels geometrically in density — essential for
peaked / heavy-tailed laws (e.g. the multivariate t) where evenly-spaced linear levels
bunch every ring at the central spike and reveal nothing of the tails.'''
gx = np.linspace(*xr, npts); gy = np.linspace(*yr, npts)
GX, GY = np.meshgrid(gx, gy)
Z = pdf_fn(np.column_stack([GX.ravel(), GY.ravel()])).reshape(npts, npts)
lv = np.geomspace(Z.max() * 0.004, Z.max() * 0.95, 11) if log_levels else levels
fig = plt.figure(figsize=(13.5, 4.7))
ax1 = fig.add_subplot(1, 2, 1, projection="3d")
ax1.plot_surface(GX, GY, Z, cmap=cmap, rstride=2, cstride=2, linewidth=0, antialiased=True, alpha=.95)
ax1.set_title(title + " · density surface z=f(x₁,x₂)", fontsize=10)
ax1.set_xlabel("x₁"); ax1.set_ylabel("x₂"); ax1.set_zlabel("f"); ax1.view_init(28, -60)
ax2 = fig.add_subplot(1, 2, 2)
ax2.contourf(GX, GY, Z, levels=lv, cmap=cmap, alpha=.9, extend="max")
ax2.contour(GX, GY, Z, levels=lv, colors="white", linewidths=.6)
if samples is not None:
ax2.scatter(samples[:, 0], samples[:, 1], s=2, alpha=.10, color="white")
ax2.set_title("iso-density contours" + (" + samples" if samples is not None else "")); ax2.set_xlabel("x₁"); ax2.set_ylabel("x₂")
ax2.set_aspect("auto"); plt.tight_layout(); plt.show()
print("helpers ready")
helpers ready
1. Multivariate Normal — the elliptical anchor¶
$\mathbf X\sim N_k(\boldsymbol\mu,\Sigma)$ has density $$f(\mathbf x)=(2\pi)^{-k/2}|\Sigma|^{-1/2}\exp\!\Big[-\tfrac12(\mathbf x-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf x-\boldsymbol\mu)\Big].$$ It depends on $\mathbf x$ only through the Mahalanobis distance, so its level sets are ellipses centred at $\boldsymbol\mu$, with axes along the eigenvectors of $\Sigma$ and lengths $\propto\sqrt{\text{eigenvalues}}$ — the correlation literally rotates and stretches the contour. We sample by the Cholesky factor, $\mathbf x=\boldsymbol\mu+L\mathbf z$ with $\Sigma=LL^\top$, $\mathbf z\sim N(0,I)$. Two facts make the Gaussian the backbone of multivariate analysis and of Gibbs sampling: every marginal is Gaussian, and every conditional is Gaussian.
mu = np.array([0.5, -0.3]); Sig = np.array([[1.5, 0.8],[0.8, 1.0]])
X = M.mvn_rng(mu, Sig, 4000, rng)
surf_contour("Bivariate Normal", lambda P: M.mvn_pdf(P, mu, Sig), (-4,5), (-4,3.5), samples=X)
print("pdf vs scipy (max abs diff):", np.max(np.abs(M.mvn_pdf(X, mu, Sig) - S.multivariate_normal(mu, Sig).pdf(X))))
print("sample mean", np.round(X.mean(0),3), "vs", mu, " | sample cov\n", np.round(np.cov(X.T),3), "vs\n", Sig)
pdf vs scipy (max abs diff): 8.326672684688674e-17 sample mean [ 0.51 -0.288] vs [ 0.5 -0.3] | sample cov [[1.532 0.811] [0.811 1.005]] vs [[1.5 0.8] [0.8 1. ]]
# Feature: correlation rotates/stretches the contour ellipse; 95% probability ellipse
from matplotlib.patches import Ellipse
fig, ax = plt.subplots(1, 3, figsize=(14, 4.2))
for k,(rho,c) in enumerate([(-0.7,BLUE),(0.0,GREEN),(0.7,RED)]):
Sg = np.array([[1.0, rho],[rho, 1.0]])
Xs = M.mvn_rng(np.zeros(2), Sg, 3000, rng)
ax[k].scatter(Xs[:,0], Xs[:,1], s=3, alpha=.18, color=c)
vals, vecs = np.linalg.eigh(Sg); ang = np.degrees(np.arctan2(*vecs[:,1][::-1]))
for nsd,al in [(np.sqrt(5.991),.9)]: # chi2_2(0.95)=5.991 -> 95% ellipse
e = Ellipse((0,0), 2*nsd*np.sqrt(vals[1]), 2*nsd*np.sqrt(vals[0]), angle=ang, fill=False, color="k", lw=2)
ax[k].add_patch(e)
ax[k].set_title(f"ρ = {rho:+.1f} (95% ellipse)"); ax[k].set_xlim(-4,4); ax[k].set_ylim(-4,4); ax[k].set_aspect("equal")
plt.tight_layout(); plt.show()
print("The correlation ρ rotates the ellipse to ±45° and elongates it; the black curve is the 95% probability contour (χ²₂).")
The correlation ρ rotates the ellipse to ±45° and elongates it; the black curve is the 95% probability contour (χ²₂).
# Marginals and conditionals of a Gaussian are Gaussian (the Gibbs-sampling property)
mu = np.array([0.5, -0.3]); Sig = np.array([[1.5, 0.8],[0.8, 1.0]])
idx, cmean, ccov = M.mvn_conditional(mu, Sig, [1], [1.5]) # x0 | x1 = 1.5
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
X = M.mvn_rng(mu, Sig, 200000, rng)
tt = np.linspace(-4, 5, 300)
ax[0].hist(X[:,0], bins=120, density=True, color="#cfe3f6", edgecolor="white", lw=.2, label="marginal samples")
ax[0].plot(tt, S.norm.pdf(tt, mu[0], np.sqrt(Sig[0,0])), color=RED, lw=2, label=f"N({mu[0]}, {Sig[0,0]}) — marginal")
ax[0].set_title("Marginal of X₁ is Gaussian"); ax[0].legend(fontsize=8); ax[0].set_xlabel("x₁")
sel = np.abs(X[:,1]-1.5) < 0.05
ax[1].hist(X[sel,0], bins=60, density=True, color="#d8ecd8", edgecolor="white", lw=.2, label="samples with x₂≈1.5")
ax[1].plot(tt, S.norm.pdf(tt, cmean[0], np.sqrt(ccov[0,0])), color=RED, lw=2, label=f"N({cmean[0]:.2f}, {ccov[0,0]:.2f}) — conditional")
ax[1].set_title("Conditional X₁ | X₂=1.5 is Gaussian"); ax[1].legend(fontsize=8); ax[1].set_xlabel("x₁")
plt.tight_layout(); plt.show()
print(f"Conditional mean {cmean[0]:.3f}, variance {ccov[0,0]:.3f} — the closed forms behind every Gaussian Gibbs update.")
Conditional mean 1.940, variance 0.860 — the closed forms behind every Gaussian Gibbs update.
2. Multivariate Student-$t$ — elliptical, but heavy-tailed¶
$\mathbf X\sim t_k(\boldsymbol\mu,\Sigma,\nu)$ shares the Gaussian's elliptical contours but replaces the exponential decay with a polynomial one, so extreme joint events are far more likely — the standard robust/heavy-tailed multivariate model (and the copula behind the $t$-copula). It is a Gaussian scaled by an independent chi-square: $\mathbf X=\boldsymbol\mu+\sqrt{\nu/w}\,L\mathbf z$, $w\sim\chi^2_\nu$. Its covariance is $\frac{\nu}{\nu-2}\Sigma$ (for $\nu>2$); as $\nu\to\infty$ it becomes the Normal. Below, the surface and contours look Gaussian near the centre but the tail (low-density contours) is fatter.
mu = np.array([0.0, 0.0]); Sig = np.array([[1.0, 0.6],[0.6, 1.0]]); nu = 4
Xt = M.mvt_rng(mu, Sig, nu, 4000, rng)
surf_contour(f"Bivariate Student-t (ν={nu})", lambda P: M.mvt_pdf(P, mu, Sig, nu), (-6,6), (-6,6), samples=Xt, cmap="magma", log_levels=True)
print("pdf vs scipy (max abs diff):", np.max(np.abs(M.mvt_pdf(Xt, mu, Sig, nu) - S.multivariate_t(mu, Sig, df=nu).pdf(Xt))))
print("sample cov", np.round(np.cov(Xt.T),3), " vs ν/(ν-2)·Σ =", np.round(nu/(nu-2)*Sig,3).tolist())
pdf vs scipy (max abs diff): 1.1102230246251565e-16 sample cov [[2.136 1.345] [1.345 2.271]] vs ν/(ν-2)·Σ = [[2.0, 1.2], [1.2, 2.0]]
# Feature: MVN vs MVt tails -- a density cross-section on a log scale
tt = np.linspace(-7, 7, 400); pts = np.column_stack([tt, np.zeros_like(tt)])
Sig = np.array([[1.0,0.6],[0.6,1.0]])
plt.figure(figsize=(9.5,4.2))
plt.semilogy(tt, M.mvn_pdf(pts, np.zeros(2), Sig), color=BLUE, lw=2, label="Multivariate Normal")
for nu_,c in [(2,RED),(6,ORANGE),(30,GREEN)]:
plt.semilogy(tt, M.mvt_pdf(pts, np.zeros(2), Sig, nu_), lw=2, color=c, label=f"MV-t (ν={nu_})")
plt.title("Cross-section f(x₁, 0): the multivariate t has far heavier joint tails"); plt.xlabel("x₁"); plt.ylabel("density (log)")
plt.legend(fontsize=8.5); plt.ylim(1e-6, 1); plt.show()
print("Near the centre the t looks Gaussian; in the tails it decays polynomially, so joint extremes are much more likely.")
Near the centre the t looks Gaussian; in the tails it decays polynomially, so joint extremes are much more likely.
3. Multivariate Skew-Normal — tilting the ellipse¶
Azzalini's multivariate skew-normal multiplies a Gaussian density by a tilting factor, $$f(\mathbf x)=2\,\phi_k(\mathbf x-\boldsymbol\xi;\Omega)\,\Phi\!\big(\boldsymbol\alpha^\top\boldsymbol\omega^{-1}(\mathbf x-\boldsymbol\xi)\big),$$ so the shape vector $\boldsymbol\alpha$ skews the distribution along a direction while keeping Gaussian-light tails. Its contours are no longer symmetric ellipses — they bunch on one side. We sample it by the Azzalini–Dalla Valle construction $\mathbf Y=\boldsymbol\delta|U_0|+\mathbf V$ (the multivariate analogue of the $\delta|Z_0|+\sqrt{1-\delta^2}Z_1$ trick from Folder 3), and confirm the marginals are univariate skew-normals with the expected asymmetry.
xi = np.array([0.0, 0.0]); Om = np.array([[1.0, 0.5],[0.5, 1.0]]); alpha = np.array([5.0, -2.0])
Xs = M.mvskewnorm_rng(xi, Om, alpha, 4000, rng)
surf_contour("Bivariate Skew-Normal (α=[5,−2])", lambda P: M.mvskewnorm_pdf(P, xi, Om, alpha), (-3.5,4), (-4,3.5), samples=Xs, cmap="cividis")
from scipy.stats import skew
print("marginal skews:", np.round(skew(M.mvskewnorm_rng(xi,Om,alpha,200000,rng), axis=0),3),
"\n— component 1 is strongly right-skewed; component 2 stays near-symmetric because the correlation")
print("redistributes the skew (the marginal shape is set by Ω̄α = [4, 0.5], not by α itself). The joint contours are clearly non-elliptical.")
marginal skews: [ 0.449 -0.007] — component 1 is strongly right-skewed; component 2 stays near-symmetric because the correlation redistributes the skew (the marginal shape is set by Ω̄α = [4, 0.5], not by α itself). The joint contours are clearly non-elliptical.
4. Multinomial & Dirichlet-Multinomial — count vectors¶
Multinomial(n, p) is the joint law of the category counts in $n$ independent draws over $K$ classes — the multivariate Binomial and the likelihood of a contingency table. Each marginal $X_j\sim\text{Binomial}(n,p_j)$, the counts are negatively correlated (they sum to $n$), and it is conjugate to the Dirichlet. We sample it by sequential conditional Binomials (stick-breaking a count).
Dirichlet-Multinomial(n, α) is the compound where the probability vector is itself random, $\mathbf p\sim\text{Dirichlet}(\alpha)$: the mixing adds over-dispersion, the count analogue of the Beta-Binomial and the standard model for correlated or clustered categorical data.
# The Multinomial pmf over the simplex of counts (K=3, n=12), and Binomial marginals
n = 12; p = np.array([0.2, 0.5, 0.3])
pts = np.array([[i, j, n-i-j] for i in range(n+1) for j in range(n+1-i)])
prob = M.multinomial_pmf(pts, n, p)
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.4))
x2 = pts[:,1] + 0.5*pts[:,2]; y2 = (np.sqrt(3)/2)*pts[:,2]
sc = ax[0].scatter(x2, y2, c=prob, s=90, cmap="viridis", edgecolor="k", lw=.3)
ax[0].plot([0,n,n/2,0],[0,0,np.sqrt(3)/2*n,0],"k-",lw=1); ax[0].axis("off"); ax[0].set_aspect("equal")
ax[0].set_title("Multinomial(12, [.2,.5,.3]) pmf on the count simplex"); plt.colorbar(sc, ax=ax[0], shrink=.8)
Xm = M.multinomial_rng(n, p, 100000, rng)
ks = np.arange(0, n+1)
ax[1].bar(ks-0.15, np.bincount(Xm[:,1], minlength=n+1)/len(Xm), 0.3, color=BLUE, label="X₂ samples")
ax[1].plot(ks, S.binom.pmf(ks, n, p[1]), "o", color=RED, ms=5, label=f"Binomial(12, {p[1]}) — marginal")
ax[1].set_title("Marginal X₂ is Binomial"); ax[1].legend(fontsize=8); ax[1].set_xlabel("count")
plt.tight_layout(); plt.show()
print("pmf vs scipy:", np.max(np.abs(M.multinomial_pmf(pts,n,p) - np.array([S.multinomial(n,p).pmf(x) for x in pts]))))
pmf vs scipy: 1.8041124150158794e-16
# Over-dispersion: Dirichlet-Multinomial vs Multinomial at the SAME mean
n = 20; alpha = np.array([2.0, 5.0, 3.0]); pbar = alpha/alpha.sum()
Xm = M.multinomial_rng(n, pbar, 80000, rng)
Xdm = M.dirmult_rng(n, alpha, 80000, rng)
ks = np.arange(0, n+1)
plt.figure(figsize=(9.5,4.2))
plt.bar(ks-0.2, np.bincount(Xm[:,1], minlength=n+1)/len(Xm), 0.4, color=BLUE, alpha=.8, label=f"Multinomial var={Xm[:,1].var():.1f}")
plt.bar(ks+0.2, np.bincount(Xdm[:,1], minlength=n+1)/len(Xdm), 0.4, color=RED, alpha=.8, label=f"Dirichlet-Multinomial var={Xdm[:,1].var():.1f}")
plt.title("Same mean, more spread: the Dirichlet-Multinomial is over-dispersed"); plt.xlabel("count of category 2"); plt.legend(fontsize=9); plt.show()
print("pmf vs scipy dirichlet_multinomial:",
np.max(np.abs(M.dirmult_pmf(np.array([[6,8,6],[4,10,6]]),n,alpha) - np.array([float(S.dirichlet_multinomial(alpha,n).pmf(x)) for x in [[6,8,6],[4,10,6]]]))))
pmf vs scipy dirichlet_multinomial: 2.42861286636753e-17
5. Gaussian Markov random fields (CAR) — structured spatial Gaussians¶
A Gaussian Markov random field is a multivariate Normal defined through its precision matrix $Q=\Sigma^{-1}$ rather than its covariance. The point is that a zero in $Q$ is a conditional independence: $Q_{ij}=0$ means $X_i\perp X_j$ given all the others. The conditional autoregressive (CAR) model sets $Q=\tau(D-\rho W)$ from a neighbourhood graph $W$ — the workhorse of disease mapping and spatial econometrics. Sampling never forms the dense $\Sigma$: with $Q=LL^\top$ we solve $L^\top\mathbf x=\mathbf z$. The parameter $\rho$ tunes spatial smoothness — near $1$, neighbouring cells become strongly coupled.
# Sampled spatial fields on a 25x25 lattice for increasing spatial dependence rho
nr = nc = 25; W = M.grid_adjacency(nr, nc)
fig, ax = plt.subplots(1, 3, figsize=(14, 4.6))
for k,rho in enumerate([0.3, 0.9, 0.99]):
Q = M.car_precision(W, 1.0, rho)
field = M.gmrf_rng(Q, 1, rng)[0].reshape(nr, nc)
im = ax[k].imshow(field, cmap="RdBu_r", interpolation="nearest")
ax[k].set_title(f"CAR field, ρ = {rho}"); ax[k].axis("off"); plt.colorbar(im, ax=ax[k], shrink=.75)
plt.tight_layout(); plt.show()
print("Higher ρ couples neighbours more strongly, producing smoother spatial fields — the CAR prior's smoothing knob.")
Higher ρ couples neighbours more strongly, producing smoother spatial fields — the CAR prior's smoothing knob.
# The precision is sparse (the Markov graph); and precision-sampling recovers Q^{-1}
nr = nc = 6; W = M.grid_adjacency(nr, nc); Q = M.car_precision(W, 1.0, 0.9)
fig, ax = plt.subplots(1, 2, figsize=(12.5, 4.4))
ax[0].imshow(Q != 0, cmap="Greys", interpolation="nearest")
ax[0].set_title("Sparsity of the precision Q (a zero = conditional independence)"); ax[0].set_xlabel("node"); ax[0].set_ylabel("node")
Xg = M.gmrf_rng(Q, 150000, rng)
ax[1].scatter(np.linalg.inv(Q).ravel(), np.cov(Xg.T).ravel(), s=8, color=BLUE, alpha=.5)
lim = [np.linalg.inv(Q).min()-.05, np.linalg.inv(Q).max()+.05]; ax[1].plot(lim, lim, "--", color=RED, lw=1.2)
ax[1].set_title("Empirical covariance vs Q⁻¹"); ax[1].set_xlabel("Q⁻¹ entries"); ax[1].set_ylabel("empirical cov entries")
plt.tight_layout(); plt.show()
print("Dense covariance from a sparse precision: the GMRF trades a full Σ for a graph — the reason spatial models scale.")
Dense covariance from a sparse precision: the GMRF trades a full Σ for a graph — the reason spatial models scale.
6. Copulas — the dependence structure, separated (reference)¶
Everything above ties the margins and the dependence together in one object. A copula pulls them apart. By Sklar's theorem, any joint CDF factors as $$F(x_1,\dots,x_k)=C\big(F_1(x_1),\dots,F_k(x_k)\big),$$ where the margins $F_i$ are arbitrary and the copula $C$ — a distribution on the unit cube with uniform margins — carries all and only the dependence. This is what lets you pair, say, a heavy-tailed margin with a skewed one and a Gaussian dependence, independently. The illustration below shows the Gaussian copula density (the dependence of a bivariate Normal, with its margins transformed to uniform): the contours reveal the dependence alone.
This is a reference section. The full treatment — Gaussian, $t$, Clayton, Gumbel and Frank copulas, tail dependence, and copula-GARCH — lives in the dedicated Copulas & Tail Dependence project of the Risk & Asset Allocation series.
# Gaussian copula density c(u,v) = φ₂(Φ⁻¹u, Φ⁻¹v; ρ) / (φ(Φ⁻¹u) φ(Φ⁻¹v))
def gaussian_copula_density(U, V, rho):
a, b = S.norm.ppf(U), S.norm.ppf(V)
num = M.mvn_pdf(np.column_stack([a.ravel(), b.ravel()]), np.zeros(2), np.array([[1,rho],[rho,1]])).reshape(U.shape)
return num / (S.norm.pdf(a) * S.norm.pdf(b))
gu = np.linspace(0.02, 0.98, 120); GU, GV = np.meshgrid(gu, gu)
fig = plt.figure(figsize=(13.5, 4.7))
ax1 = fig.add_subplot(1,2,1, projection="3d")
Zc = gaussian_copula_density(GU, GV, 0.7)
ax1.plot_surface(GU, GV, Zc, cmap="viridis", rstride=2, cstride=2, linewidth=0, alpha=.95)
ax1.set_title("Gaussian copula density (ρ=0.7)", fontsize=10); ax1.set_xlabel("u"); ax1.set_ylabel("v"); ax1.view_init(30,-60)
ax2 = fig.add_subplot(1,2,2)
ax2.contourf(GU, GV, Zc, levels=12, cmap="viridis"); ax2.contour(GU, GV, Zc, levels=12, colors="white", linewidths=.5)
ax2.set_title("copula density on the unit square (uniform margins)"); ax2.set_xlabel("u"); ax2.set_ylabel("v")
plt.tight_layout(); plt.show()
print("With the margins transformed to uniform, only the dependence remains — the copula concentrates mass on the diagonal.")
With the margins transformed to uniform, only the dependence remains — the copula concentrates mass on the diagonal.
7. Timing¶
Multivariate sampling cost is dominated by one Cholesky factorisation (shared across all draws) plus a matrix multiply — so the from-scratch samplers run at essentially library speed. The GMRF adds a triangular solve; the compound Dirichlet-Multinomial pays for a per-draw Multinomial loop.
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
N = 100000
mu2 = np.zeros(2); S2 = np.array([[1.0,0.6],[0.6,1.0]]); al = np.array([2.,5.,3.]); pbar = al/al.sum()
Wg = M.grid_adjacency(8,8); Qg = M.car_precision(Wg,1.0,0.9)
benches = [
("MV-Normal", lambda: M.mvn_rng(mu2,S2,N,rng), lambda: S.multivariate_normal(mu2,S2).rvs(N,random_state=rng), "Cholesky"),
("MV-t(5)", lambda: M.mvt_rng(mu2,S2,5,N,rng), lambda: S.multivariate_t(mu2,S2,df=5).rvs(N,random_state=rng), "Chol + chi²"),
("Multinomial", lambda: M.multinomial_rng(20,pbar,N,rng), lambda: S.multinomial(20,pbar).rvs(N,random_state=rng), "cond. Binomials"),
("GMRF(64)", lambda: M.gmrf_rng(Qg,N//10,rng), None, "precision solve"),
("Dir-Mult", lambda: M.dirmult_rng(20,al,N//20,rng), None, "compound loop"),
]
rows = []
for name, fs, sp, kind in benches:
nn = {"GMRF(64)":N//10, "Dir-Mult":N//20}.get(name, N)
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(2).to_string())
# summary text computed from the table above -- never hardcode a timing claim
cmp_ = df.dropna(subset=["ratio"])
fmt = lambda sub: ", ".join("%s (%.2fx)" % (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("\nOne Cholesky is factored once and reused for every draw, so the elliptical samplers cost a matrix")
print("multiply per variate. Against scipy — ahead: " + ahead + "; behind: " + behind + ".")
print("The GMRF adds a triangular solve against a sparse precision, and the slowest in the set is "
+ "%s at %.2f M draws/s:" % (slowest, df["from-scratch (M/s)"].min()))
print("the Dirichlet-Multinomial draws a fresh Dirichlet then a Multinomial per observation, which is the")
print("honest cost of a compound law — it cannot reuse one factorisation the way the elliptical laws do.")
from-scratch (M/s) scipy (M/s) ratio method dist MV-Normal 18.50 28.24 0.66 Cholesky MV-t(5) 7.72 13.20 0.58 Chol + chi² Multinomial 10.06 9.96 1.01 cond. Binomials GMRF(64) 0.56 NaN NaN precision solve Dir-Mult 0.02 NaN NaN compound loop One Cholesky is factored once and reused for every draw, so the elliptical samplers cost a matrix multiply per variate. Against scipy — ahead: Multinomial (1.01x); behind: MV-Normal (0.66x), MV-t(5) (0.58x). The GMRF adds a triangular solve against a sparse precision, and the slowest in the set is Dir-Mult at 0.02 M draws/s: the Dirichlet-Multinomial draws a fresh Dirichlet then a Multinomial per observation, which is the honest cost of a compound law — it cannot reuse one factorisation the way the elliptical laws do.
8. Summary¶
Six multivariate laws, each built from scratch and validated against scipy.stats — and, crucially, drawn as 3-D density surfaces and iso-density contours so the dependence structure is visible.
The themes. The multivariate Normal and $t$ are elliptical: their contours are ellipses set by $\Sigma$, sampled through the Cholesky factor, with the $t$ adding heavy joint tails via a chi-square scale — and the Gaussian's marginals-and-conditionals-are-Gaussian property is exactly what powers Gibbs sampling. The skew-normal tilts those ellipses with a shape vector. The Multinomial and Dirichlet-Multinomial are the count-vector analogues of the Binomial and Beta-Binomial, the latter over-dispersed. The Gaussian Markov random field / CAR encodes a spatial graph in a sparse precision matrix — zeros are conditional independences — giving structured, scalable spatial Gaussians. And copulas separate dependence from margins entirely, cross-linked to the dedicated Copulas project.
Next in the catalog: matrix-variate & covariance distributions — the Wishart, Inverse-Wishart, Normal-Inverse-Wishart and LKJ laws that place priors on covariance and correlation matrices themselves.