import numpy as np, pandas as pd, matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa
import compound_np as C

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)
print("helpers ready")

helpers ready

mu, phi, p = 3.0, 1.5, 1.6
Y = C.tweedie_rng(mu, phi, p, 300000, rng)
mv = C.tweedie_mean_var(mu, phi, p); lam, a, th = C.tweedie_params(mu, phi, p)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].hist(Y[Y>0], bins=80, density=True, color="#cfe3f6", edgecolor="white", lw=.3)
ax[0].axvline(0, color=RED, lw=0); ax[0].bar([0], [np.mean(Y==0)], width=0.18, color=RED, label=f"P(Y=0)={np.mean(Y==0):.3f}")
ax[0].set_title(f"Tweedie(μ={mu}, φ={phi}, p={p}): spike at 0 + continuous part"); ax[0].set_xlabel("aggregate claim Y"); ax[0].legend(fontsize=8.5); ax[0].set_xlim(-0.5, 18)
# the power-variance law: Var = phi * mu^p, slope p on log-log
mus = np.linspace(0.5, 8, 20)
for p_,c in [(1.3,BLUE),(1.6,GREEN),(1.9,RED)]:
    v = [C.tweedie_rng(m, 1.0, p_, 40000, rng).var() for m in mus]
    ax[1].loglog(mus, v, "o", color=c, ms=4, label=f"p={p_} (slope {p_})")
    ax[1].loglog(mus, 1.0*mus**p_, "-", color=c, lw=1, alpha=.6)
