Automatic random variate generation in Python

—The generation of random variates is an important tool that is required in many applications. Various software programs or packages contain generators for standard distributions like the normal, exponential or Gamma, e.g., the programming language R and the packages SciPy and NumPy in Python. However, it is not uncommon that sampling from new/non-standard distributions is required. Instead of deriving speciﬁc generators in such situations, so-called automatic or black-box methods have been developed. These allow the user to generate random variates from fairly large classes of distributions by only specifying some properties of the distributions (e.g. the density and/or cumulative distribution function). In this note, we describe the implementation of such methods from the C library UNU.RAN in the Python package SciPy and provide a brief overview of the functionality.


Introduction
The generation of random variates is an important tool that is required in many applications.Various software programs or packages contain generators for standard distributions, e.g., R ([R C21]) and SciPy ([VGO + 20]) and NumPy ([HMvdW + 20]) in Python.Standard references for these algorithms are the books [Dev86], [Dag88], [Gen03], and [Knu14].An interested reader will find many references to the vast existing literature in these works.While relying on general methods such as the rejection principle, the algorithms for well-known distributions are often specifically designed for a particular distribution.This is also the case in the module stats in SciPy that contains more than 100 distributions and the module random in NumPy with more than 30 distributions.However, there are also so-called automatic or black-box methods for sampling from large classes of distributions with a single piece of code.For such algorithms, information about the distribution such as the density, potentially together with its derivative, the cumulative distribution function (CDF), and/or the mode must be provided.See [HLD04] for a comprehensive overview of these methods.Although the development of such methods was originally motivated to generate variates from nonstandard distributions, these universal methods have advantages that make their usage attractive even for sampling from standard distributions.We mention some of the important properties (see [LH00], [HLD04], [DHL10]): • The algorithms can be used to sample from truncated distributions.
• For inversion methods, the structural properties of the underlying uniform random number generator are preserved and the numerical accuracy of the methods can be controlled by a parameter.Therefore, inversion is usually the only method applied for simulations using quasi-Monte Carlo (QMC) methods.
• Depending on the use case, one can choose between a fast setup with slow marginal generation time and vice versa.
The latter point is important depending on the use case: if a large number of samples is required for a given distribution with fixed shape parameters, a slower setup that only has to be run once can be accepted if the marginal generation times are low.If small to moderate samples sizes are required for many different shape parameters, then it is important to have a fast setup.The former situation is referred to as the fixed-parameter case and the latter as the varying parameter case.
Implementations of various methods are available in the C library UNU.RAN ( [HL07]) and in the associated R package Runuran (https://cran.r-project.org/web/packages/Runuran/index.html,[TL03]).The aim of this note is to introduce the Python implementation in the SciPy package that makes some of the key methods in UNU.RAN available to Python users in SciPy 1.8.0.These general tools can be seen as a complement to the existing specific sampling methods: they might lead to better performance in specific situations compared to the existing generators, e.g., if a very large number of samples are required for a fixed parameter of a distribution or if the implemented sampling method relies on a slow default that is based on numerical inversion of the CDF.For advanced users, they also offer various options that allow to fine-tune the generators (e.g., to control the time needed for the setup step).

Automatic algorithms in SciPy
Many of the automatic algorithms described in [HLD04] and [DHL10] are implemented in the ANSI C library, UNU.RAN (Universal Non-Uniform RANdom variate generators).Our goal was to provide a Python interface to the most important methods from UNU.RAN to generate univariate discrete and continuous non-uniform random variates.The following generators have been implemented in SciPy 1.8.0:Before describing the implementation in SciPy in Section scipy_impl, we give a short introduction to random variate generation in Section intro_rv_gen.

A very brief introduction to random variate generation
It is well-known that random variates can be generated by inversion of the CDF F of a distribution: if U is a uniform random number on (0, 1), X := F −1 (U) is distributed according to F. Unfortunately, the inverse CDF can only be expressed in closed form for very few distributions, e.g., the exponential or Cauchy distribution.If this is not the case, one needs to rely on implementations of special functions to compute the inverse CDF for standard distributions like the normal, Gamma or beta distributions or numerical methods for inverting the CDF are required.Such procedures, however, have the disadvantage that they may be slow or inaccurate, and developing fast and robust inversion algorithms such as HINV and PINV is a non-trivial task.HINV relies on Hermite interpolation of the inverse CDF and requires the CDF and PDF as an input.PINV only requires the PDF.The algorithm then computes the CDF via adaptive Gauss-Lobatto integration and an approximation of the inverse CDF using Newton's polynomial interpolation.Note that an approximation of the inverse CDF can be achieved by interpolating the points (F(x i ), x i ) for points x i in the domain of F, i.e., no evaluation of the inverse CDF is required.
For discrete distributions, F is a step-function.To compute the inverse CDF F −1 (U), the simplest idea would be to apply sequential search: if X takes values 0, 1, 2, . . .with probabilities p 0 , p 1 , p 2 , . . ., start with j = 0 and keep incrementing j until F( j) = p 0 + • • • + p j ≥ U. When the search terminates, X = j = F −1 (U).Clearly, this approach is generally very slow and more efficient methods have been developed: if X takes L distinct values, DGT realizes very fast inversion using so-called guide tables / hash tables to find the index j.In contrast DAU is not an inversion method but uses the alias method, i.e., tables are precomputed to write X as an equi-probable mixture of L twopoint distributions (the alias values).
The rejection method has been suggested in [VN51].In its simplest form, assume that f is a bounded density on Note that the accepted points (U,V ) are uniformly distributed in the region between the x-axis and the graph of the PDF.Hence, X := V has the desired distribution f .This is a special case of the general version: if f , g are two densities on an interval J such that f (x) ≤ c • g(x) for all x ∈ J and a constant c ≥ 1, sample U uniformly distributed on [0, 1] and X distributed according to g until c • U • g(X) ≤ f (X).Then X has the desired distribution f .It can be shown that the expected number of iterations before the acceptance condition is met is equal to c. Hence, the main challenge is to find hat functions g for which c is small and from which random variates can be generated efficiently.TDR solves this problem by applying a transformation T to the density such that x → T ( f (x)) is concave.A hat function can then be found by computing tangents at suitable design points.Note that by its nature any rejection method requires not always the same number of uniform variates to generate one non-uniform variate; this makes the use of QMC and of some variance reduction methods more difficult or impossible.On the other hand, rejection is often the fastest choice for the varying parameter case.
The Ratio-Of-Uniforms method (ROU, [KM77]) is another general method that relies on rejection.The underlying principle is that if (U,V ) is uniformly distributed on the set where f is a PDF with support (a, b), then X := U/V follows a distribution according to f .In general, it is not possible to sample uniform values on A f directly.However, if for finite constants u − , u + , v + , one can apply the rejection method: generate uniform values (U,V ) on the bounding rectangle R until (U,V ) ∈ A f and return X = U/V .Automatic methods relying on the ROU method such as SROU and automatic ROU ([Ley00]) need a setup step to find a suitable region S ∈ R 2 such that A f ⊂ S and such that one can generate (U,V ) uniformly on S efficiently.

Description of the SciPy interface
SciPy provides an object-oriented API to UNU.RAN's methods.To initialize a generator, two steps are required: 1) creating a distribution class and object, 2) initializing the generator itself.
In step 1, a distributions object must be created that implements required methods (e.g., pdf, cdf).This can either be a custom object or a distribution object from the classes rv_continuous or rv_discrete in SciPy.Once the generator is initialized from the distribution object, it provides a rvs method to sample random variates from the given distribution.It also provides a ppf method that approximates the inverse CDF if the initialized generator uses an inversion method.The following example illustrates how to initialize the NumericalInversePolynomial (PINV) generator for the standard normal distribution: As NumericalInversePolynomial generator uses an inversion method, it also provides a ppf method that approximates the inverse CDF: # evaluate the approximate PPF at a few points ppf = rng.ppf([0.1,0.5, 0.9]) It is also easy to sample from a truncated distribution by passing a domain argument to the constructor of the generator.For example, to sample from truncated normal distribution: # truncate the distribution by passing a # `domain`argument rng = sampling.NumericalInversePolynomial( dist, domain=(-1, 1) ) While the default options of the generators should work well in many situations, we point out that there are various parameters that the user can modify, e.g., to provide further information about the distribution (such as mode or center) or to control the numerical accuracy of the approximated PPF.(u_resolution).Details can be found in the SciPy documentation https://docs.scipy.org/doc/scipy/reference/.The above code can easily be generalized to sample from parametrized distributions using instance attributes in the distribution class.For example, to sample from the gamma distribution with shape parameter alpha, we can create the distribution class with parameters as instance attributes: In the above example, the support method is used to set the domain of the distribution.This can alternatively be done by passing a domain parameter to the constructor.In addition to continuous distribution, two UNU.RAN methods have been added in SciPy to sample from discrete distributions.In this case, the distribution can be either be represented using a probability vector (which is passed to the constructor as a Python list or NumPy array) or a Python object with the implementation of the probability mass function.In the latter case, a finite domain must be passed to the constructor or the object should implement the support method 1 .

Underlying uniform pseudo-random number generators
NumPy provides several generators for uniform pseudo-random numbers 2 .It is highly recommended to use NumPy's default random number generator np.random.PCG64 for better speed and performance, see [O'N14] and https://numpy.org/doc/stable/ 1. Support for discrete distributions with infinite domain hasn't been added yet.
reference/random/bit_generators/index.html.To change the uniform random number generator, a random_state parameter can be passed as shown in the example below: # 64-bit PCG random number generator in NumPy urng = np.random.Generator(np.random.PCG64()) # The above line can also be replaced by: # ``urng = np.random.default_rng()`#as PCG64 is the default generator starting # from NumPy 1.19.0 # change the uniform random number generator by # passing the `random_state`argument rng = sampling.NumericalInversePolynomial( dist, random_state=urng ) We also point out that the PPF of inversion methods can be applied to sequences of quasi-random numbers.SciPy provides different sequences in its QMC module (scipy.stats.qmc).
NumericalInverseHermite provides a qrvs method which generates random variates using QMC methods present in SciPy (scipy.stats.qmc)as uniform random number generators 3 .The next example illustrates how to use qrvs with a generator created directly from a SciPy distribution object.# generate quasi random numbers using the Halton # sequence as uniform variates qrvs = rng.qrvs(size=100,qmc_engine=qrng)

