Data Science
Monte Carlo Simulation
Rejection Sampling
Statistics
Python

[Statistics by Example]: Rejection Sampling

A step-by-step walk through with math and python code

Garrett Hurst

6 min readJan 25, 2025

In Data Science, simulation of a random variable is highly valuable: it is fundamental to Monte Carlo Methods that pop up in finance, queueing, reinforcement learning, and many more applications. However, the distributions that we are interested in are often complex and hard to sample. Rejection Sampling allows us to reduce the complexity of simulating a hard distribution.

Motivation:

Suppose that we have a Probability Density Function (PDF) f(x) as defined below that we would like to generate samples from:

# Note here that we did not include the integration constant divisor.
# Shortly, we will apply a little trick to `cancel this out`
def f_target(x):
    return np.exp(-x ** 4) * np.sin(10 * x) ** 2

integration_constant = scipy.integrate.quad(f_target, -np.inf, np.inf)[0]

# plot the function. 
x = np.linspace(-3, 3, 500)
f_x = f_target(x) / integration_constant
plt.title('Desired Distribution')
plt.ylabel('f(x)')
plt.xlabel('x')
plt.plot(x, f_x)
plt.show()

In a basic statistics course, it is taught that this can be done with the Inverse Transform Method, and while this can be done, it is very challenging for several reasons.

Calculating the normalization constant in the denominator takes a lot of work.
Finding the CDF also takes a lot of work
Once the CDF has been found, inverting it will also take a LOT of work.

This can be done, but it is challenging and computationally expensive. Just to hit the point home, this example is one dimensional, and as Data Scientists, we rarely care about just one dimensional data. When this problem is scaled to the n’th dimension, it gets very computationally expensive to make these calculations.

Idea:

Since we know the distribution of f(x), we can find another distribution g(x) that contains f(x). Then, if samples are taken from g(x), the ones that do not belong to f(x) can be rejected.

Algorithm:

Propose a distribution, g(x), and a constant, C, such that f(x) ≤ C g(x)for all x
Collect i.i.d. samples from the proposed distribution X ~ g(x)
Collect i.i.d. samples U ~ Uniform[0,1]
Accept the samples that meet the condition:

Selecting g(x) and C:

In order to generate samples under all parts of f(x), the g(x) that we propose must be greater than f(x) in all locations from negative infinity to positive infinity (or just the portion of the distribution that you care about sampling from). Thus, the below condition must be sufficed:

Source: Author

In our example from above, it can be seen that f(x) exponentially decays towards zero as |x| grows. Thus, a reasonable g(x) to propose would be an exponential function. We’ll use:

Source: Author

Now, we must find the constant C such that:

Recall that the normalization integral constant in the denominator of f(x) was one of the things that was challenging to compute. Since we have the freedom to choose whatever C we want to suffice this condition, let’s write our C in terms of another constant, D, that will cancel out this nasty integral.

Now, if we let C = D / Z,

This new constant, D, is much easier to solve for. From a simple inspection of f(x), we can see that since sin(10x)² takes values [0,1], it will only reduce the value of the exponential term. Thus, we can simplify the problem to:

Source: Author

Now, since x⁴ grows faster than x for all values greater than 1, we can reduce the area to search to D to the interval -1 < x < 1. Then, we can reduce this further to 0 ≤ x < 1 because the absolute value is taken on both sides. After a little rearrangement, we get:

Source: Author

Thus, by maximizing -x⁴ + |x| on the interval 0 ≤ x < 1, we can find that:

Source: Author

When plotting this, we can see that the distribution g(x) fully covers f(x). In the code below, C was calculated for plotting purposes. However, in practice only D is needed since the integral normalization constant cancels out (remember, our goal is to avoid making this integral calculation).

# scaling constants
D = np.exp(-(1/4)**(4/3) + (1/4)**(1/3))
C = D / integration_constant

# exponential distribution mirrored with abs value
def g_selected(x):
    return scipy.stats.expon.pdf(np.abs(x))

g_x = C * g_selected(x)
plt.title('Desired Distribution Under Proposal Distribution')
plt.ylabel('pdf value')
plt.xlabel('x')
plt.plot(x, f_x, label='Desired f(x)')
plt.plot(x, g_x, label='Proposal g(x)')
plt.legend()
plt.show()

Generating Samples:

Now that we have g(x), and our modified constant D all that is left to do is plug things into the algorithm from above. The below code generates a sample of 100,000 i.i.d. samples from f(x)

n = 100000
n_generated = 0
samples_f = []
while len(samples_f) < n:
    # generate U~uniform[0,1]
    uniform = np.random.uniform(0,1)
    # generate a X~Exp(|x|) (multiplying by -1 at random ror abs value in pdf)
    sample  = np.random.exponential() * np.random.choice([-1, 1])
    # count number of accepted samples to get the acceptance rate
    n_generated += 1
    # check if the sample should be rejected
    if uniform <= f_target(sample) / (D * g_selected(sample)):
        samples_f.append(sample)


print(f'Total Samples Generated: {n_generated}')
print(f'Samples Accepted: {n} ({round(n/n_generated * 100,2)} %)')

plt.title('Distribution of Samples Generated')
plt.ylabel('pdf value')
plt.xlabel('x')
plt.plot(x, f_x, label='Desired Distribution f(x)')
plt.plot(x, g_x, label='Proposal  Distribution g(x)')
plt.hist(samples_f, density=True, alpha=0.75, bins=500, label='Sample Distribution')
plt.legend()
plt.show()

Total Samples Generated: 353750
Samples Accepted: 100000 (28.27 %)

Notice that the acceptance rate for the samples generated was ~30%. This is dependent on the g(x) and C values that are selected. A ( C g(x) ) that is able to more tightly enclose f(x) will have a higher acceptance rate. However, choosing these values is a combination of an art and a science. This example could have achieved similar (or arguably better) results by choosing g(x) to be Gaussian with an optimal sigma, or by choosing g(x) to be be Uniform and cutting off the tails of the original f(x). Our goal here is to save the work, so it is up to you to choose a g(x) that is a balance between easy to simulate while maintaining an acceptable acceptance rate.