Compound & Nonparametric Distributions — the R engine¶

Statistical Distributions — the parameter becomes a distribution or a function (in R)¶

The independent R counterpart to compound_np_python.ipynb, and the final notebook of the catalog. Each object is built from scratch in base R and validated through its construction and moments (these are processes, not densities): the Tweedie (compound Poisson-Gamma) by its mean, variance and mass at zero; the Dirichlet process by its stick-breaking weights and Chinese-Restaurant cluster counts; and the Gaussian process by its kernel-defined covariance and closed-form posterior.

As in Python we draw them — sampled loss distributions, sampled random measures, sampled random functions (including a 2-D-input GP surface via persp). Self-contained; base R only.

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

norm_rng<-function(mu,s,n){ if(n==0) return(numeric(0)); 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(n==0) return(numeric(0)); 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 }

# Tweedie (compound Poisson-Gamma)
tw_params<-function(mu,phi,p) list(lam=mu^(2-p)/(phi*(2-p)), a=(2-p)/(p-1), theta=phi*(p-1)*mu^(p-1))
tw_rng<-function(mu,phi,p,n){ pr<-tw_params(mu,phi,p); N<-rpois(n,pr$lam); Nt<-sum(N)
  jumps<-if(Nt>0) pr$theta*gam_shape(pr$a,Nt) else numeric(0); Y<-numeric(n); idx<-rep(1:n,N); if(Nt>0) for(t in seq_along(idx)) Y[idx[t]]<-Y[idx[t]]+jumps[t]; Y }
tw_p0<-function(mu,phi,p) exp(-tw_params(mu,phi,p)$lam)

# GEM stick-breaking + CRP
gem<-function(alpha,K,n){ b<-1-matrix(runif(n*K),n,K)^(1/alpha); rem<-t(apply(1-b,1,cumprod)); rs<-cbind(1, rem[,-K,drop=FALSE]); b*rs }
crp<-function(alpha,n){ sizes<-numeric(0); assign<-integer(n)
  for(i in 1:n){ w<-c(sizes,alpha); w<-w/sum(w); j<-findInterval(runif(1),cumsum(w))+1
    if(j>length(sizes)) sizes<-c(sizes,1) else sizes[j]<-sizes[j]+1; assign[i]<-j }; list(assign=assign, sizes=sizes) }
dp_ecl<-function(alpha,n) sum(alpha/(alpha+(1:n)-1))

# Gaussian process
rbf_k<-function(X1,X2,ls=1,v=1){ d2<-outer(rowSums(X1^2),rowSums(X2^2),"+")-2*X1%*%t(X2); v*exp(-0.5*pmax(d2,0)/ls^2) }
mat32_k<-function(X1,X2,ls=1,v=1){ d2<-outer(rowSums(X1^2),rowSums(X2^2),"+")-2*X1%*%t(X2); s<-sqrt(3)*sqrt(pmax(d2,0))/ls; v*(1+s)*exp(-s) }
mat12_k<-function(X1,X2,ls=1,v=1){ d2<-outer(rowSums(X1^2),rowSums(X2^2),"+")-2*X1%*%t(X2); v*exp(-sqrt(pmax(d2,0))/ls) }
gp_prior<-function(X,ker,n,jit=1e-8){ m<-nrow(X); K<-ker(X,X)+jit*diag(m); Z<-matrix(norm_rng(0,1,n*m),n,m); Z%*%chol(K) }
gp_post<-function(Xtr,ytr,Xte,ker,noise=1e-6){ K<-ker(Xtr,Xtr)+noise*diag(nrow(Xtr)); Ks<-ker(Xtr,Xte); Kss<-ker(Xte,Xte)
  # solve through the Cholesky factor, never forming K^{-1} (Rasmussen & Williams Alg. 2.1)
  L<-t(chol(K)); a<-backsolve(t(L), forwardsolve(L, ytr)); V<-forwardsolve(L, Ks)
  list(mean=as.vector(t(Ks)%*%a), cov=Kss-t(V)%*%V) }
cat("from-scratch R compound & nonparametric objects ready\n")

from-scratch R compound & nonparametric objects ready

1. Compound distributions — the Tweedie (compound Poisson-Gamma)¶

$N\sim\mathrm{Poisson}(\lambda)$, $Y=\sum_{i=1}^N \mathrm{Gamma}(a,\theta)$: a point mass at 0 plus a continuous positive part; variance follows the power law $\phi\mu^p$.