ax[1].set_title("Power-variance law: Var(Y) = φ·μ^p"); ax[1].set_xlabel("mean μ (log)"); ax[1].set_ylabel("variance (log)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"mean {Y.mean():.3f} (analytic {mv['mean']}), variance {Y.var():.3f} (analytic {mv['var']:.3f}), P(Y=0) {np.mean(Y==0):.4f} (analytic {C.tweedie_p0(mu,phi,p):.4f})")
print(f"compound pieces: N~Poisson(λ={lam:.2f}), jumps~Gamma(shape={a:.2f}, scale={th:.2f}).")

mean 2.992 (analytic 3.0), variance 8.635 (analytic 8.699), P(Y=0) 0.0750 (analytic 0.0753)
compound pieces: N~Poisson(λ=2.59), jumps~Gamma(shape=0.67, scale=1.74).

# Stick-breaking weights: small alpha -> few dominant sticks, large alpha -> many small ones
fig, ax = plt.subplots(1, 3, figsize=(14.5, 4))
for k,alpha in enumerate([0.5, 3.0, 15.0]):
    W = C.gem_rng(alpha, 25, 1, rng)[0]
    ax[k].bar(np.arange(1, 26), W, color=[BLUE,GREEN,RED][k])
    ax[k].set_title(f"Stick-breaking weights, α={alpha}"); ax[k].set_xlabel("stick k"); ax[k].set_ylim(0, 0.75)
    ax[k].text(0.55, 0.8, f"E[w₁]=1/(1+α)={1/(1+alpha):.2f}", transform=ax[k].transAxes, fontsize=8.5)
plt.tight_layout(); plt.show()
W = C.gem_rng(3.0, 300, 100000, rng)
print(f"GEM(α=3): E[w₁]={W[:,0].mean():.3f} (analytic {1/4:.3f}); E[w₂]={W[:,1].mean():.3f} (analytic {3/16:.3f}); Σwₖ→{W.sum(1).mean():.4f}")

GEM(α=3): E[w₁]=0.250 (analytic 0.250); E[w₂]=0.188 (analytic 0.188); Σwₖ→1.0000

# The DP as a random discrete distribution: several draws (atoms from a Normal base measure)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
for i,c in zip(range(4), [BLUE,GREEN,RED,PURP]):
    locs, w = C.dp_random_measure(3.0, lambda K: rng.normal(0, 2, K), 300, rng)
    order = np.argsort(locs)
    ax[0].vlines(locs[order], 0, w[order]+0.002*i, color=c, alpha=.7, lw=1)
ax[0].set_title("Four DP(α=3) draws — each a random discrete distribution"); ax[0].set_xlabel("atom location θ"); ax[0].set_ylabel("weight"); ax[0].set_xlim(-6,6)
# Chinese Restaurant Process: number of clusters grows like alpha*log(n)
ns = np.unique(np.logspace(0, 3.3, 40).astype(int))
for alpha,c in [(1.0,BLUE),(5.0,RED)]:
    means = [np.mean([len(C.crp_rng(alpha, n, rng)[1]) for _ in range(30)]) for n in ns]
    ax[1].plot(ns, means, "o", color=c, ms=3, label=f"CRP simulated, α={alpha}")
    ax[1].plot(ns, [C.dp_expected_clusters(alpha, n) for n in ns], "-", color=c, lw=1.5, label=f"E[#clusters], α={alpha}")
ax[1].set_xscale("log"); ax[1].set_title("Number of clusters ≈ α·log n"); ax[1].set_xlabel("n customers (log)"); ax[1].set_ylabel("# occupied tables"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Small α → a few dominant atoms and few clusters; larger α → more, smaller atoms and more clusters (∝ α log n).")

Small α → a few dominant atoms and few clusters; larger α → more, smaller atoms and more clusters (∝ α log n).

# GP prior function draws: the kernel sets the smoothness
Xg = np.linspace(0, 10, 200)[:, None]
fig, ax = plt.subplots(1, 3, figsize=(14.5, 3.8))
for k,(name, ker, c) in enumerate([("RBF (smooth)", lambda a,b: C.rbf_kernel(a,b,1.2,1.0), BLUE),
                                   ("Matérn-3/2", lambda a,b: C.matern32_kernel(a,b,1.2,1.0), GREEN),
                                   ("Matérn-1/2 (rough)", lambda a,b: C.matern12_kernel(a,b,1.2,1.0), RED)]):
    F = C.gp_prior_sample(Xg, ker, 5, rng)
    for f in F: ax[k].plot(Xg.ravel(), f, color=c, alpha=.7, lw=1.2)
    ax[k].set_title(f"GP prior draws — {name}"); ax[k].set_xlabel("x"); ax[k].set_ylim(-3.2, 3.2)
plt.tight_layout(); plt.show()
print("Each curve is one function drawn from the GP prior. The kernel dictates the smoothness — same variance, different regularity.")

Each curve is one function drawn from the GP prior. The kernel dictates the smoothness — same variance, different regularity.

# GP regression: posterior mean, uncertainty band, and posterior function samples
def truef(x): return np.sin(x) + 0.3*x
Xtr = np.array([0.5, 1.5, 3.0, 4.5, 6.5, 8.5])[:, None]
ytr = truef(Xtr.ravel()) + rng.normal(0, 0.15, len(Xtr))
Xte = np.linspace(-0.5, 10.5, 250)[:, None]
ker = lambda a,b: C.rbf_kernel(a,b,1.3,1.2)
mean, cov = C.gp_posterior(Xtr, ytr, Xte, ker, noise=0.15**2)
sd = np.sqrt(np.clip(np.diag(cov), 0, None))
Fs = C.gp_posterior_sample(Xtr, ytr, Xte, ker, 5, rng, noise=0.15**2)
plt.figure(figsize=(10.5, 4.6))
plt.fill_between(Xte.ravel(), mean-2*sd, mean+2*sd, color="#cfe3f6", alpha=.8, label="posterior ±2 sd")
for f in Fs: plt.plot(Xte.ravel(), f, color=GREEN, lw=.8, alpha=.5)
plt.plot(Xte.ravel(), mean, color=BLUE, lw=2.2, label="posterior mean")
plt.plot(Xte.ravel(), truef(Xte.ravel()), "--", color=GREY, lw=1.3, label="true function")
plt.plot(Xtr.ravel(), ytr, "o", color=RED, ms=7, label="observations")
plt.title("Gaussian-process regression: mean, uncertainty, and posterior draws"); plt.xlabel("x"); plt.legend(fontsize=8.5); plt.show()
print("The band is tight at the data and widens in the gaps and beyond — calibrated uncertainty, the hallmark of GP regression.")

The band is tight at the data and widens in the gaps and beyond — calibrated uncertainty, the hallmark of GP regression.

# The kernel matrix (the covariance structure) and a 2-D-input GP surface (a random function of two variables)
fig = plt.figure(figsize=(13.5, 4.7))
ax1 = fig.add_subplot(1, 2, 1)
Xk = np.linspace(0, 10, 60)[:, None]; K = C.rbf_kernel(Xk, Xk, 1.2, 1.0)
im = ax1.imshow(K, cmap="viridis", extent=[0,10,10,0]); ax1.set_title("RBF kernel matrix K(x,x′) — the covariance of the function"); ax1.set_xlabel("x′"); ax1.set_ylabel("x"); plt.colorbar(im, ax=ax1, shrink=.8)
# 2-D input: sample a random surface f(x1,x2) from a GP
g = np.linspace(0, 5, 32); G1, G2 = np.meshgrid(g, g)
XX = np.column_stack([G1.ravel(), G2.ravel()])
Fs = C.gp_prior_sample(XX, lambda a,b: C.rbf_kernel(a,b,1.3,1.0), 1, rng)[0].reshape(32, 32)
ax2 = fig.add_subplot(1, 2, 2, projection="3d")
ax2.plot_surface(G1, G2, Fs, cmap="terrain", linewidth=0, antialiased=True)
ax2.set_title("One draw from a 2-D-input GP: a random surface f(x₁,x₂)", fontsize=10); ax2.set_xlabel("x₁"); ax2.set_ylabel("x₂"); ax2.view_init(35, -55)
plt.tight_layout(); plt.show()
print("Left: the kernel matrix IS the function's covariance (nearby inputs strongly correlated). Right: over a 2-D input the GP prior draws a random landscape.")

Left: the kernel matrix IS the function's covariance (nearby inputs strongly correlated). Right: over a 2-D input the GP prior draws a random landscape.

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
Xg = np.linspace(0,10,200)[:,None]
rows = [
  ("Tweedie",        best_time(lambda: C.tweedie_rng(3,1.5,1.6,100000,rng)), 100000, "compound Poisson-Gamma"),
  ("Stick-breaking", best_time(lambda: C.gem_rng(3.0,200,100000,rng)),       100000, "closed-form Beta(1,α)"),
  ("CRP (n=500)",    best_time(lambda: C.crp_rng(3.0,500,rng)),              1,       "sequential seating"),
  ("GP prior (m=200)", best_time(lambda: C.gp_prior_sample(Xg, lambda a,b: C.rbf_kernel(a,b,1.2,1.0), 10, rng)), 10, "Cholesky O(m³)"),
]
df = pd.DataFrame([(n, (cnt/t if cnt>1 else 1/t), k) for n,t,cnt,k in rows],
                  columns=["object", "draws/sec (or 1/sec)", "method"]).set_index("object")
print(df.round(1).to_string())
print("\nTweedie and stick-breaking are cheap and vectorised; the CRP is sequential; the GP is Cholesky-bound (O(m³)) —")
print("the reason scalable GP methods (inducing points, sparse kernels) exist.")

                  draws/sec (or 1/sec)                  method
object                                                        
Tweedie                      3878449.4  compound Poisson-Gamma
Stick-breaking                254364.4   closed-form Beta(1,α)
CRP (n=500)                      273.9      sequential seating
GP prior (m=200)                8378.0          Cholesky O(m³)

Tweedie and stick-breaking are cheap and vectorised; the CRP is sequential; the GP is Cholesky-bound (O(m³)) —
the reason scalable GP methods (inducing points, sparse kernels) exist.

Compound & Nonparametric Distributions¶

Statistical Distributions — when the parameter becomes a distribution or a function¶

The three ideas¶

Families covered¶

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

2. A random distribution — the Dirichlet process¶

3. A random function — the Gaussian process¶

4. Timing¶

5. Summary — and the end of the catalog¶