Extreme Value, Survival & Reliability

Statistical Distributions — the laws of maxima and of time-to-failure

Folder 4 of the Statistical Distributions catalog. Two questions drive this notebook, and they turn out to share a mathematics:

• How large is the worst case? The maximum of many observations — the biggest flood in a century, the largest insurance claim, the record heat — does not follow the distribution of a single observation. Extreme value theory says block maxima converge to one three-parameter family, the Generalized Extreme Value (GEV) law, and exceedances over a high threshold to the Generalized Pareto (GPD).
• When does it fail? The time until a component breaks, a patient relapses, or a life ends is a survival problem, described not by the density but by the hazard function $h(x)=f(x)/S(x)$ — the instantaneous failure rate given survival so far.

We build each law from scratch — PDF, CDF, RNG, moments, and (for the survival laws) the hazard — and validate against scipy.stats. The actuarial loss/income family (GB2, Burr, Dagum) rounds out the catalogue of heavy positive-support distributions.

The two unifying objects

• The GEV / GPD limit laws. Just as the Central Limit Theorem makes the Normal universal for sums, the Fisher–Tippett–Gnedenko theorem makes the GEV universal for maxima: whatever the parent distribution, its normalised block maxima converge to a GEV with shape $\xi$ — $\xi<0$ bounded (Weibull-type), $\xi=0$ light (Gumbel), $\xi>0$ heavy (Fréchet). The Pickands–Balkema–de Haan theorem is its threshold counterpart, giving the GPD.
• The hazard function. $h(x)$ and the survival $S(x)=1-F(x)$ re-express any positive law in the language of failure: a decreasing hazard is infant mortality, constant is memoryless (exponential), increasing is wear-out, and their sum is the classic bathtub curve.

Families covered