In [2]:
set.seed(1)
mu<-3; phi<-1.5; p<-1.6; Y<-tw_rng(mu,phi,p,3e5)
par(mfrow=c(1,2), mar=c(4,4,3,1))
hist(Y[Y>0], breaks=80, freq=FALSE, col="#cfe3f6", border="white", xlim=c(-0.5,18), main=sprintf("Tweedie(mu=%g,phi=%g,p=%g): spike at 0 + continuous",mu,phi,p), xlab="aggregate claim Y")
rect(-0.15,0,0.15,mean(Y==0), col=RED, border=NA); legend("topright", sprintf("P(Y=0)=%.3f",mean(Y==0)), fill=RED, bty="n", cex=.8)
mus<-seq(0.5,8,length=16)
plot(NULL, xlim=range(mus), ylim=c(0.3,80), log="xy", xlab="mean mu (log)", ylab="variance (log)", main="Power-variance law: Var = phi * mu^p")
for(a in list(c(1.3,BLUE),c(1.6,GREEN),c(1.9,RED))){ pp<-as.numeric(a[1]); v<-sapply(mus, function(m) var(tw_rng(m,1,pp,3e4))); points(mus,v,pch=19,col=a[2],cex=.7); lines(mus, mus^pp, col=a[2]) }
legend("topleft", c("p=1.3","p=1.6","p=1.9"), col=c(BLUE,GREEN,RED), lwd=2, bty="n", cex=.8)
par(mfrow=c(1,1))
cat(sprintf("mean %.3f (analytic %g); var %.3f (analytic %.3f); P(Y=0) %.4f (analytic %.4f)\n", mean(Y), mu, var(Y), phi*mu^p, mean(Y==0), tw_p0(mu,phi,p)))

mean 2.998 (analytic 3); var 8.718 (analytic 8.699); P(Y=0) 0.0754 (analytic 0.0753)
No description has been provided for this image

2. A random distribution — the Dirichlet process¶

Stick-breaking $\beta_k\sim\mathrm{Beta}(1,\alpha)$, $w_k=\beta_k\prod_{j<k}(1-\beta_j)$: a random discrete distribution. The Chinese Restaurant Process gives its cluster count, $\approx\alpha\log n$.

In [3]:
par(mfrow=c(1,3), mar=c(4,4,3,1))
for(a in list(c(0.5,BLUE),c(3,GREEN),c(15,RED))){ alpha<-as.numeric(a[1]); W<-gem(alpha,25,1)[1,]
  barplot(W, names.arg=1:25, col=a[2], ylim=c(0,0.75), main=sprintf("Stick-breaking weights, alpha=%g",alpha), xlab="stick k"); mtext(sprintf("E[w1]=1/(1+a)=%.2f",1/(1+alpha)), cex=.65) }
par(mfrow=c(1,1))
W<-gem(3,300,1e5); cat(sprintf("GEM(a=3): E[w1]=%.3f (analytic %.3f); E[w2]=%.3f (analytic %.3f); sum->%.4f\n", mean(W[,1]), 1/4, mean(W[,2]), 3/16, mean(rowSums(W))))

GEM(a=3): E[w1]=0.251 (analytic 0.250); E[w2]=0.188 (analytic 0.188); sum->1.0000
No description has been provided for this image
In [4]:
par(mfrow=c(1,2), mar=c(4,4,3,1))
plot(NULL, xlim=c(-6,6), ylim=c(0,0.6), xlab="atom location theta", ylab="weight", main="Four DP(alpha=3) draws: random discrete distributions")
cols<-c(BLUE,GREEN,RED,PURP); for(i in 1:4){ W<-gem(3,300,1)[1,]; loc<-rnorm(300,0,2); o<-order(loc); segments(loc[o],0,loc[o],W[o]+0.004*i, col=adjustcolor(cols[i],0.7)) }
ns<-unique(round(10^seq(0,3.3,length=30)))
plot(NULL, xlim=range(ns), ylim=c(0,40), log="x", xlab="n customers (log)", ylab="# clusters", main="Number of clusters ~ alpha*log n")
for(a in list(c(1,BLUE),c(5,RED))){ alpha<-as.numeric(a[1]); mm<-sapply(ns, function(n) mean(replicate(25, length(crp(alpha,n)$sizes)))); points(ns,mm,pch=19,col=a[2],cex=.6); lines(ns, sapply(ns, function(n) dp_ecl(alpha,n)), col=a[2], lwd=1.5) }
legend("topleft", c("CRP alpha=1","CRP alpha=5"), col=c(BLUE,RED), lwd=2, bty="n", cex=.8)
par(mfrow=c(1,1)); cat("Small alpha -> few dominant atoms & clusters; larger alpha -> more, smaller ones (clusters ~ alpha log n).\n")

Small alpha -> few dominant atoms & clusters; larger alpha -> more, smaller ones (clusters ~ alpha log n).
No description has been provided for this image

3. A random function — the Gaussian process¶

