Multivariate Distributions — the R engine

Statistical Distributions — joint laws and dependence, visualised in R

The independent R counterpart to multivariate_python.ipynb. Each law's density, sampling algorithm, marginals and conditionals is built from scratch in base R (using solve, chol, lgamma as primitives). As in Python, every continuous law is drawn two ways — a 3-D density surface (persp) and its iso-density contours (image + contour), the level sets along which the density is constant.

Validation uses base R's dmultinom/rmultinom (Multinomial), the independence check ($\rho=0$ Normal = product of dnorms), sample covariance recovery (MVN, MV-$t$, GMRF), and the Dirichlet-Multinomial over-dispersion. Self-contained; base R only.

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

# ---- RNG primitives ----
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 }

# ---- multivariate laws ----
mvn_pdf<-function(X,mu,Sig){ X<-matrix(X,ncol=length(mu)); k<-length(mu); Si<-solve(Sig); ld<-as.numeric(determinant(Sig,logarithm=TRUE)$modulus)
  d<-sweep(X,2,mu); quad<-rowSums((d%*%Si)*d); exp(-0.5*(k*log(2*pi)+ld+quad)) }
mvn_rng<-function(mu,Sig,n){ k<-length(mu); U<-chol(Sig); Z<-matrix(norm_rng(0,1,n*k),n,k); sweep(Z%*%U,2,mu,"+") }
mvt_pdf<-function(X,mu,Sig,nu){ X<-matrix(X,ncol=length(mu)); k<-length(mu); Si<-solve(Sig); ld<-as.numeric(determinant(Sig,logarithm=TRUE)$modulus)
  d<-sweep(X,2,mu); quad<-rowSums((d%*%Si)*d); exp(lgamma((nu+k)/2)-lgamma(nu/2)-0.5*(k*log(nu*pi)+ld)-(nu+k)/2*log1p(quad/nu)) }
mvt_rng<-function(mu,Sig,nu,n){ k<-length(mu); U<-chol(Sig); Z<-matrix(norm_rng(0,1,n*k),n,k); w<-2*gam_shape(nu/2,n); sweep(sqrt(nu/w)*(Z%*%U),2,mu,"+") }
mvn_cond<-function(mu,Sig,io,xo){ k<-length(mu); iu<-setdiff(1:k,io); Sbb<-Sig[io,io,drop=FALSE]; Sab<-Sig[iu,io,drop=FALSE]
  list(mean=mu[iu]+Sab%*%solve(Sbb)%*%(xo-mu[io]), cov=Sig[iu,iu,drop=FALSE]-Sab%*%solve(Sbb)%*%t(Sab)) }
msn_pdf<-function(X,xi,Om,al){ X<-matrix(X,ncol=length(xi)); om<-sqrt(diag(Om)); z<-sweep(X,2,xi); zz<-sweep(z,2,om,"/")
  2*mvn_pdf(X,xi,Om)*pnorm(as.numeric(zz%*%al)) }
msn_rng<-function(xi,Om,al,n){ om<-sqrt(diag(Om)); k<-length(xi); Ob<-Om/outer(om,om); den<-sqrt(1+as.numeric(t(al)%*%Ob%*%al)); delta<-as.numeric(Ob%*%al)/den
  Cv<-Ob-outer(delta,delta); Uv<-chol(Cv+1e-12*diag(k)); u0<-abs(norm_rng(0,1,n)); ZZ<-matrix(norm_rng(0,1,n*k),n,k)
  Y<-outer(u0,delta)+ZZ%*%Uv; sweep(sweep(Y,2,om,"*"),2,xi,"+") }

# multinomial (base dmultinom/rmultinom as reference)
mult_rng<-function(n,p,size){ K<-length(p); out<-matrix(0L,size,K); rem<-rep(n,size); prem<-1
  for(j in 1:(K-1)){ pj<-min(max(p[j]/prem,0),1); xj<-rbinom(size,rem,pj); out[,j]<-xj; rem<-rem-xj; prem<-prem-p[j] }; out[,K]<-rem; out }
