Hi all,

In `stats.sampling`, various methods are implemented to generate random variates based on a specified distribution (the methods are from the UNURAN C library, Tirth implemented this as part of GSoC 2022). The tools are mainly for advanced users and are not very easy to use. In a follow-up project (Numfocus SDG), I worked with one of the UNURAN authors to make these features more easily accessible: we created a config for a large number of the SciPy distributions (36 to be precise) such that `NumericalInversePolynomial` can be directly used w/o having to care about the settings. For example, to sample from the Gamma distribution, one can just write

rng = stats.sampling.FastGeneratorInversion('gamma', (1.3,))
rng.rvs(size=10)

Further supported features:

- a fast inverse CDF `ppf_fast` that is computed when the object is created
- sampling with quasi-random numbers from `stats.qmc` (`qrvs`)
- truncation: all distributions can be truncated/restricted to an interval using the keyword `domain` (so in addition to the 36 distributions, a few more like truncnorm and truncweibull_min are covered)
- relying on the usual distribution infrastructure, methods like `pdf`, `cdf` are added to the distribution object (this is important in case truncation is used, the existing infrastructure does not support the evaluation of the pdf, cdf)

For a lot of distributions, the performance of rvs is better than that of `rv_continuous` if large samples are required (for small samples, the cost of the setup step to derive `ppf_fast` is too high compared to the marginal sampling times). Below are some examples (time in ms, generation of 10**7 rvs, statistics are computed over 10 trials). For the new method, one has to take into account the setup time + the marginal sampling time:

----- beta with parameters (2.5, 3.9) -----
New method setup speed (min., mean, std): 31.3, 48.2, 17.5
New method sampling speed (min., mean, std): 469.2, 540.2, 52.1
SciPy method sampling speed (min., mean, std): 1305.7, 1467.1, 137.6

----- crystalball with parameters (0.7, 1.2) -----
New method setup speed (min., mean, std): 66.5, 88.7, 19.8
New method sampling speed (min., mean, std): 399.6, 451.8, 45.2
SciPy method sampling speed (min., mean, std): 2515.3, 2809.9, 236.1

----- gennorm with parameters (2.7,) -----
New method setup speed (min., mean, std): 12.0, 16.8, 2.0
New method sampling speed (min., mean, std): 395.4, 414.5, 17.0
SciPy method sampling speed (min., mean, std): 1798.0, 1933.1, 98.5

The functionality is not supposed to replace the existing `rv_continuous` infrastructure: the main idea is to make the UNURAN functionality more easily available and to allow truncation and faster qrvs sampling for a large number of SciPy distributions.

I also want to thank Pamphile for providing helpful advice during the development.

Best regards,

Christoph