Benchmarking
To analyze the performance of the implementation, we tested the methods applied to several standard distributions against the generators in NumPy and the original UNU.RAN C library.In addition, we selected one non-standard distribution to demonstrate that substantial reductions in the runtime can be achieved compared to other implementations.All the benchmarks were carried out using NumPy 1.22.4 and SciPy 1.8.1 running in a single core on Ubuntu 20.04.3 LTS with Intel(R) Core(TM) i7-8750H CPU (2.20GHz clock speed, 16GB RAM).We run the benchmarks with NumPy's MT19937 (Mersenne Twister) and PCG64 random number generators (np.random.MT19937 and np.random.PCG64) in Python and use NumPy's C implementation of MT19937 in the UNU.RAN C benchmarks.As explained above, the use of PCG64 is recommended, and MT19937 is only included to compare the speed of the Python implementation and the C library by relying on the same uniform number generator (i.e., differences in the performance of the uniform number generation are not taken into account).The code for all the benchmarks can be found on https://github.com/tirthasheshpatel/unuran_benchmarks.
The methods used in NumPy to generate normal, gamma, and beta random variates are: • the ziggurat algorithm ( [MT00b]) to sample from the standard normal distribution, 2. By default, NumPy's legacy random number generator, MT19937 (np.random.RandomState()) is used as the uniform random number generator for consistency with the stats module in SciPy.
Benchmarking against the normal, gamma, and beta distributions Table 1 compares the performance for the standard normal, Gamma and beta distributions.We recall that the density of the Gamma distribution with shape parameter a > 0 is given by x ∈ (0, ∞) → x a−1 e −x and the density of the beta distribution with shape parameters α, β > 0 is given by x ∈ (0, 1) where Γ(•) and B(•, •) are the Gamma and beta functions.The results are reported in Table 1.
We summarize our main observations: 1) The setup step in Python is substantially slower than in C due to expensive Python callbacks, especially for PINV and HINV.However, the time taken for the setup is low compared to the sampling time if large samples are drawn.Note that as expected, SROU has a very fast setup such that this method is suitable for the varying parameter case.
2) The sampling time in Python is slightly higher than in C for the MT19937 random number generator.If the recommended PCG64 generator is used, the sampling time in Python is slightly lower.The only exception is SROU: due to Python callbacks, the performance is substantially slower than in C.However, as the main advantage of SROU is the fast setup time, the main use case is the varying parameter case (i.e., the method is not supposed to be used to generate large samples).3) PINV, HINV, and TDR are at most about 2x slower than the specialized NumPy implementation for the normal distribution.For the Gamma and beta distribution, they even perform better for some of the chosen shape parameters.These results underline the strong performance of these black-box approaches even for standard distributions.4) While the application of PINV requires bounded densities, no issues are encountered for α = 0.05 since the unbounded part is cut off by the algorithm.However, the setup can fail for very small values of α.
Benchmarking against a non-standard distribution We benchmark the performance of PINV to sample from the generalized normal distribution ([Sub23]) whose density is given by x ∈ (−∞, ∞) → pe −|x| p 2Γ(1/p) against the method proposed in [NP09] and against the implementation in SciPy's gennorm distribution.The approach in [NP09] relies on transforming Gamma variates to the generalized normal distribution whereas SciPy relies on computing the inverse of CDF of the Gamma distribution (https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammainccinv.html).The results for different values of p are shown in Table 2.
PINV is usually about twice as fast than the specialized method and about 15-150 times faster than SciPy's implementation 4 .We also found an R package pgnorm (https: //cran.r-project.org/web/packages/pgnorm/) that implements various approaches from [KR13].In that case, PINV is usually about 70-200 times faster.This clearly shows the benefit of using a black-box algorithm.

