Example 2.8

In what follows, we will implement rejection methods for sampling continuous distributions.

Let us implement here Example 2.8 of lecture notes.

First let us describe the Beta density.

import numpy as np
import matplotlib.pyplot as plt

def beta_pdf(x, alpha, beta):
    return (x**(alpha-1) * (1-x)**(beta-1) * np.math.gamma(alpha+beta) /
            (np.math.gamma(alpha) * np.math.gamma(beta)))

Now we are ready to implement the rejection method as shown in Example 2.8.

Example 2.8. We aim at sampling from \(p(x) = \text{Beta}(x;\alpha,\beta)\) with \(\alpha = \beta = 2\). We will do this in a way that will demonstrate the Fundamental Theorem of Simulation. As noted in the example, we know that \(p^\star = 1.5\) in this case, which is the maximum of this density. In order to draw samples under the curve \(p(x)\), we will first sample uniformly in the “box” \((X', U') \sim \text{Unif}([0,1] \times [0, p^\star])\) and accept/reject samples if \(U' \leq p(X')\). Note that this sampler is just for demonstration as the rejection samplers.

alpha = 2
beta = 2
p_star = 1.5

rng = np.random.default_rng(5)

xx = np.linspace(0, 1, 1000)

n = 5000

fig, axs = plt.subplots(1, 2, figsize=(10, 4))

# for demonstration purposes, we do this in a for loop
# this can be done all at once

acc_samples = np.array([])
acc_u = np.array([])
# note that we do not have to store uniform samples
# this is again for demonstration of the FTS.

for i in range(n):
    x = rng.uniform(0, 1)
    u = rng.uniform(0, p_star)

    if u <= beta_pdf(x, alpha, beta):
        acc_u = np.append(acc_u, u)
        acc_samples = np.append(acc_samples, x)

axs[0].plot(xx, beta_pdf(xx, alpha, beta), color='k', linewidth=2)
axs[0].add_patch(plt.Rectangle((0, 0), 1, p_star, color=[0.8, 0, 0], fill=False, lw=2))
axs[0].scatter(acc_samples, acc_u, color=[0.8, 0, 0], s=12)
axs[0].set_xlim(0, 1)
axs[0].set_ylim(0, 2.1)
axs[1].cla()
axs[1].plot(xx, beta_pdf(xx, alpha, beta), color=[0.8, 0, 0], linewidth=2)
axs[1].hist(acc_samples, bins=100, density=True, color='k')
axs[1].set_xlim(0, 1)
plt.show()

<ipython-input-1-069f25697792>:5: DeprecationWarning: `np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`
  return (x**(alpha-1) * (1-x)**(beta-1) * np.math.gamma(alpha+beta) /
<ipython-input-1-069f25697792>:6: DeprecationWarning: `np.math` is a deprecated alias for the standard library `math` module (Deprecated Numpy 1.25). Replace usages of `np.math` with `math`
  (np.math.gamma(alpha) * np.math.gamma(beta)))

One can see the samples on the left that are uniformly distributed. On the right, the histogram of the accepted samples is shown (only accepted \(x\) values which are acc_samples in the code). Animating this process is too costly in a notebook but is provided in the course.