1. Extreme value building blocks — Gumbel, Weibull, Fréchet
2. The unifying maxima law — GEV, and block-maxima convergence
3. Peaks over threshold — the Generalized Pareto & mean-excess
4. The hazard function & survival — Weibull bathtub, Log-logistic, Gompertz–Makeham
5. Actuarial loss & income — GB2 and its Burr / Dagum special cases
6. Timing

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import extreme_survival as E
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(E.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. Extreme value building blocks — Gumbel, Weibull, Fréchet

The three classical types of maxima, distinguished by their tail:

Gumbel(μ, β) — the Type I law, for maxima of light-tailed parents (Normal, Exponential, Gamma). Right-skewed with a fixed skewness $\approx1.14$; its double-exponential CDF $e^{-e^{-z}}$ is the workhorse of flood and temperature extremes. Sampled by $\mu-\beta\ln(-\ln U)$.

Weibull(k, λ) — the Type III law (for maxima of bounded parents, in reversed form) and, far more commonly, the central distribution of reliability engineering. The shape $k$ sets the hazard's direction: $k<1$ decreasing (infant mortality), $k=1$ constant (the Exponential), $k>1$ increasing (wear-out). Sampled by $\lambda(-\ln U)^{1/k}$.

Fréchet(α, s, m) — the Type II law, for maxima of heavy-tailed (power-law) parents. Its own tail is heavy: only moments of order $<\alpha$ exist. Sampled by $m+s(-\ln U)^{-1/\alpha}$.

mu,beta = 1.0, 2.0; xs = np.linspace(-6, 12, 500)
panel("Gumbel(1, 2)", xs, E.gumbel_pdf(xs,mu,beta), E.gumbel_cdf(xs,mu,beta),
      E.gumbel_rng(mu,beta,150000,rng), E.gumbel_moments(mu,beta),
      S.gumbel_r.pdf(xs,mu,beta), S.gumbel_r.cdf(xs,mu,beta), sp_stats(S.gumbel_r(mu,beta)), hist_range=(-6,12))

                 analytical (formula)  scipy  empirical (RNG)
mean                            2.154  2.154            2.150
variance                        6.580  6.580            6.571
skewness                        1.140  1.140            1.153
excess kurtosis                 2.400  2.400            2.480

k,lam = 1.8, 3.0; xs = np.linspace(0, 10, 500)
panel("Weibull(k=1.8, λ=3)", xs, E.weibull_pdf(xs,k,lam), E.weibull_cdf(xs,k,lam),
      E.weibull_rng(k,lam,150000,rng), E.weibull_moments(k,lam),
      S.weibull_min.pdf(xs,k,scale=lam), S.weibull_min.cdf(xs,k,scale=lam), sp_stats(S.weibull_min(k,scale=lam)), hist_range=(0,10))

                 analytical (formula)  scipy  empirical (RNG)
mean                            2.668  2.668            2.671
variance                        2.352  2.352            2.357
skewness                        0.779  0.779            0.786
excess kurtosis                 0.557  0.557            0.609

al,s,m = 2.5, 2.0, 1.0; xs = np.linspace(1.0, 16, 500)
panel("Fréchet(α=2.5, s=2, m=1)", xs, E.frechet_pdf(xs,al,s,m), E.frechet_cdf(xs,al,s,m),
      E.frechet_rng(al,s,m,150000,rng), E.frechet_moments(al,s,m),
      S.invweibull.pdf(xs,al,loc=m,scale=s), S.invweibull.cdf(xs,al,loc=m,scale=s), sp_stats(S.invweibull(al,loc=m,scale=s)),
      hist_range=(1,16), emp_note="(α=2.5: mean & variance exist, but skewness needs α>3 and kurtosis α>4 — hence the NaNs. "
                "Note scipy still prints finite numbers there: its generic moment machinery does not check the\n existence condition, so those two values are artefacts, not a disagreement.)")

                 analytical (formula)   scipy  empirical (RNG)
mean                            3.978   3.978            3.992
variance                        9.493   9.493            8.844
skewness                          NaN  -5.396           17.782
excess kurtosis                   NaN  10.727          758.604
(α=2.5: mean & variance exist, but skewness needs α>3 and kurtosis α>4 — hence the NaNs. Note scipy still prints finite numbers there: its generic moment machinery does not check the
 existence condition, so those two values are artefacts, not a disagreement.)

2. The unifying maxima law — Generalized Extreme Value

The GEV(μ, σ, ξ) subsumes all three types through a single shape parameter ξ: $$F(x)=\exp\!\Big[-\big(1+\xi\tfrac{x-\mu}\sigma\big)^{-1/\xi}\Big],\qquad \xi\to0:\ \exp(-e^{-z}).$$

• $\xi<0$ — Weibull-type: a bounded upper tail (there is a hard maximum).
• $\xi=0$ — Gumbel-type: a light, unbounded exponential tail.
• $\xi>0$ — Fréchet-type: a heavy, power-law tail (moments fail for $\xi\ge1/k$).

This is the Fisher–Tippett–Gnedenko theorem in one formula: normalised block maxima of essentially any parent converge to a GEV. We verify that below by taking maxima of samples and watching the GEV appear. (Note scipy.genextreme uses shape $c=-\xi$.)

mu,sig,xi = 1.0, 2.0, 0.3; xs = np.linspace(-4, 16, 500)
panel("GEV(μ=1, σ=2, ξ=0.3)  Fréchet-type", xs, E.gev_pdf(xs,mu,sig,xi), E.gev_cdf(xs,mu,sig,xi),
      E.gev_rng(mu,sig,xi,150000,rng), E.gev_moments(mu,sig,xi),
      S.genextreme.pdf(xs,-xi,loc=mu,scale=sig), S.genextreme.cdf(xs,-xi,loc=mu,scale=sig), sp_stats(S.genextreme(-xi,loc=mu,scale=sig)),
      hist_range=(-4,16))

                 analytical (formula)   scipy  empirical (RNG)
mean                            2.987   2.987            3.002
variance                       23.698  23.698           23.810
skewness                       13.484  13.484            7.277
excess kurtosis                   NaN     NaN          155.824

# Feature: one family, three tail types
xs = np.linspace(-6, 12, 500)
plt.figure(figsize=(10,4.2))
for xi_,c,lab in [(-0.4,BLUE,"ξ=−0.4  (bounded, Weibull-type)"),(0.0,GREY,"ξ=0  (Gumbel-type)"),(0.4,RED,"ξ=+0.4  (heavy, Fréchet-type)")]:
    plt.plot(xs, E.gev_pdf(xs,1,2,xi_), lw=2.2, color=c, label=lab)
plt.axvline(1 - 2/(-0.4), color=BLUE, ls=":", lw=1)   # right endpoint for xi<0
plt.title("GEV(1, 2, ξ): the shape ξ sets the tail type"); plt.xlabel("x"); plt.legend(fontsize=9); plt.ylim(0,0.25); plt.show()

# Fisher–Tippett–Gnedenko in action: block maxima converge to a GEV
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
# (a) maxima of light-tailed Exponentials -> Gumbel (ξ=0)
blocks = -np.log(rng.random((60000, 80)))              # 60000 blocks of 80 Exp(1)
Mn = blocks.max(axis=1)
bn = np.log(80); Zn = Mn - bn                           # standardising: a_n=1, b_n=ln n
xx = np.linspace(-2, 8, 300)
ax[0].hist(Zn, bins=80, density=True, color="#cfe3f6", edgecolor="white", lw=.3, label="max of 80 Exponentials")
ax[0].plot(xx, E.gumbel_pdf(xx,0,1), color=RED, lw=2, label="Gumbel(0,1) limit")
ax[0].set_title("Light-tailed parent → Gumbel (ξ=0)"); ax[0].set_xlabel("standardised max"); ax[0].legend(fontsize=8)
# (b) maxima of heavy-tailed Pareto -> Fréchet (ξ>0)
al = 2.0; par = (1 - rng.random((60000, 80)))**(-1/al)  # Pareto(α=2), xm=1
Mn = par.max(axis=1); Zn = Mn / (80**(1/al))            # normalising a_n = n^{1/α}
xx = np.linspace(0.05, 8, 300)
w = np.full(len(Zn), 1/(len(Zn)*(8/80)))
ax[1].hist(Zn, bins=np.linspace(0,8,81), weights=w, color="#cfe3f6", edgecolor="white", lw=.3, label="max of 80 Pareto(2)")
ax[1].plot(xx, E.frechet_pdf(xx,al,1,0), color=RED, lw=2, label="Fréchet(α=2) limit")
ax[1].set_title("Heavy-tailed parent → Fréchet (ξ>0)"); ax[1].set_xlabel("normalised max"); ax[1].legend(fontsize=8); ax[1].set_xlim(0,8)
plt.tight_layout(); plt.show()
print("Whatever the parent, the normalised block maximum converges to a GEV — Gumbel for light tails, Fréchet for heavy.")

Whatever the parent, the normalised block maximum converges to a GEV — Gumbel for light tails, Fréchet for heavy.

3. Peaks over threshold — the Generalized Pareto

Block maxima throw away all but one point per block. The peaks-over-threshold (POT) approach instead keeps every exceedance above a high threshold $u$ — and the Pickands–Balkema–de Haan theorem says those exceedances converge to a Generalized Pareto(σ, ξ): $$F(y)=1-\big(1+\xi\tfrac{y}\sigma\big)^{-1/\xi},\qquad \xi\to0:\ 1-e^{-y/\sigma}.$$ It is the GEV's threshold twin, with the same shape $\xi$ classifying the tail, and it is the standard tool for Value-at-Risk and reinsurance. A key diagnostic is the mean-excess (mean residual life) function $e(u)=\mathbb E[X-u\mid X>u]$, which is linear in $u$ for a GPD, with slope $\xi/(1-\xi)$.

sig,xi = 1.5, 0.3; xs = np.linspace(0, 14, 500)
panel("GPD(σ=1.5, ξ=0.3)", xs, E.gpd_pdf(xs,0,sig,xi), E.gpd_cdf(xs,0,sig,xi),
      E.gpd_rng(0,sig,xi,150000,rng), E.gpd_moments(0,sig,xi),
      S.genpareto.pdf(xs,xi,scale=sig), S.genpareto.cdf(xs,xi,scale=sig), sp_stats(S.genpareto(xi,scale=sig)), hist_range=(0,14))

                 analytical (formula)   scipy  empirical (RNG)
mean                            2.143   2.143            2.140
variance                       11.480  11.480           10.891
skewness                       16.444  16.444            7.576
excess kurtosis                   NaN     NaN          145.881

# Mean-excess plot: linear for a GPD, and the POT limit theorem
def mean_excess(data, us):
    return np.array([np.mean(data[data>u]-u) if np.any(data>u) else np.nan for u in us])
data = E.gpd_rng(0, 1.5, 0.3, 200000, rng)
us = np.linspace(0, 8, 40)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot(us, mean_excess(data, us), "o", color=BLUE, ms=4, label="empirical e(u)")
ax[0].plot(us, (1.5 + 0.3*us)/(1-0.3), color=RED, lw=2, label="theoretical (σ+ξu)/(1−ξ)")
ax[0].set_title("Mean-excess is linear for a GPD"); ax[0].set_xlabel("threshold u"); ax[0].set_ylabel("E[X−u | X>u]"); ax[0].legend(fontsize=8)
# POT: exceedances of a parent (here Student-t) over a high threshold -> GPD
tdata = S.t.rvs(4, size=400000, random_state=rng)
u = np.quantile(tdata, 0.95); exc = tdata[tdata>u] - u
xx = np.linspace(0, exc.max(), 300)
sig_hat, xi_hat = 1.6, 0.25   # illustrative GPD fit params for t(4) tail
ax[1].hist(exc, bins=60, density=True, color="#cfe3f6", edgecolor="white", lw=.3, label="exceedances over 95th pct of t(4)")
ax[1].plot(xx, E.gpd_pdf(xx,0,sig_hat,xi_hat), color=RED, lw=2, label=f"GPD(σ≈{sig_hat}, ξ≈{xi_hat})")
ax[1].set_title("Peaks over threshold → GPD"); ax[1].set_xlabel("excess over threshold"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Left: the GPD's mean-excess rises linearly. Right: exceedances of a heavy-tailed parent over a high threshold are GPD-shaped.")

Left: the GPD's mean-excess rises linearly. Right: exceedances of a heavy-tailed parent over a high threshold are GPD-shaped.

4. The hazard function & survival — Weibull, Log-logistic, Gompertz–Makeham

Survival analysis re-expresses a positive distribution through the survival function $S(x)=\Pr(X>x)=1-F(x)$ and the hazard $h(x)=f(x)/S(x)$ — the failure rate given survival to $x$. The shape of $h$ tells the story:

Weibull — the hazard $h(x)=\frac{k}{\lambda}(x/\lambda)^{k-1}$ is monotone: decreasing ($k<1$, infant mortality), constant ($k=1$, the memoryless Exponential), or increasing ($k>1$, wear-out). Summing a decreasing and an increasing Weibull hazard produces the classic bathtub curve.

Log-logistic (α, β) — the density of a log-logistic random variable; its hazard is non-monotone for $\beta>1$ (rises then falls), the signature pattern for events like post-treatment mortality that peak and then subside. Closed-form quantile, so easy to simulate.

Gompertz–Makeham (λ, η, θ) — the actuarial mortality law: an exponentially increasing hazard $\eta e^{\theta x}$ (senescence) plus, in the Makeham extension, a constant $\lambda$ (age-independent accidents). Gompertz is the $\lambda=0$ case, sampled in closed form; the full Makeham has no elementary quantile and is sampled by numerical inverse-transform.

a,b = 3.0, 2.5; xs = np.linspace(0, 18, 500)
panel("Log-logistic(α=3, β=2.5)", xs, E.loglogistic_pdf(xs,a,b), E.loglogistic_cdf(xs,a,b),
      E.loglogistic_rng(a,b,150000,rng), E.loglogistic_moments(a,b),
      S.fisk.pdf(xs,b,scale=a), S.fisk.cdf(xs,b,scale=a), sp_stats(S.fisk(b,scale=a)), hist_range=(0,18))

                 analytical (formula)   scipy  empirical (RNG)
mean                            3.964   3.964            3.959
variance                       22.770  22.770           18.868
skewness                          NaN     NaN           14.280
excess kurtosis                   NaN     NaN          666.884

eta,theta = 0.5, 0.8; xs = np.linspace(0, 6, 500)
panel("Gompertz(η=0.5, θ=0.8)", xs, E.gompertz_pdf(xs,eta,theta), E.gompertz_cdf(xs,eta,theta),
      E.gompertz_rng(eta,theta,150000,rng), E.gompertz_moments(eta,theta),
      S.gompertz.pdf(xs,eta/theta,scale=1/theta), S.gompertz.cdf(xs,eta/theta,scale=1/theta), sp_stats(S.gompertz(eta/theta,scale=1/theta)),
      hist_range=(0,6))

                 analytical (formula)  scipy  empirical (RNG)
mean                            1.009  1.009            1.008
variance                        0.428  0.428            0.427
skewness                        0.524  0.524            0.521
excess kurtosis                -0.372 -0.372           -0.381

# THE signature figure: hazard shapes and the bathtub curve
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.4))
xs = np.linspace(0.02, 4, 400)
for k_,c,lab in [(0.6,BLUE,"Weibull k=0.6  (decreasing — infant mortality)"),
                 (1.0,GREY,"Weibull k=1  (constant — memoryless)"),
                 (2.5,RED,"Weibull k=2.5  (increasing — wear-out)")]:
    ax[0].plot(xs, E.weibull_hazard(xs,k_,1.5), lw=2.2, color=c, label=lab)
