Density Functions and Random Number Generators of Alpha-Stable Distributions

1Introduction¶

Heavy-tailed distributions arise naturally in many scientific and engineering applications, including communications, signal processing, finance, physics, and network modeling Nassar et al. (2008)Tsihrintzis et al. (1995)Crovella et al. (1998)Nolan (2013)Taleb (2022). In such settings, rare but large deviations occur more frequently than predicted by light-tailed models such as the Gaussian distribution. As a result, Gaussian-based modeling can significantly underestimate the probability and impact of extreme events.

Among heavy-tailed models, $\alpha$-stable distributions occupy a central theoretical position. They generalize the Gaussian distribution through the stability property: sums of independent $\alpha$-stable random variables remain $\alpha$-stable, up to location and scale changes. They also arise as limiting distributions for normalized sums of independent random variables under the generalized central limit theorem B. V. Gnedenko and A. N. Kolmogorov (1968). These properties make $\alpha$-stable laws natural candidates for modeling aggregate effects generated by many independent heavy-tailed sources.

A major practical limitation of $\alpha$-stable distributions is that their probability density functions (PDF)s generally do not admit closed-form expressions, except in the special cases of the Cauchy and Lévy distributions. Consequently, numerical tools are required for density evaluation, modeling, and simulation. Existing approaches for scalar $\alpha$-stable random variate generation, such as the Chambers–Mallows–Stuck (CMS) method Chambers et al. (1976), are well established, but they do not directly address general multivariate $\alpha$-stable distributions. In the multivariate case, available simulation methods typically focus on special spectral measures, such as isotropic or sub-Gaussian models Nolan (2013), rather than arbitrary spectral measures.

This paper presents AUB-HTP, a Python package for computational work with $\alpha$-stable distributions. The package provides scalar density computation using several numerical methods, including Zolotarev integral representations Zolotarev (1986), Skorohod-Pollard-Bergström formulas A. V. Skorohod (1961)Pollard (1946)Bergström (1952), and direct numerical inversion of the characteristic function^[1]. These methods provide complementary tradeoffs in terms of accuracy, speed, and robustness across parameter regimes. In addition, AUB-HTP implements random variate generation for scalar and multivariate $\alpha$-stable distributions. For the multivariate case, the package uses LePage series representations LePage et al. (1981), which support flexible spectral measures and perform particularly well in challenging heavy-tailed regimes, including small values of $\alpha$.

The goal of AUB-HTP is to make $\alpha$-stable modeling more accessible to computational scientists by providing reliable density evaluation, flexible simulation tools, and reproducible numerical routines within the Python ecosystem. The main contributions of this work are: (i) a unified Python implementation of several scalar $\alpha$-stable PDF evaluation methods; (ii) a flexible random variate generator for scalar and multivariate $\alpha$-stable distributions; (iii) support for general spectral measures through LePage series; and (iv) numerical benchmarks comparing accuracy, runtime, and robustness against existing implementations.

Notations¶

We adopt the following conventions. The vector $\mathbf{x}\in\mathbb{R}^d$ represents a column vector and $\mathbf{x}^T$ is its transpose. The Euclidean inner product of two vectors $\mathbf{x} = (x_1,x_2,\cdot\cdot\cdot,x_d)^T$ and $\mathbf{y} = (y_1,y_2,\cdot\cdot\cdot,y_d)^T$ is denoted by $\mathbf{x}^T\mathbf{y}=\langle\mathbf{x},\mathbf{y}\rangle=\sum_{i=1}^{d}x_iy_i$ and the $L_2$ norm of $\mathbf{x}$ by $\|\mathbf{x}\|=\sqrt{\langle\mathbf{x},\mathbf{x}\rangle}$. The unit sphere in $\mathbb{R}^d$ is denoted by $\mathbb{S}^{d-1}=\{\mathbf{x}\in\mathbb{R}^d: \|\mathbf{x}\|=1\}$. The class of all borel subsets in a metric space $S$ is represented by $\mathcal{B}(S)$. We use $\mathbf{X} \perp\!\!\!\perp \mathbf{Y}$ to indicate that the two random vectors $\mathbf{X}$ and $\mathbf{Y}$ are independent. When a random vector $\mathbf{X}$ is distributed uniformly on the unit sphere, we write $\mathbf{X}\sim\mathcal{U}(\mathbb{S}^{d-1})$. For a random variable $X$ that is Gamma distributed, we write $X\sim\text{Gamma}(\alpha,\beta)$, and we use the following parameterization for the PDF

where $\Gamma(\cdot)$ is the Gamma function.

2Preliminaries on $\alpha$-Stable Distributions¶

We adopt the definition of a multivariate stable vector given in Samorodnitsky & Taqqu (2000):

