## 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] 500
```

If 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.191538
```

In 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.267563256
```

**Important**:

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.

**Using unnecessary coercion may impose a certain computational cost, depending on the size of the vector of observations generated, for example, in a simulation study.**