ax[0].plot(xs, E.loglogistic_hazard(xs,1.5,2.5), lw=2.2, color=GREEN, ls="--", label="Log-logistic  (non-monotone hump)")
ax[0].set_title("Hazard shapes h(x) = f(x)/S(x)"); ax[0].set_xlabel("x"); ax[0].set_ylabel("hazard"); ax[0].legend(fontsize=8); ax[0].set_ylim(0,3)
# bathtub = decreasing Weibull + constant + increasing Weibull (competing risks add hazards)
xs = np.linspace(0.02, 10, 400)
h_bath = E.weibull_hazard(xs,0.5,1.2) + 0.15 + E.weibull_hazard(xs,4.0,8.0)
ax[1].plot(xs, E.weibull_hazard(xs,0.5,1.2), color=BLUE, lw=1.3, ls=":", label="infant (k<1)")
ax[1].plot(xs, 0.15+0*xs, color=GREY, lw=1.3, ls=":", label="random (constant)")
ax[1].plot(xs, E.weibull_hazard(xs,4.0,8.0), color=RED, lw=1.3, ls=":", label="wear-out (k>1)")
ax[1].plot(xs, h_bath, color=PURP, lw=2.6, label="total = bathtub curve")
ax[1].set_title("The bathtub curve: sum of three hazards"); ax[1].set_xlabel("age"); ax[1].set_ylabel("hazard"); ax[1].legend(fontsize=8); ax[1].set_ylim(0,1.2)
plt.tight_layout(); plt.show()
print("Left: the hazard's shape classifies a lifetime law. Right: infant-mortality + random + wear-out hazards add to the bathtub.")

