Matrix-Variate & Covariance Distributions — the R engine

Statistical Distributions — priors on covariance & correlation matrices, in R

The independent R counterpart to matrixvar_python.ipynb. Each law's density, sampling algorithm and mean is built from scratch in base R (using chol, solve, lgamma). As a distribution over matrices is hard to picture, we draw each sampled covariance as its ellipse — a distribution over covariance matrices becomes a cloud of ellipses — and each correlation matrix as a heatmap.

Validation uses: the $p=1$ reduction (a scalar Wishart is a scaled $\chi^2$, checked against dchisq), sample-mean recovery ($\nu V$ for the Wishart, $\Psi/(\nu-p-1)$ for the Inverse-Wishart), an independent-sampler cross-check against base R's rWishart, and the LKJ off-diagonal marginal against its analytic Beta form. Self-contained; base R only.

options(repr.plot.width=13, repr.plot.height=4.6)
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"; PURP<-"#6b46c1"

norm_rng<-function(mu,s,n){ m<-ceiling(n/2); u1<-runif(m); u2<-runif(m); r<-sqrt(-2*log(u1)); z<-c(r*cos(2*pi*u2), r*sin(2*pi*u2)); mu+s*z[1:n] }
gam_shape<-function(k,n){ if(k<1){ g<-gam_shape(k+1,n); u<-runif(n); return(g*u^(1/k)) }
  d<-k-1/3; c<-1/sqrt(9*d); res<-numeric(n); todo<-rep(TRUE,n)
  while(any(todo)){ m<-sum(todo); x<-norm_rng(0,1,m); v<-(1+c*x)^3; u<-runif(m)
    ok<- v>0 & log(u) < 0.5*x^2 + d - d*v + d*log(ifelse(v>0,v,1)); idx<-which(todo); acc<-idx[ok]; res[acc]<-d*v[ok]; todo[acc]<-FALSE }
  res }
rbeta_s<-function(a,b){ x<-gam_shape(a,1); y<-gam_shape(b,1); x/(x+y) }

mvgammaln<-function(a,p){ j<-1:p; p*(p-1)/4*log(pi)+sum(lgamma(a+(1-j)/2)) }
wishart_logpdf<-function(W,nu,V){ p<-nrow(V); ldW<-as.numeric(determinant(W,log=TRUE)$modulus); ldV<-as.numeric(determinant(V,log=TRUE)$modulus); tr<-sum(diag(solve(V,W)))
  (nu-p-1)/2*ldW - tr/2 - nu*p/2*log(2) - nu/2*ldV - mvgammaln(nu/2,p) }
wishart_rng<-function(nu,V,n){ p<-nrow(V); L<-t(chol(V)); out<-array(0,c(p,p,n))
  for(r in 1:n){ A<-matrix(0,p,p); for(i in 1:p){ A[i,i]<-sqrt(2*gam_shape((nu-i+1)/2,1)); if(i>1) for(j in 1:(i-1)) A[i,j]<-norm_rng(0,1,1) }; M<-L%*%A; out[,,r]<-M%*%t(M) }; out }
iw_rng<-function(nu,Psi,n){ Wi<-wishart_rng(nu,solve(Psi),n); out<-array(0,dim(Wi)); for(r in 1:n) out[,,r]<-solve(Wi[,,r]); out }
niw_rng<-function(mu0,kappa,Psi,nu,n){ p<-length(mu0); Sig<-iw_rng(nu,Psi,n); mus<-matrix(0,n,p); for(r in 1:n){ L<-t(chol(Sig[,,r]/kappa)); mus[r,]<-mu0+L%*%norm_rng(0,1,p) }; list(mu=mus, Sigma=Sig) }

lkj_one<-function(eta,p){ P<-matrix(0,p,p); R<-diag(p); beta<-eta+(p-1)/2
  for(k in 1:(p-1)){ beta<-beta-0.5; for(i in (k+1):p){ P[k,i]<-2*rbeta_s(beta,beta)-1; pc<-P[k,i]
    if(k>1) for(l in (k-1):1) pc<-pc*sqrt((1-P[l,i]^2)*(1-P[l,k]^2))+P[l,i]*P[l,k]
    R[k,i]<-pc; R[i,k]<-pc } }; R }