In particular, we are interested in $\alpha$-stable ^[2] random vectors.

In the scalar case ($d = 1$) –and whenever $0<\alpha<2$, the spectral measure boils down to two mass points on $s = \pm 1$, and the characteristic function in (2) is expressible in terms of 4 parameters: the stability index $0<\alpha<2$, skewness $-1\leq\beta\leq1$, scale $\gamma>0$, and location $\delta\in\mathbb{R}$. In this case, we denote an $\alpha$-stable random variable $X$ using the same notation as in Nolan Nolan (2020), $X\!\sim\!\mathcal{S}(\alpha,\beta,\gamma,\delta;k)$ where $k$ is the parameterization index (or type). We will use the symbols $\mathcal{S}_0$ and $\mathcal{S}_1$ to refer to the two common type-0 and type-1 parameterization, respectively. In the multivariate case, we denote a $d$-dimension non-degenerate $\alpha$-stable random vector $\vect{X}$ with stability index $\alpha$, spectral measure $\Lambda$ and shift vector $\boldsymbol{\mu}^0$, by $\vect{X}\sim\mathbf{S}_\alpha(\Lambda,d,\boldsymbol{\mu}^0)$. Note that Definition 2 corresponds to the generalization of the scalar case under the parameterization $k=1$.

Parameterizations of scalar $\alpha$-stable distributions¶

There two parameterizations $\mathcal{S}(\alpha, \beta, \gamma, \delta; 0)$ and $\mathcal{S}(\alpha, \beta, \gamma, \delta; 1)$ differ primarily in how they define the location (shift) of the distribution, particularly near $\alpha = 1$. For a given $\alpha$ and $\beta$, let $Z$ be a “standard" $\alpha$-stable variable corresponding to $\gamma = 1$ and $\delta = 0$.

• Parameterization 0:

  A random variable $X \sim \mathcal{S}(\alpha, \beta, \gamma, \delta; 0)$ if

  $$X \overset{d}{=} \begin{cases} \gamma \left( Z - \beta \tan\left( \frac{\pi \alpha}{2} \right) \right) + \delta, & \alpha \neq 1 \\ \gamma Z + \delta, & \alpha = 1. \end{cases}$$

  (3)

  The characteristic function of $X$ is given by

  $$\mathbb{E}[e^{iuX}] = \begin{cases} \exp\left( -|\gamma u|^\alpha \left[ 1 + \ii\beta \tan\left( \frac{\pi \alpha}{2} \right) \operatorname{sign}(u)(|\gamma u|^{1 - \alpha} - 1) \right] + \ii \delta u \right), & \alpha \neq 1 \\ \exp\left( -|\gamma u| \left[ 1 + \ii\beta \frac{2}{\pi} \operatorname{sign}(u) \log|\gamma u| \right] + \ii \delta u \right), & \alpha = 1. \end{cases}$$

  (4)

  This parameterization has the advantage of avoiding discontinuities near $\alpha = 1$ by adjusting the shift behavior. It is often preferred for numerical stability Nolan (2020).

• Parameterization 1:

  A random variable $X \sim \mathcal{S}(\alpha, \beta, \gamma, \delta; 1)$ if:

  $$X \overset{d}{=} \begin{cases} \gamma Z + \delta, & \alpha \neq 1 \\ \gamma Z + \left( \delta + \beta \frac{2}{\pi} \gamma \log \gamma \right), & \alpha = 1, \end{cases}$$

  (5)

  with a corresponding characteristic function given by

  $$\mathbb{E}[e^{iuX}] = \begin{cases} \exp\left( -|\gamma u|^\alpha \left[ 1 - \ii\beta \tan\left( \frac{\pi \alpha}{2} \right) \operatorname{sign}(u) \right] + \ii \delta u \right), & \alpha \neq 1 \\ \exp\left( -|\gamma u| \left[ 1 + \ii\beta \frac{2}{\pi} \operatorname{sign}(u) \log|u| \right] + \ii \delta u \right), & \alpha = 1. \end{cases}$$

  (6)

  Unlike type-0, this parameterization requires an explicit adjustment to the location when $\alpha = 1$.

Note that the choice of parameterization directly affects the shift and symmetry of the distribution, especially when $\beta \neq 0$ and $\alpha$ is close to 1. In both forms, the distribution can be standardized by setting $\gamma = 1$ and $\delta = 0$, in which case we use the abbreviations $\mathcal{S}(\alpha, \beta; 0)$ and $\mathcal{S}(\alpha, \beta; 1)$ respectively.

In the univariate case, our package can simulate both parameterizations $k=0,1$.

3Scalar Density Evaluation¶

