Density Functions and Random Number Generators of $α$-Stable Distributions

Abstract: Heavy-tailed distributions are increasingly found to better fit empirical data in engineering, finance, physics, network science, and related fields. Among them, $α$-stable distributions play a central role being limiting laws in the generalized central limit theorem: they are expected to be exceptionally good models whenever sums of multiple independent heavy-tailed sources are at play. Despite their theoretical importance, their practical use remains challenging: $α$-stable probability densities generally do not have closed-form expressions, and numerical evaluation and random variate generation can be difficult, especially in the multivariate setting. This paper presents AUB-HTP, a Python package for numerical computation and simulation of $α$-stable distributions. The package provides scalar density evaluation using several complementary methods, including Zolotarev-type integral representations, series formulas, and numerical inversion of characteristic functions. It also provides random variate generation for scalar and multivariate $α$-stable distributions, with support for flexible spectral measures through LePage series representations. Numerical experiments demonstrate that AUB-HTP improves the accuracy, stability, and parameter coverage of existing tools for scalar density computation, while adding new capabilities for multivariate simulation. The package is designed to support reproducible computational work involving heavy-tailed models across a broad range of scientific applications.

• The paper introduces AUB-HTP, a Python package that combines characteristic-function inversion, Zolotarev integrals, and tail series to evaluate scalar α-stable densities with a median speedup of approximately 58× over SciPy.
• The package generates univariate and multivariate α-stable random variables using CMS and truncated Lévy–LePage methods, supporting discrete, isotropic, elliptical, mixed, and arbitrary custom spectral measures.
• The paper derives an explicit truncation-error MSE bound and demonstrates broad practical utility, while identifying unresolved concerns around independent validation near α = 1, term-count limits, and future multivariate density evaluation.

AUB-HTP is a Python package for numerical evaluation of probability density functions (PDFs) and generation of random variates for $\alpha$-stable distributions, developed at the American University of Beirut. The package addresses two persistent practical obstacles in working with stable laws: the absence of closed-form densities (except for the Cauchy and Lévy cases) and the difficulty of simulating multivariate stable vectors with general spectral measures. The paper's contributions are fourfold: a unified scalar density evaluator combining several complementary numerical methods, a flexible scalar and multivariate random number generator built on LePage series, support for arbitrary spectral measures, and benchmarking against SciPy's reference implementation (2606.28530).

Background and motivation

