Create samples from trapezoidal probability density function

Question

I have defined my probability density function as follows:

def trapezoidal_pdf(x, a, b, c, d):
    const = (2 / (d + c - a - b))
    
    if a <= x < b:
        probability = const * (x - a) / (b - a)
        return probability
    elif b <= x < c:
        probability = const
        return probability
    elif c <= x < d:
        probability = const * (d - x) / (d - c)
        return probability
    else:
        return 0.0  # Outside the defined range, probability is 0

How do I get samples from the pdf so when I plot the histogram it looks like a trapezoid?

I have tried using np.random.random, but I think that's not what I want

Answer 1

ok, let me do a quick try, following https://en.wikipedia.org/wiki/Trapezoidal_distribution.

You could look at PDF (and CDF) of the distribution, and there are three distinctive parts - first triangle, middle box, last triangle. And you could sample in the first step what part you sample from, and then sample it from part distribution. Middle part is easy, just a scaled uniform, and side triangles are not hard either

Along the lines, not optimized, Python 3.11 Win x64

import numpy as np
import matplotlib.pyplot as plt

def select_part(abcd, rv: np.float64) -> int:
    a, b, c, d = abcd
    
    den = 1.0/(d + c - a - b)
    
    a1 = (b-a)*den
    a2 = 1.0 - (d - c)*den
    a3 = 1.0
    
    if rv <= a1:
        return 1
    
    if rv <= a2:
        return 2
    
    if rv <= a3:
        return 3
    
    return -1 # just in case, should never happen

def sample_trap(abcd, rv1: np.float64, rv2: np.float64) -> np.float64:
    a, b, c, d = abcd
    
    part = select_part(abcd, rv1)
    
    if part == 1:
        return a + np.sqrt(rv2)*(b - a)
    
    if part == 2:
        return b + (c-b)*rv2
    
    if part == 3:
        return d - np.sqrt(rv2)*(d - c)
    
    return -999999.0 
        
rng = np.random.default_rng(135797531)

abcd = (-2., 0., 1., 2.)

N = 1000000

res = np.empty(N)

for k in range(0, N):
    rv1 = rng.random()
    rv2 = rng.random()
    res[k] = sample_trap(abcd, rv1, rv2)

fig, axs = plt.subplots(1, 2, tight_layout=True)

axs[0].hist(res, bins=50)

and you should get something like