Bayesian Priors & Data Augmentation¶
Statistical Distributions — the computational machinery of Bayesian inference¶
Folder 5 of the Statistical Distributions catalog. The earlier folders modelled data. This one collects the distributions that make Bayesian computation work — the priors we place on constrained and scale parameters, and the latent variables that turn an intractable likelihood into a sequence of easy conjugate updates. They are less famous than the Normal or the Gamma, but they are what runs inside every modern MCMC sampler.
Each law is built from scratch — PDF, CDF, RNG, moments — and validated against scipy.stats. Two of them have no ordinary closed-form density (the Pólya-Gamma) or a Bessel-function normaliser (the GIG); for those we validate through moments and the Laplace transform, exactly as a working probabilist would.
Three jobs these distributions do¶
- Constrain a parameter. A variance must be positive; a correlation lives in $(-1,1)$; a latent utility is truncated by the observed choice. Truncated and half distributions are the priors and the augmentation variables for these.
- Mix scales for shrinkage & robustness. Writing a heavy-tailed or sparse prior as a scale mixture of Normals makes it conjugate conditional on the scale. The Inverse-Gaussian and Generalized Inverse Gaussian are the mixing laws behind the Bayesian LASSO and the whole Normal-variance-mixture family.
- Linearise a hard likelihood. The Pólya-Gamma distribution is the single most important modern augmentation: conditional on a Pólya-Gamma latent, a logistic or negative-binomial likelihood becomes exactly Gaussian, so Gibbs sampling just works. The Dirichlet is the conjugate prior on the probability simplex.
Families covered¶
- Truncated & half priors — Truncated Normal, Half-Normal, Half-Cauchy, Half-$t$
- Positive mixing distributions — Inverse-Gaussian (Wald), Generalized Inverse Gaussian
- The augmentation centrepiece — Pólya-Gamma
- The simplex prior — Dirichlet
- Augmentation in action — Normal-variance mixtures (Bayesian LASSO) and the Pólya-Gamma logistic identity
- Timing
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import bayespriors as B
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 fmt_m(d): return [d["mean"], d["var"], d.get("skew"), d.get("exkurt")]
def panel(name, xs, pdf, cdf, samples, m_fs, sp_pdf=None, sp_cdf=None, m_sp=None,
hist_range=None, show_moments=True, emp_note=None):
fig, ax = plt.subplots(1, 3, figsize=(14.5, 3.5))
ax[0].plot(xs, pdf, color=BLUE, lw=2.2, label="from-scratch")
if sp_pdf is not None: ax[0].plot(xs, sp_pdf, "--", color=RED, lw=1.3, label="scipy")
ax[0].set_title(name+" · PDF"); ax[0].legend(fontsize=8); ax[0].set_xlabel("x"); ax[0].set_ylabel("f(x)")
ax[1].plot(xs, cdf, color=BLUE, lw=2.2, label="from-scratch")
if sp_cdf is not None: ax[1].plot(xs, sp_cdf, "--", color=RED, lw=1.3, label="scipy")
ss = np.sort(samples[(samples>=xs.min())&(samples<=xs.max())])
ax[1].plot(ss, np.searchsorted(np.sort(samples), ss)/len(samples), color=GREY, lw=1.0, alpha=.55, label="empirical CDF")
ax[1].set_title("CDF (line=scratch, dashes=scipy)"); ax[1].set_ylim(-.03,1.03); ax[1].set_xlabel("x"); ax[1].legend(fontsize=8)
lo, hi = (xs.min(), xs.max()) if hist_range is None else hist_range
edges = np.linspace(lo, hi, 91); bw = edges[1]-edges[0]
w = np.full(len(samples), 1.0/(len(samples)*bw))
ax[2].hist(samples, bins=edges, weights=w, color="#cfe3f6", edgecolor="white", lw=.3, label="RNG samples")
ax[2].plot(xs, pdf, color=BLUE, lw=1.8, label="analytical PDF")
ax[2].set_title("RNG vs analytical (N=%s)" % f"{len(samples):,}"); ax[2].legend(fontsize=8); ax[2].set_xlabel("x")
plt.tight_layout(); plt.show()
if show_moments:
cols = {"analytical (formula)": fmt_m(m_fs)}
if m_sp is not None: cols["scipy"] = list(m_sp)
cols["empirical (RNG)"] = fmt_m(B.emp_moments(samples))
print(pd.DataFrame(cols, index=["mean","variance","skewness","excess kurtosis"]).round(3).to_string())
if emp_note: print(emp_note)
def sp_stats(dist): return [float(x) for x in dist.stats("mvsk")]
print("helpers ready")
helpers ready
1. Truncated & half priors — constrained parameters and scale priors¶
Truncated Normal(μ, σ; a, b) — a Normal restricted to $[a,b]$ and renormalised. It is the data-augmentation variable of the probit model (Albert–Chib): the latent utility is a Normal truncated by the sign of the observed choice. Sampled exactly by inverse-CDF between the truncation points — draw $U$, return $\Phi^{-1}\!\big(\Phi(\alpha)+U(\Phi(\beta)-\Phi(\alpha))\big)$.
Half-Normal(σ) — the Normal folded at zero, the mildest positive prior on a scale/standard-deviation.
Half-Cauchy(σ) — Gelman's recommended weakly-informative prior for a hierarchical standard deviation: it concentrates near zero yet its heavy tail permits large values when the data demand them (it has no moments). The default in much applied Bayesian work.
Half-t(σ, ν) — the general folded-$t$ scale prior, interpolating between the Half-Cauchy ($\nu=1$) and the Half-Normal ($\nu\to\infty$).
mu,sig,a,b = 0.5, 1.0, -1.0, 3.0; xs = np.linspace(-1.5, 3.5, 500)
al,be = (a-mu)/sig, (b-mu)/sig
panel("Truncated Normal(0.5, 1; [-1,3])", xs, B.truncnorm_pdf(xs,mu,sig,a,b), B.truncnorm_cdf(xs,mu,sig,a,b),
B.truncnorm_rng(mu,sig,a,b,150000,rng), B.truncnorm_moments(mu,sig,a,b),
S.truncnorm.pdf(xs,al,be,mu,sig), S.truncnorm.cdf(xs,al,be,mu,sig), sp_stats(S.truncnorm(al,be,mu,sig)),
hist_range=(a,b))
analytical (formula) scipy empirical (RNG) mean 0.621 0.621 0.619 variance 0.729 0.729 0.727 skewness 0.277 0.277 0.279 excess kurtosis -0.522 -0.522 -0.510
sig = 1.5; xs = np.linspace(0, 6, 500)
panel("Half-Normal(σ=1.5)", xs, B.halfnorm_pdf(xs,sig), B.halfnorm_cdf(xs,sig),
B.halfnorm_rng(sig,150000,rng), B.halfnorm_moments(sig),
S.halfnorm.pdf(xs,scale=sig), S.halfnorm.cdf(xs,scale=sig), sp_stats(S.halfnorm(scale=sig)), hist_range=(0,6))
analytical (formula) scipy empirical (RNG) mean 1.197 1.197 1.198 variance 0.818 0.818 0.818 skewness 0.995 0.995 1.004 excess kurtosis 0.869 0.869 0.944
s = 1.2; xs = np.linspace(0, 14, 500)
panel("Half-Cauchy(σ=1.2)", xs, B.halfcauchy_pdf(xs,s), B.halfcauchy_cdf(xs,s),
B.halfcauchy_rng(s,150000,rng), B.halfcauchy_moments(s),
S.halfcauchy.pdf(xs,scale=s), S.halfcauchy.cdf(xs,scale=s), hist_range=(0,14), show_moments=False,
emp_note="Half-Cauchy has NO moments (the heavy tail) — that heavy tail is exactly why it makes a good\n"
"weakly-informative scale prior: strong pull toward 0, but large scales are never ruled out.")
Half-Cauchy has NO moments (the heavy tail) — that heavy tail is exactly why it makes a good weakly-informative scale prior: strong pull toward 0, but large scales are never ruled out.
# Feature: the three scale priors compared (why Half-Cauchy is the default)
xs = np.linspace(0, 8, 500)
plt.figure(figsize=(10,4.2))
plt.plot(xs, B.halfnorm_pdf(xs,1), color=BLUE, lw=2.2, label="Half-Normal(1) (light tail)")
plt.plot(xs, B.halfstudent_pdf(xs,1,3), color=GREEN, lw=2.2, label="Half-t(1, ν=3)")
plt.plot(xs, B.halfcauchy_pdf(xs,1), color=RED, lw=2.2, label="Half-Cauchy(1) (heavy tail)")
plt.title("Weakly-informative scale priors: all peak at 0, but the tail sets how large σ may go")
plt.xlabel("scale parameter"); plt.ylabel("prior density"); plt.legend(fontsize=9); plt.show()
print("The Half-Cauchy keeps meaningful mass far from 0, so the data can override the prior when a large scale is needed.")
The Half-Cauchy keeps meaningful mass far from 0, so the data can override the prior when a large scale is needed.
2. Positive mixing distributions — Inverse-Gaussian & GIG¶
Heavy-tailed and sparse priors are usually written as scale mixtures of Normals: $\beta\mid\tau\sim N(0,\tau)$ with $\tau$ drawn from a positive mixing law. Two mixing laws dominate.
Inverse-Gaussian / Wald(μ, λ) — the first-passage-time law of Brownian motion, $f(x)\propto x^{-3/2}\exp\!\big(-\tfrac{\lambda(x-\mu)^2}{2\mu^2 x}\big)$. It is the conditional distribution of the local scale in the Bayesian LASSO. We sample it exactly by the Michael–Schucany–Haas transformation — a normal, a quadratic, and a coin flip.
Generalized Inverse Gaussian GIG(p, a, b) — $f(x)\propto x^{p-1}\exp\!\big(-\tfrac12(ax+b/x)\big)$, normalised by a modified Bessel function $K_p(\sqrt{ab})$. It is the umbrella mixing law: Gamma ($b\to0$), Inverse-Gamma ($a\to0$) and Inverse-Gaussian ($p=-\tfrac12$) are all special cases, and it is the conditional posterior of the scale in the normal-mixture representation of the Bayesian LASSO and the generalized-hyperbolic family. Its quantile is not elementary, so we sample by numerical inverse-transform (specialised exact samplers, e.g. Devroye 2014, also exist).
mu,lam = 1.5, 2.0; xs = np.linspace(0.02, 8, 500)
panel("Inverse-Gaussian(μ=1.5, λ=2)", xs, B.invgauss_pdf(xs,mu,lam), B.invgauss_cdf(xs,mu,lam),
B.invgauss_rng(mu,lam,150000,rng), B.invgauss_moments(mu,lam),
S.invgauss.pdf(xs,mu/lam,scale=lam), S.invgauss.cdf(xs,mu/lam,scale=lam), sp_stats(S.invgauss(mu/lam,scale=lam)),
hist_range=(0,8))
analytical (formula) scipy empirical (RNG) mean 1.500 1.500 1.503 variance 1.688 1.688 1.684 skewness 2.598 2.598 2.585 excess kurtosis 11.250 11.250 11.010
p,a,b = 0.5, 1.5, 2.0; xs = np.linspace(0.03, 9, 500)
w, sc = np.sqrt(a*b), np.sqrt(b/a)
panel("GIG(p=0.5, a=1.5, b=2)", xs, B.gig_pdf(xs,p,a,b), B.gig_cdf(xs,p,a,b),
B.gig_rng(p,a,b,150000,rng), B.gig_moments(p,a,b),
S.geninvgauss.pdf(xs,p,w,scale=sc), S.geninvgauss.cdf(xs,p,w,scale=sc), sp_stats(S.geninvgauss(p,w,scale=sc)),
hist_range=(0,9))
analytical (formula) scipy empirical (RNG) mean 1.821 1.821 1.819 variance 1.659 1.659 1.661 skewness 1.830 1.830 1.829 excess kurtosis 5.312 5.312 5.311
# Feature: the GIG umbrella — its limits are the Gamma, Inverse-Gamma and Inverse-Gaussian
xs = np.linspace(0.02, 8, 500)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot(xs, B.gig_pdf(xs,2.0,2.0,1e-6), color=BLUE, lw=3, alpha=.5, label="GIG(p=2, a=2, b→0)")
ax[0].plot(xs, S.gamma.pdf(xs,2.0,scale=1.0), color=RED, lw=1.5, ls="--", label="= Gamma(2, 1)")
ax[0].plot(xs, B.gig_pdf(xs,-0.5,2.0,2.0), color=GREEN, lw=3, alpha=.5, label="GIG(p=−½, a=2, b=2)")
ax[0].plot(xs, B.invgauss_pdf(xs,1.0,2.0), color=ORANGE, lw=1.5, ls="--", label="= Inverse-Gaussian(1, 2)")
ax[0].set_title("GIG special cases: Gamma & Inverse-Gaussian limits"); ax[0].set_xlabel("x"); ax[0].legend(fontsize=8); ax[0].set_ylim(0,1.0)
# the p knob shifts the shape
for p_,c in zip([-1.0, 0.0, 1.5, 3.0], [RED,ORANGE,GREEN,BLUE]):
ax[1].plot(xs, B.gig_pdf(xs,p_,2.0,2.0), lw=2, color=c, label=f"p={p_}")
ax[1].set_title("GIG(p, a=2, b=2): the order p reshapes the density"); ax[1].set_xlabel("x"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("GIG(b→0)=Gamma, GIG(a→0)=Inverse-Gamma, GIG(p=−½)=Inverse-Gaussian — one family behind the scale-mixture zoo.")
GIG(b→0)=Gamma, GIG(a→0)=Inverse-Gamma, GIG(p=−½)=Inverse-Gaussian — one family behind the scale-mixture zoo.
3. The augmentation centrepiece — Pólya-Gamma¶
The logistic likelihood $\dfrac{(e^\psi)^{a}}{(1+e^\psi)^{b}}$ is the obstacle to easy Bayesian inference: it is not conjugate to a Gaussian prior on $\psi$, so plain Gibbs sampling stalls. Polson, Scott & Windle (2013) removed the obstacle with the Pólya-Gamma distribution $\mathrm{PG}(b,c)$ and the identity $$\frac{(e^\psi)^{a}}{(1+e^\psi)^{b}}=2^{-b}\,e^{\kappa\psi}\int_0^\infty e^{-\omega\psi^2/2}\,p_{\mathrm{PG}(b,0)}(\omega)\,d\omega,\qquad \kappa=a-\tfrac b2 .$$ Read right-to-left: conditional on the latent $\omega\sim\mathrm{PG}$, the logistic term is a Gaussian in $\psi$ — so a logistic (or negative-binomial) regression becomes a sequence of conjugate Normal updates. $\mathrm{PG}(b,c)$ has no elementary density, but it has a clean infinite Gamma-sum representation, $$\omega \stackrel{d}{=}\frac1{2\pi^2}\sum_{k=1}^\infty \frac{G_k}{(k-\tfrac12)^2+c^2/4\pi^2},\qquad G_k\stackrel{iid}\sim\mathrm{Gamma}(b,1),$$ which we truncate to sample. We validate three ways: the $b=1$ density via its alternating (Jacobi-theta) series, the mean $\tfrac{b}{2c}\tanh\tfrac c2$, and the Laplace transform $\mathbb E[e^{-t\omega}]=\cosh^b(\tfrac c2)/\cosh^b\!\big(\sqrt{(t+c^2/2)/2}\big)$.
# PG(1,c): histogram of the Gamma-sum sampler vs the analytical series density
fig, ax = plt.subplots(1, 2, figsize=(14, 4))
for c_,col in [(0.0,BLUE),(2.0,RED)]:
om = B.polyagamma_rng(1.0, c_, 200000, rng)
xx = np.linspace(1e-3, 2.2, 400)
edges = np.linspace(0, 2.2, 80); bw = edges[1]-edges[0]; w = np.full(len(om), 1/(len(om)*bw))
ax[0].hist(om, bins=edges, weights=w, color=col, alpha=.25)
ax[0].plot(xx, B.polyagamma1_pdf(xx, c_), color=col, lw=2.2, label=f"PG(1, c={c_})")
ax[0].axvline(B.polyagamma_moments(1,c_)["mean"], color=col, ls=":", lw=1.2)
ax[0].set_title("PG(1,c): Gamma-sum sampler vs analytical series density"); ax[0].set_xlabel("ω"); ax[0].legend(fontsize=8); ax[0].set_xlim(0,2.2)
# Laplace-transform validation (the rigorous check for a density-only law)
t = np.linspace(0, 8, 100); om = B.polyagamma_rng(2.0, 1.5, 300000, rng)
elt = np.array([np.mean(np.exp(-tt*om)) for tt in t]); tlt = B.polyagamma_laplace(t, 2.0, 1.5)
ax[1].plot(t, tlt, color=RED, lw=2, label="theoretical E[e^{−tω}]")
ax[1].plot(t, elt, ".", color=BLUE, ms=4, label="empirical (Gamma-sum draws)")
ax[1].set_title("PG(2, 1.5) Laplace transform"); ax[1].set_xlabel("t"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Series density integrates to {np.trapezoid(B.polyagamma1_pdf(np.linspace(1e-3,8,200001),2.0), np.linspace(1e-3,8,200001)):.4f};",
f"max |empirical − theoretical Laplace transform| = {np.max(np.abs(elt-tlt)):.2e}.")
Series density integrates to 1.0000; max |empirical − theoretical Laplace transform| = 4.23e-04.
4. The simplex prior — Dirichlet¶
The Dirichlet(α) is the multivariate generalisation of the Beta: a distribution over probability vectors $(x_1,\dots,x_K)$ with $x_i>0,\ \sum x_i=1$. It is the conjugate prior for the multinomial/categorical — the prior on the cell probabilities of a contingency table, the topic weights of an LDA document, or the mixture weights of a finite mixture. Its density is $f(x)\propto\prod_i x_i^{\alpha_i-1}$, and it has two beautiful properties we verify: each marginal is a Beta, $x_i\sim\mathrm{Beta}(\alpha_i,\alpha_0-\alpha_i)$ with $\alpha_0=\sum_j\alpha_j$, and it is sampled simply by normalising independent Gammas, $x_i=G_i/\sum_j G_j$ with $G_i\sim\mathrm{Gamma}(\alpha_i,1)$. The concentration $\alpha_0$ controls sparsity: small $\alpha$ pushes mass to the corners (near one-hot vectors), large $\alpha$ toward the centre (near-uniform).
# Dirichlet on the 3-simplex: the concentration knob (ternary scatter)
def to_2d(P): # project the simplex onto the plane
return P[:,1] + 0.5*P[:,2], (np.sqrt(3)/2)*P[:,2]
fig, ax = plt.subplots(1, 3, figsize=(14.5, 4.3))
for k,(al,title) in enumerate([(np.array([0.4,0.4,0.4]),"α=(0.4,0.4,0.4) sparse → corners"),
(np.array([1.,1.,1.]),"α=(1,1,1) uniform on simplex"),
(np.array([6.,6.,6.]),"α=(6,6,6) concentrated → centre")]):
P = B.dirichlet_rng(al, 4000, rng); x,y = to_2d(P)
ax[k].scatter(x, y, s=3, color=[BLUE,GREEN,RED][k], alpha=.35)
ax[k].plot([0,1,0.5,0],[0,0,np.sqrt(3)/2,0], color="k", lw=1)
ax[k].set_title(title, fontsize=9.5); ax[k].axis("off"); ax[k].set_aspect("equal")
plt.tight_layout(); plt.show()
# validate: a marginal is Beta, and the density matches scipy
al = np.array([2.,3.,5.]); P = B.dirichlet_rng(al, 200000, rng)
a_i,b_i = B.dirichlet_marginal_beta(al,0)
print(f"Marginal x1 ~ Beta({a_i:.0f},{b_i:.0f}): empirical mean {P[:,0].mean():.4f} vs analytical {a_i/(a_i+b_i):.4f},",
f"var {P[:,0].var():.5f} vs {a_i*b_i/((a_i+b_i)**2*(a_i+b_i+1)):.5f}")
pts = np.array([[0.2,0.3,0.5],[0.1,0.6,0.3]])
print("pdf vs scipy.dirichlet:", np.max(np.abs(B.dirichlet_pdf(pts,al) - np.array([S.dirichlet(al).pdf(p) for p in pts]))))
Marginal x1 ~ Beta(2,8): empirical mean 0.2002 vs analytical 0.2000, var 0.01462 vs 0.01455 pdf vs scipy.dirichlet: 0.0
5. Augmentation in action — scale mixtures and the logistic identity¶
Two demonstrations that these distributions really do collapse hard problems into easy ones.
(a) Bayesian LASSO as a Normal-variance mixture. The double-exponential (Laplace) shrinkage prior is exactly a Normal whose variance is Exponentially mixed: $$\beta\mid\tau\sim N(0,\tau),\quad \tau\sim\mathrm{Exp}(\lambda^2/2)\ \Longrightarrow\ \beta\sim\mathrm{Laplace}(0,1/\lambda).$$ So the sparse LASSO prior becomes a conditionally Gaussian update — the entire reason Gibbs sampling the Bayesian LASSO is easy.
(b) The Pólya-Gamma logistic identity. We check the Polson–Scott–Windle identity numerically: the logistic function $\sigma(\psi)=e^\psi/(1+e^\psi)$ equals its Gaussian-mixture representation $\tfrac12 e^{\psi/2}\,\mathbb E_{\omega\sim\mathrm{PG}(1,0)}[e^{-\omega\psi^2/2}]$ — a Gaussian in $\psi$ averaged over the latent $\omega$.
# (a) Laplace = Normal(0, τ) with τ ~ Exponential -> the Bayesian LASSO prior
lam = 1.0
tau = -np.log(rng.random(400000)) * (2/lam**2) # Exp(rate=λ²/2) == scale 2/λ²
beta = np.sqrt(tau) * B._std_normal(400000, rng) # β | τ ~ N(0, τ)
xx = np.linspace(-6, 6, 400)
plt.figure(figsize=(9.5,4))
edges = np.linspace(-6,6,120); bw=edges[1]-edges[0]; w=np.full(len(beta),1/(len(beta)*bw))
plt.hist(beta, bins=edges, weights=w, color="#cfe3f6", edgecolor="white", lw=.3, label="N(0,τ) with τ~Exp (mixture)")
plt.plot(xx, S.laplace.pdf(xx, 0, 1/lam), color=RED, lw=2.2, label="Laplace(0, 1/λ) density")
plt.title("Bayesian LASSO: a Normal with Exponentially-mixed variance IS the Laplace prior")
plt.xlabel("β"); plt.legend(fontsize=9); plt.show()
print("The sparse Laplace prior is a conditionally-Gaussian scale mixture — which is exactly what makes it Gibbs-sampleable.")
The sparse Laplace prior is a conditionally-Gaussian scale mixture — which is exactly what makes it Gibbs-sampleable.
# (b) The Pólya-Gamma logistic identity: σ(ψ) == ½ e^{ψ/2} E_PG[e^{−ω ψ²/2}]
psi = np.linspace(-6, 6, 240)
omega = B.polyagamma_rng(1.0, 0.0, 60000, rng) # a fixed pool of PG(1,0) draws
mix = 0.5 * np.exp(psi/2) * np.array([np.mean(np.exp(-omega * p**2 / 2)) for p in psi])
sigma = np.exp(psi) / (1 + np.exp(psi))
plt.figure(figsize=(9.5,4))
plt.plot(psi, sigma, color=BLUE, lw=3, alpha=.6, label="logistic σ(ψ) = e^ψ/(1+e^ψ)")
plt.plot(psi, mix, "--", color=RED, lw=1.8, label="½ e^{ψ/2} · E_PG[e^{−ω ψ²/2}] (Gaussian mixture)")
plt.title("Polson–Scott–Windle identity: the logistic IS a Gaussian scale-mixture over Pólya-Gamma")
plt.xlabel("ψ"); plt.ylabel("value"); plt.legend(fontsize=9); plt.show()
print(f"max |logistic − Pólya-Gamma mixture| over the grid = {np.max(np.abs(sigma-mix)):.3e} (Monte-Carlo error).")
print("Conditional on ω, the logistic term is Gaussian in ψ — so logistic regression becomes conjugate Normal updates.")
max |logistic − Pólya-Gamma mixture| over the grid = 2.048e-03 (Monte-Carlo error). Conditional on ω, the logistic term is Gaussian in ψ — so logistic regression becomes conjugate Normal updates.
6. Timing: exact transforms vs numerical samplers¶
The split is familiar. The inverse-CDF (Truncated Normal), folding (Half-Normal / Half-Cauchy / Half-$t$), the Michael–Schucany–Haas Inverse-Gaussian, the Gamma-sum Pólya-Gamma and the normalised-Gamma Dirichlet are all vectorised and fast. Only the GIG pays for a numerical inverse-transform (its quantile is not elementary), and the Pólya-Gamma's cost scales with the number of series terms retained.
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 = 200000
benches = [
("Trunc-Normal", lambda: B.truncnorm_rng(0.5,1,-1,3,N,rng), lambda: S.truncnorm.rvs(-1.5,2.5,0.5,1,size=N,random_state=rng), "inverse-CDF"),
("Half-Normal", lambda: B.halfnorm_rng(1.5,N,rng), lambda: S.halfnorm.rvs(scale=1.5,size=N,random_state=rng), "|Normal|"),
("Half-Cauchy", lambda: B.halfcauchy_rng(1.2,N,rng), lambda: S.halfcauchy.rvs(scale=1.2,size=N,random_state=rng), "|Cauchy|"),
("Inv-Gaussian", lambda: B.invgauss_rng(1.5,2,N,rng), lambda: S.invgauss.rvs(0.75,scale=2,size=N,random_state=rng), "Michael–Schucany–Haas"),
("GIG", lambda: B.gig_rng(0.5,1.5,2,N,rng), lambda: S.geninvgauss.rvs(0.5,np.sqrt(3),scale=np.sqrt(4/3),size=N,random_state=rng), "numerical inverse-CDF"),
("Pólya-Gamma", lambda: B.polyagamma_rng(1.0,1.5,N,rng), None, "Gamma-sum (200 terms)"),
("Dirichlet(3)", lambda: B.dirichlet_rng(np.array([2.,3.,5.]),N,rng), None, "normalised Gammas"),
]
rows = []
for name, fs, sp, kind in benches:
tf = best_time(fs); ts = best_time(sp) if sp else np.nan
rows.append((name, N/tf/1e6, (N/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())
order = df.index; y = np.arange(len(order)); h=.38
fig, ax = plt.subplots(figsize=(11,4.6))
ax.barh(y+h/2, df["from-scratch (M/s)"], h, color=BLUE, label="from-scratch")
ax.barh(y-h/2, df["scipy (M/s)"].fillna(0), h, color=GREY, alpha=.8, label="scipy (— = no direct scipy sampler)")
ax.set_yticks(y); ax.set_yticklabels(order); ax.set_xscale("log")
ax.set_xlabel("throughput — million samples / second (log)"); ax.invert_yaxis()
ax.set_title("RNG throughput: exact transforms are fast; only the GIG needs a numerical CDF")
ax.legend(); plt.tight_layout(); plt.show()
# 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"
no_ref = ", ".join(df.index[df["ratio"].isna()])
slowest = df["from-scratch (M/s)"].idxmin()
print("\nAgainst scipy, where a direct equivalent exists — ahead: " + ahead + ".")
print("Behind: " + behind + ". Folding and the MSH transform are exact but do more arithmetic per draw than")
print("scipy's compiled versions, so being exact does not by itself buy parity.")
print("No scipy sampler exists at all for: " + no_ref + ". The slowest generator in the whole set is "
+ "%s at %.2f M draws/s;" % (slowest, df["from-scratch (M/s)"].min()))
print("the Pólya-Gamma's Gamma-sum representation costs 200 Gamma draws per variate — that is the price of the")
print("augmentation that makes a logistic likelihood conjugate, paid once per latent per iteration.")
from-scratch (M/s) scipy (M/s) ratio method dist Trunc-Normal 33.87 6.62 5.11 inverse-CDF Half-Normal 56.45 84.72 0.67 |Normal| Half-Cauchy 57.30 59.57 0.96 |Cauchy| Inv-Gaussian 24.06 43.59 0.55 Michael–Schucany–Haas GIG 7.13 14.50 0.49 numerical inverse-CDF Pólya-Gamma 0.06 NaN NaN Gamma-sum (200 terms) Dirichlet(3) 4.24 NaN NaN normalised Gammas
Against scipy, where a direct equivalent exists — ahead: Trunc-Normal (5.11x). Behind: Half-Normal (0.67x), Half-Cauchy (0.96x), Inv-Gaussian (0.55x), GIG (0.49x). Folding and the MSH transform are exact but do more arithmetic per draw than scipy's compiled versions, so being exact does not by itself buy parity. No scipy sampler exists at all for: Pólya-Gamma, Dirichlet(3). The slowest generator in the whole set is Pólya-Gamma at 0.06 M draws/s; the Pólya-Gamma's Gamma-sum representation costs 200 Gamma draws per variate — that is the price of the augmentation that makes a logistic likelihood conjugate, paid once per latent per iteration.
7. Summary¶
Eight distributions that are the plumbing of Bayesian computation, each with a from-scratch PDF / CDF / RNG / moments validated against scipy.stats — and, for the two without an ordinary closed-form density, validated through the moments and Laplace transform instead.
The roles. Constrain: the Truncated Normal (the probit augmentation variable) and the Half-Normal / Half-Cauchy / Half-$t$ scale priors, with the Half-Cauchy the weakly-informative default precisely because its heavy tail never rules out a large scale. Mix: the Inverse-Gaussian (sampled exactly by Michael–Schucany–Haas) and the GIG umbrella, the positive mixing laws that turn heavy-tailed and sparse priors into conditionally-Gaussian updates. Linearise: the Pólya-Gamma, whose Gamma-sum representation makes logistic and negative-binomial likelihoods exactly Gaussian given the latent — the identity demonstrated numerically. Simplex: the Dirichlet, the conjugate prior on probability vectors, sampled by normalising Gammas with Beta marginals. The augmentation demos close the loop: the Bayesian LASSO as a Normal-variance mixture, and the logistic as a Pólya-Gamma Gaussian mixture.
Next in the catalog: multivariate distributions — the multivariate Normal and $t$, the Dirichlet's continuous cousins, spatial GMRF/CAR priors, and a cross-link to the copula toolkit.