从梯形概率密度函数创建样本

Question

我的概率密度函数定义如下：

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

如何从 pdf 中获取样本，以便在绘制直方图时它看起来像梯形？

我尝试过使用

np.random.random

，但我认为这不是我想要的

Answer 1

好的，让我快速尝试一下，如下 https://en.wikipedia.org/wiki/Trapezoidal_distribution。

你可以看一下分布的 CDF，它有三个不同的部分 - 第一个三角形、中间的盒子、最后一个三角形。您可以在第一步中对从哪个部分进行采样，然后从部件分布中进行采样。中间部分很简单，只是一个缩放制服，侧面三角形也不难

沿着这条线，未优化，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 (b - a)*np.sqrt(rv2) + a
    
    if part == 2:
        return (c-b)*rv2 + b
    
    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)