#!/usr/bin/python3
#
# Rejection sampling

import numpy as np

# P is some PDF supported on [Dmin, Dmax] to be sampled. 
# c is a constant > max(P).
# Returns a single number drawn from P.
def rejection_sampling(P, c, Dmin, Dmax):
    while True:
        # draw x uniformly from [Dmin, Dmax):
        x = (Dmax - Dmin) * np.random.random() + Dmin
        # accept with probability P(x) / c; reject otherwise
        if P(x) / c > np.random.random():
           return x
           
# test with the PDF P(x) = sin(x)/2 on [0, pi):
if __name__ == '__main__':

    import matplotlib.pyplot as plt

    def P(x):
        if x >= 0.0 and x < np.pi:
            return 0.5 * np.sin(x)
        return 0.0
    
    N = 10000
    sample = np.zeros(N)
    for i in range(N):
        sample[i] = rejection_sampling(P, 0.5, 0.0, np.pi)
    plt.hist(sample, bins=30)
    plt.show()