mult_pmf<-function(x,n,p){ x<-matrix(x,ncol=length(p)); exp(lgamma(n+1)-rowSums(lgamma(x+1))+rowSums(sweep(x,2,log(p),"*")))*(rowSums(x)==n) }
dm_pmf<-function(x,n,al){ x<-matrix(x,ncol=length(al)); A<-sum(al); exp(lgamma(n+1)-rowSums(lgamma(x+1))+lgamma(A)-lgamma(n+A)+rowSums(sweep(lgamma(sweep(x,2,al,"+")),2,lgamma(al),"-")))*(rowSums(x)==n) }
dm_rng<-function(n,al,size){ K<-length(al); G<-matrix(0,size,K); for(i in 1:K) G[,i]<-gam_shape(al[i],size); P<-G/rowSums(G)
  out<-matrix(0L,size,K); for(r in 1:size) out[r,]<-mult_rng(n,P[r,],1); out }

# GMRF / CAR
grid_adj<-function(nr,nc){ N<-nr*nc; W<-matrix(0,N,N)
  for(i in 0:(nr-1)) for(j in 0:(nc-1)){ a<-i*nc+j+1; for(dd in list(c(-1,0),c(1,0),c(0,-1),c(0,1))){ ii<-i+dd[1]; jj<-j+dd[2]; if(ii>=0&ii<nr&jj>=0&jj<nc) W[a, ii*nc+jj+1]<-1 } }; W }
car_prec<-function(W,tau,rho){ D<-diag(rowSums(W)); tau*(D-rho*W) }
gmrf_rng<-function(Q,n){ N<-nrow(Q); U<-chol(Q); Z<-matrix(norm_rng(0,1,n*N),N,n); t(backsolve(U,Z)) }

surf_contour<-function(title, pdffun, xr, yr, np=90, log_levels=FALSE){
  gx<-seq(xr[1],xr[2],length=np); gy<-seq(yr[1],yr[2],length=np)
  G<-as.matrix(expand.grid(x=gx,y=gy)); Z<-matrix(pdffun(G), np, np)
  par(mfrow=c(1,2), mar=c(1,1,3,1))
  persp(gx,gy,Z, theta=-55, phi=28, col="#9ecae1", border=NA, shade=.4, xlab="x1", ylab="x2", zlab="f", main=paste0(title," | density surface"))
  par(mar=c(4,4,3,1)); image(gx,gy,Z, col=hcl.colors(24,"YlGnBu",rev=TRUE), xlab="x1", ylab="x2", main="iso-density contours")
  # geometric (log) density levels so a peaked / heavy-tailed law (e.g. the t) shows tail contours, not just a central blob
  lv <- if(log_levels) exp(seq(log(max(Z)*0.004), log(max(Z)*0.95), length=11)) else pretty(range(Z), 10)
  contour(gx,gy,Z, add=TRUE, col="white", drawlabels=FALSE, levels=lv)
  par(mfrow=c(1,1))
}
cat("from-scratch R multivariate distributions + surf_contour() ready\n")

from-scratch R multivariate distributions + surf_contour() ready

1. Multivariate Normal — the elliptical anchor

Density surface and iso-density ellipses; correlation rotates the ellipse; sampling by Cholesky; Gaussian marginals & conditionals.

set.seed(1)
mu<-c(0.5,-0.3); Sig<-matrix(c(1.5,0.8,0.8,1.0),2,2)
surf_contour("Bivariate Normal", function(P) mvn_pdf(P,mu,Sig), c(-4,5), c(-4,3.5))
X<-mvn_rng(mu,Sig,3e5)
cat(sprintf("rho=0 independence check: max|mvn_pdf - dnorm*dnorm| = %.2e\n",
    { g<-as.matrix(expand.grid(seq(-3,3,length=40),seq(-3,3,length=40))); max(abs(mvn_pdf(g,c(0,0),diag(2)) - dnorm(g[,1])*dnorm(g[,2]))) }))
cat("sample mean", round(colMeans(X),3), " (true 0.5,-0.3); sample cov diag", round(diag(cov(X)),3), "off", round(cov(X)[1,2],3), "(true 1.5,1.0,0.8)\n")

rho=0 independence check: max|mvn_pdf - dnorm*dnorm| = 4.16e-17

sample mean 0.497 -0.302  (true 0.5,-0.3); sample cov diag 1.497 1.004 off 0.802 (true 1.5,1.0,0.8)

