Exercise: Rejection Sampling

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm

%matplotlib inline
mu_1 = -2.5
sigma_square_1 = 4.
sigma_1 = np.sqrt(sigma_square_1)

mu_2 = 3.5
sigma_square_2 = 4.
sigma_2 = np.sqrt(sigma_square_2)
prob_1 = 0.4
size=40
rv_1 = norm(loc = mu_1, scale = sigma_1)
rv_2 = norm(loc = mu_2, scale = sigma_2)
x_ = np.arange(-14, 16, .1)

p_green = lambda x: prob_1 * rv_1.pdf(x) + (1-prob_1) * rv_2.pdf(x)
plt.plot(x_, p_green(x_) , "g-")


sigma_red = 5
mu_red = 1.
q_red = norm(loc = mu_red, scale = sigma_red)

plt.plot(x_, q_red.pdf(x_) , "r-")

_ = plt.xlabel("x")

Assume that we can't sample from the green curve $ p(x) $ , but we can sample from the red one $ q(x) $ ).

And we can compute $ p(x) $ at each position $ x $ .

Exercise: Implement rejection sampling to get a sample from $ p(x) $ by

# x-values of red points are rejected
# x-values of green point are accepted
plt.plot(x_, p_green(x_) , "g-")
plt.plot(x_, f_red(x_)  , "r-")
plt.plot(samples, y_accept, 'g.', label='Samples')
plt.plot(rejected_samples, y_reject, 'r.', label='Samples')
_ = plt.hist(samples, bins=25, density=True)
plt.plot(x_, p_green(x_) , "g-")