# What is this random number generation algorithm and why does it work (or does it)?

Found the following (Monte Carlo) RNG algorithm in a certain article. Let $f(x),f:\mathbb{R}^n \rightarrow \mathbb{R}$ be the function from which the samples are desired.

1. Draw a random point $x \in \Omega$ (uniform distribution in $\mathbb{R}^n$), where $\Omega$ is the desired region $\Omega \subset \mathbb{R}^n$.
2. Evaluate the probability $p(x) = \int_{\omega(x)} f(x′) d x′$, where $\omega(x) \subset \Omega$ is a finite but small integration region around $x$, similar (with the same volume) for each drawn $x$.
3. Draw a random number $0 \leq q \leq 1$ (uniform distribution for $q \in \mathbb{R}$). The point $x$ is included in the sample set $S$ if $q \leq p(x)/p_m$, where $p_m = \max_{x \in \Omega} p(x)$.
4. Return to 1 until $S$ is large enough.

In the source is not explained how $\omega(x)$ should be chosen and we must assume at least that $f(x)>0$ only in compact set $\Omega$ (and zero elsewhere) and that $f$ is bounded from above.

So I have not designed the algorithm. I just want to figure out if the algorithm is valid method to generate random samples from the distribution $f$ (and under what conditions).

I haven't been thus far able to rigorously prove that the algorithm works but I have made the following attempt at solution:

Assume that $f$ is some density in $\Omega = [a,b] \subset \mathbb{R}$ and for each $x\in\Omega$ let $\omega(x)=[x-\epsilon,x+\epsilon]$ for some small $\epsilon > 0$.

Then $p(x) = \int_{x-\epsilon}^{x+\epsilon}f(x')dx' \approx \int_{x-\epsilon}^{x+\epsilon}f(x)dx' = 2\epsilon f(x)$ (and this holds accurately if $f$ is close to constant near $x$) and we accept point $x$ if $q \leq \frac{2\epsilon f(x)}{p_m}$ (and we have that $q \sim \text{Uniform}([0,1])$). If we now compare this to the usual rejection sampling where the instrumental distribution is uniform i.e. $g(x)$ is uniform pdf, we see that the difference is the additional multiplication by $2\epsilon$. Thus it seems that the algorithm is valid in some approximative sense as long as $2\epsilon \leq 1$ and the boundaries do not cause issues (i.e. the density goes to zero long enough before the boundaries $x=a,x=b$). Also if $\epsilon$ is chosen very small then it is obvious that the acceptance probability will be very small making the algorithm even more inefficient than what rejection sampling might be.

So is the above algorithm some well-known variant of accept-rejection method? Could someone provide a proof or link to a proof where it is explained if the algorithm works or not.

• which article ? Commented Aug 18, 2015 at 12:43
• @Glen_b I won't share the article here and that's not necessary in any case. The above algorithm is just mentioned in it and just like I have written but I'm not sure if the algorithm is good or not. Commented Aug 18, 2015 at 13:00
• Unfortunately, people often misquote (or completely misunderstand) their sources in questions like this, and that's why it's important to include a reference to the source. Commented Aug 18, 2015 at 13:10
• @BrianBorchers ok, I can understand that. However, my question in any case is not whether the algorithm is correct or not for some specific case in some paper but whether the algorithm as stated above is a valid general RNG algorithm for any probability density (with reasonable assumptions for $f$, i.e. no point masses) or not? Commented Aug 18, 2015 at 13:59
• Did you look at TAOCP Vol 2 already? If it's not there, then come back. Commented Aug 18, 2015 at 17:46

## 1 Answer

Your description of the method is a bit odd. Your $p(x)$ is just $f(x)$ scaled by a factor of $2\epsilon$ with a bit of smoothing of $f(x)$. It isn't clear why $f(x)$ needs to be smoothed.

It also isn't clear why you want to generate a large sample of random values so that the emprical distribution of the sample matches $f(x)$. It's more common to simply generate as many samples from the probability distribution as required.

However, the method you've described is of the general class of acceptance-rejection methods. These are discussed in many textbooks on simulation and random variate generation. For example, see Averill M. Law's Simulation Modeling & Analysis textbook.

• Thank you for the observations. Unfortunately I think you have slightly misunderstood my point here. I have edited the question now. So I'm just trying to understand if the algorithm works or not. It is indeed quite similar to accept-rejection method (as I show in my hand-waving analysis) but I'm not sure. I have consulted some books (e.g. Robert and Casella: Monte Carlo Statistical Methods) but I haven't found the very same algorithm there. Unfortunately I can't access the book you linked to me. Commented Aug 18, 2015 at 18:07
• The quantity $p(x)$ in this algorithm is essentially $2\epsilon f(x)$ while $p_{m}$ is essentially $2\epsilon$ times the maximum value of $f(x)$ over the support of the distribution. Thus the ratio $2\epsilon f(x)/p_{m}$ is essentially $f(x)/f_{m}$, which is exactly the factor used in a conventional acceptance-rejection method. Your method is potentially inaccurate if $p(x)$ isn't close to $2\epsilon f(x)$ because $\epsilon$ is too large. Thus you should stick to the conventional use of $f(x)/f_{m}$. Commented Aug 18, 2015 at 19:59
• To answer your revised question- the algorithm is an odd variant of the well known acceptance-rejection method that will fail if $\epsilon$ is too large. Commented Aug 18, 2015 at 21:07
• You are right, I mixed up $p_m = \max_{x\in\Omega}p(x)$ and $\max_{x\in\Omega}f(x)$ in my analysis and that's why didn't completely figured out it myself and got confused! Anyway, this is to clear to me now, I will mark your answer as accepted. Commented Aug 18, 2015 at 21:30
• and, indeed there seems to no reason to use this algorithm instead of the usual acceptance-rejection method. It just causes extra work as one needs to deal with the choice of $\omega(x)$ and worry about the potential approximation error and the choice of $\epsilon$. Commented Aug 18, 2015 at 21:34