# correlation rotates/stretches the 95% probability ellipse
par(mfrow=c(1,3), mar=c(4,4,3,1))
th<-seq(0,2*pi,length=100)
for(rho in c(-0.7,0,0.7)){ Sg<-matrix(c(1,rho,rho,1),2,2); Xs<-mvn_rng(c(0,0),Sg,3000)
  plot(Xs, pch=20, col=adjustcolor(BLUE,0.2), cex=.5, xlim=c(-4,4), ylim=c(-4,4), asp=1, xlab="x1", ylab="x2", main=sprintf("rho = %+.1f (95%% ellipse)",rho))
   e<-eigen(Sg); A<-e$vectors%*%diag(sqrt(e$values*5.991))%*%rbind(cos(th),sin(th)); lines(A[1,],A[2,], lwd=2) }
par(mfrow=c(1,1))
cat("The correlation rotates the ellipse to +/-45 deg; the black curve is the 95% contour (chi-square_2 = 5.991).\n")

The correlation rotates the ellipse to +/-45 deg; the black curve is the 95% contour (chi-square_2 = 5.991).

# marginal & conditional of a Gaussian are Gaussian
mu<-c(0.5,-0.3); Sig<-matrix(c(1.5,0.8,0.8,1.0),2,2); X<-mvn_rng(mu,Sig,2e5)
cn<-mvn_cond(mu,Sig,2,1.5)
par(mfrow=c(1,2), mar=c(4,4,3,1)); tt<-seq(-4,5,length=300)
hist(X[,1], breaks=100, freq=FALSE, col="#cfe3f6", border="white", main="Marginal of X1 is Gaussian", xlab="x1"); lines(tt, dnorm(tt,mu[1],sqrt(Sig[1,1])), col=RED, lwd=2)
sel<-abs(X[,2]-1.5)<0.05; hist(X[sel,1], breaks=50, freq=FALSE, col="#d8ecd8", border="white", main="Conditional X1 | X2=1.5 is Gaussian", xlab="x1"); lines(tt, dnorm(tt,cn$mean,sqrt(cn$cov)), col=RED, lwd=2)
par(mfrow=c(1,1)); cat(sprintf("Conditional mean %.3f, variance %.3f -- the closed forms behind Gaussian Gibbs updates.\n", cn$mean, cn$cov))

Conditional mean 1.940, variance 0.860 -- the closed forms behind Gaussian Gibbs updates.

2. Multivariate Student-t — elliptical, heavy-tailed

Same elliptical contours as the Normal but polynomial tails; sampled as a Gaussian scaled by a chi-square; covariance $\nu/(\nu-2)\,\Sigma$.

mu<-c(0,0); Sig<-matrix(c(1,0.6,0.6,1),2,2); nu<-4
surf_contour(sprintf("Bivariate Student-t (nu=%d)",nu), function(P) mvt_pdf(P,mu,Sig,nu), c(-6,6), c(-6,6), log_levels=TRUE)
Xt<-mvt_rng(mu,Sig,6,3e5); cat("MV-t(6) sample cov\n"); print(round(cov(Xt),3)); cat("true nu/(nu-2)*Sigma =", round(6/4*Sig,3),"\n")
# tail comparison cross-section (log scale)
tt<-seq(-7,7,length=400); pts<-cbind(tt,0)
plot(tt, mvn_pdf(pts,c(0,0),Sig), type="l", log="y", lwd=2, col=BLUE, ylim=c(1e-6,1), xlab="x1", ylab="density (log)", main="Cross-section f(x1,0): MV-t has heavier joint tails")
for(a in list(c(2,RED),c(6,ORANGE),c(30,GREEN))) lines(tt, mvt_pdf(pts,c(0,0),Sig,as.numeric(a[1])), lwd=2, col=a[2])
legend("bottom", c("Normal","t(2)","t(6)","t(30)"), col=c(BLUE,RED,ORANGE,GREEN), lwd=2, bty="n", cex=.8, horiz=TRUE)

MV-t(6) sample cov

      [,1]  [,2]
[1,] 1.504 0.903
[2,] 0.903 1.506

true nu/(nu-2)*Sigma = 1.5 0.9 0.9 1.5 

3. Multivariate Skew-Normal — tilted contours

Azzalini's shape vector skews the Gaussian along a direction; the contours are no longer symmetric ellipses. Sampled by the Azzalini–Dalla Valle construction.

