Understanding the Method
In situations where we cannot use the inversion method (situations where obtaining the quantile function is not possible) and neither do we know a transformation involving a random variable from which we can generate observations, we can make use of the acceptance-rejection method.
Suppose \(X\) and \(Y\) are random variables with probability density function (pdf) or probability function (pf) \(f\) and \(g\), respectively. Furthermore, suppose there exists a constant \(c\) such that
\[\frac{f(x)}{g(y)} \leq c,\] for every value of \(x\), with \(f(x) > 0\). To use the acceptance-rejection method to generate observations of the random variable \(X\), using the algorithm below, first find a random variable \(Y\) with pdf or pf \(g\), such that it satisfies the above condition.
Important:
It is important that the chosen random variable \(Y\) is such that you can easily generate its observations. This is because the acceptance-rejection method is computationally more intensive than more direct methods such as the transformation method or the inversion method, which only requires the generation of pseudo-random numbers with a uniform distribution.
Algorithm of the Acceptance-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 go back to step 1.
Proof: Consider the discrete case, that is, \(X\) and \(Y\) are random variables with pfs \(f\) and \(g\), respectively. By step 3 of the algorithm above, we have \(\{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 Law of Total Probability, we have:
\[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-rejection method, we accept the occurrence of \(Y\) as an occurrence of \(X\) with probability \(1/c\). Moreover, by Bayes’ Theorem, we have
\[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 algorithm’s procedure is equivalent to accepting a value from \(X\) that has pf \(f\). For the continuous case, the proof is similar.
Important:
Notice that to reduce the computational cost of the method, we should choose \(c\) in such a way that we can maximize \(P(accept)\). Therefore, choosing an excessively large value of the constant \(c\) will reduce the probability of accepting an observation from \(Y\) as an observation of the random variable \(X\).
Note:
Computationally, it is convenient to consider \(Y\) as a random variable with a uniform distribution on the support of \(f\), since generating observations from a uniform distribution is straightforward on any computer. For the discrete case, considering \(Y\) with a discrete uniform distribution might be a good alternative.
Installation and loading the package
The AcceptReject package is available on CRAN and can be
installed using the following command:
install.packages("AcceptReject")
# or
install.packages("remotes")
remotes::install_github("prdm0/AcceptReject", force = TRUE)
# Load the package
library(AcceptReject)Using the accept_reject Function
Among various functions provided by the AcceptReject library, the
acceptance_rejection function implements the acceptance-rejection
method. The accept_reject() function has the following
signature:
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,
...
)Many of the arguments the user will not need to change, as the
AcceptReject::accept_reject() function already has default
values for them. However, it is important to note that the f argument is
the probability density function (pdf) or probability function (pf) of
the random variable \(X\) from which
observations are desired to be generated. The args_f
argument is a list of arguments that will be passed to the f function.
The c argument is the value of the constant c
that will be used in the acceptance-rejection method. If the user does
not provide a value for c, the accept_reject()
function will calculate the value of c that maximizes the
probability of accepting observations from \(Y\) as observations from \(X\).
Note:
accept_reject()
function can work in a parallelized manner on Unix-based operating
systems. If you use operating systems such as Linux or MacOS, you can
benefit from parallelization of the accept_reject()
function. To do this, simply set the parallel = TRUE
argument. In Windows, parallelization is not supported, and setting the
parallel = TRUE argument will have no effect.
You do not need to define the c argument when using the
accept_reject() function. By default, if
c = NULL, the accept_reject() function will
calculate the value of c that maximizes the probability of
accepting observations from \(Y\) as
observations from \(X\). However, if
you want to set a value for c, simply pass a value to the
c argument.
Details of the optimization of c:
c can be a
difficult optimization process, but generally, it is not. Therefore,
unless you have reasons to set a value for c, it is
recommended to use the default value c = NULL. For very
complicated functions, you may choose a sufficiently large
c to ensure that the method works well.
The ... argument allows changing the desired precision
in the optimize() function, used to optimize the value of
the constant c. By default,
tol = .Machine$double.eps^0.25 is used.
Examples
Below are some examples of using the accept_reject()
function to generate pseudo-random observations of discrete and
continuous random variables. It should be noted that in the case of
\(X\) being a discrete random variable,
it is necessary to provide the argument continuous = FALSE,
whereas in the case of \(X\) being
continuous (the default), you must consider
continuous = TRUE.
Generating discrete observations
As an example, let \(X \sim Poisson(\lambda
= 0.7)\). We will generate \(n =
1000\) observations of \(X\)
using the acceptance-rejection method, using the
accept_reject() function. Note that it is necessary to
provide the xlim argument. Try to set an upper limit value
for which the probability of \(X\)
assuming that value is zero or very close to zero. In this case, we
choose xlim = c(0, 20), where
dpois(x = 20, lambda = 0.7) is very close to zero
(1.6286586^{-22}).
library(AcceptReject)
#>
#> Attaching package: 'AcceptReject'
#> The following object is masked from 'package:stats':
#>
#> qqplot
# Ensuring Reproducibility
set.seed(0)
# Generating observations
data <- AcceptReject::accept_reject(
n = 1000L,
f = dpois,
continuous = FALSE,
args_f = list(lambda = 0.7),
xlim = c(0, 20),
parallel = FALSE
)
# Viewing organized output with useful information
print(data)
#>
#> ── Accept-Reject Samples ───────────────────────────────────────────────────────
#> ℹ It's not necessary, but if you want to extract the observations, use as.vector().
#> ✔ Case: discrete
#> ✔ Number of observations: 1000
#> ✔ c: 9.9317
#> ✔ Probability of acceptance (1/c): 0.1007
#> ✔ Observations: 0 1 0 1 0 1 0 0 0 0...
#> ✔ xlim = 0 20
#>
#> ────────────────────────────────────────────────────────────────────────────────
# Calculating the true probability function for each observed value
values <- unique(data)
true_prob <- dpois(values, lambda = 0.7)
# Calculating the observed probability for each value in the observations vector
obs_prob <- table(data) / length(data)
# Plotting the probabilities and observations
plot(values, true_prob, type = "p", pch = 16, col = "blue",
xlab = "x", ylab = "Probability", main = "Probability Function")
# Adding the observed probabilities
points(as.numeric(names(obs_prob)), obs_prob, pch = 16L, col = "red")
legend("topright", legend = c("True probability", "Observed probability"),
col = c("blue", "red"), pch = 16L, cex = 0.8)
grid()
Note that it is necessary to specify the nature of the random
variable from which observations are desired to be generated. In the
case of discrete variables, the argument continuous = FALSE
must be passed.
Now, consider that we want to generate observations from a random variable \(X \sim Binomial(n = 5, p = 0.7)\). Below, we will generate \(n = 2000\) observations of \(X\).
library(AcceptReject)
# Ensuring reproducibility
set.seed(0)
# Generating observations
data <- AcceptReject::accept_reject(
n = 2000L,
f = dbinom,
continuous = FALSE,
args_f = list(size = 5, prob = 0.5),
xlim = c(0, 20),
parallel = FALSE
)
# Viewing organized output with useful information
print(data)
#>
#> ── Accept-Reject Samples ───────────────────────────────────────────────────────
#> ℹ It's not necessary, but if you want to extract the observations, use as.vector().
#> ✔ Case: discrete
#> ✔ Number of observations: 2000
#> ✔ c: 6.25
#> ✔ Probability of acceptance (1/c): 0.16
#> ✔ Observations: 3 1 4 2 3 1 2 4 2 2...
#> ✔ xlim = 0 20
#>
#> ────────────────────────────────────────────────────────────────────────────────
# Calculating the true probability function for each observed value
values <- unique(data)
true_prob <- dbinom(values, size = 5, prob = 0.5)
# Calculating the observed probability for each value in the observations vector
obs_prob <- table(data) / length(data)
# Plotting the probabilities and observations
plot(values, true_prob, type = "p", pch = 16, col = "blue",
xlab = "x", ylab = "Probability", main = "Probability Function")
# Adding the observed probabilities
points(as.numeric(names(obs_prob)), obs_prob, pch = 16L, col = "red")
legend("topright", legend = c("True probability", "Observed probability"),
col = c("blue", "red"), pch = 16L, cex = 0.8)
grid()
Generating continuous observations
To expand beyond examples of generating pseudo-random observations of
discrete random variables, consider now that we want to generate
observations from a random variable \(X \sim
\mathcal{N}(\mu = 0, \sigma^2 = 1)\). We chose the normal
distribution because we are familiar with its form, but you can choose
another distribution if desired. Below, we will generate
n = 2000 observations using the acceptance-rejection
method. Note that continuous = TRUE.
library(AcceptReject)
# Ensuring reproducibility
set.seed(0)
# Generating observations
data <- AcceptReject::accept_reject(
n = 2000L,
f = dnorm,
continuous = TRUE,
args_f = list(mean = 0, sd = 1),
xlim = c(-4, 4),
parallel = FALSE
)
# Viewing organized output with useful information
print(data)
#>
#> ── Accept-Reject Samples ───────────────────────────────────────────────────────
#> ℹ It's not necessary, but if you want to extract the observations, use as.vector().
#> ✔ Case: continuous
#> ✔ Number of observations: 2000
#> ✔ c: 3.1915
#> ✔ Probability of acceptance (1/c): 0.3133
#> ✔ Observations: 1.2864 1.0329 -2.3522 -0.9272 2.1587 -0.0184 1.7409 1.2134 -2.9956 -1.8622...
#> ✔ xlim = -4 4
#>
#> ────────────────────────────────────────────────────────────────────────────────
hist(
data,
main = "Generating Gaussian observations",
xlab = "x",
probability = TRUE,
ylim = c(0, 0.4)
)
x <- seq(-4, 4, length.out = 500L)
y <- dnorm(x, mean = 0, sd = 1)
lines(x, y, col = "red", lwd = 2)
legend("topright", legend = "True density", col = "red", lwd = 2)
In the examples above, the graphs were built without using the
plot.accept_reject() function. This is just to show that
you can manipulate the returning object using the
accept_reject() function, that is, the class object
accept_reject.
However, the plot.accept_reject() function can be used
to generate graphs in a simpler way. Below, an example of how to use the
plot.accept_reject() function to generate the probability
density plot of the normal distribution. However, note that the
plot.accept_reject() function makes the plotting task
simpler and more direct. See the following example:
library(AcceptReject)
library(cowplot) # install.packages("cowplot")
# Ensuring reproducibility
set.seed(0)
simulation <- function(n){
AcceptReject::accept_reject(
n = n,
f = dnorm,
continuous = TRUE,
args_f = list(mean = 0, sd = 1),
xlim = c(-4, 4),
parallel = FALSE
)
}
# Inspecting
a <- plot(simulation(n = 250L))
b <- plot(simulation(n = 2500L))
c <- plot(simulation(n = 25000L))
d <- plot(simulation(n = 250000L))
plot_grid(a, b, c, d, nrow = 2L, labels = c("a", "b", "c", "d"))
See another example, in the discrete case:
library(AcceptReject)
library(cowplot) # install.packages("cowplot")
# Ensuring Reproducibility
set.seed(0)
simulation <- function(n){
AcceptReject::accept_reject(
n = n,
f = dpois,
continuous = FALSE,
args_f = list(lambda = 0.7),
xlim = c(0, 20),
parallel = FALSE
)
}
a <- plot(simulation(25L))
b <- plot(simulation(250L))
c <- plot(simulation(2500L))
d <- plot(simulation(25000L))
plot_grid(a, b, c, d, nrow = 2L, labels = c("a", "b", "c", "d"))
Accessing metadata
The accept_reject() function returns an object of class
accept_reject. When executing the print()
function on an object of this class, an organized output will be shown.
However, you can operate on this instance of the
accept_reject class as any atomic vector. In the example
below, notice that you can obtain a histogram with the
hist() function or check the size of the vector of
observations generated using the length() function.
library(AcceptReject)
data <- accept_reject(
n = 1000L,
f = dnorm,
continuous = TRUE,
args_f = list(mean = 0, sd = 1),
xlim = c(-4, 4)
)
# Creating a histogram
hist(data)
# Checking the size of the vector of observations
length(x)
#> [1] 500If you want to access some metadata, use the attr()
function. Check the list of attributes by doing:
library(AcceptReject)
data <- accept_reject(
n = 100L,
f = dnorm,
continuous = TRUE,
args_f = list(mean = 0, sd = 1),
xlim = c(-4, 4)
)
attributes(data)
#> $dim
#> [1] 100 1
#>
#> $class
#> [1] "accept_reject"
#>
#> $f
#> <partialised>
#> function (...)
#> f(mean = 0, sd = 1, ...)
#> <environment: 0x5fb29195ae58>
#>
#> $args_f
#> $args_f$mean
#> [1] 0
#>
#> $args_f$sd
#> [1] 1
#>
#>
#> $c
#> [1] 3.191538
#>
#> $continuous
#> [1] TRUE
#>
#> $xlim
#> [1] -4 4
# Accessing the value c
attr(data, "c")
#> [1] 3.191538In any case, it is important to highlight that, in general, you will
not need to access these attributes. The greatest interest will be in
having access to the vector of observations generated. If you want to
access the observation values directly in an atomic vector in R without
attributes, without an organized printout, simply coerce the object
using the as.vector() function, as shown in the following
example:
library(AcceptReject)
data <- accept_reject(
n = 100L,
f = dnorm,
continuous = TRUE,
args_f = list(mean = 0, sd = 1),
xlim = c(-4, 4)
)
class(data)
#> [1] "accept_reject"
print(data)
#>
#> ── Accept-Reject Samples ───────────────────────────────────────────────────────
#> ℹ It's not necessary, but if you want to extract the observations, use as.vector().
#> ✔ Case: continuous
#> ✔ Number of observations: 100
#> ✔ c: 3.1915
#> ✔ Probability of acceptance (1/c): 0.3133
#> ✔ Observations: -0.8843 1.2379 0.7305 -0.021 -0.5377 0.0052 -0.5952 1.266 -0.0175 -0.3715...
#> ✔ xlim = -4 4
#>
#> ────────────────────────────────────────────────────────────────────────────────
# Coercing the object into an atomic vector without attributes
data <- as.vector(data)
print(data)
#> [1] -0.884300686 1.237916782 0.730452033 -0.020969238 -0.537684873
#> [6] 0.005160321 -0.595158529 1.265968353 -0.017486773 -0.371506846
#> [11] 0.862292295 0.672636049 -0.008881878 -0.004335830 1.799664292
#> [16] -0.144384425 -0.375649925 1.414700381 -0.142955763 -0.512890248
#> [21] 1.214014906 -1.426491648 -0.335570900 -0.789570918 2.328796349
#> [26] -1.935499994 -0.402352029 -0.645955391 -0.619398503 -0.996748215
#> [31] 0.089246217 -0.139643574 -0.111617487 1.344497407 1.071883820
#> [36] -0.818856794 -0.987588339 0.570187641 -0.228377435 0.832923556
#> [41] -0.488165969 -0.301001403 0.672187231 1.241363853 0.958037097
#> [46] -0.005516769 -0.773214558 2.338708282 -0.246145425 0.030920289
#> [51] -0.195434039 1.434930967 -0.055245180 -0.476447430 -0.245117888
#> [56] 1.703198925 -0.185116883 -0.732147191 1.833668537 0.220142737
#> [61] 1.103287160 0.380271858 -0.957886273 -2.288235800 -1.741464285
#> [66] 0.793730455 -0.168553168 -0.575795487 -0.618257321 -1.186362244
#> [71] 0.096642327 0.651305571 -0.683857242 -0.277836451 1.187555924
#> [76] -1.789664479 -0.975050317 -0.027773783 -1.749138258 1.365316559
#> [81] -0.590321852 -0.951668080 1.729699053 -2.679250218 -1.679565914
#> [86] -0.483523101 1.079006450 -2.266770320 1.604207372 -0.337534644
#> [91] -1.033202954 -1.062941184 0.926517019 -0.769777745 0.312817041
#> [96] -1.398557486 1.643004848 0.425445648 0.238909276 -0.267563256Important:
You will not need to coerce the object of the
accept_reject class into an atomic vector with no
attributes unless you have a specific reason to do so. The object of the
accept_reject class is an atomic vector with attributes,
and you can operate on it like any atomic vector. Everything you can do
with an atomic vector, you can do with an object of the
accept_reject class.