The AUB-HTP package provides numerical evaluation of scalar $\alpha$-stable probability density functions through a SciPy-like interface. This component complements the random number generator: the generator produces samples, while the density module provides the reference curve used for model verification, visualization, and numerical comparison.

Except for Cauchy and Levy laws, $\alpha$-stable distributions do not have a closed-form density, which must therefore be computed using numerical inversion or analytic representations that are valid in specific parameter regions.

Towards fulfilling our main objective of reducing expensive quadrature calls while preserving accuracy across a wide range of $\alpha$, $\beta$, and sample values, the density evaluator combines three methods:

• direct inversion of Nolan’s $\mathcal{S}_1$ characteristic function near the mode Nolan (2020)

• Zolotarev’s integral representation in the middle range Zolotarev (1986)

• Skorohod-Pollard-Bergström’s series expansions in the tails A. V. Skorohod (1961)Wintner (1941)Pollard (1946).

3.1Input standardization¶

All density calls are first mapped to a standardized variable allowing the internal routines to work on standardized densities while still supporting user-facing location and scale parameters. More specifically, given an input sample value $x$ the code evaluates

then rescales the final density by $1/\gamma$. For the $\mathcal{S}_1$ parameterization for example, the shift is given by Nolan (2020, equation 1.5)

as can be deduced from equation (5).

Additionally, if $f_{\alpha,\beta}(x)$ is the PDF of a normalized $\alpha$-stable random variable, for negative evaluation points the code uses the symmetry relation Nolan (2020, proposition 1.1)

and most formulas only need to be evaluated on the positive half-line.

3.2Hybrid evaluation strategy¶

The main density function acts as a dispatcher. It divides the input array into regions and applies the method that is most accurate for each region.

• For small $|z|$ (see equation (7)), the package uses direct inversion of the characteristic function. While this method is slower than the series expansions, it is stable near the mode where tail formulas are not accurate.

• For intermediate values of $|z|$, the package uses Zolotarev’s integral representation. This provides a useful compromise between direct inversion and asymptotic series methods.

• For large $|z|$, the package uses Skorohod-Pollard-Bergström series’ expansions. These formulas are fast in the tails and avoid repeated numerical integration.

The switching points between methods are not fixed constants. They are instead selected from precomputed cutoff functions depending on $(\alpha,\beta)$. These cutoffs were obtained using grid-based comparison with scipy.stats.levy_stable.pdf. At runtime, the code queries the corresponding cutoff and applies Boolean masks to route each input point to the selected method.

Density formulas¶

$\bullet$ For direct inversion in the $\mathcal{S}_1$ parameterization Nolan (2020), the density is computed from the characteristic function. The standardized density is

These integrals are evaluated with scipy.integrate.quad_vec.

$\bullet$ For the middle region, the code uses Zolotarev’s integral representation Zolotarev (1986). For $\alpha\neq1$ and $x>0$,

where

and where $V(\theta)$ is evaluated following Zolotarev’s formula Zolotarev (1986). This integral is computed using the function scipy.integrate.quad.

$\bullet$ For the tail region, we use Skorohod-type series expansions A. V. Skorohod (1961)Wintner (1941)Pollard (1946) of the form

where

The number $N$ of terms used to calculate the series is appropriately chosen depending on the parameter range. In the current package, the $\alpha<1$ branch uses a larger truncation level, while the $\alpha>1$ branch requires fewer terms in the tested regimes.

3.3Validation¶

The density evaluator was tested against scipy.stats.levy_stable.pdf on grids of $\alpha$, $\beta$, and $x$. The near-mode tests used $x\in[-10,10]$, while the tail tests examined values up to $x = 100$.

The comparison recorded the following metrics:

• Maximum absolute error:

  $$E_{\max}^{\mathrm{abs}} = \max_{1 \leq i \leq n} \left| \hat{f}(x_i)-f_{\mathrm{SciPy}}(x_i) \right|.$$

  (15)

• Maximum relative error:

  $$E_{\max}^{\mathrm{rel}} = \max_{1 \leq i \leq n} \frac{ \left| \hat{f}(x_i)-f_{\mathrm{SciPy}}(x_i) \right| }{ \max\left( \left|f_{\mathrm{SciPy}}(x_i)\right|, \varepsilon \right) },$$

  (16)

  where $\varepsilon>0$ prevents division by zero.

• Mean absolute error:

  $$E_{\mathrm{mean}}^{\mathrm{abs}} = \frac{1}{n} \sum_{i=1}^{n} \left| \hat{f}(x_i)-f_{\mathrm{SciPy}}(x_i) \right|.$$

  (17)

• Runtime against SciPy, reported through the speedup factor:

  $$S = \frac{ T_{\mathrm{SciPy}} }{ T_{\mathrm{AUB\text{-}HTP}} }.$$

  (18)

