Paper: Density Functions and Random Number Generators of Alpha-Stable Distributions

[Ahmad-ElHajj]

No description provided.

[github-actions]

Curvenote Preview

Directory	Preview	Checks	Updated (UTC)
papers/alphastable_distributions	🔍 Inspect	✅ 57 checks passed (22 optional)	Aug 14, 2026, 9:37 PM

[JennEYoon]

Reviewers, here are additional information.

Reviewers -- readme: https://github.com/scipy-conference/scipy_proceedings/blob/2026/README.md#instructions-for-reviewers

Reviewers -- How to add comments on GitHub Pull Request:
https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/reviewing-changes-in-pull-requests/commenting-on-a-pull-request

[zero-to-launch]

I accept as a reviewer.

[JennEYoon]

Hi Prithvish @prithvish-doshi-17
Thank you for accepting!

[peterhys]

I can review as well.

[zero-to-launch]

Hi @Ahmad-ElHajj , I tried rendering the main.tex as a pdf in overleaf, and saw multiple errors. It would be great if you can share the pdf for review, or have some instructions to render the pdf correctly. Thank you.

[Ahmad-ElHajj]

@prithvish-doshi-17 Thank you for your interest in this work.

You're right. Latex intended for MyST will not compile on overleaf, and vice versa. We actually started on overleaf and had to make the necessary modifications to make MyST compile the latex to html.

The latex-to-typst (and subsequently latex-to-pdf) compiler seems to be even stricter. To fix it, I made some changes to the paper's .tex and added a plugin to convert case statements (\begin{cases}\end{cases}) into their typst equivalent.
Locally, $ myst build --pdf should work fine now. And here's the PDF of the preprint: https://arxiv.org/pdf/2606.28530

[Ahmad-ElHajj]

added a plugin to convert case statements (\begin{cases}\end{cases}) into their typst equivalent.

@JennEYoon It seems like GitHub Actions detected the plugin and is refusing to compile the paper for security reasons. Without this plugin, cases are silently ignored and do not get rendered on the PDF. Could you please tag the technical team for support?

[peterhys]

Hi @Ahmad-ElHajj @JennEYoon

Here is my review:

The paper presents a substantial package for alpha-stable density computation and the addition of LePage multivariate implementation and spectral measure samplers. I think this is nice work, and it could benefit the scientific Python community. The paper is well written. Some revision is needed to clarify the statements in the paper.

Comments:

1. I think the paper can benefit from discussing some other existing implementations in the introduction. In particular, I think there are several good implementations of alpha-stable distributions in R that the paper should cite. For example, libstable (by Royuela-del-Val, et al. doi.org/10.18637/jss.v078.i01), which also provides whole parameter space density computation. And maybe Swihart & Nolan, 2022 (10.32614/RJ-2022-056), alphastable Teimouri et al. (arxiv:1809.09546)
2. A key claim of the paper is that the computed density remains smooth compared to SciPy. However, a brief inspection of the code suggests boundary padding and monotonicity smoothing are used (pad_grid, remove_all_monotonicity_spikes and remove_left_monotonicity_spikes in "pdf/pdf.py"). SciPy does not perform curve-level smoothing. I think these operations should be explicitly mentioned in the paper. I think the authors should show the raw output comparison between the package and SciPy for Figure 2 instead of the smoothed version.
3. Smoothing can introduce bias, and make results grid-dependent (it looks like a fixed 10,000-point grid is used for interpolation get_alpha_stable_pdf_interpolator). A comparison with and without the smoothing needs to be shown, along with an error comparison of the smoothing effect at different grid sizes.
4. The runtime comparison should provide a detailed testing environment, versions, and precision. The reported speedup may depend on cache state in the 10,000-point interpolator get_alpha_stable_pdf_interpolator. I think that should be explicitly mentioned in the paper as well, in section 3.5.
5. The fixed 0.01 MSE target should be justified in section 5. And talk about how the number of LePage terms is determined.

Minor comments:

