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