Acceptance-Rejection Method

This function implements the acceptance-rejection method for generating random numbers from a given probability density function (pdf).

Usage

accept_reject(
  n = 1L,
  continuous = TRUE,
  f = NULL,
  args_f = NULL,
  f_base = NULL,
  random_base = NULL,
  args_f_base = NULL,
  xlim = NULL,
  c = NULL,
  parallel = FALSE,
  cores = NULL,
  warning = TRUE,
  ...
)

Arguments

n: The number of random numbers to generate.
continuous: A logical value indicating whether the pdf is continuous or discrete. Default is TRUE.
f: The probability density function (continuous = TRUE), in the continuous case or the probability mass function, in the discrete case (continuous = FALSE).
args_f: A list of arguments to be passed to the f function. It refers to the list of arguments of the target distribution.
f_base: Base probability density function (for continuous case).If f_base = NULL, a uniform distribution will be used. In the discrete case, this argument is ignored, and a uniform probability mass function will be used as the base.
random_base: Random number generation function for the base distribution passed as an argument to f_base. If random_base = NULL (default), the uniform generator will be used. In the discrete case, this argument is disregarded, and the uniform random number generator function will be used.
args_f_base: A list of arguments for the base distribution. This refers to the list of arguments that will be passed to the function f_base. It will be disregarded in the discrete case.
xlim: A vector specifying the range of values for the random numbers in the form c(min, max). Default is c(0, 100).
c: A constant value used in the acceptance-rejection method. If NULL, c will be estimated automatically.
parallel: A logical value indicating whether to use parallel processing for generating random numbers. Default is FALSE.
cores: The number of cores to be used in parallel processing. Default is NULL, i.e, all available cores.
warning: A logical value indicating whether to show warnings. Default is TRUE.
...: Additional arguments to be passed to the optimize(). With this argument, it is possible to change the tol argument of optimize(). Default is tol = .Machine$double.eps^0.25).

Value

A vector of random numbers generated using the acceptance-rejection method. The return is an object of class accept_reject, but it can be treated as an atomic vector.

Details

In situations where we cannot use the inversion method (situations where it is not possible to obtain the quantile function) and we do not know a transformation that involves a random variable from which we can generate observations, we can use the acceptance and rejection method. Suppose that $X$ and $Y$ are random variables with probability density function (pdf) or probability function (pf) $f$ and $g$, respectively. In addition, suppose that there is a constant $c$ such that

$$f(x) \leq c \cdot g(x), \quad \forall x \in \mathbb{R}.$$

for all values of $t$, with $f(t)>0$. To use the acceptance and rejection method to generate observations from the random variable $X$, using the algorithm below, first find a random variable $Y$ with pdf or pf $g$, that satisfies the above condition.

Algorithm of the Acceptance and Rejection Method:

1 - Generate an observation $y$ from a random variable $Y$ with pdf/pf $g$;

2 - Generate an observation $u$ from a random variable $U\sim \mathcal{U} (0, 1)$;

3 - If $u < \frac{f(y)}{cg(y)}$ accept $x = y$; otherwise reject $y$ as an observation of the random variable $X$ and return to step 1.

Proof: Let's consider the discrete case, that is, $X$ and $Y$ are random variables with pf's $f$ and $g$, respectively. By step 3 of the above algorithm, we have that ${accept} = {x = y} = u < \frac{f(y)}{cg(y)}$. That is,

$P(accept | Y = y) = \frac{P(accept \cap {Y = y})}{g(y)} = \frac{P(U \leq f(y)/cg(y)) \times g(y)}{g(y)} = \frac{f(y)}{cg(y)}.$

Hence, by the Total Probability Theorem, we have that:

$P(accept) = \sum_y P(accept|Y=y)\times P(Y=y) = \sum_y \frac{f(y)}{cg(y)}\times g(y) = \frac{1}{c}.$

Therefore, by the acceptance and rejection method we accept the occurrence of $Y$ as being an occurrence of $X$ with probability $1/c$. In addition, by Bayes' Theorem, we have that

$P(Y = y | accept) = \frac{P(accept|Y = y)\times g(y)}{P(accept)} = \frac{[f(y)/cg(y)] \times g(y)}{1/c} = f(y).$

The result above shows that accepting $x = y$ by the procedure of the algorithm is equivalent to accepting a value from $X$ that has pf $f$.

The argument c = NULL is the default. Thus, the function accept_reject() estimates the value of c using the optimization algorithm optimize() using the Brent method. For more details, see optimize() function. If a value of c is provided, the function accept_reject() will use this value to generate the random observations. An inappropriate choice of c can lead to low efficiency of the acceptance and rejection method.

In Unix-based operating systems, the function accept_reject() can be executed in parallel. To do this, simply set the argument parallel = TRUE. The function accept_reject() utilizes the parallel::mclapply() function to execute the acceptance and rejection method in parallel. On Windows operating systems, the code will not be parallelized even if parallel = TRUE is set.

For the continuous case, a base density function can be used, where the arguments f_base, random_base and args_f_base need to be passed. If at least one of them is NULL, the function will assume a uniform density function over the interval xlim.

For the discrete case, the arguments f_base, random_base and args_f_base should be NULL, and if they are passed, they will be disregarded, as for the discrete case, the discrete uniform distribution will always be considered as the base. Sampling from the discrete uniform distribution has shown good performance for the discrete case.

The advantage of using parallelism in Unix-based systems is relative and should be tested for each case. Significant improvement is observed when considering parallelism for very large values of n. It is advisable to conduct benchmarking studies to evaluate the efficiency of parallelism in a practical situation.

References

BISHOP, Christopher. 11.4: Slice sampling. Pattern Recognition and Machine Learning. Springer, 2006.

Brent, R. (1973) Algorithms for Minimization without Derivatives. Englewood Cliffs N.J.: Prentice-Hall.

CASELLA, George; ROBERT, Christian P.; WELLS, Martin T. Generalized accept-reject sampling schemes. Lecture Notes-Monograph Series, p. 342-347, 2004.

NEUMANN V (1951). “Various techniques used in connection with random digits.” Notes by GE Forsythe, pp. 36–38.

NEAL, Radford M. Slice sampling. The Annals of Statistics, v. 31, n. 3, p. 705-767, 2003.

Examples

set.seed(0) # setting a seed for reproducibility

x <- accept_reject(
  n = 2000L,
  f = dbinom,
  continuous = FALSE,
  args_f = list(size = 5, prob = 0.5),
  xlim = c(0, 5)
)
#> ! Warning: f(5) is 0.03125. If f is defined for x >= 5, trying a upper limit might be better.
plot(x)


y <- accept_reject(
  n = 1000L,
  f = dnorm,
  continuous = TRUE,
  args_f = list(mean = 0, sd = 1),
  xlim = c(-4, 4)
)
plot(y)