Coverage for /usr/lib/python3/dist-packages/scipy/stats/_rvs_sampling.py: 11%
27 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
2from scipy._lib._util import check_random_state
5def rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=1, c=0, random_state=None):
6 """
7 Generate random samples from a probability density function using the
8 ratio-of-uniforms method.
10 Parameters
11 ----------
12 pdf : callable
13 A function with signature `pdf(x)` that is proportional to the
14 probability density function of the distribution.
15 umax : float
16 The upper bound of the bounding rectangle in the u-direction.
17 vmin : float
18 The lower bound of the bounding rectangle in the v-direction.
19 vmax : float
20 The upper bound of the bounding rectangle in the v-direction.
21 size : int or tuple of ints, optional
22 Defining number of random variates (default is 1).
23 c : float, optional.
24 Shift parameter of ratio-of-uniforms method, see Notes. Default is 0.
25 random_state : {None, int, `numpy.random.Generator`,
26 `numpy.random.RandomState`}, optional
28 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
29 singleton is used.
30 If `seed` is an int, a new ``RandomState`` instance is used,
31 seeded with `seed`.
32 If `seed` is already a ``Generator`` or ``RandomState`` instance then
33 that instance is used.
35 Returns
36 -------
37 rvs : ndarray
38 The random variates distributed according to the probability
39 distribution defined by the pdf.
41 Notes
42 -----
43 Given a univariate probability density function `pdf` and a constant `c`,
44 define the set ``A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}``.
45 If `(U, V)` is a random vector uniformly distributed over `A`,
46 then `V/U + c` follows a distribution according to `pdf`.
48 The above result (see [1]_, [2]_) can be used to sample random variables
49 using only the pdf, i.e. no inversion of the cdf is required. Typical
50 choices of `c` are zero or the mode of `pdf`. The set `A` is a subset of
51 the rectangle ``R = [0, umax] x [vmin, vmax]`` where
53 - ``umax = sup sqrt(pdf(x))``
54 - ``vmin = inf (x - c) sqrt(pdf(x))``
55 - ``vmax = sup (x - c) sqrt(pdf(x))``
57 In particular, these values are finite if `pdf` is bounded and
58 ``x**2 * pdf(x)`` is bounded (i.e. subquadratic tails).
59 One can generate `(U, V)` uniformly on `R` and return
60 `V/U + c` if `(U, V)` are also in `A` which can be directly
61 verified.
63 The algorithm is not changed if one replaces `pdf` by k * `pdf` for any
64 constant k > 0. Thus, it is often convenient to work with a function
65 that is proportional to the probability density function by dropping
66 unnecessary normalization factors.
68 Intuitively, the method works well if `A` fills up most of the
69 enclosing rectangle such that the probability is high that `(U, V)`
70 lies in `A` whenever it lies in `R` as the number of required
71 iterations becomes too large otherwise. To be more precise, note that
72 the expected number of iterations to draw `(U, V)` uniformly
73 distributed on `R` such that `(U, V)` is also in `A` is given by
74 the ratio ``area(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf)``,
75 where `area(pdf)` is the integral of `pdf` (which is equal to one if the
76 probability density function is used but can take on other values if a
77 function proportional to the density is used). The equality holds since
78 the area of `A` is equal to 0.5 * area(pdf) (Theorem 7.1 in [1]_).
79 If the sampling fails to generate a single random variate after 50000
80 iterations (i.e. not a single draw is in `A`), an exception is raised.
82 If the bounding rectangle is not correctly specified (i.e. if it does not
83 contain `A`), the algorithm samples from a distribution different from
84 the one given by `pdf`. It is therefore recommended to perform a
85 test such as `~scipy.stats.kstest` as a check.
87 References
88 ----------
89 .. [1] L. Devroye, "Non-Uniform Random Variate Generation",
90 Springer-Verlag, 1986.
92 .. [2] W. Hoermann and J. Leydold, "Generating generalized inverse Gaussian
93 random variates", Statistics and Computing, 24(4), p. 547--557, 2014.
95 .. [3] A.J. Kinderman and J.F. Monahan, "Computer Generation of Random
96 Variables Using the Ratio of Uniform Deviates",
97 ACM Transactions on Mathematical Software, 3(3), p. 257--260, 1977.
99 Examples
100 --------
101 >>> import numpy as np
102 >>> from scipy import stats
103 >>> rng = np.random.default_rng()
105 Simulate normally distributed random variables. It is easy to compute the
106 bounding rectangle explicitly in that case. For simplicity, we drop the
107 normalization factor of the density.
109 >>> f = lambda x: np.exp(-x**2 / 2)
110 >>> v_bound = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
111 >>> umax, vmin, vmax = np.sqrt(f(0)), -v_bound, v_bound
112 >>> rvs = stats.rvs_ratio_uniforms(f, umax, vmin, vmax, size=2500,
113 ... random_state=rng)
115 The K-S test confirms that the random variates are indeed normally
116 distributed (normality is not rejected at 5% significance level):
118 >>> stats.kstest(rvs, 'norm')[1]
119 0.250634764150542
121 The exponential distribution provides another example where the bounding
122 rectangle can be determined explicitly.
124 >>> rvs = stats.rvs_ratio_uniforms(lambda x: np.exp(-x), umax=1,
125 ... vmin=0, vmax=2*np.exp(-1), size=1000,
126 ... random_state=rng)
127 >>> stats.kstest(rvs, 'expon')[1]
128 0.21121052054580314
130 """
131 if vmin >= vmax:
132 raise ValueError("vmin must be smaller than vmax.")
134 if umax <= 0:
135 raise ValueError("umax must be positive.")
137 size1d = tuple(np.atleast_1d(size))
138 N = np.prod(size1d) # number of rvs needed, reshape upon return
140 # start sampling using ratio of uniforms method
141 rng = check_random_state(random_state)
142 x = np.zeros(N)
143 simulated, i = 0, 1
145 # loop until N rvs have been generated: expected runtime is finite.
146 # to avoid infinite loop, raise exception if not a single rv has been
147 # generated after 50000 tries. even if the expected numer of iterations
148 # is 1000, the probability of this event is (1-1/1000)**50000
149 # which is of order 10e-22
150 while simulated < N:
151 k = N - simulated
152 # simulate uniform rvs on [0, umax] and [vmin, vmax]
153 u1 = umax * rng.uniform(size=k)
154 v1 = rng.uniform(vmin, vmax, size=k)
155 # apply rejection method
156 rvs = v1 / u1 + c
157 accept = (u1**2 <= pdf(rvs))
158 num_accept = np.sum(accept)
159 if num_accept > 0:
160 x[simulated:(simulated + num_accept)] = rvs[accept]
161 simulated += num_accept
163 if (simulated == 0) and (i*N >= 50000):
164 msg = ("Not a single random variate could be generated in {} "
165 "attempts. The ratio of uniforms method does not appear "
166 "to work for the provided parameters. Please check the "
167 "pdf and the bounds.".format(i*N))
168 raise RuntimeError(msg)
169 i += 1
171 return np.reshape(x, size1d)