Conclusion
The interface to UNU.RAN in SciPy provides easy access to different algorithms for non-uniform variate generation for large classes of univariate continuous and discrete distributions.We have shown that the methods are easy to use and that the algorithms perform very well both for standard and non-standard Average time taken (reported in milliseconds, unless mentioned otherwise) to sample 1 million random variates from the standard normal distribution.The mean is computed over 7 iterations.Standard deviations are not reported as they were very small (less than 1% of the mean in the large majority of cases).Note that not all methods can always be applied, e.g., TDR cannot be applied to the Gamma distribution if a < 1 since the PDF is not log-concave in that case.As NumPy uses rejection algorithms with precomputed constants, no setup time is reported.

TABLE 1
distributions.A comprehensive documentation suite, a tutorial and many examples are available at https://docs.scipy.org/doc/scipy/reference/stats.sampling.htmland https://docs.scipy.org/doc/scipy/tutorial/stats/sampling.html.Various methods have been implemented in SciPy, and if specific use cases require additional functionality from UNU.RAN, the methods can easily be added to SciPy given the flexible framework that has been developed.Another area of further development is to better integrate SciPy's QMC generators for the inversion methods.Finally, we point out that other sampling methods like Markov Chain Monte Carlo and copula methods are not part of SciPy.Relevant Python packages in that context are PyMC ([PHF10]), PyStan relying on Stan ([Tea21]), Copulas (https://sdv.dev/Copulas/)and PyCopula (https://blent-ai.github.io/pycopula/).

TABLE 2
Comparing SciPy's implementation and a specialized method against PINV to sample 1 million variates from the generalized normal distribution for different values of the parameter p.Time reported in milliseconds.The mean is computer over 7 iterations.