lkj_rng<-function(eta,p,n){ out<-array(0,c(p,p,n)); for(r in 1:n) out[,,r]<-lkj_one(eta,p); out }
lkj_rho_pdf<-function(rho,eta,p){ a<-eta+(p-2)/2; ifelse(abs(rho)<1, (1-rho^2)^(a-1)*2^(1-2*a)/exp(2*lgamma(a)-lgamma(2*a)), 0) }

cov_ellipse<-function(S,nsd=1,np=100){ e<-eigen(S); t<-seq(0,2*pi,length=np); e$vectors%*%diag(nsd*sqrt(pmax(e$values,0)))%*%rbind(cos(t),sin(t)) }
cat("from-scratch R matrix-variate distributions ready\n")

from-scratch R matrix-variate distributions ready

1. Wishart — the scatter-matrix distribution

Bartlett sampling; density $\propto|W|^{(\nu-p-1)/2}\exp(-\tfrac12\mathrm{tr}(V^{-1}W))$; mean $\nu V$; diagonal is a scaled $\chi^2$. Validated at $p=1$ vs dchisq and cross-checked against base rWishart.

set.seed(1)
V<-matrix(c(2,0.6,0.6,1),2,2); nu<-8
# p=1 reduction: Wishart(nu, v) = v * chi^2_nu  -> compare to base dchisq
ww<-seq(0.5,30,length=200); v<-2
lp<-sapply(ww, function(w) wishart_logpdf(matrix(w,1,1), nu, matrix(v,1,1)))
cat(sprintf("p=1 Wishart vs scaled chi-square (dchisq): max|pdf diff| = %.2e\n", max(abs(exp(lp) - dchisq(ww/v,nu)/v))))
W<-wishart_rng(nu,V,1.5e5); cat("from-scratch sample mean:\n"); print(round(apply(W,c(1,2),mean),3)); cat("nu*V =\n"); print(nu*V)
Wb<-rWishart(1.5e5, nu, V); cat("base rWishart sample mean (independent cross-check):\n"); print(round(apply(Wb,c(1,2),mean),3))
cat(sprintf("W[1,1]: mean %.2f (nu*V11=%g), var %.1f (2*nu*V11^2=%g)\n", mean(W[1,1,]), nu*V[1,1], var(W[1,1,]), 2*nu*V[1,1]^2))

p=1 Wishart vs scaled chi-square (dchisq): max|pdf diff| = 9.02e-17

from-scratch sample mean:

       [,1]  [,2]
[1,] 16.015 4.808
[2,]  4.808 8.004

nu*V =

     [,1] [,2]
[1,] 16.0  4.8
[2,]  4.8  8.0

base rWishart sample mean (independent cross-check):

      [,1]  [,2]
[1,] 15.99 4.810
[2,]  4.81 8.015

W[1,1]: mean 16.02 (nu*V11=16), var 63.5 (2*nu*V11^2=64)

# The Wishart as a cloud of covariance ellipses; larger nu concentrates it
par(mfrow=c(1,2), mar=c(4,4,3,1))
V<-matrix(c(2,0.6,0.6,1),2,2); W<-wishart_rng(8,V,250)
plot(NULL, xlim=c(-3,3), ylim=c(-2.5,2.5), asp=1, xlab="x1", ylab="x2", main="Wishart(nu=8): 250 sampled covariance ellipses (/nu)")
for(r in 1:250){ E<-cov_ellipse(W[,,r]/8); lines(E[1,], E[2,], col=adjustcolor(BLUE,0.10)) }
Em<-cov_ellipse(V); lines(Em[1,], Em[2,], col=RED, lwd=2.5)
plot(NULL, xlim=c(-3,3), ylim=c(-2.5,2.5), asp=1, xlab="x1", ylab="x2", main="Larger nu concentrates the Wishart")
for(a in list(c(4,RED),c(20,GREEN),c(80,BLUE))){ Wn<-wishart_rng(as.numeric(a[1]),V,60); for(r in 1:60){ E<-cov_ellipse(Wn[,,r]/as.numeric(a[1])); lines(E[1,],E[2,], col=adjustcolor(a[2],0.10)) } }
legend("topright", c("nu=4","nu=20","nu=80"), col=c(RED,GREEN,BLUE), lwd=2, bty="n", cex=.8)
par(mfrow=c(1,1)); cat("Each ellipse is one sampled covariance; red is the mean V. Spread shrinks as 1/nu.\n")