Left: the hazard's shape classifies a lifetime law. Right: infant-mortality + random + wear-out hazards add to the bathtub.

# Gompertz vs Gompertz–Makeham: the constant Makeham term lifts the whole hazard
xs = np.linspace(0, 6, 400)
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].plot(xs, E.gompertz_hazard(xs,0.3,0.9), color=BLUE, lw=2.2, label="Gompertz  η e^{θx}")
ax[0].plot(xs, E.makeham_hazard(xs,0.4,0.3,0.9), color=RED, lw=2.2, label="Makeham  λ + η e^{θx}")
ax[0].axhline(0.4, color=GREY, ls=":", lw=1.2, label="Makeham constant λ")
ax[0].set_title("Mortality hazard: Gompertz vs Gompertz–Makeham"); ax[0].set_xlabel("age"); ax[0].set_ylabel("hazard"); ax[0].legend(fontsize=8)
# survival curves
ax[1].plot(xs, 1-E.gompertz_cdf(xs,0.3,0.9), color=BLUE, lw=2.2, label="Gompertz survival")
ax[1].plot(xs, 1-E.makeham_cdf(xs,0.4,0.3,0.9), color=RED, lw=2.2, label="Makeham survival")
ax[1].set_title("Survival S(x) = 1 − F(x)"); ax[1].set_xlabel("age"); ax[1].set_ylabel("S(x)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
samp = E.makeham_rng(0.4, 0.3, 0.9, 120000, rng)
print(f"Makeham sampled by numerical inverse-transform: empirical mean {samp.mean():.3f} vs analytical {E.makeham_moments(0.4,0.3,0.9)['mean']:.3f}.")

Makeham sampled by numerical inverse-transform: empirical mean 0.939 vs analytical 0.940.

5. Actuarial loss & income — GB2 and its Burr / Dagum family

Insurance losses and incomes are positive, right-skewed and heavy-tailed. The Generalized Beta of the second kind, GB2(a, b, p, q), is the four-parameter umbrella that nests most of the standard models: $$f(x)=\frac{a\,x^{ap-1}}{b^{ap}B(p,q)\,\big(1+(x/b)^a\big)^{p+q}},\qquad x>0.$$ The shape $a$ and the two beta parameters $p,q$ independently control the peak, the left rise and the right tail (moments exist only for $-ap<k<aq$). Its special cases are the field's named laws:

• $p=1$ → Burr XII / Singh–Maddala (loss severity),
• $q=1$ → Dagum / Burr III (income distribution),
• $p=q=1$ → log-logistic.

We validate GB2 by drawing it as a transformed Beta ($X=b(W/(1-W))^{1/a}$, $W\sim\text{Beta}(p,q)$) and by checking its special cases collapse onto Burr, Dagum and log-logistic exactly.

c,k,bb = 3.0, 2.0, 2.0; xs = np.linspace(0, 12, 500)
panel("Burr XII (c=3, k=2, b=2)  [GB2 p=1]", xs, E.burr12_pdf(xs,c,k,bb), E.burr12_cdf(xs,c,k,bb),
      E.burr12_rng(c,k,bb,150000,rng), E.burr12_moments(c,k,bb),
      S.burr12.pdf(xs,c,k,scale=bb), S.burr12.cdf(xs,c,k,scale=bb), sp_stats(S.burr12(c,k,scale=bb)), hist_range=(0,12))

                 analytical (formula)  scipy  empirical (RNG)
mean                            1.612  1.612            1.610
variance                        0.625  0.625            0.627
skewness                        1.589  1.589            1.691
excess kurtosis                 7.809  7.809            9.993

a,bb,p = 3.0, 2.0, 1.5; xs = np.linspace(0, 12, 500)
panel("Dagum (a=3, b=2, p=1.5)  [GB2 q=1]", xs, E.dagum_pdf(xs,a,bb,p), E.dagum_cdf(xs,a,bb,p),
      E.dagum_rng(a,bb,p,150000,rng), E.dagum_moments(a,bb,p),
      S.burr.pdf(xs,a,p,scale=bb), S.burr.cdf(xs,a,p,scale=bb), sp_stats(S.burr(a,p,scale=bb)), hist_range=(0,12))

                 analytical (formula)  scipy  empirical (RNG)
mean                            2.875  2.875            2.871
variance                        4.824  4.824            4.518
skewness                          NaN    NaN            6.768
excess kurtosis                   NaN    NaN          119.040

# The GB2 family tree: one 4-parameter law contains the others
xs = np.linspace(0.01, 15, 500)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot(xs, E.gb2_pdf(xs,3,2,1,2), color=BLUE, lw=3, alpha=.5, label="GB2(3,2,p=1,q=2)")
ax[0].plot(xs, E.burr12_pdf(xs,3,2,2), color=RED, lw=1.6, ls="--", label="= Burr XII(3,2,2)")
ax[0].plot(xs, E.gb2_pdf(xs,3,2,1.5,1), color=GREEN, lw=3, alpha=.5, label="GB2(3,2,p=1.5,q=1)")
ax[0].plot(xs, E.dagum_pdf(xs,3,2,1.5), color=ORANGE, lw=1.6, ls="--", label="= Dagum(3,2,1.5)")
ax[0].set_title("GB2 special cases coincide with Burr & Dagum"); ax[0].set_xlabel("x"); ax[0].legend(fontsize=8)
# a full GB2 vs its heavy tail on log-log
ax[1].loglog(xs, E.gb2_pdf(xs,1.5,2,2,1.2), color=PURP, lw=2, label="GB2(1.5,2,2,1.2)")
ax[1].loglog(xs, E.gb2_pdf(xs,2.0,2,2,3.0), color=BLUE, lw=2, label="GB2(2,2,2,3)  (lighter tail)")
ax[1].set_title("GB2 tail on log–log axes (slope ≈ −a·q−1)"); ax[1].set_xlabel("x (log)"); ax[1].set_ylabel("f(x) (log)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("GB2 with p=1 IS the Burr XII; with q=1 it IS the Dagum — a single family behind the actuarial zoo.")

GB2 with p=1 IS the Burr XII; with q=1 it IS the Dagum — a single family behind the actuarial zoo.

6. Timing: inverse transforms dominate

Almost every law in this folder has a closed-form quantile — the EVT laws (Gumbel, Weibull, Fréchet, GEV, GPD) and the actuarial laws (Burr, Dagum) all invert their CDF in one expression, so their from-scratch samplers are single vectorised transforms that match or beat scipy. The GB2 draws a Beta (two Gammas); only the Gompertz–Makeham, whose quantile is transcendental, pays for a numerical inverse-transform.

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 = [
  ("Gumbel",      lambda: E.gumbel_rng(1,2,N,rng),        lambda: S.gumbel_r.rvs(1,2,size=N,random_state=rng),               "inverse transform"),
  ("Weibull",     lambda: E.weibull_rng(1.8,3,N,rng),     lambda: S.weibull_min.rvs(1.8,scale=3,size=N,random_state=rng),    "inverse transform"),
  ("Fréchet",     lambda: E.frechet_rng(2.5,2,1,N,rng),   lambda: S.invweibull.rvs(2.5,loc=1,scale=2,size=N,random_state=rng),"inverse transform"),
  ("GEV",         lambda: E.gev_rng(1,2,0.3,N,rng),       lambda: S.genextreme.rvs(-0.3,loc=1,scale=2,size=N,random_state=rng),"inverse transform"),
  ("GPD",         lambda: E.gpd_rng(0,1.5,0.3,N,rng),     lambda: S.genpareto.rvs(0.3,scale=1.5,size=N,random_state=rng),    "inverse transform"),
  ("Log-logistic",lambda: E.loglogistic_rng(3,2.5,N,rng), lambda: S.fisk.rvs(2.5,scale=3,size=N,random_state=rng),          "inverse transform"),
  ("Gompertz",    lambda: E.gompertz_rng(0.5,0.8,N,rng),  lambda: S.gompertz.rvs(0.625,scale=1.25,size=N,random_state=rng), "inverse transform"),
  ("Burr XII",    lambda: E.burr12_rng(3,2,2,N,rng),      lambda: S.burr12.rvs(3,2,scale=2,size=N,random_state=rng),        "inverse transform"),
  ("Dagum",       lambda: E.dagum_rng(3,2,1.5,N,rng),     lambda: S.burr.rvs(3,1.5,scale=2,size=N,random_state=rng),        "inverse transform"),
  ("GB2",         lambda: E.gb2_rng(1.5,2,2,3,N,rng),     None,                                                             "transformed Beta"),
  ("Makeham",     lambda: E.makeham_rng(0.4,0.3,0.9,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.4))
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 dominate; only Makeham 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"])
lo, hi = cmp_["ratio"].min(), cmp_["ratio"].max()
behind = cmp_[cmp_["ratio"] < 1].index.tolist()
slowest = df["from-scratch (M/s)"].idxmin()
print("\nEvery EVT and actuarial law here inverts its CDF in closed form, and it shows: measured against")
print("scipy's generic .rvs the from-scratch versions span %.2fx-%.2fx" % (lo, hi)
      + (", all of them at or above parity." if not behind else ", trailing only on: " + ", ".join(behind) + "."))
print("The two exceptions are GB2, which costs a Beta draw (two Gammas), and the Gompertz-Makeham, whose")
print("quantile is transcendental so it must build a CDF grid before drawing. Those are the only generators")
print("here without a closed-form inverse, and the slowest in the set is %s at %.1f M draws/s -- still fast for practical N."
      % (slowest, df["from-scratch (M/s)"].min()))

              from-scratch (M/s)  scipy (M/s)  ratio                 method
dist                                                                       
Gumbel                     57.23        47.75   1.20      inverse transform
Weibull                    51.22        36.75   1.39      inverse transform
Fréchet                    46.36        44.81   1.03      inverse transform
GEV                        36.95        21.10   1.75      inverse transform
GPD                        50.06        37.21   1.35      inverse transform
Log-logistic               54.75        47.23   1.16      inverse transform
Gompertz                   61.55        39.87   1.54      inverse transform
Burr XII                   35.27        25.72   1.37      inverse transform
Dagum                      38.28        34.54   1.11      inverse transform
GB2                         5.50          NaN    NaN       transformed Beta
Makeham                     8.29          NaN    NaN  numerical inverse-CDF

Every EVT and actuarial law here inverts its CDF in closed form, and it shows: measured against
scipy's generic .rvs the from-scratch versions span 1.03x-1.75x, all of them at or above parity.
The two exceptions are GB2, which costs a Beta draw (two Gammas), and the Gompertz-Makeham, whose
quantile is transcendental so it must build a CDF grid before drawing. Those are the only generators
here without a closed-form inverse, and the slowest in the set is GB2 at 5.5 M draws/s -- still fast for practical N.

7. Summary

Eleven laws spanning extreme value and survival/reliability, each with a from-scratch PDF / CDF / RNG / moments (and hazard) validated against scipy.stats to machine precision.

The two unifying ideas. For maxima, the GEV absorbs Gumbel ($\xi=0$), Fréchet ($\xi>0$) and Weibull ($\xi<0$), and by Fisher–Tippett–Gnedenko it is the universal limit of block maxima — its threshold twin, the GPD, is the universal limit of exceedances (Pickands–Balkema–de Haan), both demonstrated by simulation. For time-to-failure, the hazard function $h=f/S$ organises the survival laws by the shape of the failure rate: Weibull's monotone hazard (and the bathtub sum), the log-logistic's non-monotone hump, and Gompertz–Makeham's mortality law of exponential ageing plus constant accident risk. Finally, the actuarial GB2 family is the four-parameter umbrella whose $p=1$, $q=1$ and $p=q=1$ special cases are the Burr XII, Dagum and log-logistic — a single law behind the loss- and income-modelling zoo.

Next in the catalog: Bayesian priors & data augmentation — the truncated, half, Pólya-Gamma and generalized-inverse-Gaussian distributions that power modern Bayesian computation.