The tests show close agreement across most of the parameter space, but they also reveal regions where the SciPy reference exhibits irregular numerical behavior. The strongest agreement occurs away from the narrow region around $\alpha=1$, where changes in the parameterization and the presence of logarithmic terms make density evaluation more sensitive. In this region, visual diagnostics are important because the relative difference may become large when the reference density is close to zero or numerically unstable. Figure 1 illustrates this behavior near $\alpha=1$ using the maximum relative difference measured around the mode. Away from $\alpha=1$, Figure 1 shows agreement with the SciPy reference. Results for the remaining error metrics, parameter ranges, and tail evaluations, together with the complete benchmark data, are available in the Github repository^[3].

3.4Discrepancies in the SciPy reference¶

While scipy.stats.levy_stable.pdf The SciPy Community (2026) is used as the main numerical reference, we observed that it has some discrepancies for certain parameter values. These discrepancies appear mostly near the sensitive region around $\alpha=1$, and they can create misleading relative difference values.

The main issue is that the SciPy PDF sometimes shows unusual zero or near-zero regions in some places while the densities are known to be smooth. Since the maximum relative difference contains a division by the reference value, even a small disagreement in such regions can produce a very large relative difference. Therefore, some of the large relative differences in the benchmark tables do not necessarily mean that our PDF generator is inaccurate. Instead, they are likely due to the instability or discontinuity in the reference PDF.

To validate this, we generated visual diagnostic plots comparing the AUB-HTP with the levy_stable PDF values. Figure 2 shows several cases where the SciPy reference generates unexpected zero windows or sharp changes, while the PDF generated using our package remains smooth.

The heatmaps are useful for locating sensitive parameter regions, but the diagnostic plots were needed to understand whether the difference is primarily attributed to our method or the numerical reference.

3.5Runtime¶

The adopted hybrid strategy reduces runtime by avoiding direct quadrature when a faster formula is accurate enough. The largest gains occur in the tail region, where series expansions replace repeated numerical integration.

In the benchmark runs, the median speedup was approximately $58\times$, with larger speedups in tail-heavy evaluations. These gains come from three design choices:

• vectorized array operations

• region-specific method dispatch

• reuse of precomputed cutoff functions.

The speedup plot in Figure 3 shows that our method is strictly faster than scipy.stats.levy_stable.pdf for almost all tested values of $\alpha$. Note that the vertical axis is on a logarithmic scale and values between 10 and 100 represent large runtime improvements. For small values of $\alpha$, especially in the near-mode range, the speedup is very high. This is because the method selection avoids unnecessary numerical integration and uses the formulas that are faster in those regions. In the tail range, the speedup is also consistent for most values of $\alpha$, since the Skorohod-Pollard-Bergström formulas replace repeated quadrature calls. The smallest speedups occur near $\alpha=1$ which is expected since this region is numerically delicate. The formulas change around $\alpha=1$, and the characteristic function contains a logarithmic correction term. The package therefore relies more on direct inversion in this region, which is slower than the series-based methods.

4Generating $\alpha$-Stable Random Vectors¶

In the univariate case, the package implements a numerically stable version of the Chambers, Mallows, and Stuck method which is presented in Chambers et al. (1976)Nolan (2020). In the multivariate setup, the proposed random number generator is based on a theoretical foundation that we state next.

4.1The LePage series¶

The LePage series Samorodnitsky & Taqqu (2000) was first introduced in LePage et al. (1981) and allows for simulating a multivariate stable vector with a stability index $\alpha$ and a given spectral measure through the use of the following theorem Bentkus et al. (2001, Theorem 4)

Theorem 1 indicates that sampling a multivariate stable vector can be done by computing the truncated series as in equation (20). Since the vector $\mathbf{X}$ resulting from the series has a zero shift vector $\boldsymbol{\mu}^0 = {\bf 0}$ Bentkus et al. (2001, Remark 1), we simply add $\boldsymbol{\mu}^0$ after sampling $\mathbf{X}_n$ using the LePage series, when needed. The analysis of the truncation error is deferred to Appendix 7.

In the following section, we present algorithms to sample according to well-known spectral measures, namely: Isotropic, Sub-Gaussian, Discrete, and Mixed spectral measures.

4.2Common spectral measures¶

4.2.1Discrete spectral measures¶

In the case of a discrete spectral measure with $M$ mass points, we store the point masses in an array $[\mathbf{v}_1, \mathbf{v}_2, \cdots, \mathbf{v}_M]$ and draw randomly with the probability of each vector $\mathbf{v}_{i}$ being

4.2.2Isotropic spectral measures¶