Each ellipse is one sampled covariance; red is the mean V. Spread shrinks as 1/nu.

# diagonal marginal is a scaled chi-square
V<-matrix(c(1.5,0.3,0.3,1),2,2); nu<-6; W<-wishart_rng(nu,V,1.5e5)
hist(W[1,1,], breaks=90, freq=FALSE, col="#cfe3f6", border="white", main="Wishart diagonal W[1,1] is a scaled chi-square", xlab="W[1,1]")
tt<-seq(0,max(W[1,1,]),length=300); lines(tt, dchisq(tt/V[1,1],nu)/V[1,1], col=RED, lwd=2)
legend("topright", c("samples","V11 * chi^2_nu"), col=c("#cfe3f6",RED), lwd=c(6,2), bty="n")

2. Inverse-Wishart — the classical covariance prior

$W\sim\mathrm{Wishart}(\nu,\Psi^{-1})\Rightarrow W^{-1}\sim\mathrm{IW}(\nu,\Psi)$; the conjugate prior for a covariance, mean $\Psi/(\nu-p-1)$; the multivariate Inverse-Gamma.

Psi<-matrix(c(3,0.6,0.6,2),2,2); nu<-8; S<-iw_rng(nu,Psi,1.5e5)
cat("IW sample mean:\n"); print(round(apply(S,c(1,2),mean),3)); cat("Psi/(nu-p-1) =\n"); print(round(Psi/(nu-2-1),3))
par(mfrow=c(1,2), mar=c(4,4,3,1))
for(nu2 in c(6,30)){ Sc<-iw_rng(nu2, Psi*(nu2-3), 150)
  plot(NULL, xlim=c(-4,4), ylim=c(-3.5,3.5), asp=1, xlab="x1", ylab="x2", main=sprintf("Inverse-Wishart covariance prior (nu=%d)",nu2))
  for(r in 1:150){ E<-cov_ellipse(Sc[,,r]); lines(E[1,],E[2,], col=adjustcolor(BLUE,0.12)) }
  Em<-cov_ellipse(Psi); lines(Em[1,],Em[2,], col=RED, lwd=2.5) }
par(mfrow=c(1,1)); cat("Each ellipse is a covariance drawn from the prior; larger nu makes it tighter (more informative).\n")

IW sample mean:

      [,1]  [,2]
[1,] 0.599 0.119
[2,] 0.119 0.400

Psi/(nu-p-1) =

     [,1] [,2]
[1,] 0.60 0.12
[2,] 0.12 0.40

Each ellipse is a covariance drawn from the prior; larger nu makes it tighter (more informative).

3. Normal-Inverse-Wishart — the joint prior for (mu, Sigma)

$\Sigma\sim\mathrm{IW}(\nu,\Psi)$, $\mu\mid\Sigma\sim N(\mu_0,\Sigma/\kappa)$: the conjugate prior for a Gaussian's mean and covariance together.

mu0<-c(1,-0.5); Psi<-matrix(c(1,0.3,0.3,1),2,2)*4; nu<-8; d<-niw_rng(mu0,1,Psi,nu,60)
par(mfrow=c(1,2), mar=c(4,4,3,1))
plot(NULL, xlim=c(-4,6), ylim=c(-5,4), asp=1, xlab="x1", ylab="x2", main="NIW: 60 draws of (mu, Sigma)")
for(r in 1:60){ E<-cov_ellipse(d$Sigma[,,r]); lines(E[1,]+d$mu[r,1], E[2,]+d$mu[r,2], col=adjustcolor(BLUE,0.25)); points(d$mu[r,1],d$mu[r,2], pch=19, col=RED, cex=.5) }
points(mu0[1],mu0[2], pch=8, cex=2, lwd=2)
dB<-niw_rng(mu0,1,Psi,nu,6e4)
plot(dB$mu, pch=20, col=adjustcolor(BLUE,0.05), xlab="mu1", ylab="mu2", main="Marginal spread of mu (a multivariate t)"); points(mu0[1],mu0[2], pch=8, col=RED, cex=2, lwd=2)
par(mfrow=c(1,1)); cat("E[mu] =", round(colMeans(dB$mu),3), "= mu0; marginally mu is a heavy-tailed multivariate t.\n")

