import matplotlib.pyplot as plt
import numpy as np

def p(x):
    return 6*x*(1-x)

N = 1000000

np.random.seed(42)
x = np.random.uniform(0, 1, N)
y = np.random.uniform(0, 1.5, N)

z = x[y < p(x)]		# this piece of numpy magic returns a subarray of
                    # x containing those values where y[i] < p(x[i])

nb = 50
plt.hist(z, nb)
sum = len(z)
plt.plot(z, p(z)*sum/nb, 'b.', markersize=1)
plt.show()