xi<-c(0,0); Om<-matrix(c(1,0.5,0.5,1),2,2); al<-c(5,-2)
surf_contour("Bivariate Skew-Normal (alpha=[5,-2])", function(P) msn_pdf(P,xi,Om,al), c(-3.5,4), c(-4,3.5))
Xs<-msn_rng(xi,Om,al,2e5); sk<-function(x){ m<-mean(x); v<-mean((x-m)^2); mean((x-m)^3)/v^1.5 }
cat("marginal skews:", round(apply(Xs,2,sk),3), "\n-- component 1 is strongly right-skewed; component 2 stays near-symmetric because the correlation\n")
cat("redistributes the skew (the marginal shape follows Omega_bar %*% alpha = [4, 0.5], not alpha itself). Joint contours are non-elliptical.\n")

marginal skews: 0.463 -0.002 
-- component 1 is strongly right-skewed; component 2 stays near-symmetric because the correlation

redistributes the skew (the marginal shape follows Omega_bar %*% alpha = [4, 0.5], not alpha itself). Joint contours are non-elliptical.

4. Multinomial & Dirichlet-Multinomial — count vectors

Multinomial (vs base dmultinom/rmultinom); Dirichlet-Multinomial adds over-dispersion at the same mean.

n<-12; p<-c(0.2,0.5,0.3)
pts<-do.call(rbind, unlist(lapply(0:n, function(i) lapply(0:(n-i), function(j) c(i,j,n-i-j))), recursive=FALSE))
cat("Multinomial pmf vs base dmultinom: max diff =", max(abs(mult_pmf(pts,n,p) - apply(pts,1,function(x) dmultinom(x,prob=p)))), "\n")
# simplex heatmap of the pmf + Binomial marginal
par(mfrow=c(1,2), mar=c(2,2,3,1)); prob<-mult_pmf(pts,n,p)
x2<-pts[,2]+0.5*pts[,3]; y2<-(sqrt(3)/2)*pts[,3]; cols<-hcl.colors(20,"YlGnBu",rev=TRUE)[cut(prob,20)]
plot(x2,y2, pch=19, col=cols, cex=2, axes=FALSE, xlab="", ylab="", main="Multinomial(12,[.2,.5,.3]) pmf on simplex", asp=1); lines(c(0,n,n/2,0),c(0,0,sqrt(3)/2*n,0))
par(mar=c(4,4,3,1)); Xm<-mult_rng(n,p,1e5); ks<-0:n
plot(ks, tabulate(Xm[,2]+1,n+1)/nrow(Xm), type="h", lwd=4, col=BLUE, xlab="count", ylab="prob", main="Marginal X2 is Binomial"); points(ks, dbinom(ks,n,p[2]), pch=19, col=RED)
par(mfrow=c(1,1))

Multinomial pmf vs base dmultinom: max diff = 1.249001e-16 

# over-dispersion: Dirichlet-Multinomial vs Multinomial, same mean
n<-20; al<-c(2,5,3); pbar<-al/sum(al)
Xm<-mult_rng(n,pbar,8e4); Xdm<-dm_rng(n,al,8e4); ks<-0:n
plot(ks-0.2, tabulate(Xm[,2]+1,n+1)/nrow(Xm), type="h", lwd=4, col=BLUE, xlab="count of category 2", ylab="prob", main=sprintf("Over-dispersion: Multinomial var %.1f vs Dir-Mult var %.1f", var(Xm[,2]), var(Xdm[,2])))
lines(ks+0.2, tabulate(Xdm[,2]+1,n+1)/nrow(Xdm), type="h", lwd=4, col=RED)
legend("topright", c("Multinomial","Dirichlet-Multinomial"), col=c(BLUE,RED), lwd=3, bty="n", cex=.8)
# actually compute the closed form and compare, rather than printing one number
xchk<-c(6,8,6)
manual<-exp(lgamma(n+1)-sum(lgamma(xchk+1))+lgamma(sum(al))-lgamma(n+sum(al))+sum(lgamma(xchk+al)-lgamma(al)))
got<-dm_pmf(rbind(xchk),n,al)
cat(sprintf("Dir-Mult pmf at (6,8,6): from-scratch %.6f vs closed form %.6f  |diff| = %.2e\n",
            got, manual, abs(got-manual)))