Any finite set of function values is jointly Gaussian with covariance $k(x,x')$. The kernel sets smoothness (RBF smooth, Matern-1/2 rough); conditioning on data gives closed-form GP regression with uncertainty.

In [5]:
Xg<-matrix(seq(0,10,length=200),ncol=1)
par(mfrow=c(1,3), mar=c(4,4,3,1))
for(a in list(list("RBF (smooth)",function(x,y) rbf_k(x,y,1.2,1),BLUE), list("Matern-3/2",function(x,y) mat32_k(x,y,1.2,1),GREEN), list("Matern-1/2 (rough)",function(x,y) mat12_k(x,y,1.2,1),RED))){
  F<-gp_prior(Xg, a[[2]], 5); plot(NULL, xlim=c(0,10), ylim=c(-3.2,3.2), xlab="x", ylab="f(x)", main=paste("GP prior draws -", a[[1]]))
  for(i in 1:5) lines(Xg[,1], F[i,], col=adjustcolor(a[[3]],0.7), lwd=1.2) }
par(mfrow=c(1,1)); cat("Each curve is one function from the GP prior; the kernel sets the smoothness (same variance, different regularity).\n")

Each curve is one function from the GP prior; the kernel sets the smoothness (same variance, different regularity).
No description has been provided for this image
In [6]:
truef<-function(x) sin(x)+0.3*x
Xtr<-matrix(c(0.5,1.5,3,4.5,6.5,8.5),ncol=1); ytr<-truef(Xtr[,1])+rnorm(6,0,0.15)
Xte<-matrix(seq(-0.5,10.5,length=250),ncol=1); ker<-function(x,y) rbf_k(x,y,1.3,1.2)
po<-gp_post(Xtr,ytr,Xte,ker,noise=0.15^2); sd<-sqrt(pmax(diag(po$cov),0))
plot(NULL, xlim=c(-0.5,10.5), ylim=c(-2,5), xlab="x", ylab="f(x)", main="Gaussian-process regression: mean + uncertainty band")
polygon(c(Xte[,1],rev(Xte[,1])), c(po$mean-2*sd, rev(po$mean+2*sd)), col="#cfe3f6", border=NA)
lines(Xte[,1], truef(Xte[,1]), col=GREY, lwd=1.3, lty=2); lines(Xte[,1], po$mean, col=BLUE, lwd=2.2); points(Xtr[,1], ytr, pch=19, col=RED, cex=1.3)
legend("topleft", c("posterior +/-2sd","posterior mean","true function","observations"), col=c("#cfe3f6",BLUE,GREY,RED), lwd=c(6,2,2,NA), pch=c(NA,NA,NA,19), lty=c(1,1,2,NA), bty="n", cex=.8)
cat("The band is tight at the data and widens in the gaps and beyond -- calibrated GP uncertainty.\n")

The band is tight at the data and widens in the gaps and beyond -- calibrated GP uncertainty.
No description has been provided for this image
In [7]:
par(mfrow=c(1,2), mar=c(3,3,3,1))
Xk<-matrix(seq(0,10,length=60),ncol=1); K<-rbf_k(Xk,Xk,1.2,1)
image(seq(0,10,length=60), seq(0,10,length=60), K, col=hcl.colors(24,"YlGnBu",rev=TRUE), xlab="x'", ylab="x", main="RBF kernel matrix K(x,x') = covariance of f")
g<-seq(0,5,length=30); G<-as.matrix(expand.grid(x1=g,x2=g)); Fs<-matrix(gp_prior(G, function(x,y) rbf_k(x,y,1.3,1), 1)[1,], 30, 30)
persp(g,g,Fs, theta=-55, phi=32, col="#a6d854", border=NA, shade=.4, xlab="x1", ylab="x2", zlab="f", main="One draw from a 2-D-input GP: a random surface")
par(mfrow=c(1,1)); cat("Left: the kernel matrix is the function's covariance. Right: over a 2-D input the GP draws a random landscape.\n")

Left: the kernel matrix is the function's covariance. Right: over a 2-D input the GP draws a random landscape.
No description has been provided for this image

4. Summary — and the end of the catalog¶

Every compound and nonparametric object's from-scratch R implementation is validated by its construction and moments: the Tweedie by mean/variance/mass-at-zero, the Dirichlet process by its stick-breaking weights and Chinese-Restaurant cluster counts, and the Gaussian process by its kernel covariance and closed-form posterior — a cross-language corroboration of the Python results, and a random 2-D surface to match.

This closes the Statistical Distributions catalog. Eight folders, built from scratch and validated against scipy and R: discrete laws; the Gaussian & exponential family; heavy tails & skewness; extreme value & survival; Bayesian priors & augmentation; multivariate; matrix-variate covariance priors; and the compound & nonparametric processes here — an arc from a coin flip to a random function.