For an isotropic spectral measure, we sample the vectors $\{\mathbf{V}_k\}_{1 \leq k \leq n}$ uniformly on $\mathbb{S}^{d-1}$. To this end, several algorithms are available and we refer the reader to Schnabel & Janke (2022) where various such algorithms are discussed. In our computations, we adopt the method presented in Muller (1959) due to its simplicity and computational efficiency. Our package sets the spectral measure of the sphere to 1, and then scales the truncated LePage series by an appropriate factor $q$ to adjust for the desired scale $\gamma$ given by the user: When the spectral measure of the sphere is 1, it can be verified that the characteristic function is given by

To achieve a scale $\gamma$, we therefore multiply the truncated series by

4.2.3Sub-Gaussian (elliptical) spectral measures¶

A sub-Gaussian random vector is characterized by a positive-definite $d\times d$ shape matrix $\Sigma$. If $\vect{X}$ is a unit-scale isotropic multivariate stable random vector with stability index $\alpha$ then $\vect{Y}=\Sigma^{\frac{1}{2}}\vect{X}$ is a sub-Gaussian random vector with shape matrix $\Sigma$ where $\Sigma^{\frac{1}{2}}$ is the Cholesky decomposition of $\Sigma$ Nolan (2013). As such, to sample from $\mathbf{Y}_n$, one can first sample a unit-scale isotropic random vector with the same stability index $\alpha$ using the LePage series as in Section 4.2.2 and multiply it by the Cholesky decomposition $\Sigma^{\frac{1}{2}}$ of $\Sigma$. In section 7.2 in the appendix, we analyze the resulting Mean Square Error (MSE) of $\mathbf{Y}_n$.

A useful alternative would be to sample the vectors $\mathbf{V}_k$ for the sub-Gaussian case by sampling vectors $\mathbf{U}_k$ uniformly distributed on the unit sphere and multiplying them by $\Sigma^{\frac{1}{2}}$. Indeed, if $\mathbf{X}_n$ is the $n$-terms LePage series, by Theorem 1, $\mathbf{X}_n$ converges in total variation to $\mathbf{X}$, and

converges in total variation to $\mathbf{Y}=\Sigma^{\frac{1}{2}}\mathbf{X}$, which is the desired sub-Gaussian random vector.

The package first samples ${\bf V}_k$ and then computes the sum presented in (25) to sample a sub-Gaussian random vector ${\bf Y}_n$.

4.2.4Mixed spectral measures¶

Now consider the scenario where the spectral measure at hand is a mixture of $L$ different spectral measures: $\Lambda_1, \Lambda_2, \cdots, \Lambda_L$. If $w_1, w_2, \cdots, w_L$ are their respective associated weights, the mixed spectral measure $\Lambda_{\text{mix}}$ is given by

To sample from this spectral measure, we proceed as follows.

• For each weight $w_i$ we associate its corresponding probability as

  $$p_i = \frac{w_i}{\sum_{\ell=1}^L w_\ell} \qquad i=1, \cdots, L.$$

  (27)

  Then, we sort the probabilities in non-decreasing order that we assume –without loss of generality– to be $p_1, \cdots, p_L$ hereafter.

• We draw $U\sim\mathcal{U}[0,1]$,

  • If $U\leq p_1$ sample from $\Lambda_1$,

  • else if $U\leq p_1+p_2$, sample from $\Lambda_2$,

  • else iterate on the remaining spectral measures.

This algorithm achieves a sampling from the mixture of $\Lambda_1, \cdots, \Lambda_L$ with probabilities $p_1, \cdots, p_L$ respectively.

5Implementation¶

Whenever the LePage series is used, the package chooses the number of terms to achieve an MSE less than 0.01 using Theorem 2 in Appendix 7. This number is capped at 50,000 terms and a warning is issued whenever this cap is reached.

5.1Univariate random number generator¶

In the univariate case, the software samples $X \sim\mathcal{S}(\alpha,\beta,\gamma,\delta;k)$, and the type 1 parameterization is used by default.

import aub_htp as ht
alpha=1
beta=0
loc=0
scale=1
n=500000
samples=ht.alpha_stable.rvs(alpha=alpha, beta=beta, loc=loc, scale=scale, size=n)

In Figure 4, we present some examples where univariate $\alpha$-stable variables were sampled using our implementation, and plot the truncated histogram of values. Each histogram is overlayed by the pdf values which are computed using the SciPy levy_stable numerical package.

An example using the type 0 parameterization is presented in Figure 5.

import aub_htp as ht
alpha=0.1
beta=1
loc=-500
scale=10
n=500000
ht.alpha_stable.with_parameterization("S0")
samples=ht.alpha_stable.rvs(alpha=alpha, beta=beta, loc=loc, scale=scale, size=n )

5.2Multivariate random number generator¶

5.2.0.1Isotropic $\alpha$-stable random vector:¶

