Generate random values per a probability density function.

PDF_rand.py

import random
import matplotlib.pyplot as plt

a = 0  # Window minimum in x
b = 2  # Window maximum in x
c = 4  # Window maximum in y
my_df = lambda x : x ** 2


def rejection(f, x, y):
    """
    This function accepts as first argument a univariate function f. This
    function's will be computed at value x, which is the second argument. The
    third argument is y, which is a test value. If f(x) < y, then this will
    return None, otherwise it will return x.
    """
    if f(x) < y:
        return None
    else:
        return x

def random_pair(x_min, x_max, y_min, y_max):
    x_delta = x_max - x_min
    y_delta = y_max - y_min
    rand_f = lambda delta, offset : delta * random.random() + offset
    return (rand_f(x_delta,x_min), rand_f(y_delta, y_min))


def density_function_random(df, x_min, x_max, y_min, y_max):
    """
    Returns random values pairs per a univariate density function.

    Bounding window values for x (independent variable) and y (being f(x)) must
    be provided. Note that y_max should be max(f(x_min <= x <= x_max)) and y_min
    should be min(f(x_min <= x <= x_max)) or the distribution will be some
    combination of slow or perturbed.
    """
    while True:
        x, y = random_pair(x_min, x_max, y_min, y_max)
        if rejection(df, x, y) is not None:
            return x

def get_n_pdf_values(n, df):
    return [density_function_random(df, a, b, 0 ,c) for i in range(n)]


vals = get_n_pdf_values(10000, my_df)
plt.hist(vals,bins=60)
plt.show()