1. Figure 2 is too small to see the details.
2. The paper seems to exceed the 8-page limit that SciPy Proceedings requires. I am not sure if some of the general descriptions in sections 2 and 3 are necessary, and also some of the examples in section 5.
3. The code snippets are not consistent with the space before and after the equal sign. Consider wrapping the code snippets to 80 - 90 characters.
4. Eq. 31 line 4 should be E[||U||^{2 \alpha}]

[zero-to-launch]

Hi @Ahmad-ElHajj , here are my review comments:

The paper presents AUB-HTP, a Python package for computing scalar α-stable probability densities and generating univariate and multivariate α-stable random variables. Its strongest contribution is the practical combination of several established numerical methods: characteristic-function inversion near the mode, Zolotarev integration in intermediate regions, and asymptotic series in the tails. It also extends beyond standard univariate functionality by supporting multivariate simulation through LePage series with isotropic, elliptical, discrete, mixed, and custom spectral measures.

Compared with existing tools such as SciPy’s levy_stable and older stable-distribution packages, the main improvement appears to be usability and broader multivariate flexibility rather than entirely new underlying mathematics. The adaptive selection between numerical methods is useful and the reported runtime improvement is promising, especially for repeated or vectorized density evaluations. The custom spectral-measure interface is also valuable because general multivariate α-stable simulation is less accessible in the Python ecosystem than standard univariate sampling.

My main concern is the validation strategy. SciPy is used both to help determine the switching cutoffs and as the main numerical reference, while disagreements near α = 1 are attributed to weaknesses in SciPy. The authors should add at least one independent validation source, such as analytical special cases, high-precision characteristic-function inversion, or comparison with another established stable-distribution implementation. Similarly, the multivariate examples are visually convincing, but characteristic-function or projection-based tests would provide stronger evidence that the generated samples follow the intended distributions.

Overall, this is a useful and well-motivated software contribution with clear practical relevance. I would recommend acceptance after revision, mainly to strengthen the independent numerical validation, clarify the comparison with existing software, and more precisely state that the novelty lies in the hybrid implementation and flexible Python interface rather than in the classical formulas themselves.

[JennEYoon]

added a plugin to convert case statements (\begin{cases}\end{cases}) into their typst equivalent.

@JennEYoon It seems like GitHub Actions detected the plugin and is refusing to compile the paper for security reasons. Without this plugin, cases are silently ignored and do not get rendered on the PDF. Could you please tag the technical team for support?

Hi Ahmad @Ahmad-ElHajj,
Sorry I took a few days off last week. :-)

Franklin (@fwkoch ) as of last Wednesday (7/22nd) is going to review all the papers and fix compile errors if he can. I will give him till this Wednesday (7/29), and raise it with him. The preprint PDF looks great. Right now, the pretty view is for reviewers, so your preprint works well.

OK I've messaged Franklin. He will respond when he can. :-)
Jennifer, edited 7/29 Wednesday.

[JennEYoon]

Dear Ahmad @Ahmad-ElHajj,

Peter said the paper is longer than the 8 pages max. I can help you reduce it later. I will also do a final pass on all figures.

[JennEYoon]

Dear Peter @peterhys and Prithvish @zero-to-launch,
Thank you for your excellent reviews!! These are truly in the Spirit of SciPy Proceedings.
Jennifer

[Ahmad-ElHajj]

Thank you for your reviews, @peterhys and @zero-to-launch . Your efforts are truly appreciated! I'm currently coordinating with the team to formulate a modification plan and a comprehensive response.

[Ahmad-ElHajj]

Dear Peter @peterhys,

Thank you for your detailed and thorough review.

1- Thank you for pointing out these references. As requested, we will include them in the updated manuscript as part of the previous work discussion and clarify the relationship with our implementation.

Regarding the density computation, we note that:

• "libstable" covers the whole parameter space like ours, but evaluates Nolan's/Zolotarev's integral by Gauss–Kronrod quadrature everywhere, in C/C++ with a MATLAB front-end in addition to an R package. Our method reserves quadrature for the near-mode and middle regions and switches to Skorohod–Pollard–Bergström series in the tails, using (α, β)-calibrated cutoffs, all natively in Python with a SciPy compatible interface.
• As for "Swihart & Nolan", it is primarily a multivariate implementation. For the univariate PDF it uses libstableR, so it inherits libstable's quadrature approach rather than introducing a new generator.
• The "alphastable" package seems to be using the same Skorohod series as our implementation. It is closest to our tail strategy. It uses convergent asymptotic series and falls back to stabledist elsewhere. We differ by adding a third regime (direct characteristic-function inversion near the mode) and providing a SciPy-compatible Python interface.