An isotropic $\alpha$-stable random vector with characteristic function

can be generated using the AUB-HTP package as follows:

import aub_htp as ht
from aub_htp.random import IsotropicSampler
alpha = 0.5
gamma = 1
shift= [0,0]
d = 2
n=500000
sampler = IsotropicSampler(number_of_dimensions= d, alpha = alpha, gamma = gamma)
samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler = sampler,shift=shift, size = n)

where $\alpha$ is the stability index, $\gamma$ is the scale parameter and $\boldsymbol{\mu}^0$ represents the shift vector. Figure 6 depicts scatter plots of the generated isotropic $\alpha$-stable random vectors.

5.2.0.2Elliptic $\alpha$-stable random vector:¶

An elliptic $\alpha$-stable random vector with shape matrix $\Sigma$ with the following characteristic function

can also be generated using AUB-HTP. We note that the elliptic sampler takes the total mass of the unit sphere as input. If the total mass of the sphere is not specified, the software estimates it using Monte Carlo simulations (with 10⁶ samples) as follows:

where the equality is due to  Samorodnitsky & Taqqu (2000, proposition 2.5.8), and where $\mathbf{U}_1,\mathbf{U}_2,\ldots\mathbf{U}_n$ are i.i.d. samples according to $\mathbf{U}\sim\mathcal{U}(\mathbb{S}^{d-1})$. Equation (30) is justified by the application of the Weak Law of Large Numbers (WLLN) because the random variable $\|\Sigma^{\frac{1}{2}}\mathbf{U}\|^{\alpha}$ has a finite second moment. Indeed,

where equality $(a)$ is due to the Cholesky decomposition definition and in $(b)$ we used the $\Sigma$-weighted norm $\|\mathbf{x}\|^2_{\Sigma}\triangleq\mathbf{x}^T\Sigma\mathbf{x}$. Finally, the inequality in $(c)$ is due to the fact that $\|\mathbf{x}\|^2_{\Sigma}\leq\lambda_{\max}(\Sigma)\|\mathbf{x}\|^2$ where $\lambda_{\text{max}}(\Sigma)$ is the largest eigenvalue of $\Sigma$, equal to the norm of matrix $\Sigma$ Horn & Johnson (1985).

Finally, we note that the MSE resulting from estimating the total mass of the sphere using equation (30) is given by Niederreiter (1992, Theorem 1.1).

An example code is as follows.

import aub_htp as ht
from aub_htp.random import EllipticSampler
alpha = 0.5
shift= [0,0]
d = 2
n=500000
sigma1=[[2,0.8],[0.8,1.5]]
sampler =EllipticSampler(number_of_dimensions = 2, alpha = alpha, sigma = sigma1)
samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler = sampler,shift=shift, size = n)

Figure 7 shows the scatter plots of the elliptical $\alpha$-stable random vector, where we use the following shape matrices

5.2.0.3Discrete spectral measures:¶

For an $\alpha$-stable random vector that has a discrete spectral measure with point masses at positions $\vect{v}_1, \cdots, \mathbf{v}_M$ and corresponding weights $w_1, \cdots, w_M$, the corresponding characteristic function is given by

Such a vector can be generated as follows:

import aub_htp as ht
from aub_htp.random import DiscreteSampler
alpha = 0.2
shift= [0,0]
d = 2
n=500000
positions=[[1,0],[0,1],[-1,0],[0,-1]]
weights=[1,1,1,1]
sampler = DiscreteSampler(alpha = alpha, positions = positions,  weights = weights)
samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler = sampler,shift=shift, size = n)

with corresponding scatter plots shown in Figure 8.

5.2.0.4Mixed spectral measure:¶

For a mixed spectral measure, let $\Lambda_1, \cdots, \Lambda_L$ be $L$ spectral measures with associated weights $w_1, \cdots, w_L$, and consider their mixture $\Lambda_{\text{mix}} = \sum_{i=1}^L w_i \Lambda_i$. To sample according to this spectral measure, we proceed as follows:

import aub_htp as ht
from aub_htp.random import DiscreteSampler, IsotropicSampler, MixedSampler
alpha = 0.3
shift= [1000,1000]
d = 2
n=500000
gamma=1

# Discrete Spectral Measure
positions=[[1,0],[0,1],[-1,0],[0,-1]]
weights=[1,1,1,1]
sampler1 = DiscreteSampler(alpha = alpha, positions = positions,  weights = weights)

# Isotropic Spectral Measure
sampler2= IsotropicSampler(number_of_dimensions=d, alpha=alpha, gamma=gamma)

#Spectral Measures and Their Associated Weights
spectral_measures=[sampler1, sampler2]
measure_weights=[0.5,0.5]

