Heavy Tails & Skewness¶
Statistical Distributions — modelling the extremes the Gaussian misses¶
Folder 3 of the Statistical Distributions catalog. The Gaussian family of Folder 2 has light tails — extremes are astronomically rare and every shape is symmetric. Real data in finance, insurance, hydrology and networks is not like that: crashes, floods, and viral events happen far more often than a Normal allows, and losses are rarely symmetric with gains. This notebook builds the heavy-tailed and skewed laws designed for exactly those features, each with a from-scratch PDF, CDF, RNG and moments, validated against scipy.stats.
What "heavy tail" and "skew" mean precisely¶
- A heavy (fat) tail decays slower than the Normal's $e^{-x^2/2}$. Power-law tails, $\Pr(X>x)\sim x^{-\alpha}$, are the extreme case: they decay so slowly that moments cease to exist — the Cauchy has no mean, the Lévy no finite anything, a $\mathrm{Pareto}(\alpha)$ only moments of order $<\alpha$. Excess kurtosis $>0$ is the mild diagnostic; an infinite tail index is the severe one.
- Skewness is asymmetry: a longer tail on one side. The Skew-Normal and Hansen skew-$t$ add a shape parameter that tilts an otherwise symmetric law.
The unifying idea: α-stable laws¶
The Central Limit Theorem says sums of finite-variance variables converge to the Normal. Drop the finite-variance assumption and sums converge instead to an α-stable law — a two-parameter family that contains the Normal ($\alpha=2$), the Cauchy ($\alpha=1$) and the Lévy ($\alpha=\tfrac12$) as special cases. Stable laws have no closed-form density in general, so we sample them by the Chambers–Mallows–Stuck algorithm and validate through the characteristic function — a highlight of this notebook.
Families covered¶
- Symmetric, heavier than Normal — Laplace, Logistic
- No moments & power laws — Cauchy, Pareto, Lévy
- Tunable shape — Student-$t$ (location-scale), Generalized Error Distribution
- Skewness — Skew-Normal
- The unifying family — α-stable (Chambers–Mallows–Stuck) and its special cases
- A finance capstone — Hansen's skew-$t$
- Tail diagnostics — survival curves & the Hill estimator
- Timing
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import heavytails as H
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["skew"], d["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):
'''Validation figure (PDF, CDF, RNG histogram) + moments table for one continuous law.
Histogram bars are normalised by the TOTAL sample size, so out-of-range extremes shrink
the visible mass correctly (essential for heavy tails).'''
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")
xs_s = np.sort(samples[(samples>=xs.min())&(samples<=xs.max())])
ax[1].plot(xs_s, np.searchsorted(np.sort(samples), xs_s)/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(H.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. Symmetric, heavier than Normal — Laplace & Logistic¶
Laplace(μ, b) — the double exponential, $f(x)=\frac1{2b}e^{-|x-\mu|/b}$: two exponential tails glued back-to-back. Its $e^{-|x|}$ tails are heavier than the Normal's $e^{-x^2}$ (excess kurtosis $3$), and the sharp peak makes it the density behind least-absolute-deviation (robust) regression and the LASSO prior. Sampled in closed form by inverse transform.
Logistic(μ, s) — density $f(x)=\dfrac{e^{-(x-\mu)/s}}{s\left(1+e^{-(x-\mu)/s}\right)^{2}}$; the bell-shaped law with slightly heavier tails than the Normal (excess kurtosis $6/5$) whose CDF is the logistic/sigmoid function — the backbone of logistic regression. Its inverse CDF is the logit $\ln\frac{U}{1-U}$, so it too samples in one line.
mu,b = 0.5, 1.2; xs = np.linspace(-8, 8, 500)
panel("Laplace(0.5, 1.2)", xs, H.laplace_pdf(xs,mu,b), H.laplace_cdf(xs,mu,b),
H.laplace_rng(mu,b,150000,rng), H.laplace_moments(mu,b),
S.laplace.pdf(xs,mu,b), S.laplace.cdf(xs,mu,b), sp_stats(S.laplace(mu,b)), hist_range=(-8,8))
analytical (formula) scipy empirical (RNG) mean 0.50 0.50 0.497 variance 2.88 2.88 2.877 skewness 0.00 0.00 0.010 excess kurtosis 3.00 3.00 3.052
mu,s = 0.5, 1.2; xs = np.linspace(-9, 9, 500)
panel("Logistic(0.5, 1.2)", xs, H.logistic_pdf(xs,mu,s), H.logistic_cdf(xs,mu,s),
H.logistic_rng(mu,s,150000,rng), H.logistic_moments(mu,s),
S.logistic.pdf(xs,mu,s), S.logistic.cdf(xs,mu,s), sp_stats(S.logistic(mu,s)), hist_range=(-9,9))
analytical (formula) scipy empirical (RNG) mean 0.500 0.500 0.496 variance 4.737 4.737 4.756 skewness 0.000 0.000 -0.002 excess kurtosis 1.200 1.200 1.305
2. No moments & power laws — Cauchy, Pareto, Lévy¶
Cauchy(x₀, γ) — the ratio of two Normals, and the Student-$t$ with one degree of freedom, $f(x)=\dfrac{1}{\pi\gamma\left[1+\left(\frac{x-x_0}{\gamma}\right)^{2}\right]}$. Its tails are so heavy that not even the mean exists: the sample average of Cauchy draws does not converge — it is itself Cauchy, no matter how many you take. The perfect cautionary tale for "just use the mean". Its quantile is $\tan$, so it samples in closed form.
Pareto(x_m, α) — the canonical power law, $f(x)\propto x^{-(\alpha+1)}$ on $x\ge x_m$; the "80/20" law of wealth, file sizes and city populations. The tail index α governs how many moments exist: the mean needs $\alpha>1$, the variance $\alpha>2$. Sampled by inverse transform $x_m/U^{1/\alpha}$.
Lévy(μ, c) — a one-sided heavy law on $(0,\infty)$ with an $x^{-3/2}$ tail so heavy that even the mean is infinite; it is the first-passage-time law of Brownian motion and the $\alpha=\tfrac12$ stable. Sampled elegantly as $c/Z^2$ with $Z\sim\text{Normal}$.
x0,g = 0.0, 1.0; xs = np.linspace(-12, 12, 500)
panel("Cauchy(0, 1)", xs, H.cauchy_pdf(xs,x0,g), H.cauchy_cdf(xs,x0,g),
H.cauchy_rng(x0,g,150000,rng), H.cauchy_moments(x0,g),
S.cauchy.pdf(xs,x0,g), S.cauchy.cdf(xs,x0,g), hist_range=(-12,12), show_moments=False,
emp_note="Moments: NONE exist — the mean, variance, skew and kurtosis of Cauchy are all undefined,\n"
"so no analytical/empirical table is shown. (A sample mean of Cauchy draws never settles down.)")
Moments: NONE exist — the mean, variance, skew and kurtosis of Cauchy are all undefined, so no analytical/empirical table is shown. (A sample mean of Cauchy draws never settles down.)
xm,al = 1.0, 3.0; xs = np.linspace(1.0, 12, 500)
panel("Pareto(x_m=1, α=3)", xs, H.pareto_pdf(xs,xm,al), H.pareto_cdf(xs,xm,al),
H.pareto_rng(xm,al,150000,rng), H.pareto_moments(xm,al),
S.pareto.pdf(xs,al,scale=xm), S.pareto.cdf(xs,al,scale=xm), sp_stats(S.pareto(al,scale=xm)),
hist_range=(1,12),
emp_note="\n(With α=3 the mean & variance exist but skewness needs α>3 and kurtosis α>4 — hence the NaNs.\n"
"The empirical values are dominated by rare huge draws and barely converge.)")
analytical (formula) scipy empirical (RNG) mean 1.50 1.50 1.498 variance 0.75 0.75 0.693 skewness NaN NaN 9.556 excess kurtosis NaN NaN 229.203 (With α=3 the mean & variance exist but skewness needs α>3 and kurtosis α>4 — hence the NaNs. The empirical values are dominated by rare huge draws and barely converge.)
mu,c = 0.0, 1.0; xs = np.linspace(0.02, 15, 500)
panel("Lévy(0, 1)", xs, H.levy_pdf(xs,mu,c), H.levy_cdf(xs,mu,c),
H.levy_rng(mu,c,150000,rng), H.levy_moments(mu,c),
S.levy.pdf(xs,mu,c), S.levy.cdf(xs,mu,c), hist_range=(0,15), show_moments=False,
emp_note="Moments: even the MEAN is infinite (the x^{-3/2} tail). Shown only as PDF/CDF/RNG.")
Moments: even the MEAN is infinite (the x^{-3/2} tail). Shown only as PDF/CDF/RNG.
3. Tunable shape — Student-$t$ & the Generalized Error Distribution¶
Two families with a knob that dials tail-heaviness continuously — the practical workhorses for robust and financial modelling.
Student-$t$ (μ, σ, ν) — location-scale $t$, $f(x)=\dfrac{\Gamma\!\left(\frac{\nu+1}{2}\right)}{\sigma\sqrt{\nu\pi}\,\Gamma\!\left(\frac{\nu}{2}\right)}\left(1+\frac{1}{\nu}\left(\frac{x-\mu}{\sigma}\right)^{2}\right)^{-(\nu+1)/2}$: Normal-like for large $\nu$, Cauchy at $\nu=1$, with the degrees of freedom $\nu$ setting the tail weight (excess kurtosis $6/(\nu-4)$). The default robust replacement for the Normal in regression and in GARCH innovations. We revisit it here as a model (with location and scale), sampling $\mu+\sigma\,Z/\sqrt{V/\nu}$ from our own Normal and Gamma engines.
Generalized Error Distribution (μ, α, β) — also called the exponential-power or Subbotin law, $f(x)\propto e^{-(|x-\mu|/\alpha)^\beta}$. The shape $\beta$ tunes kurtosis directly: $\beta=2$ is the Normal, $\beta=1$ the Laplace, $\beta\to\infty$ the Uniform, $\beta<1$ super-heavy. It is the other standard GARCH error law. Sampled as a signed $\mathrm{Gamma}^{1/\beta}$.
mu,sig,nu = 0.5, 1.2, 6; xs = np.linspace(-9, 9, 500)
panel("Student-t (μ=0.5, σ=1.2, ν=6)", xs, H.student_t_ls_pdf(xs,mu,sig,nu), H.student_t_ls_cdf(xs,mu,sig,nu),
H.student_t_ls_rng(mu,sig,nu,150000,rng), H.student_t_ls_moments(mu,sig,nu),
S.t.pdf(xs,nu,mu,sig), S.t.cdf(xs,nu,mu,sig), sp_stats(S.t(nu,mu,sig)), hist_range=(-9,9))
analytical (formula) scipy empirical (RNG) mean 0.50 0.50 0.501 variance 2.16 2.16 2.170 skewness 0.00 0.00 -0.015 excess kurtosis 3.00 3.00 3.020
mu,a,b = 0.3, 1.4, 1.5; xs = np.linspace(-8, 8, 500)
panel("GED (μ=0.3, α=1.4, β=1.5)", xs, H.ged_pdf(xs,mu,a,b), H.ged_cdf(xs,mu,a,b),
H.ged_rng(mu,a,b,150000,rng), H.ged_moments(mu,a,b),
S.gennorm.pdf(xs,b,mu,a), S.gennorm.cdf(xs,b,mu,a), sp_stats(S.gennorm(b,mu,a)), hist_range=(-8,8))
analytical (formula) scipy empirical (RNG) mean 0.300 0.300 0.297 variance 1.447 1.447 1.445 skewness 0.000 0.000 0.003 excess kurtosis 0.762 0.762 0.732
# Feature: the GED shape knob morphs Laplace (β=1) → Normal (β=2) → near-Uniform (β large)
xs = np.linspace(-4, 4, 500)
plt.figure(figsize=(10,4.2))
for b_,c,lab in [(1.0,RED,"β=1 (Laplace)"),(1.5,ORANGE,"β=1.5"),(2.0,GREEN,"β=2 (Normal)"),(5.0,BLUE,"β=5 (→ Uniform)")]:
ex = H.ged_moments(0,1,b_)["exkurt"]
plt.plot(xs, H.ged_pdf(xs,0,1,b_), lw=2, color=c, label=f"{lab} ex.kurt {ex:+.2f}")
plt.title("Generalized Error Distribution: the shape β tunes kurtosis"); plt.xlabel("x"); plt.legend(fontsize=9); plt.show()
4. Skewness — the Skew-Normal¶
Skew-Normal(ξ, ω, α) — Azzalini's one-parameter extension of the Normal, $f(x)=\frac2\omega\,\phi\!\big(\tfrac{x-\xi}\omega\big)\,\Phi\!\big(\alpha\tfrac{x-\xi}\omega\big)$: multiply a Normal density by a tilting CDF factor. The shape $\alpha$ controls asymmetry — $\alpha=0$ recovers the Normal, $\alpha>0$ skews right, $\alpha<0$ left — while keeping the tails Gaussian-light. Skewness saturates at about $\pm0.995$. We sample it by Azzalini's construction $\delta|Z_0|+\sqrt{1-\delta^2}\,Z_1$ with $\delta=\alpha/\sqrt{1+\alpha^2}$ — a hidden-truncation representation.
xi,om,al = 0.2, 1.3, 4.0; xs = np.linspace(-4, 6, 500)
panel("Skew-Normal (ξ=0.2, ω=1.3, α=4)", xs, H.skewnorm_pdf(xs,xi,om,al), H.skewnorm_cdf(xs,xi,om,al),
H.skewnorm_rng(xi,om,al,150000,rng), H.skewnorm_moments(xi,om,al),
S.skewnorm.pdf(xs,al,xi,om), S.skewnorm.cdf(xs,al,xi,om), sp_stats(S.skewnorm(al,xi,om)), hist_range=(-4,6))
analytical (formula) scipy empirical (RNG) mean 1.206 1.206 1.208 variance 0.677 0.677 0.682 skewness 0.784 0.784 0.791 excess kurtosis 0.633 0.633 0.659
# Feature: the skew knob α (negative → left, 0 → Normal, positive → right)
xs = np.linspace(-4, 4, 500)
plt.figure(figsize=(10,4.2))
for al_,c in zip([-6,-2,0,2,6], [PURP,BLUE,GREY,ORANGE,RED]):
sk = H.skewnorm_moments(0,1,al_)["skew"]
plt.plot(xs, H.skewnorm_pdf(xs,0,1,al_), lw=2, color=c, label=f"α={al_:+d} skew {sk:+.2f}")
plt.title("Skew-Normal: the shape α tilts the Gaussian"); plt.xlabel("x"); plt.legend(fontsize=9); plt.show()
5. The unifying family — α-stable laws (Chambers–Mallows–Stuck)¶
An α-stable law $S(\alpha,\beta)$ has two shape parameters: the stability $\alpha\in(0,2]$ (smaller = heavier tails; $\Pr(|X|>x)\sim x^{-\alpha}$ for $\alpha<2$) and the skewness $\beta\in[-1,1]$. It is the only possible limit of normalised sums of i.i.d. variables — the generalised Central Limit Theorem — and it swallows three laws we have already met: $$\alpha=2:\ \text{Normal}\qquad \alpha=1,\beta=0:\ \text{Cauchy}\qquad \alpha=\tfrac12,\beta=1:\ \text{Lévy}.$$ Except in those cases the density has no closed form, so we cannot write a PDF or CDF. Instead we sample directly by the Chambers–Mallows–Stuck (1976) algorithm — a deterministic transform of one Exponential and one Uniform draw — and validate two ways: the three special cases must reproduce their known densities, and the empirical characteristic function $\hat\varphi(t)=\frac1N\sum e^{itX_j}$ must match the theoretical $\varphi(t)=\exp\!\big(-|t|^\alpha[1-i\beta\,\text{sgn}(t)\tan\tfrac{\pi\alpha}2]\big)$.
# CMS special cases: ONE sampler reproduces Normal, Cauchy and Lévy
fig, ax = plt.subplots(1, 3, figsize=(14.5, 3.6))
g = H.stable_rng(2.0, 0.0, 200000, rng) # -> Normal(0, sqrt(2))
xx = np.linspace(-6,6,300); ax[0].hist(g, bins=np.linspace(-6,6,80), density=True, color="#cfe3f6", edgecolor="white", lw=.3)
ax[0].plot(xx, S.norm.pdf(xx,0,np.sqrt(2)), color=RED, lw=2); ax[0].set_title("α=2 → Normal(0, √2)"); ax[0].set_xlabel("x")
c = H.stable_rng(1.0, 0.0, 200000, rng) # -> Cauchy
xx = np.linspace(-12,12,300); w=np.full(len(c),1/(len(c)*(24/79)))
ax[1].hist(c, bins=np.linspace(-12,12,80), weights=w, color="#cfe3f6", edgecolor="white", lw=.3)
ax[1].plot(xx, S.cauchy.pdf(xx), color=RED, lw=2); ax[1].set_title("α=1, β=0 → Cauchy"); ax[1].set_xlabel("x")
lv = H.stable_rng(0.5, 1.0, 200000, rng) # -> Lévy(0,1) directly -- no re-centring needed
xx = np.linspace(0.05,15,300)
w=np.full(len(lv),1/(len(lv)*(15/79)))
ax[2].hist(lv, bins=np.linspace(0,15,80), weights=w, color="#cfe3f6", edgecolor="white", lw=.3)
ax[2].plot(xx, S.levy.pdf(xx), color=RED, lw=2); ax[2].set_title("α=½, β=1 → Lévy"); ax[2].set_xlabel("x"); ax[2].set_xlim(0,15)
plt.tight_layout(); plt.show()
print("A single Chambers–Mallows–Stuck sampler reproduces all three named laws (histograms) against their\n"
"closed-form densities (red) — including the Lévy case at α=½, β=1, which the CMS draws match\n"
"directly, with no re-centring: this parameterisation already puts it at Lévy(0,1).")
A single Chambers–Mallows–Stuck sampler reproduces all three named laws (histograms) against their closed-form densities (red) — including the Lévy case at α=½, β=1, which the CMS draws match directly, with no re-centring: this parameterisation already puts it at Lévy(0,1).
# Characteristic-function validation for a GENERIC stable law (no closed-form density exists)
alpha, beta = 1.5, 0.5
X = H.stable_rng(alpha, beta, 500000, rng)
t = np.linspace(-3, 3, 121)
ecf = np.array([np.mean(np.exp(1j*tt*X)) for tt in t])
tcf = H.stable_cf(t, alpha, beta)
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(t, tcf.real, color=RED, lw=2, label="theoretical Re φ(t)")
ax[0].plot(t, ecf.real, ".", color=BLUE, ms=4, label="empirical (CMS draws)")
ax[0].set_title(f"α-stable characteristic function, real part (α={alpha}, β={beta})"); ax[0].set_xlabel("t"); ax[0].legend(fontsize=8)
ax[1].plot(t, tcf.imag, color=RED, lw=2, label="theoretical Im φ(t)")
ax[1].plot(t, ecf.imag, ".", color=BLUE, ms=4, label="empirical (CMS draws)")
ax[1].set_title("imaginary part (the skew β tilts it)"); ax[1].set_xlabel("t"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"max |empirical − theoretical CF| over the grid = {np.max(np.abs(ecf-tcf)):.2e} (Monte-Carlo error at N=5e5).")
print("With no density to compare, the characteristic function is the rigorous check that the CMS sampler is correct.")
max |empirical − theoretical CF| over the grid = 1.38e-03 (Monte-Carlo error at N=5e5). With no density to compare, the characteristic function is the rigorous check that the CMS sampler is correct.
6. A finance capstone — Hansen's skew-$t$¶
Financial returns are simultaneously fat-tailed and asymmetric (crashes exceed rallies). Hansen's (1994) skew-$t$ is the standard density for this — the default innovation law in a large fraction of applied GARCH models. It is standardized (mean $0$, variance $1$ by construction) with two shape parameters: the degrees of freedom $\eta>2$ (tail weight) and the skew $\lambda\in(-1,1)$. Its density stitches two scaled-$t$ pieces at a kink:
$$g(z)=b\,c\Big(1+\tfrac1{\eta-2}\big(\tfrac{bz+a}{1\mp\lambda}\big)^2\Big)^{-(\eta+1)/2},\quad \begin{cases}1-\lambda & z<-a/b\\ 1+\lambda & z\ge -a/b\end{cases}$$
with constants $a,b,c$ fixed by the mean-0/variance-1 constraints. At $\lambda=0$ it collapses to the standardized Student-$t$. It has no elementary quantile, so we sample it by general inverse-transform sampling — build the CDF by integrating the density, then invert it — a technique that works for any distribution you can only write as a density. scipy has no skew-$t$, so we validate three ways: the density integrates to $1$, the $\lambda=0$ case matches our standardized $t$, and the RNG's empirical moments match the numerically-integrated ones.
eta,lam = 6.0, 0.4
xs = np.linspace(-6, 6, 500)
pdf = H.hansen_skewt_pdf(xs, eta, lam); cdf = H.hansen_skewt_cdf(xs, eta, lam)
samp = H.hansen_skewt_rng(eta, lam, 150000, rng)
m_an = H.hansen_skewt_moments(eta, lam)
# validation: integral, and lam=0 vs standardized-t
xg = np.linspace(-40,40,400001); integ = np.trapezoid(H.hansen_skewt_pdf(xg,eta,lam), xg)
std_t = H.student_t_ls_pdf(xs, 0.0, np.sqrt((eta-2)/eta), eta)
panel("Hansen skew-t (η=6, λ=0.4)", xs, pdf, cdf, samp, m_an, hist_range=(-6,6))
print(f"\nValidation: ∫ density dx = {integ:.6f} (should be 1)")
print(f"mean = {m_an['mean']:.4f}, variance = {m_an['var']:.4f} (standardized: 0 and 1 by construction)")
plt.figure(figsize=(9.5,4))
plt.plot(xs, pdf, color=BLUE, lw=2.2, label="Hansen skew-t (η=6, λ=0.4)")
plt.plot(xs, H.hansen_skewt_pdf(xs,eta,0.0), color=GREEN, lw=1.8, ls="--", label="λ=0 (symmetric)")
plt.plot(xs, std_t, color=RED, lw=1.2, ls=":", label="standardized Student-t(6)")
plt.title("Hansen skew-t: λ tilts it; at λ=0 it IS the standardized t"); plt.xlabel("z"); plt.legend(fontsize=9); plt.show()
print("The green (λ=0) and red (standardized-t) curves coincide — the skew-t nests the symmetric t exactly.")
analytical (formula) empirical (RNG) mean 0.000 0.001 variance 1.000 1.003 skewness 1.248 1.211 excess kurtosis 5.151 4.943 Validation: ∫ density dx = 1.000000 (should be 1) mean = 0.0000, variance = 1.0000 (standardized: 0 and 1 by construction)
The green (λ=0) and red (standardized-t) curves coincide — the skew-t nests the symmetric t exactly.
What "the mean does not exist" actually looks like¶
The catalog reports undefined moments rather than printing the large finite number a truncated sum would give, and that is easy to read as a technicality about integrals. It is not. It has a consequence you can watch.
The law of large numbers — the reason averaging works at all — needs the mean to exist. Where it does, the running average of a sample settles down onto it and stays there, and more data buys precision. Where it does not, the running average never settles: it drifts, jumps whenever a new extreme arrives, and is no closer to anything after a million draws than after a hundred.
Below, the same running average is computed for three laws with progressively heavier tails. The Normal has every moment. The Pareto(α=3) has a mean and a variance but no skewness. The Cauchy has no mean at all — and the difference is visible rather than notational.
fig,ax=plt.subplots(1,3,figsize=(15,4.2),sharex=True)
rng_d=np.random.default_rng(11); N=200_000; nn=np.arange(1,N+1)
cases=[("Normal(0,1)",rng_d.normal(0,1,N),0.0,BLUE,"mean = 0, and every moment exists"),
("Pareto(alpha=3)",(1-rng_d.uniform(size=N))**(-1/3.0),3/2,GREEN,"mean = 1.5; skewness needs alpha>3, so it has none"),
("Cauchy",rng_d.standard_cauchy(N),None,RED,"NO mean: the integral does not converge")]
for a,(nm,x,truth,c,note) in zip(ax,cases):
run=np.cumsum(x)/nn
a.plot(nn,run,color=c,lw=.9)
if truth is not None:
a.axhline(truth,color="#cbd5e0",lw=1.8,zorder=0,label=f"true mean = {truth:g}")
a.legend(fontsize=8,loc="lower right")
a.set_xscale("log"); a.set_xlabel("number of draws (log scale)"); a.set_title(nm,fontsize=10)
a.axhline(0,color="k",lw=.4)
if truth is not None: # limits from the curve, so nothing is clipped
lo,hi=float(run.min()),float(run.max()); pad=0.10*(hi-lo)
a.set_ylim(lo-pad,hi+pad)
ax[0].set_ylabel("running average of the sample")
plt.suptitle("The law of large numbers needs the mean to exist",y=1.02)
plt.tight_layout(); plt.show()
for nm,x,truth,c,note in cases:
run=np.cumsum(x)/nn
last=run[-1]; swing=np.max(run[nn>1000])-np.min(run[nn>1000])
print(f" {nm:16s} average after {N:,} draws {last:>9.3f} range over the last decade of draws {swing:>8.3f} ({note})")
print("\nThe first two settle: their running averages walk onto a level and stay there, and the visible wobble shrinks")
print("as the sample grows -- which is exactly what makes an estimate from a sample worth having.")
print("The Cauchy never does. Its running average is still wandering after 200,000 draws, and it will still be")
print("wandering after 200 million, because each new record-breaking draw is large enough to move the average on its")
print("own. There is no number it is converging TO. That is what 'the mean does not exist' means, and it is why the")
print("moment tables in this catalog report undefined instead of whatever a finite sum happened to return: reporting")
print("a number here would suggest an accuracy that more data cannot deliver.")
print("\nThe Pareto makes the subtler point. Its mean exists and converges perfectly well; it is the SKEWNESS that does")
print("not exist, since that needs alpha>3 and alpha is exactly 3. A distribution can be perfectly well behaved for")
print("the moments you happen to be looking at and undefined one step further out, which is why each moment carries")
print("its own existence condition rather than the distribution carrying one verdict.")
Normal(0,1) average after 200,000 draws -0.003 range over the last decade of draws 0.034 (mean = 0, and every moment exists) Pareto(alpha=3) average after 200,000 draws 1.502 range over the last decade of draws 0.042 (mean = 1.5; skewness needs alpha>3, so it has none) Cauchy average after 200,000 draws 2.749 range over the last decade of draws 3.756 (NO mean: the integral does not converge) The first two settle: their running averages walk onto a level and stay there, and the visible wobble shrinks as the sample grows -- which is exactly what makes an estimate from a sample worth having. The Cauchy never does. Its running average is still wandering after 200,000 draws, and it will still be wandering after 200 million, because each new record-breaking draw is large enough to move the average on its own. There is no number it is converging TO. That is what 'the mean does not exist' means, and it is why the moment tables in this catalog report undefined instead of whatever a finite sum happened to return: reporting a number here would suggest an accuracy that more data cannot deliver. The Pareto makes the subtler point. Its mean exists and converges perfectly well; it is the SKEWNESS that does not exist, since that needs alpha>3 and alpha is exactly 3. A distribution can be perfectly well behaved for the moments you happen to be looking at and undefined one step further out, which is why each moment carries its own existence condition rather than the distribution carrying one verdict.
7. Tail diagnostics — survival curves & the Hill estimator¶
How do you see a heavy tail? Two standard diagnostics:
- the survival function $S(x)=\Pr(X>x)$ on a log (or log-log) scale: light tails plunge, power-law tails become straight lines whose slope is the tail index $\alpha$;
- the Hill estimator, which reads $\alpha$ off the top order statistics: $\hat\alpha_k=\big[\frac1k\sum_{i=1}^{k}\ln\frac{X_{(n-i+1)}}{X_{(n-k)}}\big]^{-1}$. A flat Hill plot signals a genuine power law and estimates its index.
# Survival curves: the ordering of tail heaviness on a log-y scale
xs = np.linspace(0, 8, 400)
plt.figure(figsize=(10,4.4))
plt.semilogy(xs, 1-S.norm.cdf(xs), color=GREY, lw=2, label="Normal (e^{-x²})")
plt.semilogy(xs, 1-H.laplace_cdf(xs,0,1), color=GREEN, lw=2, label="Laplace (e^{-x})")
plt.semilogy(xs, 1-H.student_t_ls_cdf(xs,0,1,3), color=ORANGE, lw=2, label="Student-t(3) (power law)")
plt.semilogy(xs, 1-H.cauchy_cdf(xs,0,1), color=RED, lw=2, label="Cauchy (α=1)")
plt.semilogy(xs, np.where(xs>=1, 1-H.pareto_cdf(np.maximum(xs,1),1,1.5), np.nan), color=PURP, lw=2, label="Pareto(1.5)")
plt.ylim(1e-6, 1); plt.xlabel("x"); plt.ylabel("P(X > x) (log)")
plt.title("Survival curves: heavier tails sit higher and decay slower"); plt.legend(fontsize=8.5); plt.show()
# Hill estimator recovers the Pareto tail index from data
def hill(data, kmax, kmin=10):
x = np.sort(data)[::-1] # descending
lx = np.log(x)
ks = np.arange(kmin, kmax) # k below ~10 estimates a tail index from a handful
# of order statistics, which is noise, not an estimate
return ks, np.array([1.0/np.mean(lx[:k] - lx[k]) for k in ks])
lo_all, hi_all = np.inf, -np.inf
for al_true, c in [(1.5, BLUE), (2.5, RED)]:
data = H.pareto_rng(1.0, al_true, 200000, rng)
ks, ah = hill(data, 4000)
lo_all = min(lo_all, ah.min()); hi_all = max(hi_all, ah.max())
plt.plot(ks, ah, color=c, lw=1.5, label=f"data ~ Pareto(α={al_true})")
plt.axhline(al_true, color=c, ls="--", lw=1)
pad = 0.08*(hi_all-lo_all)
plt.ylim(max(0.0, lo_all-pad), hi_all+pad) # limits from the curves, so nothing is clipped
plt.xlabel("k (number of upper order statistics)"); plt.ylabel(r"Hill estimate $\hat\alpha_k$")
plt.title("Hill plot: the flat region recovers the true tail index α"); plt.legend(); plt.show()
print("Each Hill curve settles on its true α (dashed) for moderate k — the standard tail-index diagnostic.")
print("The plot starts at k = 10 deliberately. Below that the estimate is built from a handful of order")
print("statistics and swings wildly — at k = 2 it exceeds 10 for the α = 2.5 sample — which says nothing")
print("about the tail and everything about having almost no data. The instability still visible on the left")
print("is the real message: the flat middle is what you read, and choosing where it starts is the whole")
print("difficulty of using this estimator in practice.")
Each Hill curve settles on its true α (dashed) for moderate k — the standard tail-index diagnostic. The plot starts at k = 10 deliberately. Below that the estimate is built from a handful of order statistics and swings wildly — at k = 2 it exceeds 10 for the α = 2.5 sample — which says nothing about the tail and everything about having almost no data. The instability still visible on the left is the real message: the flat middle is what you read, and choosing where it starts is the whole difficulty of using this estimator in practice.
8. Timing: closed-form vs numerical samplers¶
For a fixed number of draws, the split here is between distributions with a closed-form quantile (Laplace, Logistic, Cauchy, Pareto — sampled by a single inverse-transform expression, extremely fast and at or above scipy speed) and those sampled by a construction or numerical method (GED via a Gamma draw, Skew-Normal via two normals, the CMS stable transform, and the Hansen skew-$t$ via numerical inverse-transform). The last of these pays for building a CDF grid once — the price of sampling a distribution you can only write as a density.
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 = [
("Laplace", lambda: H.laplace_rng(0.5,1.2,N,rng), lambda: S.laplace.rvs(0.5,1.2,size=N,random_state=rng), "inverse transform"),
("Logistic", lambda: H.logistic_rng(0.5,1.2,N,rng), lambda: S.logistic.rvs(0.5,1.2,size=N,random_state=rng), "inverse (logit)"),
("Cauchy", lambda: H.cauchy_rng(0,1,N,rng), lambda: S.cauchy.rvs(size=N,random_state=rng), "inverse (tan)"),
("Pareto", lambda: H.pareto_rng(1,3,N,rng), lambda: S.pareto.rvs(3,scale=1,size=N,random_state=rng), "inverse transform"),
("Student-t", lambda: H.student_t_ls_rng(0,1,6,N,rng), lambda: S.t.rvs(6,size=N,random_state=rng), "Normal/√(χ²/ν)"),
("GED", lambda: H.ged_rng(0,1.4,1.5,N,rng), lambda: S.gennorm.rvs(1.5,scale=1.4,size=N,random_state=rng), "signed Gamma^{1/β}"),
("Skew-Normal",lambda: H.skewnorm_rng(0,1,4,N,rng), lambda: S.skewnorm.rvs(4,size=N,random_state=rng), "Azzalini (2 normals)"),
("α-stable", lambda: H.stable_rng(1.5,0.5,N,rng), None, "Chambers–Mallows–Stuck"),
("Hansen t", lambda: H.hansen_skewt_rng(6,0.4,N,rng), None, "numerical inverse-CDF"),
]
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,5))
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 scipy equivalent)")
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: closed-form inverse transforms are fastest; numerical CDF sampling is the price of a density-only law")
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"
slowest = df["from-scratch (M/s)"].idxmin()
print("\nAgainst scipy's generic .rvs, the from-scratch generators ahead of it: " + ahead + ".")
print("Behind it: " + behind + " — still the same order of magnitude, with no per-draw Python loop anywhere.")
print("The stable and Hansen skew-t have NO scipy equivalent here. The CMS transform is a closed-form algebraic")
print("map and stays fast, while the numerical inverse-CDF is the slowest generator in the set (" + slowest + "),")
print("since it integrates a density grid before it can draw. That is the cost of sampling a law you can only")
print("express as a density — and it is still comfortably fast for practical N.")
from-scratch (M/s) scipy (M/s) ratio method
dist
Laplace 34.11 52.38 0.65 inverse transform
Logistic 55.85 80.61 0.69 inverse (logit)
Cauchy 43.33 39.96 1.08 inverse (tan)
Pareto 72.86 47.43 1.54 inverse transform
Student-t 8.42 34.24 0.25 Normal/√(χ²/ν)
GED 8.34 18.45 0.45 signed Gamma^{1/β}
Skew-Normal 18.36 25.36 0.72 Azzalini (2 normals)
α-stable 11.85 NaN NaN Chambers–Mallows–Stuck
Hansen t 7.80 NaN NaN numerical inverse-CDF
Against scipy's generic .rvs, the from-scratch generators ahead of it: Cauchy (1.08x), Pareto (1.54x). Behind it: Laplace (0.65x), Logistic (0.69x), Student-t (0.25x), GED (0.45x), Skew-Normal (0.72x) — still the same order of magnitude, with no per-draw Python loop anywhere. The stable and Hansen skew-t have NO scipy equivalent here. The CMS transform is a closed-form algebraic map and stays fast, while the numerical inverse-CDF is the slowest generator in the set (Hansen t), since it integrates a density grid before it can draw. That is the cost of sampling a law you can only express as a density — and it is still comfortably fast for practical N.
9. Summary¶
Nine heavy-tailed and skewed laws, each with a from-scratch PDF / CDF / RNG / moments validated against scipy.stats to machine precision (and, for the density-only laws, against the characteristic function and the mean-0/variance-1 constraints).
The themes. Heavier-than-Normal but still light: Laplace, Logistic. So heavy that moments vanish: Cauchy (no mean), Pareto (moments only up to the tail index), Lévy (infinite mean). A continuous tail knob: Student-$t$ (df $\nu$) and the GED (shape $\beta$). Asymmetry: the Skew-Normal. The grand unification: α-stable laws, sampled by Chambers–Mallows–Stuck and validated through the characteristic function, containing the Normal, Cauchy and Lévy as the cases $\alpha=2,1,\tfrac12$. And the finance capstone, Hansen's skew-$t$, fat-tailed and skewed, sampled by numerical inverse-transform. The tail-diagnostic toolkit — survival curves and the Hill estimator — shows how to detect and measure these tails in data.
Next in the catalog: extreme value, survival & reliability — the distributions of maxima and of time-to-failure (Gumbel, Fréchet, Weibull, generalized extreme value and Pareto, plus the actuarial GB2 / Burr / Dagum family).