E[mu] = 1.007 -0.497 = mu0; marginally mu is a heavy-tailed multivariate t.

4. LKJ — the correlation-matrix prior

Density $\propto|R|^{\eta-1}$; $\eta=1$ uniform, $\eta>1$ near identity, $\eta<1$ strong correlations. Off-diagonal marginal $(\rho+1)/2\sim\mathrm{Beta}(a,a)$, $a=\eta+(p-2)/2$. Sampled by the C-vine.

par(mfrow=c(1,2), mar=c(4,4,3,1))
rr<-seq(-0.999,0.999,length=400)
plot(NULL, xlim=c(-1,1), ylim=c(0,2.6), xlab="rho", ylab="density", main="LKJ correlation density (p=2)")
i<-1; for(a in list(c(0.5,RED),c(1,GREY),c(2,GREEN),c(5,BLUE))){ eta<-as.numeric(a[1]); lines(rr, lkj_rho_pdf(rr,eta,2), col=a[2], lwd=2.2)
  if(eta>=1){ R2<-lkj_rng(eta,2,4e4); rho<-R2[1,2,]; h<-hist(rho, breaks=seq(-1,1,length=50), plot=FALSE); rect(h$breaks[-length(h$breaks)],0,h$breaks[-1],h$density, col=adjustcolor(a[2],0.12), border=NA) } }
legend("topright", c("eta=0.5","eta=1","eta=2","eta=5"), col=c(RED,GREY,GREEN,BLUE), lwd=2, bty="n", cex=.8)
plot(NULL, xlim=c(-1,1), ylim=c(0,3), xlab="rho_ij", ylab="density", main="LKJ off-diagonals (p=6): larger eta -> near identity")
for(a in list(c(1,GREY),c(2,GREEN),c(10,BLUE))){ R<-lkj_rng(as.numeric(a[1]),6,3000); off<-as.vector(apply(R,3,function(M) M[upper.tri(M)])); h<-hist(off, breaks=seq(-1,1,length=60), plot=FALSE); lines(h$mids, h$density, col=a[2], lwd=2) }
legend("topright", c("eta=1","eta=2","eta=10"), col=c(GREY,GREEN,BLUE), lwd=2, bty="n", cex=.8)
par(mfrow=c(1,1))
cat(sprintf("LKJ(2,p=2) rho marginal max|hist-analytic| = %.3f (bounded case).\n", { R2<-lkj_rng(2,2,6e4)[1,2,]; h<-hist(R2,breaks=40,plot=FALSE); max(abs(h$density-lkj_rho_pdf(h$mids,2,2))) }))

LKJ(2,p=2) rho marginal max|hist-analytic| = 0.040 (bounded case).

# sampled correlation matrices as heatmaps
par(mfrow=c(1,3), mar=c(2,2,3,1))
for(eta in c(0.5,1,10)){ R<-lkj_rng(eta,6,1)[,,1]
  image(1:6, 1:6, t(apply(R,2,rev)), col=hcl.colors(20,"Blue-Red"), zlim=c(-1,1), axes=FALSE, xlab="", ylab="", main=sprintf("one LKJ(eta=%g) draw (p=6)",eta))
  for(i in 1:6) for(j in 1:6) text(j, 7-i, sprintf("%.1f",R[i,j]), cex=.7) }
par(mfrow=c(1,1)); cat("eta=0.5 shows strong off-diagonals; eta=10 is nearly the identity.\n")

eta=0.5 shows strong off-diagonals; eta=10 is nearly the identity.

5. The modern covariance prior — LKJ + separate scales

The separation strategy: $\Sigma=\mathrm{diag}(\tau)\,R\,\mathrm{diag}(\tau)$ with $\tau$ from a scale prior (Half-Cauchy) and $R\sim\mathrm{LKJ}(\eta)$ — the Stan/PyMC default, giving scales and correlations independent, interpretable priors.