Dir-Mult pmf at (6,8,6): from-scratch 0.009687 vs closed form 0.009687  |diff| = 0.00e+00

5. Gaussian Markov random fields (CAR) — structured spatial Gaussians

A multivariate Normal defined by a sparse precision $Q=\tau(D-\rho W)$: zeros in $Q$ are conditional independences. Sampled by solving $U\mathbf x=\mathbf z$ with $Q=U^\top U$. $\rho$ tunes spatial smoothness.

nr<-nc<-25; W<-grid_adj(nr,nc)
par(mfrow=c(1,3), mar=c(2,2,3,1))
for(rho in c(0.3,0.9,0.99)){ Q<-car_prec(W,1,rho); f<-matrix(gmrf_rng(Q,1)[1,], nr, nc, byrow=TRUE)
  image(1:nr,1:nc,f, col=hcl.colors(30,"Blue-Red"), axes=FALSE, xlab="", ylab="", main=sprintf("CAR field, rho = %.2f",rho)) }
par(mfrow=c(1,1)); cat("Higher rho couples neighbours more strongly, producing smoother spatial fields.\n")

Higher rho couples neighbours more strongly, producing smoother spatial fields.

nr<-nc<-6; W<-grid_adj(nr,nc); Q<-car_prec(W,1,0.9)
par(mfrow=c(1,2), mar=c(4,4,3,1))
image(1:nrow(Q), 1:ncol(Q), t(apply(Q!=0,2,rev)), col=c("white","black"), xlab="node", ylab="node", main="Sparsity of precision Q (0 = cond. independence)")
Xg<-gmrf_rng(Q,1.5e5); plot(as.vector(solve(Q)), as.vector(cov(Xg)), pch=20, col=adjustcolor(BLUE,.5), cex=.6, xlab="Q^-1 entries", ylab="empirical cov", main="Empirical covariance vs Q^-1"); abline(0,1,col=RED,lty=2,lwd=1.2)
par(mfrow=c(1,1)); cat("A dense covariance from a sparse precision -- why spatial GMRF models scale.\n")

A dense covariance from a sparse precision -- why spatial GMRF models scale.

6. Copulas — the dependence structure, separated (reference)

By Sklar's theorem $F(x_1,x_2)=C(F_1(x_1),F_2(x_2))$: the copula $C$ carries all the dependence with uniform margins. The Gaussian copula density below is that dependence alone. Full treatment in the dedicated Copulas & Tail Dependence project.

gauss_cop<-function(U,V,rho){ a<-qnorm(U); b<-qnorm(V); mvn_pdf(cbind(a,b),c(0,0),matrix(c(1,rho,rho,1),2,2))/(dnorm(a)*dnorm(b)) }
gu<-seq(0.02,0.98,length=70); G<-as.matrix(expand.grid(u=gu,v=gu)); Z<-matrix(gauss_cop(G[,1],G[,2],0.7), 70, 70)
par(mfrow=c(1,2), mar=c(1,1,3,1)); persp(gu,gu,Z, theta=-55, phi=30, col="#a1d99b", border=NA, shade=.4, xlab="u", ylab="v", zlab="c", main="Gaussian copula density (rho=0.7)")
par(mar=c(4,4,3,1)); image(gu,gu,Z, col=hcl.colors(24,"YlGnBu",rev=TRUE), xlab="u", ylab="v", main="copula density (uniform margins)"); contour(gu,gu,Z, add=TRUE, col="white", drawlabels=FALSE)
par(mfrow=c(1,1)); cat("With margins transformed to uniform, only the dependence remains -- mass concentrates on the diagonal.\n")

With margins transformed to uniform, only the dependence remains -- mass concentrates on the diagonal.

7. Summary

Every multivariate law's from-scratch R implementation is validated against base R where possible (dmultinom), by the $\rho=0$ independence check and sample-covariance recovery (Normal, $t$, GMRF), and by the Dirichlet-Multinomial over-dispersion — and, crucially, drawn as 3-D persp surfaces and image/contour level sets. A cross-language corroboration of the Python results.

The elliptical Normal and $t$, the tilted skew-normal, the count-vector Multinomial and its over-dispersed Dirichlet-Multinomial compound, the sparse-precision spatial GMRF, and the copula view of dependence. Next in the catalog: matrix-variate & covariance distributions (Wishart, Inverse-Wishart, LKJ).