mixed_sampler=MixedSampler(spectral_measures=spectral_measures,weights=measure_weights)
samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler = mixed_sampler, shift = shift, size = n)

where spectral_measures is a list containing the spectral measures that we are mixing, and measure_weights contains the associated weight of each measure. The resulting scatter plot is shown in Figure 9a.

5.2.0.5Custom spectral measures:¶

The common spectral measures presented above are built-in into AUB-HTP. If the user wants to sample from an $\alpha$-stable random vector whose spectral measure is not built-in, this is made possible using a systematic procedure outlined in the following example:

• Say the user is interested in sampling a multivariate $\alpha$-stable random vector in dimension $d=2$ with $\alpha=0.3, \boldsymbol{\mu}^0=\mathbf{0}$, and whose spectral measure is uniform over the arcs $\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$ and $\left[\frac{3\pi}{4},\frac{5\pi}{4}\right]$, that is, with probability 0.5, we draw a point uniformly on the arc $\left[\frac{3\pi}{4},\frac{5\pi}{4}\right]$, and with probability 0.5, we draw a point uniformly on the arc $\left[-\frac{\pi}{4},\frac{\pi}{4}\right]$.

• We assign equal probabilities to both arcs and the spectral measure of each arc is half the total measure of the sphere. Generally, the probability of each region is the ratio of the spectral measure of the region to the total mass of the sphere.

• We create a class of the spectral measure sampler that we want by inheriting from the base class
  BaseSpectralMeasureSampler. In the class that the user is creating, the function which samples the vectors $\mathbf{V}$ in Theorem 1 should be provided. Moreover, the dimension of the desired alpha-stable vector and the total mass of the unit sphere should be specified. The created spectral measure sampler is then passed as an argument to the rvs() method as shown in the following code.

The steps are presented below.

import aub_htp as ht
from aub_htp.random import BaseSpectralMeasureSampler
alpha = 0.3
shift= [0,0]
d = 2
n=500000

class ButterflySampler(BaseSpectralMeasureSampler):
    def sample(self, number_of_samples: int, random_state = None):
        p = np.random.rand(number_of_samples)
        theta = np.empty(number_of_samples)

        mask = p <= 0.5
        theta[mask] = np.random.uniform(-np.pi / 4, np.pi / 4, size=mask.sum())
        theta[~mask] = np.random.uniform(
            3 * np.pi / 4,
            5 * np.pi / 4,
            size=(~mask).sum()
        )

        x = np.cos(theta)
        y = np.sin(theta)

        return np.column_stack((x, y))

    def dimensions(self) -> int:
        return 2

    def mass(self) -> float:
        return 1.0

samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler=ButterflySampler(), size = n, shift=shift)

We highlight that the mass of the sphere and the dimension should be specified when implementing the sub-class. The resulting scatter plot is shown in Figure 9b. It is important to emphasize that when $\alpha \geq 1$, the custom spectral measure should satisfy the condition in equation (22).

One could mix a custom spectral measure, with another spectral measure. For example, we mix below equally the spectral measure defined above with the discrete spectral measure with $\text{masses}=[(1,0),(0,1),(-1,0),(0,-1)]$ and $\text{weights}=[1,1,0.25,0.25],$ and the random vector has $\alpha=0.25$, $\boldsymbol{\mu}^{0 T}=(0,0)$.

import aub_htp as ht
from aub_htp.random import BaseSpectralMeasureSampler, DiscreteSampler, MixedSampler
alpha = 0.25
shift= [0,0]
d = 2
n=1000000
class ButterflySampler(BaseSpectralMeasureSampler):
    def sample(self, number_of_samples: int, random_state = None):
        p = np.random.rand(number_of_samples)
        theta = np.empty(number_of_samples)

        mask = p <= 0.5
        theta[mask] = np.random.uniform(-np.pi / 4, np.pi / 4, size=mask.sum())
        theta[~mask] = np.random.uniform(
            3 * np.pi / 4,
            5 * np.pi / 4,
            size=(~mask).sum()
        )

        x = np.cos(theta)
        y = np.sin(theta)

        return np.column_stack((x, y))

    def dimensions(self) -> int:
        return 2

    def mass(self) -> float:
        return 1.0

samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler=ButterflySampler(), size = n, shift=shift)

# Discrete Spectral Measure
positions=[[1,0],[0,1],[-1,0],[0,-1]]
weights=[1,1,0.25,0.25]
sampler1 = DiscreteSampler(alpha = alpha, positions = positions, weights = weights)

#Spectral Measures and Their Associated Weights
spectral_measures=[sampler1, ButterflySampler()]
measure_weights=[0.5,0.5]