When it comes to random number generation, we note that:

• the "libstable" package of Royuela-del-Val et al. presents a univariate RNG using the same modified CMS method as us. Moreover, they cite work by Weron 2004, which mentions that this method is the fastest and most accurate. The same paper mentions the LePage series as a simulation method and references for it, but does not implement the series nor analyze it. While the reference for LePage series mentioned in Weron, provides different bounds on the LePage series MSE, our bounds are more practical for our study. It is nevertheless important to reference these bounds in our work, which we will.
• "Swihart & Nolan" (2022) present both a univariate stable random number generator and a multivariate random vector generator taken from Teimouri et al., whose work you also mentioned.
• "Teimouri et al."" generate univariate α-stable random variables using the Chambers–Mallows–Stuck (CMS) method, as we do. Additionally, they present a method to simulate from a truncated stable distribution. For the multivariate case, they describe the generation of sub-Gaussian stable vectors using the definition given in the book by Samorodnitsky and Taqqu; this consists of sampling a Gaussian random vector and multiplying it by the square root of a skewed alpha stable random variable. Moreover, to generate a bivariate stable random vector they implement a method suggested in "R. Modarres and J. P. Nolan: A method for simulating stable random vectors." In contrast, our AUB-HTP package can simulate any spectral measure for a wide range of values of $\alpha$ and for any dimension $d$. On a final note, we will implement the subGaussian RNG using the definition given in Samorodnitsky and Taqqu which is adopted by "Teimouri et al.", and compare it to our own implementation of the subGaussian RNG which is based on the LePage series.

We will revise the manuscript to discuss these references and to clarify how the scope and implementation of our package differs and extends the available methods.

2- Thank you for pointing that out, we will explicitly mention smoothing. In the examples presented in Figure 2, we note that even after we disable the smoothing procedure, we achieve the same results and our generator remains smooth in those cases. Please refer to the newly generated and attached figure. We added the smoothing method to maintain the monotonicity of the PDF whenever numerical irregularities violate the monotonicity property. We note that through our extensive set of tests, we observed minor spikes happening away from the mode in the PDF generator. The main difference between the spikes we were observing and the irregularities of Scipy is that in Scipy, irregularities happen near the mode and the PDF is set to 0 on a continuous range whereas our PDF generator may result in a few spikes on disjoint points that can be fixed by smoothing. Finally, we believe that the presented results should be kept with the smoothed version in the manuscript since our final package, which is available to the user, presents a smoothing procedure after computing the PDF.

3 & 4- Although smoothing is part of the officially released package (currently at v1.0.9), the interpolation you're referring to is for unreleased trials we've been experimenting with for the purpose of speeding up computations of some statistics that rely on multiple computations of the PDF. We may end up scrapping these, or moving them to the (unreleased/in-dev) aub_htp.statistics submodule. To avoid confusion, we've added a notice on the development branch that points you to the latest release's branch https://github.com/AUB-HTP/AUB-HTP/tree/maintenance/v1.0. We'd be happy to rerun the tests with the experimental code for your review. However, the testing and benchmarking results displayed in the paper are the code's performance without the experimental interpolation, as they were run on the latest release's aub_htp.pdf.

5- Naturally our objective is to choose a reasonably small value for the default error while keeping the required number of terms n not excessively large for efficient running time and memory. After empirically varying the parameters, $\alpha, \text{MSE}$, and $\Lambda(\mathbb{S}^{d-1})$ and checking the required number of samples for each combination, we elected to choose a value for the error of 0.01. On the other side, please note that we do allow users to choose their own error as per https://github.com/AUB-HTP/AUB-HTP/blob/43d55b5cdd186d1fa5fd4b931a8ecfc3c01a2287/aub_htp/random/alpha_stable_sampler.py#L13. This parameter is not bubbled up to the main aub_htp.alpha_stable interface so we'll make sure it is in the revision. Internally, we infer the number of terms n from theorem 2's inequality:

M1- Thank you for pointing that out. We will update the figure.
M2- Our understanding of the proceedings guidelines we found at "https://github.com/scipy-conference/scipy_proceedings#:~:text=The%20paper%20should%20be%20at%20most%206000%20words%20including%20figures%20but%20not%20including%20references%3B%20this%20is%20about%208%20pages%20for%20the%20published%20PDF"
is that the paper should be at most 6000 words including figures but not including references; this is about 8 pages for the published PDF.
Although in our case the paper does exceed 8 pages, it is still compliant with the proceedings' guidelines as it does not exceed the actual limit of 6000 words. As suggested, we will also carefully revise Sections 2 and 3, as well as the examples in Section 5, to remove unnecessary general descriptions and repetitive statements. We will aim to shorten the manuscript while preserving the clarity of the presentation and all relevant technical details and reported results.
M3- We will run the code displayed in the paper through an auto-formatter to ensure all of the formatting issues are resolved.
M4- This has been fixed. Thank you.

[Ahmad-ElHajj]

Dear Prithvish @zero-to-launch,

Your review is highly appreciated.

1- For clarification, we confirm that Scipy's PDF generator has been compared against Nolan's STABLE program (MATLAB implementation) across the (α, β) parameter space. We were able to confirm that the two generators agree everywhere except on the set of points where Scipy itself fails. Then, we verified that those points fall outside the regions where our method uses Scipy to determine the cutoffs. Therefore, our cutoffs and boundary points are not affected by Scipy's irregularities since those points were chosen in regimes where Scipy agrees with the STABLE program. In summary, Scipy was used as a reference after we made sure it agrees with STABLE on most of the parameter space. Additionally, our density generator was indeed validated against the analytically known special cases: Gaussian, Levy and Cauch distributions.

2- Our main intent is for the plots to show the effect of changing the spectral measure; however and as pointed out, they do not provide a lot of information about their correctness. For validation we use a different tool in the manuscript and show that the mean square error (MSE) of the samples goes to 0 as $n\to\infty$. For additional RNG validation, we conducted several projection-based tests using the Kolmogorov-Smirnov (KS) statistical test with respect to Nolan's STABLE Matlab package. The AUB-HTP software was configured to generate samples with a numerical error below 0.001. In all the studied cases (isotropic and discrete spectral measures), the KS D-statistics were considerably small, indicating excellent agreement with the target distributions. The resulting p-values were large as well with the exception of few cases. This result should be interpreted with caution, since the KS test becomes highly sensitive at large sample sizes which may result in low p-values even when the observed discrepancy is practically negligible. Therefore, the small D-statistic still indicates a close fit between the generated and target distributions. We plan on highlighting some of these quantitative results in the updated version of the paper which we believe will enhance its quality.

[JennEYoon]

Dear Ahmad @Ahmad-ElHajj,
Great to see that everything is going well with the review process. A quick note from me.

You are the only author without an orcid on this paper. Somewhere on the proceedings repo, there are instructions for creating one, also see https://orcid.org/. Proceedings will make you create one before submitting this paper. And you can add this paper as your very first publication on your orcid registry! :-) Your virtual poster will also need the same orcid, and you can add that as your 2nd item!

I see you’ve just uploaded your edits to the paper. Can you post the pretty view preprint for me? I don’t have a good .tex viewer right now. Eventually Franklin will fix the Curvenote preview rendering problem. If this has already been solved, let me know where or how to view it.

Jennifer

[Ahmad-ElHajj]

Hey Jennifer @JennEYoon,
I'll look into how to create an ORCID. Thank you for pointing it out!
We still haven't updated the arXiv preprint. So I'm attaching the built PDF using the proceeding's tool as of this commit. Note that since the proceedings did not yet update it, it is using the 2025 template. Otherwise, the content is fine.
full_text.pdf

For anyone else viewing this, please note this PDF is invalidated as soon as a new commit is pushed.