$\alpha$-stable distributions are the limiting laws of the generalized central limit theorem and are closed under summation up to location and scale shifts, making them natural models for aggregates of independent heavy-tailed sources in communications, radar clutter, finance, and network traffic. The paper adopts the spectral-measure characterization of Samorodnitsky and Taqqu: a random vector is $\alpha$-stable ($0<\alpha<2$) if and only if its characteristic function is determined by a finite spectral measure $\Lambda$ on the unit sphere $\mathbb{S}^{d-1}$ and a shift vector $\boldsymbol{\mu}^0$. In the scalar case, $\Lambda$ reduces to two mass points and the distribution is parameterized by the stability index $\alpha$, skewness $\beta$, scale $\alpha$0, and location $\alpha$1, under either the type-0 (Nolan's numerically stable) or type-1 parameterization. The package supports both, which matters because the two parameterizations differ materially near $\alpha$2, where a logarithmic correction to the location appears.

Hybrid scalar density evaluation

The density evaluator is the paper's most technically involved component. Rather than relying on a single method, it dispatches each input point to one of three formulas according to a region determined by the standardized value $\alpha$3:

• Near the mode: direct numerical inversion of the $\alpha$4 characteristic function, which is stable where tail asymptotics fail but slower due to quadrature.
• Intermediate region: Zolotarev's integral representation, a compromise in both accuracy and speed.
• Tails: Skorohod–Pollard–Bergström series expansions of the form $\alpha$5, which avoid repeated integration entirely.

The switching points between regions are not fixed constants; they are precomputed as cutoff functions of $\alpha$6 via grid-based comparison against scipy.stats.levy_stable.pdf, and at runtime Boolean masks route each evaluation point to the appropriate method. Input standardization and the symmetry relation $\alpha$7 restrict all internal computation to the positive half-line of a standardized variable.

Validation against the SciPy reference on grids over $\alpha$8, $\alpha$9, and $\alpha$0 (with tail evaluations up to $\alpha$1) shows close agreement away from the neighborhood of $\alpha$2. The paper makes a notable claim here: the large relative discrepancies observed near $\alpha$3 are attributed not to errors in AUB-HTP but to numerical instability in the SciPy reference itself, which exhibits spurious zero windows and sharp discontinuities in regions where the true density is smooth. This is supported by visual diagnostic plots. The claim is credible in direction but rests on visual inspection rather than an independent ground-truth method (e.g., very high-precision quadrature), so it remains an assertion supported by evidence rather than a proof.

The runtime results are strong: the median speedup over scipy.stats.levy_stable.pdf is approximately $\alpha$4, with the largest gains in tail-heavy evaluations and for small $\alpha$5, where the series expansions replace quadrature. Speedups are smallest near $\alpha$6, where the package falls back on direct inversion because of the logarithmic correction term in the characteristic function.

Multivariate random variate generation

The univariate generator implements a numerically stable variant of the Chambers–Mallows–Stuck (CMS) method. The multivariate generator is built on the Lévy–LePage series representation: given i.i.d. vectors $\alpha$7 distributed according to the (normalized) spectral measure and i.i.d. unit-mean exponentials $\alpha$8 with $\alpha$9, the truncated series

$0<\alpha<2$0

converges in total variation to the target $0<\alpha<2$1-stable vector, for $0<\alpha<2$2 unconditionally and for $0<\alpha<2$3 provided $0<\alpha<2$4.

The package ships with samplers for several spectral-measure classes: discrete measures (categorical sampling over mass points), isotropic measures (uniform sampling on the sphere via Muller's method, with an explicit scaling factor $0<\alpha<2$5 to realize a user-specified scale $0<\alpha<2$6), sub-Gaussian/elliptical measures (via Cholesky factorization of the shape matrix $0<\alpha<2$7, with Monte Carlo estimation of the total spectral mass when unspecified, justified by a finite-second-moment argument), mixed measures (composition over weighted components), and custom measures through a user-extensible base class requiring only a sample method plus the dimension and total mass. The custom-measure interface is a genuine differentiator: existing tools largely handle only isotropic or sub-Gaussian spectral measures, whereas AUB-HTP supports arbitrary ones, subject to the zero-mean condition when $0<\alpha<2$8.

The paper derives an explicit MSE bound for the truncated LePage series (Theorem in the appendix): assuming $0<\alpha<2$9,

$\Lambda$0

The exponent $\Lambda$1 (for $\Lambda$2) implies the bound decays faster for smaller $\Lambda$3, so the method is most accurate precisely in the heavily heavy-tailed regime where other tools struggle. The sub-Gaussian case inherits the bound scaled by $\Lambda$4. In implementation, the package selects the number of terms to target an MSE below 0.01, capped at 50,000 terms with a warning when the cap is reached — an implicit acknowledgment that for some parameter combinations the theoretical term count may be impractically large.

Limitations and open questions

Several limitations are conceded or evident. The density evaluator's accuracy claims are benchmarked against SciPy rather than an independent high-precision reference, which is problematic exactly where the paper accuses SciPy of instability (near $\Lambda$5). The LePage MSE bound assumes unit-norm spectral vectors, so it does not directly cover all discrete or custom measures with non-unit masses. The 50,000-term cap means the 0.01 MSE target is not always guaranteed. The package currently provides PDF evaluation and sampling only; multivariate density computation, CDF evaluation, and parameter estimation are listed as future work, and the truncation-error analysis does not extend to derived functionals of the sampled vectors. Finally, the convergence condition $\Lambda$6 for $\Lambda$7 is left to the user to verify for custom spectral measures; whether the package enforces it is not stated.

Conclusion

AUB-HTP consolidates scalar $\Lambda$8-stable density evaluation into a hybrid, region-dispatched scheme that is substantially faster than the SciPy reference (median speedup $\Lambda$9) while exposing and diagnosing numerical defects in that reference near $\mathbb{S}^{d-1}$0. Its multivariate generator, grounded in the LePage series with an explicit finite-sample MSE bound, extends simulation capability to arbitrary spectral measures — a capability largely absent from existing Python tooling. The main open questions are the establishment of an independent ground truth for the sensitive $\mathbb{S}^{d-1}$1 region and the extension of the framework to multivariate densities, distribution functions, and estimation.