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) — 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. 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$: 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.
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):
x = np.sort(data)[::-1] # descending
lx = np.log(x)
ks = np.arange(2, kmax)
return ks, np.array([1.0/np.mean(lx[:k] - lx[k]) for k in ks])
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)
plt.plot(ks, ah, color=c, lw=1.5, label=f"data ~ Pareto(α={al_true})")
plt.axhline(al_true, color=c, ls="--", lw=1)
plt.ylim(0, 3.2); 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.")
Each Hill curve settles on its true α (dashed) for moderate k — the standard tail-index diagnostic.
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 47.81 62.33 0.77 inverse transform
Logistic 84.68 96.74 0.88 inverse (logit)
Cauchy 54.83 41.06 1.34 inverse (tan)
Pareto 73.02 59.47 1.23 inverse transform
Student-t 10.06 34.17 0.29 Normal/√(χ²/ν)
GED 9.89 17.18 0.58 signed Gamma^{1/β}
Skew-Normal 22.75 32.11 0.71 Azzalini (2 normals)
α-stable 14.34 NaN NaN Chambers–Mallows–Stuck
Hansen t 6.84 NaN NaN numerical inverse-CDF
Against scipy's generic .rvs, the from-scratch generators ahead of it: Cauchy (1.34x), Pareto (1.23x). Behind it: Laplace (0.77x), Logistic (0.88x), Student-t (0.29x), GED (0.58x), Skew-Normal (0.71x) — 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).