mixed_sampler=MixedSampler(spectral_measures=spectral_measures,weights=measure_weights)
samples = ht.multivariate_alpha_stable.rvs(alpha = alpha, spectral_measure_sampler = mixed_sampler,
shift = shift, size = n)
 

The resulting scatterplot is shown in Figure 9c.

6Conclusion¶

This paper presented AUB-HTP, a Python package for density computation and random variate generation of $\alpha$-stable distributions. The package is motivated by the practical difficulty of working with $\alpha$-stable laws, whose densities generally lack closed-form expressions and whose multivariate simulation requires careful handling of the underlying spectral measure.

AUB-HTP provides scalar density evaluation through a hybrid strategy that combines input standardization, Zolotarev-type integral representations, series formulas, and numerical inversion of characteristic functions. This design improves robustness across parameter regimes and enables systematic comparisons in terms of accuracy and runtime. The validation results show that the package provides reliable scalar density values and identify parameter regimes in which SciPy’s reference implementation appears to suffer from numerical inaccuracies.

The package also implements scalar and multivariate random variate generation. For multivariate $\alpha$-stable vectors, AUB-HTP uses LePage series representations and supports several spectral-measure models, including discrete, isotropic, sub-Gaussian/elliptical, mixed, and custom spectral measures. These features extend the range of $\alpha$-stable models that can be simulated within the Python ecosystem. The truncation-error analysis further provides guidance on the approximation induced by finite LePage-series representations.

By bringing together density evaluation, flexible simulation, validation, runtime analysis, and extensible implementation tools, AUB-HTP aims to make $\alpha$-stable modeling more accessible for computational work involving heavy-tailed data. Future developments may include multivariate PDF computation, CDF evaluation, parameter-estimation routines, further numerical acceleration, and application-oriented modules for communications, signal processing, and machine learning.

7Appendix: Error analysis for the LePage series estimate¶

In this appendix, we derive an upper bound on the mean square error of the LePage series estimate; we frequently use the following asymptotics of the gamma function from <i>NIST Digital Library of Mathematical Functions</i> (n.d.): Let $x\in\mathbb{R}$ and $a$ and $b$ be non-negative scalars. Since

then $\forall a,b \in\mathbb{R}$, there exists an $n_0$ such that

7.1The LePage series truncation error¶

The proof is as follows.

We start with the case $\alpha<1$. We have

which implies

where inequality $(a)$ is due to the Cauchy-Schwarz inequality. The equality in $(b)$ is due to the fact that $\Gamma_k\sim\text{Gamma}(k,1),n\geq \frac{2}{\alpha}$, and $\mathbb{E}\left[\Gamma_k^{-\frac{2}{\alpha}}\right]=\frac{\Gamma(k-\frac{2}{\alpha})}{\Gamma(k)}$, and inequality $(c)$ is due to (35). Inequality $(d)$ is due to the fact that

For the case $\alpha\geq1$, since $\|\mathbf{V}_k\|=1\;\forall \;k=1,\cdot\cdot\cdot,n$, equation (22) implies $\mathbb{E}[\mathbf{V}]=\mathbf{0}$. Let $\mathbf{V}_j=(v_j^{1},v_j^{2},\cdot\cdot\cdot,v_j^{d})$ then we have

where equality $(a)$ is due to the fact that $\mathbf{V}_k \perp\!\!\!\perp \Gamma_{k}$ and $(b)$ is justified since $v_k^{i} \perp\!\!\!\perp v^{i}_{\ell}$, $\forall \, k\neq\ell$ and $\mathbb{E}[v^{i}_k]=0$, $\forall \,k$. The rest of the equalities and inequalities follow by similar steps as in the case $\alpha<1$.

One can clearly see that the bound on the error decays faster for smaller values of $\alpha$. This indicates that this sampling method can appropriately handle the difficult region of small values of $\alpha$.

7.2Sub-Gaussian error¶

Sampling a multivariate sub-Gaussian $\alpha$-stable vector can be done by first computing the truncated series $\mathbf{X}_n$ for the isotropic spectral measure with unit-scale, and then applying $\mathbf{Y}_n=\Sigma^{\frac{1}{2}}\mathbf{X}_n$ where $\Sigma^{\frac{1}{2}}$ is the Cholesky decomposition of $\Sigma$– the desired shape matrix as described in Section 4.2.3. The resulting MSE in that case satisfies

where we used the fact that $\|\mathbf{x}\|^2_{\Sigma} \leq \lambda_{\max}(\Sigma)\|\mathbf{x}\|^2$, with $\lambda_{\max}(\Sigma)$ –the largest eigenvalue of $\Sigma$– being the norm of matrix $\Sigma$ Horn & Johnson (1985). Finally, one can upperbound $\mathbb{E} \left[ \|\mathbf{X}_n-\mathbf{X}\|^2 \right]$ using Theorem 2.