par(mfrow=c(1,2), mar=c(4,4,3,1))
for(eta in c(1,6)){ plot(NULL, xlim=c(-6,6), ylim=c(-6,6), asp=1, xlab="x1", ylab="x2", main=sprintf("LKJ(eta=%g) + Half-Cauchy scales",eta))
  for(r in 1:200){ R<-lkj_rng(eta,2,1)[,,1]; tau<-abs(tan(pi*(runif(2)-0.5))); S<-diag(tau)%*%R%*%diag(tau)
    if(all(eigen(S,only.values=TRUE)$values>0) && sum(diag(S))<40){ E<-cov_ellipse(S); lines(E[1,],E[2,], col=adjustcolor(BLUE,0.12)) } } }
par(mfrow=c(1,1)); cat("eta=1 admits any orientation; eta=6 keeps ellipses near axis-aligned. Half-Cauchy lets a few be large.\n")

eta=1 admits any orientation; eta=6 keeps ellipses near axis-aligned. Half-Cauchy lets a few be large.

6. Timing

Wishart / Inverse-Wishart / NIW / LKJ all assemble each matrix in a per-draw loop (a Cholesky each); the honest cost of building a structured positive-definite matrix. From-scratch is compared to base rWishart.

best <- function(f, min_t=0.05){ reps<-0L; s<-proc.time()[3]; repeat{ force(f()); reps<-reps+1L; el<-proc.time()[3]-s; if(el>=min_t && reps>=2L) break }; el/reps }
V<-matrix(c(2,0.6,0.6,1),2,2); Psi<-matrix(c(3,0.6,0.6,2),2,2)
b<-list(
 list("Wishart",   function() wishart_rng(8,V,3000),   function() rWishart(3000,8,V), 3000, "Bartlett loop"),
 list("InvWishart",function() iw_rng(8,Psi,3000),      function() NA,                3000, "invert Wishart"),
 list("LKJ(2,p=4)",function() lkj_rng(2,4,3000),       function() NA,                3000, "C-vine loop"),
 list("NIW",       function() niw_rng(c(0,0),1,Psi,8,3000), function() NA,           3000, "IW + per-draw mu"))
res<-data.frame()
for(x in b){ tf<-best(x[[2]]); has<-length(suppressWarnings(tryCatch(x[[3]](),error=function(e) NA)))>1
  tr<-if(has) best(x[[3]]) else NA; res<-rbind(res, data.frame(dist=x[[1]], scratch_Mps=x[[4]]/tf/1e6, base_Mps=if(has) x[[4]]/tr/1e6 else NA, method=x[[5]])) }
print(format(res, digits=3))
barplot(res$scratch_Mps, names.arg=res$dist, las=2, col=ifelse(is.na(res$base_Mps),ORANGE,BLUE), log="y", main="From-scratch matrix-variate RNG throughput (orange = no base equivalent)", ylab="million matrices / sec")
legend("topright", c("has base rWishart","no base equivalent"), fill=c(BLUE,ORANGE), bty="n", cex=.8)
cat("\nAll four assemble matrices in per-draw loops; base rWishart is compiled C, so it is faster on the Wishart.\n")

               dist scratch_Mps base_Mps           method
elapsed     Wishart      0.0500      3.6    Bartlett loop
elapsed1 InvWishart      0.0500       NA   invert Wishart
elapsed2 LKJ(2,p=4)      0.0122       NA      C-vine loop
elapsed3        NIW      0.0231       NA IW + per-draw mu

All four assemble matrices in per-draw loops; base rWishart is compiled C, so it is faster on the Wishart.

7. Summary

Every matrix-variate law's from-scratch R implementation is validated by the $p=1$ scaled-$\chi^2$ reduction (vs dchisq), sample-mean recovery, an independent cross-check against base rWishart, and the LKJ off-diagonal marginal — and each is drawn as covariance ellipses and correlation heatmaps. A cross-language corroboration of the Python results.

The Wishart (scatter / precision), the Inverse-Wishart (the classical covariance prior), their joint Normal-Inverse-Wishart, and the LKJ correlation prior — capped by the modern separation strategy that pairs LKJ correlations with independent scale priors. Next in the catalog: compound & nonparametric distributions (Tweedie, the Dirichlet process / stick-breaking, and the Gaussian process).