RatioUniforms#
- class scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)[source]#
Generate random samples from a probability density function using the ratio-of-uniforms method.
- Parameters:
- pdfcallable
A function with signature pdf(x) that is proportional to the probability density function of the distribution.
- umaxfloat
The upper bound of the bounding rectangle in the u-direction.
- vminfloat
The lower bound of the bounding rectangle in the v-direction.
- vmaxfloat
The upper bound of the bounding rectangle in the v-direction.
- cfloat, optional.
Shift parameter of ratio-of-uniforms method, see Notes. Default is 0.
- random_state{None, int,
numpy.random.Generator
, numpy.random.RandomState
}, optionalIf seed is None (or np.random), the
numpy.random.RandomState
singleton is used. If seed is an int, a newRandomState
instance is used, seeded with seed. If seed is already aGenerator
orRandomState
instance then that instance is used.
Methods
rvs
([size])Sampling of random variates
Notes
Given a univariate probability density function pdf and a constant c, define the set
A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}
. If(U, V)
is a random vector uniformly distributed overA
, thenV/U + c
follows a distribution according to pdf.The above result (see [1], [2]) can be used to sample random variables using only the PDF, i.e. no inversion of the CDF is required. Typical choices of c are zero or the mode of pdf. The set
A
is a subset of the rectangleR = [0, umax] x [vmin, vmax]
whereumax = sup sqrt(pdf(x))
vmin = inf (x - c) sqrt(pdf(x))
vmax = sup (x - c) sqrt(pdf(x))
In particular, these values are finite if pdf is bounded and
x**2 * pdf(x)
is bounded (i.e. subquadratic tails). One can generate(U, V)
uniformly onR
and returnV/U + c
if(U, V)
are also inA
which can be directly verified.The algorithm is not changed if one replaces pdf by k * pdf for any constant k > 0. Thus, it is often convenient to work with a function that is proportional to the probability density function by dropping unnecessary normalization factors.
Intuitively, the method works well if
A
fills up most of the enclosing rectangle such that the probability is high that(U, V)
lies inA
whenever it lies inR
as the number of required iterations becomes too large otherwise. To be more precise, note that the expected number of iterations to draw(U, V)
uniformly distributed onR
such that(U, V)
is also inA
is given by the ratioarea(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf)
, where area(pdf) is the integral of pdf (which is equal to one if the probability density function is used but can take on other values if a function proportional to the density is used). The equality holds since the area ofA
is equal to0.5 * area(pdf)
(Theorem 7.1 in [1]). If the sampling fails to generate a single random variate after 50000 iterations (i.e. not a single draw is inA
), an exception is raised.If the bounding rectangle is not correctly specified (i.e. if it does not contain
A
), the algorithm samples from a distribution different from the one given by pdf. It is therefore recommended to perform a test such askstest
as a check.References
[2]W. Hoermann and J. Leydold, “Generating generalized inverse Gaussian random variates”, Statistics and Computing, 24(4), p. 547–557, 2014.
[3]A.J. Kinderman and J.F. Monahan, “Computer Generation of Random Variables Using the Ratio of Uniform Deviates”, ACM Transactions on Mathematical Software, 3(3), p. 257–260, 1977.
Examples
>>> import numpy as np >>> from scipy import stats
>>> from scipy.stats.sampling import RatioUniforms >>> rng = np.random.default_rng()
Simulate normally distributed random variables. It is easy to compute the bounding rectangle explicitly in that case. For simplicity, we drop the normalization factor of the density.
>>> f = lambda x: np.exp(-x**2 / 2) >>> v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2) >>> umax = np.sqrt(f(0)) >>> gen = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng) >>> r = gen.rvs(size=2500)
The K-S test confirms that the random variates are indeed normally distributed (normality is not rejected at 5% significance level):
>>> stats.kstest(r, 'norm')[1] 0.250634764150542
The exponential distribution provides another example where the bounding rectangle can be determined explicitly.
>>> gen = RatioUniforms(lambda x: np.exp(-x), umax=1, vmin=0, ... vmax=2*np.exp(-1), random_state=rng) >>> r = gen.rvs(1000) >>> stats.kstest(r, 'expon')[1] 0.21121052054580314