Basic Usage¶
A tour of heavytails from first import to a fitted tail model.
Everything here runs on the standard library alone — no NumPy, no SciPy, no plotting library. The notebook is executed when the documentation is built, so every number below was produced by the code above it.
import heavytails
heavytails.__version__
'0.6.3'
Creating a distribution¶
Every family is a frozen dataclass. Parameters are validated at construction, so an invalid instance can never exist.
from heavytails import Pareto
pareto = Pareto(alpha=1.5, xm=1.0)
pareto
Pareto(alpha=1.5, xm=1.0)
from heavytails.heavy_tails import ParameterError
# A negative shape parameter is rejected immediately, not at first use.
try:
Pareto(alpha=-1.0, xm=1.0)
except ParameterError as exc:
print(f"ParameterError: {exc}")
ParameterError: Pareto requires alpha>0 and xm>0.
The distribution interface¶
Every family exposes the same five methods.
x = 2.0
print(f"pdf({x}) = {pareto.pdf(x):.6f}")
print(f"cdf({x}) = {pareto.cdf(x):.6f}")
print(f"sf({x}) = {pareto.sf(x):.6f}")
print(f"ppf(0.99) = {pareto.ppf(0.99):.6f}")
pdf(2.0) = 0.265165 cdf(2.0) = 0.646447 sf(2.0) = 0.353553 ppf(0.99) = 21.544347
Why sf exists¶
In the far tail, 1 - cdf(x) is catastrophic: once cdf(x) rounds to 1.0, the complement is exactly zero and every digit is lost. The survival function is computed directly and stays accurate.
far_out = 1e12
print(f"1 - cdf({far_out:.0e}) = {1 - pareto.cdf(far_out):.6e}")
print(f"sf({far_out:.0e}) = {pareto.sf(far_out):.6e}")
1 - cdf(1e+12) = 0.000000e+00 sf(1e+12) = 1.000000e-18
Sampling¶
Pass a seed for reproducible draws. The same seed and parameters always produce the same sequence.
samples = pareto.rvs(10_000, seed=42)
print(f"n = {len(samples)}")
print(f"min = {min(samples):.4f}")
print(f"max = {max(samples):.4f}")
print(f"median = {sorted(samples)[len(samples) // 2]:.4f}")
n = 10000 min = 1.0000 max = 490.0822 median = 1.5888
# Reproducibility: same seed, same numbers.
pareto.rvs(5, seed=7) == pareto.rvs(5, seed=7)
True
Heavy tails behave differently¶
With alpha = 1.5 the mean exists but the variance does not. A sample mean therefore wanders instead of settling — this is not a bug in the sampler, it is what an infinite variance looks like.
running = []
total = 0.0
for i, value in enumerate(samples, start=1):
total += value
if i in (100, 500, 1_000, 5_000, 10_000):
running.append((i, total / i))
theoretical_mean = 1.5 / (1.5 - 1) * 1.0 # alpha * xm / (alpha - 1)
print(f"theoretical mean = {theoretical_mean:.4f}\n")
for n, mean in running:
print(f"mean of first {n:6,d} = {mean:.4f}")
theoretical mean = 3.0000 mean of first 100 = 3.2378 mean of first 500 = 3.1839 mean of first 1,000 = 3.7109 mean of first 5,000 = 3.1126 mean of first 10,000 = 2.9227
The largest few observations dominate the sum. That single fact drives most of the practical consequences of heavy tails.
ordered = sorted(samples, reverse=True)
total = sum(ordered)
for k in (1, 10, 100):
share = sum(ordered[:k]) / total
print(
f"top {k:3d} of {len(samples):,} observations carry {share:6.2%} of the total"
)
top 1 of 10,000 observations carry 1.68% of the total top 10 of 10,000 observations carry 7.44% of the total top 100 of 10,000 observations carry 19.25% of the total
Estimating the tail index¶
The estimators return the extreme-value index gamma = 1 / alpha. Invert it to read alpha back.
from heavytails import hill_estimator, moment_estimator, pickands_estimator
gamma_hill = hill_estimator(samples, k=100)
gamma_pickands = pickands_estimator(samples, k=100)
gamma_moment, alpha_moment = moment_estimator(samples, k=100)
print("true alpha = 1.5000\n")
print(f"Hill gamma = {gamma_hill:.4f} -> alpha = {1 / gamma_hill:.4f}")
print(f"Pickands gamma = {gamma_pickands:.4f} -> alpha = {1 / gamma_pickands:.4f}")
print(f"Moment gamma = {gamma_moment:.4f} -> alpha = {alpha_moment:.4f}")
true alpha = 1.5000 Hill gamma = 0.6516 -> alpha = 1.5347 Pickands gamma = 0.8111 -> alpha = 1.2329 Moment gamma = 0.6024 -> alpha = 1.6600
Choosing k¶
Every estimator depends on how many upper order statistics it uses. Small k means low bias and high variance; large k is the reverse. The standard tool is a Hill plot: sweep k and look for a plateau.
print(f"{'k':>6} {'gamma':>8} {'alpha':>8}")
for k in (25, 50, 100, 200, 400, 800, 1600):
gamma = hill_estimator(samples, k=k)
print(f"{k:6d} {gamma:8.4f} {1 / gamma:8.4f}")
k gamma alpha
25 0.6774 1.4762
50 0.6062 1.6497
100 0.6516 1.5347
200 0.6806 1.4693
400 0.6500 1.5384
800 0.6826 1.4650
1600 0.6858 1.4580
The estimates are reasonably stable in the middle of that range, which is the plateau to read from. At the small end the variance dominates; at the large end observations from the body of the distribution bias the estimate.
See Tail Index Estimation Theory for what is going on underneath.
Diagnostics¶
heavytails.plotting returns coordinates rather than figures, so it stays dependency-free. A power-law tail is a straight line on log–log axes.
from heavytails.plotting import tail_loglog_plot
points = tail_loglog_plot(samples)
print(f"{len(points):,} (log x, log P(X > x)) pairs")
# Fit a line through the tail region by least squares to recover the slope.
tail = points[-2000:-50]
n = len(tail)
mean_x = sum(px for px, _ in tail) / n
mean_y = sum(py for _, py in tail) / n
cov = sum((px - mean_x) * (py - mean_y) for px, py in tail)
var = sum((px - mean_x) ** 2 for px, _ in tail)
slope = cov / var
print(f"slope of the log-log tail = {slope:.4f} (expected about -1.5)")
10,000 (log x, log P(X > x)) pairs slope of the log-log tail = -1.4480 (expected about -1.5)
To draw it, hand the coordinates to any plotting library:
import matplotlib.pyplot as plt
xs, ys = zip(*points)
plt.plot(xs, ys, ".", markersize=2)
plt.xlabel(r"$\log x$")
plt.ylabel(r"$\log P(X > x)$")
See Tail Diagnostics for how to read the result.
Comparing families¶
Several heavy-tailed families can look alike over a limited range. Comparing them on the same data is more honest than fitting one and declaring success.
from heavytails import BurrXII, Cauchy, LogNormal, StudentT
families = {
"Pareto(1.5, 1)": Pareto(alpha=1.5, xm=1.0),
"LogNormal(0, 1)": LogNormal(mu=0.0, sigma=1.0),
"StudentT(2)": StudentT(nu=2.0),
"Cauchy(0, 1)": Cauchy(x0=0.0, gamma=1.0),
"BurrXII(1.2, 2.5, 3)": BurrXII(c=1.2, k=2.5, s=3.0),
}
print(f"{'family':<22} {'median':>10} {'p99':>12} {'p99.9':>14}")
for name, dist in families.items():
print(
f"{name:<22} {dist.ppf(0.5):>10.4f} {dist.ppf(0.99):>12.4f} {dist.ppf(0.999):>14.4f}"
)
family median p99 p99.9 Pareto(1.5, 1) 1.5874 21.5443 100.0000 LogNormal(0, 1) 1.0000 10.2405 21.9822 StudentT(2) 0.0000 6.9646 22.3271 Cauchy(0, 1) 0.0000 31.8205 318.3088 BurrXII(1.2, 2.5, 3) 1.1593 12.0597 28.4141
The medians are comparable while the extreme quantiles differ by orders of magnitude. Choosing the wrong family barely affects the centre and completely changes the tail — which is the whole reason to model the tail explicitly.
Discrete families¶
The discrete distributions follow the same interface, with pmf in place of pdf.
from heavytails import Zipf
zipf = Zipf(s=2.0, kmax=1000)
print(f"{'k':>4} {'pmf(k)':>10}")
for k in (1, 2, 5, 10, 100):
print(f"{k:4d} {zipf.pmf(k):10.6f}")
k pmf(k) 1 0.608297 2 0.152074 5 0.024332 10 0.006083 100 0.000061
Where to go next¶
- Distributions Overview — every family and its parameters
- Tail Index Estimation — the practical estimation guide
- Tail Diagnostics — checking for a power law before fitting
- Parameter Fitting — MLE, method of moments and model comparison
- Extreme Value Theory — the framework behind it all
- CLI Reference — the same workflow from a terminal